program ftcs

c General information:

c     Explicit solution of unsteady heat conduction in a uniform bar.

c     The equation solved is
c                          u_t = kappa u_xx
c                   for 0 < x < l and 0 < t < tmax
c     where u is the temperature, l is the length of the bar and tmax is
c     the desired final time to be computed.

c     The numerical scheme is the Euler, or forward-time central-space,
c     method.

c     Copyright 1996 Leon van Dommelen
c     Version 1.0 Leon van Dommelen 12/16/96

c Usage information:

c     Program ftcs looks for its parameters in a file 'ftcs.in' of the
c     form:

c       <Value of kappa>            <comment>
c       <Value of l>                <comment>
c       <Value of tmax>             <comment>
c       <Value of jmax>             <comment>

c     where kappa, l and tmax are as above and jmax is the number of mesh
c     intervals to use along the bar.
c        For better readability, comment lines may be inserted anywhere
c     in file ftcs.in.  All lines are considered comments, unless the
c     first nonblank character is a digit, sign, or a decimal point.

c     During execution, the program will prompt for a numerical time
c     step to use.  It will list some possible choices for this time
c     step.  The program may slightly adjust the chosen time step to
c     exactly reproduce the desired ending time tmax.

c     Program output is written to output file 'ftcs.out'. Note that
c     the output subroutine out may define additional output files.

c External variables and info for compiling or changing program ftcs:

c     The next statement can be removed if your compiler does not support it:
      implicit none

c     Program ftcs requires a separate subroutine:
c           initx(l,jdim,jmax,ndim,io,x,t,u)
c     to determine the initial mesh point locations x(j).
c     Subroutine initx can be found in file:
c           ../../init/init_x.f
      external initx

c     Program ftcs requires a separate subroutine:
c           initu(kappa,l,jdim,jmax,ndim,n,x,t,io,task,u)
c     to set the initial values for the temperature u(j,0).
c     For an example subroutine initu see:
c           ../../init/exact/init_u.f
c     For a template for writing new subroutines initu see:
c           ../../init/template/init_u.f
      external initu

c     Program ftcs requires a separate subroutine:
c           bc1(kappa,l,jdim,jmax,ndim,n,x,t,io,task,u)
c     to set the left hand boundary condition.
c     For an example subroutine bc1 see:
c           ../../bc_1/diri_exa/bc_1.f
c     For a template for writing new subroutines bc1 see:
c           ../../bc_1/template/bc_1.f
      external bc1

c     Program ftcs requires a separate subroutine:
c           bc2(kappa,l,jdim,jmax,ndim,n,x,t,io,task,u)
c     to set the right hand boundary condition.
c     For an example subroutine bc2 see:
c           ../../bc_2/diri_exa/bc_2.f
c     For a template for writing new subroutines bc2 see:
c           ../../bc_2/template/bc_2.f
      external bc2

c     Program ftcs requires a separate subroutine:
c           out(kappa,l,jdim,jmax,ndim,n,x,t,u,io,task)
c     to handle output of the results.
c     For an example subroutine out see:
c           ../../out/coef_exa/out.f
c     For a template for writing new subroutines out see:
c           ../../out/template/out.f
      external out

c     Program ftcs uses the following utility routines from
c     ../../../lib/util.f:

c        exit0(module) exits the program after successful completion.
      external exit0

c        fatal(module,text1,text2,text3) kills the program after a
c     fatal error, printing the lines "text1", "text2" and "text3".
      external fatal

c        fclose(module,filnam,io) closes the I/O unit "io" on which
c     file "filnam" is open.
      external fclose

c        fcreat(module,filnam,io) creates the new output file "filnam",
c     opening it on I/O unit "io". It selects "io" if initialized to zero.
      external fcreat

c        fopen(module,filnam,io) opens the existing file "filnam" on
c     I/O unit "io". It selects "io" if initialized to zero.
      external fopen

c        iread(module,text,ioin,io,show) reads integer "iread" from I/O
c     unit "ioin". It writes it and its description "text" to I/O unit
c     "io" if nonzero, and to the screen if "show" is nonzero.
      integer iread
      external iread

c        iwrite(module,text,val,io,show) writes integer "val" and its
c     description "text" to I/O unit "io" if nonzero, and to the screen
c     if "show" is nonzero.
      external iwrite

