Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Bugfix DNRA and introduce max_limiter #413

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
+44 −39
Open

Bugfix DNRA and introduce max_limiter #413

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 17 additions & 15 deletions hamocc/mo_extNsediment.F90
View file
Open in desktop
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 @@ -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
39 changes: 20 additions & 19 deletions hamocc/mo_extNwatercol.F90
View file
Open in desktop
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 @@ -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

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
8 changes: 4 additions & 4 deletions hamocc/mo_ocprod.F90
View file
Open in desktop
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 @@ -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
Expand Down Expand Up @@ -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))
Expand Down
4 changes: 3 additions & 1 deletion hamocc/mo_param_bgc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


!********************************************************************
Expand Down Expand Up @@ -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 !
Expand Down