program   convect
	include 'convin'

c  1d simulation program using finite differences using 
c  1 - the upwind scheme
c  2 - the leapfrog scheme
c  3 - the Lax-Wendroff 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 are normalized fields for integration   
c   read input parameters and calculate runtime params:
      call pread
c   determine grid:
      call grid
      write(14) nx
      write(14) x
c   initialize time:      
      dt = c*dx/u
      nt = 0
      time = 0.0      
c   write some parameters:      
      call pwrite
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
       call bdcon
       if (mod(nt,nout) .eq. 0) call out
       do while( (nt .lt. ntmax) .and. (time .lt. 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
         call bdcon
c   write output f:      
        if (mod(nt,nout) .eq. 0) call out
       end do
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 'convin'
c   input data:
c   parameters explained in input file sim1.dat
       open(11,file='conv1.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) s,u,c
       read(11,1) nout
       read(11,3) f0
       read(11,1) init
       read(11,3) lambda,fb0(1),fb0(2),fb1(1),fb1(2)
       read(11,1) nterm
       close(11) 
       s0 = 1.0-c
       if (method .eq. 3 .or. method .eq. 5)  then
         s0 = 0.5*c
         s1 = 0.5*c*c
       endif
    1  format(i5)
    2  format(e10.2)
    3  format(f8.2)
      return
      end
c ------------ : ------------
      subroutine grid
      include 'convin'
       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 pwrite
      include 'convin'
       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)'
       write(16,2) nx,ntmax,tmax, xmin, xmax
       write(16,3) alpha,s,dt,dx,nout
       write(16,4) u,c
       write(16,5) f0,lambda,fb1(1),fb1(2),fb0(1),fb0(2)
       write(16,6) nterm
       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)'
       write(6,2) nx,ntmax,tmax, xmin, xmax
       write(6,3) alpha,s,dt,dx,nout
       write(6,4) u,c
       write(6,5) f0,lambda,fb0(1),fb0(2),fb1(1),fb1(2)
       write(6,6) nterm
    2  format(1x,'nx =',i4,'   ntmax =',i5,'     tmax =',f8.2,
     +       '  xmin =',f5.1,'  xmax =',f5.1)
    3  format(1x,'  alpha =',e10.3,' s =',f5.3,'  dt =',f8.3,
     +       '    dx =',f6.4,' nout =',i5)
    4  format(1x,' u=',f6.2,'  c =',f6.2,/)
    5  format(1x,'f0 =',f6.1,'  lambda =',f6.2,'  fxmin0 =',f5.2,
     +       ' fxmax0 =',f5.2,' fxmin =',f5.2,' fxmax =',f5.2,/)
       write(16,5) nterm
    6  format(1x,'No of terms in exact solution (nterm) =',i5)
      return
      end
c ------------ : ------------
      subroutine initcon1
      include 'convin'
       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 'convin'
       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 'convin'
       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 out
      return
      end
c ------------ : ------------
      subroutine bdcon
      include 'convin'
       f(1)  = f(nx-2)
       f(nx) = f(3)
      return
      end
c ------------ : ------------
      subroutine upwind
      include 'convin'
      integer  i
       do 10 i = 1,nx
   10   fold(i) = f(i)
       do 20 i = 2,nx-1
   20   f(i) = s0*fold(i) + c*fold(i-1)
      return
      end
c ------------ : ------------
      subroutine upleap
      include 'convin'
      integer  i
      write(6,*) 'nt:', nt, s0,c,f(2),f(3),f(4),f(5), f(6)
       do 10 i = 1,nx
   10   fold(i) = f(i)
       do 20 i = 2+mod(nt,2),nx-1,2
   20   f(i) = s0*fold(i) + c*fold(i-1) 
c   20   f(i) = fold(i) - 0.5*c*( fold(i+1)-fold(i-1) )
      write(6,*) 'nt:', nt, s0,c,f(2),f(3),f(4),f(5), f(6)
      return
      end
c ------------ : ------------
      subroutine leap
      include 'convin'
      integer  i
       do 20 i = 2+mod(nt,2),nx-1,2
   20   f(i) = f(i) - c*( f(i+1)-f(i-1) )
      return
      end
c ------------ : ------------
      subroutine laxwend
      include 'convin'
      integer  i
       do 10 i = 1,nx
   10   fold(i) = f(i)
       do 30 i = 2,nx-1
   30   f(i) = fold(i) - s0*( fold(i+1)-fold(i-1) ) 
     +            +s1*( fold(i+1)+fold(i-1)-2.*fold(i) )
      return
      end
c ------------ : ------------
      subroutine out
      include 'convin'
      integer  i
      real     fout(nx),foutex(nx)
       do 10 i = 1,nx
        foutex(i) = f0*fexact(i)
   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
       write(14) time,rms,fout,foutex
       write(*,*) 'time step:', nt
    1  format(' time = ',f5.0,'  T =', 11f7.2)
      return
      end
c ------------ : ------------
      subroutine analyt1
      include 'convin'
      integer  i
      real     sum
       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)*alpha*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
       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 'convin'
      integer  i
      real     sum
       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*alpha*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
       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 'convin'
      integer  i
      real     sum
       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)*alpha*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
       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