!
!  This fortran 90 program is used in MATH 365 - Computational
!  Fluid Dynamics at James Madison University. This program was
!  written by James Sochacki of the Department of Mathematics
!  at James Madison University as part of the NSF Grant - A Collaborative 
!  Compuatational Science Program between JMU and North Carolina Central
!  University. His info is
!  
!  Jim Sochacki
!  Department of Mathematics
!  James Madison University
!  Harrisonburg, VA   22807
!  jim@math.jmu.edu


      program burgerc_lf

! This code calculates the solution to the nonlinear burger's equation in
! conservation form using Lax-Wendroff. The initial condition is a front based
! on arctan(x).

      implicit none

      real ::  a_left,a_right,q,r,s
      real ::  dx,dt,time
      real ::  lambda1,lambda2,L,uxb,uxf
      real,dimension(:),allocatable ::  u0,u1
      integer :: count,tsteps,jmax,i,j,interval
      character*11 :: filename
      character*6 :: name

      print*,' This code computes the solution for the burger equation '
      print*,' u_t + (0.5*u^2)_x = 0 ; u(x,0) = q*arctan(s x) + r  '
      print*,' using leap frog for -L < x < L. '
      print*,'  '
       print*,' Give a_left the value of w(x,0) at -infinity '
       read*,a_left
       print*,' Give a_right the value of w(x,0) at -infinity '
       read*,a_right
       print*,' Give s the "slope" of w(0,0) '
       read*,s
       print*,' Give the time to run the model '
       read*,time
       print*,' Give the time step '
       read*, dt
       print*,' Give L '
       read*,L
       print*,' Give the grid size '
       read*, dx
       print*,' Give the snapshot interval (A snapshot will be made '
       print*,' every nth time step. You input n.) '
       read*,interval
       print*,' Give the filename where these snapshots should be '
       print*,' saved (6 characters exactly). '
       read (5,' (a6)') name
       
       jmax = nint(2*L/dx)+1
       tsteps = nint(time/dt)+1

        allocate(u0(0:jmax+1))
        allocate(u1(0:jmax+1))

        r = (a_left+a_right)/2
        q = (a_right-a_left)/(4.*atan(1.0))

        do j=0,jmax+1

          u0(j) = q*atan(s*(j*dx-L))+r
          u1(j) = u0(j)

        end do

        i = 0
        write (filename,'(a6,i5.5)') name,i

        open (6,file=filename,status='new')
                
        do j=1,jmax

          write (6,'(f14.8)') u1(j)

        end do
            
        count=0
        close (6)
!
! Start the calculations.
!

         do i=1,tsteps

           count=count+1
           if (count==interval) print*,' step = ',i

           do j=1,jmax

              uxf = (u0(j+1)+u0(j))/2
              uxb = (u0(j)+u0(j-1))/2
              u1(j) = u0(j)-(u0(j+1)**2-u0(j-1)**2)*(dt/(4*dx))+ &
                      (uxf*(u0(j+1)**2-u0(j)**2)- &
                       uxb*(u0(j)**2-u0(j-1)**2))*(dt/(2*dx))**2
           end do
!
! Create the snapshots.

          if (count==interval) then

        write (filename,'(a6,i5.5)') name,i

        open (6,file=filename,status='new')
                
        do j=1,jmax

          write (6,'(f14.8)') u1(j)

        end do
            
        count=0
        close (6)

          end if
!
!  Update the calculations in time
!
          do j=1,jmax

            u0(j)=u1(j)

          end do

        end do

        print*,' the snapshots are in file: ',name

      end program burgerc_lf