added akvaplan/antibiotic

src/models/akvaplan/CMakeLists.txt

@@ -3,6 +3,7 @@ add_library(fabm_models_akvaplan OBJECT
             tracer.F90
             plume_injection.F90
             tracer_sed.F90
+            antiparasitic.F90
            )
 
 add_dependencies(fabm_models_akvaplan fabm_base)

src/models/akvaplan/akvaplan_model_library.F90

@@ -8,6 +8,7 @@ module akvaplan_model_library
    use akvaplan_tracer
    use akvaplan_plume_injection
    use akvaplan_tracer_sed
+   use akvaplan_antiparasitic
 
    implicit none
 
@@ -31,6 +32,7 @@ subroutine create(self,name,model)
          case ('tracer');          allocate(type_tracer::model)
          case ('plume_injection'); allocate(type_plume_injection::model)
          case ('tracer_sed');      allocate(type_tracer_sed::model)
+         case ('antiparasitic');   allocate(type_antiparasitic::model)
          ! Add new models here
          case default
             call self%type_base_model_factory%create(name,model)

src/models/akvaplan/antiparasitic.F90

@@ -0,0 +1,285 @@
+#include "fabm_driver.h"
+
+module akvaplan_antiparasitic
+
+   ! Model for a dissolved compound ("antiparasitic") that partially absorbs onto particulate organic matter (POM).
+   ! Partitioning between dissolved and adsorbed phases is assumed to be instantaneous and
+   ! controlled by the organic carbon : water partition coefficient.
+   ! The compound experiences temperature-dependent degradation at phase-dependent rates.
+   ! The model can contain several classes of POM with different sinking velocities.
+   ! At the bed, both the compound and POM are incorporated into the sediment; from where they may be resuspended again.
+   ! Implementation by Jorn Bruggeman and Ricardo Torres, Plymouth Marine Laboratory
+
+   use fabm_types
+
+   implicit none
+
+   private
+
+   type, extends(type_base_model), public :: type_antiparasitic
+      ! The model's own variables
+      type (type_state_variable_id),        allocatable :: id_pom(:)
+      type (type_bottom_state_variable_id), allocatable :: id_pom_bot(:)
+      type (type_state_variable_id)                     :: id_antiparasitic
+      type (type_bottom_state_variable_id)              :: id_antiparasitic_bot
+
+      ! Environmental dependencies
+      type (type_dependency_id)            :: id_T     ! temperature
+      type (type_dependency_id)            :: id_h     ! cell thickness
+      type (type_horizontal_dependency_id) :: id_shear ! bottom shear
+
+      ! Identifiers for diagnostic variables
+      type (type_diagnostic_variable_id)            :: id_degradation_pom_flux, id_degradation_c_flux          ! pelagic degradation flux for pom and treatment (antiparasitic!)
+      type (type_horizontal_diagnostic_variable_id) :: id_resuspension_flux_pom, id_deposition_flux_pom        ! deposition and resuspension fluxes of all pom 
+      type (type_horizontal_diagnostic_variable_id) :: id_resuspension_flux_c, id_deposition_flux_c            ! deposition and resuspension fluxes of antiparasitic treatment 
+      type (type_horizontal_diagnostic_variable_id) :: id_degradation_pom_bot_flux, id_degradation_c_bot_flux  ! sediment degradation fluxes for pom and treatment (antiparasitic!)
+
+      ! Parameters
+      real(rk), allocatable :: w(:)
+      real(rk), allocatable :: tau_bot_crit(:)
+      real(rk) :: bed_por
+      real(rk) :: erate
+      real(rk) :: K_oc
+      real(rk) :: E_a
+      real(rk) :: T_ref
+      real(rk) :: k_w
+      real(rk) :: k_ads
+      real(rk) :: k_pom
+      real(rk) :: max_dt    ! maximum time step (s) to be expected - used to limit benthic-pelagic fluxes for numerical stability
+   contains
+      ! Model procedures
+      procedure :: initialize
+      procedure :: do
+      procedure :: do_bottom
+      procedure :: get_vertical_movement
+   end type
+
+   real(rk), parameter :: secs_pr_day = 86400.0_rk
+   real(rk), parameter :: R = 8.3144598_rk    ! universal gas constant (J mol-1 K-1)
+   real(rk), parameter :: Kelvin = 273.15_rk  ! offset between degrees Celsius and Kelvin
+
+contains
+
+   subroutine initialize(self, configunit)
+      class (type_antiparasitic), intent(inout), target :: self
+      integer,                    intent(in)            :: configunit
+
+      integer                         :: i, npom
+      character(len=8)                :: index
+      logical                         :: conserved
+      character(len=attribute_length) :: standard_name
+
+      ! Obtain parameter values (convert all rate constants to s-1 to match FABM's internal time unit)
+      call self%get_parameter(npom, 'npom', '', 'number of particulate organic matter classes')
+      allocate(self%id_pom(npom), self%id_pom_bot(npom), self%w(npom), self%tau_bot_crit(npom))
+      do i=1,npom
+         write (index,'(i0)') i
+         call self%get_parameter(self%w(i), 'w'//trim(index), 'm d-1', 'sinking velocity of particulate organic matter class '//trim(index), scale_factor=1.0_rk/secs_pr_day)
+         call self%register_state_variable(self%id_pom(i), 'pom'//trim(index), 'g C m-3', 'concentration of particulate organic matter class '//trim(index), initial_value=0.0_rk, vertical_movement=-self%w(i), minimum=0.0_rk)
+         call self%register_state_variable(self%id_pom_bot(i), 'pom'//trim(index)//'_bot', 'g C m-2', 'particulate organic matter class '//trim(index)//' in sediment', initial_value=0.0_rk, minimum=0.0_rk)
+      end do
+      call self%get_parameter(self%K_oc, 'K_oc', 'm3 g-1', 'organic carbon : water partition coefficient')
+
+      ! Degradation of antiparasitic and POM
+      call self%get_parameter(self%E_a, 'E_a', 'J mol-1', 'activation energy for antiparasitic decay (0: no temperature dependence)', minimum=0.0_rk, default=0.0_rk)
+      call self%get_parameter(self%T_ref, 'T_ref', 'degree Celsius', 'reference temperature', default=20.0_rk)
+      call self%get_parameter(self%k_ads, 'k_ads', 'd-1', 'degradation of adsorbed antiparasitic at reference temperature', default=0.0_rk, scale_factor=1.0_rk/secs_pr_day, minimum=0._rk)
+      call self%get_parameter(self%k_w, 'k_w', 'd-1', 'degradation of dissolved antiparasitic at reference temperature', default=0.0_rk, scale_factor=1.0_rk/secs_pr_day, minimum=0._rk)
+      call self%get_parameter(self%k_pom, 'k_pom', 'd-1', 'degradation of particulate organic matter', default=0.0_rk, scale_factor=1.0_rk/secs_pr_day, minimum=0._rk)
+
+      ! Resuspension of POM
+      call self%get_parameter(self%bed_por, 'bed_por', '-', 'bed porosity')
+      call self%get_parameter(self%erate, 'erate', 'g m-2 d-1', 'bed erodibility', default=0.0_rk, scale_factor=1.0_rk/secs_pr_day, minimum=0._rk)
+      do i=1,npom
+         write (index,'(i0)') i
+         call self%get_parameter(self%tau_bot_crit(i), 'tau_bot_crit'//trim(index), 'Pa', 'critical shear stress for particulate organic matter class  '//trim(index), default=0.0_rk, scale_factor=1.0_rk/secs_pr_day, minimum=0.0_rk)
+      end do
+
+      ! Numerics and debugging
+      call self%get_parameter(self%max_dt, 'max_dt', 's', 'maximum time step used to limit vertical exchanges', default=900.0_rk)
+      call self%get_parameter(conserved, 'conserved', '', 'treat antiparasitic and POM as conserved quantities (activates budget tracking in host)', default=.false.)
+
+      ! Register state variables for the antiparasitic compound.
+      call self%register_state_variable(self%id_antiparasitic,'c','g m-3','total antiparasitic concentration (dissolved + adbsorbed)', initial_value=0.0_rk, minimum=0.0_rk)
+      call self%register_state_variable(self%id_antiparasitic_bot,'c_bot','g m-2','total antiparasitic in sediment', minimum=0.0_rk)
+
+      if (conserved) then
+         standard_name = get_safe_name(trim(self%get_path())//'_total_ab')
+         call self%add_to_aggregate_variable(type_bulk_standard_variable(name=standard_name(2:), units='g m-3', aggregate_variable=.true.,conserved=.true.), self%id_antiparasitic)
+         call self%add_to_aggregate_variable(type_bulk_standard_variable(name=standard_name(2:), units='g m-3', aggregate_variable=.true.,conserved=.true.), self%id_antiparasitic_bot)
+
+         standard_name = get_safe_name(trim(self%get_path())//'_total_pom')
+         do i=1,npom
+           call self%add_to_aggregate_variable(type_bulk_standard_variable(name=standard_name(2:), units='g m-3', aggregate_variable=.true.,conserved=.true.), self%id_pom(i))
+           call self%add_to_aggregate_variable(type_bulk_standard_variable(name=standard_name(2:), units='g m-3', aggregate_variable=.true.,conserved=.true.), self%id_pom_bot(i))
+         end do
+     end if
+
+      ! Register environmental dependencies.
+      call self%register_dependency(self%id_T, standard_variables%temperature)
+      call self%register_dependency(self%id_shear, standard_variables%bottom_stress)
+      call self%register_dependency(self%id_h, standard_variables%cell_thickness)
+
+      ! Register diagnostic variables
+      call self%register_diagnostic_variable(self%id_degradation_pom_flux, 'pom_degradation', 'g C m-3 s-1', 'degradation of total particulate organic matter in pelagic')
+      call self%register_diagnostic_variable(self%id_degradation_c_flux, 'c_degradation', 'g m-3 s-1', 'degradation of total antiparasitic in pelagic')
+
+      call self%register_diagnostic_variable(self%id_resuspension_flux_pom, 'resuspension_flux_pom', 'g C m-2 s-1', 'total organic matter resuspension', source=source_do_bottom)
+      call self%register_diagnostic_variable(self%id_deposition_flux_pom, 'deposition_flux_pom', 'g C m-2 s-1', 'total organic matter deposition', source=source_do_bottom)
+      call self%register_diagnostic_variable(self%id_resuspension_flux_c, 'resuspension_flux_c', 'g m-2 s-1', 'antiparasitic resuspension', source=source_do_bottom)
+      call self%register_diagnostic_variable(self%id_deposition_flux_c, 'deposition_flux_c', 'g m-2 s-1', 'antiparasitic deposition', source=source_do_bottom)
+
+      call self%register_diagnostic_variable(self%id_degradation_pom_bot_flux, 'pom_bot_degradation','g C m-2 s-1', 'degradation of total organic matter in sediment', source=source_do_bottom)
+      call self%register_diagnostic_variable(self%id_degradation_c_bot_flux, 'c_bot_degradation','g m-2 s-1', 'degradation of total antiparasitic in sediment', source=source_do_bottom)
+
+   end subroutine initialize
+
+   subroutine do(self, _ARGUMENTS_DO_)
+      class (type_antiparasitic), intent(in) :: self
+      _DECLARE_ARGUMENTS_DO_
+
+      integer :: i
+      real(rk) :: pom(size(self%id_pom)), pom_tot, frac_ads, frac_wat, T, f, c, k_wat, k_ads
+
+      _LOOP_BEGIN_
+
+         do i=1,size(self%id_pom)
+            _GET_(self%id_pom(i), pom(i))
+         end do
+         pom_tot = sum(pom)
+         frac_wat = 1.0_rk / (self%K_oc * pom_tot + 1.0_rk)
+         frac_ads = 1.0_rk - frac_wat
+
+         ! Specific degradation (s-1) for adsorbed and dissolved antiparasitic
+         _GET_(self%id_T, T)
+         f = exp(self%E_a / R * (1.0_rk / (Kelvin + self%T_ref) - 1.0_rk / (Kelvin + T)))
+         k_ads = self%k_ads * f
+         k_wat = self%k_w * f
+         _GET_(self%id_antiparasitic, c)
+         _SET_ODE_(self%id_antiparasitic, -(frac_ads * k_ads + frac_wat * k_wat) * c)
+
+         ! Degradation (s-1) for POM
+         do i=1,size(self%id_pom)
+            _SET_ODE_(self%id_pom(i), -self%k_pom * pom(i))
+         end do
+
+         ! Save diagnostics for total pom and antiparasitic degradation  
+         _SET_DIAGNOSTIC_(self%id_degradation_pom_flux, self%k_pom * pom_tot) 
+         _SET_DIAGNOSTIC_(self%id_degradation_c_flux, (frac_ads * k_ads + frac_wat * k_wat) * c) 
+
+      _LOOP_END_
+
+   end subroutine do
+
+   subroutine do_bottom(self, _ARGUMENTS_DO_BOTTOM_)
+      class (type_antiparasitic), intent(in) :: self
+      _DECLARE_ARGUMENTS_DO_BOTTOM_
+
+      integer :: i
+      real(rk) :: h, pom(size(self%id_pom)), pom_bot(size(self%id_pom)), pom_tot, pom_bot_tot, frac_ads, w, c, tau_bot, erosion(size(self%id_pom)), rel_antiparasitic_erosion, T, f, k_ads, c_bot
+      real(rk) :: pom_deposition
+
+      _HORIZONTAL_LOOP_BEGIN_
+
+         _GET_(self%id_h, h)
+         do i = 1, size(self%id_pom)
+            _GET_(self%id_pom(i), pom(i))
+            _GET_HORIZONTAL_(self%id_pom_bot(i), pom_bot(i))
+         end do
+         pom_tot = sum(pom)
+         pom_bot_tot = sum(pom_bot)
+
+         ! Erosion rate (g m-2 s-1) per POM class
+         _GET_HORIZONTAL_(self%id_shear, tau_bot)
+         erosion = min(pom_bot/self%max_dt, self%erate * (1.0_rk - self%bed_por) * max(tau_bot / self%tau_bot_crit - 1.0_rk, 0.0_rk))
+
+         ! Calculate relative erosion for antiparasitic substance.
+         if (pom_bot_tot > 0) then
+            rel_antiparasitic_erosion = sum(erosion)/pom_bot_tot
+         else
+            rel_antiparasitic_erosion = 0
+         end if
+
+         ! Sedimentation, resuspension, degradation per POM class
+         do i = 1, size(self%id_pom)
+           ! From pelagic to bottom. Convention is negative out of pelagic, positive into water
+            _SET_BOTTOM_EXCHANGE_(self%id_pom(i), -min(h/self%max_dt, self%w(i)) * pom(i) + erosion(i))
+           ! From bottom to pelagic. Convention is positive into bottom, negative into water. Added degradation of pom_bot
+            _SET_BOTTOM_ODE_(self%id_pom_bot(i), min(h/self%max_dt, self%w(i)) * pom(i) - erosion(i) - self%k_pom * pom_bot(i))
+         end do
+
+         ! Sinking rate of antiparasitic is weighted average of sinking rate per antiparasitic fraction
+         ! (those adsorbed to individual POM classes, plus one non-sinking fraction for the dissolved phase)
+         frac_ads = 1.0_rk - 1.0_rk/(self%K_oc * pom_tot + 1.0_rk)
+         if (pom_tot > 0) then
+            w = min(h/self%max_dt, frac_ads * sum(pom * self%w) / pom_tot)
+         else
+            w = 0
+         end if
+
+         ! Specific degradation (s-1) for adsorbed antiparasitic in sediment
+         ! (assume dissolved fraction in pore water is negligible)
+         _GET_(self%id_T, T)
+         f = exp(self%E_a / R * (1.0_rk / (Kelvin + self%T_ref) - 1.0_rk / (Kelvin + T)))
+         k_ads = self%k_ads * f
+
+         ! Sedimentation and resuspension for antiparasitic
+         _GET_(self%id_antiparasitic, c)
+         _GET_HORIZONTAL_(self%id_antiparasitic_bot, c_bot)
+         _SET_BOTTOM_EXCHANGE_(self%id_antiparasitic, -w * c + rel_antiparasitic_erosion * c_bot)
+         _SET_BOTTOM_ODE_(self%id_antiparasitic_bot, w * c - rel_antiparasitic_erosion * c_bot - k_ads * c_bot)
+
+         ! Save diagnostic total pom bottom degradation  
+         _SET_HORIZONTAL_DIAGNOSTIC_(self%id_degradation_pom_bot_flux, self%k_pom * pom_bot_tot) 
+         ! Save diagnostic antiparasitic bottom degradation  
+         _SET_HORIZONTAL_DIAGNOSTIC_(self%id_degradation_c_bot_flux, k_ads * c_bot ) 
+
+         ! Save diagnostic total pom bottom deposition
+         pom_deposition = 0.0_rk
+         do i = 1, size(self%id_pom)
+            pom_deposition = pom_deposition + min(h/self%max_dt, self%w(i)) * pom(i) 
+         end do
+         _SET_HORIZONTAL_DIAGNOSTIC_(self%id_deposition_flux_pom, pom_deposition)
+
+         ! Save diagnostic total pom bottom resuspension
+         _SET_HORIZONTAL_DIAGNOSTIC_(self%id_resuspension_flux_pom, sum(erosion))
+
+         ! Save diagnostic antiparasitic bottom deposition  
+         _SET_HORIZONTAL_DIAGNOSTIC_(self%id_deposition_flux_c, w * c )
+
+         ! Save diagnostic antiparasitic bottom resuspension  
+         _SET_HORIZONTAL_DIAGNOSTIC_(self%id_resuspension_flux_c, rel_antiparasitic_erosion * c_bot) 
+
+       _HORIZONTAL_LOOP_END_
+
+   end subroutine do_bottom
+
+   subroutine get_vertical_movement(self, _ARGUMENTS_GET_VERTICAL_MOVEMENT_)
+      class (type_antiparasitic), intent(in) :: self
+      _DECLARE_ARGUMENTS_GET_VERTICAL_MOVEMENT_
+
+      integer :: i
+      real(rk) :: pom(size(self%id_pom)), pom_tot, frac_ads, w
+
+      _LOOP_BEGIN_
+
+         ! Sinking rate of antiparasitic is weighted average of sinking rate per antiparasitic fraction
+         ! (those adsorbed to individual POM classes, plus one non-sinking fraction for the dissolved phase)
+         ! NB we do not need to set sinking rate per POM class - that is handled by the vertical_movement argument specified during initialization
+         do i = 1, size(self%id_pom)
+            _GET_(self%id_pom(i), pom(i))
+         end do
+         pom_tot = sum(pom)
+         frac_ads = 1.0_rk - 1.0_rk/(self%K_oc * pom_tot + 1.0_rk)
+         if (pom_tot > 0) then
+            w = frac_ads * sum(pom * self%w) / pom_tot
+         else
+            w = 0
+         end if
+         _SET_VERTICAL_MOVEMENT_(self%id_antiparasitic, -w)
+
+      _LOOP_END_
+
+   end subroutine get_vertical_movement
+
+end module

testcases/fabm-akvaplan-antiparasitic.yaml

@@ -0,0 +1,49 @@
+instances:
+   ab:
+      model: akvaplan/antiparasitic
+      parameters:
+         npom: 8
+         w1: 216.0    # sinking velocity of POM class 1 (m d-1)
+         w2: 648.0    # sinking velocity of POM class 2 (m d-1)
+         w3: 1080.0   # sinking velocity of POM class 3 (m d-1)
+         w4: 1728.0   # sinking velocity of POM class 4 (m d-1)
+         w5: 3240.0   # sinking velocity of POM class 5 (m d-1)
+         w6: 6480.0   # sinking velocity of POM class 6 (m d-1)
+         w7: 7603.2   # sinking velocity of POM class 7 (m d-1)
+         w8: 10368.0  # sinking velocity of POM class 8 (m d-1)
+         K_oc: 0.005  # water : POM partition coefficient of antiparasitic (m3 g-1)
+         E_a: 0       # activation energy for antiparasitic degradation (J mol-1)
+         T_ref: 20    # reference temperature for antiparasitic degradation (degrees Celsius)
+         k_w: 0.0     # degradation rate of dissolved antiparasitic at reference temperature (d-1)
+         k_ads: 0.0     # degradation rate of adsorbed antiparasitic at reference temperature (d-1)
+         k_pom: 0.0     # degradation rate of POM (d-1)
+         erate: 0.0035  # maximum erosion rate (g m-2 d-1) - should this vary per POM class???
+         bed_por: 0.4   # sediment porosity
+         tau_bot_crit1: 0.005 # critical bottom shear stress for resuspension of POM class 1 (Pa)
+         tau_bot_crit2: 0.005 # critical bottom shear stress for resuspension of POM class 2 (Pa)
+         tau_bot_crit3: 0.005 # critical bottom shear stress for resuspension of POM class 3 (Pa)
+         tau_bot_crit4: 0.005 # critical bottom shear stress for resuspension of POM class 4 (Pa)
+         tau_bot_crit5: 0.005 # critical bottom shear stress for resuspension of POM class 5 (Pa)
+         tau_bot_crit6: 0.005 # critical bottom shear stress for resuspension of POM class 6 (Pa)
+         tau_bot_crit7: 0.005 # critical bottom shear stress for resuspension of POM class 7 (Pa)
+         tau_bot_crit8: 0.005 # critical bottom shear stress for resuspension of POM class 8 (Pa)
+      initialization:
+         pom1: 1 # pelagic concentration of POM class 1 (g m-3)
+         pom2: 1 # pelagic concentration of POM class 2 (g m-3)
+         pom3: 1 # pelagic concentration of POM class 3 (g m-3)
+         pom4: 1 # pelagic concentration of POM class 4 (g m-3)
+         pom5: 1 # pelagic concentration of POM class 5 (g m-3)
+         pom6: 1 # pelagic concentration of POM class 6 (g m-3)
+         pom7: 1 # pelagic concentration of POM class 7 (g m-3)
+         pom8: 1 # pelagic concentration of POM class 8 (g m-3)
+         pom1_bot: 0 # density of POM class 1 in sediment (g m-2)
+         pom2_bot: 0 # density of POM class 2 in sediment (g m-2)
+         pom3_bot: 0 # density of POM class 3 in sediment (g m-2)
+         pom4_bot: 0 # density of POM class 4 in sediment (g m-2)
+         pom5_bot: 0 # density of POM class 5 in sediment (g m-2)
+         pom6_bot: 0 # density of POM class 6 in sediment (g m-2)
+         pom7_bot: 0 # density of POM class 7 in sediment (g m-2)
+         pom8_bot: 0 # density of POM class 8 in sediment (g m-2)
+         c: 1 # pelagic concentration of antiparasitic (g m-3)
+         c_bot: 0 # density of antiparasitic in sediment (g m-2)
+