diff --git a/Project.toml b/Project.toml index ebbdf412a..a8c3ef66f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "OceanBioME" uuid = "a49af516-9db8-4be4-be45-1dad61c5a376" authors = ["Jago Strong-Wright ", "John R Taylor ", "Si Chen "] -version = "0.4.1" +version = "0.5.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/src/model_implementation.md b/docs/src/model_implementation.md index 648f35a7e..1d3d864ac 100644 --- a/docs/src/model_implementation.md +++ b/docs/src/model_implementation.md @@ -225,10 +225,10 @@ using OceanBioME.Boundaries.Sediments: sinking_flux import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver -@inline nitrogen_flux(grid, advection, bgc::NutrientPhytoplankton, tracers, i, j) = - sinking_flux(i, j, grid, advection, Val(:P), bgc, tracers) +@inline nitrogen_flux(i, j, k, grid, advection, bgc::NutrientPhytoplankton, tracers) = + sinking_flux(i, j, k, grid, advection, Val(:P), bgc, tracers) -@inline carbon_flux(bgc::NutrientPhytoplankton, tracers, i, j) = nitrogen_flux(bgc, tracers, i, j) * 6.56 +@inline carbon_flux(bgc::NutrientPhytoplankton, tracers, i, j) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 @inline remineralisation_receiver(::NutrientPhytoplankton) = :N ``` diff --git a/src/Boundaries/Sediments/Sediments.jl b/src/Boundaries/Sediments/Sediments.jl index 3487b8cad..e295ab63c 100644 --- a/src/Boundaries/Sediments/Sediments.jl +++ b/src/Boundaries/Sediments/Sediments.jl @@ -3,7 +3,7 @@ module Sediments export SimpleMultiG, InstantRemineralisation using KernelAbstractions -using OceanBioME: ContinuousFormBiogeochemistry +using OceanBioME: ContinuousFormBiogeochemistry, BoxModelGrid using Oceananigans using Oceananigans.Architectures: device using Oceananigans.Utils: launch! @@ -13,13 +13,16 @@ using Oceananigans.Fields: CenterField, Face using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity using Oceananigans.Grids: zspacing using Oceananigans.Operators: volume -using Oceananigans.Fields: Center +using Oceananigans.Fields: Center, ConstantField +using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, immersed_cell import Adapt: adapt_structure, adapt abstract type AbstractSediment end abstract type FlatSediment <: AbstractSediment end +@inline bottom_index(i, j, sediment) = @inbounds bottom_index_array(sediment)[i, j] + sediment_fields(::AbstractSediment) = () @inline update_sediment_tendencies!(bgc, sediment, model) = nothing @@ -43,14 +46,42 @@ end @inbounds G⁻[i, j, 1] = G⁰[i, j, 1] end +@inline nitrogen_flux() = 0 +@inline carbon_flux() = 0 +@inline remineralisation_receiver() = nothing + +@inline advection_scheme(advection, val_tracer) = advection +@inline advection_scheme(advection::NamedTuple, val_tracer::Val{T}) where T = advection[T] + +@inline function sinking_flux(i, j, k, grid, advection, val_tracer::Val{T}, bgc, tracers) where T + return - advective_tracer_flux_z(i, j, k, grid, advection_scheme(advection, val_tracer), biogeochemical_drift_velocity(bgc, val_tracer).w, tracers[T]) / + volume(i, j, k, grid, Center(), Center(), Center()) +end + +calculate_bottom_indices(::BoxModelGrid) = 1 +calculate_bottom_indices(grid) = ones(Int, size(grid)[1:2]...) + +@kernel function find_bottom_cell(grid, bottom_indices) + i, j = @index(Global, NTuple) + + Nz = size(grid, 3) -@inline nitrogen_flux(grid, adveciton, bgc, tracers, i, j) = 0 -@inline carbon_flux(grid, adveciton, bgc, tracers, i, j) = 0 -@inline remineralisation_receiver(bgc) = nothing + k_bottom = 1 -@inline sinking_flux(i, j, grid, advection, val_tracer::Val{T}, bgc, tracers) where T = - - advective_tracer_flux_z(i, j, 1, grid, advection, biogeochemical_drift_velocity(bgc, val_tracer).w, tracers[T]) / - volume(i, j, 1, grid, Center(), Center(), Center()) + while immersed_cell(i, j, k_bottom, grid.underlying_grid, grid.immersed_boundary) && k_bottom < Nz + k_bottom += 1 + end + + @inbounds bottom_indices[i, j] = k_bottom +end + +function calculate_bottom_indices(grid::ImmersedBoundaryGrid) + indices = zeros(Int, size(grid)[1:2]...) + + launch!(grid.architecture, grid, :xy, find_bottom_cell, grid, indices) + + return indices +end include("coupled_timesteppers.jl") include("simple_multi_G.jl") diff --git a/src/Boundaries/Sediments/instant_remineralization.jl b/src/Boundaries/Sediments/instant_remineralization.jl index 2c08caad3..1ea6e8d09 100644 --- a/src/Boundaries/Sediments/instant_remineralization.jl +++ b/src/Boundaries/Sediments/instant_remineralization.jl @@ -7,24 +7,27 @@ becoming N, and a fraction being perminantly burried. Burial efficiency from [RemineralisationFraction](@citet). """ -struct InstantRemineralisation{FT, F, TE} <: FlatSediment +struct InstantRemineralisation{FT, F, TE, B} <: FlatSediment burial_efficiency_constant1 :: FT burial_efficiency_constant2 :: FT burial_efficiency_half_saturaiton :: FT fields :: F tendencies :: TE + bottom_indices :: B InstantRemineralisation(burial_efficiency_constant1::FT, burial_efficiency_constant2::FT, burial_efficiency_half_saturaiton::FT, fields::F, - tendencies::TE) where {FT, F, TE} = - new{FT, F, TE}(burial_efficiency_constant1, - burial_efficiency_constant2, - burial_efficiency_half_saturaiton, - fields, - tendencies) + tendencies::TE, + bottom_indices::B) where {FT, F, TE, B} = + new{FT, F, TE, B}(burial_efficiency_constant1, + burial_efficiency_constant2, + burial_efficiency_half_saturaiton, + fields, + tendencies, + bottom_indices) end """ @@ -56,15 +59,18 @@ function InstantRemineralisation(; grid, tracer_names = (:N_storage, ) # add field slicing back ( indices = (:, :, 1)) when output writer can cope - fields = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)) - tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)), - G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names))) + fields = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)) + tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)), + G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names))) + + bottom_indices = calculate_bottom_indices(grid) return InstantRemineralisation(burial_efficiency_constant1, burial_efficiency_constant2, burial_efficiency_half_saturaiton, fields, - tendencies) + tendencies, + bottom_indices) end adapt_structure(to, sediment::InstantRemineralisation) = @@ -77,12 +83,16 @@ adapt_structure(to, sediment::InstantRemineralisation) = sediment_tracers(::InstantRemineralisation) = (:N_storage, ) sediment_fields(model::InstantRemineralisation) = (N_storage = model.fields.N_storage, ) +@inline bottom_index_array(sediment::InstantRemineralisation) = sediment.bottom_indices + @kernel function _calculate_tendencies!(sediment::InstantRemineralisation, bgc, grid, advection, tracers, timestepper) i, j = @index(Global, NTuple) - Δz = zspacing(i, j, 1, grid, Center(), Center(), Center()) + k = bottom_index(i, j, sediment) + + Δz = zspacing(i, j, k, grid, Center(), Center(), Center()) - flux = nitrogen_flux(grid, advection, bgc, tracers, i, j) * Δz + flux = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * Δz burial_efficiency = sediment.burial_efficiency_constant1 + sediment.burial_efficiency_constant2 * ((flux * 6.56) / (7 + flux * 6.56)) ^ 2 diff --git a/src/Boundaries/Sediments/simple_multi_G.jl b/src/Boundaries/Sediments/simple_multi_G.jl index ab62f856d..77e804910 100644 --- a/src/Boundaries/Sediments/simple_multi_G.jl +++ b/src/Boundaries/Sediments/simple_multi_G.jl @@ -6,7 +6,7 @@ import Base: show, summary Hold the parameters and fields for a simple "multi G" single-layer sediment model. Based on the Level 3 model described by [Soetaert2000](@citet). """ -struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE} <: FlatSediment +struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE, B} <: FlatSediment fast_decay_rate :: FT slow_decay_rate :: FT @@ -24,6 +24,7 @@ struct SimpleMultiG{FT, P1, P2, P3, P4, F, TE} <: FlatSediment fields :: F tendencies :: TE + bottom_indices :: B end """ @@ -108,22 +109,26 @@ function SimpleMultiG(; grid, tracer_names = (:C_slow, :C_fast, :N_slow, :N_fast, :C_ref, :N_ref) # add field slicing back ( indices = (:, :, 1)) when output writer can cope - fields = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)) - tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names)), - G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid;) for tracer in tracer_names))) + fields = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)) + tendencies = (Gⁿ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names)), + G⁻ = NamedTuple{tracer_names}(Tuple(CenterField(grid) for tracer in tracer_names))) + + bottom_indices = calculate_bottom_indices(grid) F = typeof(fields) TE = typeof(tendencies) - - return SimpleMultiG{FT, P1, P2, P3, P4, F, TE}(fast_decay_rate, slow_decay_rate, - fast_redfield, slow_redfield, - fast_fraction, slow_fraction, refactory_fraction, - nitrate_oxidation_params, - denitrification_params, - anoxic_params, - solid_dep_params, - fields, - tendencies) + B = typeof(bottom_indices) + + return SimpleMultiG{FT, P1, P2, P3, P4, F, TE, B}(fast_decay_rate, slow_decay_rate, + fast_redfield, slow_redfield, + fast_fraction, slow_fraction, refactory_fraction, + nitrate_oxidation_params, + denitrification_params, + anoxic_params, + solid_dep_params, + fields, + tendencies, + bottom_indices) end adapt_structure(to, sediment::SimpleMultiG) = @@ -149,16 +154,20 @@ sediment_fields(model::SimpleMultiG) = (C_slow = model.fields.C_slow, C_ref = model.fields.C_ref, N_ref = model.fields.N_ref) +@inline bottom_index_array(sediment::SimpleMultiG) = sediment.bottom_indices + @kernel function _calculate_tendencies!(sediment::SimpleMultiG, bgc, grid, advection, tracers, timestepper) i, j = @index(Global, NTuple) - Δz = zspacing(i, j, 1, grid, Center(), Center(), Center()) + k = bottom_index(i, j, sediment) + + Δz = zspacing(i, j, k, grid, Center(), Center(), Center()) @inbounds begin - carbon_deposition = carbon_flux(grid, advection, bgc, tracers, i, j) * Δz + carbon_deposition = carbon_flux(i, j, k, grid, advection, bgc, tracers) * Δz - nitrogen_deposition = nitrogen_flux(grid, advection, bgc, tracers, i, j) * Δz + nitrogen_deposition = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * Δz # rates C_min_slow = sediment.fields.C_slow[i, j, 1] * sediment.slow_decay_rate diff --git a/src/Light/2band.jl b/src/Light/2band.jl index f4b1f4a75..6f7b4bd32 100644 --- a/src/Light/2band.jl +++ b/src/Light/2band.jl @@ -105,8 +105,8 @@ function TwoBandPhotosyntheticallyActiveRadiation(; grid, field = CenterField(grid; boundary_conditions = regularize_field_boundary_conditions( - FieldBoundaryConditions(top = ValueBoundaryCondition(surface_PAR)), - grid, :PAR)) + FieldBoundaryConditions(top = ValueBoundaryCondition(surface_PAR)), + grid, :PAR)) return TwoBandPhotosyntheticallyActiveRadiation(water_red_attenuation, water_blue_attenuation, diff --git a/src/Light/Light.jl b/src/Light/Light.jl index 2a561d3b3..f34cfffad 100644 --- a/src/Light/Light.jl +++ b/src/Light/Light.jl @@ -18,6 +18,7 @@ using Oceananigans.BoundaryConditions: fill_halo_regions!, regularize_field_boundary_conditions, ContinuousBoundaryFunction using Oceananigans.Units +using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid import Adapt: adapt_structure, adapt import Base: show, summary diff --git a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl index 36cd9202e..d728ade6a 100644 --- a/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl +++ b/src/Models/AdvectedPopulations/LOBSTER/LOBSTER.jl @@ -321,6 +321,10 @@ function LOBSTER(; grid, particles::P = nothing) where {FT, LA, S, P} + if !isnothing(sediment_model) && !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 + sinking_velocities = setup_velocity_fields(sinking_speeds, grid, open_bottom) optionals = Val((carbonates, oxygen, variable_redfield)) @@ -469,19 +473,19 @@ const lobster_variable_redfield = Union{LOBSTER{<:Any, <:Any, <:Any, <:Val{(fals LOBSTER{<:Any, <:Any, <:Any, <:Val{(true, true, true)}, <:Any, <:Any}} -@inline nitrogen_flux(grid, advection, bgc::LOBSTER, tracers, i, j) = - sinking_flux(i, j, grid, advection, Val(:sPOM), bgc, tracers) + - sinking_flux(i, j, grid, advection, Val(:bPOM), bgc, tracers) +@inline nitrogen_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = + sinking_flux(i, j, k, grid, advection, Val(:sPOM), bgc, tracers) + + sinking_flux(i, j, k, grid, advection, Val(:bPOM), bgc, tracers) -@inline carbon_flux(grid, advection, bgc::LOBSTER, tracers, i, j) = nitrogen_flux(grid, advection, bgc, tracers, i, j) * bgc.organic_redfield +@inline carbon_flux(i, j, k, grid, advection, bgc::LOBSTER, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * bgc.organic_redfield -@inline nitrogen_flux(grid, advection, bgc::lobster_variable_redfield, tracers, i, j) = - sinking_flux(i, j, grid, advection, Val(:sPON), bgc, tracers) + - sinking_flux(i, j, grid, advection, Val(:bPON), bgc, tracers) +@inline nitrogen_flux(i, j, k, grid, advection, bgc::lobster_variable_redfield, tracers) = + sinking_flux(i, j, k, grid, advection, Val(:sPON), bgc, tracers) + + sinking_flux(i, j, k, grid, advection, Val(:bPON), bgc, tracers) -@inline carbon_flux(grid, advection, bgc::lobster_variable_redfield, tracers, i, j) = - sinking_flux(i, j, grid, advection, Val(:sPOC), bgc, tracers) + - sinking_flux(i, j, grid, advection, Val(:bPOC), bgc, tracers) +@inline carbon_flux(i, j, k, grid, advection, bgc::lobster_variable_redfield, tracers) = + sinking_flux(i, j, k, grid, advection, Val(:sPOC), bgc, tracers) + + sinking_flux(i, j, k, grid, advection, Val(:bPOC), bgc, tracers) @inline remineralisation_receiver(::LOBSTER) = :NH₄ end # module diff --git a/src/Models/AdvectedPopulations/NPZD.jl b/src/Models/AdvectedPopulations/NPZD.jl index 3d4981b0a..704bbb0e6 100644 --- a/src/Models/AdvectedPopulations/NPZD.jl +++ b/src/Models/AdvectedPopulations/NPZD.jl @@ -199,6 +199,10 @@ function NutrientPhytoplanktonZooplanktonDetritus(; grid, sinking_velocities = setup_velocity_fields(sinking_speeds, grid, open_bottom) + if !isnothing(sediment_model) && !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 + return NutrientPhytoplanktonZooplanktonDetritus(initial_photosynthetic_slope, base_maximum_growth, nutrient_half_saturation, @@ -339,9 +343,10 @@ adapt_structure(to, npzd::NPZD) = adapt(to, npzd.particles)) -@inline nitrogen_flux(grid, advection, bgc::NPZD, tracers, i, j) = sinking_flux(i, j, grid, advection, Val(:D), bgc, tracers) + - sinking_flux(i, j, grid, advection, Val(:P), bgc, tracers) - -@inline carbon_flux(bgc::NPZD, tracers, i, j) = nitrogen_flux(bgc, tracers, i, j) * 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) + +@inline carbon_flux(i, j, k, grid, advection, bgc::NPZD, tracers) = nitrogen_flux(i, j, k, grid, advection, bgc, tracers) * 6.56 + @inline remineralisation_receiver(::NPZD) = :N end # module diff --git a/test/test_sediments.jl b/test/test_sediments.jl index 2d6ebede1..6e48e650c 100644 --- a/test/test_sediments.jl +++ b/test/test_sediments.jl @@ -9,6 +9,8 @@ using Oceananigans.Fields: TracerFields using Oceananigans.Operators: volume, Azᶠᶜᶜ +using OceanBioME.LOBSTERModel: lobster_variable_redfield + function intercept_tendencies!(model, intercepted_tendencies) for tracer in keys(model.tracers) copyto!(intercepted_tendencies[tracer], model.timestepper.Gⁿ[tracer]) @@ -25,7 +27,14 @@ end set_defaults!(::InstantRemineralisation) = nothing -function set_defaults!(::LOBSTER, model) +set_defaults!(::LOBSTER, model) = + set!(model, P = 0.4686, Z = 0.5363, + NO₃ = 2.3103, NH₄ = 0.0010, + DOM = 0.8115, + sPOM = 0.2299, bPOM = 0.0103) + + +set_defaults!(::lobster_variable_redfield, model) = set!(model, P = 0.4686, Z = 0.5363, NO₃ = 2.3103, NH₄ = 0.0010, DIC = 2106.9, Alk = 2408.9, @@ -33,7 +42,7 @@ function set_defaults!(::LOBSTER, model) DOC = 5.3390, DON = 0.8115, sPON = 0.2299, sPOC = 1.5080, bPON = 0.0103, bPOC = 0.0781) -end + set_defaults!(::NutrientPhytoplanktonZooplanktonDetritus, model) = set!(model, N = 2.3, P = 0.4, Z = 0.5, D = 0.2) @@ -43,13 +52,21 @@ total_nitrogen(sed::SimpleMultiG) = sum(sed.fields.N_fast) + total_nitrogen(sed::InstantRemineralisation) = sum(sed.fields.N_storage) -total_nitrogen(::LOBSTER, model) = sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.DON) + sum(model.tracers.sPON) + sum(model.tracers.bPON) +total_nitrogen(::LOBSTER, model) = sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.DOM) + sum(model.tracers.sPOM) + sum(model.tracers.bPOM) +total_nitrogen(::lobster_variable_redfield, model) = sum(model.tracers.NO₃) + sum(model.tracers.NH₄) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.DON) + sum(model.tracers.sPON) + sum(model.tracers.bPON) total_nitrogen(::NutrientPhytoplanktonZooplanktonDetritus, model) = sum(model.tracers.N) + sum(model.tracers.P) + sum(model.tracers.Z) + sum(model.tracers.D) -function test_flat_sediment(grid, biogeochemistry; timestepper = :QuasiAdamsBashforth2) - model = NonhydrostaticModel(; grid, biogeochemistry, - closure = nothing, - timestepper) +function test_flat_sediment(grid, biogeochemistry, model; timestepper = :QuasiAdamsBashforth2) + model = isa(model, NonhydrostaticModel) ? model(; grid, + biogeochemistry, + closure = nothing, + timestepper, + buoyancy = nothing) : + model(; grid, + biogeochemistry, + closure = nothing, + buoyancy = nothing, + tracers = nothing) set_defaults!(model.biogeochemistry.sediment_model) @@ -88,23 +105,38 @@ display_name(::LOBSTER) = "LOBSTER" display_name(::NutrientPhytoplanktonZooplanktonDetritus) = "NPZD" display_name(::SimpleMultiG) = "Multi-G" display_name(::InstantRemineralisation) = "Instant remineralisation" +display_name(::RectilinearGrid) = "Rectilinear grid" +display_name(::LatitudeLongitudeGrid) = "Latitude longitude grid" +display_name(::ImmersedBoundaryGrid) = "Immersed boundary grid" + + +bottom_height(x, y) = -1000 + 500 * exp(- (x^2 + y^2) / 250) # a perfect hill -@testset "Sediment" begin +@testset "Sediment integration" begin for architecture in (CPU(), ) - grid = RectilinearGrid(architecture; size=(3, 3, 50), extent=(10, 10, 500)) - - for timestepper in (:QuasiAdamsBashforth2, :RungeKutta3), - sediment_model in (InstantRemineralisation(; grid), SimpleMultiG(; grid)) - for biogeochemistry in (NutrientPhytoplanktonZooplanktonDetritus(; grid, open_bottom = true, - sediment_model), - LOBSTER(; grid, - carbonates = true, oxygen = true, variable_redfield = true, - open_bottom = true, - sediment_model)) - - if !(isa(sediment_model, SimpleMultiG) && isa(biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)) - @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) on $(display_name(biogeochemistry))" - @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry))" test_flat_sediment(grid, biogeochemistry; timestepper) + grids = [RectilinearGrid(architecture; size=(3, 3, 50), extent=(10, 10, 500)), + LatitudeLongitudeGrid(architecture; size = (3, 3, 16), latitude = (0, 10), longitude = (0, 10), z = (-500, 0)), + ImmersedBoundaryGrid( + LatitudeLongitudeGrid(architecture; size = (3, 3, 16), latitude = (0, 10), longitude = (0, 10), z = (-500, 0)), + GridFittedBottom(bottom_height))] + for grid in grids + for timestepper in (:QuasiAdamsBashforth2, :RungeKutta3), + sediment_model in (InstantRemineralisation(; grid), SimpleMultiG(; grid)), + model in (NonhydrostaticModel, HydrostaticFreeSurfaceModel) + for biogeochemistry in (NutrientPhytoplanktonZooplanktonDetritus(; grid, sediment_model), + LOBSTER(; grid, + carbonates = ifelse(isa(sediment_model, SimpleMultiG), true, false), + oxygen = ifelse(isa(sediment_model, SimpleMultiG), true, false), + variable_redfield = ifelse(isa(sediment_model, SimpleMultiG), true, false), + sediment_model)) + # get rid of incompatible combinations + run = ifelse((model == NonhydrostaticModel && (isa(grid, ImmersedBoundaryGrid) || isa(grid, LatitudeLongitudeGrid))) || + (model == HydrostaticFreeSurfaceModel && timestepper == :RungeKutta3) || + (isa(sediment_model, SimpleMultiG) && isa(biogeochemistry, NutrientPhytoplanktonZooplanktonDetritus)), false, true) + if run + @info "Testing sediment on $(typeof(architecture)) with $timestepper and $(display_name(sediment_model)) in $(display_name(biogeochemistry)) on $(display_name(grid)) with $(model)" + @testset "$architecture, $timestepper, $(display_name(sediment_model)), $(display_name(biogeochemistry)), $(display_name(grid)), $(model)" test_flat_sediment(grid, biogeochemistry, model; timestepper) + end end end end