subroutine invtri (z,topbc,botbc,vmc,tdt,kmz,dzq,dzurq,dzlrq,mask,
     *                   istr,iend,aidif)
c
c=======================================================================
c                                                                    ===
c  This routine solves  the vertical diffusion equation implicitly   ===
c  using the method of inverting a tridiagonal matrix as described   ===
c  in Richtmyer and Morton (1967).   It assumes that the variables   ===
c  are  defined at grid points  and  the  top and  bottom boundary   ===
c  conditions are flux conditions. (Adapted from MOM).               ===
c                                                                    ===
c  On input:                                                         ===
c                                                                    ===
c     Z        right hand side terms (real array)                    ===
c     TOPBC    top boundary condition (real array)                   ===
c     BOTBC    bottom boundary condition (real array)                ===
c     VMC      vertical mixing coefficient (cm^2/s; real)            ===
c     TDT      twice of timestep (seconds; real array)               ===
c     KMZ      level indicator (integer array)                       ===
c     DZQ      thickness of vertical boxes (real array)              ===
c     DZURQ    upper box vertical thickness inverse (real array)     ===
c     DZLRQ    lower box vertical thickness inverse (real array)     ===
c     MASK     land/sea mask (real array)                            ===
c     ISTR     starting index (integer)                              ===
c     IEND     ending index (integer)                                ===
c     AIDIF    coefficient for implicit time differencing for        ===
c              vertical diffusion (real)                             ===
c  On output:                                                        ===
c                                                                    ===
c     Z        returned solution (real array)                        ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,istr,iend,k,kmz(imt),kz
      FLOAT
     *      aidif
      FLOAT
     *      a(imt,km),b(imt,km),botbc(imt),c(imt,km),d(imt,km),
     *      dzlrq(imt,km),dzq(imt,km),dzurq(imt,km),e(imt,0:km),
     *      f(imt,0:km),g(imt),mask(imt,km),tdt(km),topbc(imt),
     *      vmc(imt,km),z(imt,km)
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c  Set tridiagonal system.
c
#ifndef barotropic
      do 100 k=2,km
      do 100 i=istr,iend
        a(i,k)=aidif*vmc(i,k-1)*dzurq(i,k)*tdt(k)
        c(i,k)=aidif*vmc(i,k)*dzlrq(i,k)*tdt(k)
        b(i,k)=c1+a(i,k)+c(i,k)
        d(i,k)=z(i,k)
        e(i,k-1)=c0
        f(i,k-1)=c0
 100  continue
#endif
c
c  Set surface boundary conditions.
c
      k=1
      do 110 i=istr,iend
        a(i,k)=aidif*tdt(k)/dzq(i,k)
        c(i,k)=aidif*vmc(i,k)*dzlrq(i,k)*tdt(k)
        b(i,k)=c1+c(i,k)
        d(i,k)=z(i,k)
        e(i,k-1)=c0
        f(i,k-1)=c0
 110  continue
c
c  Set bottom boundary conditions.
c
      do 120 i=istr,iend
        kz=kmz(i)
        if(kz.ne.0) then
          b(i,kz)=c1+a(i,kz)
          c(i,kz)=aidif*tdt(kz)/dzq(i,kz)
          e(i,kz)=c0
          f(i,kz)=-botbc(i)
        endif
 120  continue
c
c-----------------------------------------------------------------------
c  Now invert.
c-----------------------------------------------------------------------
c
      do 130 k=km,1,-1
      do 130 i=istr,iend
        if(k.le.kmz(i)) then
          g(i)=c1/(b(i,k)-c(i,k)*e(i,k))
          e(i,k-1)=a(i,k)*g(i)
          f(i,k-1)=(d(i,k)+c(i,k)*f(i,k))*g(i)
        endif
 130  continue
c
c Put surface boundary conditions.
c
      do 140 i=istr,iend
        z(i,1)=(e(i,0)*topbc(i)+f(i,0))*mask(i,1)
 140  continue
c
#ifndef barotropic
      do 150 k=2,km
      do 150 i=istr,iend
        z(i,k)=(e(i,k-1)*z(i,k-1)+f(i,k-1))*mask(i,k)
 150  continue
#endif
      return
      end