program LAX
C
C example of Lax algorithm (numerical viscosity)
C for sollution  of 
C PDE  du/dx = - du/dt
C cf Numerical Recipes in Fortran, p. 826-7
C
C  CWJ SDSU April 2005
C
      implicit none

C......... parameters of box..............

      real L  !  box is from 0 to L
      real dx !  discretiztion of box
      integer Nx !  number of points in box

C......... parameters of time.............

      real dt !  time step
      integer Nt !  number of time steps to take
      real t
      integer n

C.....parameters of wave..................

      integer size
      parameter(size=1000)
      real u1(0:size),u2(0:size)

C.......... parameters of initial gaussian

      real xc,xw    ! center and width
      integer j

C............READ IN PARAMETERS OF BOX...........

      print*,' Enter size of box L '
      read*,L
      print*,' Enter # of points '
      read*,Nx
      dx = L/float(nx)      

C...........read in initial parameters of gaussian wave 

      print*,' Enter center, width of initial gaussian '
      read*,xc,xw

C...........fill in gaussian.........................

      do j = 1,Nx-1
         u1(j) = exp(- (dx*j - xc)**2/(2*xw**2))
      enddo
      u1(0) = 0.
      u1(Nx) = 0.0

      open(unit=1,file='initial.dat',status='unknown')
      do j = 0,Nx
        write(1,*)j*dx,u1(j)
      enddo
      close(unit=1)

C...........read in parameters of time step.......

      print*,' Enter time step '
      print*,' (suggest less than ',dx,') '
      read*,dt       

      print*,' Enter time to travel  '
      read*,t
      Nt = int(t/dt)

      do n = 1,nt/2
         call Lax_step(size,Nx,dx,dt,u1,u2)
         call Lax_step(size,Nx,dx,dt,u2,u1)
      enddo
      
      open(unit=2,file='final.dat',status = 'unknown')
      do j = 0,Nx
        write(2,*)j*dx,u1(j)
      enddo
      close(unit=2)      
      open(unit=3,file='compare.dat',status = 'unknown')
      xc = xc-(nt/2)*2*dt
      do j = 0,Nx
        write(3,*)j*dx,exp(- (dx*j - xc)**2/(2*xw**2))
      enddo
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine Lax_step(size,Nx,dx,dt,uold,unew)
C
C  does a single Lax step
C
      implicit none

      integer size
      real uold(0:size),unew(0:size)
   
      integer Nx
      real dx
      real dt

      integer j

      do j = 1,Nx-1
        unew(j) = 0.5*(uold(j+1)+uold(j-1))+
     &   dt*(uold(j+1)-uold(j-1))/(2*dx) 
      enddo      
      unew(0) = 0. ! uold(0)+dt*uold(1)/(2*dx)
      unew(Nx) =0. !  uold(nx)+dt*uold(nx-1)/(2*dx)

      return
      end