From 69da88c87ec59a51d1e2143c1f76111526ed6498 Mon Sep 17 00:00:00 2001 From: Jorn Bruggeman Date: Tue, 4 Feb 2020 21:17:44 +0000 Subject: [PATCH] added akvaplan/antibiotic --- src/models/akvaplan/CMakeLists.txt | 1 + .../akvaplan/akvaplan_model_library.F90 | 2 + src/models/akvaplan/antiparasitic.F90 | 285 ++++++++++++++++++ testcases/fabm-akvaplan-antiparasitic.yaml | 49 +++ 4 files changed, 337 insertions(+) create mode 100644 src/models/akvaplan/antiparasitic.F90 create mode 100644 testcases/fabm-akvaplan-antiparasitic.yaml diff --git a/src/models/akvaplan/CMakeLists.txt b/src/models/akvaplan/CMakeLists.txt index c0a50929..e395ca8b 100644 --- a/src/models/akvaplan/CMakeLists.txt +++ b/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) diff --git a/src/models/akvaplan/akvaplan_model_library.F90 b/src/models/akvaplan/akvaplan_model_library.F90 index a588fd9c..aa4339c7 100644 --- a/src/models/akvaplan/akvaplan_model_library.F90 +++ b/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) diff --git a/src/models/akvaplan/antiparasitic.F90 b/src/models/akvaplan/antiparasitic.F90 new file mode 100644 index 00000000..fdfe7c97 --- /dev/null +++ b/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 diff --git a/testcases/fabm-akvaplan-antiparasitic.yaml b/testcases/fabm-akvaplan-antiparasitic.yaml new file mode 100644 index 00000000..e46ece96 --- /dev/null +++ b/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) +