Merge branch 'master' into 72-nitrous-oxide

src/bacteria_docdyn.F90

@@ -14,15 +14,17 @@ module ersem_bacteria_docdyn
 
    type,extends(type_ersem_pelagic_base),public :: type_ersem_bacteria_docdyn
       ! Variables
-      type (type_state_variable_id) :: id_O3c, id_O2o, id_TA
+      type (type_state_variable_id) :: id_O3c, id_O2o, id_TA, id_N3n
       type (type_state_variable_id) :: id_R1c, id_R2c, id_R3c
       type (type_state_variable_id) :: id_R1p
       type (type_state_variable_id) :: id_R1n
-      type (type_state_variable_id) :: id_N1p,id_N4n,id_N7f
+      type (type_state_variable_id) :: id_N1p,id_N4n,id_N7f,id_N6
       type (type_dependency_id)     :: id_ETW,id_eO2mO2
       type (type_state_variable_id),allocatable,dimension(:) :: id_RPc,id_RPp,id_RPn,id_RPf
       type (type_model_id),         allocatable,dimension(:) :: id_RP
-      type (type_diagnostic_variable_id) :: id_fB1O3c, id_fB1NIn, id_fB1N1p, id_bgeff
+
+      type (type_diagnostic_variable_id) :: id_fB1O3c, id_fB1NIn, id_fB1N1p,id_bgeff,id_fdenit,id_fanox,id_freox
+
       type (type_diagnostic_variable_id) :: id_fR1B1c, id_fR2B1c, id_fR3B1c,id_fRPB1c,id_fB1R1c, id_fB1R2c, id_fB1R3c
       type (type_diagnostic_variable_id) :: id_fR1B1n,id_fB1R1n,id_fR1B1p,id_fB1R1p,id_fRPB1n,id_fRPB1p
       type (type_diagnostic_variable_id) :: id_minn,id_minp
@@ -39,6 +41,8 @@ module ersem_bacteria_docdyn
       real(rk) :: rR2B1X,rR3B1X
       real(rk),allocatable :: sRPR1(:)
       real(rk) :: frB1R3
+      real(rk) :: DeniX,reoX,omroX,omonX,chN3oX
+      integer  :: denit
 
       ! Remineralization
       real(rk) :: sR1N1X,sR1N4X
@@ -84,6 +88,22 @@ subroutine initialize(self,configunit)
       call self%get_parameter(self%qpB1cX,  'qpc',     'mmol P/mg C','maximum phosphorus to carbon ratio')
       call self%get_parameter(self%qnB1cX,  'qnc',     'mmol N/mg C','maximum nitrogen to carbon ratio')
       call self%get_parameter(self%urB1_O2X,'ur_O2',   'mmol O_2/mg C','oxygen consumed per carbon respired')
+      call self%get_parameter(self%denit,   'denit',   '-',           'denitrification switch (0: off, 1: on)',default=0)
+
+      if (self%denit == 1) then
+      call self%get_parameter(self%DeniX,   'DeniX',   '1/d',         'specific denitrification rate')
+      call self%get_parameter(self%reoX,    'reoX',    '1/d',         'specific reoxidation rate of reduction equivalents')
+      call self%get_parameter(self%omroX,   'omroX',   'mmol HS mmol O2-1', 'stoichiometric coefficient')
+      call self%get_parameter(self%omonX,   'omonX',   'mmol O2 mmol N-1',  'stoichiometric coefficient for denitrification reaction')
+      call self%get_parameter(self%chN3oX,  'chN3o',   '(mmol O_2/m^3)^3', 'Michaelis-Menten constant for cubic oxygen dependence of nitrification')
+
+      call self%register_state_dependency(self%id_N3n,'N3n','mmol N/m^3','nitrate')
+      call self%register_state_dependency(self%id_N6,'N6','mmol HS-/m^3','reduction equivalent')
+
+      call self%register_diagnostic_variable(self%id_fdenit,'fdenit','mmol N/m^3/d','denitrification', missing_value=0._rk)
+      call self%register_diagnostic_variable(self%id_fanox,'fanox',  '-',           'fanox',           missing_value=0._rk)
+      call self%register_diagnostic_variable(self%id_freox,'freox',  '-',           'freox',           missing_value=0._rk)
+      end if
 
       ! Remineralization parameters
       call self%get_parameter(self%sR1N1X,   'sR1N1',   '1/d',    'mineralisation rate of labile dissolved organic phosphorus')
@@ -189,7 +209,7 @@ subroutine do(self,_ARGUMENTS_DO_)
       real(rk) :: ETW,eO2mO2
       real(rk) :: B1c,B1n,B1p
       real(rk) :: B1cP,B1nP,B1pP
-      real(rk) :: N1pP,N4nP,R1c,R1cP,R1pP,R1nP,R2c
+      real(rk) :: N1pP,N4nP,R1c,R1cP,R1pP,R1nP,R2c,N3n,N6,O2o
       real(rk) :: qpB1c,qnB1c
       real(rk) :: etB1,eO2B1
       real(rk) :: sB1RD,sutB1,rumB1,sugB1,rugB1,rraB1,fB1O3c
