Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Hmw ct jsw/pisces #5

Merged
merged 4 commits into from
Aug 24, 2024
Merged

Hmw ct jsw/pisces #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9a34b35
fixed boundaries thing
jagoosw Aug 23, 2024
978eef9
super to sub script for par names
jagoosw Aug 23, 2024
7051614
made column script run again
jagoosw Aug 23, 2024
5860655
Merge branch 'hmw-ct-jsw/PISCES' into hmw-ct-jsw/PISCES
ciadht Aug 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 7 additions & 4 deletions src/Models/AdvectedPopulations/PISCES/PISCES.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,11 +47,13 @@ using Oceananigans.Fields: Field, TracerFields, CenterField, ZeroField, Constant
using OceanBioME.Light: TwoBandPhotosyntheticallyActiveRadiation, default_surface_PAR
using OceanBioME: setup_velocity_fields, show_sinking_velocities, Biogeochemistry, ScaleNegativeTracers
using OceanBioME.BoxModels: BoxModel
using OceanBioME.Boundaries.Sediments: sinking_flux

using Oceananigans.Biogeochemistry: AbstractContinuousFormBiogeochemistry

import OceanBioME: redfield, conserved_tracers, chlorophyll


import OceanBioME: redfield, conserved_tracers, maximum_sinking_velocity, chlorophyll


import Oceananigans.Biogeochemistry: required_biogeochemical_tracers,
required_biogeochemical_auxiliary_fields,
Expand All @@ -63,10 +65,11 @@ import OceanBioME: maximum_sinking_velocity
import Adapt: adapt_structure, adapt
import Base: show, summary

import OceanBioME.Boundaries.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers
import OceanBioME.Models.Sediments: nitrogen_flux, carbon_flux, remineralisation_receiver, sinking_tracers

struct PISCES{FT, PD, ZM, OT, W, CF, ZF, FFMLD, FFEU} <: AbstractContinuousFormBiogeochemistry


growth_rate_at_zero :: FT # add list of parameters here, assuming theyre all just numbers FT will be fine for advect_particles_kernel
growth_rate_reference_for_light_limitation :: FT
basal_respiration_rate :: FT
Expand Down Expand Up @@ -915,7 +918,7 @@ adapt_structure(to, pisces::PISCES) =
# you can updatye these if you want it to have a pretty way of showing uyou its a pisces model
summary(::PISCES{FT}) where {FT} = string("PISCES{$FT}")

show(io::IO, model::PISCES) where {FT, B, W, PD, ZM, OT} = print(io, string("Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) model")) # maybe add some more info here
show(io::IO, model::PISCES) = print(io, string("Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) model")) # maybe add some more info here

@inline maximum_sinking_velocity(bgc::PISCES) = maximum(abs, bgc.sinking_velocities.bPOM.w) # might need ot update this for wghatever the fastest sinking pareticles are

Expand Down
2 changes: 1 addition & 1 deletion src/Models/AdvectedPopulations/PISCES/calcite.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ end
return rain_ratio(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, T, PAR, zₘₓₗ, bgc)*(ηᶻ*grazing_Z(P, D, POC, T, bgc)[2]*Z+ηᴹ*grazing_M(P, D, Z, POC, T, bgc)[2]*M + 0.5*(mᴾ*concentration_limitation(P, Kₘ)*P + sh*wᴾ*P^2)) #eq76
end

#Forcing for calcite PAR₂``
#Forcing for calcite
@inline function (bgc::PISCES)(::Val{:CaCO₃}, x, y, z, t, P, D, Z, M, Pᶜʰˡ, Dᶜʰˡ, Pᶠᵉ, Dᶠᵉ, Dˢⁱ, DOC, POC, GOC, SFe, BFe, PSi, NO₃, NH₄, PO₄, Fe, Si, CaCO₃, DIC, Alk, O₂, T, zₘₓₗ, zₑᵤ, Si̅, D_dust, Ω, PAR, PAR₁, PAR₂, PAR₃)

