program   transport
	include 'transin'

c  1d simulation program for linear transport using finite differences 
c  1 - the upwind scheme
c  2 - the leapfrog scheme
c  3 - the Lax-Wendroff scheme
c  4 - Crank-Nicholson scheme
c  output of program 
c    sim1.bin - unformated (first nx,x, then each output: time, f)     
c    sim1.asci - formated (first all parameters, then: time, f) 
      open (14,file='sim1.bin',form='unformatted')
      open (16,file='sim1.asci')
c   f,fold,foold are normalized fields for integration   
c   read input parameters
      call pread
c   determine grid:
      call grid
c   calculate runtime params:
      call params
      write(14) nx
      write(14) x
c   write some parameters:      
      call pwrite
      write(14) c
      write(14) alpha
c   set initial conditions:      
       if (init .eq. 1) call initcon1
       if (init .eq. 2) call initcon2
       if (init .eq. 3) call initcon3
c   integration: 
c    intftcs - 1d heat conduction with ftcs 
c              (forward time centered space)
       nt = nt+1
       time = time + dt
       if (method .eq. 1) call upwind
       if (method .eq. 2) call upleap
       if (method .eq. 3) call laxwend
       if (method .eq. 4 .and. que .eq. 0.) call crankni
       if (method .eq. 4 .and. que .ne. 0.) call cnupwind
       call bdcon
       if (mod(nt,nout) .eq. 0) call out
       do while( (nt .le. ntmax) .and. (time .le. tmax) )
         nt = nt+1
         time = time + dt
         if (method .eq. 1) call upwind
         if (method .eq. 2) call leap
         if (method .eq. 3) call laxwend
         if (method .eq. 4 .and. que .eq. 0.) call crankni
         if (method .eq. 4 .and. que .ne. 0.) call cnupwind
         call bdcon
c   write output f:      
        if (mod(nt,nout) .eq. 0) call out
       end do
      call out
c   determine and write exact solution to order nterm:  
      if (init .eq. 1) call analyt1
      if (init .eq. 2) call analyt2
      if (init .eq. 3) call analyt3
      close(14) 
      close(16) 
      stop
      end
c ------------ : ------------
      subroutine pread
      include 'transin'
c   input data:
c   parameters explained in input file sim1.dat
       open(11,file='trans1.dat')          
       read(11,1) method
       read(11,1) ntmax
       read(11,3) tmax
       read(11,3) xmin,xmax
       read(11,2) alpha
       read(11,3) u,c,dt
       read(11,1) cordt
       read(11,3) delta,que
       read(11,1) nout
       read(11,3) f0
       read(11,1) init
       read(11,3) lambda
       read(11,1) parbound
       read(11,3) fb(1),fb(2)
       read(11,1) nterm
       close(11) 
       
       if (method.eq.4 .and. que.ne.0. .and. parbound.eq.1) then
         parbound=2
         write(6,*) '!!!! NO PERIODIC BOUNDARY CONDITION IMPLEMENTED',
     +            '      FOR THE CRANK-NICHOLSON SCHEME!!!! '
         write(6,*) '!!!! BOUNDARY CONDITION SET TO DIRICHLET!!!! '
       endif
    1  format(i5)
    2  format(e10.2)
    3  format(f8.2)
      return
      end
c ------------ : ------------
      subroutine grid
      include 'transin'
       integer i
       xsize=xmax-xmin
       dx = (xmax-xmin)/(nx-3.)
       do 10 i = 1,nx
   10    x(i)  = xmin+dx*(i-2.)
      return
      end
c ------------ : ------------
      subroutine params
      include 'transin'
c   compute parameters for schemes
       if (cordt .eq. 0) dt = c*dx/u
       if (cordt .eq. 1) c = u*dt/dx
       s = alpha*dt/dx/dx
       rcell=0.
       if (alpha .ne. 0.) rcell = u*dx/alpha
       nt = 0
       time = 0.0      
       
