!
!  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

!       This code looks at the leap frog finite difference solution 
!       to the advection diffusion equation 
!       w_t + c w_x + d w_(xx) = 0 ; -L < x < L   
!       w(x,0) = q*atan(s*x)+r 
!       q and r are determined from a_left and a_right the values of w(x,0)
!       at -infinity and infinity, respectively.

!  w0 is the array of space (x) grid values for w at the past time step
!  w1 is the array of space (x) grid values for w at the current time step
!  w2 is the array of space (x) grid values for w at the future time step
!  for the leap frog scheme. The first values for w1 are obtained using
!  the upwind finite difference scheme.

       program adv_diff_lf

       implicit none

       real, dimension(:), allocatable :: w0,w1,w2
       real :: c,d,dt,dx,L,time,q,r,s,a_left,a_right
       integer :: jmax,ja,jb,iadd,ntstps,i,j,n,count,interval
       real :: c1,c2
       character*11 :: filename
       character*6 :: name

       print*,' This code looks at the leap frog finite difference solution '
       print*,' to the advection diffusion equation '
       print*,' w_t + c w_x + d w_(xx) = 0 ; -L < x < L '  
       print*,' w(x,0) = q*atan(s*x)+r '
       
       print*,' Give c '
       read*, c
       print*,' Give d '
       read*, d
       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 is the number of x grid values. ntstps is the number of time steps
      
       jmax = nint(2*L/dx)+1
       ntstps = nint(time/dt)+1

! set up the memory for the grid space arrays.

       allocate(w0(jmax+1))
       allocate(w1(jmax+1))
       allocate(w2(jmax+1))

!   set up the difference parameters
       
           c1 = c*dt/dx
           c2 = -d*dt/(dx*dx)

!  set up the initial condition

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

       do j=1,jmax+1

         w0(j) = q*atan(s*(-L+(j-1)*dx))+r
         w1(j) = w0(j)
         w2(j) = w0(j)

       end do

!   set up the difference parameters
 
       c2 = -2.*d*dt/(dx*dx)
       c1 = c*dt/dx

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

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

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

        end do
            
        count=0
        close (6)


! The first time step is calculated using upwind

       if (c<0) then    

         do j = 1,jmax
 
           w1(j) = w0(j)+c1*(w0(j+1)-w0(j))+  &
                   c2/2.*(w0(j+1)-2.*w0(j)+w0(j-1))
         end do

       else    

         do j = 1,jmax
 
           w1(j) = w0(j)+c1*(w0(j)-w0(j-1))+  &
                   c2/2.*(w0(j+1)-2.*w0(j)+w0(j-1))
         end do

       end if

!  start the time iterations

       do n = 1,ntstps 

! use leap frog to approximate the PDE

         do j=2,jmax-1

           w2(j) = w0(j)+c1*(w1(j+1)-w1(j-1))+     &
                   c2*(w1(j+1)-2.*w1(j)+w1(j-1))

         end do

         count = count+1 

         if (count==interval) then
!  output a snapshot

            write (filename,'(a6,i5.5)') name,n
 
            open (6,file=filename,status='new')
            
            
            do j=1,jmax

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

            end do

            
           count=0
           close (6)


         end if
         
! update the time step

         do j=1,jmax
           w0(j) = w1(j)
           w1(j) = w2(j)
         end do

       end do

       close (6)

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

       end program adv_diff_lf