subroutine diag(j) c c======================================================================= c === c This routine initializes, calculates and writes out enegy and === c diagnostics. The energy/diagnostics are written to a NetCDF === c file. === c === c Calls: NCSNC, NCVID, NCVPT, NCVPT1 (NetCDF library) === c EXITUS, WRTCDF === 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 #include #include #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c logical first integer id,i,j,jj,k,ll,m,n,tindx integer count(4),start(4) FLOAT * basinare,basinvol,buoyerr,engleak,fx,giga,micro,peta,scl, * sscle,tera,tbrz,totdx,totdz,tscle,vbrz,vscle parameter (giga=c1em9,micro=c1e6,sscle=c1em10,peta=c1em15, * tera=c1em12,tscle=4.186e-15,vscle=c1em12) c save basinare,basinvol,first,tindx data first /.true./ c c======================================================================= c Begin executable code. c======================================================================= c c----------------------------------------------------------------------- c Initialize various quantities when J=1. c----------------------------------------------------------------------- c if (j.eq.1) then c c Initialize variables associated with the diagnostics of momentum. c avwv(1)=c0 avwv(2)=c0 avwv(3)=c0 smflx(1)=c0 smflx(2)=c0 do 120 ll=1,mterms engext(ll)=c0 do 100 i=1,imt zuseng(i,ll)=c0 zvseng(i,ll)=c0 100 continue do 110 k=0,km engint(k,ll)=c0 termbm(k,ll,1)=c0 termbm(k,ll,2)=c0 110 continue 120 continue do 130 jj=1,jmt do 130 k=1,km tmt(jj,k)=c0 130 continue c c Initialize variables associated with the diagnostics of tracers. c do 140 k=0,km buoyint(k)=c0 140 continue do 170 m=1,nt atwv(m)=c0 stflx(m)=c0 asst(m)=c0 do 150 k=0,km tvar(k,m)=c0 150 continue do 160 ll=1,tterms do 160 k=0,km termbt(k,ll,m)=c0 160 continue 170 continue do 180 m=1,ntmin2 do 180 jj=1,jmt do 180 ll=1,tterms ttn(ll,jj,m)=c0 180 continue c c Initialize CFL monitoring CFL parameters. c cflup =c0 cflvp =c0 cflwtp=c0 cflwup=c0 return endif c c----------------------------------------------------------------------- c Compute the zonally integrated meridional mass transport as well as c the northward transport of each tracer quantity. c----------------------------------------------------------------------- c c Compute meridional mass transport. c do 190 k=1,km vbr(k)=c0 do 190 i=2,imtm1 vbr(k)=vbr(k)+v(i,k)*dxu(i)*cs(j) if(k.eq.1) then tmt(j,1)=vbr(1)*dzvqz(i,1,0) else tmt(j,k)=tmt(j,k-1)+vbr(k)*dzvqz(i,k,0) endif 190 continue c c Compute northward transport of each tracer. c do 200 m=1,nt do 200 k=1,km tbrs(k,m)=tbrn(k,m) tbrn(k,m)=c0 200 continue c if(j.eq.2) then do 210 m=1,nt do 210 k=1,km tbrs(k,m)=c0 210 continue do 240 k=1,km totdx=c0 do 220 i=2,imtm1 totdx=totdx+dxt(i)*cst(j)*fm(i,k) do 220 m=1,nt tbrs(k,m)=tbrs(k,m)+ * t(i,k,m)*fm(i,k)*dxt(i)*cst(j)*dzqz(i,k,0) 220 continue if(totdx.ne.c0) then do 230 m=1,nt tbrs(k,m)=tbrs(k,m)/totdx 230 continue endif 240 continue endif c do 270 k=1,km totdx=c0 do 250 i=2,imtm1 totdx=totdx+dxt(i)*cst(j)*fm(i,k) do 250 m=1,nt tbrn(k,m)=tbrn(k,m)+ * tp(i,k,m)*fm(i,k)*cst(j)*dxt(i)*dzqz(i,k,1) 250 continue if(totdx.ne.c0) then do 260 m=1,nt tbrn(k,m)=tbrn(k,m)/totdx 260 continue endif 270 continue c fx=ah*dyur(j) do 280 m=1,nt do 280 k=1,km ttn(1,j,m)=ttn(1,j,m)+vbr(k)*p5*(tbrn(k,m)+tbrs(k,m)) do 280 i=2,imtm1 ttn(6,j,m)=ttn(6,j,m)+p25* * (v(i,k)*dxu(i)+v(i-1,k)*dxu(i-1))*cs(j)* * (t(i,k,m)*dzqz(i,k,0)+tp(i,k,m)*dzqz(i,k,1)) if((mixtrc.eq.2).or.(mixtrc.eq.3)) then ttn(7,j,m)=ttn(7,j,m)-fx*fm(i,k)*fmp(i,k)*(tp(i,k,m)-t(i,k,m)) * *dxt(i)*cs(j)*dzqz(i,k,0) endif 280 continue c do 300 m=1,nt do 300 i=2,imtm1 totdz=c0 vbrz=c0 tbrz=c0 do 290 k=1,km if((kmt(i).ge.k).and.(kmtp(i).ge.k)) then vbrz=vbrz+(v(i ,k)*dxu(i )*dzvqz(i ,k,0)+ * v(i-1,k)*dxu(i-1)*dzvqz(i-1,k,0)) tbrz=tbrz+(t(i,k,m)*dzqz(i,k,0)+tp(i,k,m)*dzqz(i,k,1)) totdz=totdz+dzqz(i,k,0) endif 290 continue if(totdz.ne.c0) then tbrz=tbrz/totdz ttn(3,j,m)=ttn(3,j,m)+p25*vbrz*tbrz*cs(j) ttn(5,j,m)=ttn(5,j,m)-(smf(i,1)*dxu(i)+smf(i-1,1)*dxu(i-1))* * (t(i,1,m)+tp(i,1,m)-tbrz)* * cs(j)/(c8*omega*sine(i,j)) endif 300 continue c do 310 m=1,nt ttn(2,j,m)=ttn(6,j,m)-ttn(1,j,m) ttn(4,j,m)=ttn(6,j,m)-ttn(3,j,m)-ttn(5,j,m) ttn(8,j,m)=ttn(6,j,m)+ttn(7,j,m) 310 continue c c======================================================================= c Write out volume data into energy/diagnostics NetCDF file. c======================================================================= c c Write ocean basin area and volume on first pass. c if(first) then basinare=area*cm2tom2 varid=ncvid(ncnrgid,'area',rcode) call ncvpt1(ncnrgid,varid,1,basinare,rcode) basinvol=volume*cm3tom3 varid=ncvid(ncnrgid,'volume',rcode) call ncvpt1(ncnrgid,varid,1,basinvol,rcode) first=.false. endif c c Write out time. c if(j.eq.2) then tnrgindx=((itt-1)/nnergy)+1 tindx=tnrgindx call ncvpt1(ncnrgid,tnrgid,tindx,ttsec-dtts,rcode) if(rcode.ne.0) then write(stdout,900) 'time' call exitus('DIAG') endif endif c c Write out buoyancy (gigaWatts). c if(iout(15).eq.1) then scl=cm2tom2*rho0*basinvol*giga call wrtcdf(j,buoy,imt,km,scl,buoyid,sindx,tindx,ncnrgid) if(j.eq.2) then call wrtcdf(j-1,buoy,imt,km,scl,buoyid,sindx,tindx,ncnrgid) elseif(j.eq.jmtm2) then call wrtcdf(j+1,buoy,imt,km,scl,buoyid,sindx,tindx,ncnrgid) call wrtcdf(j+2,buoy,imt,km,scl,buoyid,sindx,tindx,ncnrgid) endif endif c c Write out total kinetic energy (petaJoules). c if((iout(10).eq.1).or.(iout(10).eq.3)) then id=ketotid scl=cm2tom2*rho0*basinvol*peta call wrtcdf(j,eke(1,1,1),imt,km,scl,id,sindx,tindx,ncnrgid) if(j.eq.2) then call wrtcdf(j-1,eke(1,1,1),imt,km,scl,id,sindx,tindx,ncnrgid) elseif(j.eq.jmtm2) then call wrtcdf(j+1,eke(1,1,1),imt,km,scl,id,sindx,tindx,ncnrgid) call wrtcdf(j+2,eke(1,1,1),imt,km,scl,id,sindx,tindx,ncnrgid) endif endif c c Write time rate of change of kinetic energy per component per unit c volume (microWatts/m^3). c if((iout(10).eq.2).or.(iout(10).eq.3)) then id=ketrmid scl=cm2tom2*rho0*micro do 320 n=2,mterms-2 call wrtcdf(j,eke(1,1,n),imt,km,scl,id,n-1,tindx,ncnrgid) if(j.eq.2) then call wrtcdf(j-1,eke(1,1,n),imt,km,scl,id,n-1,tindx,ncnrgid) elseif(j.eq.jmtm2) then call wrtcdf(j+1,eke(1,1,n),imt,km,scl,id,n-1,tindx,ncnrgid) call wrtcdf(j+2,eke(1,1,n),imt,km,scl,id,n-1,tindx,ncnrgid) endif 320 continue endif c c Write time rate of change of tracers per component per unit c volume. c if(iout(20).eq.1) then do 323 m=1,nt do 323 k=1,km do 323 i=1,imt trmbtv(i,k,1,m)=trmbtv(i,k,2,m) do 323 n=3,tterms trmbtv(i,k,1,m)=trmbtv(i,k,1,m)+trmbtv(i,k,n,m) 323 continue scl=c1 do 327 m=1,nt id=trmtvid(m) do 327 n=1,tterms call wrtcdf(j,trmbtv(1,1,n,m),imt,km,scl,id,n,tindx, & ncnrgid) if(j.eq.2) then call wrtcdf(j-1,trmbtv(1,1,n,m),imt,km,scl,id,n,tindx, & ncnrgid) elseif(j.eq.jmtm2) then call wrtcdf(j+1,trmbtv(1,1,n,m),imt,km,scl,id,n,tindx, & ncnrgid) call wrtcdf(j+2,trmbtv(1,1,n,m),imt,km,scl,id,n,tindx, & ncnrgid) endif 327 continue endif c c======================================================================= c Compute globally averaged data into energy/diagnostics NetCDF file. c======================================================================= c if(j.lt.jmtm2) return c c Compute total time rate of change of external mode kinetic energy c (gigaWatts). Then, normalize by the basin volume. c do 330 jj=2,jmtm2 fx=cst(jj)*dyt(jj)/c2dtsf do 330 i=2,imtm1 engext(1)=engext(1)-p(i,jj)*ztd(i,jj)*fx*dxt(i) 330 continue scl=cm2tom2*rho0*cm3tom3*giga do 340 ll=1,mterms engext(ll)=engext(ll)*scl/volume 340 continue c c----------------------------------------------------------------------- c Integrate previously computed integrals vertically. The total c integrals are summed over K and stored into the zeroth element c of the array and the normalized by the basin volume. c----------------------------------------------------------------------- c c Average buoyancy work per unit volume (gigaWatts). c scl=cm2tom2*rho0*cm3tom3*giga do 350 k=km,1,-1 buoyint(k)=buoyint(k)*scl buoyint(0)=buoyint(0)+buoyint(k) 350 continue buoyint(0)=buoyint(0)/basinvol c c Rate of change of internal mode kinetic energy (gigaWatts). c scl=cm2tom2*rho0*cm3tom3*giga do 370 ll=1,mterms do 360 k=km,1,-1 engint(k,ll)=engint(k,ll)*scl engint(0,ll)=engint(0,ll)+engint(k,ll) 360 continue engint(0,ll)=engint(0,ll)/basinvol 370 continue c c Term balances for the time rate of change of momentum per unit volume c (momentum units * m^3/s). c scl=cm3tom3 do 390 ll=1,mterms do 380 k=km,1,-1 termbm(k,ll,1)=termbm(k,ll,1)*scl termbm(0,ll,1)=termbm(0,ll,1)+termbm(k,ll,1) termbm(k,ll,2)=termbm(k,ll,2)*scl termbm(0,ll,2)=termbm(0,ll,2)+termbm(k,ll,2) 380 continue termbm(0,ll,1)=termbm(0,ll,1)/basinvol termbm(0,ll,2)=termbm(0,ll,2)/basinvol 390 continue c c Term balances for the time rate of change of tracer per unit volume c (tracer units * m^3/s). c scl=cm3tom3 do 410 m=1,nt do 410 ll=1,tterms do 400 k=km,1,-1 termbt(0,ll,m)=termbt(k,ll,m)*scl termbt(0,ll,m)=termbt(0,ll,m)+termbt(k,ll,m) 400 continue termbt(0,ll,m)=termbt(0,ll,m)/basinvol 410 continue c c Normalize other diagnostic by the appropriate basin volume. c avwv(1)=avwv(1)/volume avwv(2)=avwv(2)/volume avwv(3)=avwv(3)/volume smflx(1)=smflx(1)/volume smflx(2)=smflx(2)/volume do 420 m=1,nt tvar(0,m)=tvar(0,m)/volume atwv(m)=atwv(m)/volume stflx(m)=stflx(m)/volume asst(m)=asst(m)/area 420 continue c c Compute buoyancy conversion error, nonlinear exchange error, and c residual (ficticious source) terms. c buoyerr=buoyint(0)-engint(0,2)-engext(2) engleak=engint(0,3)+engint(0,4)+engint(0,5)+engint(0,6)+ * engext(3)+engext(4)+engext(5)+engext(6) plicin=engint(0,1) plicex=engext(1) do 425 jj=2,mterms if ((jj.ne.10).and.(jj.ne.11)) then plicin=plicin-engint(0,jj) plicex=plicex-engext(jj) endif 425 continue c c Compute time rate of change of tracer due to convection. c do 430 m=1,nt tconv(m)=termbt(0,1,m) do 430 jj=2,tterms tconv(m)=tconv(m)-termbt(0,jj,m) 430 continue c c Scale meridional mass transport to Sverdrups. c do 440 jj=1,jmt do 440 k=1,km tmt(jj,k)=vscle*tmt(jj,k) 440 continue c c Scale northward transport of heat to petawatts and northward c transport of salt to 10^10 cm^3/sec. c do 450 jj=1,jmt do 450 ll=1,tterms ttn(ll,jj,1)=tscle*ttn(ll,jj,1) ttn(ll,jj,2)=sscle*ttn(ll,jj,2) 450 continue c c======================================================================= c Write out globally averaged data into energy/diagnostics NetCDF file. c======================================================================= c c Average buoyancy work per unit volume. c start(1)=1 count(1)=kmp1 start(2)=tindx count(2)=1 call ncvpt(ncnrgid,bavgid,start,count,buoyint(0),rcode) if(rcode.ne.0) then write(stdout,900) 'buoy_avg' call exitus('DIAG') endif c c Time rate of change of average internal mode kinetic energy c components. c start(2)=1 count(2)=kmp1 start(3)=tindx count(3)=1 do 460 n=1,mterms start(1)=n count(1)=1 call ncvpt(ncnrgid,ikeavgid,start,count,engint(0,n),rcode) if(rcode.ne.0) then write(stdout,901) 'ke_intavg',n call exitus('DIAG') endif 460 continue c c Time rate of change of average external mode kinetic energy c components. c start(1)=1 count(1)=mterms start(2)=tindx count(2)=1 call ncvpt(ncnrgid,ekeavgid,start,count,engext,rcode) if(rcode.ne.0) then write(stdout,900) 'ke_extavg' call exitus('DIAG') endif c c Term balances for the time rate of change of tracer. c start(3)=1 count(3)=kmp1 start(4)=tindx count(4)=1 do 470 m=1,nt start(1)=m count(1)=1 do 470 n=1,tterms start(2)=n count(2)=1 call ncvpt(ncnrgid,tbtrmid,start,count,termbt(0,n,m),rcode) if(rcode.ne.0) then write(stdout,902) 't_bterms',n,m call exitus('DIAG') endif 470 continue c c Term balances for the time rate of change of momentum. c start(3)=1 count(3)=kmp1 start(4)=tindx count(4)=1 do 480 m=1,2 start(1)=m count(1)=1 do 480 n=1,mterms start(2)=n count(2)=1 call ncvpt(ncnrgid,mbtrmid,start,count,termbm(0,n,m),rcode) if(rcode.ne.0) then write(stdout,903) 'm_bterms',n,m call exitus('DIAG') endif 480 continue c c Northward transport of tracer by components. c start(3)=1 count(3)=jmt start(4)=tindx count(4)=1 do 490 m=1,nt start(1)=m count(1)=1 do 490 n=1,tterms start(2)=n count(2)=1 call ncvpt(ncnrgid,ttnid,start,count,ttn(n,1,m),rcode) if(rcode.ne.0) then write(stdout,902) 'ttn',n,m call exitus('DIAG') endif 490 continue c c Meridional mass transport. c start(1)=1 count(1)=nlev start(3)=tindx count(3)=1 do 500 jj=1,jmt start(2)=jj count(2)=1 call ncvpt(ncnrgid,tmtid,start,count,tmt(jj,1),rcode) if(rcode.ne.0) then write(stdout,904) 'tmt',jj call exitus('DIAG') endif 500 continue c c Surface tracer flux. c start(1)=1 count(1)=nt start(2)=tindx count(2)=1 call ncvpt(ncnrgid,stflxid,start,count,stflx,rcode) if(rcode.ne.0) then write(stdout,900) 'stflx' call exitus('DIAG') endif c c Surface momentum flux. c start(1)=1 count(1)=2 start(2)=tindx count(2)=1 call ncvpt(ncnrgid,smflxid,start,count,smflx,rcode) if(rcode.ne.0) then write(stdout,900) 'smflx' call exitus('DIAG') endif c c Average surface tracer. c start(1)=1 count(1)=nt start(2)=tindx count(2)=1 call ncvpt(ncnrgid,asstid,start,count,asst,rcode) if(rcode.ne.0) then write(stdout,900) 'asst' call exitus('DIAG') endif c c Average tracer within volume. c start(1)=1 count(1)=nt start(2)=tindx count(2)=1 call ncvpt(ncnrgid,atwvid,start,count,atwv,rcode) if(rcode.ne.0) then write(stdout,900) 'atwv' call exitus('DIAG') endif c c Time rate of change of tracer due to convection. c start(1)=1 count(1)=nt start(2)=tindx count(2)=1 call ncvpt(ncnrgid,tconvid,start,count,tconv,rcode) if(rcode.ne.0) then write(stdout,900) 'tconv' call exitus('DIAG') endif c c Time rate of change of tracer variance. c start(2)=1 count(2)=kmp1 start(3)=tindx count(3)=1 do 510 m=1,nt start(1)=m count(1)=1 call ncvpt(ncnrgid,tvarid,start,count,tvar(0,m),rcode) if(rcode.ne.0) then write(stdout,905) 'tvar',m call exitus('DIAG') endif 510 continue c c Average velocity within volume. c start(1)=1 count(1)=3 start(2)=tindx count(2)=1 call ncvpt(ncnrgid,avwvid,start,count,avwv,rcode) if(rcode.ne.0) then write(stdout,900) 'avwv' call exitus('DIAG') endif c c Synchronize energy/diagnostics NetCDF file to disk to allow other c processes to access data immediately after it is written. c if(j.eq.jmtm2) then call ncsnc(ncnrgid,rcode) if(rcode.ne.0) then write(stdout,906) call exitus('DIAG') endif endif c c======================================================================= c Monitor velocities by showing which come closest to the CFL c condition. c======================================================================= c write(stdout,907) itt write(stdout,908) 'U ',cflum,cflup,icflu,jcflu,kcflu write(stdout,908) 'V ',cflvm,cflvp,icflv,jcflv,kcflv write(stdout,908) 'W ',cflwtm,cflwtp,icflwt,jcflwt,kcflwt write(stdout,908) 'WU',cflwum,cflwup,icflwu,jcflwu,kcflwu write(stdout,909) c 900 format(/' DIAG - error while writing variable: ',a) 901 format(/' DIAG - error while writing variable: ',a,2x, * ' at TERM = ',i2) 902 format(/' DIAG - error while writing variable: ',a,2x, * ' at TERM = ',i2,2x,' and TRACER = ',i2) 903 format(/' DIAG - error while writing variable: ',a,2x, * ' at TERM = ',i2,2x,' and COMPONENT = ',i2) 904 format(/' DIAG - error while writing variable: ',a,2x, * ' at row J = ',i4) 905 format(/' DIAG - error while writing variable: ',a,2x, * ' at TRACER = ',i2) 906 format(/' DIAG - unable to synchronize energy NetCDF to disk.') 907 format(/,' CFL Summary for timestep = ',i10,/) 908 format(' maximum ',a2,' velocity (',g10.3,') is ',f7.2,'% the', * ' CFL at (i,j,k) = (',i3,',',i3,',',i3,')') 909 format() return end