!
      subroutine trsolv(a,b,c,f,x,nlev,k1,k2,nlon,i1,i2,iprint)
      implicit none
!
! Tri-diagonal solver.
!   a(k,i)*x(k-1,i) + b(k,i)*x(k,i) + c(k,i)*x(k+1,i) = f(k,i)
! This is called from sub minor (comp.F)
! 3/05 btf: Adapted for glbmean from timegcm1.2
!
! Args:
      integer,intent(in) :: nlev,k1,k2,nlon,i1,i2,iprint
      real,dimension(nlev,nlon),intent(in) ::
     |  a, ! input coefficients
     |  b, ! input coefficients
     |  c, ! input coefficients
     |  f  ! input RHS
      real,dimension(nlev,nlon),intent(out) :: x  ! output
!
! Local:
      integer :: k,kk,i,nlevs
      real,dimension(nlev,nlon) ::
     |  w1,w2,w3  ! work arrays
!
      nlevs = k2-k1+1
!
! Lower boundary (W(K1)=B(K1):
      do i=i1,i2
        w1(k1,i) = b(k1,i)
      enddo
!
! Set up work arrays:
      do i=i1,i2
        do k=k1+1,k2
!
! W(KF+K-1)=C(K-1)/W(K-1):
          w2(k-1,i) = c(k-1,i) / w1(k-1,i)
!
! W(K)=A(K)*W(KF+K-1)
          w1(k,i) = a(k,i) * w2(k-1,i)
!
! W(K)=B(K)-W(K)
          w1(k,i) = b(k,i) - w1(k,i)
        enddo ! k=k1+1,k2

        if (iprint > 0) then
          write(6,"('trsolv: i=',i3,' nlevs=',i3,' w1=',/,(6e12.4))")
     |      i,nlevs,w1(k1:k2,i)
          write(6,"('trsolv: i=',i3,' nlevs=',i3,' w2=',/,(6e12.4))")
     |      i,nlevs,w2(k1:k2,i)
        endif
      enddo ! i=lon0,lon1
!
! Lower boundary (W(2*KF+K1)=F(K1)/W(K1)):
      do i=i1,i2
        w3(k1,i) = f(k1,i) / w1(k1,i)
      enddo ! i=i1,i2
!
      do i=i1,i2
        do k=k1+1,k2
!
! W(2*KF+K)=A(K)*W(2*KF+K-1)
          w3(k,i) = a(k,i) * w3(k-1,i)
!
! W(2*KF+K)=F(K)-W(2*KF+K)
          w3(k,i) = f(k,i) - w3(k,i)         
!
! W(2*KF+K)=W(2*KF+K)/W(K)
          w3(k,i) = w3(k,i) / w1(k,i)
        enddo ! k=k1+1,k2
!
        if (iprint > 0)
     |    write(6,"('trsolv: i=',i3,' nlevs=',i3,' w3=',/,(6e12.4))")
     |      i,nlevs,w3(k1:k2,i)
      enddo ! i=i1,i2
!
! Upper boundary (X(K2)=W(2*KF+K2)):
      do i=i1,i2
        x(k2,i) = w3(k2,i)       
      enddo ! i=i1,i2
!
! Back substitution:
      do i=i1,i2
        do kk=k1+1,k2          
          k = k1+k2-kk ! k2-1,k1,-1
!
! X(K)=W(KF+K)*X(K+1)
          x(k,i) = w2(k,i) * x(k+1,i)
! X(K)=W(2*KF+K)-X(K)
          x(k,i) = w3(k,i) - x(k,i)
        enddo ! k=k1+1,k2
        if (iprint > 0)
     |    write(6,"('trsolv: i=',i3,' nlevs=',i3,' x=',/,(6e12.4))")
     |      i,nlevs,x(k1:k2,i)
      enddo ! i=i1,i2
      end subroutine trsolv
