! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mom/diago.F subroutine diago !======================================================================= ! 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, 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 "csbc.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" include "cembm.h" include "ctavg.h" include "timeavgs.h" include "diaga.h" include "npzd.h" !----------------------------------------------------------------------- ! 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),/) !----------------------------------------------------------------------- ! print integrals for monitoring the time step !----------------------------------------------------------------------- if (tsiperts .and. eots) call ta_mom_tsi (1) if (tsits .and. ntatio .gt. 0 .and. eoseg) 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) avgtimtsi = time avgpertsi = dtts*ntatio*secday if (avgpertsi .le. 1e-6) avgpertsi = 0. tmp = 0.5 avgtimtsi = avgtimtsi - tmp*avgpertsi/365. ! 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 ! subtract out any land to ocean flux tai_cfa2o = tai_cfa2o - tai_dicwflx*yrlen*daylen*12.e-21 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) call mom_tsi_out (fname, avgpertsi, avgtimtsi, 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 &, 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) !----------------------------------------------------------------------- ! 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) !----------------------------------------------------------------------- ! 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,')') !----------------------------------------------------------------------- ! 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'/) !----------------------------------------------------------------------- ! 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)',/) !----------------------------------------------------------------------- ! 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)') !----------------------------------------------------------------------- ! save the time mean averages on the "averaging" grid !----------------------------------------------------------------------- if (timavgts .and. timavgint .gt. c0 .and. eoseg) then call avgout endif !----------------------------------------------------------------------- ! accumulate sea surface time averages !----------------------------------------------------------------------- if (timavgperts .and. eots) then nta_sscar = nta_sscar + 1 do j=2,jmtm1 do i=2,imtm1 ! sspH is really H+ until the final averaging 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 return end 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" include "npzd.h" include "sed.h" include "state.h" include "dens.h" real otsf(km) real dmsk(imt,jmt) !----------------------------------------------------------------------- ! time averaged integrated data !----------------------------------------------------------------------- if (m .eq. 0.) then ! zero ntatio = 0 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. elseif (m .eq. 1) then ! accumulate ntatio = ntatio + 1 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. 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. 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 tai_dic = tai_dic + tbar(0,idic,1) call areaavg (sbc(1,1,idicflx), dmsk, tmp) tai_dicflx = tai_dicflx + tmp tai_dicwflx = tai_dicwflx + dicwflx 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 tai_alk = tai_alk + tbar(0,ialk,1) tai_o2 = tai_o2 + tbar(0,io2,1) call areaavg (sbc(1,1,io2flx), dmsk, tmp) tai_o2flx = tai_o2flx + tmp 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) tai_no3 = tai_no3 + tbar(0,ino3,1) tai_di = tai_di + tbar(0,idiaz,1) 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 dmsk(:,:) = 1. where (kmt(:,:) .eq. 0) dmsk(:,:) = 0. ! sspH is really H+ until the final averaging 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 elseif (m .eq. 2 .and. ntatio .ne. 0) then ! average rntatio = 1./float(ntatio) tai_otmax = tai_otmax*rntatio tai_otmin = tai_otmin*rntatio tai_slh = tai_slh*rntatio 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 ! sspH is really H+ until the final averaging if (tai_sspH .gt. 0.) & tai_sspH = -alog10(tai_sspH*rntatio) tai_ssCO3 = tai_ssCO3*rntatio tai_ssOc = tai_ssOc*rntatio tai_ssOa = tai_ssOa*rntatio tai_sspCO2 = tai_sspCO2*rntatio endif return end 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