subroutine mobii #if defined O_mom && defined O_npzd ! initialize MOBI parameters implicit none include "size.h" include "npzd.h" include "calendar.h" include "coord.h" include "stdunits.h" include "scalar.h" integer k, ioun real alpha, par # if defined O_npzd_fe_limitation include "param.h" include "pconst.h" integer ib(10), ic(10), iou, m, j, i character (120) :: fname, new_file_name logical exists real c1e9 real tmpijkm(1:100,1:100,1:3,1:12) # endif namelist /npzd/ alpha, kw, kc, abio, bbio, cbio, k1n, nup &, gamma1, gbio, epsbio, nuz, nud0, LFe &, wd0, par, dtnpzd, redctn, ki, redptn, capr &, dcaco3, nupt0, redotn, jdiar, kzoo, geZ &, zprefP, zprefZ, zprefDet, zprefD, kfe, kfe_D &, nupt0_D, sgbdfac, diazntp, nup_D, dfr, dfrt &, nudon0, nudop0, eps_assim, eps_excr, eps_nfix &, eps_wcdeni, eps_bdeni0, eps_recy, hdop &, mw, mwz ! set defaults for namelist npzd alpha = 0.16 ! Initial slope P-I curve [(W/m^2)^(-1)/day] kw = 0.04 ! Light attenuation due to water [1/m] kc = 0.047 ! Light atten. by phytoplankton [1/(m*mmol/m^3)] ki = 5.0 ! Light attenuation through sea ice & snow abio = 0.6 ! a; Maximum growth rate parameter [1/day] bbio = 1.066 ! b cbio = 1.0 ! c [1/deg_C] k1n = 0.7 ! Half saturation constant for N uptake [mmol/m^3] nup = 0.03 ! Specific mortality rate (Phytoplankton) [1/day] nup_D = 0.0001 ! Specific mortaility rate (Diazotrophs [1/day] nupt0 = 0.015 ! Fast-recycling mortality rate (Phytoplankton) [1/day] nupt0_D = 0.001 ! Fast-recyling mortality rate (Diazotrophs) [1/day] gamma1 = 0.70 ! gama1; Assimilation efficiency (zpk) gbio = 0.38 ! Maximum grazing rate [1/day] epsbio = 1.6 ! Prey capture rate [(mmol/m^3)^(-2)day^(-1)] nuz = 0.06 ! Quadratic mortality (zpk) [(mmol/m^3)^(-2)day^(-1)] nud0 = 0.07 ! detritus remineralization rate [1/day] nudon0 = 2.33e-5 ! DON remineralization rate [1/day] nudop0 = 7.e-5 ! DOP remineralization rate [1/day] LFe = 1.0 ! Iron limitation wd0 = 16.0 ! Sinking speed of detritus at surface [m/day] mwz = 100000. !Depth where sinking below remains constant (cm) mw = 0.045 !Sinking speed increase with depth (s-1) par = 0.43 ! fraction of photosythetically active radiation dtnpzd = dtts/4. ! time step of biology [s] redctn = 7.1 ! C/N Redfield ratio (Schneider et al., 2003, GBC) redptn = 1./16. ! P/N Redfield ratio capr = 0.022 ! carbonate to carbon production ratio dcaco3 = 650000.0 ! remineralisation depth of calcite [cm] redotn = 10.6 ! O2/N Redfield ratio (Anderson & Sarmiento, 1994, GBC) jdiar = 0.08 ! factor reducing the growth rate of diazotrophs kzoo = 0.15 ! half sat. constant for Z grazing [mmol N m^3] geZ = 0.6 ! Zooplankton growth efficiency sgbdfac = 1.0 ! sub-grid benthic denitrification rate factor diazntp = 28. !diazotroph N:P ratio dfr = 0.08 ! phyt mortality refractory/semi-labile DOM fraction dfrt = 0.01 ! phyt fast-recy refracotyr/semi-labile DOM fraction hdop = 0.4 ! DOP growth rate handicap # if defined O_npzd_nitrogen zprefP = 0.30 ! Zooplankton preference for P zprefZ = 0.30 ! Zooplankton preference for other Z zprefDet = 0.30 ! Zooplankton preference for Detritus zprefD = 0.10 ! Zooplankton preference for diazotrophs # else zprefP = 0.35 ! Zooplankton preference for P zprefZ = 0.35 ! Zooplankton preference for other Z zprefDet = 0.30 ! Zooplankton preference for Detritus # endif # if defiend O_npzd_nitrogen_15 eps_assim = 6. eps_excr = 4. eps_nfix = 1. eps_wcdeni = 25. eps_bdeni0 = 6. eps_recy = 1. # endif # if defined O_npzd_fe_limitation kfe = 0.12 ! Half saturation constant for Fe limitation kfe_D = 0.12 ! Half saturation constant for Diaz Fe limitation # endif ! read namelist call getunit (ioun, 'control.in', 'f s r') read (ioun, npzd, end=108) 108 continue write (stdout, npzd) call relunit (ioun) ! convert units of NPZD parameters to MOM units redctn = redctn*1.e-3 redotn = redotn*1.e-3 redotp = redotn/redptn redctp = redctn/redptn redotc = redotn/redctn redntp = 1./redptn redntc = 1./redctn diazptn = 1./diazntp k1p = k1n*redptn kw = kw*1.e-2 kc = kc*1.e-2 ki = ki*1.e-2 wd0 = wd0*1.e2 alpha = alpha/daylen abio = abio/daylen nup = nup/daylen nup_D = nup_D/daylen nupt0 = nupt0/daylen nupt0_D = nupt0_D/daylen gbio = gbio/daylen epsbio = epsbio/daylen nuz = nuz/daylen nud0 = nud0/daylen nudop0 = nudop0/daylen nudon0 = nudon0/daylen tap = 2.*alpha*par ! calculate sinking speed of detritus divided by grid width do k=1,km ! linear increase wd0-200m with depth if (zt(k) .lt. mwz) then wd(k) = (wd0+mw*zt(k))/daylen/dzt(k) ! [s-1] else wd(k) = (wd0+mw*mwz)/daylen/dzt(k) ! [s-1] endif rkwz(k) = 1./(kw*dzt(k)) enddo ztt(1)=0.0 do k=1,km ztt(k+1)=(-1)*zw(k) enddo # if defined O_npzd_fe_limitation ! read in the dissolved Fe conc from O_dissolved_fe.nc fe_dissolved(:,:,:,:) = 0. ib(:) = 1 ic(:) = imtm2 ic(2) = jmtm2 ic(3) = 3 ic(4) = 12 fname = new_file_name ("O_fe_dissolved.nc") inquire (file=trim(fname), exist=exists) if (exists) then c1e9 = 1000000000 call openfile (trim(fname), iou) call getvara ('O_DISSOLVED_FE_LGM',iou, ic(1)*ic(2)*ic(3)*ic(4) &, ib, ic, tmpijkm, c1e9, c0) fe_dissolved(2:imtm1,2:jmtm1,:,:) = tmpijkm(1:imtm2 &, 1:jmtm2,:,:) ! print*,'fe_dissolved(50,50,1,1)=',fe_dissolved(50,50,1,1) do m=1,12 do k=1,3 do j=1,jmt fe_dissolved(1,j,k,m) = fe_dissolved(imtm1,j,k,m) fe_dissolved(imt,j,k,m) = fe_dissolved(2,j,k,m) enddo do i=1,imt fe_dissolved(i,1,k,m) = fe_dissolved(i,2,k,m) fe_dissolved(i,jmt,k,m) = fe_dissolved(2,j,k,m) enddo enddo enddo else print*,"Warning => Cannot find", trim(fname) endif # endif # if defined O_carbon || defined O_npzd_alk !--------------------------------------------------------------------- ! calculate variables used in calcite remineralization !--------------------------------------------------------------------- rcak(1) = -(exp(-zw(1)/dcaco3)-1.0)/dzt(1) rcab(1) = -1./dzt(1) do k=2,km rcak(k) = -(exp(-zw(k)/dcaco3))/dzt(k) & + (exp(-zw(k-1)/dcaco3))/dzt(k) rcab(k) = (exp(-zw(k-1)/dcaco3))/dzt(k) enddo # endif return end ! END mobii subroutine mobi (kmx, twodt, rctheta, dayfrac, swr, tnpzd &, t_in, po4_in # if defined O_npzd_o2 &, o2_in # endif # if defined O_carbon_13 &, s_in, dic_in, alk_in, co2_in, dic13_in # endif # if defined O_npzd_nitrogen &, no3_in, sgb_in, sgomsk_in # if defined O_npzd_nitrogen_15 &, din15_in # endif # endif # if defined O_npzd_fe_limitation &, felimit_in, felimit_D_in # endif # if defined O_sed &, sedcorgflx, sedcalflx) # endif &, src) ! main driver for Model of Ocean Biogeochemistry and Isotopes (MOBI) ! input: kmx, twodt, rctheta, dayfrac, swr, tnpzd, t_in, ! po4_in, (o2_in, s_in, dic_in, alk_in, co2_in, dic13_in, no3_in, ! sgb_in, sgomsk_in, din15_in, felimit_in, felimit_D_in) ! output: src, (sedcorgflx, sedcalflx) ! if O_save_npzd is switch on additional output (rnpp etc.) ! is in npzd.h ! copied from tracer.F (Andreas Schmittner, Oct 2, 2013) ! modified to a vertical column model implicit none integer k,kmx include "size.h" ! km maximum # of levels include "coord.h" ! zt, dzt include "grdvar.h" ! dztr = 1/dzt include "mw.h" ! ispo4 etc. include "npzd.h" real snpzd(ntnpzd), tnpzd(km,ntnpzd), src(km,nsrc) real dic_npzd_sms(km), twodt, rctheta, gl, impo, expo real npp, remi, excr, graz, morp, morpt, morz, temp, swr, dayfrac real graz_Det, graz_Z, phin, dz, prca, dprca, nud, bct, bctz real t_in(km), po4_in(km), sedrr, nudop, nudon # if defined O_npzd_extra_diagnostics real avej, avej_D, gmax, no3P, po4P, po4_D # endif # if defined O_npzd_o2 real o2_in(km), fo2, so2 # endif # if defined O_carbon_13 real s_in(km), dic_in(km), alk_in(km), co2_in, dic13_in(km) real rc13impo, rc13expo, ac13_DIC_aq real ac13_aq_POC, ac13b, prca13, rdic13, rtdic13 real pH, co2star, dco2star, pCO2, dpco2, CO3 real Omega_c, Omega_a, pCO2_old, pCO2_new real atmpres, depth # endif # if defined O_npzd_nitrogen real no3_in(km), sgb_in(km), sgomsk_in, lno3, sg_bdeni real morp_D, lntp real npp_D, graz_D, morpt_D, no3flag, wcdeni, bdeni(km), nfix(km) real din15_in(km), din15flag, eps_bdeni real rn15impo, rn15expo, uno3, rno3, bwcdeni, bbdeni # endif # if defined O_npzd_fe_limitation real felimit_D_in(km), felimit_in(km) # endif # if defined O_sed real sedcorgflx, sedcalflx # endif expo = 0.0 impo = 0.0 phin = 0.0 ! integrated phytoplankton prca = 0.0 ! integrated production of calcite sedrr = 0.0 # if defined O_save_npzd rsedrr = 0.0 rprca = 0.0 # endif # if defined O_npzd_nitrogen_15 rn15impo = 0.0 rn15expo = 0.0 # endif # if defined O_carbon_13 rc13impo = 0.0 rc13expo = 0.0 prca13 = 0.0 # endif !1111 start of main k-loop do k=1,kmx # if defined O_npzd_nitrogen_15 rn15impo = rn15expo # endif # if defined O_carbon_13 ! c13 biological fractionation atmpres = 1.0 !atm depth = zt(k)/100. !m call co2calc_SWS (t_in(k), s_in(k), dic_in(k), alk_in(k) &, co2_in, atmpres, depth, pH, co2star, dco2star &, pCO2, dpco2 &, CO3, Omega_c, Omega_a) ac13_DIC_aq = -1.0512994e-4*t_in(k)+1.011765 ! Popp et al. (1989) Am. J. Sci. ac13_aq_POC = -0.017*log10(co2star*1000.)+1.0034 ac13b = ac13_aq_POC/ac13_DIC_aq rc13impo = rc13expo # endif swr = swr*exp(-kc*phin) ! fix phin bug ! previous phyt layers already saved in swr term c phin = phin + max(tnpzd(k,3), trcmin)*dzt(k) phin = max(tnpzd(k,3), trcmin)*dzt(k) # if defined O_npzd_nitrogen & + max(tnpzd(k,9), trcmin)*dzt(k) # endif gl = swr*exp(ztt(k)*rctheta) impo = expo*dztr(k) bct = bbio**(cbio*t_in(k)) bctz = (0.5*(tanh(o2_in(k) - 8.)+1)) & *bbio**(cbio*min(t_in(k),20.)) # if defined O_npzd_o2 ! decrease remineralisation rate in oxygen minimum zone nud = nud0*(0.65+0.35*tanh(o2_in(k) - 6.0)) # else nud = nud0 # endif nudon = nudon0 nudop = nudop0 !----------------------------------------------------------------------- ! call the npzd ecosystem dynamics model !----------------------------------------------------------------------- call npzd_src (tnpzd(k,:), gl, bct, impo, dzt(k) &, dayfrac, wd(k), rkwz(k), nud, nudop, nudon &, snpzd, expo, graz, morp, morz, graz_Det, graz_Z # if defined O_save_npzd &, npp, morpt, remi, excr # if defined O_npzd_nitrogen &, npp_D, graz_D, morp_D, morpt_D &, nfix(k) # endif # if defined O_npzd_extra_diagnostics &, avej, avej_D, gmax, no3P &, po4P, po4_D # endif # endif # if defined O_npzd_fe_limitation &, felimit_in(k), felimit_D_in(k) # endif &, bctz # if defined O_npzd_nitrogen_15 &, rn15impo, rn15expo # endif # if defined O_carbon_13 &, rc13impo, rc13expo, ac13b # endif & ) ! These are source/sink terms snpzd = snpzd*rdtts expo = expo*rnbio # if defined O_npzd_nitrogen_15 rn15expo = rn15expo*rnbio # endif # if defined O_carbon_13 rc13expo = rc13expo*rnbio # endif # if defined O_save_npzd rexpo(k) = expo rgraz(k) = graz*rnbio rgraz_Det(k) = graz_Det*rnbio rgraz_Z(k) = graz_Z*rnbio rmorp(k) = morp*rnbio rmorz(k) = morz*rnbio rnpp(k) = npp*rnbio rmorpt(k) = morpt*rnbio rremi(k) = remi*rnbio rexcr(k) = excr*rnbio # if defined O_npzd_nitrogen rnpp_D(k) = npp_D*rnbio rgraz_D(k) = graz_D*rnbio rmorp_D(k) = morp_D*rnbio rmorpt_D(k) = morpt_D*rnbio rnfix(k) = nfix(k)*rnbio # endif # if defined O_npzd_extra_diagnostics ravej(k) = avej*rnbio ravej_D(k) = avej_D*rnbio rgmax(k) = gmax*rnbio rno3P(k) = no3P*rnbio rpo4P(k) = po4P*rnbio rpo4_D(k) = po4_D*rnbio # endif # endif !----------------------------------------------------------------------- ! calculate detritus at the bottom and remineralize !----------------------------------------------------------------------- if (k .eq. kmx) then sedrr = expo*dzt(k) # if defined O_sed sedcorgflx = sedcorgflx + expo*redctn*twodt*dzt(k) # endif # if defined O_save_npzd rremi(k) = rremi(k) + expo rsedrr = rsedrr + sedrr # endif snpzd(1) = snpzd(1) + redptn*expo snpzd(6) = snpzd(6) + redctn*expo # if defined O_carbon_13 snpzd(16) = snpzd(16) + rc13expo*redctn*expo # endif # if defined O_npzd_nitrogen !----------------------------------------------------------------------- ! benthic denitrification model of Bohlen et al., 2012, GBC (eq. 10) ! 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,no3_in(k)-trcmin) # if defined O_npzd_nitrogen_15 din15flag = 0.5+sign(0.5,din15_in(k)-trcmin) # else din15flag = 1. # endif lno3 = 0.5*tanh(no3_in(k)*10 - 5.0) bdeni(k) = (0.06 + 0.19*0.99**(max(o2_in(k),trcmin) & - max(no3_in(k),trcmin)))*max(expo,trcmin)*redctn & *1.e3 bdeni(k) = min(bdeni(k), expo) bdeni(k) = max(bdeni(k), 0.) bdeni(k) = bdeni(k)*(0.5 + lno3)*no3flag*din15flag & *sgomsk_in if (k .eq. 1) then bdeni(k) = 0.0 else if (k .eq. 2) then bdeni(k) = 0.08*bdeni(k) endif snpzd(7) = snpzd(7) + expo - bdeni(k) # if defined O_npzd_nitrogen_15 rno3 = max(din15_in(k),trcmin*rn15std/(1+rn15std)) & /max(no3_in(k)-din15_in(k),trcmin*rn15std/(1+rn15std)) rno3 = min(rno3, 2.*rn15std) rno3 = max(rno3, rn15std/2.) eps_bdeni = eps_bdeni0 !*exp(-2.5e-6*(zt(k)+dzt(k)/2.)) bbdeni = rno3 - eps_bdeni*rno3/1000. snpzd(10) = snpzd(10) + rn15expo*expo & - bbdeni/(1+bbdeni)*bdeni(k) # endif expo = 0. sg_bdeni = 0. else if (sgb_in(k) .gt. 0) then sedrr = sgb_in(k)*expo*dzt(k) no3flag = 0.5+sign(0.5,no3_in(k)-trcmin) # if defined O_npzd_nitrogen_15 din15flag = 0.5+sign(0.5,din15_in(k)-trcmin) # else din15flag = 1. # endif lno3 = 0.5*tanh(no3_in(k)*10 - 5.0) sg_bdeni = (0.06 + 0.19*0.99**(max(o2_in(k),trcmin) & - max(no3_in(k),trcmin))) & *max(expo*sgb_in(k),trcmin)*redctn*1.e3 sg_bdeni = min(sg_bdeni, sgb_in(k)*expo) sg_bdeni = max(sg_bdeni, 0.) c if (zt(k) .le. 50000. .and. sgb_in(k) .lt. 0.5) c & sg_bdeni = sgbdfac*sg_bdeni if (sgb_in(k) .gt. 0. .and. zt(k) .le. 24000) & sg_bdeni =min((1-sgb_in(k))*sgbdfac+1.,20.)*sg_bdeni sg_bdeni = sg_bdeni*(0.5 + lno3)*no3flag*din15flag c & *sgomsk_in c if (k .eq. 1) then c sg_bdeni = 0.0 c else if (k .eq. 2) then c sg_bdeni = 0.08*sg_bdeni c endif snpzd(1) = snpzd(1) + sgb_in(k)*expo*redptn snpzd(6) = snpzd(6) + sgb_in(k)*expo*redctn snpzd(7) = snpzd(7) + sgb_in(k)*expo - sg_bdeni # if defined O_carbon_13 snpzd(16) = snpzd(16) + rc13expo*sgb_in(k)*expo*redctn # endif # if defined O_npzd_nitrogen_15 rno3 = max(din15_in(k), trcmin*rn15std/(1+rn15std)) & /max(no3_in(k)-din15_in(k),trcmin*rn15std/(1+rn15std)) rno3 = min(rno3, 2.*rn15std) rno3 = max(rno3, rn15std/2.) eps_bdeni = eps_bdeni0 !*exp(-2.5e-6*(zt(k))) bbdeni = rno3 - eps_bdeni*rno3/1000. snpzd(10) = snpzd(10) + rn15expo*sgb_in(k)*expo & - bbdeni/(1+bbdeni)*sg_bdeni # endif expo = expo - sgb_in(k)*expo # if defined O_save_npzd rremi(k) = rremi(k) + sgb_in(k)*expo rsedrr = rsedrr + sedrr # endif else sg_bdeni = 0. endif bdeni(k) = sg_bdeni # endif endif # if defined O_npzd_nitrogen && O_save_npzd rbdeni(k) = bdeni(k) # endif !----------------------------------------------------------------------- ! set source/sink terms !----------------------------------------------------------------------- c if (k .ge. 3) then c src(k,ispo4) = snpzd(1) + 0.06e-8 !add to PO4 inventory c else src(k,ispo4) = snpzd(1) c endif src(k,isdop) = snpzd(2) src(k,isphyt) = snpzd(3) src(k,iszoop) = snpzd(4) src(k,isdetr) = snpzd(5) # if defined O_npzd_nitrogen src(k,isdic) = snpzd(6) src(k,isno3) = snpzd(7) src(k,isdon) = snpzd(8) src(k,isdiaz) = snpzd(9) # if defined O_npzd_nitrogen_15 src(k,isdin15) = snpzd(10) src(k,isdon15) = snpzd(11) src(k,isphytn15) = snpzd(12) src(k,iszoopn15) = snpzd(13) src(k,isdetrn15) = snpzd(14) src(k,isdiazn15) = snpzd(15) # endif # endif # if defined O_carbon_13 src(k,isdic13) = snpzd(16) src(k,isdoc13) = snpzd(17) src(k,isphytc13) = snpzd(18) src(k,iszoopc13) = snpzd(19) src(k,isdetrc13) = snpzd(20) # if defined O_npzd_nitrogen src(k,isdiazc13) = snpzd(21) # endif # endif ! production of calcite dprca = (morp+morz+(graz+graz_Z)*(1.-gamma1)) & *capr*redctn*rnbio prca = prca + dprca*dzt(k) ! These are sources and sinks of DIC (i.e. remin - pp) ! all are based on po4 uptake and remineralization ! dprca is a correction term # if defined O_carbon dic_npzd_sms(k) = snpzd(6) src(k,isdic) = (snpzd(6) - dprca) # if defined O_carbon_13 rtdic13 = max(dic13_in(k),trcmin*rc13std/(1+rc13std)) & / max(dic_in(k),trcmin) rtdic13 = min(rtdic13, 2.*rc13std/(1+rc13std)) rtdic13 = max(rtdic13, 0.5*rc13std/(1+rc13std)) prca13 = prca13 + dprca*dzt(k)*rtdic13 src(k,isdic13) = src(k,isdic13) - rtdic13*dprca # endif # endif # if defined O_npzd_alk src(k,isalk) = (-snpzd(6)*redntc*1.e-3 - 2.*dprca) # endif ! calculate total export to get total import for next layer expo = expo*dzt(k) enddo !1111 end of main k-loop # if defined O_save_npzd rprca = prca # endif # if defined O_npzd_o2 !2222 start of second k-loop do k=1,kmx ! limit oxygen consumption below concentrations of ! 5umol/kg as recommended in OCMIP fo2 = 0.5*tanh(o2_in(k) - 5.0) ! 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 = dic_npzd_sms(k)*redotc + nfix(k)*rnbio*1.25e-3 c so2 = dic_npzd_sms(k)*redotc src(k,iso2) = -so2*(0.5 + fo2) # if defined O_npzd_nitrogen ! add denitrification as source term for NO3 no3flag = 0.5+sign(0.5,no3_in(k)-trcmin) # if defined O_npzd_nitrogen_15 din15flag = 0.5+sign(0.5,din15_in(k)-trcmin) # else din15flag = 1. # endif lno3 = 0.5*tanh(no3_in(k) - 3.0) lntp = 0.5*tanh(no3_in(k)/(redntp*po4_in(k))*100 - 60.) ! 800 = 0.8*1000 = (elec/mol O2)/(elec/mol NO3)*(mmol/mol) wcdeni = 800.*no3flag*so2*(0.5 - fo2)*(0.5 + lno3)*din15flag wcdeni = max(wcdeni, 0.) src(k,isno3) = src(k,isno3) - wcdeni # if defined O_npzd_nitrogen_15 ! calculate isotope effect of water column denitrification uno3 = wcdeni*twodt/no3_in(k) uno3 = min(uno3, 0.999) uno3 = max(uno3, trcmin) rno3 = max(din15_in(k),trcmin*rn15std/(1+rn15std)) & / max(no3_in(k)-din15_in(k), & trcmin*rn15std/(1+rn15std)) rno3 = min(rno3, 2.*rn15std) rno3 = max(rno3, rn15std/2.) bwcdeni = rno3 + eps_wcdeni*(1-uno3)/uno3*log(1-uno3)*rno3/1000. src(k,isdin15) = src(k,isdin15) - (bwcdeni/(1+bwcdeni))*wcdeni # endif # if defined O_deni_alk ! Correct the ALK stoichiometry to account for N2 fixation src(k,isalk) = src(k,isalk) + wcdeni*1.e-3 !bdeni src(k,isalk) = src(k,isalk) + bdeni(k)*1.e-3 ! Now add account for N2 fixation (ALK production is tied to PO4 change and ! thus in the case of N2 fixation is not correct). src(k,isalk) = src(k,isalk) - nfix(k)*rnbio*1.e-3 # endif # if defined O_save_npzd rwcdeni(k) = wcdeni # endif # endif enddo !2222 end of second k-loop # endif !----------------------------------------------------------------------- ! remineralize calcite !----------------------------------------------------------------------- !3333 start of third k-loop do k=1,kmx-1 # if defined O_carbon src(k,isdic) = src(k,isdic) + prca*rcak(k) # if defined O_carbon_13 src(k,isdic13) = src(k,isdic13) + prca13*rcak(k) # endif # endif # if defined O_npzd_alk src(k,isalk) = src(k,isalk) + 2.*prca*rcak(k) # endif enddo !3333 end of third k-loop # if defined O_sed sedcalflx = sedcalflx + (prca*rcab(kmx) & - prca*rcak(kmx))*twodt*dzt(kmx) # endif # if defined O_carbon src(kmx,isdic) = src(kmx,isdic) + prca*rcab(kmx) # if defined O_carbon_13 src(kmx,isdic13) = src(kmx,isdic13) + prca13*rcab(kmx) # endif # endif # if defined O_npzd_alk src(kmx,isalk) = src(kmx,isalk) + 2.*prca*rcab(kmx) # endif return end ! END mobi subroutine npzd_src (bioin, gl, bct, impo, dzt &, dayfrac, wwd, rkw, nud, nudop, nudon &, bioout, expoout, grazout, morpout &, morzout, graz_Det_out, graz_Zout # if defined O_save_npzd &, nppout, morptout, remiout, excrout # if defined O_npzd_nitrogen &, npp_Dout, graz_Dout, morp_Dout, morpt_Dout &, nfixout # endif # if defined O_npzd_extra_diagnostics &, avej_out, avej_D_out, gmax_out, no3P_out &, po4P_out, po4_D_out # endif # endif # if defined O_npzd_fe_limitation &, felimit, felimit_D # endif &, bctz # if defined O_npzd_nitrogen_15 &, rn15impo, rn15expoout # endif # if defined O_carbon_13 &, rc13impo, rc13expoout, ac13b # endif & ) !======================================================================= ! computes source terms of the NPZD model ! initial version of code adapted from Xavier Giraud: ! Giraud et al. 2000, J Mar Res, 58, 609-630 ! original model reference: ! Oeschlies and Garcon 1999, Global Biogeochem. Cycles 13, 135-160 ! Schmittner et al. 2005, Global Biogeochem. Cycles 19, GB3004, ! doi:10.1029/2004GB002283. ! Schmittner et al. 2008, Global Biogeochem. Cycles 22, GB1013 ! ! This version was modified by David Keller and corrects the zooplankton ! grazing formulation. Note that zooplankton are now allowed to graze ! on themselves and detritus, in addition to phyt. and diazotrophs. ! The calculation of light has also been corrected. ! Keller et al. 2012, Geosci. Model Dev., 5, 1195-1220 ! ! The nitrogen isotope model (nitrogen_15) has been developed by ! Chris Somes and is described in ! Somes et al. 2010, Global Biogeochem. Cycles 24, GB4019 ! Somes et al. 2010, Geophys. Res. Lett. 37, L23605 and ! Somes et al. 2013, Biogeosc. 10, 5889-5910 ! ! The carbon isotope model (carbon_13) has been written by Andreas ! Schmittner and is documented in ! Schmittner et al. 2013 Biogeosc. 10, 5793-5816 ! Chris Somes modified the code converting from alpha to beta formulation ! ! Note that nutrient (N) represents phosphate ! ! input variables: ! ! bioin(1:5) = N,DON,P,Z,D [mmol m-3] ! bioin(6) = nitrate [mmol m-3] ! bioin(7) = DON [mmol m-3] ! bioin(8) = diazotrophs [mmol m-3] ! ! gl = 2.*light at top of grid box ! nbio = number of time steps ! dtbio = time step [s] ! bct = bbio**(cbio*temperature) ! impo = import of detritus from above [mmol m-3] ! dzt = depth of grid box [cm] ! dayfrac = day length (fraction: 0 < dayfrac < 1) ! wwd = sinking speed of detritus/dzt ! rkw = reciprical of kw*dzt(k) ! nud = remineralisation rate of detritus [s-1] ! ! output variables: ! ! bioout = change from bioin [mmol m-3] ! nppout = net primary production [mmol m-3] ! grazout = grazing [mmol m-3] ! morpout = quadratic mortality of phytoplankton [mmol m-3] ! morptout = specific mortality of phytoplankton [mmol m-3] ! morzout = mortality of zooplankton [mmol m-3] ! remiout = remineralisation [mmol m-3] ! excrout = excretion [mmol m-3] ! expoout = detrital export [mmol m-3] ! npp_Dout = NPP of diazotrophs ! graz_Dout = grazing of diazotrophs ! morp_Dout = mortality of diazotrophs ! nfixout = rate of N2 fixation ! graz_Det_out = grazing of detritus ! graz_Zout = grazing on othe zooplankton ! avej_out = light-depend phyt. growth rate ! avej_D_out = light-depend Diaz growth rate ! gmax_out = temp-depend. zoo growth rate ! no3P_out = no3 depend. phyt growth rate ! po4P_out = po4 depend. phyt growth rate ! po4_D_out = po4 depend. Diaz growth rate ! ! New grazing formulation variables and parameters ! ! The following terms determine ingestion according to a ! a Holling II curve (i.e. Michaelis Menten): ! ! Ingestion = max_graz_rate * (Ft/(Ft + kzoo)) ! ! where Ft is the weighted measure of the total food available ! and equals the sum of the different prey types times the ! preference of Z for that type of prey ! ! zprefP = Z preference for P ! zprefD = Z preference for Diaz ! zprefDet = Z preference for detritus ! zprefZ = Z preference for other Z ! kzoo = half saturation coefficienct for Z ingestion mmol N m-3 ! ing_P = zooplankton ingestion of phytoplankon ! ing_D = zooplankton ingestion of diazotrophs ! ing_Det = zooplankton ingestion of detritus ! ing_Z = zooplankton ingestion of other zooplankton ! thetaZ = Michaelis-Menten denominator ! ! felimit = Fe limitation parameter ! felmit_D = Fe limitation parameter for diazotrophs ! !======================================================================= implicit none integer n real gl, f1, biopo4, biophyt, biozoop, biodetr, jmax, u_P, g_P real morp, morpt, morz, remi, excr, expo, impo, nppout, grazout real morpout, morptout, morzout, remiout, excrout, expoout real avej_out, avej_D_out, gmax_out, no3P_out, po4P_out, po4_D_out real dzt, po4flag, phytflag, zoopflag, detrflag, wwd, rkw, gd real nupt, nud, biono3, u_D,npp_D, npp_Dout, no3flag, biodiaz, bct real diazflag, g_D,graz_D, morpt_D, jmax_D, gd_D, avej_D, no3upt_D real morpt_Dout, graz_Dout, nfixout, biop2, u1, u2, phi1, phi2 real avej, graz_Det_out, graz_Zout, thetaZ, ing_P, ing_D, dayfrac real ing_Det, ing_Z, g_Z, g_Det, graz_Z, graz_Det, gmax real no3P, po4P, po4_D, felimit, bctz, felimit_D real excrdiaz, biodop, biodon, dopflag, donflag, recy_don real recy_dop, npp, graz, biodin15, biodon15, biophytn15 real biodetrn15, biodiazn15, bassim, fcassim, bexcr, fcexcr real bnfix, fcnfix, rn15impo, rn15expoout, dig, dig_P, dig_Z real dig_Det, dig_D, excr_P, excr_Z, excr_Det, excr_D, sf, sf_P real sf_Z, sf_Det, sf_D, nr_excr_D, uno3, rno3, rzoop real rtdin15, rtphytn15, rtzoopn15, rtdetrn15, rtdiazn15, rtdon15 real din15flag, biozoopn15, morp_Dout real don15flag, phytn15flag, zoopn15flag, detrn15flag, diazn15flag real limP, limP_D, dopupt_D_flag, dopupt_D, morp_D, nupt_D real udon, rdon, brecy, fcrecy, biodic, dopupt_flag, dopupt real biodic13, biophytc13, biozoopc13, biodetrc13, biodiazc13 real dic13flag, phytc13flag, zoopc13flag, detrc13flag, diazc13flag real biodoc13, doc13flag, biodoc real rdic13, rtphytc13, rtzoopc13, rtdetrc13, rtdiazc13, ac13b real rc13impo, rc13expoout, dicflag, bc13npp, rtdic13, fcnpp real rtdoc13, nudop, nudon include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "calendar.h" include "npzd.h" real bioin(ntnpzd), bioout(ntnpzd) bioout(:) = 0.0 biopo4 = bioin(1) biodop = bioin(2) biophyt = bioin(3) biozoop = bioin(4) biodetr = bioin(5) # if defined O_npzd_nitrogen biodic = bioin(6) biono3 = bioin(7) biodon = bioin(8) biodiaz = bioin(9) # if defined O_npzd_nitrogen_15 biodin15 = bioin(10) biodon15 = bioin(11) biophytn15 = bioin(12) biozoopn15 = bioin(13) biodetrn15 = bioin(14) biodiazn15 = bioin(15) # endif # endif # if defined O_carbon_13 biodic13 = bioin(16) biodoc13 = bioin(17) biophytc13 = bioin(18) biozoopc13 = bioin(19) biodetrc13 = bioin(20) # if defined O_npzd_nitrogen biodiazc13 = bioin(21) # endif # endif ! flags prevent negative values by setting outgoing fluxes to ! zero if tracers are lower than trcmin po4flag = 0.5 + sign(0.5,biopo4 - trcmin) dopflag = 0.5 + sign(0.5,biodop - trcmin) phytflag = 0.5 + sign(0.5,biophyt - trcmin) zoopflag = 0.5 + sign(0.5,biozoop - trcmin) detrflag = 0.5 + sign(0.5,biodetr - trcmin) # if defined O_npzd_nitrogen no3flag = 0.5 + sign(0.5,biono3 - trcmin) donflag = 0.5 + sign(0.5,biodon - trcmin) diazflag = 0.5 + sign(0.5,biodiaz - trcmin) # if defined O_npzd_nitrogen_15 din15flag = 0.5 + sign(0.5, biodin15 - trcmin) don15flag = 0.5 + sign(0.5, biodon15 - trcmin) phytn15flag = 0.5 + sign(0.5, biophytn15 - trcmin) zoopn15flag = 0.5 + sign(0.5, biozoopn15 - trcmin) detrn15flag = 0.5 + sign(0.5, biodetrn15 - trcmin) diazn15flag = 0.5 + sign(0.5, biodiazn15 - trcmin) # else din15flag = 1. don15flag = 1. diazn15flag = 1. phytn15flag = 1. zoopn15flag = 1. detrn15flag = 1. # endif # endif # if defined O_carbon_13 dic13flag = 0.5 + sign(0.5,biodic13 - trcmin) doc13flag = 0.5 + sign(0.5,biodoc13 - trcmin) phytc13flag = 0.5 + sign(0.5,biophytc13 - trcmin) zoopc13flag = 0.5 + sign(0.5,biozoopc13 - trcmin) detrc13flag = 0.5 + sign(0.5,biodetrc13 - trcmin) diazc13flag = 0.5 + sign(0.5,biodiazc13 - trcmin) # endif ! limit tracers to positive values bioin(:) = max(bioin(:), trcmin) biopo4 = max(biopo4, trcmin) biodop = max(biodop, trcmin) biophyt = max(biophyt, trcmin) biozoop = max(biozoop, trcmin) biodetr = max(biodetr, trcmin) # if defined O_npzd_nitrogen biodic = max(biodic, trcmin) biono3 = max(biono3, trcmin) biodon = max(biodon, trcmin) biodiaz = max(biodiaz, trcmin) # if defined O_npzd_nitrogen_15 biodin15 = max(biodin15, trcmin) biodon15 = max(biodon15, trcmin) biophytn15 = max(biophytn15, trcmin) biozoopn15 = max(biozoopn15, trcmin) biodetrn15 = max(biodetrn15, trcmin) biodiazn15 = max(biodiazn15, trcmin) # endif # endif # if defined O_carbon_13 biodic13 = max(biodic13, trcmin) biodoc13 = max(biodoc13, trcmin) biophytc13 = max(biophytc13, trcmin) biozoopc13 = max(biozoopc13, trcmin) biodetrc13 = max(biodetrc13, trcmin) biodiazc13 = max(biodiazc13, trcmin) # endif ! photosynthesis after Evans & Parslow (1985) ! notation as in JGOFS report No. 23 p. 6 f1 = exp((-kw - kc*(biophyt+biodiaz))*dzt) # if defined O_npzd_fe_limitation ! In the following "felimit" is determined by an iron mask and ! is used to limit phytoplankton growth in HNLC regions jmax = abio*bct*felimit # else jmax = abio*bct # endif gd = jmax*dayfrac u1 = max(gl/gd,1.e-6) u2 = u1*f1 ! for the following approximation ensure that u1 < 20 phi1 = log(u1+sqrt(1.+u1**2.))-(sqrt(1.+u1**2.)-1.)/u1 phi2 = log(u2+sqrt(1.+u2**2.))-(sqrt(1.+u2**2.)-1.)/u2 avej = gd*(phi1 - phi2)/((kw+kc*(biophyt+biodiaz))*dzt) ! Make the max grazing rate a function of temperature gmax = gbio*bctz ! Note bctz, which sets an upper limit on the effects of temp on the ! grazing rate, is set in tracers.F # if defined O_npzd_nitrogen # if defined O_npzd_fe_limitation jmax_D = max(0.,abio*bct*felimit_D)*jdiar # else jmax_D = max(0.,abio*bct)*jdiar # endif gd_D = max(1.e-14,jmax_D*dayfrac) u1 = max(gl/gd_D,1.e-6) u2 = u1*f1 ! for the following approximation ensure that u1 < 20 phi1 = log(u1+sqrt(1.+u1**2.))-(sqrt(1.+u1**2.)-1.)/u1 phi2 = log(u2+sqrt(1.+u2**2.))-(sqrt(1.+u2**2.)-1.)/u2 avej_D = gd_D*(phi1 - phi2)/((kw+kc*(biophyt+biodiaz))*dzt) # endif # if defined O_npzd_nitrogen ! check grazing preferences = 1 for N case IF ((zprefP + zprefDet + zprefZ + zprefD).ne.1) THEN zprefP = 0.3 zprefZ = 0.3 zprefDet = 0.3 zprefD = 0.1 END IF # else ! check that grazing preferences = 1 for no N case IF ((zprefP + zprefDet + zprefZ).ne.1) THEN zprefP = 0.35 zprefZ = 0.35 zprefDet = 0.30 END IF # endif nupt = nupt0*bct nupt_D = nupt0_D*bct expoout = 0.0 grazout = 0.0 morpout = 0.0 morzout = 0.0 graz_Det_out = 0.0 graz_Zout = 0.0 # if defined O_npzd_nitrogen_15 rn15expoout = 0.0 # endif # if defined O_carbon_13 rc13expoout = 0.0 # endif # if defined O_save_npzd nppout = 0.0 morptout = 0.0 remiout = 0.0 excrout = 0.0 # if defined O_npzd_nitrogen npp_Dout = 0.0 graz_Dout = 0.0 morpt_Dout = 0.0 morp_Dout = 0.0 nfixout = 0.0 # endif # if defined O_npzd_extra_diagnostics avej_out = 0.0 avej_D_out = 0.0 gmax_out = 0.0 no3P_out = 0.0 po4P_out = 0.0 po4_D_out = 0.0 # endif # endif do n=1,nbio ! growth rate of phytoplankton ! consume DOP when it is more efficient if (hdop*biodop/(k1p + biodop) .gt. biopo4/(k1p + biopo4)) then dopupt_flag = 1. limP = biodop/(k1p + biodop) u_P = min(avej, jmax*limP*hdop) else dopupt_flag = 0. limP = biopo4/(k1p + biopo4) u_P = min(avej, jmax*limP) endif po4P = jmax*biopo4/(k1p + biopo4) # if defined O_npzd_nitrogen ! nitrate limitation u_P = min(u_P, jmax*biono3/(k1n + biono3)) no3P = jmax*biono3/(k1n + biono3) ! growth rate of diazotrophs smaller than other phytoplankton and ! not nitrate limited if (hdop*biodop/(k1p + biodop) .gt. biopo4/(k1p + biopo4)) then dopupt_D_flag = 1. limP_D = biodop/(k1p + biodop) u_D = min(avej_D, jmax_D*limP_D*hdop) else dopupt_D_flag = 0. limP_D = biopo4/(k1p + biopo4) u_D = min(avej_D, jmax_D*limP_D) endif po4_D = jmax_D*biopo4/(k1p + biopo4) ! Set the grazing coefficients for the N case thetaZ = zprefP*biophyt + zprefDet*biodetr + zprefZ*biozoop & + zprefD*biodiaz + kzoo ing_P = zprefP/thetaZ ing_Det = zprefDet/thetaZ ing_Z = zprefZ/thetaZ ing_D = zprefD/thetaZ # else ! If "else" then set the grazing coefficients for the no N case ! note kzoo is in terms of N so convert to P thetaZ = zprefP*biophyt + zprefDet*biodetr + zprefZ*biozoop & +kzoo*redptn ing_P = zprefP/thetaZ ing_Det = zprefDet/thetaZ ing_Z = zprefZ/thetaZ # endif npp = u_P*biophyt dopupt = npp*dopupt_flag # if defined O_npzd_nitrogen npp_D = max(0.,u_D*biodiaz) ! grazing on diazotrophs g_D = gmax*ing_D*biodiaz graz_D = g_D*biozoop morpt_D = nupt_D*biodiaz ! linear mortality morp_D = nup_D*biodiaz*biodiaz c no3upt_D = biono3/(k1n + biono3)*npp_D ! nitrate uptake no3upt_D = (0.5+0.5*tanh(biono3-5.))*npp_D ! nitrate uptake dopupt_D = npp_D*dopupt_D_flag # endif ! grazing on P g_P = gmax*ing_P*biophyt graz = g_P*biozoop ! grazing on Z g_Z = gmax*ing_Z*biozoop graz_Z = g_Z*biozoop ! grazing on Detritus g_Det = gmax*ing_Det*biodetr graz_Det = g_Det*biozoop ! morp = nup*biophyt morpt = nupt*biophyt recy_don = nudon*bct*biodon recy_dop = nudop*bct*biodop morz = nuz*biozoop*biozoop remi = nud*bct*biodetr expo = wwd*biodetr graz = graz*phytflag*zoopflag*phytn15flag*zoopn15flag graz_Z = graz_Z*zoopflag*zoopn15flag graz_Det = graz_Det*detrflag*zoopflag*detrn15flag*zoopn15flag morp = morp*phytflag*phytn15flag morpt = morpt*phytflag*phytn15flag morz = morz*zoopflag*zoopn15flag remi = remi*detrflag*detrn15flag expo = expo*detrflag*detrn15flag recy_dop = recy_dop*dopflag # if defined O_npzd_nitrogen if (dopupt_flag .eq. 1) then npp = npp*dopflag*no3flag*din15flag else npp = npp*po4flag*no3flag*din15flag endif if (dopupt_D .eq. 1) then npp_D = npp_D*dopflag else npp_D = npp_D*po4flag endif graz_D = graz_D*diazflag*zoopflag*diazn15flag*zoopn15flag morpt_D = morpt_D*diazflag*diazn15flag morp_D = morp_D*diazflag*diazn15flag no3upt_D = no3upt_D*no3flag*din15flag recy_don = recy_don*donflag*don15flag # else npp = npp*po4flag # endif # if defined O_npzd_nitrogen ! digestion of grazed material dig_P = gamma1*graz dig_Z = gamma1*graz_Z dig_Det = gamma1*graz_Det dig_D = gamma1*graz_D*(redntp/diazntp) dig = dig_P + dig_Z + dig_Det + dig_D ! excretion based on growth efficiency excr_P = gamma1*(1-geZ)*graz excr_Z = gamma1*(1-geZ)*graz_Z excr_Det = gamma1*(1-geZ)*graz_Det excr_D = gamma1*(1-geZ)*graz_D*(redntp/diazntp) excr = excr_P + excr_Z + excr_Det + excr_D ! excrete "extra" non-Redfield diazotroph N nr_excr_D = gamma1*graz_D*(1-(redntp/diazntp)) & + (1-gamma1)*graz_D*(1-(redntp/diazntp)) ! sloppy feeding sf_P = (1-gamma1)*graz sf_Z = (1-gamma1)*graz_Z sf_Det = (1-gamma1)*graz_Det sf_D = (1-gamma1)*graz_D*(redntp/diazntp) sf = sf_P + sf_Z + sf_Det + sf_D # else c excr = (gamma1-geZ)*(graz+graz_Z+graz_Det) # endif # if defined O_npzd_nitrogen_15 ! calculate isotope parameters ! See Somes et al., 2010, GBC for details/results uno3 = npp*dtbio/biono3 uno3 = min(uno3, 0.999) uno3 = max(uno3, trcmin) rno3 = biodin15/(biono3-biodin15) rno3 = min(rno3, 2*rn15std) rno3 = max(rno3, rn15std/2.) bassim = rno3 + eps_assim*(1-uno3)/uno3*log(1-uno3)*rno3/1000. fcassim = bassim/(1+bassim) udon = recy_don*dtbio/biodon udon = min(udon, 0.999) udon = max(udon, trcmin) rdon = biodon15/(biodon-biodon15) rdon = min(rdon, 2*rn15std) rdon = max(rdon, rn15std/2.) brecy = rdon + eps_recy*(1-udon)/udon*log(1-udon)*rdon/1000. fcrecy = brecy/(1+brecy) rzoop = biozoopn15/(biozoop-biozoopn15) rzoop = min(rzoop, 2.*rn15std) rzoop = max(rzoop, rn15std/2.) bexcr = rzoop - eps_excr*rzoop/1000. fcexcr = bexcr/(1+bexcr) bnfix = rn15std - eps_nfix*rn15std/1000. fcnfix = bnfix/(1+bnfix) rtdin15 = biodin15/biono3 rtdin15 = min(rtdin15, 2*rn15std/(1+rn15std)) rtdin15 = max(rtdin15, rn15std/(1+rn15std)/2.) rtdon15 = biodon15/biodon rtdon15 = min(rtdon15, 2*rn15std/(1+rn15std)) rtdon15 = max(rtdon15, rn15std/(1+rn15std)/2.) rtphytn15 = biophytn15/biophyt rtphytn15 = min(rtphytn15, 2.*rn15std/(1+rn15std)) rtphytn15 = max(rtphytn15, rn15std/(1+rn15std)/2.) rtzoopn15 = biozoopn15/biozoop rtzoopn15 = min(rtzoopn15, 2.*rn15std/(1+rn15std)) rtzoopn15 = max(rtzoopn15, rn15std/(1+rn15std)/2.) rtdetrn15 = biodetrn15/biodetr rtdetrn15 = min(rtdetrn15, 2.*rn15std/(1+rn15std)) rtdetrn15 = max(rtdetrn15, rn15std/(1+rn15std)/2.) rtdiazn15 = biodiazn15/biodiaz rtdiazn15 = min(rtdiazn15, 2.*rn15std/(1+rn15std)) rtdiazn15 = max(rtdiazn15, rn15std/(1+rn15std)/2.) # endif # if defined O_carbon_13 rdic13 = biodic13/(biodic-biodic13) rdic13 = min(rdic13, 2.*rc13std) rdic13 = max(rdic13, 0.5*rc13std) bc13npp = ac13b*rdic13 fcnpp = bc13npp/(1+bc13npp) rtdic13 = biodic13/biodic rtdic13 = min(rtdic13, 2*rc13std/(1+rc13std)) rtdic13 = max(rtdic13, 0.5*rc13std/(1+rc13std)) rtdoc13 = biodoc13/(biodon*redctn) rtdoc13 = min(rtdoc13, 2*rc13std/(1+rc13std)) rtdoc13 = max(rtdoc13, 0.5*rc13std/(1+rc13std)) rtphytc13 = biophytc13/(biophyt*redctn) rtphytc13 = min(rtphytc13, 2.*rc13std/(1+rc13std)) rtphytc13 = max(rtphytc13, 0.5*rc13std/(1+rc13std)) rtzoopc13 = biozoopc13/(biozoop*redctn) rtzoopc13 = min(rtzoopc13, 2.*rc13std/(1+rc13std)) rtzoopc13 = max(rtzoopc13, 0.5*rc13std/(1+rc13std)) rtdetrc13 = biodetrc13/(biodetr*redctn) rtdetrc13 = min(rtdetrc13, 2.*rc13std/(1+rc13std)) rtdetrc13 = max(rtdetrc13, 0.5*rc13std/(1+rc13std)) rtdiazc13 = biodiazc13/(biodiaz*redctn) rtdiazc13 = min(rtdiazc13, 2.*rc13std/(1+rc13std)) rtdiazc13 = max(rtdiazc13, 0.5*rc13std/(1+rc13std)) # endif # if defined O_npzd_nitrogen ! nutrients equation biopo4 = biopo4 + dtbio*(redptn*((1.-dfrt)*morpt + excr + remi & - (npp-dopupt)) + diazptn*(morpt_D - (npp_D-dopupt_D)) & + recy_dop) ! DOP equation biodop = biodop + dtbio*(redptn*dfr*morp + redptn*dfrt*morpt & - redptn*dopupt - diazptn*dopupt_D - recy_dop) ! phytoplankton equation biophyt = biophyt + dtbio*(npp - morp - graz - morpt) ! zooplankton equation biozoop = biozoop + dtbio*(dig - morz - graz_Z - excr) ! detritus equation biodetr = biodetr + dtbio*((1.-dfr)*morp + sf + morz - remi & - graz_Det - expo + impo + morp_D*(redntp/diazntp)) ! DIC equation biodic = biodic + dtbio*redctn*((1.-dfrt)*morpt + excr + remi & - npp + morpt_D - npp_D + recy_don + nr_excr_D & + morp_D*(1-(redntp/diazntp))) ! nitrate (NO3) equation biono3 = biono3 + dtbio*((1.-dfrt)*morpt + excr + morpt_D & + morp_D*(1-(redntp/diazntp)) + nr_excr_D + remi + recy_don & - npp - no3upt_D) ! DON equation biodon = biodon + dtbio*(dfr*morp + dfrt*morpt - recy_don) ! diazotroph equation biodiaz = biodiaz + dtbio*(npp_D - morp_D - morpt_D - graz_D) # else ! nutrients equation biopo4 = biopo4 + dtbio*redptn*(remi + excr - npp + morpt) ! phytoplankton equation biophyt = biophyt + dtbio*(npp - morp - graz - morpt) ! zooplankton equation biozoop = biozoop + dtbio*(geZ*(graz + graz_Det + graz_Z) - morz & - graz_Z) ! detritus equation biodetr = biodetr + dtbio*((1.-gamma1)*(graz + graz_Det & + graz_Z) + morp + morz - remi - graz_Det - expo + impo) # endif # if defined O_npzd_nitrogen_15 ! isotope equations biodin15 = biodin15 + dtbio*(rtphytn15*(1.-dfrt)*morpt & + fcexcr*excr + rtdiazn15*morpt_D + rtdiazn15*nr_excr_D & + rtdiazn15*morp_D*(1-(redntp/diazntp)) + rtdetrn15*remi & + fcrecy*recy_don - fcassim*npp - fcassim*no3upt_D) biodon15 = biodon15 + dtbio*(dfr*rtphytn15*morp & + dfrt*rtphytn15*morpt - fcrecy*recy_don) biophytn15 = biophytn15 + dtbio*(fcassim*npp - rtphytn15*morp & - rtphytn15*graz - rtphytn15*morpt) biozoopn15 = biozoopn15 + dtbio*(rtphytn15*dig_P & + rtzoopn15*dig_Z + rtdetrn15*dig_Det + rtdiazn15*dig_D & - rtzoopn15*morz - rtzoopn15*graz_Z - fcexcr*excr) biodetrn15 = biodetrn15 + dtbio*(rtphytn15*(1.-dfr)*morp & + rtphytn15*sf_P + rtzoopn15*sf_Z + rtdetrn15*sf_Det & + rtdiazn15*sf_D + rtzoopn15*morz - rtdetrn15*remi & - rtdetrn15*graz_Det - rtdetrn15*expo + rn15impo*impo & + rtdiazn15*morp_D*(redntp/diazntp)) biodiazn15 = biodiazn15 + dtbio*(fcnfix*(npp_D-no3upt_D) & + fcassim*no3upt_D - rtdiazn15*morp_D - rtdiazn15*graz_D & - rtdiazn15*morpt_D) # endif # if defined O_carbon_13 ! !!!!!!!!!!!!!!!!!! isotope equations biodic13 = biodic13 + dtbio*redctn*(rtphytc13*(1.-dfrt)*morpt & + rtzoopc13*excr + rtdiazc13*morpt_D + rtdiazc13*nr_excr_D & + rtdiazc13*morp_D*(1-(redntp/diazntp)) + rtdetrc13*remi & + rtdoc13*recy_don - fcnpp*npp - fcnpp*npp_D) biodoc13 = biodoc13 + dtbio*redctn*(dfr*rtphytc13*morp & + rtphytc13*dfrt*morpt - rtdoc13*recy_don) biophytc13 = biophytc13 + dtbio*redctn*(fcnpp*npp & - rtphytc13*morp - rtphytc13*graz - rtphytc13*morpt) biozoopc13 = biozoopc13 + dtbio*redctn*(rtphytc13*dig_P & + rtzoopc13*dig_Z + rtdetrc13*dig_Det + rtdiazc13*dig_D & - rtzoopc13*morz - rtzoopc13*graz_Z - rtzoopc13*excr) biodetrc13 = biodetrc13 + dtbio*redctn*(rtphytc13*(1.-dfr)*morp & + rtphytc13*sf_P + rtzoopc13*sf_Z + rtdetrc13*sf_Det & + rtdiazc13*sf_D + rtzoopc13*morz - rtdetrc13*remi & - rtdetrc13*graz_Det - rtdetrc13*expo + rc13impo*impo & + rtdiazc13*morp_D*(redntp/diazntp)) biodiazc13 = biodiazc13 + dtbio*redctn*(fcnpp*npp_D & - rtdiazc13*morp_D - rtdiazc13*graz_D - rtdiazc13*morpt_D) # endif expoout = expoout + expo # if defined O_npzd_nitrogen_15 rn15expoout = rn15expoout + rtdetrn15 # endif # if defined O_carbon_13 rc13expoout = rc13expoout + rtdetrc13 # endif grazout = grazout + graz morpout = morpout + morp morzout = morzout + morz graz_Det_out = graz_Det_out + graz_Det graz_Zout = graz_Zout + graz_Z # if defined O_save_npzd nppout = nppout + npp morptout = morptout + morpt remiout = remiout + remi excrout = excrout + excr # if defined O_npzd_nitrogen npp_Dout = npp_Dout + npp_D graz_Dout = graz_Dout + graz_D morpt_Dout = morpt_Dout + morpt_D morp_Dout = morp_Dout + morp_D nfixout = nfixout + npp_D - no3upt_D # endif # if defined O_npzd_extra_diagnostics avej_out = avej_out + avej avej_D_out = avej_D_out + avej_D gmax_out = gmax_out + gmax no3P_out = no3P_out + no3P po4P_out = po4P_out + po4P po4_D_out = po4_D_out + po4_D # endif # endif ! flags prevent negative values by setting outgoing fluxes to ! zero if tracers are lower than trcmin if (po4flag .eq. 1) po4flag = 0.5 + sign(0.5,biopo4 - trcmin) if (dopflag .eq. 1) dopflag = 0.5 + sign(0.5,biodop - trcmin) if (phytflag .eq. 1) phytflag = 0.5 + sign(0.5,biophyt - trcmin) if (zoopflag .eq. 1) zoopflag = 0.5 + sign(0.5,biozoop - trcmin) if (detrflag .eq. 1) detrflag = 0.5 + sign(0.5,biodetr - trcmin) # if defined O_npzd_nitrogen if (no3flag .eq. 1) no3flag = 0.5 + sign(0.5,biono3 - trcmin) if (donflag .eq. 1) donflag = 0.5 + sign(0.5,biodon - trcmin) if (diazflag .eq. 1) diazflag = 0.5 + sign(0.5,biodiaz - trcmin) # if defined O_npzd_nitrogen_15 if (din15flag .eq. 1) & din15flag = 0.5 + sign(0.5, biodin15 - trcmin) if (don15flag .eq. 1) & don15flag = 0.5 + sign(0.5, biodon15 - trcmin) if (phytn15flag .eq. 1) & phytn15flag = 0.5 + sign(0.5, biophytn15 - trcmin) if (zoopn15flag .eq. 1) & zoopn15flag = 0.5 + sign(0.5, biozoopn15 - trcmin) if (detrn15flag .eq. 1) & detrn15flag = 0.5 + sign(0.5, biodetrn15 - trcmin) if (diazn15flag .eq. 1) & diazn15flag = 0.5 + sign(0.5, biodiazn15 - trcmin) # endif # endif # if defined O_carbon_13 if (dic13flag .eq. 1) & dic13flag = 0.5 + sign(0.5,biodic13 - trcmin) if (doc13flag .eq. 1) & doc13flag = 0.5 + sign(0.5,biodoc13 - trcmin) if (phytc13flag .eq. 1) & phytc13flag = 0.5 + sign(0.5,biophytc13 - trcmin) if (zoopc13flag .eq. 1) & zoopc13flag = 0.5 + sign(0.5,biozoopc13 - trcmin) if (detrc13flag .eq. 1) & detrc13flag = 0.5 + sign(0.5,biodetrc13 - trcmin) # if defined O_npzd_nitrogen if (diazc13flag .eq. 1) & diazc13flag = 0.5 + sign(0.5,biodiazc13 - trcmin) # endif # endif enddo bioout(1) = biopo4 - bioin(1) bioout(2) = biodop - bioin(2) bioout(3) = biophyt - bioin(3) bioout(4) = biozoop - bioin(4) bioout(5) = biodetr - bioin(5) # if defined O_npzd_nitrogen bioout(6) = biodic - bioin(6) bioout(7) = biono3 - bioin(7) bioout(8) = biodon - bioin(8) bioout(9) = biodiaz - bioin(9) # if defined O_npzd_nitrogen_15 bioout(10) = biodin15 - bioin(10) bioout(11) = biodon15 - bioin(11) bioout(12) = biophytn15 - bioin(12) bioout(13) = biozoopn15 - bioin(13) bioout(14) = biodetrn15 - bioin(14) bioout(15) = biodiazn15 - bioin(15) # endif # endif # if defined O_carbon_13 bioout(16) = biodic13 - bioin(16) bioout(17) = biodoc13 - bioin(17) bioout(18) = biophytc13 - bioin(18) bioout(19) = biozoopc13 - bioin(19) bioout(20) = biodetrc13 - bioin(20) bioout(21) = biodiazc13 - bioin(21) # endif #endif return end #if defined O_npzd_fe_limitation subroutine fe_limit(i,j,k,felimit,felimit_D) implicit none integer i,j,k,m,fe_n,fe_jlo,fe_m,fe_k parameter (fe_n = 14) real fe_x(fe_n), fe_y(fe_n), fe_dy, fe_conc real felimit, felimit_D include "size.h" include "npzd.h" include "tmngr.h" if (k.le.3) then ! create x (time) and y (data) arrays for fe interpolation ! first zero them out just in case fe_y(:) = 0 fe_x(:) = 0 ! write in values for days 0 and 365 for interp. bounds fe_y(1) = fe_dissolved(i,j,k,1) fe_x(1) = 0 fe_y(14) = fe_dissolved(i,j,k,12) fe_x(14) = 365 ! now create the rest of the array based on the BLING ! data which is from day 16 of each month do m=2,13 fe_y(m) = fe_dissolved(i,j,k,m-1) if (m.eq.2) then fe_x(2) = 16 else if (m.eq.3) then fe_x(3) = 47 else if (m.eq.4) then fe_x(4) = 75 else if (m.eq.5) then fe_x(5) = 106 else if (m.eq.6) then fe_x(6) = 136 else if (m.eq.7) then fe_x(7) = 167 else if (m.eq.8) then fe_x(8) = 197 else if (m.eq.9) then fe_x(9) = 228 else if (m.eq.10) then fe_x(10) = 259 else if (m.eq.11) then fe_x(11) = 289 else if (m.eq.12) then fe_x(12) = 320 else if (m.eq.13) then fe_x(13) = 350 endif enddo fe_jlo = 2 fe_m = 4 ! find fe data day index call hunt (fe_x,fe_n,dayoyr,fe_jlo) ! initialize the fe data array at the right day fe_k = min(max(fe_jlo-(fe_m-1)/2,1),fe_n+1-fe_m) ! interpolate the fe data, note this does not use the whole array call polint (fe_x(fe_k),fe_y(fe_k) &, fe_m,dayoyr,fe_conc,fe_dy) ! calculate the fe limitation term felimit = fe_conc/(kfe + fe_conc) felimit_D = fe_conc/(kfe_D + fe_conc) else felimit = 1 felimit_D = 1 end if ! print*,'felimit=',felimit return end subroutine polint (xa,ya,n,x,y,dy) ! implicit none ! integer n,i,m,ns, nmax parameter (nmax = 10) ! largest anticipated value of n real dy,x,y,xa(n),ya(n) real den,dif,dift,ho,hp,w,c(nmax),d(nmax) ! ! Given arrays xa and ya, each of length n, and a given value x, this routine ! returns a value y, and an error estimate dy. If P(x) is the polynomial of ! degree N-1 such that P(xai) = yai, i=1,....,n then the returned value y=P(x) ! ! ns=1 dif=abs(x-xa(1)) ! do i=1,n dift = abs(x-xa(i)) IF (dift.lt.dif) THEN ns = i dif = dift END IF ! c(i) = ya(i) d(i) = ya(i) ! enddo ! y = ya (ns) ! this is the initial approximation to y ns = ns-1 ! do m=1,n-1 do i=1,n-m ho = xa(i) - x hp = xa (i+m) -x w = c(i+1) - d(i) den = ho - hp IF (den.eq.0.) PAUSE 'failure in polint' den = w/den d(i) = hp*den c(i) = ho*den enddo IF (2*ns.lt.n-m) THEN dy = c(ns+1) ELSE dy = d(ns) ns = ns - 1 END IF y = y + dy enddo return end ! subroutine hunt (xa,n,x,jlo) ! integer jlo,n,inc,jhi,jm real x, xa(n) logical ascnd ! ! Given an array xa(1:n), and given a value x, returns a value jlo ! such that x is between xa(jlo) and xa(jlo +1). xa(1:n) must be ! monotonic, either increasing or decreasing. jlo=0 or jlo=n is ! returned to indicate that x is out of range. jlo on input is taken ! as the initial guess for jlo on output ! ascnd = xa(n).ge.xa(1) ! True if ascending order of table, false otherwise IF (jlo.le.0.or.jlo.gt.n) THEN ! Input guess not useful. Go to bisection jlo = 0 jhi = n+1 goto 3 END IF inc = 1 ! Set the hunting increment IF (x.ge.xa(jlo).eqv.ascnd) THEN ! Hunting up: 1 jhi = jlo + inc IF (jhi.gt.n) THEN ! Done hunting, since off end of table jhi = n + 1 ELSE IF (x.ge.xa(jhi).eqv.ascnd) THEN ! Not done hunting jlo = jhi inc = inc + inc ! so double the increment goto 1 ! and try again END IF ! Done hunting, value bracketed. ELSE ! Hunt down: jhi = jlo 2 jlo = jhi-inc IF (jlo.lt.1) THEN !Done hunting, since off end of table jlo = 0 ELSE IF (x.lt.xa(jlo).eqv.ascnd) THEN ! Not done hunting jhi = jlo inc = inc + inc ! so double the increment goto 2 ! and try again END IF ! Done hunting, value bracketed END IF ! Hunt is done, so begin the final bisection phase: ! 3 IF (jhi-jlo.eq.1) THEN IF (x.eq.xa(n)) jlo = n-1 IF (x.eq.xa(1)) jlo = 1 RETURN END IF jm = (jhi + jlo)/2 IF (x.ge.xa(jm).eqv.ascnd) THEN jlo = jm ELSE jhi = jm END IF ! GOTO 3 ! end #endif