subroutine diago
#if defined O_mom
!=======================================================================
! write out diagnostic output
!=======================================================================
implicit none
integer numdsp, jrow, i, nislsp, mdscan, is, indp, ie, js, je
integer io, m, mask, k, msk, num, n, iobadt, iobads, jte, jue
integer jwte, jwue, ll, nv, ks, ke, iv, ih, lll, maxm, mloop
integer ms, me, l, len, j, jj, iou
logical defined
real uext, vext, rnum, npt, dspcrt, rgrav, scl, reltim, rvolgk
real rvolgt, rareag, avgper, fx, tnew, tsml, tbig, plicin, plicex
real buoerr, enleak, dtconv, dtfilt, taux, tauy, erru, errv
real contu, dchg, cont, tperiod, pwatts, csalt
real tmbavg, convrt, sumsto, sumdiv, sumflx, sumdif, sumsor
real sumvol, terror, serror, cmmday, cwatts, zmau, zmsmf, zmsm
real zmat, zmstf, zmst, tconv, tfilt, rnavg, timunit, tmp
character(120) :: fname, ftbt, file_stamp, new_file_name
character(32) :: nstamp
save ftbt
data ftbt /' '/
integer ntrec, nyear, nmonth, nday, nhour, nmin, nsec
real time
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "calendar.h"
include "ctmb.h"
include "coord.h"
include "cprnts.h"
include "cregin.h"
include "diag.h"
include "docnam.h"
include "emode.h"
include "grdvar.h"
include "iounit.h"
include "levind.h"
include "mw.h"
include "scalar.h"
include "stab.h"
include "state.h"
include "switch.h"
include "tmngr.h"
include "vmixc.h"
# if defined O_embm
include "cembm.h"
# endif
# if defined O_tracer_averages
include "ctavg.h"
# endif
# if defined O_save_carbon_carbonate_chem
include "timeavgs.h"
include "diaga.h"
# endif
include "npzd.h"
# if defined O_mom_tbt
include "tbt.h"
real errt(nt), tm(imt,jmt,km)
# endif
# if defined O_diagnostic_surf_height && defined O_stream_function
character(8) :: bc_symm
! psgrad= gradient of sea surface pressure (1,2) = (u,v) component
! dsp = diagnostic sea surface pressure converted to height
! divf = divergence of sea surface pressure gradients
# if !defined O_rigid_lid_surface_pressure && !defined O_implicit_free_surface
real divf(imt,jmt)
save divf
# endif
real psgrad(imt,jmt,2), dsp(imt,jmt)
save dsp, numdsp
!-----------------------------------------------------------------------
! reconstruct the diagnostic surface pressure
!-----------------------------------------------------------------------
! initialize the surface pressure fields and averaging counter
if (first) then
numdsp = 0
do jrow=1,jmt
do i=1,imt
divf(i,jrow) = c0
dsp(i,jrow) = c0
res(i,jrow) = c0
enddo
enddo
endif
if (dspperts) then
! accumulate forcing for diagnostic surface height over the
! prescribed averaging period "dspper"
! increment the counter for averaging the divergence
numdsp = numdsp + 1
! construct surface pressure gradients
! (note: uext & vext are returned but not needed)
call calc_psgrad (psgrad, uext, vext, 1, jmt, 1, imt)
do i=1,imt
psgrad(i,jmt,1) = c0
psgrad(i,jmt,2) = c0
psgrad(i,1,1) = c0
psgrad(i,1,2) = c0
enddo
call setbcx (psgrad(1,1,1), imt, jmt)
call setbcx (psgrad(1,1,2), imt, jmt)
! compute the divergence of the sea surface pressure gradients
call spforc (psgrad, dxu, dyu, csu, h, res)
! accumulate the divergence
do jrow=1,jmt
do i=1,imt
divf(i,jrow) = divf(i,jrow) + res(i,jrow)
enddo
enddo
endif
if (dspts) then
write (stdout,'(///,40x,a,/)')
& 'D I A G O N O S T I C S U R F A C E H E I G H T'
! average the divergence, zero the guess
rnum = c1/numdsp
do jrow=1,jmt
do i=1,imt
divf(i,jrow) = rnum*divf(i,jrow)
dsp(i,jrow) = c0
enddo
enddo
! initialize coefficients for the conjugate gradient solver
call spc9pt (dxu, dyu, csu, h, cf)
npt = 9
variable = 'surfpres'
bc_symm = 't even'
! number of islands must be zero for surface pressure
! dspcrt = tolerence of solution for dsp
nislsp = 0
dspcrt = tolrsp
call congr (npt, variable, bc_symm, dsp, dsp, divf, res
&, cf
&, mxscan, mdscan, dspcrt
&, imask, iperm, jperm, iofs, nislsp, nippts
&, converged, esterr)
write (stdout,'(/a,i5,a,e14.7)')
& '=> diagnostic surface pressure scans=',mdscan,' estimated err='
&, esterr
if (.not.converged) then
write (stdout,'(/,a,/)')
& '=> Warning: convergence not reached for diag surface pressure'
converged = .true.
endif
! subtract the null space from the surface pressure
call checkerboard (dsp, map)
call border (dsp, 't even')
! remove mean from the surface pressure
call zero_level (dsp, 'surf press', map, dxt, dyt, cst)
call border (dsp, 't even')
! convert surface pressure to surface elevation and write out
rgrav = c1/(rho0*grav)
do jrow=1,jmt
do i=1,imt
dsp(i,jrow) = rgrav*dsp(i,jrow)
enddo
enddo
if (iodsp .eq. stdout .or. iodsp .lt. 0) then
is = indp (slonxy, xt, imt)
ie = indp (elonxy, xt, imt)
js = indp (slatxy, yt, jmt)
je = indp (elatxy, yt, jmt)
scl = c1
write (stdout,'(/1x,a,i5,a)')
& 'Diagnostic surf height is averaged over ',numdsp,' time steps'
write (stdout,7100) 'Diagnostic surface height (cm)'
&, itt, xt(is), xt(ie), yt(js), yt(je), scl
call matrix (dsp(1,1), imt, is, ie, -js, -je, scl)
7100 format(1x,a30,1x,'ts=',i7
&, ', lon:',f6.2,' ==> ',f6.2,', lat:',f6.2,' ==> ',f6.2
&, ', scaling=',1pg10.3)
endif
if (iodsp .ne. stdout .or. iodsp .lt. 0) then
reltim = prelyr
call getunit (io, 'diag_surf.dta'
&, 'unformatted sequential append ieee')
iotext = 'read (iodsp) imt, jmt, reltim, dspper, xt, yt'
write (io) pstamp, iotext, expnam
write (io) imt, jmt, reltim, dspper, xt, yt
iotext = 'read (iodsp) ((dsp(i,j),i=1,imt),j=1,jmt)'
write (io) pstamp, iotext, expnam
call wrufio (io, dsp, imt*jmt)
call relunit (io)
write (stdout,'(a,a,f8.3,a,a,i10,a,a)')
& ' => Diagnostic surface height (averaged '
&, 'over ',dspper,' days) written '
&, 'unformatted to file diag_surf.dta on ts=', itt, ' ',stamp
endif
! reset counter for next average. zero the averaged divergence
numdsp = 0
do jrow=1,jmt
do i=1,imt
divf(i,jrow) = c0
enddo
enddo
endif
# endif
# if defined O_tracer_averages
!-----------------------------------------------------------------------
! compute tracer averages under horizontal regions and write them
! out
!-----------------------------------------------------------------------
if (tavgts) then
! initialize sums for horizontal regions & tracer averages
do m=1,nt
sumgf(m) = c0
avggf(m) = c0
sumgt(m) = c0
avggt(m) = c0
do mask=1,nhreg
sumbt(mask,m) = c0
avgbt(mask,m) = c0
avgbf(mask,m) = c0
enddo
do k=1,km
sumgk(k,m) = c0
avggk(k,m) = c0
enddo
enddo
! compute sums for tracer averages over horizontal regions
do m=1,nt
do mask=1,nhreg
sumgf(m) = sumgf(m) + sumbf(mask,m)
do k=1,km
sumbt(mask,m) = sumbt(mask,m) + sumbk(mask,k,m)
sumgk(k,m) = sumgk(k,m) + sumbk(mask,k,m)
enddo
sumgt(m) = sumgt(m) + sumbt(mask,m)
enddo
enddo
do k=1,km
if (volgk(k) .gt. c0) then
rvolgk = c1 / volgk(k)
do m=1,nt
avggk(k,m) = sumgk(k,m) * rvolgk
do mask=1,nhreg
if (volbk(mask,k) .gt. c0) then
avgbk(mask,k,m) = sumbk(mask,k,m) / volbk(mask,k)
endif
enddo
enddo
endif
enddo
rvolgt = c1 / volgt
rareag = c1 / areag
do m=1,nt
avggt(m) = sumgt(m) * rvolgt
avggf(m) = sumgf(m) * rareag
do mask=1,nhreg
if (volbt(mask) .gt. c0) then
avgbt(mask,m) = sumbt(mask,m) / volbt(mask)
endif
if (areab(mask) .gt. c0) then
avgbf(mask,m) = sumbf(mask,m) / areab(mask)
endif
enddo
enddo
! write out regional tracer means
if (iotavg .eq. stdout .or. iotavg .lt. 0) then
write (stdout,'(//,30x,a,/)') 'T R A C E R A V E R A G E S'
do m=1,nt
write(stdout,9004) trname(m), itt, stamp
write(stdout,9001) (mask,mask=1,nhreg)
do k=1,km
write(stdout,9002) k, avggk(k,m),
& (avgbk(mask,k,m),mask=1,nhreg)
enddo
write(stdout,9003) avggt(m), (avgbt(mask,m),mask=1,nhreg)
write(stdout,9014) m, avggf(m), (avgbf(msk,m),msk=1,nhreg)
enddo
endif
if (iotavg .ne. stdout .or. iotavg .lt. 0) then
reltim = prelyr
call getunit (io, 'tracer_avg.dta'
&, 'unformatted sequential append ieee')
write (stdout,*) ' => Regional tracer averages written'
&, ' unformatted to file tracer_avg.dta on ts=',itt, ' ',stamp
iotext = 'read (iotavg) reltim, nt, nhreg, km'
write (io) pstamp, iotext, expnam
write (io) reltim, nt, nhreg, km
iotext = 'read (iotavg) ((avggk(k,n),k=1,km),n=1,nt)'
write (io) pstamp, iotext, expnam
call wrufio (io, avggk, km*nt)
iotext =
& 'read (iotavg) (((avgbk(l,k,n),l=1,nhreg),k=1,km),n=1,nt)'
write (io) pstamp, iotext, expnam
call wrufio (io, avgbk, nhreg*km*nt)
iotext = 'read (iotavg) (avggt(n),n=1,nt)'
write (io) pstamp, iotext, expnam
call wrufio (io, avggt, nt)
iotext = 'read (iotavg) (avggf(n),n=1,nt)'
write (io) pstamp, iotext, expnam
call wrufio (io, avggf, nt)
iotext = 'read (iotavg) ((avgbt(l,n),l=1,nhreg),n=1,nt)'
write (io) pstamp, iotext, expnam
call wrufio (io, avgbt, nhreg*nt)
iotext = 'read (iotavg) ((avgbf(l,n),l=1,nhreg),n=1,nt)'
write (io) pstamp, iotext, expnam
call wrufio (io, avgbf, nhreg*nt)
call relunit (io)
endif
endif
9001 format(' k',' All Regions ',9(1x,i7,5x))
9002 format(1x,i4,10(1x,e12.6))
9003 format(' AVG',10(1x,e12.6))
9004 format(/' Volume Weighted Averages for ',a12,' on ts =',i10,a33)
9014 format(/' FLX',i1,10(1x,e12.6),/)
# endif
# if defined O_time_step_monitor
!-----------------------------------------------------------------------
! print integrals for monitoring the time step
!-----------------------------------------------------------------------
if (tsiperts .and. eots) call ta_mom_tsi (1)
if (tsits .and. ntatio .gt. 0) then
call ta_mom_tsi (2)
time = year0 + accel_yr0 + (relyr - accel_yr0)*accel
call rdstmp (stamp, nyear, nmonth, nday, nhour, nmin, nsec)
nyear = time
call mkstmp (nstamp, nyear, nmonth, nday, nhour, nmin, nsec)
if (iotsi .eq. stdout .or. iotsi .lt. 0) then
write (stdout,9601) itt, nstamp, tai_ek, tai_dt, tai_ds, tai_t
&, tai_s, tai_tvar, tai_svar, tai_scan
endif
call def_tsi
call def_tsi_mom (fname)
# if defined O_save_carbon_totals
! convert from umol cm-3 to Pg (redctn contained 1e-3)
tai_cocn = (tai_dic + (tai_p + tai_z + tai_d + tai_di)
& *redctn)*tcellv*12.e-21
! convert from umol cm-2 s-1 to Pg year-1
tai_cfa2o = tai_dicflx*tcella(1)*yrlen*daylen*12.e-21
# if defined O_sed && !defined O_sed_uncoupled
! subtract out any land to ocean flux
tai_cfa2o = tai_cfa2o - tai_dicwflx*yrlen*daylen*12.e-21
# endif
# else
tai_cocn = 0.
tai_cfa2o = 0.
# endif
tai_olivint =
& tai_olivine*tcella(1)*yrlen*daylen/4
& /1.e21
tai_olivint_dis =
& tai_olivine_dis*tcella(1)*yrlen*daylen/4
& /1.e21
# if !defined O_save_time_relyear0
! make output time relative to year 1
time = time - 1.
# endif
avgper = tsiper*accel
if (avgper .le. 1e-6) avgper = 0.
# if defined O_save_time_endper
tmp = 0.
# elif defined O_save_time_startper
tmp = 1.
# else
tmp = 0.5
# endif
# if defined O_units_time_years
# if defined O_calendar_360_day
time = time - tmp*avgper/360.
# elif defined O_calendar_gregorian
time = time - tmp*avgper/365.25
# else
time = time - tmp*avgper/365.
# endif
# else
# if defined O_calendar_360_day
time = time*360. - tmp*avgper
# elif defined O_calendar_gregorian
time = time*365.25 - tmp*avgper
# else
time = time*365. - tmp*avgper
# endif
# endif
call mom_tsi_out (fname, avgper, time, nstamp, tai_ek, tai_t
&, tai_s, tai_tvar, tai_svar, tai_dt, tai_ds
&, tai_scan, tai_hflx, tai_sflx, tai_dic
&, tai_dicflx, tai_alk, tai_o2, tai_o2flx
&, tai_po4, tai_p, tai_z, tai_d, tai_no3
&, tai_di, tai_c14, tai_dc14, tai_c14flx
&, tai_cfc11, tai_cfc11flx, tai_cfc12
&, tai_cfc12flx, tai_otmax, tai_otmin
&, tai_slh, tai_sspH, tai_ssCO3, tai_ssOc
&, tai_ssOa, tai_sspCO2, tai_cocn, tai_cfa2o
&, tai_olivine,tai_olivint,tai_olivine_dis
&, tai_olivint_dis,ntrec)
call ta_mom_tsi (0)
endif
9601 format (1x,'ts=',i7, 1x, a32, ', ke=', 1pe13.6,' |dT|=',1pe13.6
&, ' |dS|=',1pe13.6,' Tbar=',1pe13.6,' Sbar=',1pe13.6
&, ' Tvar=',1pe13.6,' Svar=',1pe13.6, ' scans=',1pe13.6)
# endif
# if defined O_stability_tests
!-----------------------------------------------------------------------
! show stability and CFL conditions, reynolds and peclet numbers
!-----------------------------------------------------------------------
if (stabts) then
write (stdout,'(///20x,a/15x,a,/)')
& ' S T A B I L I T Y A N A L Y S I S'
&, '(The indicated locations are the most unstable ones)'
write (stdout,*) ' longitudinal domain: ',cflons, ' to ',cflone
write (stdout,*) ' latitudinal domain: ',cflats, ' to ',cflate
write (stdout,*) ' depth domain (m) : ',cfldps*0.01
&, ' to ',cfldpe*0.01
write (stdout,'(/60x,a/)') ' CFL summary'
write (stdout
&,'(a,g10.3,a,f7.2,a,/a,i4,a,i4,a,i3,a,3x,a,f7.2,a,f7.2,a,f7.0,a)')
& ' Local U velocity (',cflum,') is ',cflup,' % of the CFL limit'
&, ' at location: (i,j,k) = (',icflu,',',jcflu,',',kcflu,')'
&, ' (lon,lat,dpt) = (',xu(icflu),',',yu(jcflu),','
&, 0.01*zt(kcflu),')'
write (stdout
&,'(a,g10.3,a,f7.2,a,/a,i4,a,i4,a,i3,a,3x,a,f7.2,a,f7.2,a,f7.0,a)')
& ' Local V velocity (',cflvm,') is ',cflvp,' % of the CFL limit'
&, ' at location: (i,j,k) = (',icflv,',',jcflv,',',kcflv,')'
&, ' (lon,lat,dpt) = (',xu(icflv),',',yu(jcflv),','
&, 0.01*zt(kcflv),')'
write (stdout
&,'(a,g10.3,a,f7.2,a,/a,i4,a,i4,a,i3,a,3x,a,f7.2,a,f7.2,a,f7.0,a)')
& ' Local adv_vbu (',cflwum,') is ',cflwup,' % of the CFL limit'
&, ' at location: (i,j,k) = (',icflwu,',',jcflwu,',',kcflwu,')'
&, ' (lon,lat,dpt) = (',xu(icflwu),',',yu(jcflwu),','
&, 0.01*zw(kcflwu),')'
write (stdout
&,'(a,g10.3,a,f7.2,a,/a,i4,a,i4,a,i3,a,3x,a,f7.2,a,f7.2,a,f7.0,a)')
& ' Local adv_vbt (',cflwtm,') is ',cflwtp,' % of the CFL limit'
&, ' at location: (i,j,k) = (',icflwt,',',jcflwt,',',kcflwt,')'
&, ' (lon,lat,dpt) = (',xu(icflwt),',',yu(jcflwt),','
&, 0.01*zw(kcflwt),')'
fx = 100.0
if (cflup .gt. fx .or. cflvp .gt. fx .or. cflwup .gt. fx .or.
& cflwtp .gt. fx) then
write (stdout,*)
& ' => Warning. CFL exceeded... computational mode exists!'
endif
write (stdout,'(/60x,a24/)') ' Reynolds number summary'
write (stdout,10300) reynx, ireynx, jreynx, kreynx
&, xu(ireynx), yu(jreynx), 0.01*zt(kreynx)
write (stdout,10310) reynu, reynmu
write (stdout,10400) reyny, ireyny, jreyny, kreyny
&, xu(ireyny), yu(jreyny), 0.01*zt(kreyny)
write (stdout,10410) reynv, reynmv
write (stdout,10500) reynz, ireynz, jreynz, kreynz
&, xu(ireynz), yu(jreynz), 0.01*zt(kreynz)
write (stdout,10510) reynw, reynmw
write (stdout,'(/60x,a22/)') ' Peclet number summary'
write (stdout,10600) peclx, ipeclx, jpeclx, kpeclx
&, xu(ipeclx), yu(jpeclx), 0.01*zt(kpeclx)
write (stdout,10610) peclu, peclmu
write (stdout,10700) pecly, ipecly, jpecly, kpecly
&, xu(ipecly), yu(jpecly), 0.01*zt(kpecly)
write (stdout,10710) peclv, peclmv
write (stdout,10800) peclz, ipeclz, jpeclz, kpeclz
&, xu(ipeclz), yu(jpeclz), 0.01*zt(kpeclz)
write (stdout,10810) peclw, peclmw
! show ficticious tracer extremums
call getunit (iostab, 'iostab', 'formatted sequential rewind')
rewind iostab
write (stdout,'(/,10x,a/,11x,a,1pe10.3,a/)')
& 'Spurious creation of local tracer extremum (if any) follow:'
&, '(where tracer exceeds local extremum by ',tdig,'*tracer)'
do num=1,1000
read (iostab,'(i4, i4, i4, i2, 3g14.7)', end=101, err=101)
& i, k, jrow, n, tnew, tsml, tbig
write (stdout,'(1x,a,i4,a,i4,a,i4,a,i2,3(a,g14.7))')
& 't(i,k,jrow,n) = t(',i,',',k,',',jrow,',',n, ') = '
&, tnew,' : local min was ',tsml, ' : local max was ', tbig
if (num .eq. 100) then
write (stdout,'(/a/)') 'Bailing out after 100 locations...'
go to 101
endif
enddo
101 continue
call relunit (iostab)
! show T and S outside allowable bounds used for density coeffs
call getunit (iobadt, 'iobadt', 'formatted sequential rewind')
rewind iobadt
write (stdout,'(/,10x,a/)')
& 'Temperatures (if any) outside allowable ranges follow:'
do num=1,1000
read (iobadt,'(i4, i4, i4, i2, 3g14.7)', end=201, err=201)
& i, k, jrow, n, tnew, tsml, tbig
write (stdout,'(1x,a,i4,a,i4,a,i4,a,i2,3(a,g14.7))')
& 't(i,k,jrow,n) = t(',i,',',k,',',jrow,',',n, ') = '
&, tnew,' : tmin is ',tsml, ' : tmax is ', tbig
if (num .eq. 100) then
write (stdout,'(/a/)') 'Bailing out after 100 locations...'
go to 201
endif
enddo
201 continue
call relunit (iobadt)
call getunit (iobads, 'iobads', 'formatted sequential rewind')
rewind iobads
write (stdout,'(/,10x,a/)')
& 'Salinities (if any) outside allowable ranges follow:'
do num=1,1000
read (iobads,'(i4, i4, i4, i2, 3g14.7)', end=301, err=301)
& i, k, jrow, n, tnew, tsml, tbig
write (stdout,'(1x,a,i4,a,i4,a,i4,a,i2,3(a,g14.7))')
& 'Error condition: t(i,k,jrow,n) = t(',i,',',k,',',jrow
&, ',',n, ') = ', tnew,' : smin is ',tsml, ' : smax is ', tbig
if (num .eq. 100) then
write (stdout,'(/a/)') 'Bailing out after 100 locations...'
go to 301
endif
enddo
301 continue
call relunit (iobads)
write (stdout,'(/60x,a/)') ' End Stability Analysis'
endif
10300 format (1x,'Maximim zonal Reynolds number is ',1pe9.2
&, ' at location: (i,j,k) = (',i4,',',i4,',',i4,'),'
&, ' (lon,lat,dpt) = (',e9.2,',',e9.2,',',e9.2,')')
10310 format (1x,' local U =',1pe9.2, ' and mixing =',e9.2)
10400 format (1x,'Maximim meridional Reynolds number is ',1pe9.2
&, ' at location: (i,j,k) = (',i4,',',i4,',',i4,'),'
&, ' (lon,lat,dpt) = (',e9.2,',',e9.2,',',e9.2,')')
10410 format (1x,' local V =',1pe9.2, ' and mixing =',e9.2)
10500 format (1x,'Maximim vertical Reynolds number is ',1pe9.2
&, ' at location: (i,j,k) = (',i4,',',i4,',',i4,'),'
&, ' (lon,lat,dpt) = (',e9.2,',',e9.2,',',e9.2,')')
10510 format (1x,' local Wu =',1pe9.2, ' and mixing =',e9.2)
10600 format (1x,'Maximim zonal Peclet number is ',1pe9.2
&, ' at location: (i,j,k) = (',i4,',',i4,',',i4,'),'
&, ' (lon,lat,dpt) = (',e9.2,',',e9.2,',',e9.2,')')
10610 format (1x,' local U =',1pe9.2, ' and mixing =',e9.2)
10700 format (1x,'Maximim meridional Peclet number is ',1pe9.2
&, ' at location: (i,j,k) = (',i4,',',i4,',',i4,'),'
&, ' (lon,lat,dpt) = (',e9.2,',',e9.2,',',e9.2,')')
10710 format (1x,' local V =',1pe9.2, ' and mixing =',e9.2)
10800 format (1x,'Maximim vertical Peclet number is ',1pe9.2
&, ' at location: (i,j,k) = (',i4,',',i4,',',i4,'),'
&, ' (lon,lat,dpt) = (',e9.2,',',e9.2,',',e9.2,')')
10810 format (1x,' local Wt =',1pe9.2, ' and mixing =',e9.2)
# endif
# if defined O_energy_analysis
!-----------------------------------------------------------------------
! add external mode component of total work done
!-----------------------------------------------------------------------
if (glents) then
call ge3 (c2dtuv)
endif
!-----------------------------------------------------------------------
! integrate previously computed energy components vertically
!-----------------------------------------------------------------------
if (glents) then
jte = 1
jue = 1
jwte = 1
jwue = 1
do k=1,km
wtlev(k,0) = c0
wulev(k,0) = c0
enddo
do jrow=2,jmt-1
do k=km,1,-1
buoy(0,1) = buoy(0,1) + buoy(k,jrow)
wtlev(k,0) = wtlev(k,0) + wtlev(k,jrow)
wulev(k,0) = wulev(k,0) + wulev(k,jrow)
enddo
do ll=1,8
do k=km,1,-1
engint(0,ll,1) = engint(0,ll,1) + engint(k,ll,jrow)
enddo
engext(ll,1) = engext(ll,1) + engext(ll,jrow)
enddo
if (abs(tcerr(jrow)) .gt. abs(tcerr(jte))) jte = jrow
if (abs(ucerr(jrow)) .gt. abs(ucerr(jte))) jue = jrow
if (abs(wtbot(jrow)) .gt. abs(wtbot(jwte))) jwte = jrow
if (abs(wubot(jrow)) .gt. abs(wubot(jwue))) jwue = jrow
enddo
buoy(0,1) = buoy(0,1)/ucellv
do ll=1,8
engint(0,ll,1) = engint(0,ll,1)/ucellv
engext(ll,1) = engext(ll,1)/ucellv
enddo
plicin = engint(0,1,1) - engint(0,2,1) - engint(0,3,1)
& - engint(0,4,1) - engint(0,5,1) - engint(0,6,1)
plicex = engext(1,1) - engext(2,1) - engext(3,1)
& - engext(4,1) - engext(5,1) - engext(6,1)
buoerr = buoy(0,1) - engint(0,6,1) - engext(6,1)
enleak = engint(0,2,1) + engint(0,3,1) + engext(2,1)
& + engext(3,1)
if (ioglen .eq. stdout .or. ioglen .lt. 0) then
write (stdout,'(///40x,a,/)') 'E N E R G Y A N A L Y S I S'
write (stdout,9100)
& 'Globally averaged work done'
&, ' for ts =', itt, stamp, ucellv, ucella(1)
write (stdout,9101) ' time rate of change ',engint(0,1,1)
&, engint(0,1,1), engext(1,1), engext(1,1)
write (stdout,9101) ' horizontal advection', engint(0,2,1)
&, engint(0,2,1), engext(2,1), engext(2,1)
write (stdout,9101) ' vertical advection ',engint(0,3,1)
&, engint(0,3,1), engext(3,1), engext(3,1)
write (stdout,9101) ' horizontal friction ',engint(0,4,1)
&, engint(0,4,1), engext(4,1), engext(4,1)
write (stdout,9101) ' vertical friction ',engint(0,5,1)
&, engint(0,5,1), engext(5,1), engext(5,1)
write (stdout,9101) ' pressure forces ',engint(0,6,1)
&, engint(0,6,1), engext(6,1), engext(6,1)
write (stdout,9101) ' ficticious sources ',plicin
&, plicin, plicex, plicex
write (stdout,9101) ' work by wind ',engint(0,7,1)
&, engint(0,7,1), engext(7,1), engext(7,1)
write (stdout,9101) ' bottom drag ',engint(0,8,1)
&, engint(0,8,1), engext(8,1), engext(8,1)
write (stdout,9110) buoy(0,1), buoy(0,1), buoerr, buoerr
&, enleak, enleak
write (stdout,9111) tcerr(jte), itcerr(jte), jtcerr(jte)
&, ktcerr(jte)
&, ucerr(jue), iucerr(jue), jucerr(jue)
&, kucerr(jue)
write (stdout,9112) wtbot(jwte), iwtbot(jwte), jwtbot(jwte)
&, kwtbot(jwte)
&, wubot(jwue), iwubot(jwue), jwubot(jwue)
&, kwubot(jwue)
write (stdout,'(/a,a,//a,a,a,a,a,/1x,a,a)')
& 'Average vertical velocity through bottom of "T" and "u" '
&, 'cells at each level','level', ' adv_vbt err '
&, ' "T" cell area ', ' adv_vbu err ',' "U" cell area'
&, '(Note: adv_vbu err only goes to zero when non-zero values on'
&, ' boundary cells are taken into account)'
do k=1,km
if (tcella(k) .ne. c0) then
wtlev(k,0) = wtlev(k,0)/tcella(k)
else
wtlev(k,0) = c0
endif
if (ucella(k) .ne. c0) then
wulev(k,0) = wulev(k,0)/ucella(k)
else
wulev(k,0) = c0
endif
write (stdout,'(i4,4(2x,e14.7))')
& k, wtlev(k,0), tcella(k), wulev(k,0), ucella(k)
enddo
endif
if (ioglen .ne. stdout .or. ioglen .lt. 0) then
reltim = prelyr
call getunit (io, 'energy_int.dta'
&, 'unformatted sequential append ieee')
write (stdout, *)
& ' ==> Global energy integrals written unformatted'
&, ' to file energy_int.dta on ts =', itt, stamp
iotext = 'read(ioglen) reltim,(engint(i),engext(i),i=1,8)'
iotext(46:) = ',plicin,plicex,buoy,buoerr,enleak'
write (io) pstamp, iotext, expnam
write (io) reltim, (engint(0,i,1),engext(i,1),i=1,8)
&, plicin, plicex, buoy(0,1), buoerr, enleak
call relunit (io)
endif
endif
9100 format(///,1x,
&/1x,a,a,i10,a/1x,'ocean volume =',e16.9,' cm**3'
&, ', ocean surface area =',e16.9,' cm**2'//' work by:',14x
&, 'internal mode external mode'/)
9101 format(a21,2(1pe15.6, ' (',z16,' hex)'))
9110 format(/' work by buoyancy forces =',1pe14.6, ' (',z16,' hex)'/
&, ' energy conversion error =',1pe14.6, ' (',z16,' hex)'/
&, ' nonlinear error =',1pe14.6, ' (',z16,' hex)'/)
9111 format(/' max "t" cell continuity error =',1pe14.6, ' at location'
&, ' (i,jrow,k) = ','(',i4,',',i4,',',i4,')'
&, /' max "u" cell continuity error =',1pe14.6, ' at location'
&, ' (i,jrow,k) = ','(',i4,',',i4,',',i4,')')
9112 format(/' max bottom "adv_vbt" (error) =',1pe14.6
&, ' at location', ' (i,jrow,k) = ','(',i4,',',i4,',',i4,')'
&, /' max bottom "adv_vbu" (slope vel) =',1pe14.6
&, ' at location',' (i,jrow,k) = ','(',i4,',',i4,',',i4,')')
# endif
# if defined O_term_balances
!-----------------------------------------------------------------------
! add the external mode part of d/dt into the momentum balance,
! the external mode part of the implicit coriolis term, and the
! surface pressure gradients into the specified volumes
!-----------------------------------------------------------------------
if (trmbts) then
call utb3
!-----------------------------------------------------------------------
! integrate previously computed term balance components vertically
!-----------------------------------------------------------------------
do n=0,numreg
if (n .gt. 0) then
nv = (n-1)/nhreg + 1
ks = llvreg(nv,1)
ke = llvreg(nv,2)
do ll=1,17
do k=ke,ks,-1
termbm(0,ll,1,n) = termbm(0,ll,1,n) + termbm(k,ll,1,n)
termbm(0,ll,2,n) = termbm(0,ll,2,n) + termbm(k,ll,2,n)
enddo
enddo
else
ks = 1
ke = km
endif
do m=1,nt
do ll=1,15
do k=ke,ks,-1
termbt(0,ll,m,n) = termbt(0,ll,m,n) + termbt(k,ll,m,n)
enddo
enddo
! construct change due to convection and filtering
dtconv = termbt(0,10,m,n) - termbt(0,9,m,n)
dtfilt = termbt(0,1,m,n) - termbt(0,10,m,n)
termbt(0,9,m,n) = dtconv
termbt(0,10,m,n) = dtfilt
enddo
enddo
! normalize integrals by appropriate volume (or area)
do n=0,numreg
do m=1,nt
if (n .le. nhreg) then
stflx(m,n) = stflx(m,n)*rareat(n)
asst(m,n) = asst(m,n)*rareat(n)
endif
do ll=1,15
termbt(0,ll,m,n) = termbt(0,ll,m,n)*rvolt(n)
enddo
enddo
enddo
do n=0,numreg
if (n .le. nhreg) then
smflx(1,0) = smflx(1,0) + smflx(1,n)
smflx(2,0) = smflx(1,0) + smflx(2,n)
smflx(1,n) = smflx(1,n)*rareau(n)
smflx(2,n) = smflx(2,n)*rareau(n)
endif
if (n .gt. 0) then
avgw(n) = avgw(n)*rvolu(n)
do ll=1,17
termbm(0,ll,1,n) = termbm(0,ll,1,n)*rvolu(n)
termbm(0,ll,2,n) = termbm(0,ll,2,n)*rvolu(n)
enddo
endif
enddo
smflx(1,0) = smflx(1,0)*rareau(0)
smflx(2,0) = smflx(2,0)*rareau(0)
if (iotrmb .eq. stdout .or. iotrmb .lt. 0) then
write (stdout,'(///,40x,a,/)') 'T E R M B A L A N C E S'
n = 0
taux = smflx(1,n)
tauy = smflx(2,n)
write (stdout,10106)
& 'All regions added together for ts ='
&, itt, stamp, volt(n), areat(n)
&, volu(n), areau(n)
write (stdout,10104)
write (stdout,10098) n, ' smf(1) = ', taux,' dynes/cm**2 '
write (stdout,10098) n, ' smf(2) = ', tauy,' dynes/cm**2 '
do m=1,nt
write (stdout,10098) n, ustf(m,1), stflx(m,n), ustf(m,2)
enddo
write (stdout,10098) n, ' tot heat = '
&, (termbt(0,15,1,n)*volt(n))
&, ' deg C * cm**3 '
write (stdout,10098) n, ' sst = ',asst(1,n)
&, ' deg C '
do n=0,numreg
if (n .eq. 1) then
write (stdout,10050)
& 'Regional averaged Momentum & Tracer Term Balances for ts = '
&, itt, stamp
endif
iv = 0
if (n .gt. 0) then
iv = (n-1)/nhreg + 1
ih = n - (iv-1)*nhreg
write (stdout,10100)
& 'Momentum terms averaged over region #'
&, n, ': ', hregnm(ih), vregnm(iv), volu(n), areau(n)
write (stdout,10104)
write (stdout,10101) n,' dU/dt = ', termbm(0,1,1,n)
&, ' dV/dt = ', termbm(0,1,2,n)
write (stdout,10101) n,' -Px = ', termbm(0,2,1,n)
&, ' -Py = ', termbm(0,2,2,n)
write (stdout,10101) n,' -surf Px= ', termbm(0,12,1,n)
&, ' -surf Py= ', termbm(0,12,2,n)
write (stdout,10101) n,' -(UU)x = ', termbm(0,3,1,n)
&, ' -(UV)x = ', termbm(0,3,2,n)
write (stdout,10101) n,' -(VU)y = ', termbm(0,4,1,n)
&, ' -(VV)y = ', termbm(0,4,2,n)
write (stdout,10101) n,' -(WU)z = ', termbm(0,5,1,n)
&, ' -(WV)z = ', termbm(0,5,2,n)
write (stdout,10101) n,'ADV_Umet = ', termbm(0,13,1,n)
&, '-ADV_Vmet= ', termbm(0,13,2,n)
write (stdout,10101) n,' DIFF_Ux= ', termbm(0,6,1,n)
&, ' DIFF_Vx= ', termbm(0,6,2,n)
write (stdout,10101) n,' DIFF_Uy= ', termbm(0,7,1,n)
&, ' DIFF_Vy= ', termbm(0,7,2,n)
write (stdout,10101) n,' DIFF_Uz= ', termbm(0,8,1,n)
&, ' DIFF_Vz= ', termbm(0,8,2,n)
write (stdout,10101) n,'DIFF_Umet= ', termbm(0,9,1,n)
&, 'DIFF_Vmet= ', termbm(0,9,2,n)
write (stdout,10101) n,' fV = ', termbm(0,10,1,n)
&, ' -fU = ', termbm(0,10,2,n)
write (stdout,10101) n,' source = ', termbm(0,11,1,n)
&, ' source = ', termbm(0,11,2,n)
erru = c0
errv = c0
do lll=2,13
erru = erru + termbm(0,lll,1,n)
errv = errv + termbm(0,lll,2,n)
enddo
write (stdout,10101) n,' error = ', termbm(0,1,1,n)-erru
&, ' error = ', termbm(0,1,2,n)-errv
write (stdout,*) ' '
write (stdout,10101) n,' -U(U)x = ', termbm(0,14,1,n)
&, ' -U(V)x = ', termbm(0,14,2,n)
write (stdout,10101) n,' -V(U)y = ', termbm(0,15,1,n)
&, ' -V(V)y = ', termbm(0,15,2,n)
write (stdout,10101) n,' -W(U)z = ', termbm(0,16,1,n)
&, ' -W(V)z = ', termbm(0,16,2,n)
! mass conservation within volume: Ux + Vy + Wz
contu = (termbm(0,3,1,n) + termbm(0,4,1,n) +
& termbm(0,5,1,n)) - (termbm(0,14,1,n) +
& termbm(0,15,1,n) + termbm(0,16,1,n))
write (stdout,10101) n,' mass err= ', contu
write (stdout,10101) n,' ubar = ', termbm(0,17,1,n)
&, ' vbar = ', termbm(0,17,2,n)
&, ' wbar = ', avgw(n)
endif
if (iv .eq. 1) then
write (stdout,10101) n,' surf Uz= ', smflx(1,n)
&, ' surf Vz= ', smflx(2,n)
endif
if (n .eq. 0) then
write (stdout,10051)
& 'Global averaged (all basins) Tracer Term Balances for ts = '
&, itt, stamp
else
write (stdout,10100)
& 'Tracer terms averaged over region #'
&, n, ': ', hregnm(ih), vregnm(iv), volt(n), areat(n)
endif
do m=1,nt
dchg = termbt(0,2,m,n) + termbt(0,3,m,n) +
& termbt(0,4,m,n) + termbt(0,5,m,n) +
& termbt(0,6,m,n) + termbt(0,7,m,n) +
& termbt(0,8,m,n) + termbt(0,9,m,n) +
& termbt(0,10,m,n)
terr(m) = termbt(0,1,m,n) - dchg
enddo
maxm = (nt-1)/7 + 1
do mloop=1,maxm
ms = (mloop-1)*7 + 1
me = min(ms + 6,nt)
write (stdout,10103) (trname(m),m=ms,me)
write (stdout,10102) n,' dT/dt = ', (termbt(0,1,m,n)
&, m=ms,me)
write (stdout,10102) n,' -(UT)x = ', (termbt(0,2,m,n)
&, m=ms,me)
write (stdout,10102) n,' -(VT)y = ', (termbt(0,3,m,n)
&, m=ms,me)
write (stdout,10102) n,' -(WT)z = ', (termbt(0,4,m,n)
&, m=ms,me)
write (stdout,10102) n,' DIFF_Tx= ', (termbt(0,5,m,n)
&, m=ms,me)
write (stdout,10102) n,' DIFF_Ty= ', (termbt(0,6,m,n)
&, m=ms,me)
write (stdout,10102) n,' DIFF_Tz= ', (termbt(0,7,m,n)
&, m=ms,me)
write (stdout,10102) n,' source = ', (termbt(0,8,m,n)
&, m=ms,me)
write (stdout,10102) n,' convct = ', (termbt(0,9,m,n)
&, m=ms,me)
write (stdout,10102) n,' filter = ', (termbt(0,10,m,n)
&, m=ms,me)
write (stdout,10102) n,' error = ', (terr(m)
&, m=ms,me)
write (stdout,*) ' '
write (stdout,10102) n,' -U(T)x = ', (termbt(0,11,m,n)
&, m=ms,me)
write (stdout,10102) n,' -V(T)y = ', (termbt(0,12,m,n)
&, m=ms,me)
write (stdout,10102) n,' -W(T)z = ', (termbt(0,13,m,n)
&, m=ms,me)
! mass conservation within volume: Ux + Vy + Wz
cont = (termbt(0,2,1,n) + termbt(0,3,1,n) +
& termbt(0,4,1,n)) - (termbt(0,11,1,n) +
& termbt(0,12,1,n) + termbt(0,13,1,n))
write (stdout,10102) n,' mass err= ', cont
write (stdout,10102) n,' chg var= ', (termbt(0,14,m,n)
&, m=ms,me)
write (stdout,10102) n,' tbar = ', (termbt(0,15,m,n)
&, m=ms,me)
if (iv .eq. 1) then
write (stdout,10102) n,' surflx = ', (stflx(m,n)
&, m=ms,me)
taux = smflx(1,n)
tauy = smflx(2,n)
write (stdout,10105) ' Regionally averaged quantities:'
write (stdout,'(1x)')
write (stdout,10098) n
&, ' smf(1) = ', taux,' dynes/cm**2 '
write (stdout,10098) n
&, ' smf(2) = ', tauy,' dynes/cm**2 '
do m=1,nt
write (stdout,10098) n
&, ustf(m,1), stflx(m,n), ustf(m,2)
enddo
if (ms .eq. 1) then
write (stdout,10098) n, ' tot heat = '
&, (termbt(0,15,1,n)*volt(n)),' deg C * cm**3 '
write (stdout,10098) n
&, ' sst = ', asst(1,n),' deg C '
endif
endif
enddo
enddo
endif
if (iotrmb .ne. stdout .or. iotrmb .lt. 0) then
reltim = prelyr
tperiod = 0.0
call getunit (io, 'term_bal.dta'
&, 'unformatted sequential append ieee')
write (stdout, *)
& ' ==> Term balances written unformatted to file term_bal.dta'
&, ' for ts =', itt, ' ',stamp
iotext =' read (iotrmb) reltim, nt, numreg, nhreg, km'
iotext(45:) =', ((ustf*15(n,i),n=1,nt),i=1,2)'
write (io) pstamp, iotext, expnam
write (io) reltim, nt, numreg, nhreg, km, ustf
iotext ='read (iotrmb) (hregnm*40(n),n=1,nhreg)'
iotext(39:)=', (vregnm*20(n),n=1,nvreg)'
write (io) pstamp, iotext, expnam
write (io) hregnm, vregnm
iotext ='read (iotrmb) (trname*12(n),n=1,nt)'
write (io) pstamp, iotext, expnam
write (io) trname
iotext ='read (iotrmb) (volu(l),l=0,numreg)'
write (io) pstamp, iotext, expnam
call wrufio (io, volu, numreg+1)
iotext ='read (iotrmb) (volt(l),l=0,numreg)'
write (io) pstamp, iotext, expnam
call wrufio (io, volt, numreg+1)
iotext ='read (iotrmb) (areau(l),l=0,numreg)'
write (io) pstamp, iotext, expnam
call wrufio (io, areau, numreg+1)
iotext ='read (iotrmb) (areat(l),l=0,numreg)'
write (io) pstamp, iotext, expnam
call wrufio (io, areat, numreg+1)
iotext ='read (iotrmb) ((smflx(i,l),i=1,2),l=0,nhreg)'
write (io) pstamp, iotext, expnam
call wrufio (io, smflx, 2*(nhreg+1))
iotext ='read (iotrmb) ((stflx(n,l),n=1,nt),l=0,nhreg)'
write (io) pstamp, iotext, expnam
call wrufio (io, stflx, nt*(nhreg+1))
iotext ='read (iotrmb) ((asst(n,l),n=1,nt),l=0,nhreg)'
write (io) pstamp, iotext, expnam
call wrufio (io, asst, nt*(nhreg+1))
iotext ='read (iotrmb) (avgw(l),l=1,numreg)'
write (io) pstamp, iotext, expnam
call wrufio (io, avgw, numreg)
iotext = 'read (iotrmb) (((termbm'
iotext(24:) = '(i,n,l),i=1,17),n=1,2),l=1,numreg)'
write (io) pstamp, iotext, expnam
write (io) (((termbm(0,i,n,l),i=1,17),n=1,2),l=1,numreg)
iotext = 'read (iotrmb) (((termbt'
iotext(24:) = '(i,n,l),i=1,15),n=1,nt),l=0,numreg)'
write (io) pstamp, iotext, expnam
write (io) (((termbt(0,i,n,l),i=1,15),n=1,nt),l=0,numreg)
call relunit (io)
endif
endif
10050 format(///,1x,/1x,a,i10,a/)
10051 format(///,1x,a,i10,a/)
10098 format(1x,'(',i4,')',a12,(1pe16.8,a15))
10100 format(/1x,a,i5,a,a,a/1x,'ocean volume =',e16.9,' cm**3'
&, ', horizontal ocean area =',e16.9,' cm**2'/)
10101 format(1x,'(',i4,')',a11,1pe15.7, 2x, a11,1pe15.7, 2x, a11
&, 1pe15.7)
10102 format(1x,'(',i4,')',a11,7(1pe15.7))
10103 format(' Region',11x,8a15)
10104 format(' Region')
10105 format(/8x,a32)
10106 format(///,1x,
&/1x,a,i10,a/1x,'ocean "t" volume =',e16.9,' cm**3 '
&, ', ocean "t" surface area =',e16.9,' cm**2 '
&/ 1x,'ocean "u" volume =',e16.9,' cm**3, '
&, ' ocean "u" surface area =',e16.9,' cm**2'/)
# endif
# if defined O_gyre_components
!-----------------------------------------------------------------------
! write out the gyre_components
!-----------------------------------------------------------------------
! convert heat transport to petawatts,
! salt transport to 10**10 cm**3/sec
if (gyrets) then
pwatts = 4.186e-15
csalt = 1.e-10
do jrow=1,jmt
do ll=1,8
ttn(ll,jrow,1)=ttn(ll,jrow,1)*pwatts
ttn(ll,jrow,2)=ttn(ll,jrow,2)*csalt
enddo
enddo
if (iogyre .eq. stdout .or. iogyre .lt. 0) then
write (stdout,'(///,40x,a,/)') 'G Y R E C O M P O N E N T S'
write (stdout,8195)
do jrow=2,jmtm2
l = jmt - jrow
write (stdout,8196) l, (ttn(i,l,1),i=1,8)
&, (ttn(i,l,2),i=1,8)
enddo
do m=0,nhreg
do jrow=1,jmt
do ll=6,8
ttn2(ll,jrow,1,m) = ttn2(ll,jrow,1,m)*4.186e-15
ttn2(ll,jrow,2,m) = ttn2(ll,jrow,2,m)*1.e-10
enddo
enddo
enddo
write (stdout,82001)
do m=0,nhreg
if (m .ne. 0) then
write (stdout,8201) hregnm(m)
else
write (stdout,8202)
endif
write (stdout,8203)
do jrow=1,jmt
write (stdout,8204) jrow,(ttn2(i,jrow,1,m),i=6,8)
enddo
enddo
write (stdout,8205)
do m=0,nhreg
if (m .ne. 0) then
write (stdout,8201) hregnm(m)
else
write (stdout,8202)
endif
write (stdout,8203)
do jrow=1,jmt
write (stdout,8204) jrow,(ttn2(i,jrow,2,m),i=6,8)
enddo
enddo
endif
if (iogyre .ne. stdout .or. iogyre .lt. 0) then
reltim = prelyr
call getunit (io, 'gyre_comp.dta'
&, 'unformatted sequential append ieee')
write (stdout,*) ' => Gyre components written'
&, ' unformatted to file gyre_comp.dta on ts=',itt, ' ',stamp
iotext = 'read (iogyre) jmt, ntmin2, reltim'
write (io) pstamp, iotext, expnam
write (io) jmt, ntmin2, reltim
iotext =
& 'read (iogyre) (((ttn(l,j,n),l=1,8),j=1,jmt),n=1,ntmin2)'
write (io) pstamp, iotext, expnam
call wrufio (io, ttn, 8*jmt*ntmin2)
iotext =
& 'read (iogyre) (((ttn2(6..8,1..jmt,1..2,0..nhreg)'
write (io) pstamp, iotext, expnam
len = 3*jmt*nt*(1+nhreg)
call wrufio (io, ttn2, len)
call relunit (io)
endif
endif
8195 format(//,30x,'Gyre Components'/
&,6x,'northward transport of heat (x10**15 watts)'
& ,22x,'northward transport of salt (x10**10 cm**3/sec)',/,
& 3x,2(3x,'x mean x eddy z mean z eddy ekman tot adv ',
& 'diffus total'))
8196 format(i4,8f8.3,1x,8f8.3)
82001 format (/,6x,'northward transport of heat (x10**15 watts)',/)
8201 format (/,22x,a40,/)
8202 format (/,22x,' Global ',/)
8203 format (8x,'total advection',2x,'total diffusion',2x,'total')
8204 format (2x,i3,3x,3(e12.6,5x))
8205 format (/,6x,'northward transport of salt (x10**10 cm**3/sec)',/)
# endif
# if defined O_meridional_overturning
!-----------------------------------------------------------------------
! write out the meridional overturning (mass)
!-----------------------------------------------------------------------
if (vmsfts) then
if (iovmsf .eq. stdout .or. iovmsf .lt. 0) then
write (stdout,'(///,40x,a,/)')
& 'M E R I D I O N A L O V E R T U R N I N G'
scl = 1.e12
write (stdout,8194)
js = indp (slatxy, yt, jmt)
je = indp (elatxy, yt, jmt)
call matrix (vmsf, jmt, js, je, 1, km, scl)
endif
if (iovmsf .ne. stdout .or. iovmsf .lt. 0) then
reltim = prelyr
call getunit (io, 'overturn.dta'
&, 'unformatted sequential append ieee')
write (stdout,*) ' => Meridional overturning of mass'
&, ' written unformatted to file overturn.dta on ts=',itt
&, ' ',stamp
iotext = 'read (iovmsf) jmt, km, reltim'
write (io) pstamp, iotext, expnam
write (io) jmt, km, reltim
iotext = 'read (iovmsf) (zw(k),k=1,km)'
write (io) pstamp, iotext, expnam
call wrufio (io, zw, km)
iotext = 'read (iovmsf) (yu(j),j=1,jmt)'
write (io) pstamp, iotext, expnam
call wrufio (io, yu, jmt)
iotext = 'read (iovmsf) ((vmsf(j,k),j=1,jmt),k=1,km)'
write (io) pstamp, iotext, expnam
call wrufio (io, vmsf, jmt*km)
call relunit (io)
endif
endif
8194 format(/,' meridional overturning stream function (sverdrups)')
# endif
# if defined O_show_external_mode
# if defined O_rigid_lid_surface_pressure || defined O_implicit_free_surface
!-----------------------------------------------------------------------
! write out the external mode velocities and surface pressure
!-----------------------------------------------------------------------
if (extts .and. eots) then
if (ioext .eq. stdout .or. ioext .lt. 0) then
write (stdout,'(///,40x,a,/)') 'E X T E R N A L M O D E'
is = indp (slonxy, xt, imt)
ie = indp (elonxy, xt, imt)
js = indp (slatxy, yt, jmt)
je = indp (elatxy, yt, jmt)
scl=1.0
write (stdout,8000) ' ubar '
&, itt, xt(is), xt(ie), yt(js), yt(je), scl
call matrix (ubar(1,1,1), imt, is, ie, -js, -je, scl)
write (stdout,8000) ' vbar '
&, itt, xt(is), xt(ie), yt(js), yt(je), scl
call matrix (ubar(1,1,2), imt, is, ie, -js, -je, scl)
scl = 0.0
write (stdout,8000) ' div ps'
&, itt, xt(is), xt(ie), yt(js), yt(je), scl
call matrix (divf(1,1), imt, is, ie, -js, -je, scl)
scl = grav
! write (stdout,8000) ' surface press(gm/cm/sec**2)'
write (stdout,8000)
& 'sea surface height (surface_press/grav) in cm'
&, itt, xt(is), xt(ie), yt(js), yt(je), scl
call matrix (ps(1,1,1), imt, is, ie, -js, -je, scl)
endif
if (ioext .ne. stdout .or. ioext .lt. 0) then
reltim = relyr
# if defined O_rigid_lid_surface_pressure
call getunit (io, 'surf_press.dta'
&, 'unformatted sequential append ieee')
# endif
# if defined O_implicit_free_surface
call getunit (io, 'ifree_surf.dta'
&, 'unformatted sequential append ieee')
# endif
iotext = 'read (ioext) imt, jmt, reltim'
write (io) stamp, iotext, expnam
write (io) imt, jmt, reltim
iotext = 'read (ioext) (xt(i),i=1,imt)'
write (io) stamp, iotext, expnam
call wrufio (io, xt, imt)
iotext = 'read (ioext) (yt(j),j=1,jmt)'
write (io) stamp, iotext, expnam
call wrufio (io, yt, jmt)
iotext = 'read (ioext) ((ps(i,j),i=1,imt),j=1,jmt)'
write (io) stamp, iotext, expnam
call wrufio (io, ps(1,1,1), imt*jmt)
write (stdout,*) ' => Surface pressure written unformatted'
&, ' to unit ',io, ' on ts=',itt, stamp
# if defined O_rigid_lid_surface_pressure
&, ' to file surf_press.dta on ts=',itt, stamp
# endif
# if defined O_implicit_free_surface
&, ' to file ifree_surf.dta on ts=',itt, stamp
# endif
call relunit (io)
endif
endif
# endif
# if defined O_stream_function
!-----------------------------------------------------------------------
! write out the stream function
!-----------------------------------------------------------------------
if (extts .and. eots) then
if (ioext .eq. stdout .or. ioext .lt. 0) then
write (stdout,'(///,40x,a,/)') 'E X T E R N A L M O D E'
is = indp (slonxy, xt, imt)
ie = indp (elonxy, xt, imt)
js = indp (slatxy, yt, jmt)
je = indp (elatxy, yt, jmt)
scl=1.e12
write (stdout,8000) ' stream function (sverdrups)'
&, itt, xt(is), xt(ie), yt(js), yt(je), scl
call matrix (psi(1,1,1), imt, is, ie, -js, -je, scl)
endif
if (ioext .ne. stdout .or. ioext .lt. 0) then
reltim = relyr
call getunit (io, 'psi.dta'
&, 'unformatted sequential append ieee')
iotext = 'read (ioext) imt, jmt, reltim'
write (io) stamp, iotext, expnam
write (io) imt, jmt, reltim
iotext = 'read (ioext) (xt(i),i=1,imt)'
write (io) stamp, iotext, expnam
call wrufio (io, xt, imt)
iotext = 'read (ioext) (yt(j),j=1,jmt)'
write (io) stamp, iotext, expnam
call wrufio (io, yt, jmt)
iotext = 'read (ioext) ((psi(i,j),i=1,imt),j=1,jmt)'
write (io) stamp, iotext, expnam
call wrufio (io, psi(1,1,1), imt*jmt)
write (stdout,*) ' => Stream function written unformatted'
&, ' to file psi.dta on ts=',itt, stamp
call relunit (io)
endif
endif
# endif
8000 format(1x,a,1x,'ts=',i7
&,', lon:',f6.2,' ==> ',f6.2,', lat:',f6.2,' ==> ',f6.2
&,', scaling=',1pg10.3)
# endif
# if defined O_meridional_tracer_budget
!-----------------------------------------------------------------------
! construct and write out the averaged components for the
! meridional tracer balance
!-----------------------------------------------------------------------
if (tmbts) then
! quantities have been weighted by volume at each cell and
! summed over all ocean longitudes and depths for each latitude.
! average quantities over "numtmb" days and do units conversion
! to petawatts for temperature and giga cm**3/sec for salt
tmbavg = numtmb * dtts
rnum = c1 / float(numtmb)
do mask=1,ntmbb
do jrow=1,jmt
smdvol(jrow,mask) = rnum*smdvol(jrow,mask)
enddo
enddo
do mask=1,ntmbb
do n=1,nt
if (n .eq. 1) then
convrt = rnum * 4.186e-15
elseif (n .eq. 2) then
convrt = rnum * 1.e-9
else
convrt = c1
endif
do jrow=1,jmt
tstor(jrow,n,mask) = convrt*tstor(jrow,n,mask)
tdiv(jrow,n,mask) = convrt*tdiv(jrow,n,mask)
tflux(jrow,n,mask) = convrt*tflux(jrow,n,mask)
tdif(jrow,n,mask) = convrt*tdif(jrow,n,mask)
tsorc(jrow,n,mask) = convrt*tsorc(jrow,n,mask)
enddo
enddo
enddo
! accumulate terms for all basins and store in basin # 0
do jrow=1,jmt
smdvol(jrow,0) = c0
do mask=1,ntmbb
smdvol(jrow,0) = smdvol(jrow,0) + smdvol(jrow,mask)
enddo
enddo
do n=1,nt
do jrow=1,jmt
tstor(jrow,n,0) = c0
tdiv(jrow,n,0) = c0
tflux(jrow,n,0) = c0
tdif(jrow,n,0) = c0
tsorc(jrow,n,0) = c0
enddo
enddo
do mask=1,ntmbb
do n=1,nt
do jrow=1,jmt
tstor(jrow,n,0) = tstor(jrow,n,0) + tstor(jrow,n,mask)
tdiv(jrow,n,0) = tdiv(jrow,n,0) + tdiv(jrow,n,mask)
tflux(jrow,n,0) = tflux(jrow,n,0) + tflux(jrow,n,mask)
tdif(jrow,n,0) = tdif(jrow,n,0) + tdif(jrow,n,mask)
tsorc(jrow,n,0) = tsorc(jrow,n,0) + tsorc(jrow,n,mask)
enddo
enddo
enddo
! write out the results
if (iotmb .eq. stdout .or. iotmb .lt. 0) then
write (stdout,'(///,20x,a,/)')
& 'M E R I D I O N A L T R A C E R B U D G E T'
do n=1,nt
! construct sum in latitude for all basins
sumsto = c0
sumdiv = c0
sumflx = c0
sumdif = c0
sumsor = c0
sumvol = c0
do jrow=1,jmt
sumsto = sumsto + tstor(jrow,n,0)
sumdiv = sumdiv + tdiv(jrow,n,0)
sumflx = sumflx + tflux(jrow,n,0)
sumdif = sumdif + tdif(jrow,n,0)
sumsor = sumsor + tsorc(jrow,n,0)
sumvol = sumvol + smdvol(jrow,0)
enddo
write (stdout,*) ' '
if (n .eq. 1) then
write (stdout,*) '(Tracer # 1 results are in petawatts)'
elseif (n .eq. 2) then
write (stdout,*)
& '(Tracer # 2 results are in giga cm**3/sec)'
else
write (stdout,*) ' '
endif
write (stdout,8050) n, tmbavg*secday, stamp
do jrow=jmt,1,-1
terror = tstor(jrow,n,0) - (tdiv(jrow,n,0) +
& tflux(jrow,n,0) + tdif(jrow,n,0) + tsorc(jrow,n,0))
write (stdout,8100) jrow, yt(jrow), tstor(jrow,n,0)
&, tdiv(jrow,n,0), tflux(jrow,n,0), tdif(jrow,n,0)
&, tsorc(jrow,n,0), terror, smdvol(jrow,0)
enddo
serror = sumsto - (sumdiv + sumflx + sumdif + sumsor)
write (stdout,8200) sumsto, sumdiv, sumflx, sumdif
&, sumsor, serror, sumvol
enddo
endif
if (iotmb .ne. stdout .or. iotmb .lt. 0) then
reltim = relyr
call getunit (io, 'tracer_bud.dta'
&, 'unformatted sequential append ieee')
write (stdout,*) ' => Meridional Tracer Budget written '
&, ' unformatted to file tracer_bud.dta on ts=',itt, ' ',stamp
iotext = 'read (iotmb) imt, jmt, nt, ntmbb, reltim, avgper'
write (io) stamp, iotext, expnam
write (io) imt, jmt, nt, ntmbb, reltim, tmbavg*secday
iotext = 'read (iotmb) (yt(j),j=1,jmt), (dyt(j),j=1,jmt)'
write (io) stamp, iotext, expnam
write (io) (yt(j),j=1,jmt), (dyt(j),j=1,jmt)
iotext =
& 'read (iotmb) (((tstor(j,n,l),j=1,jmt),n=1,nt),l=0,ntmbb)'
write (io) stamp, iotext, expnam
call wrufio (io, tstor, jmt*nt*(ntmbb+1))
iotext =
& 'read (iotmb) (((tdiv(j,n,l),j=1,jmt),n=1,nt),l=0,ntmbb)'
write (io) stamp, iotext, expnam
call wrufio (io, tdiv, jmt*nt*(ntmbb+1))
iotext =
& 'read (iotmb) (((tflux(j,n,l),j=1,jmt),n=1,nt),l=0,ntmbb)'
write (io) stamp, iotext, expnam
call wrufio (io, tflux, jmt*nt*(ntmbb+1))
iotext =
& 'read (iotmb) (((tdif(j,n,l),j=1,jmt),n=1,nt),l=0,ntmbb)'
write (io) stamp, iotext, expnam
call wrufio (io, tdif, jmt*nt*(ntmbb+1))
iotext =
& 'read (iotmb) (((tsorc(j,n,l),j=1,jmt),n=1,nt),l=0,ntmbb)'
write (io) stamp, iotext, expnam
call wrufio (io, tsorc, jmt*nt*(ntmbb+1))
iotext =
& 'read (iotmb) ((smdvol(j,l),j=1,jmt),l=0,ntmbb)'
write (io) stamp, iotext, expnam
call wrufio (io, smdvol, jmt*(ntmbb+1))
call relunit (io)
endif
! zero quantities in preparation for next average
numtmb = 0
do mask=0,ntmbb
do jrow=1,jmt
smdvol(jrow,mask) = c0
enddo
enddo
do mask=0,ntmbb
do m=1,nt
do jrow=1,jmt
tstor(jrow,m,mask) = c0
tdiv(jrow,m,mask) = c0
tflux(jrow,m,mask) = c0
tdif(jrow,m,mask) = c0
tsorc(jrow,m,mask) = c0
enddo
enddo
enddo
endif
8050 format (/1x,'Meridional Tracer Balance for tracer #',i2
&,' averaged over depth and longitude',/1x
&,' and averaged over a period of ', f10.2
&,' days ending on ', a32//
&,' jrow lat (T)t (VT)y DIFF_Tz DIFF_Ty'
&,' source error volume'/)
8100 format (1x,i4,2x,f6.2,1x,7(1pe10.3,1x))
8200 format (1x,'global sum =',1x,7(1pe10.3,1x))
# endif
# if defined O_show_zonal_mean_of_sbc
!-----------------------------------------------------------------------
! print the zonal mean boundary conditions and related items
!-----------------------------------------------------------------------
if (zmbcts) then
cmmday = 10.0*daylen
cwatts = 4.186e4
do jrow=1,jmt
if (zmau(jrow) .ne. c0) then
zmsmf(jrow,1) = zmsmf(jrow,1) / zmau(jrow)
zmsmf(jrow,2) = zmsmf(jrow,2) / zmau(jrow)
zmsm(jrow,1) = zmsm(jrow,1) / zmau(jrow)
zmsm(jrow,2) = zmsm(jrow,2) / zmau(jrow)
endif
if (zmat(jrow) .ne. c0) then
zmstf(jrow,1) = cwatts*zmstf(jrow,1)
zmstf(jrow,2) = cmmday*zmstf(jrow,2)
zmst(jrow,2) = 1.0e3*zmst(jrow,2)
do m=1,nt
zmstf(jrow,m) = zmstf(jrow,m) / zmat(jrow)
zmst(jrow,m) = zmst(jrow,m) / zmat(jrow)
enddo
endif
enddo
if (iozmbc .eq. stdout .or. iozmbc .lt. 0) then
write (stdout,'(///,27x,a,/)')
& 'Z O N A L M E A N O F S U R F A C E B. C.'
write (stdout,10200)
do jj=1,jmt
jrow = jmt + 1 - jj
write (stdout,10201) jrow, yt(jrow), zmsmf(jrow,1)
&, zmsmf(jrow,2), zmstf(jrow,1), zmstf(jrow,2)
&, zmsm(jrow,1), zmsm(jrow,2), zmst(jrow,1), zmst(jrow,2)
enddo
endif
if (iozmbc .ne. stdout .or. iozmbc .lt. 0) then
reltim = prelyr
call getunit (io, 'zmean_sbc.dta'
&, 'unformatted sequential append ieee')
write (stdout,*) ' => Zonal mean surface bc written'
&, ' unformatted to file zmean_sbc.dta on ts=',itt, ' ',stamp
iotext = 'read (iozmbc) jmt, nt, reltim'
write (io) pstamp, iotext, expnam
write (io) jmt, nt, reltim
iotext = 'read (iozmbc) (yt(j),j=1,jmt)'
write (io) pstamp, iotext, expnam
call wrufio (io, yt, jmt)
iotext = 'read (iozmbc) (yu(j),j=1,jmt)'
write (io) pstamp, iotext, expnam
call wrufio (io, yu, jmt)
iotext = 'read (iozmbc) ((zmsmf(j,n),j=1,jmt),n=1,2)'
write (io) pstamp, iotext, expnam
call wrufio (io, zmsmf, jmt*2)
iotext = 'read (iozmbc) ((zmstf(j,n),j=1,jmt),n=1,nt)'
write (io) pstamp, iotext, expnam
call wrufio (io, zmstf, jmt*nt)
iotext = 'read (iozmbc) ((zmsm(j,n),j=1,jmt),n=1,2)'
write (io) pstamp, iotext, expnam
call wrufio (io, zmsm, jmt*2)
iotext = 'read (iozmbc) ((zmst(j,n),j=1,jmt),n=1,nt)'
write (io) pstamp, iotext, expnam
call wrufio (io, zmst, jmt*nt)
call relunit (io)
endif
endif
10200 format(//
&,26x,'Zonal Mean of Surface Boundary Condition diagnostic'
&, //,1x,' row',1x,' lat ', ' Taux ',' Tauy '
&,' Heat Flux ',' Prec-Evap ',' Surface u ',' Surface v '
&,' SST ',' SSS ',/,12x
&,' dyn/cm**2 ',' dyn/cm**2 ',' watts/m**2',' mm/day '
&,' cm/sec ',' cm/sec ',' deg C ',' ppt-35 ',/)
10201 format (1x, i4, 1x, f6.2, 8(1pe11.3))
# endif
# if defined O_time_averages
!-----------------------------------------------------------------------
! save the time mean averages on the "averaging" grid
!-----------------------------------------------------------------------
if (timavgts .and. timavgint .gt. c0) then
call avgout
endif
# endif
# if defined O_xbts
!-----------------------------------------------------------------------
! add in surface pressure gradients and missing external mode
! parts of d/dt and implicit coriolis term for all XBTs
!-----------------------------------------------------------------------
if (xbtperts) then
call uxbt3
endif
!-----------------------------------------------------------------------
! save all time averaged XBT data at end of averaging period
!-----------------------------------------------------------------------
if ((xbtts .or. eorun) .and. xbtint .gt. c0) then
call xbto
endif
# endif
# if defined O_mom_tbt
if ((tbtts .or. eorun) .and. tbtint > c0 .and. ntbtts > 0) then
!-----------------------------------------------------------------------
! construct change due to convection and filtering
!-----------------------------------------------------------------------
do j=1,jmt
do i=2,imtm1
do k=1,kmt(i,j)
do n=1,nt
tconv = tbt(i,j,k,n,10) - tbt(i,j,k,n,9)
tfilt = tbt(i,j,k,n,1) - tbt(i,j,k,n,10)
tbt(i,j,k,n,9) = tconv
tbt(i,j,k,n,10) = tfilt
enddo
enddo
enddo
enddo
!-----------------------------------------------------------------------
! average the data, write it out, then zero the accumulators
!-----------------------------------------------------------------------
rnavg = c1/ntbtts
tm(:,:,:) = 0.
errt(:) = 0.
do j=1,jmt
do i=2,imtm1
do k=1,kmt(i,j)
tm(i,j,k) = 1.
do n=1,nt
do m=1,11
tbt(i,j,k,n,m) = rnavg*tbt(i,j,k,n,m)
enddo
do m=2,10
errt(n) = errt(n) + tbt(i,j,k,n,m)
enddo
errt(n) = errt(n) - tbt(i,j,k,n,1)
enddo
enddo
do n=1,nt
tbtsf(i,j,n) = rnavg*tbtsf(i,j,n)
enddo
enddo
enddo
!-----------------------------------------------------------------------
! save all time averaged data at end of the averaging period
!-----------------------------------------------------------------------
time = year0 + accel_yr0 + (relyr - accel_yr0)*accel
if (ftbt .eq. ' ') then
call rdstmp (stamp, nyear, nmonth, nday, nhour, nmin, nsec)
nyear = time
call mkstmp (nstamp, nyear, nmonth, nday, nhour, nmin, nsec)
ftbt = file_stamp ('mom_tbt',nstamp,'.nc')
ftbt = new_file_name (ftbt)
call inqdefined (ftbt, defined)
if (defined) then
call opennext (ftbt, relyr, ntrec, iou)
else
call opennew (ftbt, iou)
endif
call mom_tbt_def (ftbt, imt, jmt, km, nt, xt, yt
&, calendar, expnam, runstamp, mapt)
endif
# if !defined O_save_time_relyear0
! make output time relative to year 1
time = time - 1.
# endif
avgper = tbtper*accel
if (avgper .le. 1e-6) avgper = 0.
# if defined O_save_time_endper
tmp = 0.
# elif defined O_save_time_startper
tmp = 1.
# else
tmp = 0.5
# endif
# if defined O_units_time_years
# if defined O_calendar_360_day
time = time - tmp*avgper/360.
# elif defined O_calendar_gregorian
time = time - tmp*avgper/365.25
# else
time = time - tmp*avgper/365.
# endif
# else
# if defined O_calendar_360_day
time = time*360. - tmp*avgper
# elif defined O_calendar_gregorian
time = time*365.25 - tmp*avgper
# else
time = time*365. - tmp*avgper
# endif
# endif
is = 1
ie = imt
js = 1
je = jmt
call mom_tbt_out (ftbt, is, ie, js, je, imt, jmt, km
&, nt, xt, yt, zt, xu, yu, zw, dxt, dyt, dzt
&, dxu, dyu, dzw, avgper, time, nstamp, mapt
&, tbt, tbtsf, tlat, tlon, ulat, ulon, kmt
&, mskhr, tm)
write (stdout,'(a,i4,a,a,a,i10,a,a)')
& '=> Ocn tracer term balances #', ntrec, ' written to '
&, trim(ftbt),' on ts = ',itt, ', ', stamp
!-----------------------------------------------------------------------
! zero the accumulators
!-----------------------------------------------------------------------
do j=1,jmt
do i=1,imt
do n=1,nt
do k=1,kmt(i,j)
do m=1,11
tbt(i,j,k,n,m) = 0.
enddo
enddo
tbtsf(i,j,n) = 0.
enddo
enddo
enddo
ntbtts = 0
endif
# endif
# if defined O_save_carbon_carbonate_chem
!-----------------------------------------------------------------------
! accumulate sea surface time averages
!-----------------------------------------------------------------------
if (timavgperts .and. eots) then
nta_sscar = nta_sscar + 1
do j=2,jmtm1
do i=2,imtm1
ta_sspH(i,j) = ta_sspH(i,j) + sspH(i,j)
ta_ssCO3(i,j) = ta_ssCO3(i,j) + ssCO3(i,j)
ta_ssOc(i,j) = ta_ssOc(i,j) + ssOc(i,j)
ta_ssOa(i,j) = ta_ssOa(i,j) + ssOa(i,j)
ta_sspCO2(i,j) = ta_sspCO2(i,j) + sspCO2(i,j)
enddo
enddo
endif
# endif
#endif
return
end
#if defined O_mom && defined O_time_step_monitor
subroutine ta_mom_tsi (m)
!=======================================================================
! ocean data time integral averaging
! input:
! m = flag (0 = zero, 1 = accumulate, 2 = write)
!=======================================================================
implicit none
integer i, j, k, ib(10), ic(10), m, n
real dens, tq, sq, drodt, drods, drhodt, drhods, ddensdtdt
real ddensdtds, ddensdsds, rntatio, rnv_otsf, rnt_slh, tarea
real cstdyt, area, s, tins, tinsit, dins, sum, d, tmp, tmp1, tmp2
include "size.h"
include "param.h"
include "pconst.h"
include "stdunits.h"
include "csbc.h"
include "coord.h"
include "cregin.h"
include "grdvar.h"
include "diag.h"
include "diaga.h"
include "emode.h"
include "levind.h"
include "mw.h"
include "iounit.h"
# if defined O_npzd
include "npzd.h"
# endif
# if defined O_sed && !defined O_sed_uncoupled
include "sed.h"
# endif
include "state.h"
include "dens.h"
# if defined O_tai_otsf
real otsf(km)
# endif
real dmsk(imt,jmt)
!-----------------------------------------------------------------------
! time averaged integrated data
!-----------------------------------------------------------------------
if (m .eq. 0.) then
! zero
ntatio = 0
# if defined O_tai_otsf_from_averages
nv_otsf = 0
v_otsf(:,:) = 0.
# endif
# if defined O_tai_slh_from_averages
nt_slh = 0
t_slh(:,:,:,:) = 0.
# endif
tai_otmax = 0.
tai_otmin = 0.
tai_slh = 0.
tai_ek = 0.
tai_t = 0.
tai_s = 0.
tai_tvar = 0.
tai_svar = 0.
tai_dt = 0.
tai_ds = 0.
tai_scan = 0.
tai_hflx = 0.
tai_sflx = 0.
tai_dic = 0.
tai_dicflx = 0.
tai_dicwflx = 0.
tai_alk = 0.
tai_o2 = 0.
tai_o2flx = 0.
tai_po4 = 0.
tai_p = 0.
tai_z = 0.
tai_d = 0.
tai_no3 = 0.
tai_di = 0.
tai_c14 = 0.
tai_dc14 = 0.
tai_c14flx = 0.
tai_sspH = 0.
tai_ssCO3 = 0.
tai_ssOc = 0.
tai_ssOa = 0.
tai_sspCO2 = 0.
tai_olivine=0.
tai_olivine_dis=0.
elseif (m .eq. 1) then
! accumulate
ntatio = ntatio + 1
# if defined O_tai_otsf && !defined O_tai_otsf_from_averages
tmp1 = 0.
tmp2 = 0.
if (nv_otsf .gt. 0) then
tmp1 = -1.e20
tmp2 = 1.e20
rnv_otsf = 1./float(nv_otsf)
do j=jsot,jeot
do k=1,km
v_otsf(j,k) = v_otsf(j,k)*rnv_otsf*dzt(k)*csu(j)
enddo
! integrate up from the bottom
otsf(km) = 0.
do k=km-1,1,-1
otsf(k) = otsf(k+1) - v_otsf(j,k)
enddo
do k=ksot,keot
if (otsf(k) .gt. tmp1) tmp1 = otsf(k)
if (otsf(k) .lt. tmp2) tmp2 = otsf(k)
enddo
enddo
endif
tai_otmax = tai_otmax + tmp1
tai_otmin = tai_otmin + tmp2
nv_otsf = 0
v_otsf(:,:) = 0.
# endif
# if defined O_tai_slh && !defined O_tai_slh_from_averages
tmp = 0.
if (nt_slh .gt. 0) then
t_slh(:,:,:,:) = t_slh(:,:,:,:)/float(nt_slh)
tarea = 0.
do j=2,jmtm1
cstdyt = cst(j)*dyt(j)
do i=2,imtm1
if (kmt(i,j) .ge. 1) then
area = dxt(i)*cstdyt
tarea = tarea + area
do k=1,kmt(i,j)
s = (t_slh(i,j,k,2)*1000.) + 35.
d = zt(k)*0.01
tins = tinsit(t_slh(i,j,k,1), s, d)
dins = dens(tins-to(k), t_slh(i,j,k,2)-so(k), k)
! volume integrate (rho - rho0)/rho0
! assume rho0 = 1g/cm3 so difference is just delta rho
! dslh = sum[((rho1-rho0)/rho0)*dvol]/area
! - sum[((rho2-rho0)/rho0)*dvol}/area
! = sum[((rho1-rho2)/rho0)*dvol]/area
tmp = tmp - (dins - d_slh(i,j,k))*area*dzt(k)
enddo
endif
enddo
enddo
if (tarea .gt. 0) tmp = tmp/tarea
endif
tai_slh = tai_slh + tmp
nt_slh = 0
t_slh(:,:,:,:) = 0.
# endif
do j=2,jmtm1
do k=km,1,-1
ektot(0,1) = ektot(0,1) + ektot(k,j)
enddo
do n=1,nt
do k=km,1,-1
tbar(0,n,1) = tbar(0,n,1) + tbar(k,n,j)
travar(0,n,1) = travar(0,n,1) + travar(k,n,j)
dtabs(0,n,1) = dtabs(0,n,1) + dtabs(k,n,j)
enddo
enddo
enddo
ektot(0,1) = ektot(0,1)/ucellv
do n=1,nt
tbar(0,n,1) = tbar(0,n,1)/tcellv
travar(0,n,1) = travar(0,n,1)/tcellv - tbar(0,n,1)**2
dtabs(0,n,1) = dtabs(0,n,1)/tcellv
enddo
tai_ek = tai_ek + ektot(0,1)
tai_t = tai_t + tbar(0,itemp,1)
tai_s = tai_s + tbar(0,isalt,1)
tai_tvar = tai_tvar + travar(0,itemp,1)
tai_svar = tai_svar + travar(0,isalt,1)
tai_dt = tai_dt + dtabs(0,itemp,1)
tai_ds = tai_ds + dtabs(0,isalt,1)
tai_scan = tai_scan + mscan
dmsk(:,:) = 1.
where (kmt(:,:) .eq. 0) dmsk(:,:) = 0.
call areaavg (sbc(1,1,ihflx), dmsk, tmp)
tai_hflx = tai_hflx + tmp
call areaavg (sbc(1,1,isflx), dmsk, tmp)
tai_sflx = tai_sflx + tmp
# if defined O_carbon
tai_dic = tai_dic + tbar(0,idic,1)
call areaavg (sbc(1,1,idicflx), dmsk, tmp)
tai_dicflx = tai_dicflx + tmp
# if defined O_sed && !defined O_sed_uncoupled
tai_dicwflx = tai_dicwflx + dicwflx
# endif
# if defined O_carbon_14
tai_c14 = tai_c14 + tbar(0,ic14,1)
tai_dc14 = tai_dc14 + dc14bar/tcellv
call areaavg (sbc(1,1,ic14flx), dmsk, tmp)
tai_c14flx = tai_c14flx + tmp
# endif
# endif
# if defined O_npzd_alk
tai_alk = tai_alk + tbar(0,ialk,1)
# endif
# if defined O_npzd_o2
tai_o2 = tai_o2 + tbar(0,io2,1)
call areaavg (sbc(1,1,io2flx), dmsk, tmp)
tai_o2flx = tai_o2flx + tmp
# endif
# if defined O_npzd
tai_po4 = tai_po4 + tbar(0,ipo4,1)
tai_p = tai_p + tbar(0,iphyt,1)
tai_z = tai_z + tbar(0,izoop,1)
tai_d = tai_d + tbar(0,idetr,1)
# if defined O_npzd_nitrogen
tai_no3 = tai_no3 + tbar(0,ino3,1)
tai_di = tai_di + tbar(0,idiaz,1)
# endif
# endif
# if defined O_cfcs_data || defined O_cfcs_data_transient
tai_cfc11 = tai_cfc11 + tbar(0,icfc11,1)
call areaavg (sbc(1,1,icfc11flx), dmsk, tmp)
tai_cfc11flx = tai_cfc11flx + tmp
tai_cfc12 = tai_cfc12 + tbar(0,icfc12,1)
call areaavg (sbc(1,1,icfc12flx), dmsk, tmp)
tai_cfc12flx = tai_cfc12flx + tmp
# endif
# if defined O_save_carbon_carbonate_chem
dmsk(:,:) = 1.
where (sspH(:,:).gt.11.5 .or. sspH(:,:).le.0.5) dmsk(:,:) = 0.
call areaavg (sspH, dmsk, tmp)
tai_sspH = tai_sspH + tmp
call areaavg (ssOc, dmsk, tmp)
tai_ssOc = tai_ssOc + tmp
call areaavg (ssOa, dmsk, tmp)
tai_ssOa = tai_ssOa + tmp
call areaavg (ssCO3, dmsk, tmp)
tai_ssCO3 = tai_ssCO3 + tmp
call areaavg (sspCO2, dmsk, tmp)
tai_sspCO2 = tai_sspCO2 + tmp
# endif
call areaavg(coa_olivine(1,1),dmsk,tmp)
tai_olivine=tai_olivine+tmp
! write(*,*) tmp
call areaavg(coa_diss(1,1),dmsk,tmp)
tai_olivine_dis=tai_olivine_dis+tmp
elseif (m .eq. 2 .and. ntatio .ne. 0) then
! average
rntatio = 1./float(ntatio)
# if defined O_tai_otsf_from_averages
if (nv_otsf .gt. 0) then
rnv_otsf = 1./float(nv_otsf)
tai_otmax = -1.e20
tai_otmin = 1.e20
do j=jsot,jeot
do k=1,km
v_otsf(j,k) = v_otsf(j,k)*rnv_otsf*dzt(k)*csu(j)
enddo
! integrate up from the bottom
otsf(km) = 0.
do k=km-1,1,-1
otsf(k) = otsf(k+1) - v_otsf(j,k)
enddo
do k=ksot,keot
if (otsf(k) .gt. tai_otmax) tai_otmax = otsf(k)
if (otsf(k) .lt. tai_otmin) tai_otmin = otsf(k)
enddo
enddo
endif
# else
tai_otmax = tai_otmax*rntatio
tai_otmin = tai_otmin*rntatio
# endif
# if defined O_tai_slh_from_averages
if (nt_slh .gt. 0) then
rnt_slh = 1./float(nt_slh)
t_slh(:,:,:,:) = t_slh(:,:,:,:)*rnt_slh
tai_slh = 0.0
tarea = 0.0
do j=2,jmtm1
cstdyt = cst(j)*dyt(j)
do i=2,imtm1
if (kmt(i,j) .ge. 1) then
area = dxt(i)*cstdyt
tarea = tarea + area
do k=1,kmt(i,j)
s = (t_slh(i,j,k,2)*1000.) + 35.
d = zt(k)*0.01
tins = tinsit(t_slh(i,j,k,1), s, d)
dins = dens(tins-to(k), t_slh(i,j,k,2)-so(k), k)
! volume integrate (rho - rho0)/rho0
! assume rho0 = 1g/cm3 so difference is just delta rho
! dslh = sum[((rho1-rho0)/rho0)*dvol]/area
! - sum[((rho2-rho0)/rho0)*dvol}/area
! = sum[((rho1-rho2)/rho0)*dvol]/area
tai_slh = tai_slh - (dins - d_slh(i,j,k))*area*dzt(k)
enddo
endif
enddo
enddo
if (tarea .gt. 0) tai_slh = tai_slh/tarea
endif
# else
tai_slh = tai_slh*rntatio
# endif
tai_ek = tai_ek*rntatio
tai_t = tai_t*rntatio
tai_s = tai_s*rntatio
tai_tvar = tai_tvar*rntatio
tai_svar = tai_svar*rntatio
tai_dt = tai_dt*rntatio
tai_ds = tai_ds*rntatio
tai_scan = tai_scan*rntatio
tai_hflx = tai_hflx*rntatio
tai_sflx = tai_sflx*rntatio
tai_dic = tai_dic*rntatio
tai_dicflx = tai_dicflx*rntatio
tai_dicwflx = tai_dicwflx*rntatio
tai_alk = tai_alk*rntatio
tai_o2 = tai_o2*rntatio
tai_o2flx = tai_o2flx*rntatio
tai_po4 = tai_po4*rntatio
tai_p = tai_p*rntatio
tai_z = tai_z*rntatio
tai_d = tai_d*rntatio
tai_no3 = tai_no3*rntatio
tai_di = tai_di*rntatio
tai_c14 = tai_c14*rntatio
tai_dc14 = tai_dc14*rntatio
tai_c14flx = tai_c14flx*rntatio
tai_cfc11 = tai_cfc11*rntatio
tai_cfc11flx = tai_cfc11flx*rntatio
tai_cfc12 = tai_cfc12*rntatio
tai_cfc12flx = tai_cfc12flx*rntatio
tai_sspH = tai_sspH*rntatio
tai_ssCO3 = tai_ssCO3*rntatio
tai_ssOc = tai_ssOc*rntatio
tai_ssOa = tai_ssOa*rntatio
tai_sspCO2 = tai_sspCO2*rntatio
tai_olivine=tai_olivine*rntatio
tai_olivine_dis=tai_olivine_dis*rntatio
endif
! write(*,*) tai_olivine
return
end
#endif
#if defined O_mom && defined O_tai_slh
real function tinsit (ti, si, di)
!=======================================================================
! function to calculate insitu temperature using Newtons method
! input:
! ti = potential temperature (C)
! si = salinity (psu)
! di = depth (m)
! output:
! tinsit = institu temperataure
!=======================================================================
implicit none
integer iter
real, intent(in) :: di, si, ti
real(kind=8) d, dfdt, error, s, t, ta, tb, tol
error = 1.
iter = 0
tol = 1.e-5
t = ti
d = di
s= si
do while (abs(error) .gt. tol .and. iter .lt. 25)
iter = iter + 1
call potem(t, s, d, ta)
call potem(t+1.e-3 ,s, d, tb)
dfdt = (tb - ta)*1.e3
error = (ta - ti)/dfdt
t = t - error
enddo
if (abs(error) .gt. tol) then
print*, "=> Error: insitu temperature did not converge"
print*, " ti = ", ti
print*, " si = ", si
print*, " di = ", di
print*, " error = ", abs(error)
print*, " tinsit = ", t
print*, " Stopped in function tinsit in diago.F"
stop
endif
tinsit = t
return
end
#endif