c       delta=1./6.+c**2/12
c        que=0.5+0.25*c*c
       
       star = 0.
       if (method .eq. 1)  then
         s1 = s+c
         s2 = 1. - 2.*s - c
       endif
       if (method .eq. 2)  then
         star = 1.+2.*s
         s1   = (c+2.*s)/star
         s2   = (1.-2.*s)/star
         s3   = (-c+2.*s)/star
       endif
       if (method .eq. 3)  then
         star = s + 0.5*c*c
         s0   = -que*c/3.
         s1   = 0.5*c + star + que*c
         s2   = 1. - 2.*star - que*c
         s3   = -0.5*c + star + que*c/3.
       endif
       if (method .eq. 4)  then
         a0 = que*c/6.
         a1 = delta - 0.25*c - 0.5*s - 0.5*que*c
         a2 = 1. - 2.*delta + s + 0.5*que*c
         a3 = delta + 0.25*c - 0.5*s - 0.5*que*c/3.
         s0 = -que*c/6.
         s1 = delta + 0.25*c + 0.5*s + 0.5*que*c
         s2 = 1. - 2.*delta - s - 0.5*que*c
         s3 = delta - 0.25*c + 0.5*s + 0.5*que*c/3.
       endif
      write(6,*) delta
      return
      end
c ------------ : ------------
      subroutine pwrite
      include 'transin'
       if (method .eq. 1) write(16,*) 'Upwind scheme (explicit)'
       if (method .eq. 2) write(16,*) 'Leapfrog scheme (explicit)'
       if (method .eq. 3) write(16,*) 'Lax-Wendroff scheme (explicit)'
       if (method .eq. 4) then 
         write(16,*) 'Crank-Nicholson scheme (implicit)'
       endif
       write(16,2) nx,ntmax,tmax, xmin, xmax
       write(16,3) dt,dx,nout
       write(16,4) u,alpha,c,s,rcell
       if (method .eq. 3 .or. method .eq. 4) write(16,5) delta, que
       write(16,6) init,lambda,f0
       if (parbound .eq. 1) write(16,*) 'Periodic boundary conditions'
       if (parbound .eq. 2) write(16,*) 
     +              'dirichlet boundary conditions with:'
       if (parbound .eq. 3) write(16,*) 
     +              'neumann boundary conditions with derivatives:'
       if (parbound .ne. 1) write(16,7) fb(1),fb(2)
       write(16,8) nterm
       write(16,*) 'Stability:'
       if (method .eq. 1) write(16,9) c+2.*s, rcell, 2./(1.-c)
       if (method .eq. 2) write(16,10) c, c**2
       if (method .eq. 3) write(16,11) c**2, 2*star, rcell
       if (method .eq. 4) write(16,12) rcell
       
       
       if (method .eq. 1) write(6,*) 'Upwind scheme (explicit)'
       if (method .eq. 2) write(6,*) 'Leapfrog scheme (explicit)'
       if (method .eq. 3) write(6,*) 'Lax-Wendroff scheme (explicit)'
       if (method .eq. 4) then 
         write(6,*) 'Crank-Nicholson scheme (implicit)'
       endif
       write(6,2) nx,ntmax,tmax, xmin, xmax
       write(6,3) dt,dx,nout
       write(6,4) u,alpha,c,s,rcell
       if (method .eq. 3 .or. method .eq. 4) write(6,5) delta, que
       write(6,6) init,lambda,f0
       if (parbound .eq. 1) write(6,*) ' Periodic boundary conditions'
       if (parbound .eq. 2) write(6,*) 
     +              'dirichlet boundary conditions with:'
       if (parbound .eq. 3) write(6,*) 
     +              'neumann boundary conditions with derivatives:'
       if (parbound .ne. 1) write(6,7) fb(1),fb(2)
       write(6,8) nterm
       if (method .eq. 1) write(6,9) c+2.*s, rcell, 2./(1.-c)
       if (method .eq. 2) write(6,10) c, c**2
       if (method .eq. 3) write(6,11) c**2, 2*star, rcell
       if (method .eq. 4) write(6,12) rcell
    2  format(1x,' nx =',i4,'   ntmax =',i5,'     tmax =',f8.2,
     +       '  xmin =',f5.1,'  xmax =',f5.1)
    3  format(1x,' dt =',f8.4,'    dx =',f6.4,' nout =',i5)
    4  format(1x,' u  =',f6.2,'  alpha =',e10.3,
     +           ' c =',f6.4,' s =',f6.4,'  rcell =',e10.3)
    5  format(1x,' delta=',f8.5,'  q(upwind) =',f8.5)
    6  format(1x,' init(cond) =',i3,'  lambda =',f8.4,
     +           ' f0(norm) =',f6.1)
    7  format(1x,' f_bound_min =',f7.3,' f_bound_max =',f7.3)
    8  format(1x,' Terms in exact solution (nterm) =',i4)

    9  format(1x,' Stability: c+2.*s =',f7.3,' <=1; ',/,
     +           '  Accuracy: rcell =',e10.3,',<= 2./(1.-c) =',e10.3)
   10  format(1x,' Stability: c =',f7.3,',<=1; ',
     +           '  Accuracy: c**2 =',e10.3,',<= 1')
   11  format(1x,' Stability: c**2 =',f7.3,',<= ',f8.3,' =2*s_star <=1',
     +         /,'  Oscillations: rcell =',e10.3,',<= 2')
   12  format(1x,'  Oscillations: rcell =',e10.3,',<= 2')
      return
      end