c        rask(module,text,io) gets the number "rask" from the user,
c     using the prompt "text", and writes it to I/O unit "io" if nonzero.
      double precision rask
      external rask

c        rread(module,text,ioin,io,show) reads number "rread" from I/O
c     unit "ioin". It writes it and its description "text" to I/O unit
c     "io" if nonzero, and to the screen if "show" is nonzero.
      double precision rread
      external rread

c        rwrite(module,text,val,io,show) writes number "val" and its
c     description "text" to I/O unit "io" if nonzero, and to the screen
c     if "show" is nonzero.
      external rwrite

c        warn(module,io,text1,text2,text3) writes the warning lines
c     "text1", "text2" and "text3" to the screen, and to I/O unit "io"
c     if nonzero.
      external warn

c Program variables:

c     The heat conduction coefficient:
      double precision kappa

c     The desired maximum time to compute to:
      double precision tmax
c     The size of the time step:
      double precision dt

c     The time step number n:
      integer n
c     The maximum value of n, corresponding to time tmax:
      integer nmax

c     Program ftcs declares only storage for two time levels.  This is
c     achieved by storing using index nn=mod(n,2) instead of index n
c     itself. (Note that mod(n,2) equals 0 or 1 depending on whether n
c     is even or odd.)

c     The value of nn:
      integer nn
c     The value of nn corresponding to n+1:
      integer nnp1

c     The array of mesh times, with only storage for 2 time levels (see
c     the comments above):
      double precision t(0:1)

c     The length of the bar:
      double precision l
c     The size of the spatial step:
      double precision dx

c     The mesh point index:
      integer j
c     The maximum mesh point index:
      integer jmax
c     The dimensioned maximum mesh point index:
      integer jdim
      parameter (jdim=128)

c     The x-values of the mesh points:
      double precision x(0:jdim)
c        The values x(j) are set by subroutine initx, which will take
c     the implementation of the boundary conditions into account.  The
c     value of x(0) will be approximately zero, depending on the left
c     boundary condition, and the value of x(jmax) will be approximately
c     l, depending on the right boundary condition.

c     The temperature values of the mesh points, with only storage for
c     2 time levels (see the comments above at the declaration of nn):
      double precision u(0:jdim,0:1)

c     The Fortran I/O unit number for the output file:
      integer io

c     The Fortran I/O unit number for the input file:
      integer ioin

c     Various constants, defined to make conversion of precision simpler:
      double precision zero,one,two,six,o2
      parameter (zero=0.d0,one=1.d0,two=2.d0,six=6.d0,o2=.5d0)

c Executable statements:

c Open the output file ftcs.out on I/O unit io using fcreat of util.f:

      io=0
      call fcreat('ftcs','ftcs.out',io)

c Open the input file ftcs.in on I/O unit ioin using fopen of util.f:

      ioin=0
      call fopen('ftcs','ftcs.in',ioin)

c Read the input file:

c     Read the heat conduction coefficient:
      kappa=rread('ftcs',
     &'The heat conduction coefficient kappa',ioin,io,1)
c     Allow negative kappa (the backward heat equation):
      if(kappa.lt.zero)call warn('ftcs',io,
     &'Negative heat conduction coefficient.',' ',' ')

c     Read the length of the bar:
      l=rread('ftcs','The length of the bar l',ioin,io,1)
c     Enforce a positive bar length:
      if(l.le.zero)call fatal('ftcs',
     &'The length of the bar must be positive.',' ',' ')

c     Final time to compute:
      tmax=rread('ftcs','The maximum time to compute',ioin,io,1)
c     Do not allow backward marching:
      if(tmax.lt.zero)call fatal('ftcs',
     &'The maximum time cannot be negative.',' ',' ')
c     Allow tmax to be zero to compute the initial plane only:
      if(tmax.eq.zero)call warn('ftcs',io,
     &'Only initial data will be computed.',' ',' ')

c     Read the number of mesh intervals along the bar to use:
      jmax=iread('ftcs','The number of mesh intervals',ioin,io,1)
c     Check whether we have enough storage available:
      if(jmax.gt.jdim)then
         call iwrite('ftcs','    ftcs: '//
     &   '  Requested number of mesh intervals',jmax,io,1)
         call iwrite('ftcs','          '//
     &   'Dimensioned number of mesh intervals',jdim,io,1)
         call fatal('ftcs','The requested number of mesh intervals '//
     &   'exceeds the program limit.',' ',' ')
      endif
