#if !defined usrdiagnostic | !defined nesttime
subroutine tracer(j)
#else
subroutine tracer(j,dtsr,dtrm)
#endif
c
c=======================================================================
c ===
c TRACER computes, for one row, the NT tracers, where: ===
c ===
c J = the row number ===
c ===
c=======================================================================
c
c-----------------------------------------------------------------------
c Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <fullwd.h>
#include <scalar.h>
#include <onedim.h>
#include <fields.h>
#include <workspa.h>
#include <cvmix.h>
#include <cvbc.h>
#include <cdiag.h>
#include <hybrid.h>
#include <vertslabs.h>
#include <filtdat.h>
#include <options.h>
#ifdef bndy_rlx
# include <bndyrlx.h>
#endif
#ifdef linear_physics
# include <rhomean.h>
#endif
#if defined ext_tide & defined advtide
# include <tidesp.h>
#endif
#if defined codunlim | defined codlim
# include <cbiopnh.h>
#endif
#if defined surfpress & defined freesurf
# include <fieldsbar.h>
# include <sinfo.h>
# include <vertical.h>
#endif
c
c-----------------------------------------------------------------------
c Define local data.
c-----------------------------------------------------------------------
c
integer i,iu,j,k,kz,m,sofar
FLOAT
* boxvol,boxfac,fx,fxa,fxb,ssst,sstf
#ifdef explicitvmix
# ifndef barotropic
integer ks,n
FLOAT
* factor
# endif
# else
FLOAT
* rc2dt,twodt(km)
#endif
#if defined linear_physics & !defined bioMcGillic & !defined bioFasham & !defined bioAnder & !defined bioDuse
FLOAT
* der1,derkm,f,val
FLOAT
* d2(mprof),wk(mprof),zvec(imt)
c
parameter (der1=c1e30,derkm=c1e30)
#endif
#if defined usrdiagnostic & defined nesttime
FLOAT
* ddtrm(2),ddtsr(2),dtrm(2),dtsr(2)
# ifdef iforttime
& ,dum
FLOAT
& dtime
# endif
#endif
#if defined surfpress & defined freesurf
FLOAT
* wfree(imt,kmp1)
#endif
c
c=======================================================================
c Begin executable code.
c=======================================================================
c
#ifdef frozentrc
c=======================================================================
c Omit timesteping of tracers during initialization.
c (density is maintained constant for NTDGN timesteps).
c=======================================================================
c
if(itt.le.ntdgn) then
do 100 m=1,nt
do 100 k=1,km
do 100 i=1,imt
ta(i,k,m)=t(i,k,m)
100 continue
return
endif
c
#endif
#if defined usrdiagnostic & defined nesttime
c-----------------------------------------------------------------------
c Initialize time counters.
c-----------------------------------------------------------------------
c
do m = 1,2
dtrm(m) = c0
dtsr(m) = c0
end do
c
#endif
c=======================================================================
c Begin introductory section, preparing various ======================
c arrays for the computation of the tracers ======================
c=======================================================================
c
c-----------------------------------------------------------------------
c Compute the advective coefficients FUW at west face of T grid box
c and FVN at the north face of T grid box.
c-----------------------------------------------------------------------
c
fxa=cstr(j)*dytr(j)
fxb=fxa*cs(j)
do 110 k=1,km
do 110 i=2,imt
iu = min(i,imu)
#if (!defined ext_tide | !defined advtide) & (!defined surfpress | !defined freesurf)
fuw(i,k)=(u (i-1,k)*dyu(j )*dzvqz(i-1,k,0)+
* um(i-1,k)*dyu(j-1)*dzvqz(i-1,k,1))*fxa
fvn(i,k)=(v (iu ,k)*dxuq(i ,k)*dzvqz(iu ,k,0)+
* v (i-1,k)*dxuq(i-1,k)*dzvqz(i-1,k,0))*fxb*dxt4rq(i,k)
110 continue
#elif defined surfpress & defined freesurf
fuw(i,k)=( (ucl(i-1,k)*dzvqz(i-1,k,0)+ dzvqzb(i-1,k,0)*
& (ubaro(i-1,j)+ubarob(i-1,j))*p5)*dyu(j ) +
* (uclm(i-1,k)*dzvqz(i-1,k,1)+ dzvqzb(i-1,k,1)*
& (ubaro(i-1,j-1)+ubarob(i-1,j-1))*p5)*dyu(j-1) )*fxa
fvn(i,k)=( (vcl(iu,k)*dzvqz(iu ,k,0)+ dzvqzb(i ,k,0)*
& (vbaro(iu,j)+vbarob(iu,j))*p5)*dxuq(i ,k) +
* (vcl(i-1,k)*dzvqz(i-1,k,0)+ dzvqzb(i-1,k,0)*
& (vbaro(i-1,j)+vbarob(i-1,j))*p5)*dxuq(i-1,k) )*fxb*
* dxt4rq(i,k)
110 continue
#else
fuw(i,k)=(u (i-1,k)*dyu(j )*dzvqz(i-1,k,0)*ustretch(i-1,j)+
* um(i-1,k)*dyu(j-1)*dzvqz(i-1,k,1)*ustretch(i-1,j-1))*fxa
fvn(i,k)=(v (iu ,k)*dxuq(i ,k)*dzvqz(iu ,k,0)*ustretch(iu,j)+
* v (i-1,k)*dxuq(i-1,k)*dzvqz(i-1,k,0)*ustretch(i-1,j))*
* fxb*dxt4rq(i,k)
110 continue
c
c Add tidal contributions.
c
do 113 k=1,km
do 113 i=2,imt
iu = min(i,imu)
fuwtd(i,k)=
* (utide (i-1,k)*dyu(j )*dzvqz(i-1,k,0)*ustretch(i-1,j)+
* utidem(i-1,k)*dyu(j-1)*dzvqz(i-1,k,1)*ustretch(i-1,j-1))*fxa
fvntd(i,k)=
* (vtide (iu ,k)*dxuq(i ,k)*dzvqz(iu ,k,0)*ustretch(iu,j)+
* vtide (i-1,k)*dxuq(i-1,k)*dzvqz(i-1,k,0)*ustretch(i-1,j))
* *fxb*dxt4rq(i,k)
# ifdef advtide0
fuwtd(i,k)=fuwtd(i,k)*sadv
fvntd(i,k)=fvntd(i,k)*sadv
# endif
113 continue
#endif
#ifdef cyclic
c
c Set Cyclic boundary conditions.
c
do 117 k=1,km
fuw(imt,k)=fuw(2,k)
fuw(1,k)=fuw(imtm1,k)
#if defined ext_tide & defined advtide
fuwtd(imt,k)=fuwtd(2,k)
fuwtd(1,k)=fuwtd(imtm1,k)
#endif
117 continue
#endif
#if defined surfpress | ( defined ext_tide & defined advtide )
c-----------------------------------------------------------------------
c Compute "omega" vertical velocity in T columns.
c-----------------------------------------------------------------------
c
c Set W at the bottom. Kinematic boundary condition
c
do 120 i=2,imtm1
kz=km
w(i,kz+1)=c0
# if defined freesurf
wfree(i,kz+1)=c0
# endif
120 continue
c
c Compute change of W between levels.
c
# ifdef freesurf
fx=c1/c2dtps
c
# endif
do 140 k=1,km
do 140 i=2,imtm1
# if defined ext_tide & defined advtide
w(i,k)=c2*((fuw(i+1,k)+fuwtd(i+1,k)
* -fuw (i,k)-fuwtd(i,k))*dxt4rq(i,k)
* +(fvn(i ,k)+fvntd(i ,k)
* -fvst(i,k)-fvstdt(i,k)) )
# else
w(i,k)=c2*((fuw(i+1,k)-fuw (i,k))*dxt4rq(i,k)
* + (fvn(i ,k)-fvst(i,k)) )
# if defined freesurf
wfree(i,k) = detat(i,j)*fx*dzt(i,j,k)*htav(i,j)
w(i,k) = w(i,k) + wfree(i,k)
# endif
# endif
140 continue
c
c Integrate upwards from the bottom
c
do 150 k=km,1,-1
do 150 i=2,imtm1
w(i,k)=-w(i,k)+w(i,k+1)
# if defined freesurf
wfree(i,k) = wfree(i,k) + wfree(i,k+1)
# endif
150 continue
#else
c
c-----------------------------------------------------------------------
c Compute "omega" vertical velocity in T columns.
c-----------------------------------------------------------------------
c
c Set "omega" vertical velocity at the surface to zero (rigid-lid).
c
do 120 i=2,imtm1
w(i,1)=c0
120 continue
c
c Compute change of W between levels.
c
do 140 k=1,km
do 140 i=2,imtm1
w(i,k+1)=c2*((fuw(i+1,k)-fuw (i,k))*dxt4rq(i,k)
* +(fvn(i ,k)-fvst(i,k)) )
140 continue
c
c Integrate downward from the surface.
c
do 150 k=1,km
do 150 i=2,imtm1
w(i,k+1)=w(i,k)+w(i,k+1)
150 continue
#endif
c
c-----------------------------------------------------------------------
c Compute standard vertical velocity from the omega vertical
c velocity.
c-----------------------------------------------------------------------
c
if((diagts.or.wrtts).and.eots) then
do 170 k=1,km
wvelt(1,k)=c0
wvelt(imt,k)=c0
do 170 i=2,imtm1
wvelt(i,k)=p5*(w(i,k+1)+w(i,k))-
* (u(i-1,k)+u(i,k)+um(i-1,k)+um(i,k))
* *cstr(j)*dxt4rq(i,k)*p5*
* ( (vdepth(i ,k,jrn)-
* vdepth(i-1,k,jrn) )+
* (vdepth(i ,k,jrs) -
* vdepth(i-1,k,jrs) ) )-
* (v(i-1,k)+v(i,k)+vm(i-1,k)+vm(i,k))
* *dyt4r(j)*p5*
* ( (vdepth(i ,k,jrn)-
* vdepth(i ,k,jrs)) +
* (vdepth(i-1,k,jrn)-
* vdepth(i-1,k,jrs)) )
#if defined surfpress & defined freesurf
* + p5*(wfree(i,k+1)+wfree(i,k))
#endif
170 continue
#ifdef cyclic
do 175 k=1,km
wvelt(1,k)=wvelt(imtm1,k)
wvelt(imt,k)=wvelt(2,k)
175 continue
#endif
endif
c
c=======================================================================
c End introductory section ===========================================
c=======================================================================
c
c=======================================================================
c Begin computation of the tracers. =====================
c the new values "TA", will first be loaded with =====================
c the time rate of change, and then updated. =====================
c=======================================================================
c
#ifdef bioDuse
c
c-----------------------------------------------------------------------
c Calculate primary productivity for use in biosource.
c-----------------------------------------------------------------------
c
call priprod
#endif
do 350 m=1,nt
c
c-----------------------------------------------------------------------
c Calculate quantities for the computation of vertical diffusion.
c-----------------------------------------------------------------------
c
#ifndef barotropic
do 180 k=1,kmm1
do 180 i=2,imtm1
vtf(i,k,m)=vdc(i,k)*(tb(i,k,m)-tb(i,k+1,m))/dzzqz(i,k+1,0)
180 continue
#endif
c
c Set the K=0 elements of VTF to reflect surface tracer flux and set
c the K=KZ elements of VTF to reflect insulation condition.
c
do 190 i=2,imtm1
kz=kmt(i)
vtf(i,0,m)=stf(i,m)
vtf(i,kz,m)=btf(i,m)
190 continue
#if !defined notadvt | defined linear_physics
c
c-----------------------------------------------------------------------
c Compute total advection of tracers.
c-----------------------------------------------------------------------
# if !defined notadvt & !defined linear_physics
c
c Compute zonal advection of tracer.
c
do 200 k=1,km
do 200 i=2,imtm1
# if !defined ext_tide | !defined advtide
UTx(i,k)=(fuw(i ,k)*(t(i ,k,m)+t(i-1,k,m))-
* fuw(i+1,k)*(t(i+1,k,m)+t(i ,k,m)))
* *dxt4rq(i,k)/dzqz(i,k,0)
# else
UTx(i,k)=((fuw(i ,k)+fuwtd(i ,k))*(t(i ,k,m)+t(i-1,k,m))-
* (fuw(i+1,k)+fuwtd(i+1,k))*(t(i+1,k,m)+t(i ,k,m)))
* *dxt4rq(i,k)/(dzqz(i,k,0)*tstretch(i,j))
# endif
200 continue
c
c Compute meridional advection of tracer.
c
do 210 k=1,km
do 210 i=2,imtm1
# if !defined ext_tide | !defined advtide
VTy(i,k)=( fvst(i,k)*(t (i,k,m)+tm(i,k,m))
* -fvn (i,k)*(tp(i,k,m)+t (i,k,m)))
* /dzqz(i,k,0)
# else
VTy(i,k)=( (fvst(i,k)+fvstdt(i,k))*(t (i,k,m)+tm(i,k,m))
* -(fvn (i,k)+fvntd (i,k))*(tp(i,k,m)+t (i,k,m)))
* /(dzqz(i,k,0)*tstretch(i,j))
# endif
210 continue
# endif
c
c Compute vertical advection of tracer.
c
# ifndef linear_physics
# ifndef barotropic
do 220 k=2,km
do 220 i=2,imtm1
tempb(i,k)=w(i,k)*(t(i,k-1,m)+t(i,k,m))
220 continue
# endif
do 230 i=2,imtm1
tempb(i, 1)=w(i, 1)*t(i, 1,m)
tempb(i,kmp1)=w(i,kmp1)*t(i,km,m)
230 continue
# elif !defined bioMcGillic & !defined bioFasham & !defined bioAnder & !defined bioDuse
if (iflag(2).eq.1) then
call spline(tinit(1,nt+1),tinit(1,m),nprof,der1,derkm,d2,wk)
endif
do 220 i = 2, imtm1
zvec(i) = c0
if (iflag(2).eq.1) then
call splint(tinit(1,nt+1),tinit(1,m),d2,nprof,zvec(i),val,f)
else
call lintrp(nprof,tinit(1,nt+1),tinit(1,m),1,zvec(i),val)
end if
tempb(i,1) = w(i,1)*c2*val
220 continue
do 230 k = 2, kmp1
do 230 i = 2, imtm1
zvec(i) = zvec(i) + dzqz(i,k-1,0)
if (iflag(2).eq.1) then
call splint(tinit(1,nt+1),tinit(1,m),d2,nprof,zvec(i),val,f)
else
call lintrp(nprof,tinit(1,nt+1),tinit(1,m),1,zvec(i),val)
end if
tempb(i,k) = w(i,k)*c2*val
230 continue
# endif
c
do 240 k=1,km
do 240 i=2,imtm1
WTz(i,k)=(tempb(i,k+1)-tempb(i,k))*dz2rqz(i,k,0)
240 continue
#endif
c
c-----------------------------------------------------------------------
c Compute horizontal diffusion of tracers (evaluate at TAU-1 timestep).
c-----------------------------------------------------------------------
c
if((mixtrc.eq.2).or.(mixtrc.eq.3)) then
if(mixtrc.eq.2) then
call lapt_depth(j,tbm(1,1,m),tb(1,1,m),tbp(1,1,m),fmm,fm,
* fmp,Txx,Tyy)
else
call lapt_lev(j,km,tbm(1,1,m),tb(1,1,m),tbp(1,1,m),fmm,fm,
* fmp,Txx,Tyy)
endif
do 250 k=1,km
do 250 i=2,imtm1
Txx(i,k)=ah*Txx(i,k)
Tyy(i,k)=ah*Tyy(i,k)
250 continue
endif
c
c-----------------------------------------------------------------------
c Compute vertical diffusion of tracers.
c-----------------------------------------------------------------------
c
do 260 k=1,km
do 260 i=2,imtm1
Tzz(i,k)=(vtf(i,k-1,m)-vtf(i,k,m))/dzqz(i,k,0)
260 continue
c
c-----------------------------------------------------------------------
c Compute source term, Tsrc.
c-----------------------------------------------------------------------
c
c Initialize tracer source terms.
c
#if defined usrdiagnostic & defined nesttime
# ifndef iforttime
call dtime (ddtrm)
# else
dum = dtime(ddtrm)
# endif
dtrm(1) = dtrm(1) + ddtrm(1)
dtrm(2) = dtrm(2) + ddtrm(2)
#endif
do 262 k=1,km
do 262 i=1,imt
Tsrc(i,k)=c0
262 continue
#if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
if(m.gt.2) then
call biosource(j,m)
endif
#endif
#ifdef pttrcsrc
call tsource(j,m)
#endif
#if defined usrdiagnostic & defined nesttime
# ifndef iforttime
call dtime (ddtsr)
# else
dum = dtime(ddtsr)
# endif
dtsr(1) = dtsr(1) + ddtsr(1)
dtsr(2) = dtsr(2) + ddtsr(2)
#endif
c
#ifdef bndy_rlx
c-----------------------------------------------------------------------
c Compute boundary relaxation source terms.
c-----------------------------------------------------------------------
c
do 265 k=1,km
do 265 i=2,imtm1
if(itt.eq.1)then
t_0(i,j,k,m)=tb(i,k,m)
endif
Tbrlx(i,k)=tfacbrlx(i,j)*(t_0(i,j,k,m)-tb(i,k,m))
265 continue
c
#endif
c-----------------------------------------------------------------------
c Calculate the new tracer quantities allowing for implicit treatment
c of vertical diffusion. Reset land points to zero.
c-----------------------------------------------------------------------
c
do 270 k=1,km
do 270 i=2,imtm1
#if defined bioadjvert | defined bioadjloc
if(m.le.2)then
ta(i,k,m)=t(i,k,m)
else
#endif
#if !defined bndy_rlx | !defined imp_bnd_rlx
ta(i,k,m)=fm(i,k)*(tb(i,k,m)+c2dtts*(
# else
ta(i,k,m)=fm(i,k)*(tb(i,k,m)+
* c2dtts/(c1+p5*c2dtts*tfacbrlx(i,j))*(
#endif
#if !defined bioadjvert & !defined bioadjloc
# if !defined notadvt & !defined linear_physics
* UTx(i,k)+VTy(i,k)
# endif
#endif
#ifndef bioadjloc
#if !defined notadvt | defined linear_physics
* +WTz(i,k)
#endif
#endif
#if !defined bioadjvert & !defined bioadjloc
* +Txx(i,k)+Tyy(i,k)
#endif
#ifdef bioadjloc
c note that vertical diffusion has *not* (by design) been shut off
#endif
#ifdef explicitvmix
* +Tzz(i,k)
#else
* +Tzz(i,k)*(c1-aidif)
#endif
#ifdef bndy_rlx
* +Tbrlx(i,k)
#endif
* +Tsrc(i,k)))
#if defined bioadjvert | defined bioadjloc
endif
#endif
270 continue
#ifndef explicitvmix
c
c-----------------------------------------------------------------------
c Solve vertical diffusion implicitly.
c-----------------------------------------------------------------------
c
c Store terms to compute implicit vertical mixing on diagnostic
c timesteps.
c
if(diagts.and.eots) then
do 280 k=1,km
do 280 i=2,imtm1
tempa(i,k)=ta(i,k,m)
280 continue
endif
do 290 k=1,km
twodt(k)=c2dtts
290 continue
# if (!defined codunlim & !defined codlim) | defined codvmix
c
c Add in the implicit vertical diffusion.
c
# else
c Add in the implicit vertical diffusion (except for cod).
c
if (m.ne.icod) then
# endif
call invtri (ta(1,1,m),stf(1,m),btf(1,m),vdc,twodt,kmt,
* dzqz(1,1,0),dzturq,dztlrq,fm,2,imtm1,aidif)
# if (defined codunlim | defined codlim) & !defined codvmix
end if
# endif
c
c Compute residual implicit vertical diffusion.
c
if(diagts.and.eots) then
do 300 k=1,km
rc2dt=c1/twodt(k)
do 300 i=2,imtm1
tempa(i,k)=rc2dt*(ta(i,k,m)-tempa(i,k))
300 continue
endif
#endif
c
c-----------------------------------------------------------------------
c Do analysis of tracers on diagnostic timesteps.
c-----------------------------------------------------------------------
c
if(diagts.and.eots) then
fx=cst(j)*dyt(j)
c
c Compute buoyancy: total energy exchange between potential and
c kinetic.
c
if(m.eq.1) then
fxa=p5*grav
do 320 i=2,imtm1
fxb=fx*dxt(i)
kz=kmt(i)
if(kz.ge.2) then
do 310 k=2,kz
buoy(i,k)=-wvelt(i,k-1)*(rhos(i,k-1)+rhos(i,k))*fxa
buoyint(k)=buoyint(k)+buoy(i,k)*fxb*dzqz(i,k-1,0)
310 continue
endif
320 continue
endif
c
c Compute term balance components for time rate of change of tracers
c per unit volume.
c
do 340 i=2,imtm1
kz=kmt(i)
if(kz.ne.0) then
fxa=fx*dxt(i)
sstf=fxa*stf(i,m)
ssst=fxa*t(i,1,m)
asst(m)=asst(m)+ssst
stflx(m)=stflx(m)+sstf
do 330 k=1,kz
boxvol=fxa*dzqz(i,k,0)
trmbtv(i,k,2,m)=fm(i,k)*UTx(i,k)
termbt(k,2,m)=termbt(k,2,m)+UTx(i,k)*boxvol
trmbtv(i,k,3,m)=fm(i,k)*VTy(i,k)
termbt(k,3,m)=termbt(k,3,m)+VTy(i,k)*boxvol
trmbtv(i,k,4,m)=fm(i,k)*WTz(i,k)
termbt(k,4,m)=termbt(k,4,m)+WTz(i,k)*boxvol
if((mixtrc.eq.2).or.(mixtrc.eq.3)) then
trmbtv(i,k,5,m)=fm(i,k)*Txx(i,k)
termbt(k,5,m)=termbt(k,5,m)+Txx(i,k)*boxvol
trmbtv(i,k,6,m)=fm(i,k)*Tyy(i,k)
termbt(k,6,m)=termbt(k,6,m)+Tyy(i,k)*boxvol
endif
#ifdef explicitvmix
trmbtv(i,k,7,m)=fm(i,k)*Tzz(i,k)
termbt(k,7,m)=termbt(k,7,m)+Tzz(i,k)*boxvol
#else
trmbtv(i,k,7,m)=fm(i,k)*(Tzz(i,k)*(c1-aidif)+
* tempa(i,k))
termbt(k,7,m)=termbt(k,7,m)+(Tzz(i,k)*(c1-aidif)+
* tempa(i,k))*boxvol
#endif
trmbtv(i,k,8,m)=fm(i,k)*Tsrc(i,k)
termbt(k,8,m)=termbt(k,8,m)+Tsrc(i,k)*boxvol
sofar=8
#ifdef bndy_rlx
sofar=sofar+1
trmbtv(i,k,sofar,m)=fm(i,k)*(Tbrlx(i,k)+
* (ta(i,k,m)-tb(i,k,m))*p5*c2dtts*tfacbrlx(i,j))
termbt(k,sofar,m)=termbt(k,sofar,m)+boxvol*
* trmbtv(i,k,sofar,m)
#endif
atwv(m)=atwv(m)+t(i,k,m)
330 continue
endif
340 continue
endif
c
350 continue
#if defined explicitvmix & !defined barotropic
c
c-----------------------------------------------------------------------
c Convectively adjust water column if gravitationally unstable.
c-----------------------------------------------------------------------
c
c Perform NCON passes through convection loop:
c
c KS=1: compare level 1 to 2; 3 to 4; etc. and adjust if necessary.
c KS=2: compare level 2 to 3; 4 to 5; etc. and adjust if necessary.
c
if(ncon.gt.0) then
do 390 n=1,ncon
do 390 ks=1,2
c
c 1st, find density for entire slab for stability determination.
c
call state(ta(1,1,1),ta(1,1,2),tdepth(1,1,jrs),tempa)
c
c 2nd, for each tracer, mix adjoining levels if unstable.
c
do 380 m=1,nt
do 370 k=ks,kmm1,2
do 370 i=2,imtm1
if(kmt(i).gt.0) then
if(tempa(i,k).gt.tempa(i,k+1)) then
factor=dzqz(i,k,0)*dzz2rqz(i,k,0)
ta(i,k,m)=factor*ta(i,k,m)+(1.0-factor)*ta(i,k+1,m)
ta(i,k+1,m)=ta(i,k,m)
endif
endif
370 continue
380 continue
390 continue
endif
#endif
c
c-----------------------------------------------------------------------
c Integrate total changes in tracers and squared tracer on diagnostic
c timesteps over the regional volume.
c-----------------------------------------------------------------------
c
if(diagts.and.eots) then
fx=cst(j)*dyt(j)/c2dtts
do 400 m=1,nt
do 400 k=1,km
do 400 i=2,imtm1
boxfac=fm(i,k)*fx*dxtq(i,k)*dzqz(i,k,0)
termbt(k,1,m)=termbt(k,1,m)+(ta(i,k,m)-tb(i,k,m))*boxfac
tvar(k,m)=tvar(k,m)+(ta(i,k,m)**2-tb(i,k,m)**2)*boxfac
400 continue
endif
c
c-----------------------------------------------------------------------
c Accumulate integrated absolute changes in tracers over the regional
c ocean volume.
c-----------------------------------------------------------------------
c
if(mod(itt,ntsi).eq.0) then
fx=cst(j)*dyt(j)/dtts
do 430 m=1,nt
do 430 k=1,km
do 410 i=1,imt
boxfac=fm(i,k)*fx*dxtq(i,k)*dzqz(i,k,0)
tempb(i,k)=abs(ta(i,k,m)-t(i,k,m))*boxfac
410 continue
do 420 i=1,imt
dtabs(k,m,j)=dtabs(k,m,j)+tempb(i,k)
420 continue
430 continue
endif
c
c=======================================================================
c End computation of the tracers =====================================
c=======================================================================
c
c-----------------------------------------------------------------------
c Transfer quantities computed to the north of the present row
c to be defined to the south in the computation of the next row
c-----------------------------------------------------------------------
c
fx=cst(j)*dyt(j)*cstr(j+1)*dytr(j+1)
do 440 k=1,km
do 440 i=1,imt
fvst(i,k)=fvn(i,k)*fx
#if defined ext_tide & defined advtide
fvstdt(i,k)=fvntd(i,k)*fx
#endif
440 continue
#ifdef cyclic
c
c Set Cyclic boundary conditions on newly computed tracers.
c
do 450 m=1,nt
do 450 k=1,km
ta(1 ,k,m)=ta(imtm1,k,m)
ta(imt,k,m)=ta(2 ,k,m)
450 continue
#endif
#if defined usrdiagnostic & defined nesttime
# ifndef iforttime
call dtime (ddtrm)
# else
dum = dtime(ddtrm)
# endif
dtrm(1) = dtrm(1) + ddtrm(1)
dtrm(2) = dtrm(2) + ddtrm(2)
#endif
return
end