! source file: /Users/oschlies/UVIC/master/testcase/updates/diago.F subroutine diago !======================================================================= ! write out diagnostic output !======================================================================= include "param.h" include "calendar.h" include "ctmb.h" include "coord.h" include "cprnts.h" include "cregin.h" include "ctavg.h" include "diag.h" include "docnam.h" include "emode.h" include "grdvar.h" include "iounit.h" include "levind.h" include "scalar.h" include "stab.h" include "state.h" include "switch.h" include "task_on.h" include "tmngr.h" include "vmixc.h" character(120) :: fname, ftbt, file_stamp save ftbt data ftbt /' '/ integer ntrec real time !----------------------------------------------------------------------- ! 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) then call ta_mom_tsi (2) if (iotsi .eq. stdout .or. iotsi .lt. 0) then write (stdout,9601) itt, stamp, tai_ek, tai_dt, tai_ds, tai_t &, tai_s, tai_tvar, tai_svar, tai_scan endif if (iotsi .ne. stdout .or. iotsi .lt. 0) then if (iotsi .ne. stdout) then avgper = dtts*ntatio/daylen time = relyr + year0 call def_tsi call def_tsi_mom (fname) call mom_tsi_out (fname, avgper, time, stamp, 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_n2o, tai_n2oflx &, tai_abioto2, tai_abioto2flx &, tai_prefo2 &, tai_prefpo4 &, tai_prefno3 &, 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, ntrec) endif endif 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 1 !----------------------------------------------------------------------- ! 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) then call avgout endif return end subroutine ta_mom_tsi (m) !======================================================================= ! ocean data time integral averaging ! input: ! m = flag (0 = zero, 1 = accumulate, 2 = write) ! based on code by: M. Eby !======================================================================= ! implicit none include "param.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" integer i, j, jrow, k, kmax, ib(10), ic(10) include "state.h" include "dens.h" real d, dens, dins, rntatio, s, sum, tmp, tins, tinsit real C2K, p001, p035, dmsk(imt,jmt) real otsf(km) character(120) :: fname, new_file_name logical exists real, allocatable :: d_slh(:,:,:) save d_slh C2K = 273.15 p001 = 0.001 p035 = 0.035 if (.not. allocated (d_slh)) then fname = new_file_name ("slh_ref.nc") inquire (file=trim(fname), exist=exists) if (exists) then allocate ( d_slh(2:imtm1,2:jmtm1,km) ) d_slh(:,:,:) = 0. ib(:) = 1 ic(:) = imtm2 ic(2) = jmtm2 ic(3) = km ic(4) = 1 call openfile (fname, iou) call getvara ('temp', iou, imtm2*jmtm2*km, ib, ic &, t_slh(2,2,1,1), c1, -C2K) call getvara ('salinity', iou, imtm2*jmtm2*km, ib, ic &, t_slh(2,2,1,2), p001, -p035) call closefile (iou) do j=2,jmtm1 do i=2,imtm1 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) d_slh(i,j,k) = dens(tins-to(k), t_slh(i,j,k,2)-so(k), k) enddo enddo enddo endif endif !----------------------------------------------------------------------- ! time averaged integrated data !----------------------------------------------------------------------- if (m .eq. 0.) then ! zero ntatio = 0 nt_slh = 0 t_slh(:,:,:,:) = 0. nv_otsf = 0 v_otsf(:,:) = 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_alk = 0. tai_o2 = 0. tai_o2flx = 0. tai_n2o = 0. tai_n2oflx = 0. tai_abioto2 = 0. tai_abioto2flx = 0. tai_prefo2 = 0. tai_prefpo4 = 0. tai_prefno3 = 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. elseif (m .eq. 1) then ! accumulate ntatio = ntatio + 1 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_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_n2o = tai_n2o + tbar(0,in2o,1) call areaavg (sbc(1,1,in2oflx), dmsk, tmp) tai_n2oflx = tai_n2oflx + tmp tai_abioto2 = tai_abioto2 + tbar(0,iabioto2,1) call areaavg (sbc(1,1,iabioto2flx), dmsk, tmp) tai_abioto2flx = tai_abioto2flx + tmp tai_prefo2 = tai_prefo2 + tbar(0,iprefo2,1) tai_prefpo4 = tai_prefpo4 + tbar(0,iprefpo4,1) tai_prefno3 = tai_prefno3 + tbar(0,iprefno3,1) 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) elseif (m .eq. 2 .and. ntatio .ne. 0) then ! average rntatio = 1./float(ntatio) 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 if (nt_slh .gt. 0) then rnt_slh = 1./float(nt_slh) t_slh(:,:,:,:) = t_slh(:,:,:,:)*rnt_slh if (.not. allocated (d_slh)) then ! if not allocated get reference from the first density allocate ( d_slh(imt,jmt,km) ) do j=2,jmtm1 do i=2,imtm1 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) d_slh(i,j,k) = dens(tins-to(k),t_slh(i,j,k,2)-so(k),k) enddo enddo enddo else ! if reference is allocated calculate relative sea level height tai_slh = 0.0 tarea = 0.0 do j=2,jmtm1 cosdyt = cst(j)*dyt(j) do i=2,imtm1 if (kmt(i,j) .ge. 1) then area = dxt(i)*cosdyt 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) 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 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_alk = tai_alk*rntatio tai_o2 = tai_o2*rntatio tai_o2flx = tai_o2flx*rntatio tai_n2o = tai_n2o*rntatio tai_n2oflx = tai_n2oflx*rntatio tai_abioto2 = tai_abioto2*rntatio tai_abioto2flx = tai_abioto2flx*rntatio tai_prefo2 = tai_prefo2*rntatio tai_prefpo4 = tai_prefpo4*rntatio tai_prefno3 = tai_prefno3*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 endif return end real function tinsit (ti, si, di) !======================================================================= ! function to calculate insitu temperature using Newtons method ! based on code by e.c.wiebe ! input: ! ti = potential temperature (C) ! si = salinity (psu) ! di = depth (m) ! output: ! tinsit = institu temperataure ! based on code by: M. Eby !======================================================================= 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