c ------------ : ------------
      subroutine initcon1
      include 'transin'
       integer i
       real    lhalf
       lhalf=0.5*lambda
       do 10 i = 1,nx-1
   10     f(i) = 0.
       do 20 i = 2,nx-1
   20    if (x(i) .ge. -lhalf .and. x(i).le. lhalf)
     +            f(i) = cos(pi*x(i)/lambda)
       do 25 i=1,nterm
         kf(i)= (pi/2.+(i-1.)*pi)/xmax
         arg1=pi/lambda+kf(i)
         arg2=pi/lambda-kf(i)
         if (arg2 .ne. 0) cf(i)=1./xmax*( 
     +            (sin(arg1*lhalf)/arg1)+(sin(arg2*lhalf)/arg2) )
   25 continue
      do 30 i = 2,nx-1
        fexact(i) = 0
        do 30 l = 1,nterm
        fexact(i)=fexact(i)+cf(l)*cos(kf(l)*x(i))
   30 continue
       call bdcon
       call out
      return
      end
c ------------ : ------------
      subroutine initcon2
      include 'transin'
       integer  i,nwave
       nwave=xsize/(2*lambda)
       kwave=2.*pi*nwave/xsize
       do 10 i = 2,nx-1
         f(i) = sin(kwave*x(i))
         fexact(i) = f(i)
   10  continue
       call bdcon
       call out
      return
      end
c ------------ : ------------
      subroutine initcon3
      include 'transin'
       integer i
       real    lhalf
       lhalf=0.5*lambda
       do 10 i = 1,nx-1
   10     f(i) = 0.
       do 20 i = 2,nx-1
   20    if (x(i) .ge. -lhalf .and. x(i).le. lhalf)
     +            f(i) = 1.
       do 25 i=1,nterm
         kf(i)= (pi/2.+(i-1.)*pi)/xmax
         cf(i)=2./xmax*sin( kf(i)*lhalf )/kf(i)
   25 continue
      do 30 i = 2,nx-1
        fexact(i) = 0
        do 30 l = 1,nterm
        fexact(i)=fexact(i)+cf(l)*cos(kf(l)*x(i))
   30 continue
       call bdcon
       call bdcon
       call out
      return
      end
c ------------ : ------------
      subroutine bdcon
      include 'transin'
       if (parbound .eq. 1) then
        f(1)  = f(nx-2)
        f(nx) = f(3)
       endif
       if (parbound .eq. 2) then
        f(2)  = fb(1)
        f(nx-1) = fb(2)
        f(1)  = fb(1)
        f(nx) = fb(2)
       endif
       if (parbound .eq. 3) then
        f(1)  = f(3)-2.*dx*fb(1)
        f(nx) = f(nx-2)+2.*dx*fb(2)
       endif
      return
      end
c ------------ : ------------
      subroutine upwind
      include 'transin'
      integer  i
       do 10 i = 1,nx
   10   fold(i) = f(i)
       do 20 i = 2,nx-1
   20   f(i) = s1*fold(i-1) + s2*fold(i) + s*fold(i+1)
      return
      end
