From 4df7c8ea6da728828bd4f9641fbf00ce160703f9 Mon Sep 17 00:00:00 2001 From: ciadht Date: Sat, 24 Aug 2024 14:48:42 +0100 Subject: [PATCH 1/5] Changed PAR names and added chlorophyll function --- Manifest.toml | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 8b45b3531..8c9034d63 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,8 +1,8 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.2" +julia_version = "1.10.4" manifest_format = "2.0" -project_hash = "bd1de93833dafc8fa6de30ae101c862f7fa5877e" +project_hash = "37e11a3bdd396c973917441db34ded05b97faa85" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] @@ -171,7 +171,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.1.0+0" +version = "1.1.1+0" [[deps.CompositionsBase]] git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" @@ -443,20 +443,16 @@ uuid = "9c1d0b0a-7046-5b2e-a33f-ea22f176ac7e" version = "0.2.1+0" [[deps.KernelAbstractions]] -deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] -git-tree-sha1 = "35ceea58aa34ad08b1ae00a52622c62d1cfb8ce2" +deps = ["Adapt", "Atomix", "InteractiveUtils", "LinearAlgebra", "MacroTools", "PrecompileTools", "Requires", "SparseArrays", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"] +git-tree-sha1 = "0fac59881e91c7233a9b0d47f4b7d9432e534f0f" uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -version = "0.9.24" +version = "0.9.23" [deps.KernelAbstractions.extensions] EnzymeExt = "EnzymeCore" - LinearAlgebraExt = "LinearAlgebra" - SparseArraysExt = "SparseArrays" [deps.KernelAbstractions.weakdeps] EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" - LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" - SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[deps.LLVM]] deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"] @@ -673,9 +669,9 @@ version = "1.2.0" [[deps.Oceananigans]] deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"] -git-tree-sha1 = "12cdd648cc342e52686515811c69bc4bc452a94c" +git-tree-sha1 = "ed415deb1d80c2a15c9b8bf0c7df04c6dec9d6c3" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" -version = "0.91.9" +version = "0.91.8" [deps.Oceananigans.extensions] OceananigansEnzymeExt = "Enzyme" From e72c87bebcc31a9fedf27b4c52e3d187ba110cfa Mon Sep 17 00:00:00 2001 From: ciadht Date: Mon, 26 Aug 2024 15:45:48 +0100 Subject: [PATCH 2/5] Fixed bugs in columnPISCES --- .../AdvectedPopulations/PISCES/PISCES.jl | 2 +- .../AdvectedPopulations/PISCES/calcite.jl | 2 +- validation/PISCES/columnPISCES.jl | 35 +++++++------------ 3 files changed, 15 insertions(+), 24 deletions(-) diff --git a/src/Models/AdvectedPopulations/PISCES/PISCES.jl b/src/Models/AdvectedPopulations/PISCES/PISCES.jl index d73ba10b0..d0ae6ddb7 100644 --- a/src/Models/AdvectedPopulations/PISCES/PISCES.jl +++ b/src/Models/AdvectedPopulations/PISCES/PISCES.jl @@ -950,5 +950,5 @@ include("zooplankton.jl") @inline sinking_tracers(::PISCES) = (:POC, :GOC, :SFe, :BFe, :PSi, :CaCO₃) # please list them here -@inline chlorophyll(model) = model.tracers.Pᶜʰˡ + model.tracers.Dᶜʰˡ +@inline chlorophyll(bgc::PISCES, model) = model.tracers.Pᶜʰˡ + model.tracers.Dᶜʰˡ end # module diff --git a/src/Models/AdvectedPopulations/PISCES/calcite.jl b/src/Models/AdvectedPopulations/PISCES/calcite.jl index 0dea836a4..dd0ae6f09 100644 --- a/src/Models/AdvectedPopulations/PISCES/calcite.jl +++ b/src/Models/AdvectedPopulations/PISCES/calcite.jl @@ -39,7 +39,7 @@ end wᴾ = bgc.min_quadratic_mortality_of_phytoplankton ηᶻ = bgc.fraction_of_calcite_not_dissolving_in_guts.Z ηᴹ = bgc.fraction_of_calcite_not_dissolving_in_guts.M - sh = shear_rate(z, zₘₓₗ)PAR₃ + sh = shear_rate(z, zₘₓₗ) 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 diff --git a/validation/PISCES/columnPISCES.jl b/validation/PISCES/columnPISCES.jl index 53c74b1dc..b50321ab9 100644 --- a/validation/PISCES/columnPISCES.jl +++ b/validation/PISCES/columnPISCES.jl @@ -23,7 +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⁰(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 * (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) @@ -39,28 +39,14 @@ nothing #hide @inline euphotic(x, y, z, t) = - 50.0 nothing #hide -PAR_func(x, y, z, t) = PAR⁰(x,y,t)*exp(z/10) # Modify the PAR based on the nominal depth and exponential decay -PAR_func1(x, y, z, t) = 1/3*PAR⁰(x,y,t)*exp(z/10) -PAR_func2(x, y, z, t) = 1/3*PAR⁰(x,y,t)*exp(z/10) -PAR_func3(x, y, z, t) = 1/3*PAR⁰(x,y,t)*exp(z/10) - -yearly_maximum_silicate = ConstantField(1) -dust_deposition = ConstantField(0) -carbonate_sat_ratio = ConstantField(0) - -#w_GOC(z) = 30/day + (200/day - 30/day)*(max(0, abs(z)-100))/(5000) +#w_GOC(z) = 30/day + (200/day - 30/day)*(max(0, abs(z)-abs(zₘₓₗ)))/(5000) w_GOC = 30/day w_POC = 2.0/day grid = RectilinearGrid(size = (1, 1, 100), extent = (20meters, 20meters, 400meters)) 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) zₘₓₗ = FunctionField{Center, Center, Center}(MLD, grid; clock) zₑᵤ = FunctionField{Center, Center, Center}(euphotic, grid; clock) @@ -77,7 +63,7 @@ zₑᵤ = FunctionField{Center, Center, Center}(euphotic, grid; clock) biogeochemistry = PISCES(; grid, 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ₑᵤ) + mixed_layer_depth = zₘₓₗ, euphotic_layer_depth = zₑᵤ) CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition() @@ -89,17 +75,19 @@ model = NonhydrostaticModel(; grid, clock, closure = ScalarDiffusivity(VerticallyImplicitTimeDiscretization(), κ = κₜ), biogeochemistry, - boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), )) - + boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux), ), + auxiliary_fields = (; S) + ) + @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.) +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.001, DIC = 2139.0, Alk = 2366.0, O₂ = 237.0, T = funT) #Using Copernicus Data (26.665, 14.) # ## Simulation # Next we setup a simulation and add some callbacks that: # - Show the progress of the simulation # - Store the model and particles output -simulation = Simulation(model, Δt = 90minutes, stop_time = 2years) +simulation = Simulation(model, Δt = 50minutes, stop_time = 2years) progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", iteration(sim), @@ -179,7 +167,7 @@ carbon_export = zeros(length(times)) using Oceananigans.Biogeochemistry: biogeochemical_drift_velocity for (i, t) in enumerate(times) - air_sea_CO₂_flux[i] = CO₂_flux.condition.func(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) + #air_sea_CO₂_flux[i] = CO₂_flux.condition.func(0.0, 0.0, t, DIC[1, 1, grid.Nz, i], Alk[1, 1, grid.Nz, i], temp(1, 1, 0, t), 35) carbon_export[i] = (POC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:POC)).w[1, 1, grid.Nz-20] + GOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(model.biogeochemistry, Val(:GOC)).w[1, 1, grid.Nz-20]) * redfield(Val(:GOC), model.biogeochemistry) end @@ -299,4 +287,7 @@ lines!(axfDIC, times / days, air_sea_CO₂_flux / 1e3 * CO₂_molar_mass * year, lines!(axfDIC, times / days, carbon_export / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export") Legend(fig[7, 2], axfDIC, framevisible = false) +axs = [] +push!(axs, Axis(fig[7,3], xlabel = "Time (days)", title = "Mixed Layer Depth (m)")) +lines!(axs[end], (0:1day:2years)/days, x -> MLD(0, 0, 0, x*days), linewidth = 3) fig From ff959df6ea12f9ee89f39e33d68320db31e6ec49 Mon Sep 17 00:00:00 2001 From: ciadht Date: Wed, 4 Sep 2024 12:09:53 +0100 Subject: [PATCH 3/5] Commenting changes --- validation/PISCES/columnPISCES.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/validation/PISCES/columnPISCES.jl b/validation/PISCES/columnPISCES.jl index b50321ab9..bfc645505 100644 --- a/validation/PISCES/columnPISCES.jl +++ b/validation/PISCES/columnPISCES.jl @@ -21,7 +21,7 @@ const year = years = 365days 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) +# Setting up idealised functions for PAR and diffusivity (details here can be ignored but these are typical of the North Atlantic), temperaeture and euphotic layer @inline PAR⁰(x,y,t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2 @@ -33,13 +33,12 @@ nothing #hide @inline κₜ(x, y, z, t) = (1e-2 * (1 + tanh((z - MLD(x, y, z, t)) / 10)) / 2 + 1e-4) -#@inline temp(x, y, z, t) = (2.4 * cos(t * 2π / year + 50days) + 10)*exp(z/50) -@inline temp(x, y, z, t) = 14*exp(z/10) +@inline temp(x, y, z, t) = (2.4 * cos(t * 2π / year + 50days) + 10)*exp(z/10) @inline euphotic(x, y, z, t) = - 50.0 nothing #hide - +#The commented equation is the correct form of w_GOC, but have not figured out how to implement this #w_GOC(z) = 30/day + (200/day - 30/day)*(max(0, abs(z)-abs(zₘₓₗ)))/(5000) w_GOC = 30/day w_POC = 2.0/day @@ -51,9 +50,6 @@ 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. - # ## Model # First we define the biogeochemical model including carbonate chemistry (for which we also define temperature (``T``) and salinity (``S``) fields) @@ -287,6 +283,7 @@ lines!(axfDIC, times / days, air_sea_CO₂_flux / 1e3 * CO₂_molar_mass * year, lines!(axfDIC, times / days, carbon_export / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export") Legend(fig[7, 2], axfDIC, framevisible = false) +#Plotting a graph of Mixed Layer Depth axs = [] push!(axs, Axis(fig[7,3], xlabel = "Time (days)", title = "Mixed Layer Depth (m)")) lines!(axs[end], (0:1day:2years)/days, x -> MLD(0, 0, 0, x*days), linewidth = 3) From a322b21efc1a3553b87a9efa769d126518d028b7 Mon Sep 17 00:00:00 2001 From: ciadht Date: Wed, 4 Sep 2024 12:11:10 +0100 Subject: [PATCH 4/5] Calcite comment --- validation/PISCES/columnPISCES.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/validation/PISCES/columnPISCES.jl b/validation/PISCES/columnPISCES.jl index bfc645505..41bccb5b8 100644 --- a/validation/PISCES/columnPISCES.jl +++ b/validation/PISCES/columnPISCES.jl @@ -77,7 +77,7 @@ model = NonhydrostaticModel(; grid, @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.202, NH₄ = 0.25*6.202, PO₄ = 0.8722, Fe = 1.256, Si = 7.313, CaCO₃ = 0.001, DIC = 2139.0, Alk = 2366.0, O₂ = 237.0, T = funT) #Using Copernicus Data (26.665, 14.) +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.001, DIC = 2139.0, Alk = 2366.0, O₂ = 237.0, T = funT) #Using Copernicus Data (26.665, 14.), Calcite is not correct, but this is to see it on the graphs # ## Simulation # Next we setup a simulation and add some callbacks that: # - Show the progress of the simulation From a0be5646aa05ec5f300ee0503036a7f075f669a1 Mon Sep 17 00:00:00 2001 From: ciadht Date: Wed, 4 Sep 2024 12:30:07 +0100 Subject: [PATCH 5/5] Commenting BoxPISCES --- validation/PISCES/boxPISCES.jl | 28 +++++++++++----------------- validation/PISCES/columnPISCES.jl | 1 + 2 files changed, 12 insertions(+), 17 deletions(-) diff --git a/validation/PISCES/boxPISCES.jl b/validation/PISCES/boxPISCES.jl index 495c5a206..5ce2b0517 100644 --- a/validation/PISCES/boxPISCES.jl +++ b/validation/PISCES/boxPISCES.jl @@ -18,21 +18,17 @@ using Oceananigans.Fields: FunctionField const year = years = 365day nothing #hide -grid = RectilinearGrid( topology = (Flat, Flat, Flat), size = (), z = -100) +grid = RectilinearGrid( topology = (Flat, Flat, Flat), size = (), z = 0) clock = Clock(time = zero(grid)) # This is forced by a prescribed time-dependent photosynthetically available radiation (PAR) PAR_func(t) = 60 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2 MLD(t) = - 100 EU(t) = - 50 - # specify the nominal depth of the box for the PAR profile -# Modify the PAR based on the nominal depth and exponential decay -#PAR_func(t) = 18.0 # Modify the PAR based on the nominal depth and exponential decay -#PAR_func(t) = 18.0 PAR_func1(t) = PAR_func(t)/3 PAR_func2(t) = PAR_func(t)/3 PAR_func3(t)= PAR_func(t)/3 - +#Each PAR component must sum to the initial PAR function at the surface, and we are only considering the 0 layer PAR₁ = FunctionField{Center, Center, Center}(PAR_func1, grid; clock) PAR₂ = FunctionField{Center, Center, Center}(PAR_func2, grid; clock) @@ -42,8 +38,8 @@ 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.) @@ -54,12 +50,8 @@ simulation.output_writers[:fields] = JLD2OutputWriter(model, model.fields; filen prog(sim) = @info "$(prettytime(time(sim))) in $(prettytime(simulation.run_wall_time))" +#NaN checker function, could be removed if confident of model stability function non_zero_fields!(model) - #for tracer in model.fields - # if tracer != model.fields[:T] - # tracer = tracer[1,1,1] - # getproperty(model, tracer) = max(0.0, getproperty(model, tracer)) - #end @inbounds for (idx, field) in enumerate(model.fields) if isnan(field[1,1,1]) throw("$(keys(model.fields)[idx]) has gone NaN") @@ -68,9 +60,7 @@ function non_zero_fields!(model) end end - #end return nothing - end simulation.callbacks[:progress] = Callback(prog, IterationInterval(100)) @@ -97,11 +87,13 @@ for (name, tracer) in pairs(timeseries) push!(axs, Axis(fig[floor(Int, idx/4), Int(idx%4)], ylabel = "$name", xlabel = "years", xticks=(0:40))) lines!(axs[end], times / year, tracer, linewidth = 3) end + +#Plotting the function of PAR push!(axs, Axis(fig[6,1], ylabel = "PAR", xlabel = "years", xticks=(0:40))) lines!(axs[end], (0:10day:5years) / year, x -> PAR_func(x * year), linewidth = 3) -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") +#Below is merely a tester to check that values of things are conserved, these can be removed in the test files +fi = length(timeseries.P) 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] @@ -119,4 +111,6 @@ println("Iron at start = ", Iron_at_start, " Iron at end = ", Iron_at_end) println("Silicon at start = ", Silicon_at_start, " Silicon at end = ", Silicon_at_End) println("Phosphates at start = ", Phosphates_at_start, " Phosphates at end = ", Phosphates_at_end) println("Nitrogen at start = ", Nitrogen_at_start, " Nitrogen at end = ", Nitrogen_at_end) + +#Display the figure fig \ No newline at end of file diff --git a/validation/PISCES/columnPISCES.jl b/validation/PISCES/columnPISCES.jl index 41bccb5b8..de718fb49 100644 --- a/validation/PISCES/columnPISCES.jl +++ b/validation/PISCES/columnPISCES.jl @@ -94,6 +94,7 @@ progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(10day) ) +#NaN Checker function. Could be removed to improve speed, if confident of model stability function non_zero_fields!(model) @inbounds for (idx, tracer) in enumerate(model.tracers) for i in 1:50