program   sim1d
	include 'simin'

c  1d simulation program using finite differences
c  the particular scheme is specified in the integration routine
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   input parameters:
      call pread
c      s=1./6.
c      d=1.-1./12./s
      write(*,*) 'd1:',d
c   determine grid:
      call grid
      write(14) nx
      write(14) x
c   initialize time:      
      dt = dx*dx*s/alpha
      nt = 0
      time = 0.0      
      write(14) s,dt
      write(*,*) 'Timestep dt:', dt,'  Output every ',nout*dt,'  times'
c   write some parameters:      
      call pwrite
      write(*,*) 'd2:',d
c   set initial conditions:      
      call initcon
c   integration: 
c    intftcs - 1d heat conduction with ftcs 
c              (forward time centered space)
        nt = nt+1
        time = time + dt
        call intftcs
        call bdcon
        if (mod(nt,nout) .eq. 0) call out
      do while( (nt .lt. ntmax) .and. (time .lt. (tmax-0.001*dt)) )
        nt = nt+1
        time = time + dt
        write(*,*) 'time', time
        call int2l
        call bdcon
c   write output f:      
        if (mod(nt,nout) .eq. 0) call out
      end do
      write(*,*) 'd3:',d
c   determine and write exact solution to order nterm:  
      call analyt
      close(14) 
      close(16) 
      stop
      end
c ------------ : ------------
      subroutine pread
      include 'simin'
c   input data:
c   parameters explained in input file sim1.dat
       open(11,file='sim1.dat')          
       read(11,1) ntmax
       read(11,3) tmax
       read(11,3) xmin,xmax
       read(11,2) alpha
       read(11,3) s
       read(11,3) d
       read(11,1) nout
       read(11,3) f0,fini,fb0(1),fb0(2),fb1(1),fb1(2)
       read(11,1) nterm
       close(11) 
    1  format(i5)
    2  format(e10.2)
    3  format(f8.2)
      return
      end
c ------------ : ------------
      subroutine grid
      include 'simin'
       integer i
       dx = (xmax-xmin)/(nx-1.)
       do 10 i = 1,nx
   10    x(i)  = xmin+dx*(i-1.)
      return
      end
c ------------ : ------------
      subroutine pwrite
      include 'simin'
       write(16,1) 
       write(16,2) nx,ntmax,tmax, xmin, xmax
       write(16,3) s,d,alpha,dt,dx,nout
       write(16,4) f0,fini,fb1(1),fb1(2),fb0(1),fb0(2)
    1  format(1x,'Two level scheme (explicit)',//)
    2  format(1x,'nx =',i4,'   ntmax =',i5,'     tmax =',f8.2,
     +       '  xmin =',f5.1,'  xmax =',f5.1)
    3  format(1x,' s =',f5.3,' d =',f5.3,'  alpha =',e10.3,
     +           '  dt =',f8.3,'    dx =',f6.4,' nout =',i5)
    4  format(1x,'f0 =',f6.1,'  fini =',f6.2,'  fxmin =',f5.2,
     +       ' fxmax =',f5.2,' fxmin0 =',f5.2,' fxmax0 =',f5.2)
       write(16,5) nterm
    5  format(1x,'No of terms in exact solution (nterm) =',i5/)
      return
      end
c ------------ : ------------
      subroutine initcon
      include 'simin'
        integer i
       do 10 i = 1,nx-1
  10     f(i) = fini
       call bdcon
       call out
      return
      end
c ------------ : ------------
      subroutine bdcon
      include 'simin'
       f(1)  = fb1(1)
       f(nx) = fb1(2)
c       if(nt .eq. 1) f(1)  = 2.*fb0(1)
c       if(nt .eq. 1) f(nx) = 2.*fb0(2)
       if(nt .eq. 0) f(1)  = fb0(1)
       if(nt .eq. 0) f(nx) = fb0(2)
      return
      end
c ------------ : ------------
      subroutine intftcs
      include 'simin'
      integer  i
      real     s0
       s0 = 1.0-2.*s
       do 10 i = 1,nx
   10   fold(i) = f(i)
       do 20 i = 2,nx-1
   20   f(i) = s0*fold(i) + s*( fold(i+1)+fold(i-1) )
      return
      end
c ------------ : ------------
      subroutine int2l
      include 'simin'
      integer  i
      real     s00,s10,s01,s11
       s10 = 2.*s*(1.+d)/3.
       s00 = 4./3.-2.*s10
       s11 = -2./3.*s*d
       s01 = -1./3.-2.*s11
       do 10 i = 1,nx
   10   foold(i) = fold(i)
       do 20 i = 1,nx
   20   fold(i) = f(i)
       do 30 i = 2,nx-1
   30   f(i) = s00*fold(i) + s10*( fold(i+1)+fold(i-1) ) 
     +            + s01*foold(i) + s11*( foold(i+1)+foold(i-1) ) 
      return
      end
c ------------ : ------------
      subroutine out
      include 'simin'
      integer  i
      real     fout(nx)
       do 10 i = 1,nx
   10   fout(i) = f0*f(i)
       write(16,1) time
       write(16,2) (fout(i),i=1,nx)
    1  format(' Time = ',f7.0)
    2  format('   T:', 11f7.2)
       write(14) time,fout
       write(*,*) 'time step:', nt
      return
      end
c ------------ : ------------
      subroutine analyt
      include 'simin'
      integer  i,m
      real     dam,arg1,arg2,sum,rms
       sum = 0.0
       do 20 i = 1,nx
         fexact(i) = f0
         do 10 m = 1,nterm
           dam = 2.*m-1.
           arg1 = dam*pi*x(i) 
           arg2 = -alpha*dam*dam*pi*pi*time
c   limit argument of exp(dtm)
           if (arg2 .lt. -87.) arg2 = -87.0
           fexact(i) = fexact(i) - 4.*f0/dam/pi*sin(arg1)*exp(arg2)
   10    continue
         sum = sum + ( fexact(i)-f0*f(i) )**2
   20  continue
c   rms is the rms error
       rms = sqrt(sum/nx)
       write(16,1) nterm
       write(16,2) time
       write(16,3) (fexact(i),i=1,nx)
       write(16,4) rms
       write(14) time,fexact
    1  format(1x,/'No of terms in exact solution (nterm) =',i5)
    2  format(' time = ',f7.0)
    3  format('  T:', 11f7.2)
    4  format(1x,' rms dif = ', e11.4,//)
      return
      end