program advdif

      parameter( mx=400 )

      common uL, Re, pi

      dimension x(mx), dx(mx),
     &          wd1m(mx), wd10(mx), wd1p(mx),
     &          wd2m(mx), wd20(mx), wd2p(mx),
     &          u(mx), u_in(mx), u_ss(mx), u_exact(mx), a(mx), rhs(mx)

c   *** Set parameters.  The problem is solved in non-dimensional form
c   *** and these parameters are used only to scale the output.

      pi = acos(-1.0)
      xL = 1.0      ! domain length
      uL = 2.0      ! boundary value, u(xL,t)
      cw = 1.0      ! wave speed

      open(unit=1,file='frames.dat',form='unformatted')
      open(unit=2,file='final.dat')
      open(unit=3,file='error.dat')

      print *, 'enter Nx, StretchingFactor, Re, and i_fd_type'
      read(5,*) Nx, sigma, Re, i_fd_type

      if( Nx+1 .gt. mx ) then
         print *, ' '
         print *, 'Error, Nx exceeds the parameter mx'
         print *, 'Either choose a smaller Nx or enlarge',
     &           ' mx and recompile'
         print *, ' '
         stop
      end if

      if( Re .gt. 70 ) then
         print *, ' '
         print *, 'Warning, The exact solution suffers from round off',
     &           ' error for Re>70.'
         print *, 'The numerical solution is fine, but the spike in',
     &           ' the error history at small'
         print *, 'times is due to the inability to compute the exact',
     &           ' solution accurately.'
         print *, ' '
      end if

      vis = 1.0/Re
      R_delta = Re/float(Nx)

c   *** Compute the sine series coefficients for the exact solution.
c   *** An increasingly large number of modes are required at high Reynolds number.
c   *** The following curve fit expression attempts to do this.

      N_modes = nint( exp( Re/30.0 + 2.4 ) )

      if( N_modes+1 .gt. mx ) then
         print *, ' '
         print *, 'Error, ', N_modes, ' modes are required for the',
     &           ' exact solution and this exceeds the parameter mx'
         print *, 'Either choose a smaller Reynolds number or enlarge',
     &           ' mx and recompile'
         print *, ' '
         stop
      end if

      call sine_coef( x, u, u_in, u_ss, a, N_modes )

c   *** Compute the mesh

      dx_uniform = 1.0/float(Nx)

      if( sigma .eq. 0.0 ) then
         dx0 = dx_uniform
      else
         dx0 = sigma/((1.0+sigma)**float(Nx)-1.0)
      end if

      x( 1) = 0.0
      dx(1) = dx0
      dx_max = dx(1)
      dx_min = dx(1)
      do i=2, Nx+1
         x( i) = x(i-1) + dx(i-1)
         dx(i) = (1.0+sigma)*dx(i-1)
         dx_min = min( dx(i), dx_min )
         dx_max = max( dx(i), dx_max )
      end do

c   *** Determine the maximum allowable CFL number and the time step

      if( i_fd_type .eq. 0 ) then
         CFL = min( 2.0/R_delta, 0.5*R_delta )
      else
         CFL = R_delta/(2.0 + R_delta)
      end if

      print *, ' '
      print *, 'R_delta, CFL = ', R_delta, CFL
      print *, ' '
      print *, 
     &'If the simulation blows up, you may need to lower the time step.'
      print *, 'Enter a new CFL number, or a negative number to keep',
     &        ' the value shown above.'
      read(5,*) ans
      if( ans .gt. 0.0 ) CFL=ans

c   *** Determine the time step.  The factor dx_min/dx_max attempts to
c   *** account for the reduction in stability on a stretched mesh

      dt = CFL*dx_uniform*(dx_min/dx_max)

c   *** Set the weights for the first derivative

      if( i_fd_type .eq. 1 ) then   ! first order forward
         do i=2, Nx
            wd1m(i) =  0.0
            wd10(i) = -1.0/dx(i)
            wd1p(i) =  1.0/dx(i)
         end do
      else if( i_fd_type .eq.  0 ) then   ! second order central
         do i=2, Nx
            w1 = dx(i-1)/( (dx(i)+dx(i-1))*dx(i  ) )
            w2 = dx(i  )/( (dx(i)+dx(i-1))*dx(i-1) )
            wd1m(i) = -w2
            wd10(i) =  w2 - w1
            wd1p(i) =  w1
         end do
      else if( i_fd_type .eq. -1 ) then   ! first order backward
         do i=2, Nx
            wd1m(i) = -1.0/dx(i-1)
            wd10(i) =  1.0/dx(i-1)
            wd1p(i) =  0.0
         end do
      else
         print *, 'bad input for i_fd_type'
         stop
      end if

c   *** Set the weights for the second derivative

      do i=2, Nx
         w1 = 1.0/( 0.5*(dx(i)+dx(i-1))*dx(i  ) )
         w2 = 1.0/( 0.5*(dx(i)+dx(i-1))*dx(i-1) )
         wd2m(i) =   w2
         wd20(i) = -(w1 + w2)
         wd2p(i) =   w1
      end do

