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 #199

Merged
merged 3 commits into from
Aug 6, 2024
Merged

Hmw ct jsw/pisces #199

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a8e26de
added Carbon Checke
ciadht Aug 6, 2024
e7e3a04
Added sh to equations
hannahmw1 Aug 6, 2024
db7c3ec
added parameter scaling
ciadht Aug 6, 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
6 changes: 3 additions & 3 deletions src/Models/AdvectedPopulations/PISCES/PISCES.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -625,8 +625,8 @@ function PISCES(; grid, # finally the function
min_half_saturation_const_for_iron_uptake :: PD = (P = 1.0, D = 3.0), #nmolFeL⁻¹
size_ratio_of_phytoplankton :: PD = (P = 3.0, D = 3.0),
optimal_SiC_uptake_ratio_of_diatoms :: FT = 0.159, #molSi/(mol C)
optimal_iron_quota :: PD = (P = 7.0, D = 7.0), #mmolFe/(mol C)
max_iron_quota :: PD = (P = 40.0, D = 40.0), #molFe/(mol C)
optimal_iron_quota :: PD = (P = 7.0e-3, D = 7.0e-3), #mmolFe/(mol C)
max_iron_quota :: PD = (P = 40.0e-3, D = 40.0e-3), #molFe/(mol C)
phytoplankton_mortality_rate :: PD = (P = 0.01/day, D = 0.01/day),
min_quadratic_mortality_of_phytoplankton :: FT = 0.01 / day, #1/(d mol C)
max_quadratic_mortality_of_diatoms :: FT = 0.03 / day, #1/(d mol C)
Expand Down Expand Up @@ -715,7 +715,7 @@ function PISCES(; grid, # finally the function
dissolution_rate_of_silicon :: FT = 1.0,
coefficient_of_bacterial_uptake_of_iron_in_POC :: FT = 0.5,
coefficient_of_bacterial_uptake_of_iron_in_GOC :: FT = 0.5,
max_FeC_ratio_of_bacteria :: FT = 10.0e-6, #or 6
max_FeC_ratio_of_bacteria :: FT = 10.0e-3, #or 6
Fe_half_saturation_const_for_PLACEHOLDER :: FT = 2.5e-10, #or 2.5e-10 #not sure what this should be called
proportion_of_sinking_grazed_shells :: ZM = (Z = 0.3, M = 0.3), # 0.3 for both? not sure

Expand Down
12 changes: 2 additions & 10 deletions src/Models/AdvectedPopulations/PISCES/POC_and_GOC.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ end
gₚₒ_FFᴹ = g_FF*(bₘ^T)*wₚₒ*POC #29a
Φ = get_Φ(POC, GOC, sh, bgc)

return σᶻ*∑gᶻ*Z + 0.5*mᴰ*K_mondo(D, Kₘ)*D + rᶻ*(b_Z^T)*(K_mondo(Z, Kₘ) + 3*ΔO₂(O₂, bgc))*Z + mᶻ*(b_Z^T)*Z^2 + (1 - 0.5*R_CaCO₃)*(mᴾ*K_mondo(P, Kₘ)*P + wᴾ*P^2) + λₚₒ¹*GOC + Φ₁ᴰᴼᶜ + Φ₃ᴰᴼᶜ - (gₚₒᴹ + gₚₒ_FFᴹ)*M - gₚₒᶻ*Z - λₚₒ¹*POC - Φ #37
return σᶻ*∑gᶻ*Z + 0.5*mᴰ*K_mondo(D, Kₘ)*D + rᶻ*(b_Z^T)*(K_mondo(Z, Kₘ) + 3*ΔO₂(O₂, bgc))*Z + mᶻ*(b_Z^T)*Z^2 + (1 - 0.5*R_CaCO₃)*(mᴾ*K_mondo(P, Kₘ)*P + sh*wᴾ*P^2) + λₚₒ¹*GOC + Φ₁ᴰᴼᶜ + Φ₃ᴰᴼᶜ - (gₚₒᴹ + gₚₒ_FFᴹ)*M - gₚₒᶻ*Z - λₚₒ¹*POC - Φ #37
end

@inline function (bgc::PISCES)(::Val{:GOC}, 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³)
Expand Down Expand Up @@ -68,14 +68,6 @@ end
gₚₒ_FFᴹ = g_FF*bₘ^T*wₚₒ*POC
g_GOC_FFᴹ = g_FF*bₘ^T*w_GOC*GOC #29b
λₚₒ¹ = λ¹(T, O₂, bgc)

#println("Sum of positive terms is ", σᴹ*(∑gᴹ + ∑g_FFᴹ)*M + rᴹ*bₘ^T*K_mondo(M, Kₘ)*M + Pᵤₚᴹ + 0.5*R_CaCO₃*(mᴾ*K_mondo(P, Kₘ)*P + wᴾ*P^2) + 0.5*mᴰ*K_mondo(D, Kₘ)*D^3*wᴰ + Φ + Φ₂ᴰᴼᶜ)
#println("term a1 = $(σᴹ*(∑gᴹ)*M), term a2 = $(σᴹ*(gₚₒ_FFᴹ)*M), term a3 = $(σᴹ*g_GOC_FFᴹ*M), term b = $(rᴹ*bₘ^T*K_mondo(M, Kₘ)*M), term c = $(Pᵤₚᴹ), term d = $(0.5*R_CaCO₃*(mᴾ*K_mondo(P, Kₘ)*P + wᴾ*P^2)), term e = $(0.5*mᴰ*K_mondo(D, Kₘ)*D^2*wᴰ), term f = $(Φ), term g = $(Φ₂ᴰᴼᶜ)")
#println("Sum of negative terms is ", g_GOC_FFᴹ*M + λₚₒ¹*GOC)
#println("term 1 = $(g_GOC_FFᴹ*M), term 2 = $(λₚₒ¹*GOC)")
#println("-g_FF = $(g_FF), bₘ^T = $(bₘ^T), w_GOC = $(w_GOC), GOC = $(GOC) -")
#println("Total change is ", σᴹ*(∑gᴹ + ∑g_FFᴹ)*M + rᴹ*bₘ^T*K_mondo(M, Kₘ)*M + Pᵤₚᴹ + 0.5*R_CaCO₃*(mᴾ*K_mondo(P, Kₘ)*P + wᴾ*P^2) + 0.5*mᴰ*K_mondo(D, Kₘ)*D^2*wᴰ + Φ + Φ₂ᴰᴼᶜ - g_GOC_FFᴹ*M - λₚₒ¹*GOC )
#println("-------------------------------------")

return σᴹ*(∑gᴹ + ∑g_FFᴹ)*M + rᴹ*bₘ^T*(K_mondo(M, Kₘ) + 3*ΔO₂(O₂, bgc))*M + Pᵤₚᴹ + 0.5*R_CaCO₃*(mᴾ*K_mondo(P, Kₘ)*P + wᴾ*P^2) + 0.5*mᴰ*K_mondo(D, Kₘ)*D + D^2*wᴰ + Φ + Φ₂ᴰᴼᶜ - g_GOC_FFᴹ*M - λₚₒ¹*GOC #40 assumed that there is a typo in the D^2 instead of D^3
return σᴹ*(∑gᴹ + ∑g_FFᴹ)*M + rᴹ*bₘ^T*(K_mondo(M, Kₘ) + 3*ΔO₂(O₂, bgc))*M + Pᵤₚᴹ + 0.5*R_CaCO₃*(mᴾ*K_mondo(P, Kₘ)*P + sh*wᴾ*P^2) + 0.5*mᴰ*K_mondo(D, Kₘ)*D + sh*D^2*wᴰ + Φ + Φ₂ᴰᴼᶜ - g_GOC_FFᴹ*M - λₚₒ¹*GOC #40 assumed that there is a typo in the D^2 instead of D^3
end
15 changes: 10 additions & 5 deletions validation/PISCES/boxPISCES.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,13 +40,13 @@ PAR³ = FunctionField{Center, Center, Center}(PAR_func3, 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³))),
model = BoxModel(; biogeochemistry = PISCES(; grid, light_attenuation_model = PrescribedPhotosyntheticallyActiveRadiation((; PAR, PAR¹, PAR², PAR³)), flux_feeding_rate = 2.0e-3),
clock)
#set!(model, P = 0.0, D = 1.0, Z = 2.0, M = 3.0, Pᶜʰˡ= 4.0, Dᶜʰˡ = 5.0, Pᶠᵉ = 6.0, Dᶠᵉ = 7.0, Dˢⁱ = 8.0, DOC = 9.0, POC = 10.0, GOC = 11.0, SFe = 12.0, BFe = 13.0, PSi = 14.0, NO₃ = 15.0, NH₄ = 16.0, PO₄ = 17.0, Fe = 18.0, Si = 19.0, CaCO₃ = 20.0, DIC = 21.0, Alk = 22.0, O₂ = 23.0, T = 14.0)
set!(model, NO₃ = 4.0, NH₄ = 0.1, P = 4.26, D = 4.26, Z = 0.426, M = 0.426, Pᶠᵉ = 3.50, Dᶠᵉ = 3.50, Pᶜʰˡ = .1, Dᶜʰˡ = .01, Dˢⁱ = 0.525, Fe = 0.00082410, O₂ = 264.0, Si = 4.557, Alk = 2360.0, PO₄ = 0.8114, DIC = 2000.0, CaCO₃ = 0.0001, T = 14.0)
#set!(model,P = 3.963367728460601, D = 3.763831823528108, Z = 0.620887631503286, M = 0.4911996116700677, Pᶜʰˡ = 0.1263393104069646, Dᶜʰˡ = 0.0966272698878372, Pᶠᵉ = 2.916749891527781, Dᶠᵉ = 2.6966762460922764, Dˢⁱ = 0.5250058442518801, DOC = 5.492834645446811e-5, POC = 0.00010816947467085888, GOC = 1.541376629008023, SFe = 6.94778354330689e-5, BFe = 1.3780182342394662, PSi = 0.138718322180627, NO₃ = 3.862629483089866, NH₄ = 0.10480738012675432, PO₄ = 0.8031309301476024, Fe = 0.00024547654218086575, Si = 4.413896794698411, CaCO₃ = 0.011644257272404535, DIC = 1998.9796292207268, Alk = 2360.118267032333, O₂ = 265.37453137881016, T = 14.0)

