! source file: /gxfs_work1/geomar/smomw113/UVic_ESCM/2.9_dk_0814/updates/11_valley/source/mom/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 !======================================================================= implicit none character(120) :: fname, new_file_name integer istrt, iend, i, k, j, ip, kr, jq, n, jp, jrow, iou, js integer je, limit, joff, is, ie, kmx, m, kb, idiag, index integer it(10), iu(10), ib(10), ic(10), nfnpzd, mfnpzd, mxfnpzd integer id_time, id_xt, id_yt, id_zt integer fe_jlo, fe_m, fe_k, fe_n parameter (fe_n = 14) logical inqvardef, exists real rctheta, declin, gl, impo, expo, npp, time real remi, excr, graz, morp, morpt, morz, temp, swr, dayfrac real graz_Det, graz_Z, avej, avej_D, gmax, no3P, po4P, po4_D real phin, dz, prca, dprca, nud, bct, tap, fo2, so2, ai, hi, hs real npp_D, graz_D, morp_D, no3flag, deni, nfix, felimit real t_i, t_j, dz_t2r, dz_tr, dz_wtr, dx_t2r, dx_tr, dy_t2r real dy_tr, adv_tx, adv_ty, adv_tz, adv_txiso, adv_tyiso real adv_tziso, diff_tx, diff_ty, diff_tz, zmax, cont, drho real drhom1, wt, ahbi_cstr, ahbi_csu_dyur, gamma, rrstd, fy, fyz real bctz, felimit_D real Paulmier_a, Paulmier_z, Paulmier_R0, o2thresh, no3thresh # if defined TA_pipes real hpipe, wpipe # endif # if defined O_ocean_upwelling_pipes real hpipe, wpipe cwk real kpipemin real t_in, s_in, dic_in, ta_in, co2_in, pHlo, pHhi, sit_in &, pt_in, atmpres, zero, pH, co2star, dco2star,pCO2, dpco2 &, CO3, Omega_c, Omega_a, pCO2_old, kpmax, pCO2_new # if defined O_pipes_diagnostics real pCO2_orig # endif # if defined O_pipes_no_prenutrient_upward_pumping integer m_no3r,m_po4r,m1 # endif # endif real expofe, impofe, feorgads, remife, thetamax, deffe, fecol real thetachl, chl, chl_D, feprime, fesed, bfe real ofehydro cwk Ben's calibrated Fe mask runs use always subgridscale para, hence I move ! sgb from 1 above to here, with my add option real sgb include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "coord.h" include "cregin.h" include "csbc.h" include "emode.h" include "grdvar.h" include "hmixc.h" include "levind.h" include "mw.h" include "scalar.h" include "switch.h" include "timeavgs.h" include "tmngr.h" include "vmixc.h" include "diaga.h" include "ice.h" # if defined O_ocean_upwelling_pipes include "cembm.h" # endif include "atm.h" include "npzd.h" parameter (istrt=2, iend=imt-1) real twodt(km) real snpzd(ntnpzd), tnpzd(ntnpzd) real src(imt,km,jsmw:jemw,nsrc) include "isopyc.h" include "fdift.h" ! o2thresh = threshold for oxic metabolism, mmol/m3 ! This is used with fo2 to limit oxygen consumption below concentrations of 5umol/kg as recommended in OCMIP. ! With _no_sulphidic, I now compute fo2 also before npzd_src. o2thresh = 5. !----------------------------------------------------------------------- ! 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.*par ! alpha is iron dependent and calculated in npzd_src do j=js,je jrow = j + joff do i=is,ie if (kmt(i,jrow) .gt. 0) then ai = aice(i,jrow,2) hi = hice(i,jrow,2) hs = hsno(i,jrow,2) ! 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*(1. + ai*(exp(-ki*(hi + hs)) - 1.)) expo = 0.0 impo = 0.0 expofe = 0.0 impofe = 0.0 phin = 0.0 ! integrated phytoplankton prca = 0.0 ! integrated production of calcite kmx = min(kmt(i,jrow), kpzd) do k=1,kmx !----------------------------------------------------------------------- ! limit tracers to positive values !----------------------------------------------------------------------- ! { cwk 29.03.22 ! UVic version differ about the best approach where to limit tracers to ! positive values. In other UVic version this is not done here, but within ! npzd_src.F. This is related to the two (or more, parameter dependent) ! bio-loops in npzd_src.f. I stick to David's original / my version for ! long term consistency of model output. ! cwk 29.03.22 } 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) tnpzd(7) = max(t(i,k,j,idfe,taum1), trcmin) tnpzd(8) = max(t(i,k,j,idetrfe,taum1), trcmin) swr = swr*exp(-kc*phin) phin = (tnpzd(6) + tnpzd(2))*dzt(k) gl = tap*swr*exp(ztt(k)*rctheta) impo = expo*dztr(k) impofe = expofe*dztr(k) bct = bbio**(cbio*t(i,k,j,itemp,taum1)) if (t(i,k,j,itemp,taum1).gt.20) then bctz = (0.5*(tanh(t(i,k,j,io2,taum1)*1000. - 8.)+1)) & *bbio**(cbio*20) else bctz = (0.5*(tanh(t(i,k,j,io2,taum1)*1000. - 8.)+1)) & *bct endif !{ cwk 29.03.22 Currently I adopt this formulation only for the _BenPara option. ! I would need to double check otherwise, how to deal in detail with my _npzd_h2s ! code option ! cwk 29.03.22 } ! Decrease remineralisation rate in oxygen minimum zone and when there is no ! nitrate. if (t(i,k,j,ino3,taum1).le.0.) then nud = 0 else nud = nud0*(0.65+0.35*tanh(t(i,k,j,io2,taum1) & *1000.-6.)) endif !----------------------------------------------------------------------- ! call the npzd model !----------------------------------------------------------------------- call npzd_src (tnpzd, nbio(k), dtbio(k), gl, bct, impo &, dzt(k), dayfrac, wd(k), rkwz(k), nud &, snpzd, expo, graz, morp, morz &, graz_Det, graz_Z &, npp, morpt, remi, excr &, npp_D, graz_D, morp_D, nfix &, bctz &, expofe, impofe, remife &, t(i,k,j,io2,taum1) &, chl, thetachl &, chl_D &, pipeactive(i,j) & ) ! These are source/sink terms snpzd(1:4) = snpzd(1:4)*rdtts(k) snpzd(5:6) = snpzd(5:6)*rdtts(k) snpzd(7:8) = snpzd(7:8)*rdtts(k) expofe=expofe*rnbio(k) expo = expo*rnbio(k) rexpo(i,k,j) = expo rgraz(i,k,j) = graz*rnbio(k) rgraz_Det(i,k,j) = graz_Det*rnbio(k) rgraz_Z(i,k,j) = graz_Z*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) rexpofe(i,k,j) = expofe rremife(i,k,j) = remife*rnbio(k) rchl(i,k,j) = chl*rnbio(k) rchl_D(i,k,j) = chl_D*rnbio(k) ta_wd(k)=wd(k)*dzt(k)*86400./100. !----------------------------------------------------------------------- ! calculate detritus at the bottom and remineralize !----------------------------------------------------------------------- sgb = sg_bathy(i,j,k) cwk sgb is zero if no fraction of the cell in that k-level is sediment cwk for k=1, sgb is 1 on land, 0 in the open ocean cwk when it is between 0 and 1, this gives the fraction of sediment in that cell on that k-level cwk for k=19, sgb is always 0, if (sgb .gt. 0) then rremi(i,k,j) = rremi(i,k,j) + expo*sgb snpzd(1) = snpzd(1) + redptn*expo*sgb snpzd(5) = snpzd(5) + expo*sgb fesed=fetopsed*bct*redptn*expo*sgb bfe=fesed snpzd(7) = snpzd(7) + fesed if (t(i,k,j,io2,taum1)*1e3 .lt. o2min) then snpzd(7) = snpzd(7) + expofe*sgb bfe=bfe+expofe*sgb expofe = expofe - sgb * expofe endif rremife(i,k,j) = rremife(i,k,j) + bfe expo = expo - sgb * expo endif ! of (sgb .gt. 0) !----------------------------------------------------------------------- ! set major source/sink 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) if (k.gt.1) then src(i,k,j,ispo4r) = src(i,k,j,ispo4) src(i,k,j,isno3r) = src(i,k,j,isno3) endif if ( k .eq. 1 ) then snpzd(7) = snpzd(7) + ! comes in mol/m2/s gets converted to mmol/m3/s & sbc(i,j,idfeadep)*1000./(dzt(1)/100.) !dzt is in cm endif cwyao add hydrothermal flux of iron snpzd(7) = snpzd(7) & + fehydro(i,j,k) cwyao dzt is in cm rfehydro(i,k,j)=fehydro(i,j,k) src(i,k,j,isdfe) = snpzd(7) src(i,k,j,isdetrfe) = snpzd(8) ! production of calcite dprca = (morp+morz+(graz+graz_Z)*(1.-gamma1)) & *capr & *redctn*rnbio(k) prca = prca + dprca*dzt(k) src(i,k,j,isdic) = ( & snpzd(1)*redctp & - dprca) if (k.gt.1) & src(i,k,j,isdicr) = & (snpzd(1)*redctp) src(i,k,j,isalk) = ( & -snpzd(1)*redntp*1.e-3 & - 2.*dprca) if (k.gt.1) & src(i,k,j,isalkr) = & (-1)*snpzd(1)*redntp*1.e-3 !----------------------------------------------------------------------- ! 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_rgraz_Z(i,k,jrow) = ta_rgraz_Z(i,k,jrow) & + rgraz_Z(i,k,j) ta_rgraz_Det(i,k,jrow) = ta_rgraz_Det(i,k,jrow) & + rgraz_Det(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) ta_rfehydro(i,k,jrow)=ta_rfehydro(i,k,jrow) & +rfehydro(i,k,j) ta_rchl(i,k,jrow) = ta_rchl(i,k,jrow) & + rchl(i,k,j) ta_rchl_D(i,k,jrow) = ta_rchl_D(i,k,jrow) & + rchl_D(i,k,j) # if defined O_oceanmask_out ta_oceanmask(i,k,jrow) = ta_oceanmask(i,k,jrow) & + oceanmask(i,j,k) # endif endif ! calculate total export to get total import for next layer expo = expo*dzt(k) expofe = expofe*dzt(k) enddo ! k-loop which started in 420 !----------------------------------------------------------------------- ! set o2 sinks-minus-sources !----------------------------------------------------------------------- kmx = kmt(i,jrow) do k=1,kmx fo2 = 0.5*tanh(t(i,k,j,io2,taum1)*1000. - o2thresh) ! sink of oxygen (so2) ! 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) ! last endif: of _respi_data & _respi_data_transient (l 749) ! last endif: of defined _respi_data & defined _respi_data_repyr if (k.gt.1) & src(i,k,j,iso2r) = src(i,k,j,iso2) !----------------------------------------------------------------------- ! calculate denitrification, its effect on alkalinity and ideal. tracers !----------------------------------------------------------------------- ! 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) deni = max(0., 800.*no3flag*so2*(0.5 - fo2)) src(i,k,j,isno3) = src(i,k,j,isno3) - deni rdeni(i,k,jrow) = deni if (k.gt.1) & src(i,k,j,isno3d) = (-1.) * deni !----------------------------------------------------------------------- ! alkalinity corrections: denitrification, n2-fixation, sulphate reduction, etc !----------------------------------------------------------------------- ! Correct the alkalinty effect of denitrification based on the stoichiometric model ! formulated in Paulmier et al. 2009 BG cwk In the following I added more comments in order to better understand what is done here. cwk In particularly, it was confusing that the code deviates from equation given in Paulmier. cwk However, the computed R0 values nevertheless agree! This is documented below. cwk For some unexplained reason the term Paulmier_a is here the C:N ratio instead of the C:P ratio Paulmier_a = 1.e3*redctn cwk Excess hydrogen, z, for complete respiration with NO3 as N-endproduct; Eq 27 in Paulmier cwk z = 4 * rO2:P - 4a -8d, with rO2:P = R0 + 2d; a is the C:P ratio of organic matter, d is the N:P ratio cwk For some unexplained reason Paulmier_z uses the C:N ratio instead of the C:P ratio?? Paulmier_z = 4*1.e3*redotp - 4.*Paulmier_a - 8.*redntp cwk Compute R0 in the absence of knowing the hydrogen and oxygen contend of organic matter, using cwk equation 26 of Paulmier cwk Given typical UVic parameters, redctn=6.625, redctp=106, redntp =16, redotp=160 cwk equation 27 gives z=640 - 424 - 128 = 88 cwk equation Paulmier_z = 640 - 26.5 - 128 = 485.5 Paulmier_R0 = Paulmier_a + 0.25*Paulmier_z cwk R0 based on z from equ 27 gives: 106 + 22 = 128 cwk R0 based on Paulmier_z = 6.625 + 121.375 = 128 src(i,k,j,isalk) = src(i,k,j,isalk) & + no3flag*src(i,k,j,ispo4)*(0.5-fo2) cwk Paulmier et al equation from their table 1 is: 4/5 R0 + 3/5d -1 cwk The term -1 reflects the alkalinity effect of P (which is ignored in Uvic) cwk The term (3/5+1)d instead of 3/5d backcorrects the already considers nitrate effect on alk cwk (src(i,k,j,isalk) = -snpzd(1)*redntp*1.e-3), see above & *(4./5.*Paulmier_R0 + (3./5. +1.)*redntp) * 1.e-3 cwk gives 0.8*128 + (0.6+1)*16 = 102.4 + 25.6 = 128 ! Now account for N2 fixation (ALK production is tied to PO4 change and ! thus in the case of N2 fixation was not correct without this fix). src(i,k,j,isalk) = src(i,k,j,isalk) & - rnfix(i,k,j)*1.e-3 if (k.gt.1) & src(i,k,j,isalkr) = src(i,k,j,isalkr) ! next two lines include the effect of denitrification & + no3flag*src(i,k,j,ispo4)*(0.5-fo2) & *(4./5.*Paulmier_R0 + (3./5. +1.)*redntp) * 1.e-3 ! next line corrects for the effect of N2-fix (hardly important for k.gt.1 & - rnfix(i,k,j)*1.e-3 ! last endif: of _n2o (line 1345) ! last endif: of _npzd_nitrogen (l 875) enddo ! of do over k (line 725), beginning of block: set o2 sinks-minus-sources ! last endif: of _npzd_o2 (l 725) !----------------------------------------------------------------------- ! remineralize calcite !----------------------------------------------------------------------- kmx = kmt(i,jrow) do k=1,kmx-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) if (k.eq.1) & src(i,k,j,isalks) = 0. if (k.gt.1) & src(i,k,j,isalks) = & 2.*prca*rcak(k) enddo src(i,kmx,j,isdic) = src(i,kmx,j,isdic) & + prca*rcab(kmx) src(i,kmx,j,isalk) = src(i,kmx,j,isalk) & + 2.*prca*rcab(kmx) if (kmx.gt.1) & src(i,kmx,j,isalks) = 0. & + 2.*prca*rcab(kmx) ! wk, 20.11.19, bottom water calcite degradation error corrected ! wk, 08.01.20, additional issue with bottom water calcite dissolution corrected ! The following line numbers are outdated due to new code additions. ! Also with all the assimilation code there are several instances where the respective ! variables are modified. This always needs a close inspection by checking code/tracer.f ! for src(i,k,j,iso2r) see code in line 1218ff ! for src(i,k,j,ispo4r) see code in line 879ff ! for src(i,k,j,isno3r) see code in line 880ff ! for src(i,k,j,isdicr) see code in line 907, 1101, 1131, 1169 ! for src(i,k,j,isalkr) see code in line 926, 1105, 1135, 1172, 1334 ! for src(i,k,j,isalks) see code below in line 656ff ! for src(i,k,j,isno3d) see code in line 614ff !----------------------------------------------------------------------- ! accumulate time averages for full depth variables !----------------------------------------------------------------------- if (timavgperts .and. .not. euler2) then kmx = kmt(i,jrow) expo = prca ta_rprocal(i,jrow) = ta_rprocal(i,jrow) + prca do k=1,kmx expo = expo*dztr(k) 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) cwk 14.07.22 next two lines be aware that below the variable expo is reasigned (i.e. f. CaCO3 export) cwk this is very old ('dirty', I think) code expo = expo - prca*rcak(k) ta_rexpocal(i,k,jrow) = ta_rexpocal(i,k,jrow) + expo expo = expo*dzt(k) ta_rdeni(i,k,jrow) = ta_rdeni(i,k,jrow) + rdeni(i,k,j) ta_rremife(i,k,jrow) = ta_rremife(i,k,jrow) & + rremife(i,k,j) ta_rexpofe(i,k,jrow) = ta_rexpofe(i,k,jrow) & + rexpofe(i,k,j) enddo endif ! last endif: of _time_averages && _save_npzd (line 1175) endif enddo enddo ! last endif: of _npzd (line 404 ! last endif: of _carbon_fnpzd (line 1225) cwk 4.01.19 I moved the code for age change of the ideal age tracer to here cwk i.e. outside 1 cwk 12.03.19, however, I did not include the two outer do j, do i loops! cwk added now !----------------------------------------------------------------------- ! increment ideal age !----------------------------------------------------------------------- do j=js,je jrow = j + joff do i=is,ie if (kmt(i,jrow) .gt. 0) then kmx = kmt(i,jrow) do k=1,kmx src(i,k,j,isidealage) = 1.0 / 3.6e3 /24. enddo endif enddo enddo !----------------------------------------------------------------------- cwk 4.01.09 BGC sms are done by now, compute the transport terms: !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! solve for one tracer at a time ! n = 1 => temperature ! n = 2 => salinity ! n > 2 => other tracers (if applicable) !----------------------------------------------------------------------- do n=1,nt !----------------------------------------------------------------------- ! calculate advective tracer flux !----------------------------------------------------------------------- call adv_flux (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 ! last endif: of _isopycmix ! last endif: of !_biharmonic || _bryan_lewis_horizontal ! lase endif of _consthmix !----------------------------------------------------------------------- ! 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 # if defined O_ocean_upwelling_pipes !----------------------------------------------------------------------- ! artificial ocean upwelling via pipes !----------------------------------------------------------------------- ! This is the pipe option in the 2.8 CE, originally from Oschlies et al. 2010, ! as modified by Keller et al., 2014. WK made additional changes to this code ! in 2021, 2022 etc. do j=js,je do i=istrt,iend hpipe = 0. cwk pipespeed=-1. wpipe = pipespeed/86400 cwk kpipemin=2 cwk kpipemax=8 pipeactive(i,j)=0. # if defined O_pipes_diagnostics delpCO2_soft(i,j) = 0. delpCO2_solu(i,j) = 0. # endif # if defined O_pipes_at_kmax kpipe(i,j) = kpipemax # else kpipe(i,j) = 1 # endif # if defined O_pipes_everywhere ! if (kmt(i,j).ge.kpipemin.and.t(i,1,j,ipo4,tau).lt.1000) # else if (kmt(i,j).ge.kpipemin # if defined O_oceanmask & defined O_pipes_mask & !defined O_pipes_mask_const cwk Experimental code, 08.02.21: The oceanmask has been read-in in setmom.f cwk However, it is only used to prescribe the area, where pipes are applied, not the specific kpipe cwk 19.10.22 Note, that the second condition is needed since oceanmask has now a high positive missing_value & .and.nint(oceanmask(i,j,1)).ge.kpipemin & .and.nint(oceanmask (i,j,1)).le.kpipemax & .and.t(i,1,j,ipo4,tau).lt.pipeoligo # else & .and.t(i,1,j,ipo4,tau).lt.pipeoligo # endif & ) then t_in = t(i,1,j,1,tau) !degree Celsius s_in = 1000.*t(i,1,j,2,tau) + 35. !psu dic_in = t(i,1,j,idic,tau) !mol/m^3 ta_in = t(i,1,j,ialk,tau) !eq/m^3 # if defined O_carbon_co2_2d co2_in = at(i,j,2,ico2) # else co2_in = co2ccn # endif pHlo = 6. pHhi = 10. sit_in = 7.6875e-03 !mol/m^3 pt_in = 0.5125e-3 !mol/m^3 atmpres = 1.0 !atm zero = 0. call co2calc_SWS (t_in, s_in, dic_in, ta_in, co2_in &, pt_in, sit_in, atmpres, zero, pHlo, pHhi, pH &, co2star, dco2star, pCO2, dpco2, CO3 &, Omega_c, Omega_a) pCO2_old = pCO2 # if defined O_pipes_diagnostics pCO2_orig = pCO2 # endif kpmax=min(INT(kpipemax),kmt(i,j)) # if defined O_pipes_at_kmax kpipe(i,j) = kpipemax # else kpipe(i,j) = 1 # endif do k=2,kpmax t_in = t(i,1,j,1,tau) s_in = 1000.*t(i,k,j,2,tau) + 35. dic_in = t(i,k,j,idic,tau) & - (t(i,k,j,ipo4,tau)-t(i,1,j,ipo4,tau))*16. # if defined O_pipes_rmodcondi & *redctnefac*6.7/1000. # else & *6.7/1000. # endif ta_in = t(i,k,j,ialk,tau) call co2calc_SWS (t_in, s_in, dic_in, ta_in, co2_in &, pt_in, sit_in, atmpres, zero, pHlo, pHhi, pH &, co2star, dco2star, pCO2, dpco2, CO3 &, Omega_c, Omega_a) pCO2_new = pCO2 if (pCO2_new.lt.pCO2_old) then pCO2_old = pCO2_new # endif ! last endif: !_pipes_everywhere # if defined O_pipes_at_kmax if(kmt(i,j).gt.kpipemax)then kpipe(i,j) = kpipemax else kpipe(i,j) = max(2,kmt(i,j)) endif # else # if defined O_pipes_everywhere kpmax=min(kpipemax,kmt(i,j)) do k=2,kpmax # endif if(kmt(i,j).ge.k.and.kmt(i,j).ge.kpipemin)then kpipe(i,j) = k endif # endif # if defined O_ocean_upwelling_pipes_with_fe kpipe_fe(i,j)=kpipe(i,j) # endif # if defined O_pipes_diagnostics delpCO2_kpipe(i,j)=pCO2_new - pCO2_orig # endif # if defined O_pipes_everywhere # else # if defined O_pipes_at_kmin cwk 9.02.21: using the _at_kmin option to break here: pCO2_new was found to lower than pCO2_old EXIT # endif endif # endif enddo # if defined O_oceanmask & defined O_pipes_mask & defined O_pipes_mask_const cwk The selection of the most appropriate kpipe(i,j) done above is being ignored. cwk Instead, I read in a prestribed time-constant kpipe value from the oceanmask variable. if ((nint(oceanmask(i,j,1).ge.kpipemin)).and. & (nint(oceanmask(i,j,1).le.kpipemax))) kpipe(i,j) = nint(oceanmask(i,j,1)) # if defined O_ocean_upwelling_pipes_with_fe kpipe_fe(i,j)=kpipe(i,j) # endif # endif if (kpipe(i,j).ge.2) then # if defined O_pipes_diagnostics cwk 05.02.21 experimental code by wk cwk all property definitions (e.g. where I use k or 1) used here are subject to discussion cwk delpCO2 due to soft tissue pump cwk for t_in, s_in, dic_in, ta_in: k=1 for this purpose cwk for pumped po4, I only consider po4pre t_in = t(i,1,j,1,tau) s_in = 1000.*t(i,1,j,2,tau) + 35. dic_in = t(i,1,j,idic,tau) & - (t(i,kpipe(i,j),j,iprepo4,tau)-t(i,1,j,ipo4,tau))*16.*6.7/1000. ta_in = t(i,1,j,ialk,tau) call co2calc_SWS (t_in, s_in, dic_in, ta_in, co2_in &, pt_in, sit_in, atmpres, zero, pHlo, pHhi, pH &, co2star, dco2star, pCO2, dpco2, CO3 &, Omega_c, Omega_a) ckw delpCO2_ is not del against the atm, but pCO2_orig delpCO2_soft(i,j) = pCO2 - pCO2_orig cwk delpCO2 due to solubility pump cwk for s_in, dic_in, ta_in I use k=kpipe(i,j) cwk for alk, dic, I only consider to 'pump' the preformed components t_in = t(i,1,j,1,tau) s_in = 1000.*t(i,kpipe(i,j),j,2,tau) + 35. dic_in = t(i,kpipe(i,j),j,ipredic,tau) ta_in = t(i,kpipe(i,j),j,iprealk,tau) call co2calc_SWS (t_in, s_in, dic_in, ta_in, co2_in &, pt_in, sit_in, atmpres, zero, pHlo, pHhi, pH &, co2star, dco2star, pCO2, dpco2, CO3 &, Omega_c, Omega_a) delpCO2_solu(i,j) = pCO2 - pCO2_orig # endif pipeactive(i,j)=1. # if defined O_pipes_no_heat_exchange if(n.gt.1) then # endif # if defined O_pipes_no_nutrient_exchange cwk 24.09.19 (?), allow for pipes that do not pump nutrients if ((trim(mapt(n)) .ne. 'no3') & .and. (trim(mapt(n)) .ne. 'po4') & .and. (trim(mapt(n)) .ne. 'preno3') & .and. (trim(mapt(n)) .ne. 'prepo4') & .and. (trim(mapt(n)) .ne. 'no3r') & .and. (trim(mapt(n)) .ne. 'po4r') & .and. (trim(mapt(n)) .ne. 'no3d') & ) then # endif # if defined O_pipes_no_alk_exchange if ((trim(mapt(n)) .ne. 'alk') & .and. (trim(mapt(n)) .ne. 'prealk') & .and. (trim(mapt(n)) .ne. 'alkr') & .and. (trim(mapt(n)) .ne. 'alks') & ) then # endif # if defined O_pipes_no_dic_exchange if ((trim(mapt(n)) .ne. 'dic') & .and. (trim(mapt(n)) .ne. 'predic') & .and. (trim(mapt(n)) .ne. 'dicr') & ) then # endif # if defined O_pipes_no_prenutrient_upward_pumping cwk Make sure nitrate and phosphate are treated differently, cwk i.e. not no3 (po4), but no3r (po4r) are being pumped UP cwk Note: this option should be used only in oligotrophic runs, cwk I think. (wk 25.11.21) cwk Note2: dont combile this option with, e.g. _no_alk_exchange cwk this may not work. if (trim(mapt(n)) .eq. 'no3') then do m1=1,nt if (trim(mapt(n)) .eq. 'no3r') then m_no3r = m1 endif enddo cwk conditional added added 26.11.21) if ( t(i,kpipe(i,j),j,m_no3r,taum1).gt. & t(i,1,j,n,taum1) ) then source(i,1,j) = source(i,1,j) & -wpipe*(t(i,kpipe(i,j),j,m_no3r,taum1) & -t(i,1,j,n,taum1))/dzt(1) endif else if (trim(mapt(n)) .eq. 'po4') then do m1=1,nt if (trim(mapt(n)) .eq. 'po4r') then m_po4r = m1 endif enddo cwk conditional added added 26.11.21) if ( t(i,kpipe(i,j),j,m_po4r,taum1).gt. & t(i,1,j,n,taum1) ) then source(i,1,j) = source(i,1,j) & -wpipe*(t(i,kpipe(i,j),j,m_po4r,taum1) & -t(i,1,j,n,taum1))/dzt(1) endif else cwk This is for all other tracers, the normal code up & down source(i,1,j) = source(i,1,j) & -wpipe*(t(i,kpipe(i,j),j,n,taum1) & -t(i,1,j,n,taum1))/dzt(1) endif cwk Note that the downward pumping is not affected by do k=2,kpipe(i,j) source(i,k,j) = source(i,k,j) & -wpipe*(t(i,k-1,j,n,taum1) & -t(i,k,j,n,taum1))/dzt(k) enddo # else cwk The next 8 lines is the original code, which will be cwk compiled if no specific condition is given source(i,1,j) = source(i,1,j) & -wpipe*(t(i,kpipe(i,j),j,n,taum1) & -t(i,1,j,n,taum1))/dzt(1) do k=2,kpipe(i,j) source(i,k,j) = source(i,k,j) & -wpipe*(t(i,k-1,j,n,taum1) & -t(i,k,j,n,taum1))/dzt(k) enddo # endif # if defined O_pipes_no_dic_exchange endif # endif # if defined O_pipes_no_alk_exchange endif # endif # if defined O_pipes_no_nutrient_exchange endif # endif # if defined O_pipes_no_heat_exchange endif # endif # if defined O_pipes_everywhere # else endif # endif else # if defined O_pipes_at_kmax cwk 03.02.21 Here I correct a bug in the original code, which affected the cwk diagnostic O_zpipe, which however could also affect the run cwk time model behaviour when the Fe pipe is on. cwk We are in eutropic waters; make sure kpipe is no longer wrongly set to 8 (kpipemax) kpipe(i,j) = 1 # endif endif enddo enddo # endif !----------------------------------------------------------------------- ! Acid Downwelling via Pipes (mjuerchott, 13.10.2024) !----------------------------------------------------------------------- ! This part of the code does surface ocean alkalinity enhancement, ! adds an equal amount of acid to depth and pumps "carier-water" to depth # if defined TA_pipes do j=js,je do i=istrt,iend hpipe = 0. cwk pipespeed=-0.000503. ! 1000m, k=8: -0.0004 (Tyka et al., 2022 -> 0.015 Sv) ! 2000m, k=10: -0.00042 ! 3000m, k=13: -0.000503 wpipe = pipespeed/86400 cwk kpipemax=13 pipeactive(i,j)=0. if (kmt(i,j).ge.kpipemax) then pipeactive(i,j)=1. kpipe(i,j) = kpipemax # if defined TA_OAE_surface if ((trim(mapt(n)) .eq. 'alk') & .or. (trim(mapt(n)) .eq. 'prealk') & .or. (trim(mapt(n)) .eq. 'alkr') & .or. (trim(mapt(n)) .eq. 'alks') & ) then source(i,1,j) = source(i,1,j) & +0.0000000006154 endif ! 1000m, k=8: +0.00000000048 ! 2000m, k=10: +0.000000000513 ! 3000m, k=13: +0.0000000006154 # endif # if defined TA_acid_to_depth if ((trim(mapt(n)) .eq. 'alk') & .or. (trim(mapt(n)) .eq. 'prealk') & .or. (trim(mapt(n)) .eq. 'alkr') & .or. (trim(mapt(n)) .eq. 'alks') & ) then source(i,kpipemax,j) = source(i,kpipemax,j) & -0.00000000007494 endif ! 1000m, k=8: -0.0000000000922 ! 2000m, k=10: -0.00000000008005 ! 3000m, k=13: -0.00000000007494 # endif # if defined TA_downwelling_pipes if ((trim(mapt(n)) .ne. 'alk') & .or. (trim(mapt(n)) .ne. 'prealk') & .or. (trim(mapt(n)) .ne. 'alkr') & .or. (trim(mapt(n)) .ne. 'alks') & ) then source(i,kpipemax,j) = source(i,kpipemax,j) & -wpipe*(t(i,pipeactive(i,j),j,n,taum1) & -t(i,kpipemax,j,n,taum1))/dzt(kpipemax) do k=1,kpipe(i,j)-1 source(i,k,j) = source(i,k,j) & -wpipe*(t(i,k+1,j,n,taum1) & -t(i,k,j,n,taum1))/dzt(k) enddo endif # endif endif enddo enddo # else do j=js,je do i=istrt,iend pipeactive(i,j)=0.0 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) enddo enddo enddo ! 2nd: add in portion of vertical diffusion handled implicitly call ivdift (joff, js, je, istrt, iend, n, twodt) 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 !----------------------------------------------------------------------- enddo ! of do n=1,nt (line 1400), i.e. 'solve for one tracer at a time' !----------------------------------------------------------------------- cwk !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! 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) ! ---------------------------------------------------------------------- ! restore t(iidealage,taup1) to 0.0 in the surface layer ! ! ---------------------------------------------------------------------- do j=js,je ! I dont really use the next line ..... but the output looks fine jrow = j + joff do i=istrt,iend t(i,1,j,iidealage,taup1)= 0.0 enddo enddo ! ---------------------------------------------------------------------- ! restore the remineralised tracer in the surface layer to zero ! ! ---------------------------------------------------------------------- twodt(1) = c2dtts*dtxcel(1) do j=js,je do i=istrt,iend t(i,1,j,io2r,taup1) = 0.0 t(i,1,j,ipo4r,taup1) = 0.0 t(i,1,j,ino3r,taup1) = 0.0 t(i,1,j,idicr,taup1) = 0.0 t(i,1,j,ialkr,taup1) = 0.0 t(i,1,j,ialks,taup1) = 0.0 t(i,1,j,ino3d,taup1) = 0.0 enddo enddo ! ---------------------------------------------------------------------- ! restore preformed tracers in the surface layer to the concentration ! of the respective reference tracer ! ---------------------------------------------------------------------- do j=js,je do i=istrt,iend t(i,1,j,ipreo2,taup1)= t(i,1,j,io2,taup1) t(i,1,j,iprepo4,taup1)= t(i,1,j,ipo4,taup1) t(i,1,j,ipreno3,taup1)= t(i,1,j,ino3,taup1) t(i,1,j,ipredic,taup1)= t(i,1,j,idic,taup1) t(i,1,j,iprealk,taup1)= t(i,1,j,ialk,taup1) enddo enddo !----------------------------------------------------------------------- ! 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, issalk, ialk) call asbct (joff, js, je, istrt, iend, isso2, io2) call asbct (joff, js, je, istrt, iend, isso2abiot &, io2abiot) call asbct (joff, js, je, istrt, iend, isspo4, ipo4) call asbct (joff, js, je, istrt, iend, issno3, ino3) call asbct (joff, js, je, istrt, iend, issdfe, idfe) !----------------------------------------------------------------------- ! calculate diagnostic pipe depth for the _ocean_upwelling_pipes CE option !----------------------------------------------------------------------- if (timavgperts .and. .not. euler2) then do j=js,je jrow = j + joff do i=istrt,iend if (kpipe(i,jrow).ge.2) then ta_kpipe(i,jrow) = ta_kpipe(i,jrow) + zt(kpipe(i,j)) ta_pipeactive(i,jrow) =ta_pipeactive(i,jrow) & + pipeactive(i,j) ta_kpipek(i,jrow) = ta_kpipek(i,jrow) & + kpipe(i,j) endif enddo enddo endif return end subroutine get_co2ccn(co2ccn_out) 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 !----------------------------------------------------------------------- implicit none integer i, k, j, ip, kr, jq, n, jp, jrow, js, je, joff, is, ie integer mask, m real t_i, t_j, dz_t2r, dz_tr, dz_wtr, dx_t2r, dx_tr, dy_t2r real dy_tr, adv_tx, adv_ty, adv_tz, adv_txiso, adv_tyiso real adv_tziso, diff_tx, diff_ty, diff_tz, dtdx, dtdy, dtdz real r2dt, cosdyt, fx, darea, boxar, rtwodt, sumdx, delx real sumdxr, dxdy, dxdydz include "size.h" include "param.h" include "pconst.h" include "stdunits.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 "levind.h" include "mw.h" include "scalar.h" include "switch.h" include "vmixc.h" real temp1(imt,km), temp2(imt,km), temp3(imt,km) real twodt(km) include "isopyc.h" include "fdift.h" !----------------------------------------------------------------------- ! diagnostic: integrate |d(tracer)/dt| and tracer variance on "tau" ! globally !----------------------------------------------------------------------- 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) !----------------------------------------------------------------------- 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 !----------------------------------------------------------------------- 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. !----------------------------------------------------------------------- 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) !----------------------------------------------------------------------- implicit none integer idiag, j, js, je, k, i, joff, iocv, jrow, is, ie real rdt, reltim, period include "size.h" include "param.h" include "pconst.h" include "stdunits.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 !----------------------------------------------------------------------- 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 !----------------------------------------------------------------------- implicit none integer isbc, itr, j, js, je, jrow, joff, i, is, ie real rts include "size.h" include "param.h" include "pconst.h" include "stdunits.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 !----------------------------------------------------------------------- implicit none integer j, js, je, k, i, is, ie, n, joff real rc2dt include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "levind.h" include "mw.h" include "switch.h" include "vmixc.h" real 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 subroutine swflux0 (joff, js, je, is, ie, source) return end !