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 #include #include #include #include #include #include #include #include #include #ifdef shapiro # include # include #endif #include #ifdef oias # include #endif #include #include #include #include #include #ifdef surfpress # include # include #endif #if defined ldrifters & defined rmdenbar # include #endif #if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse # include #endif #ifdef ext_tide # include #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