#set!(model, NO₃ = 4.0, NH₄ = 0.1, P = 4.26, D = .426, Z = 0.426, M = 0.426, Pᶠᵉ = 3.50, Dᶠᵉ = 3.50, Pᶜʰˡ = .1, Dᶜʰˡ = .01, Dˢⁱ = 0.525, Fe = 0.00082410, O₂ = 264.0, Si = 4.557, Alk = 2360.0, PO₄ = 0.8114, DIC = 2000.0, CaCO₃ = 0.0001, T = 14.0)
set!(model, P = 3.4969507431335574, D = 3.045411522412732, Z = 0.06632989062825102, M = 0.2013934196815033, Pᶜʰˡ = 0.1649236496634812, Dᶜʰˡ = 0.1744075407918407, Pᶠᵉ = 0.629885753182277, Dᶠᵉ = 0.5811469701875243, Dˢⁱ = 0.6798967635852701, DOC = 0.00038765439868301653, POC = 0.0015526645806093746, GOC = 15.59902032438512, SFe = 0.0003450130201159402, BFe = 3.7215124823095893, PSi = 0.5404372397384611, NO₃ = 2.828750058046986, NH₄ = 0.5246932481313445, PO₄ = 0.7647402066361438, Fe = 2.240993210202666e-5, Si = 3.7394628646278316, CaCO₃ = 0.00933485588389991, DIC = 1994.2982703537236, Alk = 2361.5639068719747, O₂ = 272.5482524012977, T = 14.0)

