#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 #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifdef bndy_rlx # include #endif #ifdef linear_physics # include #endif #if defined ext_tide & defined advtide # include #endif #if defined codunlim | defined codlim # include #endif #if defined surfpress & defined freesurf # include # include # include #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