diff --git a/src/Models/AdvectedPopulations/PISCES/PISCES.jl b/src/Models/AdvectedPopulations/PISCES/PISCES.jl index 4b980ba5d..be53b7b07 100644 --- a/src/Models/AdvectedPopulations/PISCES/PISCES.jl +++ b/src/Models/AdvectedPopulations/PISCES/PISCES.jl @@ -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) @@ -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 diff --git a/src/Models/AdvectedPopulations/PISCES/POC_and_GOC.jl b/src/Models/AdvectedPopulations/PISCES/POC_and_GOC.jl index c06632d45..5466557e7 100644 --- a/src/Models/AdvectedPopulations/PISCES/POC_and_GOC.jl +++ b/src/Models/AdvectedPopulations/PISCES/POC_and_GOC.jl @@ -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³) @@ -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 \ No newline at end of file diff --git a/validation/PISCES/boxPISCES.jl b/validation/PISCES/boxPISCES.jl index d9b16cd61..a840880e4 100644 --- a/validation/PISCES/boxPISCES.jl +++ b/validation/PISCES/boxPISCES.jl @@ -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) @@ -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 \ No newline at end of file