c ------------ : ------------
      subroutine upleap
      include 'transin'
      integer  i
       do 10 i = 1,nx
   10   fold(i) = f(i)
       do 20 i = 2,nx-1
   20   f(i) = (s+c)*fold(i-1) + (1.-2.*s-c)*fold(i) + s*fold(i+1) 
      return
      end
c ------------ : ------------
      subroutine leap
      include 'transin'
      integer  i
       do 10 i = 1,nx
   10   foold(i) = fold(i)
       do 15 i = 1,nx
   15   fold(i) = f(i)
       do 20 i = 2,nx-1
   20   f(i) = s1*fold(i-1) + s2*foold(i) + s3*fold(i+1)
      return
      end
c ------------ : ------------
      subroutine laxwend
      include 'transin'
      integer  i
       do 10 i = 1,nx
   10   fold(i) = f(i)
       if (que .eq. 0) then
        do 30 i = 2,nx-1
   30    f(i) = s1*fold(i-1) + s2*fold(i) + s3*fold(i+1)
       else
        do 35 i = 3,nx-1
   35    f(i) = s0*fold(i-2) + s1*fold(i-1) + s2*fold(i) + s3*fold(i+1)
        if (parbound .eq. 1) 
     +   f(2) = s0*fold(nx-3) + s1*fold(1) + s2*fold(2) + s3*fold(3)
        if (parbound .eq. 2) 
     +   f(2) = s0*fb(1) + s1*fold(1) + s2*fold(2) + s3*fold(3)
        if (parbound .eq. 3) 
     +   f(2) = s0*( f(2)-2.*dx*fb(1) ) 
     +                        + s1*fold(1) + s2*fold(2) + s3*fold(3)
       endif
      return
      end
c ------------ : ------------
      subroutine crankni
      include 'transin'
      integer  j
      real     a(5,nx),r(nx),td(nx)
       do 10 j = 2,nx-1
        jm=j-1
        a(1,jm)=0.
        a(2,jm)=a1
        a(3,jm)=a2
        a(4,jm)=a3
        r(jm) = s1*f(jm)+s2*f(j)+s3*f(j+1)
   10  continue
       r(1)=r(1)-a(2,1)*f(1)
       a(2,1)=0.
       a(4,nx-2)=0.
       
       call banfac(a,nx-2,1)
       call bansol(r,td,a,nx-2,1)

       do 30 i = 2,nx-1
   30   f(i) = td(i-1)
      return
      end
c ------------ : ------------
      subroutine cnupwind
      include 'transin'
      integer  j
      real     a(5,nx),r(nx),td(nx)
       do 10 j = 2,nx-1
        jm=j-1
        a(1,jm)=a0
        a(2,jm)=a1
        a(3,jm)=a2
        a(4,jm)=a3
        r(jm) = s1*f(jm)+s2*f(j)+s3*f(j+1)
   10  continue
       r(1)=r(1)-a(2,1)*f(1)
       if (que .ne. 0.) then
        do 20 j = 2,nx-2
   20     r(j) = r(j)+s0*f(j-1)
        r(1)=r(1)-a(1,1)*f(1)
        r(2)=r(2)-a(1,2)*f(1)
c  reduce to tridiagonal:   
        do 25 j = 3,nx-2
         jm = j-1
         dum=a(1,j)/a(2,jm)
         a(2,j) = a(2,j) - a(3,jm)*dum
         a(3,j) = a(3,j) - a(4,jm)*dum
         a(1,j) = 0.
         r(j) = r(j) - r(jm)*dum
   25   continue
        a(1,1) = 0.
        a(1,2) = 0.
       endif
       a(2,1)=0.
       a(4,nx-2)=0.
       
       call banfac(a,nx-2,1)
       call bansol(r,td,a,nx-2,1)

       do 30 i = 2,nx-1
   30   f(i) = td(i-1)
      return
      end
c ------------ : ------------
      subroutine out
      include 'transin'
      integer  i
      real     fout(nx),fff(nx)
       do 10 i = 1,nx
   10   fout(i) = f0*f(i)
       write(16,1) time, (fout(i),i=1,nx)
       if (init .eq. 1) call analyt1
       if (init .eq. 2) call analyt2
       if (init .eq. 3) call analyt3
       do 20 i = 1,nx
   20   fff(i) = f0*fexact(i)
       write(14) time,rms,fout,fff
       write(16,*) 'time step:', nt
    1  format(' time = ',f5.0,'  T =', 11f7.2)
      return
      end