simulation = Simulation(model; Δt = 5minutes, stop_time = 7hours)
simulation = Simulation(model; Δt = 5minutes, stop_time = 40minutes)

simulation.output_writers[:fields] = JLD2OutputWriter(model, model.fields; filename = "box.jld2", schedule = TimeInterval(5minutes), overwrite_existing = true)

Expand Down Expand Up @@ -100,4 +100,9 @@ fi = length(timeseries.P)

println("P = $(timeseries.P[fi]), D = $(timeseries.D[fi]), Z = $(timeseries.Z[fi]), M = $(timeseries.M[fi]), Pᶜʰˡ = $(timeseries.Pᶜʰˡ[fi]), Dᶜʰˡ = $(timeseries.Dᶜʰˡ[fi]), Pᶠᵉ = $(timeseries.Pᶠᵉ[fi]), Dᶠᵉ = $(timeseries.Dᶠᵉ[fi]), Dˢⁱ = $(timeseries.Dˢⁱ[fi]), DOC = $(timeseries.DOC[fi]), POC = $(timeseries.POC[fi]), GOC = $(timeseries.GOC[fi]), SFe = $(timeseries.SFe[fi]), BFe = $(timeseries.BFe[fi]), PSi = $(timeseries.PSi[fi]), NO₃ = $(timeseries.NO₃[fi]), NH₄ = $(timeseries.NH₄[fi]), PO₄ = $(timeseries.PO₄[fi]), Fe = $(timeseries.Fe[fi]), Si = $(timeseries.Si[fi]), CaCO₃ = $(timeseries.CaCO₃[fi]), DIC = $(timeseries.DIC[fi]), Alk = $(timeseries.Alk[fi]), O₂ = $(timeseries.O₂[fi]), T = 14.0")

Carbon_at_start = timeseries.P[1] + timeseries.D[1] + timeseries.Z[1] + timeseries.M[1] + timeseries.DOC[1] + timeseries.POC[1] + timeseries.GOC[1] + timeseries.DIC[1] + timeseries.CaCO₃[1]
Carbon_at_end = timeseries.P[fi] + timeseries.D[fi] + timeseries.Z[fi] + timeseries.M[fi] + timeseries.DOC[fi] + timeseries.POC[fi] + timeseries.GOC[fi] + timeseries.DIC[fi] + timeseries.CaCO₃[fi]

println("Carbon at start = ", Carbon_at_start, "Carbon at end = ", Carbon_at_end)

fig
Loading