! source file: /gxfs_work1/fs1/work-geomar6/smomw258/Kvale_etal_2020_GMD/model_output_for_paper/updates/tracer.F subroutine tracer (joff, js, je, is, ie) !#if defined 1 !======================================================================= ! 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, 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_loc, dz, prca, dprca, nud, bct, tap, fo2, so2, ai real 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 fe_dy, fe_conc, fe_x(fe_n), fe_y(fe_n), bctz, felimit_D real Paulmier_a, Paulmier_z, Paulmier_R0 real expofe, impofe, feorgads, remife, thetamax, deffe, fecol real thetachl, chl, chl_D, feprime, fesed, bfe, sgb real sil_in, feorgads_ca real remi_B,expo_B, impo_B, graz_Det_B real npp_Diat, graz_Diat, morp_Diat, morpt_Diat,felimit_Diat real deffe_Diat, thetamax_Diat real dprop, prop, expo_op, rivsil, op_diss, bsi real npp_C, graz_C, morp_C, morpt_C,felimit_C real deffe_C, thetamax_C real bdeni, din15flag, lno3, sg_bdeni real temp_in, d_in,alk_in, sal_in, expocaco3,old_diss real dissl,calpro,impocaco3,calatt real expo_out, impo_out, caco3in_loc, dissk1 real atmpres1, pHlo1, pHhi1, pH1, p_in real co2star1, dco2star1, pCO21, dpco21, CO31 real omegaca, omegaar, c_in,depth1 include "cembm.h" 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" include "atm.h" include "npzd.h" parameter (istrt=2, iend=imt-1) real twodt(km) real snpzd(ntnpzd), tnpzd(ntnpzd), gbio real settle, constant, temp_2(km) real bctd(imt,km,jsmw:jemw) real src(imt,km,jsmw:jemw,nsrc) include "isopyc.h" include "fdift.h" !----------------------------------------------------------------------- ! 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 added in npzd_src bsi = 0.0 rivsil = 0.0 !river Si deposition mol si/s?? silwflx = 0.0 !Opal flux through bottom umol si/s bctd(:,:,:) = 0.0 !temperature scaling of microbial dissolution 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 expo_B = 0.0 ! detritus from coccs and zoop impo_B = 0.0 ! detritus from coccs and zoop remi_B = 0.0 ! detritus shifting ballasted and free pools expo_op = 0.0 ! export of opal umol Si/cm^2/s dprop = 0.0 ! umol Si/cm^3/s prop = 0.0 ! integrated production of opal mol Si/m^2/s op_diss = 0.0 ! opal dissolution rate umol/cm^2/s ?? caco3in_loc = 0.0 old_diss = 0.0 calpro = 0.0 calatt = 0.0 dissl = 0.0 impocaco3 = 0.0 expocaco3 = 0.0 expo_out = 0.0 impo_out = 0.0 dissk1 = 0.0 expofe = 0.0 impofe = 0.0 phin_loc = 0.0 ! integrated phytoplankton prca = 0.0 ! integrated production of calcite sg_bdeni = 0. bdeni = 0. kmx = min(kmt(i,jrow), kpzd) do k=1,kmx !----------------------------------------------------------------------- ! Initialize tracers. !----------------------------------------------------------------------- src(i,k,j,:)=0.0 tnpzd(ibion) = t(i,k,j,ipo4,taum1) tnpzd(ibiop) = t(i,k,j,iphyt,taum1) tnpzd(ibioz) = t(i,k,j,izoop,taum1) tnpzd(ibiod) = t(i,k,j,idetr,taum1) tnpzd(ibiono3) = t(i,k,j,ino3,taum1) tnpzd(ibiodiaz) = t(i,k,j,idiaz,taum1) tnpzd(ibiodfe) = t(i,k,j,idfe,taum1) tnpzd(ibiodetrfe) = t(i,k,j,idetrfe,taum1) tnpzd(ibioc) = t(i,k,j,icocc,taum1) tnpzd(ibiodiat) = t(i,k,j,idiat,taum1) tnpzd(ibiocaco3) = t(i,k,j,icaco3,taum1) ! mmol/m3 tnpzd(ibiod_B) = t(i,k,j,idetr_B,taum1) swr = swr*exp(-(kc*phin_loc+kc_c*caco3in_loc)) ! phin_loc has been removed from the following code to fix a bug phin_loc = (max(tnpzd(ibiop),trcmin) & + max(tnpzd(ibiodiaz),trcmin) & + max(tnpzd(ibioc),trcmin) & + max(tnpzd(ibiodiat),trcmin) & )*dzt(k) caco3in_loc = max(tnpzd(ibiocaco3),trcmin)*dzt(k) impocaco3 = expocaco3*dztr(k) impo_B = expo_B*dztr(k) !ballasted detritus gl = tap*swr*exp(ztt(k)*rctheta) impo = expo*dztr(k) impofe = expofe*dztr(k) bct = bbio**(cbio*t(i,k,j,itemp,taum1)) bctd(i,k,j) = bct !store it for using below 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 ! 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 !----------------------------------------------------------------------- ! set co2 concentration or emissions by tracking average co2 !----------------------------------------------------------------------- call co2ccndata !----------------------------------------------------------------------- ! calculate delCO3 !----------------------------------------------------------------------- sal_in = 1000.0*t(i,k,j,isalt,taum1) + 35.0 sal_in = min(45.,max(sal_in,0.)) temp_in = t(i,k,j,itemp,taum1) temp_in = min(35.,max(temp_in,-2.)) d_in = max(t(i,k,j,idic,taum1),trcmin) alk_in = max(t(i,k,j,ialk,taum1),trcmin) pHlo1 = 6. pHhi1 = 10. sil_in = max(t(i,k,j,isil,taum1),trcmin) !silicate mol/m^3 p_in = 0.5125e-3 !phos mol/m^3 atmpres1 = 1.0 !atm depth1 = dzt(k)/100 !depth in meters c_in = co2ccn call co2calc_SWS (temp_in, sal_in,d_in,alk_in,c_in, p_in &, sil_in, atmpres1, depth1, pHlo1, pHhi1,pH1 &, co2star1, dco2star1, pCO21, dpco21, CO31 &, omegaca, omegaar) romca(i,k,j) = omegaca ! assume KK that Ca ion is constant (S&G 2006 pg 366) rco3(i,k,j) = CO31 !mol/m3 rco3_sat(i,k,j) = CO31/omegaca rdel_sat(i,k,j) = rco3(i,k,j)- rco3_sat(i,k,j) !mol/m3 dissk0 = (1-(CO31-CO31/omegaca))/(kcal+(CO31-CO31/omegaca)) ! dissk0 and dissk1 dimensionless if(dissk0.ge.0) then dissk1 = min(c1,dissk0) else dissk1 = 0.0 endif ! if(omegaca-1.lt.0) then ! dissk1 = dissk0*(1-omegaca) !Andreas S. 11 d^-1 from Gehlen 2006 ! else ! dissk1 = 0. ! endif sil_in = max(t(i,k,j,isil,taum1)*1e3,trcmin) !convert to mmol m-3 for npzd_src! k1si = 0.8 + 7.2 * (sil_in/(30.+sil_in)) !Aumont et al (2003) in mmol m-3 ! repeated down below unfortunately, needs tidying up sipr = 0. ! From HAMOCC5 (Aumont et al, GBC, 2003) values 0-0.5, mostly ~0.1 sipr = sipr0 * min(1.,sil_in/k1si) & * (4.-3.*min(1., & max(t(i,k,j,idfe,taum1),trcmin) &/kfe_Diat)) gbio = gbio0 !----------------------------------------------------------------------- ! call the npzd model !----------------------------------------------------------------------- call npzd_src (tnpzd, nbio(k), dtbio(k), gl, bct, impo &, dzt(k), dayfrac, wd(k), gbio &, rkwz(k), nud &, snpzd, expo, graz, morp, morz, graz_Det &, graz_Z &, npp, morpt, remi, excr &, impo_B, expo_B, remi_B, graz_Det_B &, npp_D, graz_D, morp_D, nfix &, npp_C, morpt_C, graz_C, morp_C &, thetamax_C, deffe_C &, npp_Diat, morpt_Diat, graz_Diat, morp_Diat &, thetamax_Diat, deffe_Diat &, sil_in &, impocaco3,wc(k),expocaco3,dissl &, calpro &, calatt,dissk1 &, avej, avej_D, gmax, no3P, po4P, po4_D &, bctz &, expofe, impofe, remife &, t(i,k,j,io2,taum1) & ) ! These are source/sink terms snpzd(ibion) = snpzd(ibion)*rdtts(k) snpzd(ibiop) = snpzd(ibiop)*rdtts(k) snpzd(ibioz) = snpzd(ibioz)*rdtts(k) snpzd(ibiod) = snpzd(ibiod)*rdtts(k) snpzd(ibiono3) = snpzd(ibiono3)*rdtts(k) snpzd(ibiodiaz) = snpzd(ibiodiaz)*rdtts(k) snpzd(ibioc) = snpzd(ibioc)*rdtts(k) snpzd(ibiodiat) = snpzd(ibiodiat)*rdtts(k) expocaco3 = expocaco3*rnbio(k) snpzd(ibiocaco3) = snpzd(ibiocaco3)*rdtts(k) snpzd(ibiod_B) = snpzd(ibiod_B)*rdtts(k) expo_B = expo_B*rnbio(k) !ballast detritus snpzd(ibiodfe) = snpzd(ibiodfe)*rdtts(k) snpzd(ibiodetrfe) = snpzd(ibiodetrfe)*rdtts(k) expofe=expofe*rnbio(k) expo = expo*rnbio(k) 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) rgraz_Det_B(i,k,j) = graz_Det_B*rnbio(k) rremi_B(i,k,j) = remi_B*rdtts(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) rcalpro(i,k,j) = calpro*rnbio(k) rcalatt(i,k,j) = calatt*rnbio(k) rdissl(i,k,j) = dissl*rdtts(k) rnpp_Diat(i,k,j) = npp_Diat*rnbio(k) rgraz_Diat(i,k,j) = graz_Diat*rnbio(k) rmorp_Diat(i,k,j) = morp_Diat*rnbio(k) rmorpt_Diat(i,k,j) = morpt_Diat*rnbio(k) rnpp_C(i,k,j) = npp_C*rnbio(k) rgraz_C(i,k,j) = graz_C*rnbio(k) rmorp_C(i,k,j) = morp_C*rnbio(k) rmorpt_C(i,k,j) = morpt_C*rnbio(k) rremife(i,k,j) = remife*rnbio(k) ravej(i,k,j) = avej*rnbio(k) ravej_D(i,k,j) = avej_D*rnbio(k) rgmax(i,k,j) = gmax*rnbio(k) rno3P(i,k,j) = no3P*rnbio(k) rpo4P(i,k,j) = po4P*rnbio(k) rpo4_D(i,k,j) = po4_D*rnbio(k) !----------------------------------------------------------------------- ! calculate detritus at the bottom and remineralize !----------------------------------------------------------------------- sgb = sg_bathy(i,j,k) if (sgb .gt. 0) then rremi(i,k,j) = rremi(i,k,j) + (expo+expo_B)*sgb snpzd(ibion) = snpzd(ibion) + redptn*(expo+expo_B & )*sgb snpzd(ibiono3) = snpzd(ibiono3) + (expo+expo_B & )*sgb !----------------------------------------------------------------------- ! benthic denitrification model of Bohlen et al., 2012, GBC ! NO3 is removed out of bottom water nitrate. ! See Somes et al., 2012, BGS for additional details/results !----------------------------------------------------------------------- ! limit denitrification as nitrate approaches 0 uM no3flag = 0.5+sign(0.5,t(i,k,j,ino3,taum1)-trcmin) din15flag = 1.0 lno3 = 0.5*tanh(t(i,k,j,ino3,taum1)*10 - 5.0) sg_bdeni = (0.06 + 0.19*0.99 & **(max(t(i,k,j,io2,taum1)*1000,trcmin) & - max(t(i,k,j,ino3,taum1),trcmin))) & *max((expo & + expo_B & )*sgb,trcmin)*6.625 sg_bdeni = min(sg_bdeni, sgb*( & expo_B + & expo)) sg_bdeni = max(sg_bdeni, 0.) sg_bdeni = sg_bdeni*(0.5 + lno3)*no3flag*din15flag snpzd(ibiono3) = snpzd(ibiono3) - sg_bdeni src(i,k,j,isdic) = src(i,k,j,isdic)+expocaco3*1.e-3*sgb src(i,k,j,isalk) = src(i,k,j,isalk) & + 2.*expocaco3*1.e-3*sgb !---------------------------- ! sediment transfer function for iron !---------------------------- fesed=fetopsed*bct*redptn*(expo+expo_B & )*sgb bfe=fesed snpzd(ibiodfe) = snpzd(ibiodfe) + fesed if (t(i,k,j,io2,taum1)*1e3 .lt. o2min) then snpzd(ibiodfe) = snpzd(ibiodfe) + expofe*sgb bfe=bfe+expofe*sgb expofe = expofe - sgb * expofe endif rremife(i,k,j) = rremife(i,k,j) + bfe expo = expo - sgb * expo expo_B = expo_B - sgb * expo_B !---------------------------------------------------------------------- ! calculate calcite at the bottom and dissolve !--------------------------------------------------------------------- expocaco3 = expocaco3 - sgb*expocaco3 endif !end of bathy condition rexpo(i,k,j) = expo rexpo_B(i,k,j) = expo_B rexpocaco3(i,k,j) = expocaco3 rexpofe(i,k,j) = expofe if (k .eq. kmt(i,jrow)) then rremi(i,k,j) = rremi(i,k,j) + expo + expo_B snpzd(ibion) = snpzd(ibion) + redptn*(expo+expo_B) !nitrogen snpzd(ibiono3) = snpzd(ibiono3) + (expo+expo_B) !----------------------------------------------------------------------- ! benthic denitrification model of Bohlen et al., 2012, GBC ! NO3 is removed out of bottom water nitrate. ! See Somes et al., 2012, BGS for additional details/results !----------------------------------------------------------------------- ! limit denitrification as nitrate approaches 0 uM no3flag = 0.5+sign(0.5,t(i,k,j,ino3,taum1)-trcmin) din15flag = 1.0 lno3 = 0.5*tanh(t(i,k,j,ino3,taum1)*10 - 5.0) bdeni = (0.06 + 0.19*0.99 & **(max(t(i,k,j,io2,taum1)*1000,trcmin) & - max(t(i,k,j,ino3,taum1),trcmin))) & *max((expo & + expo_B & ),trcmin)*6.625 bdeni = min(bdeni,( & expo_B + & expo)) bdeni = max(bdeni, 0.) bdeni = bdeni*(0.5 + lno3)*no3flag*din15flag snpzd(ibiono3) = snpzd(ibiono3) - bdeni fesed=fetopsed*bct*redptn*(expo+expo_B & ) bfe=fesed snpzd(ibiodfe) = snpzd(ibiodfe) + fesed if (t(i,k,j,io2,taum1)*1e3 .lt. o2min) then snpzd(ibiodfe) = snpzd(ibiodfe) + expofe bfe=bfe+expofe endif rremife(i,k,j) = rremife(i,k,j) + bfe !---------------------------------------------------------------------- ! calculate calcite at the bottom and dissolve !--------------------------------------------------------------------- endif !end of bottom box loop !----------------------------------------------------------------------- ! set source/sink terms !----------------------------------------------------------------------- src(i,k,j,ispo4) = snpzd(ibion) src(i,k,j,isphyt) = snpzd(ibiop) src(i,k,j,iszoop) = snpzd(ibioz) src(i,k,j,isdetr) = snpzd(ibiod) src(i,k,j,isno3) = snpzd(ibiono3) src(i,k,j,isdiaz) = snpzd(ibiodiaz) src(i,k,j,iscocc) = snpzd(ibioc) src(i,k,j,isdiat) = snpzd(ibiodiat) src(i,k,j,iscaco3) = snpzd(ibiocaco3) src(i,k,j,isdetr_B) = snpzd(ibiod_B) src(i,k,j,isdfe) = snpzd(ibiodfe) + fe_hydr(i,j,k) src(i,k,j,isdetrfe) = snpzd(ibiodetrfe) !----------------------------------------------------------------------- ! production, export, and dissolution of opal !----------------------------------------------------------------------- ! sipr = 0. !# if defined 1 && (defined 1 || defined O_fe_limitation) ! From HAMOCC5 (Aumont et al, GBC, 2003) values 0-0.5, mostly ~0.1 ! sipr = sipr0 * min(1.,sil_in/k1si) !# if defined O_fe_limitation || defined 1 ! & * (4.-3.*min(1., !# if !defined 1 ! & fe_conc !# else ! & max(t(i,k,j,idfe,taum1),trcmin) !# endif ! &/kfe_Diat)) !# endif !# else ! sipr = sipr0 !# endif ! production of opal. dprop = (morp_Diat + graz_Diat*(1.-gamma1)) & *sipr*redctn*rnbio(k) !umol Si cm-3 s-1 prop = prop + dprop*dzt(k) !umol Si cm-2 s-1 src(i,k,j,issil) = -dprop + si_hydr(i,j,k)*si_h_sol !umol Si cm-3 s-1 silwflx = silwflx - si_hydr(i,j,k)*si_h_sol & *dzt(k)*dxt(i)*dyt(j)*cst(j) !umol Si/s !----------------------------------------------------------------------- ! production of calcite !----------------------------------------------------------------------- ! Still calculated but now only used if kk_caco3tr is not implemented! dprca = (capr*morp_C+capr*morz & +(capr*graz_C+capr*graz_Z)*(1.-gamma1)) & *redctn*rnbio(k) prca = prca + dprca*dzt(k) !mol C/m3/s ! These are sources and sinks of DIC (i.e. remin - pp) ! all are based on po4 uptake and remineralization ! dprca is a correction term ! calculated below ! calculated below ! print*,'KK3 tr', src(i,k,j,:) !----------------------------------------------------------------------- ! 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_C(i,k,jrow) = ta_rnpp_C(i,k,jrow) & + rnpp_C(i,k,j) ta_rgraz_C(i,k,jrow) = ta_rgraz_C(i,k,jrow) & + rgraz_C(i,k,j) ta_rmorp_C(i,k,jrow) = ta_rmorp_C(i,k,jrow) & + rmorp_C(i,k,j) ta_rmorpt_C(i,k,jrow)= ta_rmorpt_C(i,k,jrow) & + rmorpt_C(i,k,j) ta_rnpp_Diat(i,k,jrow) = ta_rnpp_Diat(i,k,jrow) & + rnpp_Diat(i,k,j) ta_rgraz_Diat(i,k,jrow) = ta_rgraz_Diat(i,k,jrow) & + rgraz_Diat(i,k,j) ta_rmorp_Diat(i,k,jrow) = ta_rmorp_Diat(i,k,jrow) & + rmorp_Diat(i,k,j) ta_rmorpt_Diat(i,k,jrow)= ta_rmorpt_Diat(i,k,jrow) & + rmorpt_Diat(i,k,j) ta_romca(i,k,jrow) = ta_romca(i,k,jrow) + romca(i,k,j) ta_rco3(i,k,jrow) = ta_rco3(i,k,jrow) + rco3(i,k,j) ta_rco3_sat(i,k,jrow) = ta_rco3_sat(i,k,jrow) & + rco3_sat(i,k,j) ta_rdel_sat(i,k,jrow) = ta_rdel_sat(i,k,jrow) & + rdel_sat(i,k,j) ta_rcalatt(i,k,jrow) = ta_rcalatt(i,k,jrow) & + rcalatt(i,k,j) ta_rprocal(i,k,jrow) = ta_rprocal(i,k,jrow) & + rcalpro(i,k,j) ta_rgraz_Det_B(i,k,jrow) = ta_rgraz_Det_B(i,k,jrow) & + rgraz_Det_B(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_ravej(i,k,jrow) = ta_ravej(i,k,jrow) & + ravej(i,k,j) ta_ravej_D(i,k,jrow) = ta_ravej_D(i,k,jrow) & + ravej_D(i,k,j) ta_rgmax(i,k,jrow) = ta_rgmax(i,k,jrow) & + rgmax(i,k,j) ta_rno3P(i,k,jrow) = ta_rno3P(i,k,jrow) & + rno3P(i,k,j) ta_rpo4P(i,k,jrow) = ta_rpo4P(i,k,jrow) & + rpo4P(i,k,j) ta_rpo4_D(i,k,jrow) = ta_rpo4_D(i,k,jrow) & + rpo4_D(i,k,j) endif ! calculate total export to get total import for next layer expo = expo*dzt(k) expo_B = expo_B*dzt(k) !ballast detritus expocaco3 = expocaco3*dzt(k) expofe = expofe*dzt(k) enddo kmx = kmt(i,jrow) do k=1,kmx ! 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 ! 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) deni = max(0., 800.*no3flag*so2*(0.5 - fo2)) src(i,k,j,isno3) = src(i,k,j,isno3) - deni ! Correct the ALK stoichiometry to account for N2 fixation ! New stoichiometric model parameters formulated as in Paulmier et al. 2009 BG 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 ! 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 + deni*1.e-3 rdeni(i,k,jrow) = deni & + sg_bdeni + bdeni enddo !----------------------------------------------------------------------- ! remineralize implicit calcite and opal !----------------------------------------------------------------------- kmx = kmt(i,jrow) do k=1,kmx-1 src(i,k,j,isdic) = src(i,k,j,isdic) & + src(i,k,j,ispo4)*redctp !mmol/m3/s & + rdissl(i,k,j)*1.e-3 & - rcalpro(i,k,j)*1.e-3 & - rcalatt(i,k,j)*1.e-3 old_diss = prca*rcak(k) src(i,k,j,isalk) = src(i,k,j,isalk) & - src(i,k,j,ispo4)*redntp*1.e-3 & + 2.*rdissl(i,k,j)*1.e-3 & - 2.*rcalpro(i,k,j)*1.e-3 & - 2.*rcalatt(i,k,j)*1.e-3 enddo src(i,kmx,j,isdic) = src(i,kmx,j,isdic) & + src(i,kmx,j,ispo4)*redctp & + rdissl(i,kmx,j)*1.e-3 & - rcalpro(i,kmx,j)*1.e-3 & - rcalatt(i,kmx,j)*1.e-3 & + rexpocaco3(i,kmx,j)*1.e-3 old_diss = prca*rcab(kmx) src(i,kmx,j,isalk) = src(i,kmx,j,isalk) & - src(i,kmx,j,ispo4)*redntp*1.e-3 & + 2.*rdissl(i,kmx,j)*1.e-3 & - 2.*rcalpro(i,kmx,j)*1.e-3 & - 2.*rcalatt(i,kmx,j)*1.e-3 & + 2.*rexpocaco3(i,kmx,j)*1.e-3 ! Test some dissolution parameters, a default exponential from Aumont et al. (2003) ! can be found in setmom.F expo_op = prop !flux kmx = kmt(i,jrow) do k=1,kmx temp_2(k) = t(i,k,j,itemp,taum1)+273.15 !convert to deg K ! (1) From Gnanadesikan 1999 !lower rates than HAMOCC ! might over-estimate dissl in deep ocean (discussed in Ridgwell et al 2002) ! ropk(k) = dzt(k)*si_dis/ws0*exp(-11481./temp_2(k)) !unitless theoretical dissl ! (2) From PlankTOM10, includes oxygen. Similar to HAMOCC except it always takes the minimum=1 and you end up with high dissolution.Plus there is no evidence of an oxygen dependency! ! ropk(k) = min(1.,si_dis*exp(-11200./temp_2(k))) ! & *avail_oxy(k)/ws0*dzt(k) ! (3) Let's include the microbial loop instead. I made this up. ropk(k) = dzt(k)*si_dis/ws0*bctd(i,k,j) !exp(-11481./temp_2(k)) ! if(rmorpt_Diat(i,k,j).gt.0)then ! ropk(k) = ropk(k)+ropk(k)*bctd(i,k,j)!*0.5!0.1 ! ropk(k)=ropk(k)*bctd(i,k,j) ! endif if(ropk(k).gt.1) ropk(k) = 1. op_diss = expo_op*ropk(k) !umol Si cm-2 s-1 expo_op = expo_op - op_diss src(i,k,j,issil) = src(i,k,j,issil) + op_diss*dztr(k) !umol Si cm-3 s-1 sgb = sg_bathy(i,j,k) if (sgb .gt. 0) then ! Use a limit (2 mmol Si /m^2/day) and ratios from Sarmiento & Gruber (2006) ! Buried Si is lost from the system!!! if (expo_op .le. 0.2/24/60/60) then bsi = 0.05*expo_op*sgb !buried sil umol/cm^2/s else bsi = 0.3*expo_op*sgb endif expo_op = expo_op - bsi silwflx = silwflx + bsi*dxt(i)*dyt(j)*cst(j) if (timavgperts .and. .not. euler2) then ta_rbur_sil(i,k,jrow) = ta_rbur_sil(i,k,jrow) + bsi*dztr(k) !umol/cm3/s endif endif if (timavgperts .and. .not. euler2) then ta_rdisopl(i,k,jrow) = ta_rdisopl(i,k,jrow) + op_diss*dztr(k) if (k .ne. kmt(i,jrow)) then ta_rexpoopl(i,k,jrow) = ta_rexpoopl(i,k,jrow) +expo_op*dztr(k) !umol/cm3/s else ta_rdisopl(i,k,jrow) = ta_rdisopl(i,k,jrow)+ expo_op*dztr(kmx) endif endif if (k .eq. kmt(i,jrow)) then src(i,kmx,j,issil) = src(i,kmx,j,issil) + expo_op*dztr(kmx) ! put the leftovers back into the ocean ! if (timavgperts .and. .not. euler2) then ! ta_rdisopl(i,k,jrow) = ta_rdisopl(i,k,jrow)+ expo_op*dztr(kmx) ! endif endif enddo ! K loop ! Here I should add the age change of the ideal age tracer kmx = kmt(i,jrow) do k=1,kmx src(i,k,j,isidealage) = 1.0 / 3.6e3 /24. enddo !----------------------------------------------------------------------- ! accumulate time averages for full depth variables !----------------------------------------------------------------------- if (timavgperts .and. .not. euler2) then kmx = kmt(i,jrow) expo = prca ta_rproopl(i,jrow) = ta_rproopl(i,jrow) + prop do k=1,kmx expo = expo*dztr(k) old_diss = old_diss*dztr(k) ta_rold_diss(i,k,jrow) = ta_rold_diss(i,k,jrow) & + old_diss old_diss = old_diss*dzt(k) ta_rremi_B(i,k,jrow) = ta_rremi_B(i,k,jrow) & + rremi_B(i,k,j) ta_rexpo_B(i,k,jrow) = ta_rexpo_B(i,k,jrow) & + rexpo_B(i,k,j) !ballast detritus 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) expo = expo - prca*rcak(k) ta_rexpocal(i,k,jrow) = ta_rexpocal(i,k,jrow) + expo expo = expo*dzt(k) ta_rdissl(i,k,jrow) = ta_rdissl(i,k,jrow) & + rdissl(i,k,j) ta_rexpocaco3(i,k,jrow) = & ta_rexpocaco3(i,k,jrow) + rexpocaco3(i,k,j) ! ta_rexpoopl(i,k,jrow) = ta_rexpoopl(i,k,jrow) +expo_op*dztr(k) !umol/cm3/s ! ta_rdisopl(i,k,jrow) = ta_rdisopl(i,k,jrow) + op_diss*dztr(k) ! ta_rbur_sil(i,k,jrow) = ta_rbur_sil(i,k,jrow) + bsi*dztr(k) !umol/cm3/s 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 endif enddo enddo rivsil = (silwflx-globalsilwflx)/globaldisch !scale Si over total discharge do j=js,je jrow = j + joff do i=is,ie if (kmt(i,jrow) .gt. 0) then src(i,1,j,issil) = src(i,1,j,issil) & + rivsil*disch(i,j)*dztr(1) if (timavgperts .and. .not. euler2) then ta_rrivsil(i,jrow) = ta_rrivsil(i,jrow) & + rivsil*disch(i,j) endif 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 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 !----------------------------------------------------------------------- ! 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) 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 !----------------------------------------------------------------------- ! 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) to 0.0 in the surface layer ! ---------------------------------------------------------------------- do j=js,je jrow = j + joff do i=istrt,iend t(i,1,j,iidealage,taup1)= 0.0 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, 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, isssil, isil) 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 delta carbon 14 !----------------------------------------------------------------------- if (tsiperts .or. timavgperts) 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 !----------------------------------------------------------------------- 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 "npzd.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 !