subroutine tribc(ndim,n,v0,vn,a,b,c,d,
     & u0,u1,u2,dstore,w1,w2,w3,w4,w5,
     & u,ier,err,umax)

c General information:

c     Subroutine tribc solves a tridiagonal system of equations for
c     the n+1 unknowns u(0), u(1), ..., u(n-1), u(n) with nontridiagonal
c     boundary conditions.

c     The system to be solved must consist of n+1 equations of the form:
c
c        sum (i = 0 to n)    v0(i) u(i)     = d(0)   first equation
c     a(i) u(i-1) + b(i) u(i) + c(i) u(i+1) = d(i)   for i = 1, 2, ..., n-1
c        sum (i = 0 to n)    vn(i) u(i)     = d(n)   final equation

c     Subroutine tribc uses the Thomas algorithm. This method does not
c     use pivoting; it may crash or be inaccurate if the tridiagonal
c     system is not strictly diagonally dominant or symmetric positive
c     definite.

c     Always check that the error code ier is zero after using tribc.

c     Subroutine tribc will crash if ndim < 0 or n < 0 or n > ndim.

c     Copyright 1996 Leon van Dommelen
c     Version 1.0 Leon van Dommelen 1/7/97

c Computational procedure:

c     Subroutine tribc works by solving the intermediate tridiagonal
c     equations using three different Dirichlet end conditions. Then it
c     combines the three solutions to also satisfy the first and final
c     equations.

c Arguments:

c     Avoid typos:
      implicit none

c     Input: declared maximum equation number (0 <= ndim):
      integer ndim

c     Input: actual last equation number (0 <= n <= ndim):
      integer n

c     Input: the coefficients of the first equation:
      double precision v0(0:ndim)

c     Input: the coefficients of the final equation:
      double precision vn(0:ndim)

c     Input: the 'subdiagonal' coefficients a(i) in the equations:
      double precision a(0:ndim)
c        Storage location a(n) must exist; don't use a shortened vector.

c     Input: the 'diagonal' coefficients b(i) in the equations:
      double precision b(0:ndim)
c        Locations b(0) and b(n) must exist; don't use a shortened vector.

c     Input: the 'superdiagonal' coefficients c(i) in the equations:
      double precision c(0:ndim)
c        Location c(0) must exist; don't use a shortened vector.

c     Input: the right hand sides d(i) of the equations:
      double precision d(0:ndim)

c     Work: three solutions that satisfy the tridiagonal equations
c     but not the first or final equation:
      double precision u0(0:ndim),u1(0:ndim),u2(0:ndim)

c     Work: storage area to safeguard the values of vector d:
      double precision dstore(0:ndim)

c     Work: work arrays for use by tridia:
      double precision
     & w1(0:ndim),w2(0:ndim),w3(0:ndim),w4(0:ndim),w5(0:ndim)

c     Output: the solution values u(i):
      double precision u(0:ndim)

c     Output: an error code:
      integer ier
c        Code ier is set equal to the sum of:
c     1: if there is only one unknown;
c     2: if the equations are effectively singular;
c     4: if the subroutine crashed: no solution was found.
c        Always check that ier is zero after using tribc.

c     Output: the estimated maximum error in the coefficients of u:
      double precision err
c        This error is based on the round-off error in the coefficients
c     of the equations.  It ignores the possibility of underflow.
c        A small value of err implies that the solution will be accurate
c     as long as underflow plays no role. However, the value of err is
c     typically very pessimistic.

c     Output: the maximum absolute coefficient of u:
      double precision umax

c External variables:

c     Subroutine tribc requires a separate subroutine:
c           tridia(ndim,n,a,b,c,d,w1,w2,w3,w4,w5,u,ier,err,umax)
c     to solve a tridiagonal system of equations of the form
c                     b(0) u(0) + c(0) u(1)   = d(0)  first equation;
c       a(j) u(j-1) + b(j) u(j) + c(j) u(j+1) = d(j)  for j = 0, 1, ..., n-1;
c       a(n) u(n-1) + b(n) u(n)               = d(n)  final equation;
c     for u(0), u(1), ... u(n). Arguments w1 through w5 are work space
c     and ier is an error code which should be less than two after tridia.
c     Subroutine tridia stores its estimated error in err and the maximum
c     absolute value of the u(i) in umax.
c     Subroutine tridia can be found in:
c           ../../../lib/tridia.f
      external tridia

