subroutine trisolve ( a, b, c, x, r, ndim)
c
c    abcd:tritest.trisolve <------------- inverts the tridiagonal matrix
c                                                            march, 2002
c
c    written by: A Busy Code Developer
c    modified 1: 
c
c  PURPOSE:  This routine solves the tridiagonal set of equations:
c       _                                 _   _    _       _    _
c      | b(1)   c(1)                       | | x(1) |     | r(1) |
c      | a(2)   b(2)   c(2)                | | x(2) |     | r(2) |
c      |        a(3)   b(3)   c(3)         | | x(3) |  =  | r(3) |
c      |              . . .                | |  ... |     |  ... |
c      |              a(n-1) b(n-1) c(n-1) | |x(n-1)|     |r(n-1)|
c      |_                     a(n)   b(n) _| |_x(n)_|     |_r(n)_|
c
c
c  INPUT VARIABLES:
c    a           sub-diagonal of the tridiagonal matrix
c    b               diagonal of the tridiagonal matrix
c    c         super-diagonal of the tridiagonal matrix
c    r         source term
c    ndim      dimension of a, b, c, r, and x
c
c  OUTPUT VARIABLES:
c    x         solution to tridiagional system of equations
c
c  EXTERNALS: [none]
c
c-----------------------------------------------------------------------
c
       implicit      none
c
       integer       nmax
       parameter     ( nmax = 100 )
c
       integer       ndim    , n
c
       real*8        a       (ndim), b       (ndim), c       (ndim)
     1             , r       (ndim), x       (ndim)
       real*8        beta    (nmax), rho     (nmax)
c
c-----------------------------------------------------------------------
c
       if (ndim .gt. nmax) then
         write ( 6, 2010 ) nmax, ndim
         stop
       endif
       if (b(1) .eq. 0.0d0) then
         write ( 6, 2020 )
         stop
       endif
c
c      Gauss-eliminate the sub-diagonal.
c
       beta(1) = b(1)
       rho (1) = r(1)
       do 10 n=2,ndim
         beta(n) = b(n) - a(n) * c  (n-1) / beta(n-1)
         rho (n) = r(n) - a(n) * rho(n-1) / beta(n-1)
         if (beta(n) .eq. 0.0d0) then
           write ( 6, 2030 ) n
           stop
         endif
10     continue
c
c      Perform the back substitutions
c
       x(ndim) = rho(ndim) / beta(ndim)
       do 20 n=ndim-1,1,-1
         x(n) = ( rho(n) - c(n) * x(n+1) ) / beta(n)
20     continue
c
c-----------------------------------------------------------------------
c----------------------- Write format statements -----------------------
c-----------------------------------------------------------------------
c
2010   format('TRISOLVE: Recompile with either ndim < ',i4
     1       ,', or nmax > ',i4,'.')
2020   format('TRISOLVE: Diagonal element of reduced matrix (n=   1) '
     1       ,'is zero: no unique ',/
     2       ,'TRISOLVE: solution exists.  Abort!')
2030   format('TRISOLVE: Diagonal element of reduced matrix (n=',i4
     1       ,') is zero: no unique ',/
     2       ,'TRISOLVE: solution exists.  Abort!')
c
       return
       end
c