c     Less than one interval is not supported:
      if(jmax.lt.1)call fatal('ftcs',
     &'The number of mesh intervals should be at least one.',' ',' ')

c     Done with input
      call fclose('ftcs','ftcs.in',ioin)

c Initialize the program variables:

c     Let subroutine initx create the mesh point positions x(j):
      call initx(l,jdim,jmax,1,io,x,t,u)

c     Compute the mesh size:
      dx=x(1)-x(0)
      call rwrite('ftcs','Mesh size dx',dx,io,0)

c     Let rask read the time step from the keyboard:
      if(kappa.ne.zero)then
         call rwrite('ftcs',
     &   'The time step for 4th order accuracy',dx**2/(six*kappa),io,1)
         call rwrite('ftcs',
     &   '        The largest stable time step',dx**2/(two*kappa),io,1)
      endif
 400  dt=one
      if(tmax.gt.zero)dt=rask('ftcs','the time step (0 exits)',io)
      if(dt.eq.zero)call fatal('ftcs',
     &'Execution aborted by the user.',' ',' ')
      if(dt.lt.zero)then
         call warn('ftcs',io,
     &   'The time step must be positive.',
     &   'The entered negative value has been refused.',' ')
         goto 400
      endif

c     Compute the rounded number of time steps:
      nmax=tmax/dt+o2
      if(tmax.le.zero)nmax=0
      if(tmax.gt.zero)nmax=max(nmax,1)
      call iwrite('ftcs','The total number of time steps',nmax,io,1)
c     Adjust the time step slightly to give an integer amount of steps:
      if(nmax.gt.0)then
         dt=tmax/nmax
         call rwrite('ftcs','The actual time step to be used',dt,io,1)
      endif

c     Output the amount of work to be done:
      call iwrite('ftcs',
     &'The total number of points to be computed',jmax*nmax,io,1)

c     Start at time zero:
      t(0)=zero

c     Make sure that the various subroutines do not mistakenly think the
c     other level 1 has meaningful information for an earlier time:
      t(1)=t(0)+dt

c     Allow subroutine init to set the initial data u(j,0):
      call initu(kappa,l,jdim,jmax,1,0,x,t,io,0,u)

c     Allow subroutine bc1 to set the boundary condition at j=0:
      call bc1(kappa,l,jdim,jmax,1,0,x,t,io,0,u)

c     Allow subroutine bc2 to set the boundary condition at j=jmax:
      call bc2(kappa,l,jdim,jmax,1,0,x,t,io,0,u)

c     Perform output:
      call out(kappa,l,jdim,jmax,1,0,x,t,u,io,0)

c March the solution in time to the final time tmax:

c     For n= 0 to nmax-1, compute the time plane n+1:

      do 5000 n=0,nmax-1

c        Use the Fortran mod (remainder) function to reduce n and n+1
c        to the only two available storage levels in array u:
         nn=mod(n,2)
         nnp1=mod(n+1,2)

c        Find the time at step n+1:
         t(nnp1)=(n+1+zero)/(nmax+zero)*tmax
c           Written this way, the last time will be exactly tmax despite
c        any round-off errors.  So if the last requested output time
c        equals tmax, it will not be missed because of round-off.

c        Compute all internal temperatures using the finite difference formula:
         do 2200 j=1,jmax-1
            u(j,nnp1)=u(j,nn)+
     &       kappa*dt/dx**2*(u(j+1,nn)-two*u(j,nn)+u(j-1,nn))
 2200    continue

c        Allow subroutine bc1 to set the boundary condition at j=0:
         call bc1(kappa,l,jdim,jmax,1,nnp1,x,t,io,0,u)

c        Allow subroutine bc2 to set the boundary condition at j=jmax:
         call bc2(kappa,l,jdim,jmax,1,nnp1,x,t,io,0,u)

c        Allow subroutine out to perform any required output:
         call out(kappa,l,jdim,jmax,1,nnp1,x,t,u,io,0)

 5000 continue

c Exit:

c     Close the output file:
      call fclose('ftcs','ftcs.out',io)

c     Exit:
      call exit0('ftcs')

c     For picky compilers:
      stop
      end