program adipr

c General information:

c     Program adipr solves the two-dimensional convection equation
c
c         u_t + ax u_x + ay u_y =0
c
c     for
c
c         0 <= x <= xmax  0 <= y <= ymax  0 <= t <= tmax.

c     For simplicity it is assumed that the convection velocity
c     components ax and ay are positive, so that the boundary conditions
c     apply at x = 0 and at y = 0.  Reverse x or y to handle negative ax
c     or ay.

c     Program adipr produces the output file "adipr.out". Any existing file of
c     that name will be deleted.

c     Copyright 1997 Leon van Dommelen
c     Version 1.0 Leon van Dommelen 1/10/97

c Numerical scheme:

c     Program adipr uses the two-stage Peaceman-Rachford ADI method:
c
c      (1 + 0.5 dt ax d/dx) u-tilde = ( 1 - 0.5 dt ay d/dy) u^n
c      (1 + 0.5 dt ay d/dy) u^n+1   = ( 1 - 0.5 dt ax d/dx) u-tilde
c
c     where u^n is the dependent variable u at the last known time,
c     u^n+1 is the dependent variable at the next time to be computed;
c     dt is the numerical time step; u-tilde a purely computational
c     "intermediate unknown", and the derivatives d/dx and d/dy are
c     evaluated using standard second order central differences.
c     This method is second order accurate in space and time.

c     A boundary condition for u-tilde can be obtained by substracting
c     the two Peaceman-Rachford formulae from each other, to give:
c
c      2 u-tilde = (1 - 0.5 dt ay d/dy) u^n + (1 + 0.5 dt ay d/dy) u^n+1

c     Because of the central differencing, the scheme requires
c     artificial boundary conditions at x = xmax and at y = ymax.

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

c     Program adipr requires a separate subroutine
c           ini(jdim,ldim,io,
c               ax,ay,xmax,ymax,tmax,jmax,lmax,nmax,x,y)
c     to initialize the various variables.

c     Program adipr requires a separate subroutine
c           ic(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,io,u)
c     to set the initial condition.

c     Program adipr requires a separate subroutine
c           bcx(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,io,
c               abcx,bbcx)
c     to return the boundary condition at x = 0 as equations of the form
c           sum (j = 0 to jmax) abcx(j) u(j,l) = bbcx(l).

c     Program adipr requires a separate subroutine
c           bcy(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,io,
c               abcy,bbcy)
c     to return the boundary condition at y = 0 as equations of the form
c           sum (l = 0 to lmax) abcy(l) u(j,l) = bbcy(j).

c     Program adipr requires a separate subroutine
c           bcxa(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,io,
c                abcxa,bbcxa)
c     to return an artificial boundary condition at x = xmax as equations
c     of the form
c           sum (j = 0 to jmax) abcxa(j) u(j,l) = bbcxa(l).

c     Program adipr requires a separate subroutine
c           bcya(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,io,
c                abcya,bbcya)
c     to return an artificial boundary condition at y = ymax as equations
c     of the form
c           sum (l = 0 to lmax) abcya(l) u(j,l) = bbcya(j).

c     Program adipr requires a separate subroutine
c           out(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,u,io)
c     to output the solution at time t.

c     Program adipr requires a separate subroutine:
c           tribc(ndim,n,abc,abca,a1,a2,a3,b,
c                 w1,w2,w3,w4,w5,w6,w7,w8,w9,
c                 u,ier,err,umax)
c     to solve an almost tridiagonal system of n+1 equations of the form
c           sum (j = 0 to n) abc(j)  u(j)       = b(0)  first equation;
c      a1(j) u(j-1) + a2(j) u(j) + a3(j) u(j+1) = b(j)  for j = 1, ..., n-1;
c           sum (j = 0 to n) abca(j) u(j)       = b(n)  final equation;
c     for u(0), u(1), ... u(n). Arguments w1 through w9 are work space;
c     ier is an error code which should be less than two after tribc;
c     err is an error estimate and umax is the maximum |u(i)| value.
      external tribc

c Program variables:

c     Avoid typos:
      implicit none

c     Components of the convection velocity in the x and y directions:
      double precision ax,ay

c     Size of the domain in the x- and y-directions:
      double precision xmax,ymax

c     The maximum time to compute:
      double precision tmax

c     Index in the x-direction with its maximum and dimensioned value:
      integer j,jmax,jdim
      parameter (jdim=128)