c   *** Determine the number of time steps and the skip factor for animation

      t_final  = 5.0
      n_frames = 100
      dt_print = t_final/float(n_frames)

      nt = nint(t_final/dt)
      i_skip = nint(max(dt_print/dt, 1.0))
      n_frames = nt/i_skip + 1

c   *** Initialize the solution and write it to the frames file

      t = 0.0

      call initial(      x, Nx, u    )
      call steady_state( x, Nx, u_ss )
      call exact(        x, a, u_ss, t, Nx, N_modes, u_exact )

      write(1) Nx+1, xL, sigma, cw, Re, CFL, i_fd_type, 
     &         nt, t_final, dt, dt_print, i_skip, n_frames

      n = 0
      i_frame = 1
      write(1) i_frame, n, t
      write(1) (u(      i),i=1,Nx+1)
      write(1) (u_exact(i),i=1,Nx+1)
      write(3,10) xL*t/cw, 0.0, 0.0

c   *** Advance the solution in time

      do n=1, nt

         do i=2, Nx
            adv =     ( wd1m(i)*u(i-1) + wd10(i)*u(i) + wd1p(i)*u(i+1) )
            dif = vis*( wd2m(i)*u(i-1) + wd20(i)*u(i) + wd2p(i)*u(i+1) )
            rhs(i) = dif - adv
         end do

         do i=2, Nx
            u(i) = u(i) + dt*rhs(i)
         end do

         t = t + dt

         if( mod(n,i_skip) .eq. 0 .or. n .eq. nt ) then

            call exact( x, a, u_ss, t, Nx, N_modes, u_exact )

            sum1 = 0.0
            sum2 = 0.0
            do i=2, Nx
               sum1 = sum1 + ( u(i) - u_exact(i) )
               sum2 = sum2 + ( u(i) - u_exact(i) )**2
            end do
            avg = sum1/float(Nx-1)
            rms = sqrt( sum2/float(Nx-1) )

            write(3,10) xL*t/cw, avg, rms
10          format(1p,4e12.4)

            i_frame = i_frame + 1
            write(1) i_frame, n, t
            write(1) (u(      i),i=1,Nx+1)
            write(1) (u_exact(i),i=1,Nx+1)

         end if

      end do

      do i=1, Nx+1
         write(2,10) xL*x(i), u(i), u_exact(i)
      end do

      stop
      end


      subroutine initial( x, Nx, u )

c   *** Initial condition

      common uL, Re, pi

      real x(Nx+1), u(Nx+1)

      do i=1, Nx+1
         u(i) = uL*x(i)
      end do

      return
      end


      subroutine steady_state( x, Nx, u_ss )

c   *** Exact steady-state solution 

      common uL, Re, pi

      real x(Nx+1), u_ss(Nx+1)

      do i=1, Nx+1
         u_ss(i) = uL*( exp(Re*x(i)) - 1.0 )/( exp(Re) - 1.0 )
      end do

      return
      end


      subroutine sine_coef( x, u, u_in, u_ss, a, N_modes )

c   *** Compute the sine coefficients needed for the exact solution

      common uL, Re, pi

      real x(N_modes+1), u(N_modes+1), u_in(N_modes+1), u_ss(N_modes+1),
     &     a(N_modes)

c   *** Sample the function to be transformed on a uniform mesh

      dx0 = 1.0/float(N_modes)
      do i=1, N_modes+1
         x(i) = float(i-1)*dx0
      end do

      call initial(      x, N_modes, u_in )
      call steady_state( x, N_modes, u_ss )

      do i=1, N_modes+1
         u(i) = exp( -0.5*Re*x(i) )*( u_in(i) - u_ss(i) )
      end do

c   *** Compute the sine transform

      fac = 1.0/float(N_modes)
      do j=1, N_modes-1
         sum1 = 0.0
         do i=2, N_modes
            sum1 = sum1 + u(i)*sin( float(j)*pi*x(i) )
         end do
         a(j) = 2.0*sum1*fac
      end do
      a(N_modes  ) = 0.0
      a(N_modes+1) = 0.0

c      do j=1, N_modes+1
c         write(10,10) x(j), u(j), a(j)
c10       format(1p,3e12.4)
c      end do

      return
      end


      subroutine exact( x, a, u_ss, t, Nx, N_modes, u_exact )

c   *** Exact solution for any time

      common uL, Re, pi

      real x(Nx+1), a(N_modes), u_ss(Nx+1), u_exact(Nx+1)

      do i=1, Nx+1
         sum1 = 0.0
         do j=1, N_modes
            wave_no = float(j)*pi
            sum1 = sum1 + exp(-wave_no**2*t/Re)*a(j)*sin(wave_no*x(i))
         end do
         u_exact(i) = exp( 0.5*Re*(x(i)-0.5*t) )*sum1 + u_ss(i)
      end do

      return
      end