c===========================================================
c     wave: Solves wave equation:
c
c        u(x,t)_tt = u_xx 
c
c     on x = [0..1], t > 0 with initial conditions
c
c        u(x,0)   = exp(-((x-0.5)/0.1)^2)
c        u_t(x,0) = 0
c
c     and boundary conditions 
c
c        u(0,t) = u(1,t) = 0
c
c     Solution is obtained using method of lines, with 
c     O(h^2) approximation for u_xx, and LSODA to integrate 
c     resulting set of ODEs.
c
c     Output is in form suitable for surface-plotting via 
c     gnuplot.
c
c     Program also uses 'xvs' interface to generate 
c     .sdf files which can subsequently be visualized using
c     'xvs' visualization server.  See links in course 
c     Software page for more details.
c===========================================================
      program        wave

      implicit       none

      integer        iargc,   i4arg
      real*8         r8arg

c-----------------------------------------------------------
c     Include common block for communication with fcn
c-----------------------------------------------------------
      include       'fcn.inc'

c-----------------------------------------------------------
c     Command-line arguments
c-----------------------------------------------------------
      integer        xlevel,     olevel
      real*8         tfinal,     dtout,     tol

c-----------------------------------------------------------
c     Storage for coordinates of spatial mesh and approx.
c     solution.
c-----------------------------------------------------------
      real*8         xmin,          xmax
      parameter    ( xmin = 0.0d0,  xmax = 1.0d0 )
      integer        nxmax
      parameter    ( nxmax = 32 769 )
      real*8         x(nxmax),      y(2 * nxmax)

c-----------------------------------------------------------
c     LSODA declarations
c-----------------------------------------------------------
      external       fcn,        jac
      integer        neq
      real*8         t,          tout
      real*8         rtol,       atol
      integer        itol
      integer        itask,      istate,     iopt
      integer        lrw
      parameter    ( lrw = 22 + 2 * nxmax * 16 )
      real*8         rwork(lrw)
      integer        liw
      parameter    ( liw = 20 + 2 * nxmax )
      integer        iwork(liw)
      integer        jt

c-----------------------------------------------------------
c     Locals.
c-----------------------------------------------------------
      real*8         h
      integer        j,          nx,         stride  

c-----------------------------------------------------------
c     This function, defined in the p410f library, returns
c     its integer argument as a character string.
c-----------------------------------------------------------
      character*2    itoc

c-----------------------------------------------------------
c     Argument parsing and checking.
c-----------------------------------------------------------
      if( iargc() .ne. 5 ) go to 900
      xlevel = i4arg(1,-1)
      olevel = i4arg(2,-1)
      tfinal = r8arg(3,-1.0d0)
      dtout  = r8arg(4,-1.0d0)
      tol    = r8arg(5,-1.0d0)
      if( xlevel .lt. 1  .or.  olevel .lt. 1  .or.  
     &    olevel .gt. xlevel  .or. tfinal .lt. 0.0d0  .or.
     &    tol .lt. 0.0d0 ) go to 900

c-----------------------------------------------------------
c     Set up mesh, compute output stride, and initialize 
c     mesh coordinates and solution. 
c-----------------------------------------------------------
      nx = 2**xlevel + 1
      if( nx .gt. nxmax ) then
         write(0,*) 'wave: Requested nx = ', nx, 
     &              'exceeds maximum ', nxmax
         stop
      end if
      stride = 2**(xlevel - olevel)
      h   = (xmax - xmin) / (nx - 1)
      hm2 = 1.0d0 / (h * h)
      x(1)    = xmin
      y(1)    = 0.0d0
      y(1+nx) = 0.0d0
      do j = 2 , nx -1
         x(j)    = x(j-1) + h
         y(j)    = exp( -((x(j) - 0.5d0) / 0.1d0)**2 )
         y(j+nx) = 0.0d0
      end do
      x(nx)    = 1.0d0
      y(nx)    = 0.0d0
      y(nx+nx) = 0.0d0

c-----------------------------------------------------------
c     Set LSODA parameters 
c-----------------------------------------------------------
      neq    = 2 * nx 
      itol   = 1         ! Indicates that 'atol' is scalar
      rtol   = tol       ! Use same relative and absolute
      atol   = tol       ! tolerances.
      itask  = 1         ! Normal computation
      iopt   = 0         ! Indicates no optional inputs
      jt     = 2         ! Jacobian type 

c-----------------------------------------------------------
c     Output initial solution.
c-----------------------------------------------------------
      t = 0.0d0
      call xvs('u'//itoc(xlevel),t,x,y,nx)
      call gnuout(y,x,nx,t,stride)
c-----------------------------------------------------------
c     Integrate the approximate solution of the PDE 
c     using LSODA.
c-----------------------------------------------------------
      istate = 1
      do while( t .lt. tfinal ) 
         tout = t + dtout
c-----------------------------------------------------------
c        Call lsoda to integrate system on [t ... tout]
c-----------------------------------------------------------
         call lsoda(fcn,neq,y,t,tout,
     &              itol,rtol,atol,itask,
     &              istate,iopt,rwork,lrw,iwork,liw,jac,jt)
c-----------------------------------------------------------
c        Check return code and exit with error message if 
c        there was trouble.
c-----------------------------------------------------------
         if( istate .lt. 0 ) go to 950
c-----------------------------------------------------------
c        Output solution.
c-----------------------------------------------------------
         call xvs('u'//itoc(xlevel),t,x,y,nx)
         call gnuout(y,x,nx,t,stride)
      end do

      stop

 900  continue
         write(0,*) 'usage: wave <xlevel> <olevel> '//
     &              '<tfinal> <dtout> <tol>' 
      stop

 950  continue
         write(0,*) 'wave: Exiting due to LSODA failure'
         write(0,*) 'wave: Interval ', t, t + dtout
         write(0,*) 'wave: LSODA return code ', istate
      stop

      end

c===========================================================
c     Output to standard out for subsequent plotting via 
c     gnuplot.
c===========================================================
      subroutine gnuout(u,x,nx,t,stride)
         implicit       none

         integer        nx,           stride
         real*8         u(nx),        x(nx),        t

         integer        j
         
         do j = 1 , nx , stride
            write(*,*) t, x(j), u(j)
         end do
         write(*,*)

         return

      end