c     Index in the y-direction with its maximum and dimensioned value:
      integer l,lmax,ldim
      parameter (ldim=64)

c     Time step number and its maximum value:
      integer n,nmax

c     Mesh point x- and y-positions
      double precision x(0:jdim),y(0:ldim)

c     Current time:
      double precision t

c     Spatial mesh size in the x direction:
      double precision dx

c     Spatial mesh size in the y direction:
      double precision dy

c     The time step:
      double precision dt

c     The values of the dependent variable for the last computed time:
      double precision u(0:jdim,0:ldim)

c     The intermediate unknowns u-tilde in the ADI method:
      double precision utilde(0:jdim,0:ldim)

c     Storage for the equations during the first stage of the
c     ADI method listed above under "Numerical method":

c     The coefficients and the right hand sides of the tridiagonal
c     equations in the x-direction for u-tilde for the interior points:
      double precision a1x(0:jdim),a2x(0:jdim),a3x(0:jdim),bx(0:jdim)

c     The coefficients and the right hand sides of the equations
c     describing the boundary condition at x = 0 for u as returned by
c     subroutine bcx:
      double precision abcx(0:jdim),bbcx(0:ldim)

c     The right hand sides of the equations describing the boundary
c     condition at x = 0 for the intermediate unknown u-tilde:
      double precision bbcxt(0:ldim)

c     The coefficients and right hand sides of the equations describing
c     the artificial boundary condition at x = xmax:
      double precision abcxa(0:jdim),bbcxa(0:ldim)

c     The right hand sides of the equations describing the artificial
c     boundary condition at x = xmax for the unknown u-tilde:
      double precision bbcxat(0:ldim)

c     The solution vector of the equations as found by subroutine tribc:
      double precision solx(0:jdim)

c     Work vectors used by tribc while solving the equations:
      double precision wx1(0:jdim),wx2(0:jdim),wx3(0:jdim),wx4(0:jdim),
     & wx5(0:jdim),wx6(0:jdim),wx7(0:jdim),wx8(0:jdim),wx9(0:jdim)

c     Storage for the equations during the second stage of the
c     ADI method listed above under "Numerical method":

c     The coefficients and the right hand sides of the tridiagonal
c     equations in the y-direction for u for the interior points:
      double precision a1y(0:ldim),a2y(0:ldim),a3y(0:ldim),by(0:ldim)

c     The coefficients and the right hand sides of the equations
c     describing the boundary condition at y = 0 as returned by bcy:
      double precision abcy(0:ldim),bbcy(0:jdim)

c     The coefficients and right hand sides of the equations describing
c     the artificial boundary condition at y = ymax:
      double precision abcya(0:ldim),bbcya(0:jdim)

c     The solution vector of the equations as found by subroutine tribc:
      double precision soly(0:ldim)

c     Work vectors used by tribc while solving the equations:
      double precision wy1(0:ldim),wy2(0:ldim),wy3(0:ldim),wy4(0:ldim),
     & wy5(0:ldim),wy6(0:ldim),wy7(0:ldim),wy8(0:ldim),wy9(0:ldim)

c     Error code returned by subroutine tribc:
      integer ier
c     Estimated error made by tribc:
      double precision err,errmax
c     Maximum u value found by subroutine tribc:
      double precision umax

c     Whether an old output file exists:
      logical exists

c     Constants, defined to simplify changing precision:
      double precision zero,one,two,o2
      parameter (zero=0.d0,one=1.d0,two=2.d0,o2=0.5d0)

c Initialize the variables:

c     Create the output file after deleting any existing old file:
      inquire(file='adipr.out',exist=exists)
      if(exists)then
         open(1,file='adipr.out',status='old')
         close(1,status='delete')
      endif
      open(1,file='adipr.out',status='new')

c     Let subroutine ini initialize the computational variables:
      call ini(jdim,ldim,1,
     & ax,ay,xmax,ymax,tmax,jmax,lmax,nmax,x,y)

c     Check for the minimum mesh size:
      if(jmax.lt.1 .or. lmax.lt.1)then
         print*,'*** adipr: the mesh size must be at least 1x1'
         stop
      endif

c     Set the numerical time step:
      dt=one
      if(nmax.gt.0)dt=tmax/nmax

c     Set the horizontal mesh size:
      dx=x(1)-x(0)