c ------------ : ------------
      subroutine analyt1
      include 'transin'
      integer  i
      real     sum
      alph1=0.0
       do 10 i = 1,nx
   10   xr(i) = x(i)-u*time
       ishift=0
       do while (xr(nx-2)+ishift*xsize .lt. xmin) 
         ishift=ishift+1
       end do
       do 15 i = 2,nx-1
         xr(i)=xr(i)+ishift*xsize
         if (xr(i) .lt. xmin) xr(i)=xr(i)+xsize
   15  continue
       do 30 i = 2,nx-1
         fexact(i) = 0
         do 30 l = 1,nterm
         fexact(i)=fexact(i)+cf(l)*exp(-kf(l)*kf(l)*alph1*time)
     +                      *cos(kf(l)*xr(i))
  30   continue
       sum=0.
       do 40 i = 2,nx-1
  40    sum = sum + ( fexact(i)-f0*f(i) )**2
       if (method .eq. 2) then 
         do 50 i = 2-mod(nt,2),nx,2
   50        sum = sum + ( fexact(i)-f0*f(i) )**2
       endif
c   rms is the rms error
       rms = sqrt(sum/nx)
       if (method .eq. 5)  rms = sqrt(2.*sum/(nx+1))
       write(16,1) nterm
       write(16,2) time, (fexact(i),i=1,nx)
       write(16,3) rms
