subroutine tridia(jdim,jmax,diags,b)

c     Subroutine to solve a tridiagonal system of the form
c       diags(1,j)*u(j-1) + diags(2,j)*u(j) + diags(3,j)*u(j+1) = b(j)
c     for j = 0, 1, 2, 3, ... jmax (where diags(1,0) = diags(3,jmax) = 0.)

c     The vector b is altered by tridia and contains the values of the
c     solution vector u after tridia.

c     The matrix diags is also altered:
c     therefore, to solve the same system, but with different right hand side
c     b, again, restore the matrix diags.

c     no pivoting is performed.

c parameters

c     declared and actual dimension
      integer jdim,jmax

c     coefficients and right hand sides of the equations
      real diags(3,0:jdim),b(0:jdim)

c local variables

c     equation index
      integer j

c     amount of the previous equation to substract
      real fac

c exeutable statements

c     reduce the matrix to bidiagonal form:

      do 10 j=1,jmax

c        amount of the previous equation to substract to eliminate diags(1,j)
         fac=diags(1,j)/diags(2,j-1)

c        substract this amount
         diags(2,j)=diags(2,j)-fac*diags(3,j-1)
         b(j)=b(j)-fac*b(j-1)

 10   continue

c     solve the bidiagonal system and store the solution in b(j):

      b(jmax)=b(jmax)/diags(2,jmax)
      do 20 j=jmax-1,0,-1
         b(j)=(b(j)-diags(3,j)*b(j+1))/diags(2,j)
 20   continue

      return
      end