return (production_of_sinking_calcite(P, PO₄, NO₃, NH₄, Pᶜʰˡ, Pᶠᵉ, Fe, D, Z, M, POC, T, PAR, zₘₓₗ, z, bgc)
Expand Down
11 changes: 6 additions & 5 deletions validation/PISCES/boxPISCES.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,16 +33,17 @@ PAR_func1(t) = PAR_func(t)/3
PAR_func2(t) = PAR_func(t)/3
PAR_func3(t)= PAR_func(t)/3

PAR = FunctionField{Center, Center, Center}(PAR_func, grid; clock)
PAR¹ = FunctionField{Center, Center, Center}(PAR_func1, grid; clock)
PAR² = FunctionField{Center, Center, Center}(PAR_func2, grid; clock)
PAR³ = FunctionField{Center, Center, Center}(PAR_func3, grid; clock)

PAR₁ = FunctionField{Center, Center, Center}(PAR_func1, grid; clock)
PAR₂ = FunctionField{Center, Center, Center}(PAR_func2, grid; clock)
PAR₃ = FunctionField{Center, Center, Center}(PAR_func3, grid; clock)
zₘₓₗ = FunctionField{Center, Center, Center}(MLD, grid; clock)
zₑᵤ = FunctionField{Center, Center, Center}(EU, grid; clock)
nothing #hide

# Set up the model. Here, first specify the biogeochemical model, followed by initial conditions and the start and end times
model = BoxModel(; biogeochemistry = PISCES(; grid, light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation((; PAR, PAR¹, PAR², PAR³)), mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ),
model = BoxModel(; biogeochemistry = PISCES(; grid, light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation((; PAR, PAR₁, PAR₂, PAR₃)), mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ),

clock)

set!(model, P = 6.95, D = 6.95, Z = 0.695, M = 0.695, Pᶜʰˡ = 1.671, Dᶜʰˡ = 1.671, Pᶠᵉ =7e-6 * 1e9 / 1e6 * 6.95, Dᶠᵉ = 7e-6 * 1e9 / 1e6 * 6.95, Dˢⁱ = 1.162767, DOC = 0.0, POC = 0.0, GOC = 0.0, SFe = 7e-6 * 1e9 / 1e6 *1.256, BFe =7e-6 * 1e9 / 1e6 *1.256, NO₃ = 6.202, NH₄ = 0.25*6.202, PO₄ = 0.8722, Fe = 1.256, Si = 7.313, CaCO₃ = 0.29679635764590534, DIC = 2139.0, Alk = 2366.0, O₂ = 237.0, T = 14.0) #Using Copernicus Data (26.665, 14.)
Expand Down
25 changes: 13 additions & 12 deletions validation/PISCES/columnPISCES.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
# ## Model setup
# We load the packages and choose the default LOBSTER parameter set
using OceanBioME, Oceananigans, Printf
using OceanBioME.SLatissimaModel: SLatissima
using Oceananigans.Fields: FunctionField, ConstantField
using Oceananigans.Units

Expand All @@ -24,8 +23,7 @@ nothing #hide
# ## Surface PAR and turbulent vertical diffusivity based on idealised mixed layer depth
# Setting up idealised functions for PAR and diffusivity (details here can be ignored but these are typical of the North Atlantic)

@inline PAR⁰(x, y, t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2
#@inline PAR⁰(x,y,t) = 60.0
@inline PAR⁰(t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2

@inline H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0)

Expand Down Expand Up @@ -59,11 +57,13 @@ grid = RectilinearGrid(size = (1, 1, 100), extent = (20meters, 20meters, 400mete
clock = Clock(; time = 0.0)

PAR = FunctionField{Center, Center, Center}(PAR_func, grid; clock)
PAR¹ = FunctionField{Center, Center, Center}(PAR_func1, grid; clock)
PAR² = FunctionField{Center, Center, Center}(PAR_func2, grid; clock)
PAR³ = FunctionField{Center, Center, Center}(PAR_func3, grid; clock)

PAR₁ = FunctionField{Center, Center, Center}(PAR_func1, grid; clock)
PAR₂ = FunctionField{Center, Center, Center}(PAR_func2, grid; clock)
PAR₃ = FunctionField{Center, Center, Center}(PAR_func3, grid; clock)
zₘₓₗ = FunctionField{Center, Center, Center}(MLD, grid; clock)
zₑᵤ = FunctionField{Center, Center, Center}(euphotic, grid; clock)

#ff = FunctionField{Nothing, Nothing, Face}(w_GOC, grid)
# ## Grid
# Define the grid.
Expand All @@ -75,10 +75,11 @@ zₑᵤ = FunctionField{Center, Center, Center}(euphotic, grid; clock)
# and then setup the Oceananigans model with the boundary condition for the DIC based on the air-sea CO₂ flux.

biogeochemistry = PISCES(; grid,
light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation((; PAR, PAR¹, PAR², PAR³)), sinking_speeds = (; POC = w_POC, SFe = w_POC, GOC = w_GOC, BFe = w_GOC, PSi = w_GOC, CaCO₃ = w_GOC), mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ)

light_attenuation_model = MultiBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR = PAR⁰),
sinking_speeds = (; POC = w_POC, SFe = w_POC, GOC = w_GOC, BFe = w_GOC, PSi = w_GOC, CaCO₃ = w_GOC),
mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ)

CO₂_flux = GasExchange(; gas = :CO₂)
CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition()

funT = FunctionField{Center, Center, Center}(temp, grid; clock)
S = ConstantField(35)
Expand All @@ -88,9 +89,9 @@ model = NonhydrostaticModel(; grid,
clock,
closure = ScalarDiffusivity(VerticallyImplicitTimeDiscretization(), κ = κₜ),
biogeochemistry,
boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), ),
auxiliary_fields = (; S, zₘₓₗ, zₑᵤ, Si̅ = yearly_maximum_silicate, D_dust = dust_deposition, Ω = carbonate_sat_ratio, PAR, PAR¹, PAR², PAR³)
)
boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), ))


@info "Setting initial values..."
set!(model, P = 6.95, D = 6.95, Z = 0.695, M = 0.695, Pᶜʰˡ = 1.671, Dᶜʰˡ = 1.671, Pᶠᵉ =7e-6 * 1e9 / 1e6 * 6.95, Dᶠᵉ = 7e-6 * 1e9 / 1e6 * 6.95, Dˢⁱ = 1.162767, DOC = 0.0, POC = 0.0, GOC = 0.0, SFe = 7e-6 * 1e9 / 1e6 *1.256, BFe =7e-6 * 1e9 / 1e6 *1.256, NO₃ = 6.202e-2, NH₄ = 0.25*6.202e-2, PO₄ = 0.8722, Fe = 1.256, Si = 7.313, CaCO₃ = 0.29679635764590534, DIC = 2139.0, Alk = 2366.0, O₂ = 237.0, T = funT) #Using Copernicus Data (26.665, 14.)
# ## Simulation
Expand Down