c       write(6,3) rms
    1  format(1x,'No of terms in exact solution (nterm) =',i5)
    2  format(/,' time = ',f5.0,'  T =', 11f7.2)
    3  format(1x,' rms dif = ', e11.4,//)
      return
      end
c ------------ : ------------
      subroutine analyt2
      include 'transin'
      integer  i
      real     sum
      alph1=0.0
       do 10 i = 1,nx
   10   xr(i) = x(i)-u*time
       ishift=0
       do while (xr(nx-2)+ishift*xsize .lt. xmin) 
         ishift=ishift+1
       end do
       do 15 i = 2,nx-1
         xr(i)=xr(i)+ishift*xsize
         if (xr(i) .lt. xmin) xr(i)=xr(i)+xsize
         fexact(i) = exp(-kwave*kwave*alph1*time)*sin(kwave*xr(i))
   15  continue
       sum=0.
       do 40 i = 2,nx-1
  40    sum = sum + ( fexact(i)-f0*f(i) )**2
       if (method .eq. 2) then 
         do 50 i = 2-mod(nt,2),nx,2
   50        sum = sum + ( fexact(i)-f0*f(i) )**2
       endif
c   rms is the rms error
       rms = sqrt(sum/nx)
       if (method .eq. 5)  rms = sqrt(2.*sum/(nx+1))
       write(16,1) nterm
       write(16,2) time, (fexact(i),i=1,nx)
       write(16,3) rms
c       write(6,3) rms
    1  format(1x,'No of terms in exact solution (nterm) =',i5)
    2  format(/,' time = ',f5.0,'  T =', 11f7.2)
    3  format(1x,' rms dif = ', e11.4,//)
      return
      end
c ------------ : ------------
      subroutine analyt3
      include 'transin'
      integer  i
      real     sum
      alph1=0.0
       do 10 i = 1,nx
   10   xr(i) = x(i)-u*time
       ishift=0
       do while (xr(nx-2)+ishift*xsize .lt. xmin) 
         ishift=ishift+1
       end do
       do 15 i = 2,nx-1
         xr(i)=xr(i)+ishift*xsize
         if (xr(i) .lt. xmin) xr(i)=xr(i)+xsize
   15  continue
       do 30 i = 2,nx-1
         fexact(i) = 0
         do 30 l = 1,nterm
         fexact(i)=fexact(i)+cf(l)*exp(-kf(l)*kf(l)*alph1*time)
     +                      *cos(kf(l)*xr(i))
  30   continue
       sum=0.
       do 40 i = 2,nx-1
  40    sum = sum + ( fexact(i)-f0*f(i) )**2
       if (method .eq. 2) then 
         do 50 i = 2-mod(nt,2),nx,2
   50        sum = sum + ( fexact(i)-f0*f(i) )**2
       endif
c   rms is the rms error
       rms = sqrt(sum/nx)
       if (method .eq. 5)  rms = sqrt(2.*sum/(nx+1))
       write(16,1) nterm
       write(16,2) time, (fexact(i),i=1,nx)
       write(16,3) rms
c       write(6,3) rms
    1  format(1x,'No of terms in exact solution (nterm) =',i5)
    2  format(/,' time = ',f5.0,'  T =', 11f7.2)
    3  format(1x,' rms dif = ', e11.4,//)
      return
      end
c ------------ : ------------
      subroutine banfac(b,n,int)
      include 'transin'
c
c  factorizes band matrix arising from linear or quadratic elements
c  into l.u
c  
      real     b(5,nx)
      
c
c  int = 1, linear elements = tridiagonal system
c  int = 2, qudratic elements = pentadiagonal system
c  assumes that first equation formed at midside node
c  
      if (int .eq. 1) then
       np = n-1
       do 1 j = 1,np
        jp = j+1
        b(2,jp) = b(2,jp)/b(3,j)
        b(3,jp) = b(3,jp) - b(2,jp)*b(4,j)
    1  continue
       return
      else
       nh = n/2
       do 5 i = 1,2
        js = 3-i
        do 4 j = js,nh
         ja = 2*(j-1)
c  i = 1, first pass, reduction to tridiagonal
c  i = 2, second pass, reduction to upper triangular
         if (i .eq. 1) then
          jb = ja+2
          b(1,jb) = b(1,jb)/b(2,jb-1)
          b(2,jb) = b(2,jb) - b(1,jb)*b(3,jb-1)
          b(3,jb) = b(3,jb) - b(1,jb)*b(4,jb-1)
         else
          jb = ja+3
          b(2,jb-1) = b(2,jb-1)/b(3,jb-2)
          b(3,jb-1) = b(3,jb-1) - b(2,jb-1)*b(4,jb-2)
          if (jb .le. n) then
           b(2,jb) = b(2,jb)/b(3,jb-1)
           b(3,jb) = b(3,jb) - b(2,jb)*b(4,jb-1)
           b(4,jb) = b(4,jb) - b(2,jb)*b(5,jb-1)
          endif
         endif
    4   continue
    5  continue
      endif
      return
      end
c ------------ : ------------
      subroutine bansol(r,xx,b,n,int)
      include 'transin'
c
c  uses l.u factorization to solve for x, given r
c  
      real     r(nx),xx(nx),b(5,nx)

c
c  int = 1, tridiagonal system
c  int = 2, pentadiagonal system
c   assumes that first equation formed at midside node
c   ibc=0 if last equation formed at midside node
c   ibc=1 if last equation formed at corner node
c  
      if (int .eq. 1) then
       np = n-1
       do 1 j = 1,np
        jp = j+1
        r(jp) = r(jp) - b(2,jp)*r(j)
    1  continue
       xx(n)=r(n)/b(3,n)
       do 2 j = 1,np
        ja = n - j
        xx(ja) = ( r(ja) - b(4,ja)*xx(ja+1) )/b(3,ja)
    2  continue
       return
      else
    3  ibc = 1
       nh = n/2
       if (2*nh .eq. n) r(n+1) = 0.
       if (2*nh .eq. n) b(2,n+1) = 0.
       do 6 i = 1,2
        do 5 j = 1,nh
         ja = 2*j
         do 4 k = 1,i
          jb = ja -1 + k
          r(jb) = r(jb) - b(i,jb)*r(jb-1)
    4    continue
    5   continue
    6  continue
       nen = nh -ibc
       xx(n) = r(n)/b(3,n)
       if ( ibc .eq. 1) xx(n-1) = ( r(n-1) - b(4,n-1)*xx(n) )/b(3,n-1)
       do 7 j = 1,nen
        ja = n -2*j + 1 - ibc
        xx(ja) = ( r(ja) - b(4,ja)*xx(ja+1) - b(5,ja)*xx(ja+2) )/b(3,ja)
        xx(ja-1) = ( r(ja-1) - b(4,ja-1)*xx(ja) )/b(3,ja-1)
    7  continue
      endif
      return
      end