! source file: /Users/oschlies/UVIC/master/testcase/updates/tracer.F subroutine tracer (joff, js, je, is, ie) !======================================================================= ! compute tracers at "tau+1" for rows js through je in the MW. ! input: ! joff = offset relating "j" in the MW to latitude "jrow" ! js = starting row in the MW ! je = ending row in the MW ! is = starting longitude index in the MW ! ie = ending longitude index in the MW ! based on code by: R. C. Pacanowski !======================================================================= include "param.h" parameter (istrt=2, iend=imt-1) include "accel.h" include "coord.h" include "cregin.h" include "csbc.h" include "emode.h" include "grdvar.h" include "hmixc.h" include "isopyc.h" include "levind.h" include "mw.h" include "scalar.h" include "switch.h" include "timeavgs.h" include "tmngr.h" include "vmixc.h" dimension twodt(km) include "fdift.h" include "diaga.h" include "ice.h" include "atm.h" include "npzd.h" integer index real snpzd(ntnpzd), tnpzd(ntnpzd) real rctheta, declin, gl, impo, expo, npp real remi, excr, graz, morp, morpt, morz, temp, swr, dayfrac real phin, dz, prca, dprca, nud, bct, tap, fo2, so2 real funo2, sn2o real npp_D, graz_D, morp_D, no3flag, deni, nfix real so4flag real fy, fyz real src(imt,km,jsmw:jmw,nsrc) src(:,:,:,:) = 0. !----------------------------------------------------------------------- ! bail out if starting row exceeds ending row !----------------------------------------------------------------------- if (js .gt. je) return !----------------------------------------------------------------------- ! limit the longitude indices based on those from the argument list ! Note: this is currently bypassed. istrt and iend are set as ! parameters to optimize performance !----------------------------------------------------------------------- ! istrt = max(2,is) ! iend = min(imt-1,ie) !----------------------------------------------------------------------- ! build coefficients to minimize advection and diffusion computation !----------------------------------------------------------------------- limit = min(je+1+joff,jmt) - joff do j=js,limit jrow = j + joff do i=istrt-1,iend cstdxtr(i,j) = cstr(jrow)*dxtr(i) cstdxt2r(i,j) = cstr(jrow)*dxtr(i)*p5 cstdxur(i,j) = cstr(jrow)*dxur(i) ah_cstdxur(i,j) = diff_cet*cstr(jrow)*dxur(i) enddo enddo !----------------------------------------------------------------------- ! calculation of biological interactions !----------------------------------------------------------------------- declin = sin((mod(relyr,1.) - 0.22)*2.*pi)*0.4 ! declination do k=1,km twodt(k) = c2dtts*dtxcel(k) nbio(k) = twodt(k)/dtnpzd dtbio(k) = twodt(k)/nbio(k) rdtts(k) = 1./twodt(k) rnbio(k) = 1./nbio(k) enddo tap = 2.*alpha*par do j=js,je jrow = j + joff do i=is,ie if (kmt(i,jrow) .gt. 0) then ! calculate day fraction and incoming solar ! angle of incidence = lat - declin, refraction index = 1.33 rctheta = max(-1.5, min(1.5, tlat(i,jrow)/radian - declin)) rctheta = kw/sqrt(1. - (1. - cos(rctheta)**2.)/1.33**2.) dayfrac = min( 1., -tan(tlat(i,jrow)/radian)*tan(declin)) dayfrac = max(1e-12, acos(max(-1., dayfrac))/pi) swr = dnswr(i,jrow)*1e-3 ! convert to W/m^2 & *(1.0 + aice(i,jrow,2) ! attenuation by sea ice & snow & *(exp(-ki*(hice(i,jrow,2) + hsno(i,jrow,2))) -1.)) expo = 0.0 impo = 0.0 phin = 0.0 ! integrated phytoplankton prca = 0.0 ! integrated production of calcite kmax = min(kmt(i,jrow), kpzd) do k=1,kmax !----------------------------------------------------------------------- ! limit tracers to positive values !----------------------------------------------------------------------- tnpzd(1) = max(t(i,k,j,ipo4,taum1), trcmin) tnpzd(2) = max(t(i,k,j,iphyt,taum1), trcmin) tnpzd(3) = max(t(i,k,j,izoop,taum1), trcmin) tnpzd(4) = max(t(i,k,j,idetr,taum1), trcmin) tnpzd(5) = max(t(i,k,j,ino3,taum1), trcmin) tnpzd(6) = max(t(i,k,j,idiaz,taum1), trcmin) swr = swr*exp(-kc*phin) phin = phin + tnpzd(2)*dzt(k) gl = tap*swr*exp(ztt(k)*rctheta) impo = expo*dztr(k) bct = bbio**(cbio*t(i,k,j,itemp,taum1)) ! decrease remineralisation rate in oxygen minimum zone nud = nud0*(0.65+0.35*tanh(t(i,k,j,io2,taum1)*1000.-6.)) cao 11.3.09 old decrease from 1.0 to 0.3 (factor of 3.3) cao new decrease from 1.0 to 0.2 (factor of 5 as said in Schmittner et al.2008) c nud = nud0*(0.60+0.40*tanh(t(i,k,j,io2,taum1)*1000.-6.)) c nud = nud0*(0.75+0.25*tanh(t(i,k,j,io2,taum1)*1000.-6.)) c nud = nud0 !----------------------------------------------------------------------- ! call the npzd model !----------------------------------------------------------------------- fo2 = 0.5*tanh(t(i,k,j,io2,taum1)*1000. - 5.) cao so4flag = 0.5+sign(0.5,t(i,k,j,ino3,taum1)-trcmin) so4flag = 0.5+sign(0.5,t(i,k,j,ino3,taum1)-.5) if (fo2.gt.0.) so4flag = 1. c print*,i,j,k,so4flag,fo2 call npzd_src (tnpzd, nbio(k), dtbio(k), gl, bct, impo &, dzt(k), dayfrac, wd(k), rkwz(k), nud &, snpzd, expo, graz, morp, morz &, npp, morpt, remi, excr &, npp_D, graz_D, morp_D, nfix &, so4flag & ) snpzd(1:4) = snpzd(1:4)*rdtts(k) snpzd(5:6) = snpzd(5:6)*rdtts(k) expo = expo*rnbio(k) rexpo(i,k,j) = expo rgraz(i,k,j) = graz*rnbio(k) rmorp(i,k,j) = morp*rnbio(k) rmorz(i,k,j) = morz*rnbio(k) rnpp(i,k,j) = npp*rnbio(k) rmorpt(i,k,j) = morpt*rnbio(k) rremi(i,k,j) = remi*rnbio(k) rexcr(i,k,j) = excr*rnbio(k) rnpp_D(i,k,j) = npp_D*rnbio(k) rgraz_D(i,k,j) = graz_D*rnbio(k) rmorp_D(i,k,j) = morp_D*rnbio(k) rnfix(i,k,j) = nfix*rnbio(k) !----------------------------------------------------------------------- ! calculate detritus at the bottom !----------------------------------------------------------------------- if (k .eq. kmt(i,jrow)) then snpzd(4) = snpzd(4) + expo expo = 0. rexpo(i,k,j) = 0. elseif (k .eq. kpzd .and. expo .ne. 0.) then expo = expo*dzt(kpzd) do m=kpzd+1,kmt(i,jrow) cao in the standard configuration, kpzd=km, so this loop is never executed. cao a.o., 4.2.10 expo = expo*dztr(m) ! if below the npzd model move all biology to export src(i,m,j,isphyt) = -rdtts(m)*t(i,m,j,iphyt,taum1) t(i,m,j,iphyt,taup1) = 0. src(i,m,j,iszoop) = -rdtts(m)*t(i,m,j,izoop,taum1) t(i,m,j,izoop,taup1) = 0. src(i,m,j,isdetr) = -rdtts(m)*t(i,m,j,idetr,taum1) t(i,m,j,idetr,taup1) = 0. src(i,m,j,isdiaz) = -rdtts(m)*t(i,m,j,idiaz,taum1) t(i,m,j,idiaz,taup1) = 0. expo = expo - src(i,m,j,isphyt) - src(i,m,j,iszoop) & - src(i,m,j,isdetr) - src(i,m,j,isdiaz) if (m .eq. kmt(i,jrow)) then ! remineralize all detritus if at the bottom src(i,m,j,ispo4) = redptn*expo rremi(i,m,j) = expo rexpo(i,m,j) = 0. else ! calculate the proportion of detritus remineralized nud = nud0*bbio**(cbio*t(i,m,j,itemp,taum1)) ! decrease remineralisation in oxygen minimum zone & *(0.65+0.35*tanh(t(i,m,j,io2,taum1)*1000.-6.)) c & *(0.60+0.40*tanh(t(i,m,j,io2,taum1)*1000.-6.)) c & *(0.75+0.25*tanh(t(i,m,j,io2,taum1)*1000.-6.)) remi = expo*nud/(nud + wd(m)) src(i,m,j,ispo4) = redptn*remi rremi(i,m,j) = remi rexpo(i,m,j) = expo - remi ! calculate the detritus exported to the next level expo = (expo - remi)*dzt(m) endif src(i,m,j,isno3) = src(i,m,j,ispo4)*redntp src(i,m,j,isdic) = src(i,m,j,ispo4)*redctp src(i,m,j,isalk) = -src(i,m,j,ispo4)*redntp*1.e-3 enddo endif !----------------------------------------------------------------------- ! set source terms !----------------------------------------------------------------------- src(i,k,j,ispo4) = snpzd(1) src(i,k,j,isphyt) = snpzd(2) src(i,k,j,iszoop) = snpzd(3) src(i,k,j,isdetr) = snpzd(4) src(i,k,j,isno3) = snpzd(5) src(i,k,j,isdiaz) = snpzd(6) sst = t(i,1,j,1,tau) sss = 1000.*t(i,1,j,2,tau) + 35. dic_in = t(i,1,j,idic,tau) ta_in = t(i,1,j,ialk,tau) pt_in = 0.5125e-3 !mol/m^3, from gasbc.F, instead of t(i,1,j,ipo4,tau) sit_in = 7.6875e-03 !mol/m^3 phlo = 6. phhi = 10. atmpres = 1.0 co2_in = co2ccn call co2calc_SWS (sst, sss, dic_in, ta_in, pt_in, sit_in &, phlo, phhi, ph, co2_in, atmpres, co2star &, dco2star, pCO2surf, dpco2) pco2z0(i,j) = pCO2surf ! production of calcite dprca = (morp+morz+graz*(1.-gamma1))*capr*redctn*rnbio(k) prca = prca + dprca*dzt(k) src(i,k,j,isdic) = (snpzd(1)*redctp - dprca) src(i,k,j,isalk) = (-snpzd(1)*redntp*1.e-3 - 2.*dprca) !----------------------------------------------------------------------- ! accumulate time averages !----------------------------------------------------------------------- if (timavgperts .and. .not. euler2) then ta_rnpp(i,k,jrow) = ta_rnpp(i,k,jrow) + rnpp(i,k,j) ta_rgraz(i,k,jrow) = ta_rgraz(i,k,jrow) + rgraz(i,k,j) ta_rmorp(i,k,jrow) = ta_rmorp(i,k,jrow) + rmorp(i,k,j) ta_rmorpt(i,k,jrow)= ta_rmorpt(i,k,jrow) + rmorpt(i,k,j) ta_rmorz(i,k,jrow) = ta_rmorz(i,k,jrow) + rmorz(i,k,j) ta_rexcr(i,k,jrow) = ta_rexcr(i,k,jrow) + rexcr(i,k,j) ta_rnpp_D(i,k,jrow) = ta_rnpp_D(i,k,jrow) & + rnpp_D(i,k,j) ta_rgraz_D(i,k,jrow) = ta_rgraz_D(i,k,jrow) & + rgraz_D(i,k,j) ta_rmorp_D(i,k,jrow) = ta_rmorp_D(i,k,jrow) & + rmorp_D(i,k,j) ta_rnfix(i,k,jrow) = ta_rnfix(i,k,jrow) + rnfix(i,k,j) endif ! calculate total export to get total import for next layer expo = expo*dzt(k) enddo kmax = kmt(i,jrow) do k=1,kmax ! limit oxygen consumption below concentrations of ! 5umol/kg as recommended in OCMIP fo2 = 0.5*tanh(t(i,k,j,io2,taum1)*1000. - 5.) ! sink of oxygen !cao 22.5.12 O2 is needed to generate the equivalent of NO3 from N2 during N2 fixation ! 0.5 H2O + 0.5 N2+1.25O2 -> HNO3 ! note that so2 is -dO2/dt so2 = src(i,k,j,ispo4)*redotp & + rnfix(i,k,j)*1.25e-3 src(i,k,j,iso2) = -so2*(0.5 + fo2) ! add denitrification as source term for NO3 no3flag = 0.5+sign(0.5,t(i,k,j,ino3,taum1)-trcmin) ! 800 = 0.8*1000 = (elec/mol O2)/(elec/mol NO3)*(mmol/mol) !cao 26.8.07 leave old factor of 0.8 ... # if defined AO_run ! 667 = 0.667*1000 ! deni = 667.*no3flag*so2*(0.5 - fo2) ! ! ao, 30.5.12 The factor 4/5=0.8 is correct: deni describes only the left-hand side of the equation, ! i.e. the equivalent O2 demand, which is 4/5 the O2. ! In the set-up of the UVic model, HNO3 is added to the right hand side of the denitrification ! equation (i.e. src(i,j,k,ino3). See eq. 28 of Paulmier et al. (2009) Biogeosciences deni = 800.*no3flag*so2*(0.5 - fo2) src(i,k,j,isno3) = src(i,k,j,isno3) - deni cao all new parameters formulated in terms of Aurelien Paulmier et al.'s 2009 BG paper and cao model stoichimetric variables. cao a.o. 22.6.12 not yet tested (only version with hardwired cao src(i,k,j,isalk) = src(i,k,j,isalk) cao & + no3flag*src(i,k,j,ispo4)*(0.5-fo2)*(119.68+redntp)*1.e-3 cao had been tested Paulmier_a = 1.e3*redctn Paulmier_z = 4.*1.e3*redotp - 4.*Paulmier_a - 8.*redntp Paulmier_R0 = Paulmier_a + 0.25*Paulmier_z src(i,k,j,isalk) = src(i,k,j,isalk) & + no3flag*src(i,k,j,ispo4)*(0.5-fo2) & * (4./5.*Paulmier_R0 + (3./5.+1.)*redntp) * 1.e-3 c & * (119.68+redntp)*1.e-3 cao using Aurelien's BG paper with (4/5a+1/5b-2/5c) = (4/5R0+3/5d-1), with R0=a+1/4z cao a=C:P = 7*16 = 112 cao z= 4*-O2:P(of remin.) - 4*a - 8*d = 4*10.6*16 - 4*112 - 8*16 = 102.4 cao d=N:P = 16 cao ==> R0 = 112+25.6 = 137.6 cao ==> (4/5R0+3/5d-1) = 110.08 + 9.6 -1 = 118.68 cao STILL: use 119.68 instead of 118.68 because Alkalinity change by PO4 change is never cao accounted for in this model. cao and add redntp of aerobic remineralisation wrongly aplied above. cao cao also add N2 fixation effect! (alkalinity production is tied to po4 uptake above, cao in the case of nitrogen fixation, this is not correct and has been substracted here). src(i,k,j,isalk) = src(i,k,j,isalk) & - rnfix(i,k,j)*1.e-3 rdeni(i,k,jrow) = deni ! use Suntharalingam etal. (2000) approach OX.5 ! the value for beta is derived from ASnn_co2emit_n AD 2010 oxygen field ! note: t(i,k,j,io2,taum1) is in mol/m3 funo2 = exp((-.1)*(t(i,k,j,io2,taum1)*1000. - 1.)/1.) ! use remineralisation rate (so2) rather than oxygen consumption rate ! (this matters only in the range 1-5mmolO2/m3, though...) sn2o = (0.5e-4 + 0.0065*funo2)*so2*(0.5 + fo2) src(i,k,j,isn2o) = sn2o ! add in N2O consumption in anoxic areas: oxygen is taken from N2O ! (factor 2 because per mol there is twice as many N2O moles as O2 moles ! are needed) ! use 4mmol O2/m3 as threshold (3.6 is the minimum reached by the model) ! the threshold 4mmol O2/m3 is also suggested by Nevison et al., GBC 2003 if (t(i,k,j,io2,taum1).lt.0.004) & src(i,k,j,isn2o) = (-2.)*so2 ! skipping equation below yields LOWER N2O inventories, because so2 can ! be negative in surface layers with positive NCP (i.e., oxygen production) if (k.le.2) src(i,k,j,isn2o) = c0 src(i,k,j,isabioto2) = c0 enddo !----------------------------------------------------------------------- ! accumulate time averages for full depth variables !----------------------------------------------------------------------- kmax = kmt(i,jrow) do k=1,kmax if (timavgperts .and. .not. euler2) then ta_rremi(i,k,jrow) = ta_rremi(i,k,jrow) + rremi(i,k,j) ta_rexpo(i,k,jrow) = ta_rexpo(i,k,jrow) + rexpo(i,k,j) ta_rdeni(i,k,jrow) = ta_rdeni(i,k,jrow) + rdeni(i,k,j) endif enddo !----------------------------------------------------------------------- ! remineralize calcite !----------------------------------------------------------------------- kmax = kmt(i,jrow) do k=1,kmax-1 src(i,k,j,isdic) = src(i,k,j,isdic) + prca*rcak(k) src(i,k,j,isalk) = src(i,k,j,isalk) + 2.*prca*rcak(k) enddo src(i,kmax,j,isdic) = src(i,kmax,j,isdic) + prca*rcab(kmax) src(i,kmax,j,isalk) = src(i,kmax,j,isalk)+2.*prca*rcab(kmax) endif enddo enddo !----------------------------------------------------------------------- ! set source for c14 !----------------------------------------------------------------------- do j=js,je jrow = j + joff do i=is,ie if (kmt(i,jrow) .gt. 0) then do k=1,kmt(i,jrow) src(i,k,j,isc14) = src(i,k,j,isdic)*rstd & - 3.836e-12*t(i,k,j,ic14,taum1) enddo endif enddo enddo !----------------------------------------------------------------------- ! solve for one tracer at a time ! n = 1 => temperature ! n = 2 => salinity ! n > 2 => other tracers (if applicable) !----------------------------------------------------------------------- do n=1,nt !----------------------------------------------------------------------- ! calculate 2* FCT tracer flux !----------------------------------------------------------------------- call adv_flux_fct(joff, js, je, is, ie, n) !----------------------------------------------------------------------- ! calculate diffusive flux across eastern and northern faces ! of "T" cells due to various parameterizations for diffusion. !----------------------------------------------------------------------- ! diffusive flux on eastern face of "T" cells do j=js,je do k=1,km do i=istrt-1,iend diff_fe(i,k,j) = & ah_cstdxur(i,j)* & (t(i+1,k,j,n,taum1) - t(i,k,j,n,taum1)) enddo enddo enddo ! diffusive flux on northern face of "T" cells ! (background for isopycnal mixing) do j=js-1,je jrow = j + joff do k=1,km do i=istrt,iend diff_fn(i,k,j) = & diff_cnt* & csu_dyur(jrow)*(t(i,k,j+1,n,taum1) - t(i,k,j,n,taum1)) enddo enddo enddo !----------------------------------------------------------------------- ! calculate diffusive flux across bottom face of "T" cells !----------------------------------------------------------------------- do j=js,je do k=1,km-1 do i=istrt,iend diff_fb(i,k,j) = diff_cbt(i,k,j)*dzwr(k)* & (t(i,k,j,n,taum1) - t(i,k+1,j,n,taum1)) enddo enddo enddo !----------------------------------------------------------------------- ! compute isopycnal diffusive flux through east, north, ! and bottom faces of T cells. !----------------------------------------------------------------------- call isoflux (joff, js, je, is, ie, n) !----------------------------------------------------------------------- ! set surface and bottom vert b.c. on "T" cells for diffusion ! and advection. for isopycnal diffusion, set adiabatic boundary ! conditions. ! note: the b.c. at adv_fb(i,k=bottom,j) is set by the above code. ! However, it is not set when k=km so it is set below. ! adv_fb(i,km,j) is always zero (to within roundoff). !----------------------------------------------------------------------- do j=js,je jrow = j + joff do i=istrt,iend kb = kmt(i,jrow) diff_fb(i,0,j) = stf(i,j,n) diff_fb(i,kb,j) = btf(i,j,n) adv_fb(i,0,j) = adv_vbt(i,0,j)*(t(i,1,j,n,tau) + & t(i,1,j,n,tau)) adv_fb(i,km,j) = adv_vbt(i,km,j)*t(i,km,j,n,tau) enddo enddo !----------------------------------------------------------------------- ! set source term for "T" cells !----------------------------------------------------------------------- source(:,:,:) = c0 if (itrc(n) .ne. 0) then do j=js,je do k=1,km do i=istrt,iend source(i,k,j) = src(i,k,j,itrc(n)) enddo enddo enddo endif !----------------------------------------------------------------------- ! solve for "tau+1" tracer using statement functions to represent ! each component of the calculation !----------------------------------------------------------------------- ! 1st: solve using all components which are treated explicitly do j=js,je jrow = j + joff do k=1,km twodt(k) = c2dtts*dtxcel(k) do i=istrt,iend t(i,k,j,n,taup1) = t(i,k,j,n,taum1) + twodt(k)*( & DIFF_Tx(i,k,j) + DIFF_Ty(i,k,j,jrow,n) + DIFF_Tz(i,k,j) & - ADV_Tx(i,k,j) - ADV_Ty(i,k,j,jrow,n) - ADV_Tz(i,k,j) & + source(i,k,j) & )*tmask(i,k,j) if (n.eq.in2o) & t(i,k,j,n,taup1) = max(1.e-12,t(i,k,j,n,taup1)) enddo enddo enddo ! 2nd: add in portion of vertical diffusion handled implicitly call ivdift (joff, js, je, istrt, iend, n, twodt) if (itt.eq.3533201.and.n.eq.iabioto2 & .or.itt.eq.3533202.and.n.eq.iabioto2) then print*,itt,iabioto2,'initialize abioto2' do j=js,je do k=1,km do i=istrt,iend sst=t(i,k,j,1,taup1) ! degree Celsius sss=1000.*t(i,k,j,2,taup1) + 35. ! oxygen saturation concentration [mol/m^3] f1 = alog((298.15 - sst)/(273.15 + sst)) f2 = f1*f1 f3 = f2*f1 f4 = f3*f1 f5 = f4*f1 o2sat = exp (2.00907 + 3.22014*f1 + 4.05010*f2 & + 4.94457*f3 - 2.56847E-1*f4 + 3.88767*f5 & + sss*(-6.24523e-3 - 7.37614e-3*f1 - 1.03410e-2*f2 & - 8.17083E-3*f3) - 4.88682E-7*sss*sss) ! Convert from ml/l to mol/m^3 o2sat = o2sat/22391.6*1000.0 t(i,k,j,n,taup1) = o2sat enddo enddo print*,j,t(100,1,j,iabioto2,taup1),t(100,1,j,iabioto2,tau) & ,t(100,1,j,iabioto2,taum1) enddo endif if (n.eq.iprefo2) then do j=js,je do i=istrt,iend t(i,1,j,n,taup1)=t(i,1,j,io2,taup1) enddo enddo endif if (n.eq.iprefpo4) then do j=js,je do i=istrt,iend t(i,1,j,n,taup1)=t(i,1,j,ipo4,taup1) enddo enddo endif if (n.eq.iprefno3) then do j=js,je do i=istrt,iend t(i,1,j,n,taup1)=t(i,1,j,ino3,taup1) enddo enddo endif do j=js,je call setbcx (t(1,1,j,n,taup1), imt, km) enddo !----------------------------------------------------------------------- ! construct diagnostics associated with tracer "n" !----------------------------------------------------------------------- call diagt1 (joff, js, je, istrt, iend, n, twodt) !----------------------------------------------------------------------- ! end of tracer component "n" loop !----------------------------------------------------------------------- cao 26.4.08 end of ifdef c# if defined stat_alk c endif c# endif enddo !----------------------------------------------------------------------- ! explicit convection: adjust column if gravitationally unstable !----------------------------------------------------------------------- call convct2 (t(1,1,1,1,taup1), joff, js, je, is, ie, kmt) do j=js,je do n=1,nt call setbcx (t(1,1,j,n,taup1), imt, km) enddo enddo if (timavgperts .and. eots) then if (joff .eq. 0) nta_conv = nta_conv + 1 do j=js,je jrow = j + joff do i=istrt,iend ta_totalk(i,jrow) = ta_totalk(i,jrow) + totalk(i,j) ta_vdepth(i,jrow) = ta_vdepth(i,jrow) + vdepth(i,j) ta_pe(i,jrow) = ta_pe(i,jrow) + pe(i,j) enddo enddo endif !----------------------------------------------------------------------- ! construct diagnostics after convection !----------------------------------------------------------------------- idiag = 10 call diagt2 (joff, js, je, istrt, iend, idiag) !----------------------------------------------------------------------- ! filter tracers at high latitudes !----------------------------------------------------------------------- if (istrt .eq. 2 .and. iend .eq. imt-1) then call filt (joff, js, je) else write (stdout,'(a)') & 'Error: filtering requires is=2 and ie=imt-1 in tracer' stop '=>tracer' endif do n=1,nt do j=js,je call setbcx (t(1,1,j,n,taup1), imt, km) enddo enddo !----------------------------------------------------------------------- ! construct diagnostics after filtering (for total dT/dt) !----------------------------------------------------------------------- idiag = 1 call diagt2 (joff, js, je, istrt, iend, idiag) !----------------------------------------------------------------------- ! if needed, construct the Atmos S.B.C.(surface boundary conditions) ! averaged over this segment ! eg: SST and possibly SSS !----------------------------------------------------------------------- call asbct (joff, js, je, istrt, iend, isst, itemp) call asbct (joff, js, je, istrt, iend, isss, isalt) call asbct (joff, js, je, istrt, iend, issdic, idic) call asbct (joff, js, je, istrt, iend, issc14, ic14) call asbct (joff, js, je, istrt, iend, issalk, ialk) call asbct (joff, js, je, istrt, iend, isso2, io2) call asbct (joff, js, je, istrt, iend, issn2o, in2o) call asbct (joff, js, je, istrt, iend, issabioto2, iabioto2) call asbct (joff, js, je, istrt, iend, isspo4, ipo4) call asbct (joff, js, je, istrt, iend, issno3, ino3) !----------------------------------------------------------------------- ! calculate diagnostic delta carbon 14 !----------------------------------------------------------------------- if (tsiperts .or. timavgperts .or. snapts) then rrstd = 1000./rstd do j=js,je jrow = j + joff do k=1,km do i=istrt,iend dc14(i,k,j) = (rrstd*t(i,k,j,ic14,taup1) & /(t(i,k,j,idic,taup1) + epsln) - 1000.) & *tmask(i,k,j) enddo enddo enddo endif if (tsiperts .and. eots) then if (js+joff .eq. 2) dc14bar = 0. do j=js,je jrow = j + joff fy = cst(jrow)*dyt(jrow) do k=1,km fyz = fy*dzt(k) do i=istrt,iend dc14bar = dc14bar + dc14(i,k,j)*dxt(i)*fyz*tmask(i,k,j) enddo enddo enddo endif if (timavgperts .and. .not. euler2) then do j=js,je jrow = j + joff do k=1,km do i=istrt,iend ta_dc14(i,k,jrow) = ta_dc14(i,k,jrow) + dc14(i,k,j) enddo enddo enddo endif return end subroutine diagt1 (joff, js, je, is, ie, n, twodt) !----------------------------------------------------------------------- ! construct diagnostics associated with tracer component "n" ! input: ! joff = offset relating "j" in the MW to latitude "jrow" ! js = starting row in the MW ! je = ending row in the MW ! is = starting longitude index in the MW ! ie = ending longitude index in the MW ! n = (1,2) = (u,v) velocity component ! twodt = (2*dtts,dtts) on (leapfrog,mixing) time steps ! based on code by: R. C. Pacanowski !----------------------------------------------------------------------- include "param.h" include "accel.h" include "coord.h" include "cregin.h" include "csbc.h" include "ctavg.h" include "diag.h" include "diaga.h" include "emode.h" include "grdvar.h" include "hmixc.h" include "isopyc.h" include "levind.h" include "mw.h" include "scalar.h" include "switch.h" include "vmixc.h" dimension temp1(imt,km), temp2(imt,km), temp3(imt,km) dimension twodt(km) include "fdift.h" !----------------------------------------------------------------------- ! diagnostic: integrate |d(tracer)/dt| and tracer variance on "tau" ! globally ! based on code by: R. C. Pacanowski ! (based on diagnostic by M. Cox) !----------------------------------------------------------------------- if (tsiperts .and. eots) then do j=js,je jrow = j + joff r2dt = c1/c2dtts cosdyt = cst(jrow)*dyt(jrow) do k=1,km fx = r2dt/dtxcel(k) do i=is,ie darea = dzt(k)*dxt(i)*cosdyt*tmask(i,k,j) temp3(i,k) = t(i,k,j,n,tau)*darea temp1(i,k) = t(i,k,j,n,tau)**2*darea temp2(i,k) = abs(t(i,k,j,n,taup1)-t(i,k,j,n,taum1))* & darea*fx enddo do i=is,ie tbar(k,n,jrow) = tbar(k,n,jrow) + temp3(i,k) travar(k,n,jrow) = travar(k,n,jrow) + temp1(i,k) dtabs(k,n,jrow) = dtabs(k,n,jrow) + temp2(i,k) enddo enddo enddo endif !----------------------------------------------------------------------- ! diagnostic: accumulate tracers for averages under horizontal ! regions (use units of meters, rather than cm) ! based on code by: K. Dixon !----------------------------------------------------------------------- if (tavgts .and. eots) then do j=js,je jrow = j + joff do i=is,ie mask = mskhr(i,jrow) if (mask .ne. 0) then boxar = cst(jrow)*dxt(i)*dyt(jrow)*tmask(i,1,j)*0.0001 sumbf(mask,n) = sumbf(mask,n) + stf(i,j,n)*boxar do k=1,km sumbk(mask,k,n) = sumbk(mask,k,n) + t(i,k,j,n,tau) & *boxar*dzt(k)*tmask(i,k,j)*0.01 enddo endif enddo enddo endif !----------------------------------------------------------------------- ! diagnostic: compute the northward transport components of ! each tracer ! based on code by: R. C. Pacanowski ! (based on diagnostic by K. Bryan) !----------------------------------------------------------------------- if (gyrets .and. eots) call gyre (joff, js, je, is, ie, n) !----------------------------------------------------------------------- ! diagnostic: integrate r.h.s. terms in the tracer equations ! over specified regional volumes. ! based on code by: R. C. Pacanowski !----------------------------------------------------------------------- if (trmbts .and. eots) call ttb1 (joff, js, je, is, ie, n) return end subroutine diagt2 (joff, js, je, is, ie, idiag) !----------------------------------------------------------------------- ! construct d(tracer)/dt diagnostics ! input: ! joff = offset relating "j" in the MW to latitude "jrow" ! js = starting row in the MW ! je = ending row in the MW ! is = starting longitude index in the MW ! ie = ending longitude index in the MW ! idiag = 1 => total tracer change ! idiag = 10 => change of tracer due to filtering(also convection) !----------------------------------------------------------------------- include "param.h" include "coord.h" include "diaga.h" include "iounit.h" include "mw.h" include "scalar.h" include "switch.h" include "tmngr.h" include "timeavgs.h" !----------------------------------------------------------------------- ! diagnostic: integrate d/dt(tracer) over specified regional volumes ! after convection and filtering ! based on code by: R. C. Pacanowski !----------------------------------------------------------------------- if (trmbts .and. eots) call ttb2 (joff, js, je, is, ie, idiag) return end subroutine asbct (joff, js, je, is, ie, isbc, itr) !----------------------------------------------------------------------- ! construct the Atmos S.B.C. (surface boundary conditions) ! input: ! joff = offset relating "j" in the MW to latitude "jrow" ! js = starting row in the MW ! je = ending row in the MW ! is = starting longitude index in the MW ! ie = ending longitude index in the MW ! isbc = index for sbc ! itr = index for tracer ! based on code by: R. C. Pacanowski !----------------------------------------------------------------------- include "param.h" include "csbc.h" include "levind.h" include "mw.h" include "scalar.h" include "switch.h" ! initialize the Atmos S.B.C. at the start of each ocean segment ! (do not alter values in land) if (isbc .le. 0 .or. itr .le. 0) return if (eots .and. osegs) then do j=js,je jrow = j + joff do i=is,ie if (kmt(i,jrow) .ne. 0) sbc(i,jrow,isbc) = c0 enddo enddo endif ! accumulate surface tracers for the Atmos S.B.C. every time step if (eots) then do j=js,je jrow = j + joff do i=is,ie sbc(i,jrow,isbc) = sbc(i,jrow,isbc)+t(i,1,j,itr,taup1) enddo enddo endif ! average the surface tracers for the Atmos S.B.C. at the end of ! each ocean segment. (do not alter values in land) if (eots .and. osege) then rts = c1/ntspos do j=js,je jrow = j + joff do i=is,ie if (kmt(i,jrow) .ne. 0) & sbc(i,jrow,isbc) = rts*sbc(i,jrow,isbc) enddo enddo endif return end subroutine ivdift (joff, js, je, is, ie, n, twodt) !----------------------------------------------------------------------- ! solve vertical diffusion of tracers implicitly ! input: ! joff = offset relating "j" in the MW to latitude "jrow" ! js = starting row in the MW ! je = ending row in the MW ! is = starting longitude index in the MW ! ie = ending longitude index in the MW ! n = tracer component ! twodt = (2*dtts, dtts) on (leapfrog, mixing) time steps ! based on code by: R. C. Pacanowski !----------------------------------------------------------------------- include "param.h" include "levind.h" include "mw.h" include "switch.h" include "vmixc.h" dimension twodt(km) ! store terms to compute implicit vertical mixing on ! diagnostic time steps if (trmbts .and. eots) then do j=js,je do k=1,km do i=is,ie zzi(i,k,j) = t(i,k,j,n,taup1) enddo enddo enddo endif call invtri (t(1,1,1,n,taup1), stf(1,1,n), btf(1,1,n) &, diff_cbt(1,1,jsmw), twodt, kmt, tmask(1,1,1), is, ie &, joff, js, je) ! compute residual implicit vertical mixing if (trmbts .and. eots) then do j=js,je do k=1,km rc2dt = c1/twodt(k) do i=is,ie zzi(i,k,j) = rc2dt*(t(i,k,j,n,taup1) - zzi(i,k,j)) enddo enddo enddo endif return end