Skip to content

Commit

Permalink
Bugfix DNRA and introduce max_limiter (#413)
Browse files Browse the repository at this point in the history 
* Bugfix DNRA and introduce max_limiter
  • Loading branch information
@jmaerz @TomasTorsvik
jmaerz authored and TomasTorsvik committed Oct 23, 2024
1 parent b1d38c3 commit 3c132d5
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 41 deletions.
34 changes: 18 additions & 16 deletions hamocc/mo_extNsediment.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -91,7 +91,7 @@ module mo_extNsediment
ised_remin_sulf = 13, &
n_seddiag = 13

real :: eps = 1.e-25
real :: eps = epsilon(1.)
real :: minlim = 1.e-9

contains
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
41 changes: 21 additions & 20 deletions hamocc/mo_extNwatercol.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -73,8 +72,7 @@ module mo_extNwatercol
! public functions
public :: nitrification,denit_NO3_to_NO2,anammox,denit_dnra,extN_inv_check

real :: eps = 1.e-25
real :: minlim = 1.e-9
real :: eps = epsilon(1.)

contains

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
7 changes: 3 additions & 4 deletions hamocc/mo_ocprod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -333,14 +333,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+epsilon(1.))
! 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))
Expand Down
4 changes: 3 additions & 1 deletion hamocc/mo_param_bgc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,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


!********************************************************************
Expand Down Expand Up @@ -151,6 +152,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 be 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 !
Expand Down

0 comments on commit 3c132d5

Please sign in to comment.