c     Set the vertical mesh size:
      dy=y(1)-y(0)

c     Initialize the time step:
      n=0

c     Initialize the time:
      t=zero

c     Let subroutine ic set the initial condition on u:
      call ic(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,1,u)

c     Get the boundary condition at x = 0 from subroutine bcx:
      call bcx(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,1,
     & abcx,bbcx)

c     Get the artificial boundary condition at x = xmax from subroutine bcxa:
      call bcxa(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,1,
     & abcxa,bbcxa)

c     Get the boundary condition at y = 0 from subroutine bcy:
      call bcy(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,1,
     & abcy,bbcy)

c     Get the artificial boundary condition at y = ymax from subroutine bcya:
      call bcya(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,1,
     & abcya,bbcya)

c     Set the y = 0 and ymax boundary values according to the end conditions:
      do 200 j=1,jmax-1
         u(j,0)=bbcy(j)
         do 150 l=1,lmax
            u(j,0)=u(j,0)-abcy(l)*u(j,l)
 150     continue
         u(j,0)=u(j,0)/abcy(0)
         u(j,lmax)=bbcya(j)
         do 160 l=0,lmax-1
            u(j,lmax)=u(j,lmax)-abcya(l)*u(j,l)
 160     continue
         u(j,lmax)=u(j,lmax)/abcya(lmax)
 200  continue

c     Set the x = 0 and xmax boundary values according to the end conditions:
      do 300 l=0,lmax
         u(0,l)=bbcx(l)
         do 250 j=1,jmax
            u(0,l)=u(0,l)-abcx(j)*u(j,l)
 250     continue
         u(0,l)=u(0,l)/abcx(0)
         u(jmax,l)=bbcxa(l)
         do 260 j=0,jmax-1
            u(jmax,l)=u(jmax,l)-abcxa(j)*u(j,l)
 260     continue
         u(jmax,l)=u(jmax,l)/abcxa(jmax)
 300  continue

c     Let subroutine out output the initial plane:
      call out(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,u,1)

c     Initialize the maximum error in subroutine tribc alone:
      errmax=zero

c Start of the computation of times greater than zero:

c     For n = 0, 1, ..., nmax-1, compute the time-plane n+1
      do 5000 n=0,nmax-1

c        Find the next time:
         t=(n+1+zero)/nmax*tmax
c           The addition of zero assures that the first division is a
c        floating point one.  The order of the terms ensures that the
c        last time plane is exactly time tmax, regardless of round-off
c        errors.  Otherwise subroutine out might not output it.

c     Determination of the intermediate unknown u-tilde in the ADI
c     procedure listed above under "Numerical method":

c        The boundary condition at x = 0 for u-tilde is of similar form
c        as the one for u, namely,
c           sum (j = 0 to jmax) abcx(j) utilde(j,l) = bbcxt(l),
c        and the artificial boundary condition at x = xmax is
c           sum (j = 0 to jmax) abcxa(j) utilde(j,l) = bbcxat(l).
c        From the expression for u-tilde in terms of u at time steps n
c        and n+1 as given under "Numerical method", the right hand side
c        bbcxt is found in terms of the right hand sides bbcx for u as:
c           bbcxt = 0.5 (1 - 0.5 dt ay d/dy) bbcx^n
c                   0.5 (1 + 0.5 dt ay d/dy) bbcx^n+1
c        and similar for bbcxat.

c        Contribution of time step n to the right hand sides for u-tilde:
         do 2100 l=1,lmax-1
            bbcxt(l)=
     &       o2*(bbcx(l)-o2*dt/(two*dy)*(bbcx(l+1)-bbcx(l-1)))
            bbcxat(l)=
     &       o2*(bbcxa(l)-o2*dt/(two*dy)*(bbcxa(l+1)-bbcxa(l-1)))
 2100    continue

c        Get the new boundary condition at x = 0:
         call bcx(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,1,
     &    abcx,bbcx)

c        Get the artificial boundary condition at x = xmax:
         call bcxa(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,1,
     &    abcxa,bbcxa)

c        Contribution of time step n+1 to the right hand sides for u-tilde:
         do 2200 l=1,lmax-1
            bbcxt(l)=bbcxt(l)+
     &       o2*(bbcx(l)+o2*dt/(two*dy)*(bbcx(l+1)-bbcx(l-1)))
            bbcxat(l)=bbcxat(l)+
     &       o2*(bbcxa(l)+o2*dt/(two*dy)*(bbcxa(l+1)-bbcxa(l-1)))
 2200    continue

