c---------------------------------------------------
c  Fortran Problem 3
c---------------------------------------------------

      program          wave1d
      
      implicit         none

c=============================================================
c   Variables   ::
c-------------------------------------------------------------
c
c   h_x ,  h_t  :  specify the value of the mesh spacing
c
c   maxj        :  maximum # of spatial points to consider
c
c   n,  nswap    :  n = number of spatial points actually using,  
c   np1, nm1         nswap, np1 and nm1 are dummies used in update loop
c
c   u           :  u is solution array
c
c   dt_dx       :  Courant number
c
c   left_amp    :  Initial Amplitudes of Gaussian wavepackets
c   rite_amp
c
c   level       :  level of grid, olevel is output level
c   olevel               
c
c   ncross      :  crossing time for wave ( equals one in this case)
c
c   x           :  spatial mesh storage
c   xmin,xmax   :  min and max x; 0 < x < 1 in this case
c
c   nx, nt      :   the number of x and t points
c
c-------------------------------------------------------------
c==============================================================

c----------------------------------
c   Parsed variables
c----------------------------------
      integer    level,     i4arg,   i4_never, 
     &           ncross,    olevel,  iargc
      real*8     dt_dx,     r8arg,   r8_never,
     &           left_amp,  rite_amp
c----------------------------------
c   Initialize Gaussian variables
c----------------------------------
      integer    maxj,      nx,      i
      parameter (maxj = 120 000)

      real*8     xmin,      xmax,    h_x

      real*8     xO,        del 
      parameter (xmin = 0.0d0, xmax = 1.0d0)

      real*8     x(maxj),   l(maxj),  dl(maxj),
     &           ddl(maxj), r(maxj),  dr(maxj), ddr(maxj)
c---------------------------------
c   Third Order initial data
c---------------------------------
      real*8     u(maxj,2), h_t
c---------------------------------
c   Time-step Loop variables
c---------------------------------
      integer    n,   nm1,    np1,   j,  t,
     &           nt,  ofreq,  nswap
      real*8     time
      integer    out

c---------------------------------
c   End of Variable Section
c---------------------------------

c==================================
c   Parse command line for arguments
c==================================
      if (iargc() .ne. 6) then
      write(0,*)'usage: wave1d <level>  <dt/dx>  <ncross> ',
     &' <lt_amp>  <rt_amp>  <olevel>'
      stop
      end if

      level = i4arg (1, i4_never) 
      if (level .eq. i4_never) then
        write(0,*) ' Error:1'
      stop
      end if

      dt_dx    = r8arg (2, r8_never)
      if (dt_dx .eq. r8_never) then
        write(0,*) ' Error:2'
      stop
      end if

      ncross   = i4arg (3, i4_never)
      if (ncross .eq. i4_never) then
        write(0,*) ' Error:3'
      stop
      end if

      left_amp = r8arg (4, r8_never)
      if (left_amp .eq. r8_never) then
        write(0,*) ' Error:4'
      stop
      end if

      rite_amp = r8arg (5, r8_never)
      if (rite_amp .eq. r8_never) then
        write(0,*) ' Error:5'
      stop
      end if

      olevel   = i4arg (6, i4_never)      
      if (olevel  .eq. i4_never) then
        write(0,*) ' Error:6'
      stop
      end if


c===================================
c   Initialize data -- In this case use 2 gaussians
c===================================
        
c---------------------------------
c   Set up mesh
c--------------------------------- 
      h_x =  (xmax - xmin) / (2**level)
      nx   =  2**level +1
      do i = 1, nx
         x(i) = xmin  +  (i-1) * h_x
      end do

c---------------------------------
c   Some initial args for subroutine
c---------------------------------
      xO  = 0.5d0
      del = 0.1d0

      call dvgaussian (l, dl, ddl, x, nx, left_amp, xO, del)
      call dvgaussian (r, dr, ddr, x, nx, rite_amp, xO, del)
c---------------------------------
c   Set up initial data - 3rd order - t=0 and t=1
c---------------------------------
      h_t = (dt_dx * h_x)
      write(0,*) ' through dvvgaussian'    
      do i = 1, nx
          u(i, 1) = l(i) + r(i) 
          u(i, 2) = l(i) + r(i) + (h_t) * (dl(i)-dr(i)) + 
     &              0.5d0 * (h_t**2) * (ddl(i) + ddr(i))
      end do 

c=======================================
c   Begin Time Step Loop.  NOTE:  This loop was written by Dr. Matthew Choptuik
c=======================================
      n   = 2
      nm1 = 1
      np1 = 1
      ofreq = 2**(level - olevel)

      call vsxynt('wave', 0.0d0, x, u(1,np1), nx)
      call gnuout ( u(1,np1), x, nx, 0.0d0, ofreq )

      if (level .eq. olevel) then
      write(0,*) 'Olevel is equal to level'
          call vsxynt('wave',h_t, x, u(1,np1), nx)
          call gnuout ( u(1,np1), x, nx, h_t, ofreq )
      end if 

c-------------------------------------
c   I'm not sure if this is correct def. for nt.  But i think so.
c-------------------------------------
      nt  = ncross*nx
      write(0,*) 'nt = ', nt
      write(0,*) ' at time step loop'
      do t = 2 , nt
           u(1, np1) = 0.0d0
           do j = 2 , nx - 1
             u(j, np1) = 2.0d0 * u(j, n) - u(j, nm1) + 
     &       (dt_dx**2) * (u(j+1,n) - 2.0d0 * u(j,n) + u(j-1,n))
           end do
      u(nx,np1) = 0.0d0
c=======================================
c    Periodic Output of data
c=======================================
      time  = h_t * t
      write(0,*) 'time = ', time
c---------------------
c  Define time and out
c--------------------
      out = mod(t,ofreq)
 
      if (out .eq. 0) then
          call vsxynt('wave', time, x, u(1,np1), nx)
          call gnuout ( u(1,np1), x, nx, time, ofreq )
      end if

      nswap =  np1
      np1   =  n 
      n     =  nswap
      nm1   =  np1

      end do
   
      end