subroutine diago !======================================================================= ! write out diagnostic output !======================================================================= #include "param.h" #include "calendar.h" #include "ctmb.h" #include "coord.h" #include "cprnts.h" #include "cregin.h" #if defined tracer_averages # include "ctavg.h" #endif #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" #if defined stream_function && defined debug_extmode dimension ext(imt,jmt,2) #endif #if defined uvic_tbt # include "tbt.h" real errt(nt), tm(imt,jmt,km) #endif character(120) :: fname, ftbt, file_stamp save ftbt data ftbt /' '/ integer ntrec real time #if defined diagnostic_surf_height && defined stream_function # include "mw.h" 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 dimension psgrad(imt,jmt,2), divf(imt,jmt) dimension dsp(imt,jmt) save divf, 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 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 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) 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_dic13, tai_dic13flx, tai_phytc13 &, tai_zoopc13, tai_detrc13, tai_diazc13 &, 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, 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) #endif #if defined 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 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 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 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 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 show_external_mode # if defined rigid_lid_surface_pressure || defined 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 rigid_lid_surface_pressure call getunit (io, 'surf_press.dta' &, 'unformatted sequential append ieee') # endif # if defined 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 rigid_lid_surface_pressure &, ' to file surf_press.dta on ts=',itt, stamp # endif # if defined implicit_free_surface &, ' to file ifree_surf.dta on ts=',itt, stamp # endif call relunit (io) endif endif # endif # if defined 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) # if defined debug_extmode ! construct ubar and vbar for comparison with ! rigid_lid_surface_pressure scl=1.0 do jrow=1,jmt-1 j = jrow do i=2,imt-1 diag1 = psi(i+1,jrow+1,1) - psi(i ,jrow,1) diag0 = psi(i ,jrow+1,1) - psi(i+1,jrow,1) ext(i,j,1) = -(diag1+diag0)*dyu2r(jrow)*hr(i,jrow) ext(i,j,2) = (diag1-diag0)*dxu2r(i)*hr(i,jrow)*csur(jrow) enddo enddo call border (ext(1,1,1), 'u even') call border (ext(1,1,2), 'u odd') write (stdout,8000) ' ubar ' &, itt, xt(is), xt(ie), yt(js), yt(je), scl call matrix (ext(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 (ext(1,1,2), imt, is, ie, -js, -je, scl) # endif 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 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 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 time_averages !----------------------------------------------------------------------- ! save the time mean averages on the "averaging" grid !----------------------------------------------------------------------- if (timavgts .and. timavgint .gt. c0) then call avgout endif #endif #if defined 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 uvic_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 print*, ' tbt error: ', errt !----------------------------------------------------------------------- ! save all time averaged data at end of the averaging period !----------------------------------------------------------------------- # if defined uvic_multiple_files ntrec = 0 ftbt = file_stamp ('tbt_mom',stamp,'.nc') # else ntrec = 1 if (ftbt .eq. ' ') then ntrec = 0 ftbt = file_stamp ('tbt_mom',stamp,'.nc') endif # endif avgper = dtts*ntbtts/daylen is = 1 ie = imt js = 1 je = jmt time = relyr + year0 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, ntrec, timunit, expnam &, avgper, time, stamp, 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 return end #if defined time_step_monitor 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" # if defined uvic_npzd # include "npzd.h" # endif 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) # if defined uvic_tai_otsf real otsf(km) # endif # if defined uvic_tai_slh 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 # endif !----------------------------------------------------------------------- ! time averaged integrated data !----------------------------------------------------------------------- if (m .eq. 0.) then ! zero ntatio = 0 # if defined uvic_tai_slh nt_slh = 0 t_slh(:,:,:,:) = 0. # endif # if defined uvic_tai_otsf nv_otsf = 0 v_otsf(:,:) = 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_dic13 = 0. tai_dic13flx = 0. tai_phytc13 = 0. tai_zoopc13 = 0. tai_detrc13 = 0. tai_diazc13 = 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. 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 # if defined uvic_carbon tai_dic = tai_dic + tbar(0,idic,1) call areaavg (sbc(1,1,idicflx), dmsk, tmp) tai_dicflx = tai_dicflx + tmp # if defined uvic_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 # if defined osu_c13 tai_dic13 = tai_dic13 + tbar(0,idic13,1) tai_phytc13 = tai_phytc13 + tbar(0,iphytc13,1) tai_zoopc13 = tai_zoopc13 + tbar(0,izoopc13,1) tai_detrc13 = tai_detrc13 + tbar(0,idetrc13,1) tai_diazc13 = tai_diazc13 + tbar(0,idiazc13,1) c tai_dc13 = tai_dc13 + dc13bar/tcellv call areaavg (sbc(1,1,idic13flx), dmsk, tmp) tai_dic13flx = tai_dic13flx + tmp # endif # endif # if defined uvic_alk tai_alk = tai_alk + tbar(0,ialk,1) # endif # if defined uvic_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 uvic_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 uvic_nitrogen tai_no3 = tai_no3 + tbar(0,ino3,1) tai_di = tai_di + tbar(0,idiaz,1) # endif # endif # if defined uvic_cfc11 tai_cfc11 = tai_cfc11 + tbar(0,icfc11,1) call areaavg (sbc(1,1,icfc11flx), dmsk, tmp) tai_cfc11flx = tai_cfc11flx + tmp # endif # if defined uvic_cfc12 tai_cfc12 = tai_cfc12 + tbar(0,icfc12,1) call areaavg (sbc(1,1,icfc12flx), dmsk, tmp) tai_cfc12flx = tai_cfc12flx + tmp # endif elseif (m .eq. 2 .and. ntatio .ne. 0) then ! average rntatio = 1./float(ntatio) # if defined uvic_tai_otsf 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 # endif # if defined uvic_tai_slh 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 # 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_dic13 = tai_dic13*rntatio tai_dic13flx = tai_dic13flx*rntatio tai_phytc13 = tai_phytc13*rntatio tai_zoopc13 = tai_zoopc13*rntatio tai_detrc13 = tai_detrc13*rntatio tai_diazc13 = tai_diazc13*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 endif return end #endif #if defined uvic_tai_slh 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 #endif