c        Now find u-tilde at the interior points:

         do 3000 l=1,lmax-1

c           Find the equations for u-tilde at the interior points:

            do 2500 j=1,jmax-1

c              The subdiagonal coefficient (see "Numerical scheme" above):
               a1x(j)=-o2*dt*ax/(two*dx)
c              The diagonal coefficient (see "Numerical scheme" above):
               a2x(j)=one
c              The superdiagonal coefficient (see "Numerical scheme" above):
               a3x(j)=-a1x(j)
c              The right hand side (see "Numerical scheme" above):
               bx(j)=u(j,l)-o2*dt*ay/(two*dy)*(u(j,l+1)-u(j,l-1))

 2500       continue

c           Set the right-hand sides corresponding to the end conditions:
            bx(0)=bbcxt(l)
            bx(jmax)=bbcxat(l)

c           Hand the system of equation to subroutine tribc to solve:
            call tribc(jdim,jmax,abcx,abcxa,a1x,a2x,a3x,bx,
     &       wx1,wx2,wx3,wx4,wx5,wx6,wx7,wx8,wx9,
     &       solx,ier,err,umax)
            if(ier.ne.0)then
               print*,'*** adipr: fatal error code from tribc:',
     &          ier,err,umax
               stop
            endif
            errmax=max(errmax,err)

c           Copy the solution into u-tilde:
            do 2900 j=0,jmax
               utilde(j,l)=solx(j)
 2900       continue

 3000    continue

c        Get the new boundary condition at y = 0 from subroutine bcy:
         call bcy(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,1,
     &    abcy,bbcy)

c        Get the artificial boundary condition at y = ymax:
         call bcya(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,1,
     &    abcya,bbcya)

c        Compute the new dependent variable along the interior mesh lines:

         do 4000 j=1,jmax-1

c           Find the equations for u at the interior points:

            do 3500 l=1,lmax-1

c              The subdiagonal coefficient (see "Numerical scheme" above):
               a1y(l)=-o2*dt*ay/(two*dy)
c              The diagonal coefficient (see "Numerical scheme" above):
               a2y(l)=one
c              The superdiagonal coefficient (see "Numerical scheme" above):
               a3y(l)=-a1y(l)
c              The right hand side (see "Numerical scheme" above):
               by(l)=utilde(j,l)
     &          -o2*dt*ax/(two*dx)*(utilde(j+1,l)-utilde(j-1,l))

 3500    continue

c           Set the right-hand sides corresponding to the end conditions:
            by(0)=bbcy(j)
            by(lmax)=bbcya(j)

c           Hand the system of equation to subroutine tribc to solve:
            call tribc(ldim,lmax,abcy,abcya,a1y,a2y,a3y,by,
     &       wy1,wy2,wy3,wy4,wy5,wy6,wy7,wy8,wy9,
     &       soly,ier,err,umax)
            if(ier.ne.0)then
               print*,'*** adipr: fatal error code from tribc:',
     &          ier,err,umax
               stop
            endif
            errmax=max(errmax,err)

c           Copy the solution into u:
            do 3900 l=0,lmax
               u(j,l)=soly(l)
 3900       continue

 4000    continue

c        Find the new dependent variable on the boundaries:

         do 4100 l=0,lmax
            u(0,l)=bbcx(l)
            do 4050 j=1,jmax
               u(0,l)=u(0,l)-abcx(j)*u(j,l)
 4050       continue
            u(0,l)=u(0,l)/abcx(0)
            u(jmax,l)=bbcxa(l)
            do 4060 j=0,jmax-1
               u(jmax,l)=u(jmax,l)-abcxa(j)*u(j,l)
 4060       continue
            u(jmax,l)=u(jmax,l)/abcxa(jmax)
 4100    continue

c        Allow subroutine out to perform output:
         call out(ax,ay,xmax,ymax,jdim,jmax,ldim,lmax,x,y,t,u,1)

c        Go compute the next plane:

 5000 continue

 9000 print*,'Maximum error in solving the equations alone:',errmax
      write(1,9010)errmax
 9010 format('Maximum error in solving the equations alone:',e11.2)

      stop
      end