c     Subroutine tribc requires a separate subroutine:
c           satis2(ndim,n,u0,r10,r20,u1,r11,r21,u2,r12,r22,
C                  erru,errr,u,ier,err)
c     to combine three vectors u0, u1, and u2 so that two equations
c     become satisfied. Here r10, r11, and r12 are the residuals in the
c     first equation for u0, u1, and u2 respectively; r20, r21, and r22
c     are the residuals in the second equation; erru is the estimated
c     maximum absolute error in the vectors ui; errr is the estimated
c     maximum absolute error in the residuals, u is the combined
c     solution; ier is an error which should ne nozero after satis2; and
c     err is the estimated maximum absolute error in u.
c     Subroutine satis2 can be found in:
c           ../../../lib/satis2.f
      external satis2

c Local variables:

c     An index to keep track whether tribc has been initialized:
      integer init

c     The maximum norms of vectore v0 and vn:
      double precision v0max,vnmax
c     The L1 norms of vectors v0 and vn:
      double precision v0sum,vnsum

c     The maximum norms of u0, u1, and u2:
      double precision u0max,u1max,u2max
c     The L1 norms of u0, u1, and u2:
      double precision u0sum,u1sum,u2sum
c     The maximum errors in u0, u1, and u2:
      double precision du0,du1,du2

c     The residuals in the first equation for u0,u1, and u2:
      double precision r0u0,r0u1,r0u2
c     The errors in these residuals:
      double precision dr0u0,dr0u1,dr0u2
c     The residuals in the final equation for u0,u1, and u2:
      double precision rnu0,rnu1,rnu2
c     The errors in these residuals:
      double precision drnu0,drnu1,drnu2

c     Temporary storage for b(0), c(0), a(n), and b(n):
      double precision b0stor,c0stor,anstor,bnstor
c     Error code returned by tridia or satis2
      integer ier2

c     Index over the equations:
      integer i

c     The machine epsilon (the relative error in real numbers):
      double precision eps
      save eps

c     Numerical codes for the various problems that tribc may encounter:
      integer onedim,singul,crash
      parameter (onedim=1,singul=2,crash=4)

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

c Data statements:

      data init/0/

c Executable statements:

c     Check the arguments:
      if(ndim.lt.0 .or. n.lt.0 .or. n.gt.ndim)then
         print*,'*** tribc: ',
     &   'Invalid arguments issued to subroutine tribc.'
         print*,'    Correct the program code.'
         stop
      endif

c     Initialize the first time around:
      if(init.eq.0)then

c        Find the machine epsilon as the absolute error in the value 1:
         eps=two
 10      eps=eps/two
         if(one+eps.ne.one)goto 10

c        Mark as initialized:
         init=1
      endif

c     Initialize the error code:
      ier=zero
c     Cannot satisfy two or more equations using a single unknown:
      if(n.le.0)ier=onedim

c     Find the maximum and L1 norms of the vectors v0 and vn:
      v0max=zero
      v0sum=zero
      vnmax=zero
      vnsum=zero
      do 50 i=0,n
         v0max=max(v0max,abs(v0(i)))
         v0sum=v0sum+abs(v0(i))
         vnmax=max(vnmax,abs(vn(i)))
         vnsum=vnsum+abs(vn(i))
 50   continue

c     We will need to play around with the right hand side vector d,
c     b(0), c(0), a(n), and b(n), so we store their values for now:

      do 60 i=0,n
         dstore(i)=d(i)
 60   continue
      b0stor=b(0)
      c0stor=c(0)
      anstor=a(n)
      bnstor=b(n)

c     Find a solution u0 to the intermediate tridiagonal equations:

c     Replace the first equation by a diagonal homogeneous one:
      b(0)=one
      c(0)=zero
      d(0)=zero

c     Replace the final equation by a diagonal homogeneous one:
      a(n)=zero
      b(n)=one
      d(n)=zero

c     Let subroutine tridia solve this system:
      call tridia(ndim,n,a,b,c,d,w1,w2,w3,w4,w5,u0,ier2,du0,u0max)
      if(ier2.ge.crash)then
         ier=ier+crash
         return
      endif
      if(ier2.ge.singul)ier=ier-ier/singul*singul+singul

c     Find the L1 norm of vector u0:
      u0sum=zero
      do 120 i=0,n
         u0sum=u0sum+abs(u0(i))
 120  continue

c     Find the residual in the first equation:
      r0u0=dstore(0)
      do 130 i=0,n
         r0u0=r0u0-v0(i)*u0(i)
 130  continue
c     Estimated maximum error in the residual:
      dr0u0=v0sum*du0+eps*(v0max*u0sum+abs(dstore(0)))

c     Find the residual in the final equation:
      rnu0=dstore(n)
      do 140 i=0,n
         rnu0=rnu0-vn(i)*u0(i)
 140  continue
c     Estimated maximum error in the residual:
      drnu0=vnsum*du0+eps*(vnmax*u0sum+abs(dstore(n)))

c     Find another solution u1 to the intermediate tridiagonal equations:

c     Replace the first equation by a diagonal inhomogeneous one:
      b(0)=one
      c(0)=zero
      d(0)=one

