subroutine step
c
c=======================================================================
c ===
c STEP is called once per timestep. It initializes various ===
c quantities, bootstraps the basic row by row computation ===
c of prognostic variables, manages the I/O for the latter, ===
c and performs various analysis procedures on the ===
c progressing solution. ===
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 <fieldsbar.h>
#include <workspa.h>
#include <bndata.h>
#ifdef shapiro
# include <voldat.h>
# include <filtdat.h>
#endif
#include <options.h>
#ifdef oias
# include <oiopts.h>
#endif
#include <iounits.h>
#include <hybrid.h>
#include <vertslabs.h>
#include <extra.h>
#include <cdiag.h>
#ifdef surfpress
# include <workspb.h>
# include <sinfo.h>
#endif
#if defined ldrifters & defined rmdenbar
# include <rhomean.h>
#endif
#if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
# include <cbiopnh.h>
#endif
#ifdef ext_tide
# include <tidesp.h>
#endif
c
c-----------------------------------------------------------------------
c Define local data.
c-----------------------------------------------------------------------
c
integer i,im1,ip1,j,k,m,n,ndiskx
#ifdef shapiro
integer nn
# if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
integer ip
# endif
#endif
#ifdef surfpress
integer isave
#endif
#if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
integer nnegtrc
FLOAT
* negtrc, postrc
#endif
#ifdef peprf
logical wrtprf
#endif
FLOAT
* boot(imt,lbc),boot_save(imt,lbc),fx
#if defined surfpress | defined mfluxavg
FLOAT
* vcc,vcw,vcn,vcnw,vce,vcne,vedgesouth
integer jss
#endif
#if !defined surfpress & !defined mfluxavg
FLOAT
* vbaredge(imt)
#endif
#if defined surfpress
FLOAT
*uf(imt,jmt),vf(imt,jmt),fac
#endif
#if defined usrdiagnostic & defined nesttime
FLOAT
& dcpg(2),dco(2),dlp(2),dtrm(2),dtsr(2),tcpg(2),tco(2),tlp(2),
& tmsc(2),tmsc2(2),tmsc3(2),trlx(2),tsn1(2),tsn2(2),ttrm(2),
& ttsr(2)
# ifdef iforttime
& ,dum
FLOAT
& dtime
# endif
#endif
#ifdef cod_ing
integer ntimes,norder
#endif
c
#if defined usrdiagnostic & defined nesttime
data dcpg,dco,dlp,dtrm,dtsr,tcpg,tco,tlp,tmsc,
& tmsc2,tmsc3,trlx,tsn1,tsn2,ttrm,ttsr /32*c0/
c
#endif
c=======================================================================
c Begin executable code.
c=======================================================================
c
c=======================================================================
c Begin section for the initialization of ============================
c various quantities on each timestep ============================
c=======================================================================
c
c-----------------------------------------------------------------------
c Update timestep counter and total elapsed time.
c-----------------------------------------------------------------------
c
itt=itt+1
ttsec=ttsec+dtts
c
c-----------------------------------------------------------------------
c Set switches to indicate mixing timestep, to activate diagnostics, to
c print single line information, and to write out progressing solution.
c-----------------------------------------------------------------------
c
mixts=.false.
mxpas2=.false.
#if defined oias & !defined surfpress
if((mod(itt,nmix).eq.1).or.
* ((noi.gt.0).and.(ittoi.eq.2))) mixts=.true.
#else
if(mod(itt,nmix).eq.1) mixts=.true.
#endif
eots=.true.
if(mixts.and.eb) eots=.false.
diagts=.false.
if(mod(itt-1,nnergy).eq.0) diagts=.true.
prntsi=.false.
if(mod(itt,ntsi).eq.0) prntsi=.true.
wrtts=.false.
if(mod(itt-1,ntsout).eq.0) wrtts=.true.
c
#ifndef surfpress
c-----------------------------------------------------------------------
c If appropriate, initialize barotropic velocity field.
c-----------------------------------------------------------------------
c
# ifdef oias
if(itt.eq.1.or.ittoi.eq.2) call baro_vel
# else
if(itt.eq.1) call baro_vel
# endif
c
#endif
#if defined surfpress & defined freesurf & defined rmdenbar
c-----------------------------------------------------------------------
c Reset background mean density
c-----------------------------------------------------------------------
c
call meanrho
#endif
c
c-----------------------------------------------------------------------
c If appropriate, read in new boundary information.
c-----------------------------------------------------------------------
c
call infld(0)
c
c-----------------------------------------------------------------------
c Update permuting disc I/O units.
c-----------------------------------------------------------------------
c
ndiskb=mod(itt+2,3)+1
ndisk =mod(itt+0,3)+1
ndiska=mod(itt+1,3)+1
c
c-----------------------------------------------------------------------
c Adjust various quantities for mixing timestep.
c-----------------------------------------------------------------------
c
mxp=0
if(mixts) then
mix=1
c2dtts=dtts
c2dtuv=dtuv
c2dtsf=dtsf
do 10 j=1,jmt
do 10 i=1,imt
#if !defined surfpress | !defined freesurf
pb(i,j)=p(i,j)
#endif
ubarob(i,j)=ubaro(i,j)
vbarob(i,j)=vbaro(i,j)
10 continue
else
mix=0
c2dtts=c2*dtts
c2dtuv=c2*dtuv
c2dtsf=c2*dtsf
endif
#if defined surfpress & defined freesurf
if (itt.eq.1) c2dtps=c2dtsf
#endif
#ifdef oias
c
c If this is the end of an assimilation cycle, deactivate assimilation
c switch ITTOI.
c
if((noi.gt.0).and.(ittoi.eq.2)) then
ittoi=0
write(stdout,*) ' End of assimilation cycle = ',icoi-1
endif
#endif
c
20 continue
c
c-----------------------------------------------------------------------
c Establish over dimensioned arrays for vectorization.
c-----------------------------------------------------------------------
c
do 30 k=1,km
do 30 i=1,imt
dxtq (i,k)=dxt (i)
dxt4rq(i,k)=dxt4r(i)
dxuq (i,k)=dxu (i)
dxu2rq(i,k)=dxu2r(i)
30 continue
c
#if defined ext_tide
c-----------------------------------------------------------------------
c Compute tidal velocities along row J=1
c-----------------------------------------------------------------------
c
call get_tide (1)
c
# if defined advtide
c-----------------------------------------------------------------------
c Compute change in box vert. thickness due to free surf. displacement
c-----------------------------------------------------------------------
c
call reset_t_thickness
c
# endif
#endif
c-----------------------------------------------------------------------
c Get depths and vertical thicknesses at tracer points (rows=1,2)
c and velocity points (row=1).
c-----------------------------------------------------------------------
c
call setvert(1)
c
c-----------------------------------------------------------------------
c Queue up disk reads for this timestep.
c-----------------------------------------------------------------------
c
do 40 j=1,jmtm1
call ofind(labs(ndiskb),nslab,(j-1)*nslab+1)
call ofind(labs(ndisk),nslab,(j-1)*nslab+1)
40 continue
c
c-----------------------------------------------------------------------
c Initialize variables associated with energy diagnostics.
c-----------------------------------------------------------------------
c
if(prntsi) then
do 60 j=1,jmt
do 60 k=0,km
ektot(k,j)=c0
do 50 m=1,nt
dtabs(k,m,j)=c0
50 continue
60 continue
endif
if(diagts.and.eots.and.wnetcdf) call diag(1)
c
c=======================================================================
c End of section for initialization ==================================
c=======================================================================
c
c=======================================================================
c Begin a bootstrap procedure to prepare for the =====================
c row-by-row computation of prognostic variables =====================
c=======================================================================
c
c-----------------------------------------------------------------------
c Fetch data for row 2 from the disc.
c-----------------------------------------------------------------------
c
call oget(labs(ndiskb),nslab, +1,tb)
call oget(labs(ndisk ),nslab, +1,t )
call oget(labs(ndiskb),nslab,nslab+1,tbp)
call oget(labs(ndisk ),nslab,nslab+1,tp )
c
c-----------------------------------------------------------------------
c Switch slab incidental data into correct slab after read in.
c-----------------------------------------------------------------------
c
if(mod(itt,2)+mxp.ne.1) then
do 70 i=1,imt
bcon(i,1)=fkmup(i)
fkmup(i)=fkmtp(i)
fkmtp(i)=bcon(i,1)
boot(i,1)=fkmu(i)
fkmu(i)=fkmt(i)
fkmt(i)=boot(i,1)
70 continue
do 80 i=1,imt
bcon(i,2)=wsyp(i)
wsyp(i)=wsxp(i)
wsxp(i)=bcon(i,2)
boot(i,2)=wsy(i)
wsy(i)=wsx(i)
wsx(i)=boot(i,2)
80 continue
do 83 i=1,imt
boot_save(i,1)=bcon(i,1)
boot_save(i,2)=bcon(i,2)
83 continue
else
do 87 i=1,imt
boot_save(i,1)=fkmu(i)
boot_save(i,2)=wsy(i)
87 continue
endif
c
c-----------------------------------------------------------------------
c Convert maximum level indicators to integers.
c-----------------------------------------------------------------------
c
do 90 i=1,imt
kmtp(i)=nint(fkmtp(i))
kmup(i)=nint(fkmup(i))
90 continue
c
c-----------------------------------------------------------------------
c Move TAU-1 data to TAU level on a mixing timestep.
c-----------------------------------------------------------------------
c
if(mixts) then
do 100 m=1,nt
do 100 k=1,km
do 100 i=1,imt
tbp(i,k,m)=tp(i,k,m)
tb(i,k,m)=t(i,k,m)
100 continue
do 110 k=1,km
do 110 i=1,imt
ubp(i,k)=up(i,k)
vbp(i,k)=vp(i,k)
ub(i,k)=u(i,k)
vb(i,k)=v(i,k)
110 continue
endif
#ifdef ldrifters
c
c-----------------------------------------------------------------------
c Load first two slabs of density into volume storage (for drifters).
c-----------------------------------------------------------------------
c
call state(t(1,1,1),t(1,1,2),tdepth(1,1,jrs),rhos)
call state(tp(1,1,1),tp(1,1,2),tdepth(1,1,jrn),rhon)
# if defined rmdenbar
do 120 k=1,km
do 120 i=1,imt
rhos(i,k)=rhos(i,k)-rhobar(i,1,k)
rhon(i,k)=rhon(i,k)-rhobar(i,2,k)
120 continue
# endif
call load_sig(1,rhos)
call load_sig(2,rhon)
#endif
c
c-----------------------------------------------------------------------
c Initialize arrays for first calls to CLINIC and TRACER.
c-----------------------------------------------------------------------
c
do 130 k=1,km
do 130 i=1,imt
fvst(i,k)=c0
rhos(i,k)=c0
rhon(i,k)=c0
fmm (i,k)=c0
fm (i,k)=c0
gm (i,k)=c0
Umet(i,k)=c0
Uxx (i,k)=c0
Uyy (i,k)=c0
Vmet(i,k)=c0
Vxx (i,k)=c0
Vyy (i,k)=c0
Txx (i,k)=c0
Tyy (i,k)=c0
130 continue
c
c-----------------------------------------------------------------------
c Construct mask array for row J=2 tracers.
c-----------------------------------------------------------------------
c
do 140 k=1,km
do 140 i=1,imt
if(kmtp(i).ge.kar(k)) then
fmp(i,k)=c1
else
fmp(i,k)=c0
endif
140 continue
c
#ifndef surfpress
c-----------------------------------------------------------------------
c Set vorticity computation arrays at southern wall.
c-----------------------------------------------------------------------
c
do 150 i=1,imt
zus(i)=c0
zvs(i)=c0
150 continue
c
#else
c-----------------------------------------------------------------------
c Set roll vertically averaged momentum forcings.
c-----------------------------------------------------------------------
c
if(itt.gt.1) then
do 145 j=1,jmt
do 145 i=1,imt
zumb(i,j)=zum(i,j)
zvmb(i,j)=zvm(i,j)
145 continue
endif
c
do 150 j=1,jmt
do 150 i=1,imt
zum(i,j)=c0
zvm(i,j)=c0
150 continue
c
#endif
c-----------------------------------------------------------------------
c Save internal mode velocities for row J=2.
c-----------------------------------------------------------------------
c
do 160 k=1,km
do 160 i=1,imt
ucl(i,k)=u(i,k)
vcl(i,k)=v(i,k)
uclb(i,k)=ub(i,k)
vclb(i,k)=vb(i,k)
uclp(i,k)=up(i,k)
vclp(i,k)=vp(i,k)
uclbp(i,k)=ubp(i,k)
vclbp(i,k)=vbp(i,k)
160 continue
c
c-----------------------------------------------------------------------
c Compute advective coefficient for south face of row J=2 U,V boxes.
c-----------------------------------------------------------------------
c
#if defined surfpress | defined mfluxavg
fx=dyu2r(2)*csr(2)
do 231 k=1,km
do 231 i=2,imt-1
jss=1
ip1=min(i+1,imu)
# if !defined surfpress | !defined freesurf
vcc=dzvqz(i,k,0)*(vcl(i,k)+vbaro(i,jss))*cs(jss)
* *min(c1,fkmq(i ,jss ))
vcw=dzvqz(i-1,k,0)*(vcl(i-1,k)+vbaro(i-1,jss))*cs(jss)
* *min(c1,fkmq(i-1,jss ))
vcn=xzvqz(i,k)*(vclp(i,k)+vbaro(i,jss+1))*cs(jss+1)
* *min(c1,fkmq(i ,jss+1))
vcnw=xzvqz(i-1,k)*(vclp(i-1,k)+vbaro(i-1,jss+1))*cs(jss+1)
* *min(c1,fkmq(i-1,jss+1))
vce=dzvqz(ip1,k,0)*(vcl(ip1,k)+vbaro(ip1,jss))*cs(jss)
* *min(c1,fkmq(ip1,jss ))
vcne=xzvqz(ip1,k)*(vclp(ip1,k)+vbaro(ip1,jss+1))*cs(jss+1)
* *min(c1,fkmq(ip1,jss+1))
# else
vcc=(dzvqz(i,k,0)*vcl(i,k)+p5*dzvqzb(i,k,0)*(vbaro(i,jss)+
* vbarob(i,jss)))*cs(jss)*min(c1,fkmq(i ,jss ))
vcw=(dzvqz(i-1,k,0)*vcl(i-1,k)+dzvqzb(i-1,k,0)*(vbaro(i-1,jss)+
* vbarob(i-1,jss))*p5)*cs(jss)*min(c1,fkmq(i-1,jss ))
vcn=(xzvqz(i,k)*vclp(i,k)+p5*xzvqzb(i,k)*(vbaro(i,jss+1)+
* vbarob(i,jss+1)))*cs(jss+1)*min(c1,fkmq(i ,jss+1))
vcnw=(xzvqz(i-1,k)*vclp(i-1,k)+xzvqzb(i-1,k)*(vbaro(i-1,jss+1)+
* vbarob(i-1,jss+1))*p5)*cs(jss+1)*min(c1,fkmq(i-1,jss+1))
vce=(dzvqz(ip1,k,0)*vcl(ip1,k)+dzvqzb(i+1,k,0)*(vbaro(ip1,jss)+
* vbarob(ip1,jss))*p5)*cs(jss)*min(c1,fkmq(i+1,jss ))
vcne=(xzvqz(ip1,k)*vclp(ip1,k)+xzvqzb(i+1,k)*(vbaro(ip1,jss+1)+
* vbarob(ip1,jss+1))*p5)*cs(jss+1)*min(c1,fkmq(i+1,jss+1))
# endif
c
vedgesouth=(vcnw+vcne+vce+vcw+c2*(vcc+vcn))/c8
c
fvsu(i,k)=vedgesouth*fx
c
231 continue
# ifdef cyclic
do 232 k=1,km
fvsu(1,k)=fvsu(imtm1,k)
fvsu(imt,k)=fvsu(2,k)
232 continue
# else
do 232 k=1,km
fvsu(1,k)=fvsu(2,k) + dxt(1)/dxt(2)*(fvsu(2,k)-fvsu(3,k))
232 continue
# endif
#else
fx=dyu2r(2)*csr(2)*cst(2)
do 231 k=1,km
do 231 i=1,imt
ip1=min(i+1,imt)
# ifndef baroedge
vbaredge(i)=c2*(p(ip1,2)-p(i,2))
* *min(c1,fkmu(i ),fkmup(i))*dxur(i)/
* (hdv(i,2)+hdv(i,1))*cstr(2)
# else
vbaredge(i)=p5*(vbaro(i,1)+vbaro(i,2))
# endif
231 continue
# ifdef cyclic
vbaredge(1 )=vbaredge(imtm1)
vbaredge(imt)=vbaredge(2 )
# endif
c
do 232 k=1,km
do 232 i=1,imt
fvsu(i,k)=((vp(i,k)*xzvqz(i,k)+v(i,k)*dzvqz(i,k,0))*p5+
* p5*(xzvqz(i,k)+dzvqz(i,k,0))*vbaredge(i))*fx
c
232 continue
#endif
c
c-----------------------------------------------------------------------
c Add internal mode to external mode for row J=1.
c-----------------------------------------------------------------------
c
do 190 k=1,km
do 190 i=1,imu
ub(i,k)=ub(i,k)+ubarob(i,1)
vb(i,k)=vb(i,k)+vbarob(i,1)
u(i,k)=u(i,k)+ubaro(i,1)
v(i,k)=v(i,k)+vbaro(i,1)
190 continue
c
#ifndef surfpress
c-----------------------------------------------------------------------
c Construct flux for S, T on the southern faces of T-S boxes
c-----------------------------------------------------------------------
c
fx=cstr(2)*dytr(2)*cs(1)
do 195 k=1,km
do 195 i=1,imu
im1=max(i-1,1)
# if !defined ext_tide | !defined advtide
fvst(i,k)=(v(i,k)*dxuq(i,k)*dzvqz(i,k,0)+
* v(im1,k)*dxuq(im1,k)*dzvqz(im1,k,0))*
* fx*dxt4rq(i,k)
195 continue
# else
fvst(i,k)=(v(i,k)*dxuq(i,k)*dzvqz(i,k,0)*ustretch(i,1)+
* v(im1,k)*dxuq(im1,k)*dzvqz(im1,k,0)*ustretch(im1,1))*
* fx*dxt4rq(i,k)
195 continue
c
c-----------------------------------------------------------------------
c Compute tidal contributions for row J=1.
c-----------------------------------------------------------------------
c
do 200 k=1,km
do 200 i=1,imu
im1=max(i-1,1)
fvstdt(i,k)=
* (vtide (i ,k)*dxuq(i ,k)*dzvqz(i ,k,0)*ustretch(i,1)+
* vtide (im1,k)*dxuq(im1,k)*dzvqz(im1,k,0)*ustretch(im1,1))*
* fx*dxt4rq(i,k)
# ifdef advtide0
fvstdt(i,k)=fvstdt(i,k)*sadv
# endif
200 continue
# ifdef cyclic
do 201 k=1,km
fvstdt(1,k)=fvstdt(imtm1,k)
fvstdt(imt,k)=fvstdt(2,k)
201 continue
# endif
# endif
c
#else
c-----------------------------------------------------------------------
c Construct flux for S, T on the southern faces of T-S boxes
c-----------------------------------------------------------------------
c
fx=cstr(2)*dytr(2)*cs(1)
do 200 k=1,km
do 200 i=1,imu
im1=max(i-1,1)
# if !defined freesurf
fvst(i,k)=(v(i,k)*dxuq(i,k)*dzvqz(i,k,0)+
* v(im1,k)*dxuq(im1,k)*dzvqz(im1,k,0))*
* fx*dxt4rq(i,k)
# else
fvst(i,k)=((vcl(i,k)*dzvqz(i,k,0) + dzvqzb(i,k,0)*
* (vbaro(i,1)+vbarob(i,1))*p5)*dxuq(i,k) +
* (vcl(im1,k)*dzvqz(im1,k,0) + dzvqzb(im1,k,0)*
* (vbaro(im1,1)+vbarob(im1,1))*p5)*dxuq(im1,k) )*
* fx*dxt4rq(i,k)
# endif
200 continue
c
#endif
#ifdef cyclic
do 205 k=1,km
fvst(1,k)=fvst(imtm1,k)
fvst(imt,k)=fvst(2,k)
205 continue
c
#endif
c-----------------------------------------------------------------------
c Add external mode to internal mode for row J=2 (ocean points only).
c-----------------------------------------------------------------------
c
do 230 k=1,km
do 230 i=1,imu
if(kmup(i).ge.kar(k)) then
ubp(i,k)=ubp(i,k)+ubarob(i,2)
vbp(i,k)=vbp(i,k)+vbarob(i,2)
up (i,k)=up (i,k)+ubaro (i,2)
vp (i,k)=vp (i,k)+vbaro (i,2)
endif
230 continue
c
c-----------------------------------------------------------------------
c Set vertical boundary conditions (surface and bottom) for momentum
c and tracer at row J=1.
c-----------------------------------------------------------------------
c
#ifndef analytical
call setvbc(1)
#else
call anavbc(1)
#endif
c
c-----------------------------------------------------------------------
c Set vertical mixing coefficients for row J=1.
c-----------------------------------------------------------------------
c
#ifdef ppvmix
call ppmix(1)
#else
call cnvmix(1)
#endif
#ifdef peprf
c
c-----------------------------------------------------------------------
c For southernmost slab, write any requested profiles in ASCII format.
c-----------------------------------------------------------------------
c
call time_prf (itt,wrtprf)
if (wrtprf) call slab_prf (itt,1)
#endif
c
c=======================================================================
c End of bootstrap procedure =========================================
c=======================================================================
c
c=======================================================================
c Begin row-by-row computation of prognostic variables ===============
c=======================================================================
#if defined usrdiagnostic & defined nesttime
c
do j = 1, 2
tcpg(j) = c0
tco(j) = c0
tlp(j) = c0
ttrm(j) = c0
ttsr(j) = c0
tsn1(j) = c0
tsn2(j) = c0
enddo
# ifndef iforttime
call dtime (tmsc)
# else
dum = dtime(tmsc)
# endif
#endif
c
do 450 j=2,jmtm2
c
c-----------------------------------------------------------------------
c Set boundary conditions and then save slabs
c-----------------------------------------------------------------------
c
if(j.gt.2) then
#if !defined cyclic & !defined nest2larger
c
c Set western and eastern lateral boundary conditons on UA,VA,TA,P,ZTD
c for row=J-1.
c
call boundary(j-1,1)
#endif
#ifndef surfpress
c
c Reset southern boundary condition for vorticity tendency
c
if ((j.eq.5).and.(iopt(4).eq.3)) then
call zrobc_ori(j,south)
elseif ((j.eq.5).and.(iopt(4).eq.7).and.qcanpo) then
call zrobc_por(j,south)
endif
#endif
c
c Save data for row=J-1.
c
#ifdef shapiro
call osav(j-1)
#else
# if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
c
c Check for negative tracers, print warning if found.
c
if ((biopos.eq.3).or.(iopt(5).ne.0)) then
do 270 m=3,nt
nnegtrc=0
negtrc=c0
postrc=c0
do 260 k=1,km
do 260 i=1,imt
if(ta(i,k,m).lt.c0) then
nnegtrc=nnegtrc+1
negtrc=negtrc+ta(i,k,m)
else
postrc=postrc+ta(i,k,m)
endif
260 continue
if (nnegtrc.gt.0) then
write(stdout,900)nnegtrc,m
write(stdout,910)negtrc,postrc
if(biopos.eq.3) call exitus('STEP')
endif
270 continue
endif
c
c Insure non-negative biological tracers of row=J-1.
c
if(biopos.eq.1) then
do 280 m=3,nt
do 280 k=1,km
do 280 i=1,imt
ta(i,k,m)=max(c0,ta(i,k,m))
280 continue
endif
# endif
# if defined nest2larger | defined nest2smaller | defined AsselinFilt | defined surfpress
call osav(j-1)
# else
call oput(labs(ndiska),nslab,(j-2)*nslab+1,ta)
# endif
#endif
endif
c
if(j.eq.3) then
c
#ifndef nest2larger
c Set southern lateral boundary conditions on UA,VA,TA,P,ZTD.
c
call setvert(1)
call boundary(1,2)
call setvert(j-1)
c
#endif
c Set incidental data for row 1 and save present values of bcon in boot
c
do 285 i=1,imt
boot(i,1)=bcon(i,1)
boot(i,2)=bcon(i,2)
bcon(i,1)=boot_save(i,1)
bcon(i,2)=boot_save(i,2)
285 continue
c
c Save data for row=1.
c
#ifdef shapiro
call osav(1)
#else
# if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
c
c Check for negative tracers, print warning if found.
c
if ((biopos.eq.3).or.(iopt(5).ne.0)) then
do 300 m=3,nt
nnegtrc=0
negtrc=c0
postrc=c0
do 290 k=1,km
do 290 i=1,imt
if(ta(i,k,m).lt.c0) then
nnegtrc=nnegtrc+1
negtrc=negtrc+ta(i,k,m)
else
postrc=postrc+ta(i,k,m)
endif
290 continue
if (nnegtrc.gt.0) then
write(stdout,900)nnegtrc,m
write(stdout,910)negtrc,postrc
if(biopos.eq.3) call exitus('STEP')
endif
300 continue
endif
c
c Insure non-negative biological tracers for row=2.
c
if(biopos.eq.1) then
do 310 m=3,nt
do 310 k=1,km
do 310 i=1,imt
ta(i,k,m)=max(c0,ta(i,k,m))
310 continue
endif
# endif
# if defined nest2larger | defined nest2smaller | defined AsselinFilt | defined surfpress
call osav(1)
# else
call oput(labs(ndiska),nslab,1,ta)
# endif
#endif
c
c Reset bcon from values saved in boot
c
do 315 i=1,imt
bcon(i,1)=boot(i,1)
bcon(i,2)=boot(i,2)
315 continue
endif
c
c-----------------------------------------------------------------------
c Move all slab data down one row.
c-----------------------------------------------------------------------
c
do 320 n=1,nt
do 320 k=1,km
do 320 i=1,imt
tbm(i,k,n)=tb (i,k,n)
tm (i,k,n)=t (i,k,n)
tb (i,k,n)=tbp(i,k,n)
t (i,k,n)=tp (i,k,n)
320 continue
do 330 k=1,km
do 330 i=1,imt
ubm(i,k)=ub (i,k)
vbm(i,k)=vb (i,k)
um (i,k)=u (i,k)
vm (i,k)=v (i,k)
ub (i,k)=ubp(i,k)
vb (i,k)=vbp(i,k)
u (i,k)=up (i,k)
v (i,k)=vp (i,k)
#ifdef ext_tide
utidem(i,k)=utide (i,k)
vtidem(i,k)=vtide (i,k)
utide (i,k)=utidep(i,k)
vtide (i,k)=vtidep(i,k)
#endif
330 continue
do 340 i=1,imt
fkmum(i)=fkmu (i)
wsym (i)=wsy (i)
fkmtm(i)=fkmt (i)
wsxm (i)=wsx (i)
fkmu (i)=fkmup(i)
wsy (i)=wsyp (i)
fkmt (i)=fkmtp(i)
wsx (i)=wsxp (i)
#ifdef ext_tide
btide (i)=btidep(i)
#endif
340 continue
c
c-----------------------------------------------------------------------
c Complete read in of J+1 slab.
c-----------------------------------------------------------------------
c
call oget(labs(ndiskb),nslab,j*nslab+1,tbp)
call oget(labs(ndisk ),nslab,j*nslab+1,tp )
c
c-----------------------------------------------------------------------
c Switch slab incidental data into correct slab after read in.
c-----------------------------------------------------------------------
c
if(mod(itt,2)+mxp.ne.1) then
do 370 i=1,imt
bcon(i,1)=fkmup(i)
fkmup(i)=fkmtp(i)
fkmtp(i)=bcon(i,1)
bcon(i,2)=wsyp(i)
wsyp(i)=wsxp(i)
wsxp(i)=bcon(i,2)
370 continue
endif
c
c-----------------------------------------------------------------------
c Shift maximum level indicators down one row and set J+1 values.
c-----------------------------------------------------------------------
c
do 380 i=1,imt
kmt (i)=kmtp (i)
kmu (i)=kmup (i)
kmtp(i)=nint(fkmtp(i))
kmup(i)=nint(fkmup(i))
if(j.eq.jmtm2) kmtp(i)=kmt(i)
if(j.eq.jmtm2) kmup(i)=kmu(i)
380 continue
c
c-----------------------------------------------------------------------
c Move TAU-1 data to TAU level on a mixing timestep.
c-----------------------------------------------------------------------
c
if(mixts) then
do 390 m=1,nt
do 390 k=1,km
do 390 i=1,imt
tbp(i,k,m)=tp(i,k,m)
390 continue
do 400 k=1,km
do 400 i=1,imt
ubp(i,k)=up(i,k)
vbp(i,k)=vp(i,k)
400 continue
endif
c
c-----------------------------------------------------------------------
c Shift masks down one row and compute new masks.
c-----------------------------------------------------------------------
c
do 410 k=1,km
do 410 i=1,imt
fmm(i,k)=fm (i,k)
fm (i,k)=fmp(i,k)
410 continue
do 420 k=1,km
do 420 i=1,imt
if(kmtp(i).ge.kar(k)) then
fmp(i,k)=c1
else
fmp(i,k)=c0
endif
if(kmu(i).ge.kar(k)) then
gm(i,k)=c1
else
gm(i,k)=c0
endif
420 continue
c
c-----------------------------------------------------------------------
c Get depths and vertical thicknesses at the tracer points (row=J,J+1)
c and velocity points (row=J,J-1).
c-----------------------------------------------------------------------
c
call setvert(j)
c
c-----------------------------------------------------------------------
c Set vertical boundary conditions (surface and bottom) for momentum
c and tracers.
c-----------------------------------------------------------------------
c
#ifndef analytical
call setvbc(j)
#else
call anavbc(j)
#endif
#ifdef ext_tide
c-----------------------------------------------------------------------
c Compute tidal velocities along row J
c-----------------------------------------------------------------------
c
call get_tide (j)
c
#endif
c
c-----------------------------------------------------------------------
c Set vertical mixing coefficients.
c-----------------------------------------------------------------------
c
#ifdef ppvmix
call ppmix(j)
#else
call cnvmix(j)
#endif
c
c-----------------------------------------------------------------------
c Monitor velocities at the current time level.
c-----------------------------------------------------------------------
c Too Early (incorrect w,wu) PJH
c call cfl(j)
c
c-----------------------------------------------------------------------
c Call the main computational routines to update the row.
c-----------------------------------------------------------------------
c
#if !defined usrdiagnostic | !defined nesttime
call clinic(j)
call tracer(j)
#else
# ifndef iforttime
call dtime (dlp)
# else
dum = dtime(dlp)
# endif
call clinic(j,dcpg,dco)
call tracer(j,dtsr,dtrm)
do n = 1, 2
tcpg(n) = tcpg(n) + dcpg(n)
tco(n) = tco(n) + dco(n)
tlp(n) = tlp(n) + dlp(n)
ttrm(n) = ttrm(n) + dtrm(n)
ttsr(n) = ttsr(n) + dtsr(n)
enddo
#endif
c
c-----------------------------------------------------------------------
c Monitor velocities at the current time level. This call occurs
c after the calls to clinic & tracer since the correct vertical
c velocities are not available until then.
c-----------------------------------------------------------------------
c
call cfl(j)
c
c-----------------------------------------------------------------------
c Calculate diagnostics on diagnostic timesteps. Only if number of
c output levels NLEV is greater than zero.
c-----------------------------------------------------------------------
c
if(diagts.and.eots.and.wnetcdf) call diag(j)
c
#if defined analytical
c-----------------------------------------------------------------------
c Report other User dependent diagnostics.
c-----------------------------------------------------------------------
c
if(eots) call anadiag(j)
c
# elif defined usrdiagnostic & !defined nesttime
c-----------------------------------------------------------------------
c Report other User dependent diagnostics.
c-----------------------------------------------------------------------
c
if(eots) call userdiag(j)
c
#endif
c-----------------------------------------------------------------------
c Write out the progressing solution at specified rows on energy TSTEP.
c Only if number of output levels NLEV is greater than zero.
c-----------------------------------------------------------------------
c
if(wrtts.and.eots.and.wnetcdf) call cdfout(j)
c
#ifdef peprf
c-----------------------------------------------------------------------
c For current slab, write any requested profiles in ASCII format.
c-----------------------------------------------------------------------
c
if (wrtprf) call slab_prf (itt,j)
c
#endif
c-----------------------------------------------------------------------
c Put slab incidental data into correct slab for writeout.
c-----------------------------------------------------------------------
c
if(mod(itt,2).eq.0) then
do 430 i=1,imt
bcon(i,1)=fkmt(i)
bcon(i,2)=wsx(i)
430 continue
else
do 440 i=1,imt
bcon(i,1)=fkmu(i)
bcon(i,2)=wsy(i)
440 continue
endif
#ifdef ldrifters
c
c-----------------------------------------------------------------------
c Load density for slab J+1 into volume storage for drifters.
c-----------------------------------------------------------------------
c
call load_sig(j+1,rhon)
#endif
450 continue
c
#if defined usrdiagnostic & defined nesttime
# ifndef iforttime
call dtime (dlp)
# else
dum = dtime(dlp)
# endif
tlp(1) = tlp(1) + dlp(1)
tlp(2) = tlp(2) + dlp(2)
c
#endif
c=======================================================================
c End row-by-row computation =========================================
c=======================================================================
c
c-----------------------------------------------------------------------
c Modify UA,VA,TA at I=1,I=IMT (IMTM1) boundaries and
c save newly computed data from the final row.
c-----------------------------------------------------------------------
#if !defined cyclic & !defined nest2larger
c
c Set western and eastern lateral boundary conditons on UA,VA,TA,P,ZTD
c for row=JMTM2.
c
call boundary(jmtm2,1)
#endif
c
c Save data for row=JMTM2
c
#ifdef shapiro
call osav(jmtm2)
#else
# if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
c
c Check for negative tracers, print warning if found.
c
if ((biopos.eq.3).or.(iopt(5).ne.0)) then
do 470 m=3,nt
nnegtrc=0
negtrc=c0
postrc=c0
do 460 k=1,km
do 460 i=1,imt
if(ta(i,k,m).lt.c0) then
nnegtrc=nnegtrc+1
negtrc=negtrc+ta(i,k,m)
else
postrc=postrc+ta(i,k,m)
endif
460 continue
if (nnegtrc.gt.0) then
write(stdout,900)nnegtrc,m
write(stdout,910)negtrc,postrc
if(biopos.eq.3) call exitus('STEP')
endif
470 continue
endif
c
c Insure non-negative biological tracers for row=JMTM2.
c
if(biopos.eq.1) then
do 480 m=3,nt
do 480 k=1,km
do 480 i=1,imt
ta(i,k,m)=max(c0,ta(i,k,m))
480 continue
endif
# endif
# if defined nest2larger | defined nest2smaller | defined AsselinFilt | defined surfpress
call osav(jmtm2)
# else
call oput(labs(ndiska),nslab,(jmt-3)*nslab+1,ta)
# endif
#endif
c
#ifndef nest2larger
c Set northern lateral boundary conditons on UA,VA,TA,P,ZTD.
c
call setvert(jmtm1)
call boundary(jmtm1,4)
c
#endif
c Save data for row=JMTM1.
c
#ifdef shapiro
call osav(jmtm1)
#else
# if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
c
c Check for negative tracers, print warning if found.
c
if ((biopos.eq.3).or.(iopt(5).ne.0)) then
do 500 m=3,nt
nnegtrc=0
negtrc=c0
postrc=c0
do 490 k=1,km
do 490 i=1,imt
if(ta(i,k,m).lt.c0) then
nnegtrc=nnegtrc+1
negtrc=negtrc+ta(i,k,m)
else
postrc=postrc+ta(i,k,m)
endif
490 continue
if (nnegtrc.gt.0) then
write(stdout,900)nnegtrc,m
write(stdout,910)negtrc,postrc
if(biopos.eq.3) call exitus('STEP')
endif
500 continue
endif
c
c Insure non-negative biological tracers for row=JMT-1
c
if(biopos.eq.1) then
do 510 m=3,nt
do 510 k=1,km
do 510 i=1,imt
ta(i,k,m)=max(c0,ta(i,k,m))
510 continue
endif
# endif
# if defined nest2larger | defined nest2smaller | defined AsselinFilt | defined surfpress
call osav(jmtm1)
# else
call oput(labs(ndiska),nslab,(jmt-2)*nslab+1,ta)
# endif
#endif
#ifdef peprf
c
c-----------------------------------------------------------------------
c For northernmost slabs, write any requested profiles in ASCII format.
c-----------------------------------------------------------------------
c
if (wrtprf) call slab_prf (itt,jmtm1)
if (wrtprf) call slab_prf (itt,jmt)
c
#endif
#ifdef surfpress
c-----------------------------------------------------------------------
c Define initial values of vertically averaged baroclinic forcings.
c-----------------------------------------------------------------------
c
if (itt.eq.1) then
isave = iopt(4)
iopt(4) = 5
call boundary_zuvm
iopt(4) = isave
c
do 515 j=1,jmtm1
do 515 i=1,imu
zumb(i,j)=zum(i,j)
zvmb(i,j)=zvm(i,j)
515 continue
endif
c
#endif
#if defined nest2larger | defined nest2smaller
c-----------------------------------------------------------------------
# ifndef surfpress
# if !defined nest_ext2lrgr & !defined nest_ext2smlr
c Send and/or receive boundary conditions for barotropic vorticity.
# elif defined nest_ext2lrgr & defined nest_ext2smlr
c Send and/or receive transport for fine grid
# elif defined nest_ext2lrgr
c Receive transport from coarse grid. Send barotropic vorticity
c boundary conditions to fine grid.
# else
c Receive barotropic vorticity boundary conditions from coarse grid.
c Send transport to fine grid.
# endif
# else
# if !defined nest_ext2lrgr & !defined nest_ext2smlr
c Send and/or receive boundary conditions for barotropic divergence.
# elif defined nest_ext2lrgr & defined nest_ext2smlr
c Send and/or receive pressure for fine grid
# elif defined nest_ext2lrgr
c Receive pressure from coarse grid. Send barotropic divergence
c boundary conditions to fine grid.
# else
c Receive barotropic divergence boundary conditions from coarse grid.
c Send pressure to fine grid.
# endif
# endif
c Send and/or receive interior data from smaller grids to larger.
# ifdef nest2larger
c Receive boundary conditions from larger grid.
# endif
c-----------------------------------------------------------------------
c
# if !defined usrdiagnostic | !defined nesttime
call nest_t_align (itt)
call nest_interior (itt)
# ifdef nest2larger
call nest_rec_bc
# endif
# else
# ifndef iforttime
call dtime (tmsc2)
# else
dum = dtime(tmsc2)
# endif
call nest_t_align (itt,tsn1)
call nest_interior (itt,tsn1,tsn2)
# ifdef nest2larger
call nest_rec_bc (tsn1,tsn2)
# endif
# endif
c
#endif
#ifdef surfpress
c---------------------------------------------------------------------
c Set open boundary conditions for Pressure
c---------------------------------------------------------------------
c
# if defined freesurf & defined ext_tide & !defined mixtideonly
call bdy_tide_press
c
if(iopt(3).eq.7) then
# else
if(iopt(3).ge.3) then
# endif
call boundary_press
endif
c
c---------------------------------------------------------------------
c Extend the vertically averaged baroclinic forcing to the boundaries.
c---------------------------------------------------------------------
c
call boundary_zuvm
c
c-------------------------------------------------------------------
c Solve for the augmented (hat) Barotropic velocities (used in the
c splitting) if using the surface pressure formulation. Also apply
c boundary conditions to them and filter them.
c-------------------------------------------------------------------
c
call UV_hat
c
c------------------------------------------------------------------
c Get the RHS for the surface pressure Elliptic equation.
c------------------------------------------------------------------
c
call h_force(uf,vf,fac)
c
# if defined freesurf
call div_fld(uf,vf,ztd,hdvcf,fac)
# else
call div_fld(uf,vf,ztd,hdvc,fac)
# endif
c
#endif
#ifdef shapiro
c-----------------------------------------------------------------------
c If applicable, Shapiro filter the 2-d horizontal velocity and tracers
c fields.
c-----------------------------------------------------------------------
c
c ICNTM, ICNTH = counter (for frequency of application)
c NFRQM, NFRQH = frequency with which filter is applied
c NTIMM, NTIMH = number of times filter is applied per time step
c NORDM, NORDH = order of the filter
c
c Shapiro filter internal mode velocity.
c
if(mixvel.eq.1) then
icntv=icntv+1
if(icntv.ne.nfrqv) goto 530
if(nordv.eq.0) goto 530
icntv=0
do 520 nn=1,ntimv
do 520 k=1,km
call shap_lev(xu(1,k),imu,jmtm1,vgrid,nordv)
call shap_lev(xv(1,k),imu,jmtm1,vgrid,nordv)
520 continue
530 continue
endif
c
c Shapiro filter tracers.
c
if(mixtrc.eq.1) then
icntt=icntt+1
if(icntt.ne.nfrqt) goto 550
if(nordt.eq.0) goto 550
icntt=0
# if defined shapmean & defined pressbias
call shap_mean(0)
# endif
# ifndef cod_ing
do 540 nn=1,ntimt
do 540 m=1,nt
do 540 k=1,km
# if defined coast & defined coastedge
call set_edges(xt(1,k,m),tgrid)
# endif
call shap_lev(xt(1,k,m),imt,jmtm1,tgrid,nordt)
# if defined coast & defined coastedge
call set_edges(xt(1,k,m),tgrid)
# endif
540 continue
# else
do 540 m=1,nt
do 540 k=1,km
if(m.eq.icod) then
ntimes=10
if(k.gt.10)ntimes=ntimes+(k-9)
norder=2
else
ntimes=ntimt
norder=nordt
endif
do 540 nn=1,ntimes
# if defined coast & defined coastedge
call set_edges(xt(1,k,m),tgrid)
# endif
call shap_lev(xt(1,k,m),imt,jmtm1,tgrid,norder)
# if defined coast & defined coastedge
call set_edges(xt(1,k,m),tgrid)
# endif
540 continue
# endif
# if defined shapmean & defined pressbias
call shap_mean(1)
# endif
550 continue
endif
# if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
c
c Check for negative tracers, print warning if found.
c
if ((biopos.eq.3).or.(iopt(5).ne.0)) then
do 570 m=3,nt
nnegtrc=0
negtrc=c0
postrc=c0
do 560 k=1,km
do 560 j=1,jmtm1
do 560 i=1,imt
ip=i+(j-1)*imt
if(xt(ip,k,m).lt.c0) then
nnegtrc=nnegtrc+1
negtrc=negtrc+xt(ip,k,m)
else
postrc=postrc+xt(ip,k,m)
endif
560 continue
if (nnegtrc.gt.0) then
write(stdout,900)nnegtrc,m
write(stdout,910)negtrc,postrc
if(biopos.eq.3) call exitus('STEP')
endif
570 continue
endif
c
c Insure non-negative biological tracers.
c
if(biopos.eq.1) then
do 580 m=3,nt
do 580 k=1,km
do 580 j=1,jmtm1
do 580 i=1,imt
ip=i+(j-1)*imt
xt(ip,k,m)=max(c0,xt(ip,k,m))
580 continue
endif
# endif
#endif
c
c Put filtered data into disk.
c
#if !(defined AsselinFilt | defined surfpress) & (defined shapiro | defined nest2larger | defined nest2smaller)
call okeep(labs(ndiska))
#elif defined AsselinFilt | defined surfpress
call okeep(labs(ndiska),labs(ndisk ),labs(ndiskb))
#endif
c-----------------------------------------------------------------------
c Write out line integral analysis.
c-----------------------------------------------------------------------
c
if(prntsi.and.eots) then
do 610 j=2,jmtm2
do 590 k=km,1,-1
ektot(0,1)=ektot(0,1)+ektot(k,j)
590 continue
do 600 m=1,nt
do 600 k=km,1,-1
dtabs(0,m,1)=dtabs(0,m,1)+dtabs(k,m,j)
600 continue
610 continue
ektot(0,1)=ektot(0,1)/volume
do 620 m=1,nt
dtabs(0,m,1)=dtabs(0,m,1)/volume
620 continue
write(stdout,630) itt,ektot(0,1),dtabs(0,1,1),dtabs(0,2,1)
630 format(1x,'TS = ',i7,2x,'KE = ',1pe13.6,2x,'Dtemp = ',1pe13.6,
* 2x,'Dsalt = ',1pe13.6)
endif
c
#if defined usrdiagnostic & defined nesttime
# ifndef iforttime
call dtime (tmsc3)
# else
dum = dtime(tmsc3)
# endif
#endif
#if (!defined nest2larger | !defined nest_ext2lrgr) & !defined bioadjloc
c-----------------------------------------------------------------------
# ifndef surfpress
c Solve for stream function at advanced time level
# else
c Solve for the surface pressure at advanced time level
# endif
c-----------------------------------------------------------------------
c
# ifndef surfpress
# ifndef pcg5
call relax
# else
call streamfctn
# endif
# else
call pressure
# endif
c
#endif
c
c----------------------------------------------------------------------
c Update the Barotropic velocities and store the previous-time values.
c----------------------------------------------------------------------
c
#ifndef surfpress
call baro_vel
#else
call external_vel
#endif
c
#if defined surfpress & defined freesurf
c----------------------------------------------------------------------
c Update vertical grid parameters
c----------------------------------------------------------------------
c
call expand_depth
c
c Complete barotropic advection effects.
c
call reset_advect (labs(ndiska))
#endif
c
#if defined usrdiagnostic & defined nesttime
# ifndef iforttime
call dtime (trlx)
# else
dum = dtime(trlx)
# endif
c
#endif
#ifdef nest2smaller
c-----------------------------------------------------------------------
c Send boundary conditions to smaller grid.
c-----------------------------------------------------------------------
c
# if !defined usrdiagnostic | !defined nesttime
call nest_snd_bc (itt)
# else
call nest_snd_bc (itt,tsn1,tsn2)
# endif
c
#endif
#if defined usrdiagnostic & defined nesttime
c-----------------------------------------------------------------------
c Report other User dependent diagnostics.
c-----------------------------------------------------------------------
c
call userdiag (itt,tcpg,tco,ttrm,ttsr,tlp,tmsc,tmsc2,tmsc3,trlx,
& tsn1,tsn2)
c
#endif
c-----------------------------------------------------------------------
c If this is the end of the 1st pass of an euler backward timestep,
c set the input disc units so that the proper levels are fetched on
c the next pass. The output for the 2nd pass will be placed on the
c unused unit ("NDISKB") and transferred to the proper unit ("NDISKA")
c later. Return to the top of "STEP" to do the 2nd pass.
c-----------------------------------------------------------------------
c
if(mixts.and.eb) then
mix=0
mxp=1
eots=.true.
mixts=.false.
mxpas2=.true.
ndiskx=ndiskb
ndiskb=ndisk
ndisk=ndiska
ndiska=ndiskx
go to 20
endif
c
c-----------------------------------------------------------------------
c If this is the end of the 2nd pass of an Euler backward timestep,
c transfer the data written temporarily to "NDISKX" to its final
c destination (the original "NDISKA").
c-----------------------------------------------------------------------
c
if(mxpas2) then
ndiska=ndisk
ndisk=ndiskb
do 640 j=2,jmtm1
call oget(labs(ndiskx),nslab,(j-1)*nslab+1,ta)
call oput(labs(ndiska),nslab,(j-1)*nslab+1,ta)
640 continue
endif
c
c-----------------------------------------------------------------------
c For purposes of recovering from the disc after an abnormal stop,
c normally inactive disc units are brought up to date here.
c-----------------------------------------------------------------------
c
call oput(kflds,nwds,(ndiska-1)*nwds+1,p)
call oput(kontrl,20,1,itt)
c
#if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
900 format(/' WARNING - found negative biological tracer: ',i9,
* ' instances, tracer ',i2)
910 format(/' net negatives: ',1pe13.6,' net positive: ',1pe13.6)
#endif
return
end