Merge pull request #88 from pmlmodelling/merge-dev-to-public

Merge dev to public

src/CMakeLists.txt

@@ -24,6 +24,7 @@ add_library(fabm_models_ersem OBJECT
             carbonate.F90
             oxygen.F90
             nitrification.F90
+            nitrous_oxide.F90
             benthic_column.F90
             benthic_column_dissolved_matter.F90
             benthic_column_particulate_matter.F90

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/ersem_model_library.F90

@@ -14,6 +14,7 @@ module ersem_model_library
    use ersem_bacteria
    use ersem_bacteria_docdyn
    use ersem_nitrification
+   use ersem_nitrous_oxide
    use ersem_light
    use ersem_light_iop
    use ersem_light_iop_ady
@@ -65,6 +66,7 @@ subroutine create(self,name,model)
          case ('bacteria');                                allocate(type_ersem_bacteria::model)
          case ('bacteria_docdyn');                         allocate(type_ersem_bacteria_docdyn::model)
          case ('nitrification');                           allocate(type_ersem_nitrification::model)
+         case ('nitrous_oxide');                           allocate(type_ersem_nitrous_oxide::model)
          case ('light');                                   allocate(type_ersem_light::model)
          case ('light_iop');                               allocate(type_ersem_light_iop::model)
          case ('light_iop_ady');                           allocate(type_ersem_light_iop_ady::model)

src/nitrification.F90

@@ -15,14 +15,15 @@ module ersem_nitrification
    type,extends(type_ersem_pelagic_base),public :: type_ersem_nitrification
       ! Variables
       type (type_state_variable_id) :: id_O2o,id_TA
-      type (type_state_variable_id) :: id_N3n,id_N4n
+      type (type_state_variable_id) :: id_N3n,id_N4n,id_O5n
       type (type_dependency_id)     :: id_ETW,id_phx
-      type (type_diagnostic_variable_id)     :: id_nitrification
+      type (type_diagnostic_variable_id)     :: id_nitrification, id_fN4O5n, id_fN4N3n
 
       ! Parameters
       real(rk) :: q10
-      real(rk) :: sN4N3X,chN3oX,chN4nX
+      real(rk) :: sN4N3X,chN3oX,chN4nX, N4O5minX
       integer  :: ISWphx
+      logical  :: ISWn2o
    contains
 !     Model procedures
       procedure :: initialize
@@ -50,15 +51,20 @@ subroutine initialize(self,configunit)
       call self%get_parameter(self%ISWphx,'ISWph','',                'pH impact on nitrification (0: off, 1: on)')
       call self%get_parameter(self%sN4N3X,'sN4N3','1/d',             'specific nitrification rate')
       call self%get_parameter(self%chN3oX,'chN3o','(mmol O_2/m^3)^3','Michaelis-Menten constant for cubic oxygen dependence of nitrification')
+      call self%get_parameter(self%ISWn2o,'ISWn2o','',               'activate n2o production', default = .false.)
+      if (self%ISWn2o) call self%get_parameter(self%N4O5minX,'N4O5minX','-','minimal fraction of N2O production')
       call self%get_parameter(self%chN4nX,'chN4n','(mmol N/m^3)^3','Michaelis-Menten constant for cubic ammonium dependence of nitrification', default=0.0_rk)
 
       ! Register diagnostic variables
       call self%register_diagnostic_variable(self%id_nitrification,"rate","mmol/m3/d","rate")
+      call self%register_diagnostic_variable(self%id_fN4N3n,'fN4N3n','mmol N/m^3/d','nitrification',output=output_time_step_averaged)
+      if (self%ISWn2o) call self%register_diagnostic_variable(self%id_fN4O5n,'fN4O5n','mmol N/m^3/d','N2O production',output=output_time_step_averaged)
 
       ! Register links to nutrient and oxygen pools.
       call self%register_state_dependency(self%id_N3n,'N3n','mmol N/m^3',  'nitrate')
       call self%register_state_dependency(self%id_N4n,'N4n','mmol N/m^3',  'ammonium')
       call self%register_state_dependency(self%id_O2o,'O2o','mmol O_2/m^3','oxygen')
+      if (self%ISWn2o) call self%register_state_dependency(self%id_O5n,'O5n','mmol N/m^3','nitrous oxide')
       call self%register_state_dependency(self%id_TA,standard_variables%alkalinity_expressed_as_mole_equivalent)
 
       ! Register environmental dependencies (temperature, pH)
@@ -73,7 +79,7 @@ subroutine do(self,_ARGUMENTS_DO_)
       _DECLARE_ARGUMENTS_DO_
 
       real(rk) :: ETW,phx,O2o,N4n,N4nP
-      real(rk) :: etB1,o2state,n4state,Fph,fN4N3n
+      real(rk) :: etB1,o2state,n4state,Fph,fN4N3n,fN4O5n
 
       ! Leave spatial loops (if any)
       _LOOP_BEGIN_
@@ -101,6 +107,22 @@ subroutine do(self,_ARGUMENTS_DO_)
             Fph = min(2._rk,max(0._rk,0.6111_rk*phx-3.8889_rk))
             fN4N3n = fN4N3n * Fph
          end if
+         ! N2O production
+         ! A fixed quota of N is assumed to go to N2O during nitrification
+         ! Depending on O2 limitation this quota may increase up to 20 time
+         ! (Cadispoti, 2010). Note that N2O is given as mmol of N m-3 (i.e. 1/2
+         ! mmole N2O m-2)!  (Luca, July 2016)
+
+         if (self%ISWn2o) then
+          fN4O5n=self%N4O5minX*fN4N3n*min(20._rk,(1._rk/o2state))
+         _SET_ODE_(self%id_O5n,fN4O5n)
+         _SET_ODE_(self%id_N3n, -fN4O5n)
+         _SET_DIAGNOSTIC_(self%id_fN4O5n,fN4O5n)
+         end if
+
+         _SET_DIAGNOSTIC_(self%id_fN4N3n,fN4N3n)
+
+
 
          _SET_ODE_(self%id_N3n, + fN4N3n)
          _SET_ODE_(self%id_N4n, - fN4N3n)

src/nitrous_oxide.F90

@@ -0,0 +1,144 @@
+#include "fabm_driver.h"
+
+!------------------------------------------------------------------------
+! This module calculates saturation concentrations and air-sea
+! exchange of nitrous oxide. Production of nitrous oxide is implemented
+! in nitrification module (nitrification.F90).
+! For implementation and validation of nitrous oxide on the North-West 
+! European Shelf see Lessin et al. (2020): doi.org/10.1029/2019JG005613
+!------------------------------------------------------------------------
+
+module ersem_nitrous_oxide
+
+   use fabm_types
+   use fabm_particle
+
+   use ersem_shared
+   use ersem_pelagic_base
+
+   implicit none
+
+   private
+
+   type,extends(type_ersem_pelagic_base),public :: type_ersem_nitrous_oxide
+      ! Variables
+      type (type_dependency_id)                     :: id_ETW,id_X1X
+      type (type_horizontal_dependency_id)          :: id_wnd,id_pN2Oa
+      type (type_horizontal_diagnostic_variable_id) :: id_airsea
+      type (type_diagnostic_variable_id)            :: id_satp
+
+      integer :: iswN2O
+
+   contains
+!     Model procedures
+      procedure :: initialize
+      procedure :: do_surface
+      procedure :: do
+   end type
+
+contains
+
+   subroutine initialize(self,configunit)
+!
+! !DESCRIPTION:
+!
+! !INPUT PARAMETERS:
+      class (type_ersem_nitrous_oxide),intent(inout),target :: self
+      integer,                         intent(in)           :: configunit
+!
+!EOP
+!-----------------------------------------------------------------------
+!BOC
+      call self%initialize_ersem_base(sedimentation=.false.)
+
+      call self%add_constituent('n',0.0_rk)
+
+      call self%register_dependency(self%id_ETW,standard_variables%temperature)
+      call self%register_dependency(self%id_X1X,standard_variables%practical_salinity)
+      call self%register_dependency(self%id_wnd,standard_variables%wind_speed)
+      call self%register_dependency(self%id_pN2Oa,partial_pressure_of_n2o)
+      call self%register_diagnostic_variable(self%id_satp,'satp','%','nitrous oxide % saturation')
+      call self%register_diagnostic_variable(self%id_airsea,'airsea','mmol N2O/m^2/d','airsea flux of N2O',source=source_do_surface)
+      call self%get_parameter(self%iswN2O,'iswN2O','','air-sea flux switch (0: off, 1: on)',default=1)
+
+   end subroutine initialize
+
+   subroutine do(self,_ARGUMENTS_DO_)
+      class (type_ersem_nitrous_oxide), intent(in) :: self
+      _DECLARE_ARGUMENTS_DO_
+
+      real(rk) :: O5n,ETW,X1X
+      real(rk) :: koN2O,pN2Oa
+
+       _LOOP_BEGIN_
+          _GET_(self%id_n,O5n)
+          _GET_(self%id_ETW,ETW)
+          _GET_(self%id_X1X,X1X)
+          _GET_HORIZONTAL_(self%id_pN2Oa,pN2Oa)
+
+         koN2O = n2o_transfer_coefficient(self,ETW,X1X)
+
+       _SET_DIAGNOSTIC_(self%id_satp,(100._rk*O5n*0.5_rk/(KoN2O*pN2Oa*1.e-9*1.e6)))
+      _LOOP_END_
+   end subroutine
+
+   subroutine do_surface(self,_ARGUMENTS_DO_SURFACE_)
+      class (type_ersem_nitrous_oxide), intent(in) :: self
+      _DECLARE_ARGUMENTS_DO_SURFACE_
+
+      real(rk) :: O5n,ETW,X1X,wnd, fwind
+      real(rk) :: koN2O,sc,n2oflux,pN2Oa
+
+      _HORIZONTAL_LOOP_BEGIN_
+         _GET_(self%id_n,O5n)
+         _GET_(self%id_ETW,ETW)
+         _GET_(self%id_X1X,X1X)
+         _GET_HORIZONTAL_(self%id_wnd,wnd)
+         _GET_HORIZONTAL_(self%id_pN2Oa,pN2Oa)
+
+! Schmidt number for N2O (Wanninkhof, 1992)
+
+        sc=2301._rk-151._rk*ETW+4.74_rk*ETW**2._rk-0.06_rk*ETW**3
+        fwind =  0.39_rk * wnd**2 *(sc/660._rk)**(-0.5_rk)
+        fwind=fwind*24._rk/100._rk   ! convert to m/day
+
+        koN2O = n2o_transfer_coefficient(self,ETW,X1X)
+
+!        PNOatm is given in natm, so it is here converted to atm  by means of a 1.-9 factor. As KoNO is
+!        given in mol L-1 atm -1 (Weiss and Price, 1980), a conversion factor of
+!        1.e6 is used to have the final units of mmol m-3. N0 is given in mmol
+!        of N so it is multiplied by 0.5 to have mmol of N2O. (Luca, July 2016)
+
+        if (self%iswN2O .eq. 1) then
+         n2oflux = fwind*(KoN2O*pN2Oa*1.e-9*1.e6-O5n*0.5_rk)
+        else
+         n2oflux=0._rk
+        endif
+
+         _SET_HORIZONTAL_DIAGNOSTIC_(self%id_airsea,n2oflux)
+
+         _SET_SURFACE_EXCHANGE_(self%id_n,n2oflux*2._rk)
+
+      _HORIZONTAL_LOOP_END_
+   end subroutine do_surface
+
+   function n2o_transfer_coefficient(self,ETW,X1X) result(koN2O)
+       class (type_ersem_nitrous_oxide), intent(in) :: self
+       real(rk),                         intent(in) :: ETW,X1X
+       real(rk)                                     :: koN2O
+       real(rk)           :: tk,tk100
+       real(rk),parameter :: A1 = -62.7062_rk
+       real(rk),parameter :: A2 = 97.3066_rk
+       real(rk),parameter :: A3 = 24.1406_rk
+       real(rk),parameter :: B1 = -0.05842_rk
+       real(rk),parameter :: B2 = 0.033193_rk
+       real(rk),parameter :: B3 = -0.0051_rk
+
+       TK=ETW+273.15_rk
+       TK100=TK/100._rk
+
+        koN2O = exp(A1+A2/tk100 + A3 * log(tk100) + &
+        &       X1X * (B1 + B2 * tk100 + B3 * tk100 ** 2._rk))
+    end function
+
+end module

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
@@ -142,6 +142,7 @@ subroutine initialize(self,configunit)
              call self%register_diagnostic_variable(self%id_kflux_f,'kflux_f','nmolFe/m^3/d','1st order kinetic flux of iron',source=source_do)
           endif
       endif
+      if (index(composition,'h')/=0) call self%add_constituent('h',0.0_rk)
 
    end subroutine
 
@@ -207,6 +208,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')
 
@@ -44,6 +45,7 @@ module ersem_shared
    type (type_horizontal_standard_variable),parameter :: depth_of_bottom_interface_of_layer_1 = type_horizontal_standard_variable(name='depth_of_bottom_interface_of_layer_1',units='m')
    type (type_horizontal_standard_variable),parameter :: depth_of_bottom_interface_of_layer_2 = type_horizontal_standard_variable(name='depth_of_bottom_interface_of_layer_2',units='m')
    type (type_horizontal_standard_variable),parameter :: pelagic_benthic_transfer_constant = type_horizontal_standard_variable(name='pelagic_benthic_transfer_constant',units='d/m')
+   type (type_horizontal_standard_variable),parameter :: partial_pressure_of_n2o = type_horizontal_standard_variable(name='partial_pressure_of_n2o',units='natm')
    type (type_horizontal_standard_variable),parameter :: sediment_erosion = type_horizontal_standard_variable(name='sediment_erosion',units='m/d')
 
    ! Aggregate absorption and backscatter.