c     Replace the final equation by a diagonal homogeneous one:
      a(n)=zero
      b(n)=one
      d(n)=zero

c     Make the intermediate equations homogeneous:
      do 200 i=1,n-1
         d(i)=zero
 200  continue

c     Let subroutine tridia solve this system:
      call tridia(ndim,n,a,b,c,d,w1,w2,w3,w4,w5,u1,ier2,du1,u1max)
      if(ier2.ge.crash)then
         ier=ier+crash
         return
      endif
      if(ier2.ge.singul)ier=ier-ier/singul*singul+singul

c     Create a solution to the inhomogeneous tridiagonal equations
c     by adding a suitable multiple of the solution u1 to u0:

      if(u0max.gt.zero)then
c        The resulting vector u1:
         do 210 i=0,n
            u1(i)=u0(i)+two*u0max/u1max*u1(i)
 210     continue
c        The estimated error in this vector:
         du1=du0+two*u0max/u1max*du1
      endif

c     Find the maximum and L1 norm of the resulting vector u1:
      u1max=zero
      u1sum=zero
      do 220 i=0,n
         u1max=max(u1max,abs(u1(i)))
         u1sum=u1sum+abs(u1(i))
 220  continue

c     Find the residual in the first equation:
      r0u1=dstore(0)
      do 230 i=0,n
         r0u1=r0u1-v0(i)*u1(i)
 230  continue
c     Estimated maximum error in the residual:
      dr0u1=v0sum*du1+eps*(v0max*u1sum+abs(dstore(0)))

c     Find the residual in the final equation:
      rnu1=dstore(n)
      do 240 i=0,n
         rnu1=rnu1-vn(i)*u1(i)
 240  continue
c     Estimated maximum error in the residual:
      drnu1=vnsum*du1+eps*(vnmax*u1sum+abs(dstore(n)))

c     Find another solution u2 to the intermediate tridiagonal equations:

c     Replace the first equation by a diagonal homogeneous one:
      b(0)=one
      c(0)=zero
      d(0)=zero

c     Replace the final equation by a diagonal inhomogeneous one:
      a(n)=zero
      b(n)=one
      d(n)=one

c     Make the intermediate equations homogeneous:
      do 300 i=1,n-1
         d(i)=zero
 300  continue

c     Let subroutine tridia solve this system:
      call tridia(ndim,n,a,b,c,d,w1,w2,w3,w4,w5,u2,ier2,du2,u2max)
      if(ier2.ge.crash)then
         ier=ier+crash
         return
      endif
      if(ier2.ge.singul)ier=ier-ier/singul*singul+singul

c     Create a solution to the inhomogeneous tridiagonal equations
c     by adding a suitable multiple of the solution u2 to u0:

      if(u0max.gt.zero)then
c        The resulting vector u2:
         do 310 i=0,n
            u2(i)=u0(i)+two*u0max/u2max*u2(i)
 310     continue
c        The estimated error in this vector:
         du2=du0+two*u0max/u2max*du2
      endif

c     Find the maximum and L1 norm of the resulting vector u2:
      u2max=zero
      u2sum=zero
      do 320 i=0,n
         u2max=max(u2max,abs(u2(i)))
         u2sum=u2sum+abs(u2(i))
 320  continue

c     Find the residual in the first equation:
      r0u2=dstore(0)
      do 330 i=0,n
         r0u2=r0u2-v0(i)*u2(i)
 330  continue
c     Estimated maximum error in the residual:
      dr0u2=v0sum*du2+eps*(v0max*u2sum+abs(dstore(0)))

c     Find the residual in the final equation:
      rnu2=dstore(n)
      do 340 i=0,n
         rnu2=rnu2-vn(i)*u2(i)
 340  continue
c     Estimated maximum error in the residual:
      drnu2=vnsum*du2+eps*(vnmax*u2sum+abs(dstore(n)))

c     Let satis2 combine the three solutions to satisfy the first and
c     final equations:

      call satis2(ndim,n,u0,r0u0,rnu0,u1,r0u1,rnu1,u2,r0u2,rnu2,
     & max(du0,du1,du2),max(dr0u0,drnu0,dr0u1,drnu1,dr0u2,drnu2),
     * u,ier2,err)
      if(ier2.ge.crash)then
         ier=ier+crash
         return
      endif
      if(ier2.ge.singul)ier=ier-ier/singul*singul+singul

c     Find the maximum norm of u:
      umax=zero
      do 810 i=0,n
         umax=max(umax,abs(u(i)))
 810  continue

c     Restore the equation system:
      do 820 i=0,n
         d(i)=dstore(i)
 820  continue
      b(0)=b0stor
      c(0)=c0stor
      a(n)=anstor
      b(n)=bnstor

c Exit:

      return
      end