@@ -201,6 +221,7 @@ subroutine do(self,_ARGUMENTS_DO_)
       real(rk) :: fB1R1c
       real(rk) :: totsubst
       real(rk) :: CORROX
+      real(rk) :: fdenit,denitpot,deniteff,o2state,freox,fanox
       integer  :: iRP
       real(rk),dimension(self%nRP) :: RPc,RPcP,RPnP,RPpP
       real(rk),dimension(self%nRP) :: fRPB1c,fRPB1p,fRPB1n
@@ -239,6 +260,15 @@ subroutine do(self,_ARGUMENTS_DO_)
 
          qpB1c = B1p/B1c
          qnB1c = B1n/B1c
+
+        if (self%denit == 1) then   
+            _GET_(self%id_N3n,N3n)
+            _GET_(self%id_O2o,O2o)
+            O2o = max(0.0_rk,O2o)
+            o2state = O2o**3/(O2o**3 + self%chN3oX)
+            _GET_(self%id_N6,N6)
+        end if
+
 !..Temperature effect on pelagic bacteria:
 
          etB1 = max(0.0_rk,self%q10B1X**((ETW-10._rk)/10._rk) - self%q10B1X**((ETW-32._rk)/3._rk))
@@ -294,6 +324,25 @@ subroutine do(self,_ARGUMENTS_DO_)
          fB1R3c=self%frB1R3*rraB1
          fB1RDc = fB1R1c + fB1R2c + fB1R3c
 
+! Denitrification implemented as in Sankar et al. (2008), doi.org/10.1016/j.ecolmodel.2018.01.016
+     if (self%denit == 1) then      
+        denitpot = self%DeniX * N3n !eq.6
+        deniteff = max(0._rk, self%urB1_O2X * (1._rk-o2state) * fB1O3c / self%omonX) !eq.7
+        fdenit = min(denitpot, deniteff) !eq.5
+
+        _SET_DIAGNOSTIC_(self%id_fdenit,fdenit)
+        _SET_ODE_(self%id_N3n, -fdenit)
+
+! Reduced sulfur formation corresponds to eq.9 in Sankar et al. (2008)
+        fanox = self%omrox * (self%urB1_O2X * (1._rk-o2state) * fB1O3c - self%omonX * fdenit)
+        freox = self%reoX * etB1 * o2state * N6
+
+        _SET_DIAGNOSTIC_(self%id_fanox,freox)
+        _SET_DIAGNOSTIC_(self%id_freox,fanox)
+
+        _SET_ODE_(self%id_N6, fanox - freox)
+     end if 
+
 !..net bacterial production
 
          netb1 = rugB1 - fB1o3c - fB1RDc

src/pelagic_base.F90

@@ -13,7 +13,7 @@ module ersem_pelagic_base
    private
 
    type,extends(type_particle_model),public :: type_ersem_pelagic_base
-      type (type_state_variable_id)                 :: id_c,id_n,id_p,id_f,id_s,id_chl
+      type (type_state_variable_id)                 :: id_c,id_n,id_p,id_f,id_s,id_chl,id_h
       type (type_horizontal_dependency_id)          :: id_bedstress,id_wdepth
       type (type_dependency_id)                     :: id_dens
       type (type_horizontal_diagnostic_variable_id) :: id_w_bot
@@ -92,6 +92,7 @@ subroutine initialize(self,configunit)
          call self%add_constituent('s',0.0_rk,s0)
       end if
       if (index(composition,'f')/=0) call self%add_constituent('f',0.0_rk)
+      if (index(composition,'h')/=0) call self%add_constituent('h',0.0_rk)
 
    end subroutine
 
@@ -157,6 +158,8 @@ subroutine add_constituent(self,name,initial_value,background_value,qn,qp)
          if (use_iron) call register(self%id_f,'f','umol Fe','iron',standard_variables%total_iron,self%qxf,self%id_fdep,self%id_targetf)
       case ('chl')
          call register(self%id_chl,'Chl','mg','chlorophyll a',total_chlorophyll)
+      case ('h')
+         call register(self%id_h,'h','mmol','reduction equivalent',total_h)
       case default
          call self%fatal_error('add_constituent','Unknown constituent "'//trim(name)//'".')
       end select

src/shared.F90

@@ -26,6 +26,7 @@ module ersem_shared
 
    ! Aggregate diagnostics for e.g., carbon budgets.
    type (type_bulk_standard_variable),parameter :: total_chlorophyll = type_bulk_standard_variable(name='total_chlorophyll',units='mg/m^3',aggregate_variable=.true.)
+   type (type_bulk_standard_variable),parameter :: total_h = type_bulk_standard_variable(name='total_h',units='mmol/m^3',aggregate_variable=.true.)
    type (type_bulk_standard_variable),parameter :: total_calcite_in_biota = type_bulk_standard_variable(name='total_calcite_in_biota',units='mg C/m^3',aggregate_variable=.true.)
    type (type_bulk_standard_variable),parameter :: secchi_depth = type_bulk_standard_variable(name='secchi_depth',units='m')