From dfae26118196b43cc42b50accf6b5119058b170e Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 14 Oct 2024 15:52:49 +0100 Subject: [PATCH 1/8] removed old file --- .../PISCES/dissolved_organic_matter.jl | 137 ------------------ 1 file changed, 137 deletions(-) delete mode 100644 src/Models/AdvectedPopulations/PISCES/dissolved_organic_matter.jl diff --git a/src/Models/AdvectedPopulations/PISCES/dissolved_organic_matter.jl b/src/Models/AdvectedPopulations/PISCES/dissolved_organic_matter.jl deleted file mode 100644 index 7c3b13c14..000000000 --- a/src/Models/AdvectedPopulations/PISCES/dissolved_organic_matter.jl +++ /dev/null @@ -1,137 +0,0 @@ -""" - DissolvedOrganicMatter - -Parameterisation of dissolved organic matter which depends on a bacterial -concentration derived from the concentration of zooplankton. -""" -@kwdef struct DissolvedOrganicMatter{FT, AP} - remineralisation_rate :: FT = 0.3/day # 1 / s - microzooplankton_bacteria_concentration :: FT = 0.7 # - mesozooplankton_bacteria_concentration :: FT = 1.4 # - maximum_bacteria_concentration :: FT = 4.0 # mmol C / m³ - bacteria_concentration_depth_exponent :: FT = 0.684 # - reference_bacteria_concentration :: FT = 1.0 # mmol C / m³ - temperature_sensetivity :: FT = 1.066 # - doc_half_saturation_for_bacterial_activity :: FT = 417.0 # mmol C / m³ - nitrate_half_saturation_for_bacterial_activity :: FT = 0.03 # mmol N / m³ - ammonia_half_saturation_for_bacterial_activity :: FT = 0.003 # mmol N / m³ - phosphate_half_saturation_for_bacterial_activity :: FT = 0.003 # mmol P / m³ - iron_half_saturation_for_bacterial_activity :: FT = 0.01 # μmol Fe / m³ -# (1 / (mmol C / m³), 1 / (mmol C / m³), 1 / (mmol C / m³), 1 / (mmol C / m³) / s, 1 / (mmol C / m³) / s) - aggregation_parameters :: AP = (0.37, 102, 3530, 5095, 114) .* (10^-6 / day) - maximum_iron_ratio_in_bacteria :: FT = 0.06 # μmol Fe / mmol C - iron_half_saturation_for_bacteria :: FT = 0.3 # μmol Fe / m³ - maximum_bacterial_growth_rate :: FT = 0.6 / day # 1 / s -end - -@inline function (dom::DissolvedOrganicMatter)(::Val{:DOC}, bgc, - x, y, z, t, - P, D, Z, M, - PChl, DChl, PFe, DFe, DSi, - DOC, POC, GOC, - SFe, BFe, PSi, - NO₃, NH₄, PO₄, Fe, Si, - CaCO₃, DIC, Alk, - O₂, T, S, - zₘₓₗ, zₑᵤ, Si′, Ω, κ, mixed_layer_PAR, wPOC, wGOC, PAR, PAR₁, PAR₂, PAR₃) - - nanophytoplankton_exudation = dissolved_exudate(bgc.nanophytoplankton, bgc, y, t, P, PChl, PFe, NO₃, NH₄, PO₄, Fe, Si, T, Si′, zₘₓₗ, zₑᵤ, κ, PAR₁, PAR₂, PAR₃) - diatom_exudation = dissolved_exudate(bgc.diatoms, bgc, y, t, D, DChl, DFe, NO₃, NH₄, PO₄, Fe, Si, T, Si′, zₘₓₗ, zₑᵤ, κ, PAR₁, PAR₂, PAR₃) - - phytoplankton_exudation = nanophytoplankton_exudation + diatom_exudation - - particulate_degredation = specific_degredation_rate(bgc.particulate_organic_matter, bgc, O₂, T) * POC - - respiration_product = dissolved_upper_trophic_respiration_product(bgc.mesozooplankton, M, T) - - microzooplankton_grazing_waste = specific_dissolved_grazing_waste(bgc.microzooplankton, P, D, PFe, DFe, Z, POC, GOC, SFe, T, wPOC, wGOC) * Z - mesozooplankton_grazing_waste = specific_dissolved_grazing_waste(bgc.mesozooplankton, P, D, PFe, DFe, Z, POC, GOC, SFe, T, wPOC, wGOC) * M - - grazing_waste = microzooplankton_grazing_waste + mesozooplankton_grazing_waste - - degredation = bacterial_degradation(dom, z, Z, M, DOC, NO₃, NH₄, PO₄, Fe, T, zₘₓₗ, zₑᵤ) - - aggregation_to_particles, = aggregation(dom, bgc, z, DOC, POC, GOC, zₘₓₗ) - - return phytoplankton_exudation + particulate_degredation + respiration_product + grazing_waste - degredation - aggregation_to_particles -end - -@inline function bacteria_concentration(dom::DissolvedOrganicMatter, z, Z, M, zₘₓₗ, zₑᵤ) - bZ = dom.microzooplankton_bacteria_concentration - bM = dom.mesozooplankton_bacteria_concentration - a = dom.bacteria_concentration_depth_exponent - - zₘ = min(zₘₓₗ, zₑᵤ) - - surface_bacteria = min(4, bZ * Z + bM * M) - - depth_factor = (zₘ / z) ^ a - - return ifelse(z >= zₘ, 1, depth_factor) * surface_bacteria -end - -@inline function bacteria_activity(dom::DissolvedOrganicMatter, DOC, NO₃, NH₄, PO₄, Fe) - K_DOC = dom.doc_half_saturation_for_bacterial_activity - K_NO₃ = dom.nitrate_half_saturation_for_bacterial_activity - K_NH₄ = dom.ammonia_half_saturation_for_bacterial_activity - K_PO₄ = dom.phosphate_half_saturation_for_bacterial_activity - K_Fe = dom.iron_half_saturation_for_bacterial_activity - - DOC_limit = DOC / (DOC + K_DOC) - - L_N = (K_NO₃ * NH₄ + K_NH₄ * NO₃) / (K_NO₃ * K_NH₄ + K_NO₃ * NH₄ + K_NH₄ * NO₃) - - L_PO₄ = PO₄ / (PO₄ + K_PO₄) - - L_Fe = Fe / (Fe + K_Fe) - - # assuming typo in paper otherwise it doesn't make sense to formulate L_NH₄ like this - limiting_quota = min(L_N, L_PO₄, L_Fe) - - return limiting_quota * DOC_limit -end - -@inline function bacterial_degradation(dom::DissolvedOrganicMatter, z, Z, M, DOC, NO₃, NH₄, PO₄, Fe, T, zₘₓₗ, zₑᵤ) - Bact_ref = dom.reference_bacteria_concentration - b = dom.temperature_sensetivity - λ = dom.remineralisation_rate - - f = b^T - - Bact = bacteria_concentration(dom, z, Z, M, zₘₓₗ, zₑᵤ) - - LBact = bacteria_activity(dom, DOC, NO₃, NH₄, PO₄, Fe) - - return λ * f * LBact * Bact / Bact_ref * DOC # differes from Aumont 2015 since the dimensions don't make sense -end - -@inline function oxic_remineralisation(dom::DissolvedOrganicMatter, bgc, z, Z, M, DOC, NO₃, NH₄, PO₄, Fe, O₂, T, zₘₓₗ, zₑᵤ) - ΔO₂ = anoxia_factor(bgc, O₂) - - degredation = bacterial_degradation(dom, z, Z, M, DOC, NO₃, NH₄, PO₄, Fe, T, zₘₓₗ, zₑᵤ) - - return (1 - ΔO₂) * degredation -end - -@inline function denitrifcation(dom::DissolvedOrganicMatter, bgc, z, Z, M, DOC, NO₃, NH₄, PO₄, Fe, O₂, T, zₘₓₗ, zₑᵤ) - ΔO₂ = anoxia_factor(bgc, O₂) - - degredation = bacterial_degradation(dom, z, Z, M, DOC, NO₃, NH₄, PO₄, Fe, T, zₘₓₗ, zₑᵤ) - - return ΔO₂ * degredation -end - -@inline function aggregation(dom::DissolvedOrganicMatter, bgc, z, DOC, POC, GOC, zₘₓₗ) - a₁, a₂, a₃, a₄, a₅ = dom.aggregation_parameters - - backgroound_shear = bgc.background_shear - mixed_layer_shear = bgc.mixed_layer_shear - - shear = ifelse(z < zₘₓₗ, backgroound_shear, mixed_layer_shear) - - Φ₁ = shear * (a₁ * DOC + a₂ * POC) * DOC - Φ₂ = shear * (a₃ * GOC) * DOC - Φ₃ = (a₄ * POC + a₅ * DOC) * DOC - - return Φ₁ + Φ₂ + Φ₃, Φ₁, Φ₂, Φ₃ -end From df4a3a6d0f7276b4070ea6a12585510d0a898144 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 14 Oct 2024 17:45:16 +0100 Subject: [PATCH 2/8] attempted to automatically set float types in BGC models, carbon chem, gas exchange, and particles core --- src/Light/2band.jl | 4 +- src/Light/Light.jl | 2 +- src/Light/compute_euphotic_depth.jl | 4 +- src/Light/multi_band.jl | 12 +- .../AdvectedPopulations/LOBSTER/LOBSTER.jl | 9 +- src/Models/AdvectedPopulations/NPZD.jl | 11 +- .../AdvectedPopulations/PISCES/PISCES.jl | 123 ++--- .../dissolved_organic_carbon.jl | 31 +- .../two_size_class.jl | 79 ++- .../PISCES/phytoplankton/mixed_mondo.jl | 93 +++- .../phytoplankton/mixed_mondo_nano_diatoms.jl | 49 +- .../phytoplankton/nutrient_limitation.jl | 4 +- .../PISCES/zooplankton/defaults.jl | 38 +- .../zooplankton/food_quality_dependant.jl | 70 ++- .../CarbonChemistry/carbon_chemistry.jl | 68 ++- .../CarbonChemistry/equilibrium_constants.jl | 509 +++++++++++------- src/Models/GasExchange/GasExchange.jl | 71 +-- src/Models/GasExchange/gas_solubility.jl | 16 +- .../GasExchange/gas_transfer_velocity.jl | 50 +- .../GasExchange/generic_parameterisations.jl | 4 +- src/Models/GasExchange/schmidt_number.jl | 12 +- src/Models/GasExchange/surface_values.jl | 2 + src/Models/Individuals/SugarKelp/SugarKelp.jl | 3 + src/Particles/Particles.jl | 13 +- 24 files changed, 801 insertions(+), 476 deletions(-) diff --git a/src/Light/2band.jl b/src/Light/2band.jl index 778626a20..442653e7e 100644 --- a/src/Light/2band.jl +++ b/src/Light/2band.jl @@ -74,7 +74,7 @@ struct TwoBandPhotosyntheticallyActiveRadiation{FT, F, SPAR} end """ - TwoBandPhotosyntheticallyActiveRadiation(; grid, + TwoBandPhotosyntheticallyActiveRadiation(; grid::AbstractGrid{FT}, water_red_attenuation::FT = 0.225, # 1/m water_blue_attenuation::FT = 0.0232, # 1/m chlorophyll_red_attenuation::FT = 0.037, # 1/(m * (mgChl/m³) ^ eʳ) @@ -94,7 +94,7 @@ Keyword Arguments which should be `f(x, y, t)` where `x` and `y` are the native coordinates (i.e. meters for rectilinear grids and latitude/longitude as appropriate) """ -function TwoBandPhotosyntheticallyActiveRadiation(; grid, +function TwoBandPhotosyntheticallyActiveRadiation(; grid::AbstractGrid{FT}, water_red_attenuation::FT = 0.225, # 1/m water_blue_attenuation::FT = 0.0232, # 1/m chlorophyll_red_attenuation::FT = 0.037, # 1/(m * (mgChl/m³) ^ eʳ) diff --git a/src/Light/Light.jl b/src/Light/Light.jl index 04e982112..316fc4e58 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -13,7 +13,7 @@ using KernelAbstractions, Oceananigans.Units using Oceananigans.Architectures: device, architecture, on_architecture using Oceananigans.Utils: launch! using Oceananigans: Center, Face, fields -using Oceananigans.Grids: node, znodes, znode +using Oceananigans.Grids: node, znodes, znode, AbstractGrid using Oceananigans.Fields: CenterField, TracerFields, location using Oceananigans.BoundaryConditions: fill_halo_regions!, ValueBoundaryCondition, diff --git a/src/Light/compute_euphotic_depth.jl b/src/Light/compute_euphotic_depth.jl index 8652d9938..118e302c5 100644 --- a/src/Light/compute_euphotic_depth.jl +++ b/src/Light/compute_euphotic_depth.jl @@ -1,11 +1,11 @@ using Oceananigans.Fields: ConstantField, ZeroField -@kernel function _compute_euphotic_depth!(euphotic_depth, PAR, grid, cutoff) +@kernel function _compute_euphotic_depth!(euphotic_depth, PAR, grid::AbstractGrid{FT}, cutoff) where FT i, j = @index(Global, NTuple) surface_PAR = @inbounds (PAR[i, j, grid.Nz] + PAR[i, j, grid.Nz + 1])/2 - @inbounds euphotic_depth[i, j, 1] = -Inf + @inbounds euphotic_depth[i, j, 1] = convert(FT, -Inf) for k in grid.Nz-1:-1:1 PARₖ = @inbounds PAR[i, j, k] diff --git a/src/Light/multi_band.jl b/src/Light/multi_band.jl index bb3f55fe2..3955fc39e 100644 --- a/src/Light/multi_band.jl +++ b/src/Light/multi_band.jl @@ -37,7 +37,7 @@ struct MultiBandPhotosyntheticallyActiveRadiation{T, F, FN, K, E, C, SPAR, SPARD end """ - MultiBandPhotosyntheticallyActiveRadiation(; grid, + MultiBandPhotosyntheticallyActiveRadiation(; grid::AbstractGrid{FT}, bands = ((400, 500), (500, 600), (600, 700)), #nm base_bands = MOREL_λ, base_water_attenuation_coefficient = MOREL_kʷ, @@ -65,7 +65,7 @@ Keyword Arguments and latitude/longitude as appropriate) """ -function MultiBandPhotosyntheticallyActiveRadiation(; grid, +function MultiBandPhotosyntheticallyActiveRadiation(; grid::AbstractGrid{FT}, bands = ((400, 500), (500, 600), (600, 700)), #nm base_bands = MOREL_λ, base_water_attenuation_coefficient = MOREL_kʷ, @@ -73,7 +73,7 @@ function MultiBandPhotosyntheticallyActiveRadiation(; grid, base_chlorophyll_attenuation_coefficient = MOREL_χ, field_names = ntuple(n->par_symbol(n), Val(length(bands))), surface_PAR = default_surface_PAR, - surface_PAR_division = fill(1 / length(bands), length(bands))) + surface_PAR_division = fill(1 / length(bands), length(bands))) where FT Nbands = length(bands) kʷ = zeros(eltype(grid), Nbands) @@ -108,9 +108,9 @@ function MultiBandPhotosyntheticallyActiveRadiation(; grid, return MultiBandPhotosyntheticallyActiveRadiation(total_PAR, fields, tuple(field_names...), - kʷ, - e, - χ, + convert.(FT, kʷ), + convert.(FT, e), + convert.(FT, χ), surface_PAR, surface_PAR_division) end diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index e2a1df8c2..2914bdbe4 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -45,6 +45,7 @@ export LOBSTER using Oceananigans.Units using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField +using Oceananigans.Grids: AbstractGrid using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, default_surface_PAR using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, ScaleNegativeTracers @@ -169,7 +170,7 @@ struct LOBSTER{FT, B, W} <: AbstractContinuousFormBiogeochemistry end """ - LOBSTER(; grid, + LOBSTER(; grid::AbstractGrid{FT}, phytoplankton_preference::FT = 0.5, maximum_grazing_rate::FT = 9.26e-6, # 1/s grazing_half_saturation::FT = 1.0, # mmol N/m³ @@ -254,7 +255,7 @@ LOBSTER{Float64} with carbonates ❌, oxygen ❌, variable Redfield ratio ❌ an Modifiers: Nothing ``` """ -function LOBSTER(; grid, +function LOBSTER(; grid::AbstractGrid{FT}, phytoplankton_preference::FT = 0.5, maximum_grazing_rate::FT = 9.26e-6, # 1/s grazing_half_saturation::FT = 1.0, # mmol N/m³ @@ -457,9 +458,9 @@ const VariableRedfieldLobster = Union{LOBSTER{<:Any, <:Val{(false, false, true)} @inline redfield(::Val{:P}, bgc::LOBSTER) = (1 + bgc.organic_carbon_calcate_ratio) * bgc.phytoplankton_redfield @inline redfield(::Val{:Z}, bgc::LOBSTER) = bgc.phytoplankton_redfield -@inline redfield(::Union{Val{:NO₃}, Val{:NH₄}, Val{:Alk}, Val{:O₂}}, bgc::LOBSTER) = 0 +@inline redfield(::Union{Val{:NO₃}, Val{:NH₄}, Val{:Alk}, Val{:O₂}}, bgc::LOBSTER{FT}) where FT = convert(FT, 0) @inline redfield(::Union{Val{:sPOM}, Val{:bPOM}, Val{:DOM}}, bgc::LOBSTER) = bgc.organic_redfield -@inline redfield(::Union{Val{:sPOC}, Val{:bPOC}, Val{:DOC}, Val{:DIC}}, bgc::LOBSTER) = 1 +@inline redfield(::Union{Val{:sPOC}, Val{:bPOC}, Val{:DOC}, Val{:DIC}}, bgc::LOBSTER{FT}) where FT = convert(FT, 1) @inline redfield(i, j, k, ::Val{:sPON}, bgc::VariableRedfieldLobster, tracers) = @inbounds tracers.sPOC[i, j, k] / tracers.sPON[i, j, k] @inline redfield(i, j, k, ::Val{:bPON}, bgc::VariableRedfieldLobster, tracers) = @inbounds tracers.bPOC[i, j, k] / tracers.bPON[i, j, k] diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 96706e57d..74ab50586 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -21,6 +21,7 @@ using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry using Oceananigans.Units using Oceananigans.Fields: ZeroField +using Oceananigans.Grids: AbstractGrid using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, default_surface_PAR using OceanBioME: setup_velocity_fields, show_sinking_velocities @@ -93,7 +94,7 @@ struct NutrientPhytoplanktonZooplanktonDetritus{FT, W} <: AbstractContinuousForm end """ - NutrientPhytoplanktonZooplanktonDetritus(; grid, + NutrientPhytoplanktonZooplanktonDetritus(; grid::AbstractGrid{FT}, initial_photosynthetic_slope::FT = 0.1953 / day, # 1/(W/m²)/s base_maximum_growth::FT = 0.6989 / day, # 1/s nutrient_half_saturation::FT = 2.3868, # mmol N/m³ @@ -154,7 +155,7 @@ NutrientPhytoplanktonZooplanktonDetritus{Float64} model, with (:P, :D) sinking Modifiers: Nothing ``` """ -function NutrientPhytoplanktonZooplanktonDetritus(; grid, +function NutrientPhytoplanktonZooplanktonDetritus(; grid::AbstractGrid{FT}, initial_photosynthetic_slope::FT = 0.1953 / day, # 1/(W/m²)/s base_maximum_growth::FT = 0.6989 / day, # 1/s nutrient_half_saturation::FT = 2.3868, # mmol N/m³ @@ -316,8 +317,8 @@ adapt_structure(to, npzd::NPZD) = adapt(to, npzd.sinking_velocities)) @inline redfield(i, j, k, val_tracer_name, bgc::NPZD, tracers) = redfield(val_tracer_name, bgc) -@inline redfield(::Union{Val{:N}}, bgc::NPZD) = 0 -@inline redfield(::Union{Val{:P}, Val{:Z}, Val{:D}}, bgc::NPZD) = 6.56 +@inline redfield(::Union{Val{:N}}, bgc::NPZD{FT}) where FT = convert(FT, 0) +@inline redfield(::Union{Val{:P}, Val{:Z}, Val{:D}}, bgc::NPZD{FT}) where FT = convert(FT, 6.56) @inline nitrogen_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = sinking_flux(i, j, k, grid, advection, Val(:D), bgc, tracers) + sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers) @@ -329,6 +330,6 @@ adapt_structure(to, npzd::NPZD) = @inline conserved_tracers(::NPZD) = (:N, :P, :Z, :D) @inline sinking_tracers(bgc::NPZD) = keys(bgc.sinking_velocities) -@inline chlorophyll(bgc::NPZD, model) = 1.31 * model.tracers.P +@inline chlorophyll(bgc::NPZD{FT}, model) where FT = convert(FT, 1.31) * model.tracers.P end # module diff --git a/src/Models/AdvectedPopulations/PISCES/PISCES.jl b/src/Models/AdvectedPopulations/PISCES/PISCES.jl index 254921787..e78134bc3 100644 --- a/src/Models/AdvectedPopulations/PISCES/PISCES.jl +++ b/src/Models/AdvectedPopulations/PISCES/PISCES.jl @@ -36,7 +36,7 @@ using OceanBioME.Models.CarbonChemistryModel: CarbonChemistry using Oceananigans.Biogeochemistry: AbstractBiogeochemistry using Oceananigans.Fields: set! -using Oceananigans.Grids: φnodes, RectilinearGrid +using Oceananigans.Grids: φnodes, RectilinearGrid, AbstractGrid import OceanBioME: redfield, conserved_tracers, maximum_sinking_velocity, chlorophyll @@ -174,7 +174,7 @@ include("adapts.jl") include("show_methods.jl") """ - PISCES(; grid, + PISCES(; grid::AbstractGrid{FT}, phytoplankton = MixedMondoNanoAndDiatoms(), zooplankton = MicroAndMesoZooplankton(), dissolved_organic_matter = DissolvedOrganicCarbon(), @@ -284,67 +284,70 @@ the classes to a single `phytoplankton` if more classes are required (see `OceanBioME.Models.PISCESModel` docstring). Similarly, if a more generic `particulate_organic_matter` was desired a way to specify arbitary tracers for arguments would be required. """ -function PISCES(; grid, - phytoplankton = MixedMondoNanoAndDiatoms(), - zooplankton = MicroAndMesoZooplankton(), - dissolved_organic_matter = DissolvedOrganicCarbon(), - particulate_organic_matter = TwoCompartementCarbonIronParticles(), - - nitrogen = NitrateAmmonia(), - iron = SimpleIron(), - silicate = Silicate(), - oxygen = Oxygen(), - phosphate = Phosphate(), - - inorganic_carbon = InorganicCarbon(), - - # from Aumount 2005 rather than 2015 since it doesn't work the other way around - first_anoxia_thresehold = 6.0, - second_anoxia_thresehold = 1.0, - - nitrogen_redfield_ratio = 16/122, - phosphate_redfield_ratio = 1/122, - - mixed_layer_shear = 1.0, - background_shear = 0.01, - - latitude = PrescribedLatitude(45), - day_length = CBMDayLength(), - - mixed_layer_depth = Field{Center, Center, Nothing}(grid), - euphotic_depth = Field{Center, Center, Nothing}(grid), - - silicate_climatology = ConstantField(7.5), - - mean_mixed_layer_vertical_diffusivity = Field{Center, Center, Nothing}(grid), - mean_mixed_layer_light = Field{Center, Center, Nothing}(grid), - - carbon_chemistry = CarbonChemistry(), - calcite_saturation = CenterField(grid), - - surface_photosynthetically_active_radiation = default_surface_PAR, - - light_attenuation = - MultiBandPhotosyntheticallyActiveRadiation(; grid, - surface_PAR = surface_photosynthetically_active_radiation), - - sinking_speeds = (POC = 2/day, - # might be more efficient to just precompute this - GOC = Field(KernelFunctionOperation{Center, Center, Face}(DepthDependantSinkingSpeed(), - grid, - mixed_layer_depth, - euphotic_depth))), - open_bottom = true, - - scale_negatives = false, - invalid_fill_value = NaN, - - sediment = nothing, - particles = nothing, - modifiers = nothing) +function PISCES(FT = Float64; + grid, + phytoplankton = MixedMondoNanoAndDiatoms(FT), + zooplankton = MicroAndMesoZooplankton(FT), + dissolved_organic_matter = DissolvedOrganicCarbon(FT), + particulate_organic_matter = TwoCompartementCarbonIronParticles(FT), + + nitrogen = NitrateAmmonia{FT}(), + iron = SimpleIron{FT}(), + silicate = Silicate(), + oxygen = Oxygen{FT}(), + phosphate = Phosphate(), + + inorganic_carbon = InorganicCarbon(), + + # from Aumount 2005 rather than 2015 since it doesn't work the other way around + first_anoxia_thresehold = convert(FT, 6.0), + second_anoxia_thresehold = convert(FT, 1.0), + + nitrogen_redfield_ratio = convert(FT, 16/122), + phosphate_redfield_ratio = convert(FT, 1/122), + + mixed_layer_shear = convert(FT, 1.0), + background_shear = convert(FT, 0.01), + + latitude = PrescribedLatitude{FT}(45), + day_length = CBMDayLength(), + + mixed_layer_depth = Field{Center, Center, Nothing}(grid), + euphotic_depth = Field{Center, Center, Nothing}(grid), + + silicate_climatology = ConstantField(convert(FT, 7.5)), + + mean_mixed_layer_vertical_diffusivity = Field{Center, Center, Nothing}(grid), + mean_mixed_layer_light = Field{Center, Center, Nothing}(grid), + + carbon_chemistry = CarbonChemistry(FT), + calcite_saturation = CenterField(grid), + + surface_photosynthetically_active_radiation = default_surface_PAR, + + light_attenuation = + MultiBandPhotosyntheticallyActiveRadiation(; grid, + surface_PAR = surface_photosynthetically_active_radiation), + + sinking_speeds = (POC = convert(FT, 2/day), + # might be more efficient to just precompute this + GOC = Field(KernelFunctionOperation{Center, Center, Face}(DepthDependantSinkingSpeed{FT}(), + grid, + mixed_layer_depth, + euphotic_depth))), + open_bottom = true, + + scale_negatives = false, + invalid_fill_value = NaN, + + sediment = nothing, + particles = nothing, + modifiers = nothing) @warn "This implementation of PISCES is in early development and has not yet been validated against the operational version" + eltype(grid) == FT || @warn "The float type for parameters ($FT) does not match the grids float type ($(eltype(grid))). This will cause type instability." + if !isnothing(sediment) && !open_bottom @warn "You have specified a sediment model but not `open_bottom` which will not work as the tracer will settle in the bottom cell" end diff --git a/src/Models/AdvectedPopulations/PISCES/dissolved_organic_matter/dissolved_organic_carbon.jl b/src/Models/AdvectedPopulations/PISCES/dissolved_organic_matter/dissolved_organic_carbon.jl index 551e2960d..45f1e8f03 100644 --- a/src/Models/AdvectedPopulations/PISCES/dissolved_organic_matter/dissolved_organic_carbon.jl +++ b/src/Models/AdvectedPopulations/PISCES/dissolved_organic_matter/dissolved_organic_carbon.jl @@ -6,13 +6,32 @@ using Oceananigans.Grids: znode, Center Parameterisation of dissolved organic matter which depends on a bacterial concentration. """ -@kwdef struct DissolvedOrganicCarbon{FT, AP} - remineralisation_rate :: FT = 0.3/day # 1 / s - bacteria_concentration_depth_exponent :: FT = 0.684 # - reference_bacteria_concentration :: FT = 1.0 # mmol C / m³ - temperature_sensetivity :: FT = 1.066 # +struct DissolvedOrganicCarbon{FT, AP} + remineralisation_rate :: FT # 1 / s + bacteria_concentration_depth_exponent :: FT # + reference_bacteria_concentration :: FT # mmol C / m³ + temperature_sensetivity :: FT # # (1 / (mmol C / m³), 1 / (mmol C / m³), 1 / (mmol C / m³), 1 / (mmol C / m³) / s, 1 / (mmol C / m³) / s) - aggregation_parameters :: AP = (0.37, 102, 3530, 5095, 114) .* (10^-6 / day) + aggregation_parameters :: AP + + function DissolvedOrganicCarbon(FT = Float64; + remineralisation_rate = 0.3/day, # 1 / s + bacteria_concentration_depth_exponent = 0.684, # + reference_bacteria_concentration = 1.0, # mmol C / m³ + temperature_sensetivity = 1.066, # +# (1 / (mmol C / m³), 1 / (mmol C / m³), 1 / (mmol C / m³), 1 / (mmol C / m³) / s, 1 / (mmol C / m³) / s) + aggregation_parameters = (0.37, 102, 3530, 5095, 114) .* (10^-6 / day)) + + aggregation_parameters = convert.(FT, aggregation_parameters) + + AP = typeof(aggregation_parameters) + + return new{FT, AP}(convert(FT, remineralisation_rate), + convert(FT, bacteria_concentration_depth_exponent), + convert(FT, reference_bacteria_concentration), + convert(FT, temperature_sensetivity), + aggregation_parameters) + end end required_biogeochemical_tracers(::DissolvedOrganicCarbon) = tuple(:DOC) diff --git a/src/Models/AdvectedPopulations/PISCES/particulate_organic_matter/two_size_class.jl b/src/Models/AdvectedPopulations/PISCES/particulate_organic_matter/two_size_class.jl index 23ef232e1..f7a434fd3 100644 --- a/src/Models/AdvectedPopulations/PISCES/particulate_organic_matter/two_size_class.jl +++ b/src/Models/AdvectedPopulations/PISCES/particulate_organic_matter/two_size_class.jl @@ -14,28 +14,73 @@ and large carbon classes, `SFe` and `BFe` for the small and ̶l̶a̶r̶g̶e̶ b compartements, and `PSi` for the ̶l̶a̶r̶g̶e̶ particulate silicon (*not* the phytoplankton silicon). """ -@kwdef struct TwoCompartementCarbonIronParticles{FT, AP} - temperature_sensetivity :: FT = 1.066 # - base_breakdown_rate :: FT = 0.025 / day # 1 / s +struct TwoCompartementCarbonIronParticles{FT, AP} + temperature_sensetivity :: FT # + base_breakdown_rate :: FT # 1 / s # (1 / (mmol C / m³), 1 / (mmol C / m³), 1 / (mmol C / m³) / s, 1 / (mmol C / m³) / s) - aggregation_parameters :: AP = (25.9, 4452, 3.3, 47.1) .* (10^-6 / day) - minimum_iron_scavenging_rate :: FT = 3e-5/day # 1 / s - load_specific_iron_scavenging_rate :: FT = 0.005/day # 1 / (mmol C / m³) / s - bacterial_iron_uptake_efficiency :: FT = 0.16 # - small_fraction_of_bacterially_consumed_iron :: FT = 0.12 / bacterial_iron_uptake_efficiency - large_fraction_of_bacterially_consumed_iron :: FT = 0.04 / bacterial_iron_uptake_efficiency - base_liable_silicate_fraction :: FT = 0.5 # - fast_dissolution_rate_of_silicate :: FT = 0.025/day # 1 / s - slow_dissolution_rate_of_silicate :: FT = 0.003/day # 1 / s + aggregation_parameters :: AP + minimum_iron_scavenging_rate :: FT # 1 / s + load_specific_iron_scavenging_rate :: FT # 1 / (mmol C / m³) / s + bacterial_iron_uptake_efficiency :: FT # + small_fraction_of_bacterially_consumed_iron :: FT + large_fraction_of_bacterially_consumed_iron :: FT + base_liable_silicate_fraction :: FT # + fast_dissolution_rate_of_silicate :: FT # 1 / s + slow_dissolution_rate_of_silicate :: FT # 1 / s # calcite - base_calcite_dissolution_rate :: FT = 0.197 / day # 1 / s - calcite_dissolution_exponent :: FT = 1.0 # + base_calcite_dissolution_rate :: FT # 1 / s + calcite_dissolution_exponent :: FT # # iron in particles - maximum_iron_ratio_in_bacteria :: FT = 0.06 # μmol Fe / mmol C - iron_half_saturation_for_bacteria :: FT = 0.3 # μmol Fe / m³ - maximum_bacterial_growth_rate :: FT = 0.6 / day # 1 / s + maximum_iron_ratio_in_bacteria :: FT # μmol Fe / mmol C + iron_half_saturation_for_bacteria :: FT # μmol Fe / m³ + maximum_bacterial_growth_rate :: FT # 1 / s + + function TwoCompartementCarbonIronParticles(FT = Float64; + temperature_sensetivity = 1.066, # + base_breakdown_rate = 0.025 / day, # 1 / s +# (1 / (mmol C / m³), 1 / (mmol C / m³), 1 / (mmol C / m³) / s, 1 / (mmol C / m³) / s) + aggregation_parameters = (25.9, 4452, 3.3, 47.1) .* (10^-6 / day), + minimum_iron_scavenging_rate = 3e-5/day, # 1 / s + load_specific_iron_scavenging_rate = 0.005/day, # 1 / (mmol C / m³) / s + bacterial_iron_uptake_efficiency = 0.16, # + small_fraction_of_bacterially_consumed_iron = 0.12 / bacterial_iron_uptake_efficiency, + large_fraction_of_bacterially_consumed_iron = 0.04 / bacterial_iron_uptake_efficiency, + base_liable_silicate_fraction = 0.5, # + fast_dissolution_rate_of_silicate = 0.025/day, # 1 / s + slow_dissolution_rate_of_silicate = 0.003/day, # 1 / s + + # calcite + base_calcite_dissolution_rate = 0.197 / day, # 1 / s + calcite_dissolution_exponent = 1.0, # + + # iron in particles + maximum_iron_ratio_in_bacteria = 0.06, # μmol Fe / mmol C + iron_half_saturation_for_bacteria = 0.3, # μmol Fe / m³ + maximum_bacterial_growth_rate = 0.6 / day) # 1 / s + + aggregation_parameters = convert.(FT, aggregation_parameters) + + AP = typeof(aggregation_parameters) + + return new{FT, AP}(convert(FT, temperature_sensetivity), + convert(FT, base_breakdown_rate), + aggregation_parameters, + convert(FT, minimum_iron_scavenging_rate), + convert(FT, load_specific_iron_scavenging_rate), + convert(FT, bacterial_iron_uptake_efficiency), + convert(FT, small_fraction_of_bacterially_consumed_iron), + convert(FT, large_fraction_of_bacterially_consumed_iron), + convert(FT, base_liable_silicate_fraction), + convert(FT, fast_dissolution_rate_of_silicate), + convert(FT, slow_dissolution_rate_of_silicate), + convert(FT, base_calcite_dissolution_rate), + convert(FT, calcite_dissolution_exponent), + convert(FT, maximum_iron_ratio_in_bacteria), + convert(FT, iron_half_saturation_for_bacteria), + convert(FT, maximum_bacterial_growth_rate)) + end end const TwoCompartementPOCPISCES = PISCES{<:Any, <:Any, <:Any, <:TwoCompartementCarbonIronParticles} diff --git a/src/Models/AdvectedPopulations/PISCES/phytoplankton/mixed_mondo.jl b/src/Models/AdvectedPopulations/PISCES/phytoplankton/mixed_mondo.jl index 56807ea6c..6ab901b56 100644 --- a/src/Models/AdvectedPopulations/PISCES/phytoplankton/mixed_mondo.jl +++ b/src/Models/AdvectedPopulations/PISCES/phytoplankton/mixed_mondo.jl @@ -20,35 +20,76 @@ either `NutrientLimitedProduction` or `GrowthRespirationLimitedProduction`, which represent the typical and `newprod` versions of PISCES. """ -@kwdef struct MixedMondo{GR, NL, FT} +struct MixedMondo{GR, NL, FT} growth_rate :: GR nutrient_limitation :: NL - exudated_fracton :: FT = 0.05 # - - blue_light_absorption :: FT # - green_light_absorption :: FT # - red_light_absorption :: FT # - - mortality_half_saturation :: FT = 0.2 # mmol C / m³ - linear_mortality_rate :: FT = 0.01 / day # 1 / s - - base_quadratic_mortality :: FT = 0.01 / day # 1 / s / (mmol C / m³) - maximum_quadratic_mortality :: FT # 1 / s / (mmol C / m³) - zero for nanophytoplankton - - minimum_chlorophyll_ratio :: FT = 0.0033 # mg Chl / mg C - maximum_chlorophyll_ratio :: FT # mg Chl / mg C - - maximum_iron_ratio :: FT = 0.06 # μmol Fe / mmol C - - silicate_half_saturation :: FT = 2.0 # mmol Si / m³ - enhanced_silicate_half_saturation :: FT = 20.9 # mmol Si / m³ - optimal_silicate_ratio :: FT = 0.159 # mmol Si / mmol C - - half_saturation_for_iron_uptake :: FT # μmol Fe / m³ - - threshold_for_size_dependency :: FT = 1.0 # mmol C / m³ - size_ratio :: FT = 3.0 # + exudated_fracton :: FT # + + blue_light_absorption :: FT # + green_light_absorption :: FT # + red_light_absorption :: FT # + + mortality_half_saturation :: FT # mmol C / m³ + linear_mortality_rate :: FT # 1 / s + + base_quadratic_mortality :: FT # 1 / s / (mmol C / m³) + maximum_quadratic_mortality :: FT # 1 / s / (mmol C / m³) - zero for nanophytoplankton + + minimum_chlorophyll_ratio :: FT # mg Chl / mg C + maximum_chlorophyll_ratio :: FT # mg Chl / mg C + + maximum_iron_ratio :: FT # μmol Fe / mmol C + + silicate_half_saturation :: FT # mmol Si / m³ + enhanced_silicate_half_saturation :: FT # mmol Si / m³ + optimal_silicate_ratio :: FT # mmol Si / mmol C + + half_saturation_for_iron_uptake :: FT # μmol Fe / m³ + + threshold_for_size_dependency :: FT # mmol C / m³ + size_ratio :: FT # + + function MixedMondo(FT = Float64; + growth_rate::GR, + nutrient_limitation::NL, + exudated_fracton = 0.05, # + + blue_light_absorption, # + green_light_absorption, # + red_light_absorption, # + + mortality_half_saturation = 0.2, # mmol C / m³ + linear_mortality_rate = 0.01 / day, # 1 / s + + base_quadratic_mortality = 0.01 / day, # 1 / s / (mmol C / m³) + maximum_quadratic_mortality, # 1 / s / (mmol C / m³) - zero for nanophytoplankton + + minimum_chlorophyll_ratio = 0.0033, # mg Chl / mg C + maximum_chlorophyll_ratio, # mg Chl / mg C + + maximum_iron_ratio = 0.06, # μmol Fe / mmol C + + silicate_half_saturation = 2.0, # mmol Si / m³ + enhanced_silicate_half_saturation = 20.9, # mmol Si / m³ + optimal_silicate_ratio = 0.159, # mmol Si / mmol C + + half_saturation_for_iron_uptake, # μmol Fe / m³ + + threshold_for_size_dependency = 1.0, # mmol C / m³ + size_ratio = 3.0) where {GR, NL} # + + return new{GR, NL, FT}(growth_rate, nutrient_limitation, + exudated_fracton, + blue_light_absorption, green_light_absorption, red_light_absorption, + mortality_half_saturation, linear_mortality_rate, + base_quadratic_mortality, maximum_quadratic_mortality, + minimum_chlorophyll_ratio, maximum_chlorophyll_ratio, + maximum_iron_ratio, + silicate_half_saturation, enhanced_silicate_half_saturation, optimal_silicate_ratio, + half_saturation_for_iron_uptake, + threshold_for_size_dependency, size_ratio) + end end required_biogeochemical_tracers(phyto::MixedMondo, base) = diff --git a/src/Models/AdvectedPopulations/PISCES/phytoplankton/mixed_mondo_nano_diatoms.jl b/src/Models/AdvectedPopulations/PISCES/phytoplankton/mixed_mondo_nano_diatoms.jl index 833160801..73b66fd48 100644 --- a/src/Models/AdvectedPopulations/PISCES/phytoplankton/mixed_mondo_nano_diatoms.jl +++ b/src/Models/AdvectedPopulations/PISCES/phytoplankton/mixed_mondo_nano_diatoms.jl @@ -1,27 +1,28 @@ -function MixedMondoNanoAndDiatoms(; nano = MixedMondo(growth_rate = GrowthRespirationLimitedProduction(dark_tollerance = 3days), - nutrient_limitation = - NitrogenIronPhosphateSilicateLimitation(minimum_ammonium_half_saturation = 0.013, - minimum_nitrate_half_saturation = 0.13, - minimum_phosphate_half_saturation = 0.8, - silicate_limited = false), - blue_light_absorption = 2.1, - green_light_absorption = 0.42, - red_light_absorption = 0.4, - maximum_quadratic_mortality = 0.0, - maximum_chlorophyll_ratio = 0.033, - half_saturation_for_iron_uptake = 1.0), - diatoms = MixedMondo(growth_rate = GrowthRespirationLimitedProduction(dark_tollerance = 4days), - nutrient_limitation = - NitrogenIronPhosphateSilicateLimitation(minimum_ammonium_half_saturation = 0.039, - minimum_nitrate_half_saturation = 0.39, - minimum_phosphate_half_saturation = 2.4, - silicate_limited = true), - blue_light_absorption = 1.6, - green_light_absorption = 0.69, - red_light_absorption = 0.7, - maximum_quadratic_mortality = 0.03/day, - maximum_chlorophyll_ratio = 0.05, - half_saturation_for_iron_uptake = 3.0)) +function MixedMondoNanoAndDiatoms(FT = Float64; + nano = MixedMondo(FT; growth_rate = GrowthRespirationLimitedProduction{FT}(dark_tollerance = 3days), + nutrient_limitation = + NitrogenIronPhosphateSilicateLimitation{FT}(minimum_ammonium_half_saturation = 0.013, + minimum_nitrate_half_saturation = 0.13, + minimum_phosphate_half_saturation = 0.8, + silicate_limited = false), + blue_light_absorption = convert(FT, 2.1), + green_light_absorption = convert(FT, 0.42), + red_light_absorption = convert(FT, 0.4), + maximum_quadratic_mortality = convert(FT, 0.0), + maximum_chlorophyll_ratio = convert(FT, 0.033), + half_saturation_for_iron_uptake = convert(FT, 1.0)), + diatoms = MixedMondo(FT; growth_rate = GrowthRespirationLimitedProduction{FT}(dark_tollerance = 4days), + nutrient_limitation = + NitrogenIronPhosphateSilicateLimitation{FT}(minimum_ammonium_half_saturation = 0.039, + minimum_nitrate_half_saturation = 0.39, + minimum_phosphate_half_saturation = 2.4, + silicate_limited = true), + blue_light_absorption = convert(FT, 1.6), + green_light_absorption = convert(FT, 0.69), + red_light_absorption = convert(FT, 0.7), + maximum_quadratic_mortality = convert(FT, 0.03/day), + maximum_chlorophyll_ratio = convert(FT, 0.05), + half_saturation_for_iron_uptake = convert(FT, 3.0))) return NanoAndDiatoms(; nano, diatoms) end diff --git a/src/Models/AdvectedPopulations/PISCES/phytoplankton/nutrient_limitation.jl b/src/Models/AdvectedPopulations/PISCES/phytoplankton/nutrient_limitation.jl index 92bc393ab..9dc12ac51 100644 --- a/src/Models/AdvectedPopulations/PISCES/phytoplankton/nutrient_limitation.jl +++ b/src/Models/AdvectedPopulations/PISCES/phytoplankton/nutrient_limitation.jl @@ -7,12 +7,12 @@ iron (Fe), phosphate PO₄, and (optionally) silicate (Si) availability. Silicate limitation may be turned off (e.g. for nanophytoplankton) by setting `silicate_limited=false`. """ -@kwdef struct NitrogenIronPhosphateSilicateLimitation{FT, BT} +@kwdef struct NitrogenIronPhosphateSilicateLimitation{FT} minimum_ammonium_half_saturation :: FT # mmol N / m³ minimum_nitrate_half_saturation :: FT # mmol N / m³ minimum_phosphate_half_saturation :: FT # mmol P / m³ optimal_iron_quota :: FT = 0.007 # μmol Fe / mmol C - silicate_limited :: BT # Bool + silicate_limited :: Bool # Bool minimum_silicate_half_saturation :: FT = 1.0 # mmol Si / m³ silicate_half_saturation_parameter :: FT = 16.6 # mmol Si / m³ end diff --git a/src/Models/AdvectedPopulations/PISCES/zooplankton/defaults.jl b/src/Models/AdvectedPopulations/PISCES/zooplankton/defaults.jl index 64b7698c9..111e199d2 100644 --- a/src/Models/AdvectedPopulations/PISCES/zooplankton/defaults.jl +++ b/src/Models/AdvectedPopulations/PISCES/zooplankton/defaults.jl @@ -1,23 +1,23 @@ # this file sets up the default configuration of Z and M which graze on P, D, (Z, ) and POC -function MicroAndMesoZooplankton(; - micro = QualityDependantZooplankton(maximum_grazing_rate = 3/day, - food_preferences = (P = 1.0, D = 0.5, POC = 0.1, Z = 0), - quadratic_mortality = 0.004/day, - linear_mortality = 0.03/day, - minimum_growth_efficiency = 0.3, - maximum_flux_feeding_rate = 0.0, - undissolved_calcite_fraction = 0.5, - iron_ratio = 0.01), - meso = QualityDependantZooplankton(maximum_grazing_rate = 0.75/day, - food_preferences = (P = 0.3, D = 1.0, POC = 0.3, Z = 1.0), - quadratic_mortality = 0.03/day, - linear_mortality = 0.005/day, - minimum_growth_efficiency = 0.35, - maximum_flux_feeding_rate = 2e3 / 1e6, - undissolved_calcite_fraction = 0.75, - iron_ratio = 0.015)) - - return MicroAndMeso(; micro, meso) +function MicroAndMesoZooplankton(FT = Float64; + micro::MU = QualityDependantZooplankton(FT; maximum_grazing_rate = 3/day, + food_preferences = (P = 1.0, D = 0.5, POC = 0.1, Z = 0), + quadratic_mortality = 0.004/day, + linear_mortality = 0.03/day, + minimum_growth_efficiency = 0.3, + maximum_flux_feeding_rate = 0.0, + undissolved_calcite_fraction = 0.5, + iron_ratio = 0.01), + meso::ME = QualityDependantZooplankton(FT; maximum_grazing_rate = 0.75/day, + food_preferences = (P = 0.3, D = 1.0, POC = 0.3, Z = 1.0), + quadratic_mortality = 0.03/day, + linear_mortality = 0.005/day, + minimum_growth_efficiency = 0.35, + maximum_flux_feeding_rate = 2e3 / 1e6, + undissolved_calcite_fraction = 0.75, + iron_ratio = 0.015)) where {MU, ME} + + return MicroAndMeso{MU, ME, FT}(; micro, meso) end @inline concentration(::Val{:P}, i, j, k, fields) = @inbounds fields.P[i, j, k] diff --git a/src/Models/AdvectedPopulations/PISCES/zooplankton/food_quality_dependant.jl b/src/Models/AdvectedPopulations/PISCES/zooplankton/food_quality_dependant.jl index 2b705ea9c..f2ff84a7e 100644 --- a/src/Models/AdvectedPopulations/PISCES/zooplankton/food_quality_dependant.jl +++ b/src/Models/AdvectedPopulations/PISCES/zooplankton/food_quality_dependant.jl @@ -8,31 +8,71 @@ particulates (POC and GOC). This model assumes a fixed ratio for all other elements (i.e. N, P, Fe). """ -@kwdef struct QualityDependantZooplankton{FT, FP} - temperature_sensetivity :: FT = 1.079 # - maximum_grazing_rate :: FT # 1 / s +struct QualityDependantZooplankton{FT, FP} + temperature_sensetivity :: FT # + maximum_grazing_rate :: FT # 1 / s food_preferences :: FP - food_threshold_concentration :: FT = 0.3 # mmol C / m³ - specific_food_thresehold_concentration :: FT = 0.001 # mmol C / m³ + food_threshold_concentration :: FT # mmol C / m³ + specific_food_thresehold_concentration :: FT # mmol C / m³ - grazing_half_saturation :: FT = 20.0 # mmol C / m³ + grazing_half_saturation :: FT # mmol C / m³ - maximum_flux_feeding_rate :: FT # m / (mmol C / m³) + maximum_flux_feeding_rate :: FT # m / (mmol C / m³) - iron_ratio :: FT # μmol Fe / mmol C + iron_ratio :: FT # μmol Fe / mmol C - minimum_growth_efficiency :: FT # - non_assililated_fraction :: FT = 0.3 # + minimum_growth_efficiency :: FT # + non_assililated_fraction :: FT # - mortality_half_saturation :: FT = 0.2 # mmol C / m³ - quadratic_mortality :: FT # 1 / (mmol C / m³) / s - linear_mortality :: FT # 1 / s + mortality_half_saturation :: FT # mmol C / m³ + quadratic_mortality :: FT # 1 / (mmol C / m³) / s + linear_mortality :: FT # 1 / s # this should be called inorganic excretion factor - dissolved_excretion_fraction :: FT = 0.6 # - undissolved_calcite_fraction :: FT # + dissolved_excretion_fraction :: FT # + undissolved_calcite_fraction :: FT # + + function QualityDependantZooplankton(FT = Float64; + temperature_sensetivity = 1.079, # + maximum_grazing_rate, # 1 / s + + food_preferences, + + food_threshold_concentration = 0.3, # mmol C / m³ + specific_food_thresehold_concentration = 0.001, # mmol C / m³ + + grazing_half_saturation = 20.0, # mmol C / m³ + + maximum_flux_feeding_rate, # m / (mmol C / m³) + + iron_ratio, # μmol Fe / mmol C + + minimum_growth_efficiency, # + non_assililated_fraction = 0.3, # + + mortality_half_saturation = 0.2, # mmol C / m³ + quadratic_mortality, # 1 / (mmol C / m³) / s + linear_mortality, # 1 / s + + # this should be called inorganic excretion factor + dissolved_excretion_fraction = 0.6, # + undissolved_calcite_fraction) # + + food_preferences = NamedTuple{keys(food_preferences)}(map(fp -> convert(FT, fp), food_preferences)) + + FP = typeof(food_preferences) + + return new{FT, FP}(temperature_sensetivity, maximum_grazing_rate, + food_preferences, + food_threshold_concentration, specific_food_thresehold_concentration, + grazing_half_saturation, maximum_flux_feeding_rate, + iron_ratio, + minimum_growth_efficiency, non_assililated_fraction, + mortality_half_saturation, quadratic_mortality, linear_mortality, + dissolved_excretion_fraction, undissolved_calcite_fraction) + end end required_biogeochemical_tracers(::QualityDependantZooplankton, name_base) = tuple(name_base) diff --git a/src/Models/CarbonChemistry/carbon_chemistry.jl b/src/Models/CarbonChemistry/carbon_chemistry.jl index 48eeeedf9..2aa9e9175 100644 --- a/src/Models/CarbonChemistry/carbon_chemistry.jl +++ b/src/Models/CarbonChemistry/carbon_chemistry.jl @@ -1,18 +1,33 @@ using Roots using OceanBioME.Models: teos10_polynomial_approximation +struct CarbonChemistry{P0, PC, PB, PS, PF, PP, PSi, PW, IS, PKS, PRho} + ionic_strength :: IS + solubility :: P0 + carbonic_acid :: PC + boric_acid :: PB + water :: PW + sulfate :: PS + fluoride :: PF + phosphoric_acid :: PP + silicic_acid :: PSi + calcite_solubility :: PKS + density_function :: PRho +end + """ - CarbonChemistry(; ionic_strength = IonicStrength(), - solubility = K0(), - carbonic_acid = (K1 = K1(), K2 = K2()), - boric_acid = KB(), - water = KW(), - sulfate = KS(; ionic_strength), - fluoride = KF(; ionic_strength), - phosphoric_acid = (KP1 = KP1(), KP2 = KP2(), KP3 = KP3()), - silicic_acid = KSi(; ionic_strength), - calcite_solubility = KSP_calcite(), - density_function = teos10_polynomial_approximation) + CarbonChemistry(FT = Float64; + ionic_strength = IonicStrength(), + solubility = K0(), + carbonic_acid = (K1 = K1(), K2 = K2()), + boric_acid = KB(), + water = KW(), + sulfate = KS(; ionic_strength), + fluoride = KF(; ionic_strength), + phosphoric_acid = (KP1 = KP1(), KP2 = KP2(), KP3 = KP3()), + silicic_acid = KSi(; ionic_strength), + calcite_solubility = KSP_calcite(), + density_function = teos10_polynomial_approximation) Carbon chemistry model capable of solving for sea water pCO₂ from DIC and total alkalinity or DIC and pH. @@ -42,20 +57,27 @@ julia> pCO₂_higher_pH = carbon_chemistry(; DIC = 2000, T = 10, S = 35, pH = 7. ``` """ -@kwdef struct CarbonChemistry{P0, PC, PB, PS, PF, PP, PSi, PW, IS, PKS, PRho} - ionic_strength :: IS = IonicStrength() - solubility :: P0 = K0() - carbonic_acid :: PC = (K1 = K1(), K2 = K2()) - boric_acid :: PB = KB() - water :: PW = KW() - sulfate :: PS = KS(; ionic_strength) - fluoride :: PF = KF(; ionic_strength, sulfate_constant = sulfate) - phosphoric_acid :: PP = (KP1 = KP1(), KP2 = KP2(), KP3 = KP3()) - silicic_acid :: PSi = KSi(; ionic_strength) - calcite_solubility :: PKS = KSP_calcite() - density_function :: PRho = teos10_polynomial_approximation + + + +function CarbonChemistry(FT = Float64; + ionic_strength = IonicStrength{FT}(), + solubility = K0{FT}(), + carbonic_acid = (K1 = K1(FT), K2 = K2(FT)), + boric_acid = KB(FT), + water = KW(FT), + sulfate = KS(FT; ionic_strength), + fluoride = KF(FT; ionic_strength), + phosphoric_acid = (KP1 = KP1(FT), KP2 = KP2(FT), KP3 = KP3(FT)), + silicic_acid = KSi(FT; ionic_strength), + calcite_solubility = KSP_calcite(FT), + density_function = teos10_polynomial_approximation) # the denisity function *is* going to cause type instability but I can't see a way to fix it + + return CarbonChemistry(ionic_strength, solubility, carbonic_acid, boric_acid, water, + sulfate, fluoride, phosphoric_acid, silicic_acid, calcite_solubility, density_function) end + """ alkalinity_residual(H, p) diff --git a/src/Models/CarbonChemistry/equilibrium_constants.jl b/src/Models/CarbonChemistry/equilibrium_constants.jl index 07e3856b7..b4864dcbf 100644 --- a/src/Models/CarbonChemistry/equilibrium_constants.jl +++ b/src/Models/CarbonChemistry/equilibrium_constants.jl @@ -1,20 +1,29 @@ """ - PressureCorrection(a₀, a₁, a₂, + PressureCorrection(FT=Float64; + a₀, a₁, a₂, b₀, b₁, b₂, - R) + R = 83.14472) Parameterisation for the pressure effect on thermodynamic constants. Form from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341). """ -@kwdef struct PressureCorrection{FT} +struct PressureCorrection{FT} a₀ :: FT a₁ :: FT a₂ :: FT b₀ :: FT b₁ :: FT - R :: FT = 83.14472 + R :: FT + + function PressureCorrection(FT = Float64; + a₀, a₁, a₂, + b₀, b₁, + R = 83.14472) + + return new{FT}(a₀, a₁, a₂, b₀, b₁, R) + end end @inline function (pc::PressureCorrection)(Tk, P) @@ -28,7 +37,7 @@ end return exp((-ΔV + 0.5 * Δκ * P) * P / RT) end -@inline (pc::PressureCorrection)(T, ::Nothing) = 1 +@inline (pc::PressureCorrection{FT})(T, ::Nothing) where FT = convert(FT, 1) summary(::IO, ::PressureCorrection) = string("Equilibrium constant pressure correction") show(io::IO, pc::PressureCorrection) = print(io, "Equilibrium constant pressure correction\n", @@ -75,13 +84,13 @@ show(io::IO, k0::K0) = print(io, "Solubility constant\n", " ln(k₀/k°) = $(k0.constant) + $(k0.inverse_T) / T + $(k0.log_T) (log(T) - log(100)) + $(k0.T²) T² + ($(k0.S) + $(k0.ST) T + $(k0.ST²) T²)S") """ - K1(; constant = 61.2172, - inverse_T = -3633.86, - log_T = -9.67770, - S = 0.011555, - S² = -0.0001152, - pressure_correction = - PressureCorrection(; a₀=-25.50, a₁=0.1271, a₂=0.0, b₀=-0.00308, b₁=0.0000877)) + K1(FT = Float64; + constant = 61.2172, + inverse_T = -3633.86, + log_T = -9.67770, + S = 0.011555, + S² = -0.0001152, + pressure_correction = PressureCorrection(FT; a₀=-25.50, a₁=0.1271, a₂=0.0, b₀=-0.00308, b₁=0.0000877)) Parameterisation for aquious carbon dioxide - bicarbonate dissociation equilibrium constant. @@ -91,15 +100,25 @@ Parameterisation for aquious carbon dioxide - bicarbonate dissociation equilibri Default values from Lueker et al. (2000, Mar. Chem., 70: 105–119). """ -@kwdef struct K1{FT, PC} - constant :: FT = 61.2172 - inverse_T :: FT = -3633.86 - log_T :: FT = -9.67770 - S :: FT = 0.011555 - S² :: FT = -0.0001152 +struct K1{FT, PC} + constant :: FT + inverse_T :: FT + log_T :: FT + S :: FT + S² :: FT + + pressure_correction :: PC - pressure_correction :: PC = - PressureCorrection(; a₀=-25.50, a₁=0.1271, a₂=0.0, b₀=-0.00308, b₁=0.0000877) + function K1(FT = Float64; + constant = 61.2172, + inverse_T = -3633.86, + log_T = -9.67770, + S = 0.011555, + S² = -0.0001152, + pressure_correction::PC = PressureCorrection(FT; a₀=-25.50, a₁=0.1271, a₂=0.0, b₀=-0.00308, b₁=0.0000877)) where PC + + return new{FT, PC}(constant, inverse_T, log_T, S, S², pressure_correction) + end end @inline (c::K1)(T, S; P = nothing) = @@ -111,13 +130,13 @@ show(io::IO, k1::K1) = print(io, "First carbon dioxide dissociation constant\n", " log₁₀(k₁/k°) = $(k1.constant) + $(k1.inverse_T) / T + $(k1.log_T) log(T) + $(k1.S) S + $(k1.S²) S²") """ - K2(; constant = -25.9290, - inverse_T = -471.78, - log_T = 3.16967, - S = 0.01781, - S² = -0.0001122, - pressure_correction = - PressureCorrection(; a₀=-15.82, a₁=-0.0219, a₂=0.0, b₀=0.00113, b₁=-0.0001475)) + K2(FT = Float64; + constant = -25.9290, + inverse_T = -471.78, + log_T = 3.16967, + S = 0.01781, + S² = -0.0001122, + pressure_correction = PressureCorrection(FT; a₀=-15.82, a₁=-0.0219, a₂=0.0, b₀=0.00113, b₁=-0.0001475)) Parameterisation for bicarbonate dissociation equilibrium constant. @@ -127,15 +146,25 @@ Parameterisation for bicarbonate dissociation equilibrium constant. Default values from Lueker et al. (2000, Mar. Chem., 70: 105–119). """ -@kwdef struct K2{FT, PC} - constant :: FT = -25.9290 - inverse_T :: FT = -471.78 - log_T :: FT = 3.16967 - S :: FT = 0.01781 - S² :: FT = -0.0001122 +struct K2{FT, PC} + constant :: FT + inverse_T :: FT + log_T :: FT + S :: FT + S² :: FT + + pressure_correction :: PC - pressure_correction :: PC = - PressureCorrection(; a₀=-15.82, a₁=-0.0219, a₂=0.0, b₀=0.00113, b₁=-0.0001475) + function K2(FT = Float64; + constant = -25.9290, + inverse_T = -471.78, + log_T = 3.16967, + S = 0.01781, + S² = -0.0001122, + pressure_correction::PC = PressureCorrection(FT; a₀=-15.82, a₁=-0.0219, a₂=0.0, b₀=0.00113, b₁=-0.0001475)) where PC + + return new{FT, PC}(constant, inverse_T, log_T, S, S², pressure_correction) + end end @inline (c::K2)(T, S; P = nothing) = @@ -147,20 +176,20 @@ show(io::IO, k2::K2) = print(io, "Second carbon dioxide dissociation constant\n" " log₁₀(k₂/k°) = $(k2.constant) + $(k2.inverse_T) / T + $(k2.log_T) log(T) + $(k2.S) S + $(k2.S²) S²") """ - KB(; constant = 148.0248, - inverse_T = -8966.90, - inverse_T_sqrt_S = -2890.53, - inverse_T_S = -77.942, - inverse_T_sqrt_S³ = 1.728, - inverse_T_S² = -0.0996, - sqrt_S = 137.1942, - S = 1.62142, - log_T = -24.4344, - log_T_sqrt_S = -25.085, - S_log_T = -0.2474, - T_sqrt_S = 0.053105, - pressure_correction = - PressureCorrection(; a₀=-29.48, a₁=0.1622, a₂=-0.0026080, b₀=-0.00284, b₁=0.0)) + KB(FT = Float64; + constant = 148.0248, + inverse_T = -8966.90, + inverse_T_sqrt_S = -2890.53, + inverse_T_S = -77.942, + inverse_T_sqrt_S³ = 1.728, + inverse_T_S² = -0.0996, + sqrt_S = 137.1942, + S = 1.62142, + log_T = -24.4344, + log_T_sqrt_S = -25.085, + S_log_T = -0.2474, + T_sqrt_S = 0.053105, + pressure_correction = PressureCorrection(FT; a₀=-29.48, a₁=0.1622, a₂=-0.0026080, b₀=-0.00284, b₁=0.0)) Parameterisation for boric acid equilibrium with water. @@ -170,22 +199,45 @@ Parameterisation for boric acid equilibrium with water. Default values from Dickson (1990, Deep-Sea Res., 37, 755–766). """ -@kwdef struct KB{FT, PC} - constant :: FT = 148.0248 - inverse_T :: FT = -8966.90 - inverse_T_sqrt_S :: FT = -2890.53 - inverse_T_S :: FT = -77.942 - inverse_T_sqrt_S³ :: FT = 1.728 - inverse_T_S² :: FT = -0.0996 - sqrt_S :: FT = 137.1942 - S :: FT = 1.62142 - log_T :: FT = -24.4344 - log_T_sqrt_S :: FT = -25.085 - S_log_T :: FT = -0.2474 - T_sqrt_S :: FT = 0.053105 - - pressure_correction :: PC = - PressureCorrection(; a₀=-29.48, a₁=0.1622, a₂=-0.0026080, b₀=-0.00284, b₁=0.0) +struct KB{FT, PC} + constant :: FT + inverse_T :: FT + inverse_T_sqrt_S :: FT + inverse_T_S :: FT + inverse_T_sqrt_S³ :: FT + inverse_T_S² :: FT + sqrt_S :: FT + S :: FT + log_T :: FT + log_T_sqrt_S :: FT + S_log_T :: FT + T_sqrt_S :: FT + + pressure_correction :: PC + + function KB(FT = Float64; + constant = 148.0248, + inverse_T = -8966.90, + inverse_T_sqrt_S = -2890.53, + inverse_T_S = -77.942, + inverse_T_sqrt_S³ = 1.728, + inverse_T_S² = -0.0996, + sqrt_S = 137.1942, + S = 1.62142, + log_T = -24.4344, + log_T_sqrt_S = -25.085, + S_log_T = -0.2474, + T_sqrt_S = 0.053105, + pressure_correction::PC = PressureCorrection(FT; a₀=-29.48, a₁=0.1622, a₂=-0.0026080, b₀=-0.00284, b₁=0.0)) where PC + + return new{FT, PC}(constant, inverse_T, + inverse_T_sqrt_S, inverse_T_S, + inverse_T_sqrt_S³, inverse_T_S², + sqrt_S, S, + log_T, log_T_sqrt_S, + S_log_T, T_sqrt_S, + pressure_correction) + end end @inline (c::KB)(T, S; P = nothing) = @@ -224,17 +276,32 @@ Parameterisation for water dissociation equilibrium constant. Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). """ -@kwdef struct KW{FT, PC} - constant :: FT = 148.9652 - inverse_T :: FT = -13847.26 - log_T :: FT = -23.6521 - sqrt_S :: FT = -5.977 - inverse_T_sqrt_S :: FT = 118.67 - log_T_sqrt_S :: FT = 1.0495 - S :: FT = -0.01615 - - pressure_correction :: PC = - PressureCorrection(; a₀=-20.02, a₁=0.1119, a₂=-0.001409, b₀=-0.00513, b₁=0.0000794) +struct KW{FT, PC} + constant :: FT + inverse_T :: FT + log_T :: FT + sqrt_S :: FT + inverse_T_sqrt_S :: FT + log_T_sqrt_S :: FT + S :: FT + + pressure_correction :: PC + + function KW(FT = Float64; + constant = 148.9652, + inverse_T = -13847.26, + log_T = -23.6521, + sqrt_S = -5.977, + inverse_T_sqrt_S = 118.67, + log_T_sqrt_S = 1.0495, + S = -0.01615, + pressure_correction::PC = PressureCorrection(FT; a₀=-20.02, a₁=0.1119, a₂=-0.001409, b₀=-0.00513, b₁=0.0000794)) where PC + + return new{FT, PC}(constant, inverse_T, + log_T, sqrt_S, + inverse_T_sqrt_S, log_T_sqrt_S, S, + pressure_correction) + end end @inline (c::KW)(T, S; P = nothing) = @@ -296,24 +363,48 @@ Parameterisation for bisulfate dissociation equilibrium constant. Default values from Dickson (1990, Chem. Thermodyn., 22, 113–127). """ -@kwdef struct KS{IS, FT, PC} - ionic_strength :: IS = IonicStrength() +struct KS{IS, FT, PC} + ionic_strength :: IS - constant :: FT = 141.328 - inverse_T :: FT = -4276.1 - log_T :: FT = -23.093 - sqrt_Is :: FT = 324.57 - inverse_T_sqrt_Is :: FT = -13856.0 - log_T_sqrt_Is :: FT = -47.986 - Is :: FT = -771.54 - inverse_T_Is :: FT = 35474.0 - log_T_Is :: FT = 114.723 - inverse_T_sqrt_Is³ :: FT = -2698.0 - inverse_T_Is² :: FT = 1776.0 - log_S :: FT = -0.001005 - - pressure_correction :: PC = - PressureCorrection(; a₀=-18.03, a₁=0.0466, a₂=0.000316, b₀=-0.00453, b₁=0.00009) + constant :: FT + inverse_T :: FT + log_T :: FT + sqrt_Is :: FT + inverse_T_sqrt_Is :: FT + log_T_sqrt_Is :: FT + Is :: FT + inverse_T_Is :: FT + log_T_Is :: FT + inverse_T_sqrt_Is³ :: FT + inverse_T_Is² :: FT + log_S :: FT + + pressure_correction :: PC + + function KS(FT = Float64; + ionic_strength::IS = IonicStrength{FT}(), + constant = 141.328, + inverse_T = -4276.1, + log_T = -23.093, + sqrt_Is = 324.57, + inverse_T_sqrt_Is = -13856.0, + log_T_sqrt_Is = -47.986, + Is = -771.54, + inverse_T_Is = 35474.0, + log_T_Is = 114.723, + inverse_T_sqrt_Is³ = -2698.0, + inverse_T_Is² = 1776.0, + log_S = -0.001005, + pressure_correction::PC = + PressureCorrection(FT; a₀=-18.03, a₁=0.0466, a₂=0.000316, b₀=-0.00453, b₁=0.00009)) where {IS, PC} + + return new{IS, FT, PC}(ionic_strength, + constant, inverse_T, log_T, + sqrt_Is, inverse_T_sqrt_Is, log_T_sqrt_Is, + Is, inverse_T_Is, log_T_Is, inverse_T_sqrt_Is³, + inverse_T_Is², log_S, + pressure_correction) + end end @inline (c::KS)(T, S, Is = c.ionic_strength(S); P = nothing) = @@ -357,18 +448,34 @@ Parameterisation for hydrogen fluoride dissociation equilibrium constant. Default values from Perez and Fraga (1987, Mar. Chem., 21, 161–168). """ -@kwdef struct KF{IS, KS, FT, PC} - ionic_strength :: IS = IonicStrength() - sulfate_constant :: KS = KS(; ionic_strength) +struct KF{IS, KS, FT, PC} + ionic_strength :: IS + sulfate_constant :: KS - constant :: FT = -9.68 - inverse_T :: FT = 874.0 - sqrt_S :: FT = 0.111 - log_S :: FT = 0.0 - log_S_KS :: FT = 0.0 - - pressure_correction :: PC = - PressureCorrection(; a₀=-9.78, a₁=-0.0090, a₂=-0.000942, b₀=-0.00391, b₁=0.000054) + constant :: FT + inverse_T :: FT + sqrt_S :: FT + log_S :: FT + log_S_KS :: FT + + pressure_correction :: PC + + function KF(FT = Float64; + ionic_strength::IS = IonicStrength{FT}(), + sulfate_constant::KS = KS(FT; ionic_strength), + + constant = -9.68, + inverse_T = 874.0, + sqrt_S = 0.111, + log_S = 0.0, + log_S_KS = 0.0, + + pressure_correction::PC = + PressureCorrection(FT; a₀=-9.78, a₁=-0.0090, a₂=-0.000942, b₀=-0.00391, b₁=0.000054)) where {IS, KS, PC} + + return new{IS, KS, FT, PC}(ionic_strength, sulfate_constant, constant, inverse_T, sqrt_S, log_S, log_S_KS, pressure_correction) + end + end @inline (c::KF)(T, S, Is = c.ionic_strength(S), KS = c.sulfate_constant(T, S, Is); P = nothing) = @@ -448,23 +555,24 @@ Instance of `KP` returning the first phosphocic acid equilibrium constant. Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). """ -KP1(; constant = 115.525, - inverse_T = -4576.752, - log_T = - 18.453, - sqrt_S = 0.69171, - inverse_T_sqrt_S = -106.736, - S = -0.01844, - inverse_T_S = -0.65643, - pressure_correction = - PressureCorrection(; a₀=-14.51, a₁=0.1211, a₂=-0.000321, b₀=-0.00267, b₁=0.0000427)) = - KP(constant, - inverse_T, - log_T, - sqrt_S, - inverse_T_sqrt_S, - S, - inverse_T_S, - pressure_correction) +KP1(FT = Float64; + constant = 115.525, + inverse_T = -4576.752, + log_T = - 18.453, + sqrt_S = 0.69171, + inverse_T_sqrt_S = -106.736, + S = -0.01844, + inverse_T_S = -0.65643, + pressure_correction::PC = + PressureCorrection(FT; a₀=-14.51, a₁=0.1211, a₂=-0.000321, b₀=-0.00267, b₁=0.0000427)) where PC = + KP{FT, PC}(constant, + inverse_T, + log_T, + sqrt_S, + inverse_T_sqrt_S, + S, + inverse_T_S, + pressure_correction) """ KP2(; constant = 172.0883, @@ -485,23 +593,24 @@ Instance of `KP` returning the second phosphocic acid equilibrium constant. Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). """ -KP2(; constant = 172.0883, - inverse_T = -8814.715, - log_T = -27.927, - sqrt_S = 1.3566, - inverse_T_sqrt_S = -160.340, - S = -0.05778, - inverse_T_S = 0.37335, - pressure_correction = - PressureCorrection(; a₀=-23.12, a₁=0.1758, a₂=-0.002647, b₀=-0.00515, b₁=0.00009)) = - KP(constant, - inverse_T, - log_T, - sqrt_S, - inverse_T_sqrt_S, - S, - inverse_T_S, - pressure_correction) +KP2(FT = Float64; + constant = 172.0883, + inverse_T = -8814.715, + log_T = -27.927, + sqrt_S = 1.3566, + inverse_T_sqrt_S = -160.340, + S = -0.05778, + inverse_T_S = 0.37335, + pressure_correction::PC = + PressureCorrection(FT; a₀=-23.12, a₁=0.1758, a₂=-0.002647, b₀=-0.00515, b₁=0.00009)) where PC = + KP{FT, PC}(constant, + inverse_T, + log_T, + sqrt_S, + inverse_T_sqrt_S, + S, + inverse_T_S, + pressure_correction) """ KP3(; constant = -18.141, @@ -522,23 +631,24 @@ Instance of `KP` returning the third phosphocic acid equilibrium constant. Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). """ -KP3(; constant = - 18.141, - inverse_T = -3070.75, - log_T = 0.0, - sqrt_S = 2.81197, - inverse_T_sqrt_S = 17.27039, - S = -0.09984, - inverse_T_S = -44.99486, - pressure_correction = - PressureCorrection(; a₀=-26.57, a₁=0.2020, a₂=-0.0030420, b₀=-0.00408, b₁=0.0000714)) = - KP(constant, - inverse_T, - log_T, - sqrt_S, - inverse_T_sqrt_S, - S, - inverse_T_S, - pressure_correction) +KP3(FT = Float64; + constant = - 18.141, + inverse_T = -3070.75, + log_T = 0.0, + sqrt_S = 2.81197, + inverse_T_sqrt_S = 17.27039, + S = -0.09984, + inverse_T_S = -44.99486, + pressure_correction::PC = + PressureCorrection(FT; a₀=-26.57, a₁=0.2020, a₂=-0.0030420, b₀=-0.00408, b₁=0.0000714)) where PC = + KP{FT, PC}(constant, + inverse_T, + log_T, + sqrt_S, + inverse_T_sqrt_S, + S, + inverse_T_S, + pressure_correction) """ KSi(; ionic_strength = IonicStrength(), @@ -561,19 +671,36 @@ Parameterisation for silicic acid dissociation equilibrium constant. Default values from Millero (1995, Geochim. Cosmochim. Acta, 59, 661–677). """ -@kwdef struct KSi{IS, FT} - ionic_strength :: IS = IonicStrength() - - constant :: FT = 117.385 - inverse_T :: FT = -8904.2 - log_T :: FT = -19.334 - sqrt_Is :: FT = 3.5913 - inverse_T_sqrt_Is :: FT = -458.79 - Is :: FT = -1.5998 - inverse_T_Is :: FT = 188.74 - Is² :: FT = 0.07871 - inverse_T_Is² :: FT = -12.1652 - log_S :: FT = -0.001005 +struct KSi{IS, FT} + ionic_strength :: IS + + constant :: FT + inverse_T :: FT + log_T :: FT + sqrt_Is :: FT + inverse_T_sqrt_Is :: FT + Is :: FT + inverse_T_Is :: FT + Is² :: FT + inverse_T_Is² :: FT + log_S :: FT + + function KSi(FT = Float64; + ionic_strength::IS = IonicStrength{FT}(), + constant = 117.385, + inverse_T = -8904.2, + log_T = -19.334, + sqrt_Is = 3.5913, + inverse_T_sqrt_Is = -458.79, + Is = -1.5998, + inverse_T_Is = 188.74, + Is² = 0.07871, + inverse_T_Is² = -12.1652, + log_S = -0.001005) where IS + + return new{IS, FT}(ionic_strength, constant, inverse_T, log_T, sqrt_Is, inverse_T_sqrt_Is, Is, + inverse_T_Is, Is², inverse_T_Is², log_S) + end end @inline (c::KSi)(T, S, Is = c.ionic_strength(S); P = nothing) = @@ -659,17 +786,18 @@ Instance of `KSP` returning calcite solubility. Default values from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341). """ -KSP_calcite(; therm_constant = -171.9065, - therm_T = -0.077993, - therm_inverse_T = 2839.319, - therm_log_T = 71.595, - sea_sqrt_S = -0.77712, - sea_T_sqrt_S = 0.0028426, - sea_inverse_T_sqrt_S = 178.34, - sea_S = -0.07711, - sea_S_sqrt_S³ = 0.0041249, - pressure_correction = - PressureCorrection(; a₀=-48.76, a₁=0.5304, a₂=-0.0, b₀=-0.01176, b₁=0.0003692)) = +KSP_calcite(FT = Float64; + therm_constant = -171.9065, + therm_T = -0.077993, + therm_inverse_T = 2839.319, + therm_log_T = 71.595, + sea_sqrt_S = -0.77712, + sea_T_sqrt_S = 0.0028426, + sea_inverse_T_sqrt_S = 178.34, + sea_S = -0.07711, + sea_S_sqrt_S³ = 0.0041249, + pressure_correction = + PressureCorrection(FT; a₀=-48.76, a₁=0.5304, a₂=-0.0, b₀=-0.01176, b₁=0.0003692)) = KSP(therm_constant, therm_T, therm_inverse_T, @@ -698,17 +826,18 @@ Instance of `KSP` returning calcite solubility. Default values from Millero, F. J. (2007, Chemical Reviews, 107(2), 308–341). """ -KSP_aragonite(; therm_constant = -171.945, - therm_T = -0.077993, - therm_inverse_T = 2903.293, - therm_log_T = 71.595, - sea_sqrt_S = -0.068393, - sea_T_sqrt_S = 0.0017276, - sea_inverse_T_sqrt_S = 88.135, - sea_S = -0.10018, - sea_S_sqrt_S³ = 0.0059415, - pressure_correction = - PressureCorrection(; a₀=-45.96, a₁=0.5304, a₂=-0.0, b₀=-0.01176, b₁=0.0003692)) = +KSP_aragonite(FT = Float64; + therm_constant = -171.945, + therm_T = -0.077993, + therm_inverse_T = 2903.293, + therm_log_T = 71.595, + sea_sqrt_S = -0.068393, + sea_T_sqrt_S = 0.0017276, + sea_inverse_T_sqrt_S = 88.135, + sea_S = -0.10018, + sea_S_sqrt_S³ = 0.0059415, + pressure_correction = + PressureCorrection(FT; a₀=-45.96, a₁=0.5304, a₂=-0.0, b₀=-0.01176, b₁=0.0003692)) = KSP(therm_constant, therm_T, therm_inverse_T, diff --git a/src/Models/GasExchange/GasExchange.jl b/src/Models/GasExchange/GasExchange.jl index 76a4b98a7..c4e4fc36e 100644 --- a/src/Models/GasExchange/GasExchange.jl +++ b/src/Models/GasExchange/GasExchange.jl @@ -65,14 +65,15 @@ or a `CarbonDioxideConcentration` which diagnoses the partial pressure of CO₂ `transfer_velocity` should be a function of the form `k(u₁₀, T)`. """ -function GasExchangeBoundaryCondition(; water_concentration, - air_concentration, - transfer_velocity, - wind_speed, - discrete_form = false) +function GasExchangeBoundaryCondition(FT = Float64; + water_concentration, + air_concentration, + transfer_velocity, + wind_speed, + discrete_form = false) - wind_speed = normalise_surface_function(wind_speed; discrete_form) - air_concentration = normalise_surface_function(air_concentration; discrete_form) + wind_speed = normalise_surface_function(wind_speed; discrete_form, FT) + air_concentration = normalise_surface_function(air_concentration; discrete_form, FT) exchange_function = GasExchange(wind_speed, transfer_velocity, water_concentration, air_concentration) @@ -80,13 +81,14 @@ function GasExchangeBoundaryCondition(; water_concentration, end """ - CarbonDioxideGasExchangeBoundaryCondition(; carbon_chemistry = CarbonChemistry(), - transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = CarbonDioxidePolynomialSchmidtNumber()), - air_concentration = 413, # ppmv - wind_speed = 2, - water_concentration = nothing, - silicate_and_phosphate_names = nothing, - kwargs...) + CarbonDioxideGasExchangeBoundaryCondition(FT = Float64; + carbon_chemistry = CarbonChemistry(FT), + transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = CarbonDioxidePolynomialSchmidtNumber(FT)), + air_concentration = 413, # ppmv + wind_speed = 2, + water_concentration = nothing, + silicate_and_phosphate_names = nothing, + kwargs...) Returns a `FluxBoundaryCondition` for the gas exchange between carbon dioxide dissolved in the water specified by the `carbon_chemisty` model, and `air_concentration` with `transfer_velocity` (see @@ -99,13 +101,14 @@ and phosphate tracers, or a `NamedTuple` of values for the `carbon_chemistry` m Note: The model always requires `T`, `S`, `DIC`, and `Alk` to be present in the model. """ -function CarbonDioxideGasExchangeBoundaryCondition(; carbon_chemistry = CarbonChemistry(), - transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = CarbonDioxidePolynomialSchmidtNumber()), - air_concentration = 413, # ppmv - wind_speed = 2, - water_concentration = nothing, - silicate_and_phosphate_names = nothing, - kwargs...) +function CarbonDioxideGasExchangeBoundaryCondition(FT = Float64; + carbon_chemistry = CarbonChemistry(FT), + transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = CarbonDioxidePolynomialSchmidtNumber(FT)), + air_concentration = 413, # ppmv + wind_speed = 2, + water_concentration = nothing, + silicate_and_phosphate_names = nothing, + kwargs...) if isnothing(water_concentration) water_concentration = CarbonDioxideConcentration(; carbon_chemistry, silicate_and_phosphate_names) @@ -113,15 +116,16 @@ function CarbonDioxideGasExchangeBoundaryCondition(; carbon_chemistry = CarbonCh @warn "Make sure that the `carbon_chemistry` $(carbon_chemistry) is the same as that in `water_concentration` $(water_concentration) (or set it to `nothing`)" end - return GasExchangeBoundaryCondition(; water_concentration, air_concentration, transfer_velocity, wind_speed, kwargs...) + return GasExchangeBoundaryCondition(FT; water_concentration, air_concentration, transfer_velocity, wind_speed, kwargs...) end """ - OxygenGasExchangeBoundaryCondition(; transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = OxygenPolynomialSchmidtNumber()), - water_concentration = OxygenConcentration(), - air_concentration = 9352.7, # mmolO₂/m³ - wind_speed = 2, - kwagrs...) + OxygenGasExchangeBoundaryCondition(FT = Float64; + transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = OxygenPolynomialSchmidtNumber(FT)), + water_concentration = OxygenConcentration(), + air_concentration = 9352.7, # mmolO₂/m³ + wind_speed = 2, + kwagrs...) Returns a `FluxBoundaryCondition` for the gas exchange between oxygen dissolved in the water specified by the the `OxygenConcentration` in the base model, and `air_concentration` with `transfer_velocity` @@ -129,11 +133,12 @@ specified by the the `OxygenConcentration` in the base model, and `air_concentra `kwargs` are passed on to `GasExchangeBoundaryCondition`. """ -OxygenGasExchangeBoundaryCondition(; transfer_velocity = SchmidtScaledTransferVelocity(; schmidt_number = OxygenPolynomialSchmidtNumber()), - water_concentration = OxygenConcentration(), - air_concentration = PartiallySolubleGas(9352.7, OxygenSolubility()), # mmolO₂/m³ - wind_speed = 2, - kwargs...) = - GasExchangeBoundaryCondition(; water_concentration, air_concentration, transfer_velocity, wind_speed, kwargs...) +OxygenGasExchangeBoundaryCondition(FT = Float64; + transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = OxygenPolynomialSchmidtNumber(FT)), + water_concentration = OxygenConcentration(), + air_concentration = PartiallySolubleGas(FT; air_concentration = 9352.7, solubility = OxygenSolubility(FT)), # mmolO₂/m³ + wind_speed = 2, + kwargs...) = + GasExchangeBoundaryCondition(FT; water_concentration, air_concentration, transfer_velocity, wind_speed, kwargs...) end # module \ No newline at end of file diff --git a/src/Models/GasExchange/gas_solubility.jl b/src/Models/GasExchange/gas_solubility.jl index 4b2bd526d..a213468a8 100644 --- a/src/Models/GasExchange/gas_solubility.jl +++ b/src/Models/GasExchange/gas_solubility.jl @@ -7,8 +7,18 @@ where ``\alpha`` is the Ostwald solubility coeffieient and ``C_a`` is the concen struct PartiallySolubleGas{AC, S} air_concentration :: AC solubility :: S -end + function PartiallySolubleGas(FT = Float64; + air_concentration, + solubility :: S) where S + + air_concentration = normalise_surface_function(air_concentration; FT) + + AC = typeof(air_concentration) + + return new{AC, S}(air_concentration, solubility) + end +end @inline surface_value(gs::PartiallySolubleGas, i, j, grid, clock, model_fields) = surface_value(gs.air_concentration, i, j, grid, clock) * surface_value(gs.solubility, i, j, grid, clock, model_fields) @@ -38,5 +48,5 @@ function surface_value(sol::Wanninkhof92Solubility, i, j, grid, clock, model_fie return β / Tk end -OxygenSolubility(; A1 = -58.3877, A2 = 85.8079, A3 = 23.8439, B1 = -0.034892, B2 = 0.015578, B3 = -0.0019387) = - Wanninkhof92Solubility(A1, A2, A3, B1, B2, B3) \ No newline at end of file +OxygenSolubility(FT = Float64; A1 = -58.3877, A2 = 85.8079, A3 = 23.8439, B1 = -0.034892, B2 = 0.015578, B3 = -0.0019387) = + Wanninkhof92Solubility{FT}(A1, A2, A3, B1, B2, B3) \ No newline at end of file diff --git a/src/Models/GasExchange/gas_transfer_velocity.jl b/src/Models/GasExchange/gas_transfer_velocity.jl index e97db5a97..2f42e21d8 100644 --- a/src/Models/GasExchange/gas_transfer_velocity.jl +++ b/src/Models/GasExchange/gas_transfer_velocity.jl @@ -39,95 +39,95 @@ show(io::IO, k::SchmidtScaledTransferVelocity{KB, SC}) where {KB, SC} = " k₆₆₀(u₁₀) : $(nameof(KB))") """ - Wanninkhof99(scale_factor = 0.0283 / hour / 100) + Wanninkhof99(FT = Float64; scale_factor = 0.0283 / hour / 100) Cubic k₆₆₀ parameterisation of Wanninkhof & McGillis (1999) suitable for short term, in situ wind products. """ -Wanninkhof99(scale_factor = 0.0283 / hour / 100) = PolynomialParameterisation{3}((0, 0, 0, scale_factor)) +Wanninkhof99(FT = Float64; scale_factor = 0.0283 / hour / 100) = PolynomialParameterisation{3}(FT; coefficients = (0, 0, 0, scale_factor)) """ - Ho06(scale_factor = 0.266 / hour / 100) + Ho06(FT = Float64; scale_factor = 0.266 / hour / 100) Quadratic k₆₆₀ parameterisation of Ho et al. (2006) suitable for the QuickSCAT satellite and short-term steady wind product. """ -Ho06(scale_factor = 0.266 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) +Ho06(FT = Float64; scale_factor = 0.266 / hour / 100) = PolynomialParameterisation{2}(FT; coefficients = (0, 0, scale_factor)) """ - Nightingale00(linear = 0.333 / hour / 100, quadratic = 0.222 / hour / 100) + Nightingale00(FT = Float64; linear = 0.333 / hour / 100, quadratic = 0.222 / hour / 100) Cubic k₆₆₀ parameterisation of Nightingale et al. (2000) suitable for short term, in situ wind products (?). """ -Nightingale00(linear = 0.333 / hour / 100, quadratic = 0.222 / hour / 100) = - PolynomialParameterisation{2}((0, linear, quadratic)) +Nightingale00(FT = Float64; linear = 0.333 / hour / 100, quadratic = 0.222 / hour / 100) = + PolynomialParameterisation{2}(FT; coefficients = (0, linear, quadratic)) """ - McGillis01(constant = 3.3 / hour / 100, cubic = 0.026 / hour / 100) + McGillis01(FT = Float64; constant = 3.3 / hour / 100, cubic = 0.026 / hour / 100) Cubic k₆₆₀ parameterisation of McGillis et al. (2001) suitable for short term, in situ wind products. """ -McGillis01(constant = 3.3 / hour / 100, cubic = 0.026 / hour / 100) = - PolynomialParameterisation{3}((constant, 0, 0, cubic)) +McGillis01(FT = Float64; constant = 3.3 / hour / 100, cubic = 0.026 / hour / 100) = + PolynomialParameterisation{3}(FT; coefficients = (constant, 0, 0, cubic)) """ - Sweeny07(scale_factor = 0.27 / hour / 100) + Sweeny07(FT = Float64; scale_factor = 0.27 / hour / 100) Quadratic k₆₆₀ parameterisation of Sweeny et al. (2007) suitable for the NCEP/NCAR reanalysis 1 product """ -Sweeny07(scale_factor = 0.27 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) +Sweeny07(FT = Float64; scale_factor = 0.27 / hour / 100) = PolynomialParameterisation{2}(FT; coefficients = (0, 0, scale_factor)) """ - Wanninkhof09(constant = 3 / hour / 100, linear = 0.1 / hour / 100, quadratic = 0.064 / hour / 100, cubic = 0.011 / hour / 100) + Wanninkhof09(FT = Float64; constant = 3 / hour / 100, linear = 0.1 / hour / 100, quadratic = 0.064 / hour / 100, cubic = 0.011 / hour / 100) Cubic k₆₆₀ parameterisation of Wanninkhof et al (2009) suitable for the Cross-Calibrated Multi-Platform (CCMP) Winds product """ -Wanninkhof09(constant = 3 / hour / 100, linear = 0.1 / hour / 100, quadratic = 0.064 / hour / 100, cubic = 0.011 / hour / 100) = - PolynomialParameterisation{3}((constant, linear, quadratic, cubic)) +Wanninkhof09(FT = Float64; constant = 3 / hour / 100, linear = 0.1 / hour / 100, quadratic = 0.064 / hour / 100, cubic = 0.011 / hour / 100) = + PolynomialParameterisation{3}(FT; coefficients = (constant, linear, quadratic, cubic)) """ - Wanninkhof14(scale_factor = 0.251 / hour / 100) + Wanninkhof14(FT = Float64; scale_factor = 0.251 / hour / 100) Quadratic k₆₆₀ parameterisation of Wanninkhof et al (2014) suitable for the Cross-Calibrated Multi-Platform (CCMP) Winds product """ -Wanninkhof14(scale_factor = 0.251 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) +Wanninkhof14(FT = Float64; scale_factor = 0.251 / hour / 100) = PolynomialParameterisation{2}(FT; coefficients = (0, 0, scale_factor)) """ - ERA5(scale_factor = 0.270875 / hour / 100) + ERA5(FT = Float64; scale_factor = 0.270875 / hour / 100) Quadratic k₆₆₀ parameterisation calibrated to give 16.5 cm/hr global average (reccomended Naegler, 2009) for the ERA5 wind product by SeaFlux/Luke Gregor et al. (2023). """ -ERA5(scale_factor = 0.270875 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) +ERA5(FT = Float64; scale_factor = 0.270875 / hour / 100) = PolynomialParameterisation{2}(FT; coefficients = (0, 0, scale_factor)) """ - JRA55(scale_factor = 0.2601975 / hour / 100) + JRA55(FT = Float64; scale_factor = 0.2601975 / hour / 100) Quadratic k₆₆₀ parameterisation calibrated to give 16.5 cm/hr global average (reccomended Naegler, 2009) for the JRA55 wind product by SeaFlux/Luke Gregor et al. (2023). """ -JRA55(scale_factor = 0.2601975 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) +JRA55(FT = Float64; scale_factor = 0.2601975 / hour / 100) = PolynomialParameterisation{2}(FT; coefficients = (0, 0, scale_factor)) """ - NCEP1(scale_factor = 0.2866424 / hour / 100) + NCEP1(FT = Float64; scale_factor = 0.2866424 / hour / 100) Quadratic k₆₆₀ parameterisation calibrated to give 16.5 cm/hr global average (reccomended Naegler, 2009) for the NCEP1 wind product by SeaFlux/Luke Gregor et al. (2023). """ -NCEP1(scale_factor = 0.2866424 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) +NCEP1(FT = Float64; scale_factor = 0.2866424 / hour / 100) = PolynomialParameterisation{2}(FT; coefficients = (0, 0, scale_factor)) """ - CCMP2(scale_factor = 0.256789 / hour / 100) + CCMP2(FT = Float64; scale_factor = 0.256789 / hour / 100) Quadratic k₆₆₀ parameterisation calibrated to give 16.5 cm/hr global average (reccomended Naegler, 2009) for the CCMP2 wind product by SeaFlux/Luke Gregor et al. (2023). """ -CCMP2(scale_factor = 0.256789 / hour / 100) = PolynomialParameterisation{2}((0, 0, scale_factor)) +CCMP2(FT = Float64; scale_factor = 0.256789 / hour / 100) = PolynomialParameterisation{2}(FT; coefficients = (0, 0, scale_factor)) end # module \ No newline at end of file diff --git a/src/Models/GasExchange/generic_parameterisations.jl b/src/Models/GasExchange/generic_parameterisations.jl index 7d4c66f0e..567b0dd3e 100644 --- a/src/Models/GasExchange/generic_parameterisations.jl +++ b/src/Models/GasExchange/generic_parameterisations.jl @@ -2,10 +2,12 @@ struct PolynomialParameterisation{N, C} coefficients :: C end -function PolynomialParameterisation{N}(coefficients) where N +function PolynomialParameterisation{N}(FT = Float64; coefficients) where N length(coefficients) == N + 1 || throw(ArgumentError("You must provide N+1 coefficients for an order N polynomial")) + coefficients = convert.(FT, coefficients) + return PolynomialParameterisation{N, typeof(coefficients)}(coefficients) end diff --git a/src/Models/GasExchange/schmidt_number.jl b/src/Models/GasExchange/schmidt_number.jl index 6cf14fcdd..16cd7d20c 100644 --- a/src/Models/GasExchange/schmidt_number.jl +++ b/src/Models/GasExchange/schmidt_number.jl @@ -1,16 +1,16 @@ """ - CarbonDioxidePolynomialSchmidtNumber(; a = 2116.8, b = -136.25, c = 4.7353, d = -0.092307, e = 0.0007555) + CarbonDioxidePolynomialSchmidtNumber(FT = Float64; a = 2116.8, b = -136.25, c = 4.7353, d = -0.092307, e = 0.0007555) Schmidt number parameterisation Wanninkhof, 2014 for sea water """ -CarbonDioxidePolynomialSchmidtNumber(; a = 2116.8, b = -136.25, c = 4.7353, d = -0.092307, e = 0.0007555) = - PolynomialParameterisation{4}((a, b, c, d, e)) +CarbonDioxidePolynomialSchmidtNumber(FT = Float64; a = 2116.8, b = -136.25, c = 4.7353, d = -0.092307, e = 0.0007555) = + PolynomialParameterisation{4}(FT; coefficients = (a, b, c, d, e)) """ - OxygenPolynomialSchmidtNumber(; a = 1953.4, b = - 128.0, c = 3.9918, d = -0.050091) + OxygenPolynomialSchmidtNumber(FT = Float64; a = 1953.4, b = - 128.0, c = 3.9918, d = -0.050091) Schmidt number parameterisation Wanninkhof, 2014 for sea water """ -OxygenPolynomialSchmidtNumber(; a = 1920.4, b = -135.6, c = 5.2122, d = -0.10939, e = 0.00093777) = - PolynomialParameterisation{4}((a, b, c, d, e)) \ No newline at end of file +OxygenPolynomialSchmidtNumber(FT = Float64; a = 1920.4, b = -135.6, c = 5.2122, d = -0.10939, e = 0.00093777) = + PolynomialParameterisation{4}(FT; coefficients = (a, b, c, d, e)) \ No newline at end of file diff --git a/src/Models/GasExchange/surface_values.jl b/src/Models/GasExchange/surface_values.jl index 7bc950239..55b0a6e6a 100644 --- a/src/Models/GasExchange/surface_values.jl +++ b/src/Models/GasExchange/surface_values.jl @@ -29,4 +29,6 @@ end # fallback normalise_surface_function(f; kwargs...) = f +normalise_surface_function(f::Number; FT, kwargs...) = convert(FT, f) + normalise_surface_function(f::Function; discrete_form = false) = discrete_form ? DiscreteSurfaceFuncton(f) : ContinuousSurfaceFunction(f) \ No newline at end of file diff --git a/src/Models/Individuals/SugarKelp/SugarKelp.jl b/src/Models/Individuals/SugarKelp/SugarKelp.jl index dfced410f..7e6ad2a3e 100644 --- a/src/Models/Individuals/SugarKelp/SugarKelp.jl +++ b/src/Models/Individuals/SugarKelp/SugarKelp.jl @@ -32,6 +32,9 @@ import OceanBioME.Particles: required_particle_fields, required_tracers, coupled SugarKelp Defines the parameters for `SugarKelp` biogeochemistry. + +TODO: make a function to make it so this can be automatically converted into different float types as that wouldn't work at the moment +but there is too many parameters """ @kwdef struct SugarKelp{FT, TL} temperature_limit :: TL = LinearOptimalTemperatureRange() # TODO: split up more parameterisations like this diff --git a/src/Particles/Particles.jl b/src/Particles/Particles.jl index e1d733c8d..a67a7f6f2 100644 --- a/src/Particles/Particles.jl +++ b/src/Particles/Particles.jl @@ -10,6 +10,7 @@ using Oceananigans: NonhydrostaticModel, HydrostaticFreeSurfaceModel using OceanBioME: DiscreteBiogeochemistry, ContinuousBiogeochemistry using Oceananigans.Architectures: architecture, on_architecture +using Oceananigans.Grids: AbstractGrid import Oceananigans.Architectures: architecture import Oceananigans.Biogeochemistry: update_tendencies! @@ -75,22 +76,22 @@ Particles can also have a `scalefactor` which scales their tracer interaction (e.g. to mimic the particle representing multiple particles). """ function BiogeochemicalParticles(number; - grid, + grid::AbstractGrid{FT}, biogeochemistry, advection = LagrangianAdvection(), timestepper = ForwardEuler, field_interpolation = NearestPoint(), - scalefactors = ones(number)) + scalefactors = ones(number)) where FT arch = architecture(grid) particle_fields = required_particle_fields(biogeochemistry) - fields = NamedTuple{particle_fields}(ntuple(n->on_architecture(arch, zeros(number)), Val(length(particle_fields)))) + fields = NamedTuple{particle_fields}(ntuple(n->on_architecture(arch, zeros(FT, number)), Val(length(particle_fields)))) - x = on_architecture(arch, zeros(number)) - y = on_architecture(arch, zeros(number)) - z = on_architecture(arch, zeros(number)) + x = on_architecture(arch, zeros(FT, number)) + y = on_architecture(arch, zeros(FT, number)) + z = on_architecture(arch, zeros(FT, number)) timestepper = timestepper(particle_fields, number, arch) From e82625fd615e349c233267a7f3a3d27fee0c3015 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 14 Oct 2024 17:52:34 +0100 Subject: [PATCH 3/8] issue --- src/Models/GasExchange/surface_values.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Models/GasExchange/surface_values.jl b/src/Models/GasExchange/surface_values.jl index 55b0a6e6a..5bf72d36a 100644 --- a/src/Models/GasExchange/surface_values.jl +++ b/src/Models/GasExchange/surface_values.jl @@ -31,4 +31,4 @@ normalise_surface_function(f; kwargs...) = f normalise_surface_function(f::Number; FT, kwargs...) = convert(FT, f) -normalise_surface_function(f::Function; discrete_form = false) = discrete_form ? DiscreteSurfaceFuncton(f) : ContinuousSurfaceFunction(f) \ No newline at end of file +normalise_surface_function(f::Function; discrete_form = false, kwargs...) = discrete_form ? DiscreteSurfaceFuncton(f) : ContinuousSurfaceFunction(f) \ No newline at end of file From 68f84a8063af97d7460c8ca6d16e8f1a8296d05f Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 15 Oct 2024 16:27:20 +0100 Subject: [PATCH 4/8] bump patch --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5c4689eac..6fcea1efd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OceanBioME" uuid = "a49af516-9db8-4be4-be45-1dad61c5a376" authors = ["Jago Strong-Wright and contributors"] -version = "0.13.0" +version = "0.13.1" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" From 3e0de9fd8838d155bac39dbb7c8f1a5521fec071 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 16 Oct 2024 11:21:35 +0100 Subject: [PATCH 5/8] Removed positional arg --- .../AdvectedPopulations/PISCES/PISCES.jl | 118 +++++++++--------- 1 file changed, 59 insertions(+), 59 deletions(-) diff --git a/src/Models/AdvectedPopulations/PISCES/PISCES.jl b/src/Models/AdvectedPopulations/PISCES/PISCES.jl index e78134bc3..656fb3a98 100644 --- a/src/Models/AdvectedPopulations/PISCES/PISCES.jl +++ b/src/Models/AdvectedPopulations/PISCES/PISCES.jl @@ -284,65 +284,65 @@ the classes to a single `phytoplankton` if more classes are required (see `OceanBioME.Models.PISCESModel` docstring). Similarly, if a more generic `particulate_organic_matter` was desired a way to specify arbitary tracers for arguments would be required. """ -function PISCES(FT = Float64; - grid, - phytoplankton = MixedMondoNanoAndDiatoms(FT), - zooplankton = MicroAndMesoZooplankton(FT), - dissolved_organic_matter = DissolvedOrganicCarbon(FT), - particulate_organic_matter = TwoCompartementCarbonIronParticles(FT), - - nitrogen = NitrateAmmonia{FT}(), - iron = SimpleIron{FT}(), - silicate = Silicate(), - oxygen = Oxygen{FT}(), - phosphate = Phosphate(), - - inorganic_carbon = InorganicCarbon(), - - # from Aumount 2005 rather than 2015 since it doesn't work the other way around - first_anoxia_thresehold = convert(FT, 6.0), - second_anoxia_thresehold = convert(FT, 1.0), - - nitrogen_redfield_ratio = convert(FT, 16/122), - phosphate_redfield_ratio = convert(FT, 1/122), - - mixed_layer_shear = convert(FT, 1.0), - background_shear = convert(FT, 0.01), - - latitude = PrescribedLatitude{FT}(45), - day_length = CBMDayLength(), - - mixed_layer_depth = Field{Center, Center, Nothing}(grid), - euphotic_depth = Field{Center, Center, Nothing}(grid), - - silicate_climatology = ConstantField(convert(FT, 7.5)), - - mean_mixed_layer_vertical_diffusivity = Field{Center, Center, Nothing}(grid), - mean_mixed_layer_light = Field{Center, Center, Nothing}(grid), - - carbon_chemistry = CarbonChemistry(FT), - calcite_saturation = CenterField(grid), - - surface_photosynthetically_active_radiation = default_surface_PAR, - - light_attenuation = - MultiBandPhotosyntheticallyActiveRadiation(; grid, - surface_PAR = surface_photosynthetically_active_radiation), - - sinking_speeds = (POC = convert(FT, 2/day), - # might be more efficient to just precompute this - GOC = Field(KernelFunctionOperation{Center, Center, Face}(DepthDependantSinkingSpeed{FT}(), - grid, - mixed_layer_depth, - euphotic_depth))), - open_bottom = true, - - scale_negatives = false, - invalid_fill_value = NaN, - - sediment = nothing, - particles = nothing, - modifiers = nothing) +function PISCES(; grid, + FT = eltype(grid), + phytoplankton = MixedMondoNanoAndDiatoms(FT), + zooplankton = MicroAndMesoZooplankton(FT), + dissolved_organic_matter = DissolvedOrganicCarbon(FT), + particulate_organic_matter = TwoCompartementCarbonIronParticles(FT), + + nitrogen = NitrateAmmonia{FT}(), + iron = SimpleIron{FT}(), + silicate = Silicate(), + oxygen = Oxygen{FT}(), + phosphate = Phosphate(), + + inorganic_carbon = InorganicCarbon(), + + # from Aumount 2005 rather than 2015 since it doesn't work the other way around + first_anoxia_thresehold = convert(FT, 6.0), + second_anoxia_thresehold = convert(FT, 1.0), + + nitrogen_redfield_ratio = convert(FT, 16/122), + phosphate_redfield_ratio = convert(FT, 1/122), + + mixed_layer_shear = convert(FT, 1.0), + background_shear = convert(FT, 0.01), + + latitude = PrescribedLatitude{FT}(45), + day_length = CBMDayLength(), + + mixed_layer_depth = Field{Center, Center, Nothing}(grid), + euphotic_depth = Field{Center, Center, Nothing}(grid), + + silicate_climatology = ConstantField(convert(FT, 7.5)), + + mean_mixed_layer_vertical_diffusivity = Field{Center, Center, Nothing}(grid), + mean_mixed_layer_light = Field{Center, Center, Nothing}(grid), + + carbon_chemistry = CarbonChemistry(FT), + calcite_saturation = CenterField(grid), + + surface_photosynthetically_active_radiation = default_surface_PAR, + + light_attenuation = + MultiBandPhotosyntheticallyActiveRadiation(; grid, + surface_PAR = surface_photosynthetically_active_radiation), + + sinking_speeds = (POC = convert(FT, 2/day), + # might be more efficient to just precompute this + GOC = Field(KernelFunctionOperation{Center, Center, Face}(DepthDependantSinkingSpeed{FT}(), + grid, + mixed_layer_depth, + euphotic_depth))), + open_bottom = true, + + scale_negatives = false, + invalid_fill_value = NaN, + + sediment = nothing, + particles = nothing, + modifiers = nothing) @warn "This implementation of PISCES is in early development and has not yet been validated against the operational version" From 687a93db81c70dbfc6b70dc03d104fcd893dd3bc Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 12 Nov 2024 11:23:47 +0000 Subject: [PATCH 6/8] fixed sugar kelp --- src/Models/Individuals/SugarKelp/SugarKelp.jl | 176 ++++++++++++------ src/Models/Individuals/SugarKelp/equations.jl | 2 +- src/Models/Individuals/SugarKelp/show.jl | 2 +- 3 files changed, 124 insertions(+), 56 deletions(-) diff --git a/src/Models/Individuals/SugarKelp/SugarKelp.jl b/src/Models/Individuals/SugarKelp/SugarKelp.jl index 7e6ad2a3e..ca459b121 100644 --- a/src/Models/Individuals/SugarKelp/SugarKelp.jl +++ b/src/Models/Individuals/SugarKelp/SugarKelp.jl @@ -21,69 +21,137 @@ export SugarKelp, SugarKelpParticles using Roots +using Oceananigans.Grids: AbstractGrid using Oceananigans.Units using OceanBioME.Particles: BiogeochemicalParticles import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers -# disgsuting number of parameters """ SugarKelp Defines the parameters for `SugarKelp` biogeochemistry. - -TODO: make a function to make it so this can be automatically converted into different float types as that wouldn't work at the moment -but there is too many parameters """ -@kwdef struct SugarKelp{FT, TL} - temperature_limit :: TL = LinearOptimalTemperatureRange() # TODO: split up more parameterisations like this - growth_rate_adjustment :: FT = 4.5 - photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60) - minimum_carbon_reserve :: FT = 0.01 - structural_carbon :: FT = 0.2 - exudation :: FT = 0.5 - erosion_exponent :: FT = 0.22 - base_erosion_rate :: FT = 10^-6 - saturation_irradiance :: FT = 90 * day/ (10 ^ 6) - structural_dry_weight_per_area :: FT = 0.5 - structural_dry_to_wet_weight :: FT = 0.0785 - carbon_reserve_per_carbon :: FT = 2.1213 - nitrogen_reserve_per_nitrogen :: FT = 2.72 - minimum_nitrogen_reserve :: FT = 0.0126 - maximum_nitrogen_reserve :: FT = 0.0216 - growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2 - maximum_specific_growth_rate :: FT = 0.18 - structural_nitrogen :: FT = 0.0146 - photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24 - photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24 - photosynthesis_ref_temp_1 :: FT = 285.0 - photosynthesis_ref_temp_2 :: FT = 288.0 - photoperiod_1 :: FT = 0.85 - photoperiod_2 :: FT = 0.3 - respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24 - respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24 - respiration_ref_temp_1 :: FT = 285.0 - respiration_ref_temp_2 :: FT = 290.0 - photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1) - photosynthesis_low_temp :: FT = 271.0 - photosynthesis_high_temp :: FT = 296.0 - photosynthesis_high_arrhenius_temp :: FT = 1414.87 - photosynthesis_low_arrhenius_temp :: FT = 4547.89 - respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1) - current_speed_for_0p65_uptake :: FT = 0.03 - nitrate_half_saturation :: FT = 4.0 - ammonia_half_saturation :: FT = 1.3 - maximum_nitrate_uptake :: FT = 10 / structural_dry_weight_per_area * 24 * 14 / (10^6) # / dm^2 / h to - maximum_ammonia_uptake :: FT = 12 / structural_dry_weight_per_area * 24 * 14 / (10^6) - current_1 :: FT = 0.72 - current_2 :: FT = 0.28 - current_3 :: FT = 0.045 - base_activity_respiration_rate :: FT = 1.11e-4 * 24 - base_basal_respiration_rate :: FT = 5.57e-5 * 24 - exudation_redfield_ratio :: FT = Inf - adapted_latitude :: FT = 57.5 +struct SugarKelp{FT, TL} + temperature_limit :: TL + growth_rate_adjustment :: FT + photosynthetic_efficiency :: FT + minimum_carbon_reserve :: FT + structural_carbon :: FT + exudation :: FT + erosion_exponent :: FT + base_erosion_rate :: FT + saturation_irradiance :: FT + structural_dry_weight_per_area :: FT + structural_dry_to_wet_weight :: FT + carbon_reserve_per_carbon :: FT + nitrogen_reserve_per_nitrogen :: FT + minimum_nitrogen_reserve :: FT + maximum_nitrogen_reserve :: FT + growth_adjustment_2 :: FT + growth_adjustment_1 :: FT + maximum_specific_growth_rate :: FT + structural_nitrogen :: FT + photosynthesis_at_ref_temp_1 :: FT + photosynthesis_at_ref_temp_2 :: FT + photosynthesis_ref_temp_1 :: FT + photosynthesis_ref_temp_2 :: FT + photoperiod_1 :: FT + photoperiod_2 :: FT + respiration_at_ref_temp_1 :: FT + respiration_at_ref_temp_2 :: FT + respiration_ref_temp_1 :: FT + respiration_ref_temp_2 :: FT + photosynthesis_arrhenius_temp :: FT + photosynthesis_low_temp :: FT + photosynthesis_high_temp :: FT + photosynthesis_high_arrhenius_temp :: FT + photosynthesis_low_arrhenius_temp :: FT + respiration_arrhenius_temp :: FT + current_speed_for_0p65_uptake :: FT + nitrate_half_saturation :: FT + ammonia_half_saturation :: FT + maximum_nitrate_uptake :: FT + maximum_ammonia_uptake :: FT + current_1 :: FT + current_2 :: FT + current_3 :: FT + base_activity_respiration_rate :: FT + base_basal_respiration_rate :: FT + exudation_redfield_ratio :: FT + adapted_latitude :: FT + + function SugarKelp(FT = Float64; + temperature_limit::TL = LinearOptimalTemperatureRange{FT}(), + growth_rate_adjustment = 4.5, + photosynthetic_efficiency = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60), + minimum_carbon_reserve = 0.01, + structural_carbon = 0.2, + exudation = 0.5, + erosion_exponent = 0.22, + base_erosion_rate = 10^-6, + saturation_irradiance = 90 * day/ (10 ^ 6), + structural_dry_weight_per_area = 0.5, + structural_dry_to_wet_weight = 0.0785, + carbon_reserve_per_carbon = 2.1213, + nitrogen_reserve_per_nitrogen = 2.72, + minimum_nitrogen_reserve = 0.0126, + maximum_nitrogen_reserve = 0.0216, + growth_adjustment_2 = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)), + growth_adjustment_1 = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2, + maximum_specific_growth_rate = 0.18, + structural_nitrogen = 0.0146, + photosynthesis_at_ref_temp_1 = 1.22e-3 * 24, + photosynthesis_at_ref_temp_2 = 1.3e-3 * 24, + photosynthesis_ref_temp_1 = 285.0, + photosynthesis_ref_temp_2 = 288.0, + photoperiod_1 = 0.85, + photoperiod_2 = 0.3, + respiration_at_ref_temp_1 = 2.785e-4 * 24, + respiration_at_ref_temp_2 = 5.429e-4 * 24, + respiration_ref_temp_1 = 285.0, + respiration_ref_temp_2 = 290.0, + photosynthesis_arrhenius_temp = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1), + photosynthesis_low_temp = 271.0, + photosynthesis_high_temp = 296.0, + photosynthesis_high_arrhenius_temp = 1414.87, + photosynthesis_low_arrhenius_temp = 4547.89, + respiration_arrhenius_temp = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1), + current_speed_for_0p65_uptake = 0.03, + nitrate_half_saturation = 4.0, + ammonia_half_saturation = 1.3, + maximum_nitrate_uptake = 10 / structural_dry_weight_per_area * 24 * 14 / (10^6), + maximum_ammonia_uptake = 12 / structural_dry_weight_per_area * 24 * 14 / (10^6), + current_1 = 0.72, + current_2 = 0.28, + current_3 = 0.045, + base_activity_respiration_rate = 1.11e-4 * 24, + base_basal_respiration_rate = 5.57e-5 * 24, + exudation_redfield_ratio = Inf, + adapted_latitude = 57.5) where TL + + return new{FT, TL}(temperature_limit, + growth_rate_adjustment, photosynthetic_efficiency, + minimum_carbon_reserve, structural_carbon, + exudation, erosion_exponent, base_erosion_rate, + saturation_irradiance, + structural_dry_weight_per_area, structural_dry_to_wet_weight, + carbon_reserve_per_carbon, nitrogen_reserve_per_nitrogen, + minimum_nitrogen_reserve, maximum_nitrogen_reserve, + growth_adjustment_2, growth_adjustment_1, maximum_specific_growth_rate, + structural_nitrogen, + photosynthesis_at_ref_temp_1, photosynthesis_at_ref_temp_2, + photosynthesis_ref_temp_1, photosynthesis_ref_temp_2, photoperiod_1, photoperiod_2, + respiration_at_ref_temp_1, respiration_at_ref_temp_2, respiration_ref_temp_1, respiration_ref_temp_2, + photosynthesis_arrhenius_temp, photosynthesis_low_temp, photosynthesis_high_temp, + photosynthesis_high_arrhenius_temp, photosynthesis_low_arrhenius_temp, + respiration_arrhenius_temp, + current_speed_for_0p65_uptake, nitrate_half_saturation, ammonia_half_saturation, + maximum_nitrate_uptake, maximum_ammonia_uptake, current_1, current_2, current_3, + base_activity_respiration_rate, base_basal_respiration_rate, exudation_redfield_ratio, + adapted_latitude) + end end """ @@ -92,8 +160,8 @@ end Sets up `n` sugar kelp `BiogeochemicalParticles` with default parameters except those specified in `kelp_parameters`. `kwagrs` are passed onto `BiogeochemicalParticles`. """ -SugarKelpParticles(n; grid, kelp_parameters = NamedTuple(), kwargs...) = - BiogeochemicalParticles(n; grid, biogeochemistry = SugarKelp(; kelp_parameters...), kwargs...) +SugarKelpParticles(n; grid::AbstractGrid{FT}, kelp_parameters = NamedTuple(), kwargs...) where FT = + BiogeochemicalParticles(n; grid, biogeochemistry = SugarKelp(FT; kelp_parameters...), kwargs...) @inline required_particle_fields(::SugarKelp) = (:A, :N, :C) @inline required_tracers(::SugarKelp) = (:u, :v, :w, :T, :NO₃, :NH₄, :PAR) diff --git a/src/Models/Individuals/SugarKelp/equations.jl b/src/Models/Individuals/SugarKelp/equations.jl index a59d1c32e..417cec8f5 100644 --- a/src/Models/Individuals/SugarKelp/equations.jl +++ b/src/Models/Individuals/SugarKelp/equations.jl @@ -243,4 +243,4 @@ normed_day_length_change(φ, n) = (day_length(φ, n) - day_length(φ, n-1)) / (d δ = asin(sin(λ * π / 180) * sin(23.44 * π / 180)) ω = (sin(-0.83 * π / 180) * sin(φ * π / 180) * sin(δ)) / (cos(φ * π / 180) * cos(δ)) return ω / 180 -end \ No newline at end of file +end diff --git a/src/Models/Individuals/SugarKelp/show.jl b/src/Models/Individuals/SugarKelp/show.jl index 0a8f68b35..264d22831 100644 --- a/src/Models/Individuals/SugarKelp/show.jl +++ b/src/Models/Individuals/SugarKelp/show.jl @@ -1,6 +1,6 @@ import Base: summary, show -summary(::SugarKelp{FT}) where FT = "SugarKelp{FT} biogeochemistry" +summary(::SugarKelp{FT}) where FT = "SugarKelp{$FT} biogeochemistry" show(io::IO, kelp::SugarKelp) = print(io, summary(kelp), From 4141d93d298c1afb8d52b00d84d57b500fa9561a Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 12 Nov 2024 12:54:26 +0000 Subject: [PATCH 7/8] fix docstring --- docs/src/model_components/individuals/slatissima.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/model_components/individuals/slatissima.md b/docs/src/model_components/individuals/slatissima.md index fe95bbff9..8a85d2ab5 100644 --- a/docs/src/model_components/individuals/slatissima.md +++ b/docs/src/model_components/individuals/slatissima.md @@ -14,7 +14,7 @@ using OceanBioME kelp_bgc = SugarKelp() # output -SugarKelp{FT} biogeochemistry (Broch & Slagstad, 2012) tracking the `N`itrogen and `C`arbon in a frond of `A`rea +SugarKelp{Float64} biogeochemistry (Broch & Slagstad, 2012) tracking the `N`itrogen and `C`arbon in a frond of `A`rea ``` which can be put into `BiogeochemicalParticles`, or you can directly manifest particles: ```jldoctest @@ -24,7 +24,7 @@ grid = RectilinearGrid(size = (1, 1, 1), extent = (1, 1, 1)); particles = SugarKelpParticles(10; grid) # output -10 BiogeochemicalParticles with SugarKelp{FT} biogeochemistry: +10 BiogeochemicalParticles with SugarKelp{Float64} biogeochemistry: ├── fields: (:A, :N, :C) └── coupled tracers: (:NO₃, :NH₄, :DIC, :O₂, :DOC, :DON, :bPOC, :bPON) From dbce30e132b1adec3c1daec87836d1d4514e33da Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 13 Nov 2024 19:01:48 +0000 Subject: [PATCH 8/8] bump minor version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 56cb17ddf..ab8dbbfda 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OceanBioME" uuid = "a49af516-9db8-4be4-be45-1dad61c5a376" authors = ["Jago Strong-Wright and contributors"] -version = "0.13.2" +version = "0.13.3" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"