From eb67a7a16db95486daa8abacaa03977008851827 Mon Sep 17 00:00:00 2001 From: joeran maerz Date: Fri, 18 Oct 2024 20:46:10 +0200 Subject: [PATCH] Bugfix DNRA and introduce max_limiter --- hamocc/mo_extNsediment.F90 | 32 ++++++++++++++++--------------- hamocc/mo_extNwatercol.F90 | 39 +++++++++++++++++++------------------- hamocc/mo_ocprod.F90 | 8 ++++---- hamocc/mo_param_bgc.F90 | 4 +++- 4 files changed, 44 insertions(+), 39 deletions(-) diff --git a/hamocc/mo_extNsediment.F90 b/hamocc/mo_extNsediment.F90 index 546ae10e..42871395 100644 --- a/hamocc/mo_extNsediment.F90 +++ b/hamocc/mo_extNsediment.F90 @@ -56,8 +56,8 @@ module mo_extNsediment & bkanh4nitr_sed,bkamoxn2o_sed,bkyamox_sed,rano2nitr_sed,q10ano2nitr_sed,& & Trefano2nitr_sed,bkoxnitr_sed,bkano2nitr_sed,n2omaxy_sed,n2oybeta_sed, & & NOB2AOAy_sed,bn2o_sed,mufn2o_sed,POM_remin_q10_sed, POM_remin_Tref_sed,& - & bkox_drempoc_sed - use mo_control_bgc, only: io_stdo_bgc,dtb + & bkox_drempoc_sed,max_limiter + use mo_control_bgc, only: io_stdo_bgc use mo_sedmnt, only: powtra,sedlay,porsol,porwat implicit none @@ -195,16 +195,17 @@ subroutine sed_nitrification(j,kpie,kpje,kpke,kbnd,ptho,omask,ex_ddic,ex_dalk) ! Account for potential earlier changes in DIC and alkalinity in finiding the minimum totd = max(0., & - & min(totd, & - & powtra(i,j,k,ipownh4)/(amoxfrac + fdetnitr*nitrfrac + eps), & ! ammonium - & (powtra(i,j,k,ipowaic) + ex_ddic(i,k)) & + & min(totd, & + & max_limiter*powtra(i,j,k,ipownh4)/(amoxfrac + fdetnitr*nitrfrac + eps), & ! ammonium + & max_limiter*(powtra(i,j,k,ipowaic) + ex_ddic(i,k)) & & /(rc2n*(fdetamox*amoxfrac + fdetnitr*nitrfrac) & & + eps), & ! CO2 - & powtra(i,j,k,ipowaph)/(rnoi*(fdetamox*amoxfrac+fdetnitr*nitrfrac) + eps), & ! PO4 - & powtra(i,j,k,ipowaox) & + & max_limiter*powtra(i,j,k,ipowaph) & + & /(rnoi*(fdetamox*amoxfrac+fdetnitr*nitrfrac) + eps), & ! PO4 + & max_limiter*powtra(i,j,k,ipowaox) & & /((1.5*fno2 + fn2o - ro2nnit*fdetamox)*amoxfrac & & + (0.5 - ro2nnit*fdetnitr)*nitrfrac + eps), & ! O2 - & (powtra(i,j,k,ipowaal) + ex_dalk(i,k)) & + & max_limiter*(powtra(i,j,k,ipowaal) + ex_dalk(i,k)) & & /((2.*fno2 + fn2o + rnm1*rnoi*fdetamox)*amoxfrac & & + (rnm1*rnoi*fdetnitr)*nitrfrac + eps))) ! alkalinity amox = amoxfrac*totd @@ -260,7 +261,8 @@ subroutine sed_denit_NO3_to_NO2(j,kpie,kpje,kpke,kbnd,ptho,omask,ex_ddic,ex_dalk ano3new = powtra(i,j,k,ipowno3)/(1. + rano3denit_sed*Tdep*O2inhib*nutlim) - ano3denit = max(0.,min(powtra(i,j,k,ipowno3) - ano3new, sedlay(i,j,k,issso12)*rnoxp*s2w)) + ano3denit = max(0.,min(powtra(i,j,k,ipowno3) - ano3new, & + & max_limiter*sedlay(i,j,k,issso12)*rnoxp*s2w)) powtra(i,j,k,ipowno3) = powtra(i,j,k,ipowno3) - ano3denit powtra(i,j,k,ipowno2) = powtra(i,j,k,ipowno2) + ano3denit @@ -307,11 +309,11 @@ subroutine sed_anammox(j,kpie,kpje,kpke,kbnd,ptho,omask,ex_ddic,ex_dalk) ano2new = powtra(i,j,k,ipowno2)/(1. + rano2anmx_sed*Tdep*O2inhib*nut1lim*nut2lim) ! Account for former changes in DIC and alkalinity - ano2anmx = max(0.,min(powtra(i,j,k,ipowno2) - ano2new, & - powtra(i,j,k,ipownh4)*rno2anmx*rnh4anmxi, & - (powtra(i,j,k,ipowaic)+ex_ddic(i,k))*rno2anmx/rcar, & - powtra(i,j,k,ipowaph)*rno2anmx, & - (powtra(i,j,k,ipowaal)+ex_dalk(i,k))*rno2anmx/rnm1)) + ano2anmx = max(0.,min(max_limiter*powtra(i,j,k,ipowno2) - ano2new, & + max_limiter*powtra(i,j,k,ipownh4)*rno2anmx*rnh4anmxi, & + max_limiter*(powtra(i,j,k,ipowaic)+ex_ddic(i,k))*rno2anmx/rcar, & + max_limiter*powtra(i,j,k,ipowaph)*rno2anmx, & + max_limiter*(powtra(i,j,k,ipowaal)+ex_dalk(i,k))*rno2anmx/rnm1)) powtra(i,j,k,ipowno2) = powtra(i,j,k,ipowno2) - ano2anmx powtra(i,j,k,ipownh4) = powtra(i,j,k,ipownh4) - ano2anmx*rnh4anmx*rno2anmxi @@ -398,7 +400,7 @@ subroutine sed_denit_DNRA(j,kpie,kpje,kpke,kbnd,ptho,omask,ex_ddic,ex_dalk) fdetano2denit = rnoxpi*ano2denit/(potddet + eps) fdetan2odenit = rnoxpi*an2odenit/(potddet + eps) fdetdnra = 1. - fdetano2denit - fdetan2odenit - potddet = max(0.,min(potddet,powtra(i,j,k,issso12)*s2w)) + potddet = max(0.,min(potddet,max_limiter*sedlay(i,j,k,issso12)*s2w)) ! change of NO2 and N2O in N units ano2denit = fdetano2denit*rnoxp*potddet diff --git a/hamocc/mo_extNwatercol.F90 b/hamocc/mo_extNwatercol.F90 index 3995c98d..7b8a7846 100644 --- a/hamocc/mo_extNwatercol.F90 +++ b/hamocc/mo_extNwatercol.F90 @@ -48,7 +48,6 @@ module mo_extNwatercol !**************************************************************** use mo_vgrid, only: dp_min use mod_xc, only: mnproc - use mo_control_bgc, only: dtb use mo_param1_bgc, only: ialkali,ianh4,iano2,ian2o,iano3,idet,igasnit,iiron,ioxygen,iphosph, & & isco212 use mo_carbch, only: ocetra @@ -63,7 +62,7 @@ module mo_extNwatercol & n2oybeta,NOB2AOAy,bn2o,mufn2o, & & rc2n,ro2nnit,rnoxp,rnoxpi,rno2anmx,rno2anmxi,rnh4anmx, & & rnh4anmxi,rno2dnra,rno2dnrai,rnh4dnra,rnh4dnrai,rnm1, & - & bkphyanh4,bkphyano3,bkphosph,bkiron,ro2utammo + & bkphyanh4,bkphyano3,bkphosph,bkiron,ro2utammo,max_limiter use mo_biomod, only: nitr_NH4,nitr_NO2,nitr_N2O_prod,nitr_NH4_OM,nitr_NO2_OM,denit_NO3, & & denit_NO2,denit_N2O,DNRA_NO2,anmx_N2_prod,anmx_OM_prod implicit none @@ -74,7 +73,6 @@ module mo_extNwatercol public :: nitrification,denit_NO3_to_NO2,anammox,denit_dnra,extN_inv_check real :: eps = 1.e-25 - real :: minlim = 1.e-9 contains @@ -164,15 +162,17 @@ subroutine nitrification(kpie,kpje,kpke,kbnd,pddpo,omask,ptho) totd = max(0., & & min(totd, & - & ocetra(i,j,k,ianh4)/(amoxfrac + fdetnitr*nitrfrac + eps), & ! ammonium - & ocetra(i,j,k,isco212)/(rc2n*(fdetamox*amoxfrac + fdetnitr*nitrfrac) +eps),& ! CO2 - & ocetra(i,j,k,iphosph)/(rnoi*(fdetamox*amoxfrac + fdetnitr*nitrfrac) +eps),& ! PO4 - & ocetra(i,j,k,iiron)/(riron*rnoi*(fdetamox*amoxfrac + fdetnitr*nitrfrac) & - & + eps), & ! Fe - & ocetra(i,j,k,ioxygen) & + & max_limiter*ocetra(i,j,k,ianh4)/(amoxfrac + fdetnitr*nitrfrac + eps), & ! ammonium + & max_limiter*ocetra(i,j,k,isco212)/ & + & (rc2n*(fdetamox*amoxfrac + fdetnitr*nitrfrac) +eps),& ! CO2 + & max_limiter*ocetra(i,j,k,iphosph)/ & + & (rnoi*(fdetamox*amoxfrac + fdetnitr*nitrfrac) +eps),& ! PO4 + & max_limiter*ocetra(i,j,k,iiron)/ & + & (riron*rnoi*(fdetamox*amoxfrac + fdetnitr*nitrfrac) + eps), & ! Fe + & max_limiter*ocetra(i,j,k,ioxygen) & & /((1.5*fno2 + fn2o - ro2nnit*fdetamox)*amoxfrac & + (0.5 - ro2nnit*fdetnitr)*nitrfrac + eps), & ! O2 - & ocetra(i,j,k,ialkali) & + & max_limiter*ocetra(i,j,k,ialkali) & & /((2.*fno2 + fn2o + rnm1*rnoi*fdetamox)*amoxfrac & & + (rnm1*rnoi*fdetnitr)*nitrfrac + eps))) ! alkalinity amox = amoxfrac*totd @@ -229,10 +229,11 @@ subroutine denit_NO3_to_NO2(kpie,kpje,kpke,kbnd,pddpo,omask,ptho) Tdep = q10ano3denit**((temp-Trefano3denit)/10.) O2inhib = 1. - tanh(sc_ano3denit*ocetra(i,j,k,ioxygen)) nutlim = ocetra(i,j,k,iano3)/(ocetra(i,j,k,iano3) + bkano3denit) - + ano3new = ocetra(i,j,k,iano3)/(1. + rano3denit*Tdep*O2inhib*nutlim) - ano3denit = max(0.,min(ocetra(i,j,k,iano3) - ano3new, ocetra(i,j,k,idet)*rnoxp)) + ano3denit = max(0.,min(ocetra(i,j,k,iano3) - ano3new, & + max_limiter*ocetra(i,j,k,idet)*rnoxp)) ocetra(i,j,k,iano3) = ocetra(i,j,k,iano3) - ano3denit ocetra(i,j,k,iano2) = ocetra(i,j,k,iano2) + ano3denit @@ -284,11 +285,11 @@ subroutine anammox(kpie,kpje,kpke,kbnd,pddpo,omask,ptho) ano2new = ocetra(i,j,k,iano2)/(1. + rano2anmx*Tdep*O2inhib*nut1lim*nut2lim) ano2anmx = max(0.,min(ocetra(i,j,k,iano2) - ano2new, & - ocetra(i,j,k,ianh4)*rno2anmx*rnh4anmxi, & - ocetra(i,j,k,isco212)*rno2anmx/rcar, & - ocetra(i,j,k,iphosph)*rno2anmx, & - ocetra(i,j,k,iiron)*rno2anmx/riron, & - ocetra(i,j,k,ialkali)*rno2anmx/rnm1)) + max_limiter*ocetra(i,j,k,ianh4)*rno2anmx*rnh4anmxi, & + max_limiter*ocetra(i,j,k,isco212)*rno2anmx/rcar, & + max_limiter*ocetra(i,j,k,iphosph)*rno2anmx, & + max_limiter*ocetra(i,j,k,iiron)*rno2anmx/riron, & + max_limiter*ocetra(i,j,k,ialkali)*rno2anmx/rnm1)) ocetra(i,j,k,iano2) = ocetra(i,j,k,iano2) - ano2anmx ocetra(i,j,k,ianh4) = ocetra(i,j,k,ianh4) - ano2anmx*rnh4anmx*rno2anmxi @@ -387,7 +388,7 @@ subroutine denit_dnra(kpie,kpje,kpke,kbnd,pddpo,omask,ptho) fdetano2denit = rnoxpi*ano2denit/(potddet + eps) fdetan2odenit = rnoxpi*an2odenit/(potddet + eps) fdetdnra = 1. - fdetano2denit - fdetan2odenit - potddet = max(0.,min(potddet,ocetra(i,j,k,idet))) + potddet = max(0.,min(potddet,max_limiter*ocetra(i,j,k,idet))) ! change of NO2 and N2O in N units ano2denit = fdetano2denit*rnoxp*potddet @@ -426,7 +427,7 @@ end subroutine denit_dnra !================================================================================================================================== subroutine extN_inv_check(kpie,kpje,kpke,pdlxp,pdlyp,pddpo,omask,inv_message) use mo_inventory_bgc, only: inventory_bgc - use mo_control_bgc, only: io_stdo_bgc,dtb,use_PBGC_OCNP_TIMESTEP + use mo_control_bgc, only: io_stdo_bgc,use_PBGC_OCNP_TIMESTEP implicit none ! provide inventory calculation for extended nitrogen cycle diff --git a/hamocc/mo_ocprod.F90 b/hamocc/mo_ocprod.F90 index 90856f94..8224962d 100644 --- a/hamocc/mo_ocprod.F90 +++ b/hamocc/mo_ocprod.F90 @@ -77,7 +77,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph,psao,pp dmsp1,dmsp2,dmsp3,dmsp4,dmsp5,dmsp6,dms_gamma, & fbro1,fbro2,atten_f,atten_c,atten_uv,atten_w,bkopal,bkphy,bkzoo, & POM_remin_q10,POM_remin_Tref,opal_remin_q10,opal_remin_Tref, & - bkphyanh4,bkphyano3,bkphosph,bkiron,ro2utammo + bkphyanh4,bkphyano3,bkphosph,bkiron,ro2utammo,max_limiter use mo_biomod, only: bsiflx0100,bsiflx0500,bsiflx1000,bsiflx2000,bsiflx4000,bsiflx_bot, & calflx0100,calflx0500,calflx1000,calflx2000,calflx4000,calflx_bot, & carflx0100,carflx0500,carflx1000,carflx2000,carflx4000,carflx_bot, & @@ -137,6 +137,7 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph,psao,pp real :: wpocd,wcald,wopald,wdustd,dagg real :: wcal,wdust,wopal,wpoc real :: o2lim ! O2 limitation of ammonification (POC remin) + real, parameter :: tiny_num = 1e-25 ! sedbypass real :: florca,flcaca,flsil ! cisonew @@ -333,14 +334,13 @@ subroutine ocprod(kpie,kpje,kpke,kbnd,pdlxp,pdlyp,pddpo,omask,ptho,pi_ph,psao,pp nlim = ano3up_inh*ocetra(i,j,k,iano3)/(ocetra(i,j,k,iano3) + bkphyano3) + anh4lim grlim = min(nutlim,nlim) ! growth limitation - nh4uptfrac = 1. - if(nlim .gt. 1.e-18) nh4uptfrac = anh4lim/nlim + nh4uptfrac = anh4lim/(nlim+tiny_num) ! re-check avnut - can sum N avail exceed indiv. contrib? avanut = max(0.,min(ocetra(i,j,k,iphosph), ocetra(i,j,k,iiron)/riron, & & rnoi*((1.-nh4uptfrac)*ocetra(i,j,k,iano3) + nh4uptfrac*ocetra(i,j,k,ianh4)))) xn = avphy/(1. - pho*grlim) ! phytoplankton growth - phosy = max(0.,min(xn-avphy,avanut)) ! limit PP growth to available nutr. + phosy = max(0.,min(xn-avphy,max_limiter*avanut)) ! limit PP growth to available nutr. else avanut = max(0.,min(ocetra(i,j,k,iphosph),rnoi*ocetra(i,j,k,iano3))) avanfe = max(0.,min(avanut,ocetra(i,j,k,iiron)/riron)) diff --git a/hamocc/mo_param_bgc.F90 b/hamocc/mo_param_bgc.F90 index 8deac006..4690c1a4 100644 --- a/hamocc/mo_param_bgc.F90 +++ b/hamocc/mo_param_bgc.F90 @@ -121,7 +121,8 @@ module mo_param_bgc & bkoxamox_sed,bkanh4nitr_sed,bkamoxn2o_sed,bkyamox_sed, & & rano2nitr_sed,q10ano2nitr_sed,Trefano2nitr_sed,bkoxnitr_sed, & & bkano2nitr_sed,n2omaxy_sed,n2oybeta_sed,NOB2AOAy_sed,bn2o_sed, & - & mufn2o_sed,POM_remin_q10_sed, POM_remin_Tref_sed,bkox_drempoc_sed + & mufn2o_sed,POM_remin_q10_sed, POM_remin_Tref_sed,bkox_drempoc_sed, & + & max_limiter !******************************************************************** @@ -152,6 +153,7 @@ module mo_param_bgc real, parameter :: c14_t_half = 5700.*365. ! Half life of 14C [days] ! Extended nitrogen cycle + real, parameter :: max_limiter = 0.9999 ! maximum in concentrations that can consumed at once real, parameter :: rc2n = rcar/rnit ! iHAMOCC C:N ratio real, parameter :: ro2utammo = 140. ! Oxygen utilization per mol detritus during ammonification real, parameter :: ro2nnit = ro2utammo/rnit !