!
!  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 upwind 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 
!      The values for 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 future time step
!  for the upwind finite difference scheme.
      
       program adv_diff_upwind

       implicit none

       real, dimension(:), allocatable :: 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 upwind 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(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

         w1(j) = q*atan(s*(-L+(j-1)*dx))+r
         w2(j) = w1(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)') w2(j)

        end do
            
        count=0
        close (6)


!  start the time iterations

       do n = 1,ntstps
         
! upwind on c to approximate the solution to the PDE

         if (c>0) then

           do j=2,jmax-1

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

           end do

         else

           do j=2,jmax-1

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

           end do
               
         end if

         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
           w1(j) = w2(j)
         end do

       end do

       close (6)

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

       end program adv_diff_upwind