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