From 6575e4d1fbe977406f0dd9436ccfa90de7dab1c3 Mon Sep 17 00:00:00 2001 From: sichen Date: Thu, 18 Jan 2024 15:19:02 +0000 Subject: [PATCH 1/5] add subpolarnew example folder --- subpolarnew/kelp450_11_new/kelp.jl | 215 ++++++++++++++++++ subpolarnew/kelp450_11_new/kelp450_11.jl | 172 ++++++++++++++ subpolarnew/kelp450_11_new/kelp450_11_new.jl | 189 +++++++++++++++ subpolarnew/kelp450_11_new/kelp450_11_new2.jl | 189 +++++++++++++++ .../kelp450_11_new/kelp_old_example.jl | 206 +++++++++++++++++ test_stretch/readme.rtf | 17 ++ 6 files changed, 988 insertions(+) create mode 100644 subpolarnew/kelp450_11_new/kelp.jl create mode 100644 subpolarnew/kelp450_11_new/kelp450_11.jl create mode 100644 subpolarnew/kelp450_11_new/kelp450_11_new.jl create mode 100644 subpolarnew/kelp450_11_new/kelp450_11_new2.jl create mode 100644 subpolarnew/kelp450_11_new/kelp_old_example.jl create mode 100644 test_stretch/readme.rtf diff --git a/subpolarnew/kelp450_11_new/kelp.jl b/subpolarnew/kelp450_11_new/kelp.jl new file mode 100644 index 000000000..925530d66 --- /dev/null +++ b/subpolarnew/kelp450_11_new/kelp.jl @@ -0,0 +1,215 @@ +# # Simple active particle example +# Here, we setup a simple 1D column example with the [LOBSTER](@ref LOBSTER) biogeochemical model and active particles modelling the growth of sugar kelp. +# This example demonstrates: +# - How to setup OceanBioME's biogeochemical models +# - How to add biologically active particles which interact with the biodeochemical model +# - How to visualise results + +# This is forced by idealised mixing layer depth and surface photosynthetically available radiation (PAR) which are setup first. + +# ## Install dependencies +# First we check we have the dependencies installed +# ```julia +# using Pkg +# pkg "add OceanBioME, Oceananigans, CairoMakie, JLD2" +# ``` + +# ## Model setup +# First load the required packages +using OceanBioME, Oceananigans, Printf +using Oceananigans.Fields: FunctionField, ConstantField +using Oceananigans.Units +using Oceananigans.Architectures: arch_array + +const year = years = 365days # just for these idealised cases +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 H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0) + +@inline fmld1(t) = H(t, 50days, year) * (1 / (1 + exp(-(t - 100days) / 5days))) * (1 / (1 + exp((t - 330days) / 25days))) + +@inline MLD(t) = - (10 + 340 * (1 - fmld1(year - eps(year)) * exp(-mod(t, year) / 25days) - fmld1(mod(t, year)))) + +@inline κₜ(x, y, z, t) = 1e-2 * (1 + tanh((z - MLD(t))/10)) / 2 + 1e-4 + +@inline temp(x, y, z, t) = 2.4 * cos(t * 2π / year + 50day) + 10 +nothing #hide + +architecture = CPU() + +# ## Grid and PAR field +# Define the grid and an extra Oceananigans' field that the PAR will be stored in +Lx, Ly = 20meters, 20meters +grid = RectilinearGrid(architecture, size=(1, 1, 50), extent=(Lx, Ly, 200)) + +# Specify the boundary conditions for DIC and O₂ based on the air-sea CO₂ and O₂ flux +CO₂_flux = GasExchange(; gas = :CO₂) + +clock = Clock(; time = 0.0) +T = FunctionField{Center, Center, Center}(temp, grid; clock) +S = ConstantField(35) + +# ## Kelp Particle setup +@info "Setting up kelp particles" + +n = 5 # number of kelp fronds +z₀ = [-21:5:-1;] * 1.0 # depth of kelp fronds + +particles = SLatissima(; architecture, + x = arch_array(architecture, ones(n) * Lx / 2), + y = arch_array(architecture, ones(n) * Ly / 2), + z = arch_array(architecture, z₀), + A = arch_array(architecture, ones(n) * 10.0), + latitude = 57.5, + scalefactor = 500.0) + +# ## Setup BGC model +biogeochemistry = LOBSTER(; grid, + surface_photosynthetically_active_radiation = PAR⁰, + carbonates = true, + variable_redfield = true, + scale_negatives = true, + particles) + +model = NonhydrostaticModel(; grid, + clock, + closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), + biogeochemistry, + auxiliary_fields = (; T, S)) + +set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2409.0) + +# ## Simulation +# Next we setup the simulation along with some callbacks that: +# - Show the progress of the simulation +# - Store the model and particles output +# - Prevent the tracers from going negative from numerical error (see discussion of this in the [positivity preservation](@ref pos-preservation) implementation page) + +simulation = Simulation(model, Δt = 3minutes, stop_time = 100days) + +progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) + +simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(10days)) + +filename = "kelp" +simulation.output_writers[:profiles] = JLD2OutputWriter(model, model.tracers, + filename = "$filename.jld2", + schedule = TimeInterval(1day), + overwrite_existing = true) + +simulation.output_writers[:particles] = JLD2OutputWriter(model, (; particles), + filename = "$(filename)_particles.jld2", + schedule = TimeInterval(1day), + overwrite_existing = true) + +nothing #hide + +# ## Run! +# Finally we run the simulation +run!(simulation) + +# Now we can visualise the results with some post processing to diagnose the air-sea CO₂ flux - hopefully this looks different to the example without kelp! + + P = FieldTimeSeries("$filename.jld2", "P") + NO₃ = FieldTimeSeries("$filename.jld2", "NO₃") + Z = FieldTimeSeries("$filename.jld2", "Z") +sPOC = FieldTimeSeries("$filename.jld2", "sPOC") +bPOC = FieldTimeSeries("$filename.jld2", "bPOC") + DIC = FieldTimeSeries("$filename.jld2", "DIC") + Alk = FieldTimeSeries("$filename.jld2", "Alk") + +x, y, z = nodes(P) +times = P.times +nothing #hide + +# We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and +# the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth +# that corresponds to `k = grid.Nz - 20`. +air_sea_CO₂_flux = zeros(length(times)) +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) + carbon_export[i] = sPOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:sPOC)).w[1, 1, grid.Nz-20] + + bPOC[1, 1, grid.Nz-20, i] * biogeochemical_drift_velocity(biogeochemistry, Val(:bPOC)).w[1, 1, grid.Nz-20] +end + +# Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. + +using CairoMakie + +fig = Figure(resolution = (1000, 1500), fontsize = 20) + +axis_kwargs = (xlabel = "Time (days)", ylabel = "z (m)", limits = ((0, times[end] / days), (-85meters, 0))) + +axP = Axis(fig[1, 1]; title = "Phytoplankton concentration (mmol N/m³)", axis_kwargs...) +hmP = heatmap!(times / days, z, interior(P, 1, 1, :, :)', colormap = :batlow) +Colorbar(fig[1, 2], hmP) + +axNO₃ = Axis(fig[2, 1]; title = "Nitrate concentration (mmol N/m³)", axis_kwargs...) +hmNO₃ = heatmap!(times / days, z, interior(NO₃, 1, 1, :, :)', colormap = :batlow) +Colorbar(fig[2, 2], hmNO₃) + +axZ = Axis(fig[3, 1]; title = "Zooplankton concentration (mmol N/m³)", axis_kwargs...) +hmZ = heatmap!(times / days, z, interior(Z, 1, 1, :, :)', colormap = :batlow) +Colorbar(fig[3, 2], hmZ) + +axD = Axis(fig[4, 1]; title = "Detritus concentration (mmol C/m³)", axis_kwargs...) +hmD = heatmap!(times / days, z, interior(sPOC, 1, 1, :, :)' .+ interior(bPOC, 1, 1, :, :)', colormap = :batlow) +Colorbar(fig[4, 2], hmD) + +CO₂_molar_mass = (12 + 2 * 16) * 1e-3 # kg / mol + +axfDIC = Axis(fig[5, 1], xlabel = "Time (days)", ylabel = "Flux (kgCO₂/m²/year)", + title = "Air-sea CO₂ flux and Sinking", limits = ((0, times[end] / days), nothing)) +lines!(axfDIC, times / days, air_sea_CO₂_flux / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Air-sea flux") +lines!(axfDIC, times / days, carbon_export / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export") +Legend(fig[5, 2], axfDIC, framevisible = false) + +fig + +# We can also have a look at how the kelp particles evolve +using JLD2 + +file = jldopen("$(filename)_particles.jld2") + +iterations = keys(file["timeseries/t"]) + +A, N, C = ntuple(n -> zeros(5, length(iterations)), 3) + +times = zeros(length(iterations)) + +for (i, iter) in enumerate(iterations) + particles = file["timeseries/particles/$iter"] + A[:, i] = particles.A + N[:, i] = particles.N + C[:, i] = particles.C + + times[i] = file["timeseries/t/$iter"] +end + +fig = Figure(resolution = (1000, 800), fontsize = 20) + +axis_kwargs = (xlabel = "Time (days)", limits = ((0, times[end] / days), nothing)) + +ax1 = Axis(fig[1, 1]; ylabel = "Frond area (dm²)", axis_kwargs...) +[lines!(ax1, times / day, A[n, :], linewidth = 3) for n in 1:5] + +ax2 = Axis(fig[2, 1]; ylabel = "Total Carbon (gC)", axis_kwargs...) +[lines!(ax2, times / day, (@. A * (N + particles.structural_nitrogen) * particles.structural_dry_weight_per_area)[n, :], linewidth = 3) for n in 1:5] + +ax3 = Axis(fig[3, 1]; ylabel = "Total Nitrogen (gN)", axis_kwargs...) +[lines!(ax3, times / day, (@. A * (C + particles.structural_carbon) * particles.structural_dry_weight_per_area)[n, :], linewidth = 3) for n in 1:5] + +fig diff --git a/subpolarnew/kelp450_11_new/kelp450_11.jl b/subpolarnew/kelp450_11_new/kelp450_11.jl new file mode 100644 index 000000000..ca1ec3b60 --- /dev/null +++ b/subpolarnew/kelp450_11_new/kelp450_11.jl @@ -0,0 +1,172 @@ +# # Simple active particle example +# In this example we will setup a simple 1D column with the [LOBSTER](@ref LOBSTER) biogeochemical model and active particles modelling the growth of sugar kelp. This demonstraits: +# - How to setup OceanBioME's biogeochemical models +# - How to setup light attenuation +# - How to add biologically active particles which interact with the biodeochemical model +# - How to include optional tracer sets (carbonate chemistry and oxygen) +# - How to visulise results + +# This is forced by idealised mixing layer depth and surface photosynthetically available radiation (PAR) which are setup first + +# ## Install dependencies +# First we will check we have the dependencies installed +# ```julia +# using Pkg +# pkg"add OceanBioME, Oceananigans, Printf, GLMakie" +# ``` + +# ## Model setup +# We load the packages and choose the default LOBSTER parameter set +using OceanBioME, Oceananigans,Printf +using Oceananigans.Units: second, minute, minutes, hour, hours, day, days, year, years +using Interpolations +using JLD2 + +# ## 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) + +##PAR⁰(x, y, t) = 60*(1-cos((t+15days)*2π/(365days)))*(1 /(1 +0.2*exp(-((mod(t, 365days)-200days)/50days)^2))) .+ 2 +##H(t, t₀, t₁) = ifelse(t₀ keys(pp["timeseries"]) 15-element Vector{String}: "NO₃" "NH₄" "P" "Z" "sPON" "bPON" "DON" "DIC" "Alk" "O₂" "sPOC" "bPOC" "DOC" "PAR" "t"# + +# ## Simulation +# Next we setup the simulation along with some callbacks that: +# - Couples the particles to the biodeochemical model +# - Update the PAR field from the surface PAR and phytoplankton concentration +# - Show the progress of the simulation +# - Store the model and particles output +# - Prevent the tracers from going negative from numerical error (see discussion of this in the [positivity preservation](@ref pos-preservation) implimentation page) + +simulation = Simulation(model, Δt=3.5minutes, stop_time=duration) + +simulation.callbacks[:couple_particles] = Callback(Particles.infinitesimal_particle_field_coupling!; callsite = TendencyCallsite()) + +progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) +simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100)) + +filename = "kelp450_11" +simulation.output_writers[:profiles] = JLD2OutputWriter(model, merge(model.tracers, model.auxiliary_fields), filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing=true) +simulation.output_writers[:particles] = JLD2OutputWriter(model, (particles=model.particles, ), filename = "$(filename)_particles.jld2", schedule = TimeInterval(1day), overwrite_existing = true) + + +#simulation.callbacks[:timestep] = Callback(update_timestep!, IterationInterval(1), (c_forcing=0.1, c_adv=0.5, c_diff=0.5, w = 200/day, relaxation=0.95), TimeStepCallsite()) + +scale_negative_tracers = ScaleNegativeTracers(tracers = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON)) +simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) +plankton_redfield = model.biogeochemistry.phytoplankton_redfield +scale_negative_carbon_tracers = ScaleNegativeTracers(tracers = (:P, :Z, :DOC, :sPOC, :bPOC, :DIC), scalefactors = (P = plankton_redfield, Z = plankton_redfield, DOC = 1, sPOC = 1, bPOC = 1, DIC = 1)) +simulation.callbacks[:neg2] = Callback(scale_negative_carbon_tracers; callsite = UpdateStateCallsite()) +#wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 2.0, min_change = 0.5, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) +#simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) + +#simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=SpecifiedTimes([i*years for i=1:12]), prefix=filename) #prefix="kelp_01" +# ## Run! +# Finally we run the simulation +run!(simulation) + + + + + + + + + + + + + + + + + + diff --git a/subpolarnew/kelp450_11_new/kelp450_11_new.jl b/subpolarnew/kelp450_11_new/kelp450_11_new.jl new file mode 100644 index 000000000..966f0702d --- /dev/null +++ b/subpolarnew/kelp450_11_new/kelp450_11_new.jl @@ -0,0 +1,189 @@ +# # Simple active particle example +# In this example we will setup a simple 1D column with the [LOBSTER](@ref LOBSTER) biogeochemical model and active particles modelling the growth of sugar kelp. This demonstraits: +# - How to setup OceanBioME's biogeochemical models +# - How to setup light attenuation +# - How to add biologically active particles which interact with the biodeochemical model +# - How to include optional tracer sets (carbonate chemistry and oxygen) +# - How to visulise results + +# This is forced by idealised mixing layer depth and surface photosynthetically available radiation (PAR) which are setup first + +# ## Install dependencies +# First we will check we have the dependencies installed +# ```julia +# using Pkg +# pkg"add OceanBioME, Oceananigans, Printf, GLMakie" +# ``` + +# ## Model setup +# We load the packages and choose the default LOBSTER parameter set +using OceanBioME, Oceananigans,Printf +using Oceananigans.Units: second, minute, minutes, hour, hours, day, days#, year, years +using Interpolations +using JLD2 +using Oceananigans.Architectures: arch_array #####changed +const year = years = 365days #####changed + +# ## 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) + +##PAR⁰(x, y, t) = 60*(1-cos((t+15days)*2π/(365days)))*(1 /(1 +0.2*exp(-((mod(t, 365days)-200days)/50days)^2))) .+ 2 +##H(t, t₀, t₁) = ifelse(t₀ keys(pp["timeseries"]) 15-element Vector{String}: "NO₃" "NH₄" "P" "Z" "sPON" "bPON" "DON" "DIC" "Alk" "O₂" "sPOC" "bPOC" "DOC" "PAR" "t"# + +# ## Simulation +# Next we setup the simulation along with some callbacks that: +# - Couples the particles to the biodeochemical model +# - Update the PAR field from the surface PAR and phytoplankton concentration +# - Show the progress of the simulation +# - Store the model and particles output +# - Prevent the tracers from going negative from numerical error (see discussion of this in the [positivity preservation](@ref pos-preservation) implimentation page) + +simulation = Simulation(model, Δt=3.5minutes, stop_time=duration) + +# simulation.callbacks[:couple_particles] = Callback(Particles.infinitesimal_particle_field_coupling!; callsite = TendencyCallsite()) + +progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) +simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100)) + +filename = "kelp450_11_new" +simulation.output_writers[:profiles] = JLD2OutputWriter(model, merge(model.tracers, model.auxiliary_fields), filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing=true) +simulation.output_writers[:particles] = JLD2OutputWriter(model, (; kelp_particles), filename = "$(filename)_particles.jld2", schedule = TimeInterval(1day), overwrite_existing = true) + + +#simulation.callbacks[:timestep] = Callback(update_timestep!, IterationInterval(1), (c_forcing=0.1, c_adv=0.5, c_diff=0.5, w = 200/day, relaxation=0.95), TimeStepCallsite()) + +scale_negative_tracers = ScaleNegativeTracers(;model, tracers = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON)) #changed +simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) + +# plankton_redfield = model.biogeochemistry.phytoplankton_redfield +# scale_negative_carbon_tracers = ScaleNegativeTracers(tracers = (:P, :Z, :DOC, :sPOC, :bPOC, :DIC), scalefactors = (P = plankton_redfield, Z = plankton_redfield, DOC = 1, sPOC = 1, bPOC = 1, DIC = 1)) +# simulation.callbacks[:neg2] = Callback(scale_negative_carbon_tracers; callsite = UpdateStateCallsite()) + + +#wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 2.0, min_change = 0.5, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) +#simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) + +#simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=SpecifiedTimes([i*years for i=1:12]), prefix=filename) #prefix="kelp_01" +# ## Run! +# Finally we run the simulation +run!(simulation) + + + + + + + + + + + + + + + + + + diff --git a/subpolarnew/kelp450_11_new/kelp450_11_new2.jl b/subpolarnew/kelp450_11_new/kelp450_11_new2.jl new file mode 100644 index 000000000..966f0702d --- /dev/null +++ b/subpolarnew/kelp450_11_new/kelp450_11_new2.jl @@ -0,0 +1,189 @@ +# # Simple active particle example +# In this example we will setup a simple 1D column with the [LOBSTER](@ref LOBSTER) biogeochemical model and active particles modelling the growth of sugar kelp. This demonstraits: +# - How to setup OceanBioME's biogeochemical models +# - How to setup light attenuation +# - How to add biologically active particles which interact with the biodeochemical model +# - How to include optional tracer sets (carbonate chemistry and oxygen) +# - How to visulise results + +# This is forced by idealised mixing layer depth and surface photosynthetically available radiation (PAR) which are setup first + +# ## Install dependencies +# First we will check we have the dependencies installed +# ```julia +# using Pkg +# pkg"add OceanBioME, Oceananigans, Printf, GLMakie" +# ``` + +# ## Model setup +# We load the packages and choose the default LOBSTER parameter set +using OceanBioME, Oceananigans,Printf +using Oceananigans.Units: second, minute, minutes, hour, hours, day, days#, year, years +using Interpolations +using JLD2 +using Oceananigans.Architectures: arch_array #####changed +const year = years = 365days #####changed + +# ## 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) + +##PAR⁰(x, y, t) = 60*(1-cos((t+15days)*2π/(365days)))*(1 /(1 +0.2*exp(-((mod(t, 365days)-200days)/50days)^2))) .+ 2 +##H(t, t₀, t₁) = ifelse(t₀ keys(pp["timeseries"]) 15-element Vector{String}: "NO₃" "NH₄" "P" "Z" "sPON" "bPON" "DON" "DIC" "Alk" "O₂" "sPOC" "bPOC" "DOC" "PAR" "t"# + +# ## Simulation +# Next we setup the simulation along with some callbacks that: +# - Couples the particles to the biodeochemical model +# - Update the PAR field from the surface PAR and phytoplankton concentration +# - Show the progress of the simulation +# - Store the model and particles output +# - Prevent the tracers from going negative from numerical error (see discussion of this in the [positivity preservation](@ref pos-preservation) implimentation page) + +simulation = Simulation(model, Δt=3.5minutes, stop_time=duration) + +# simulation.callbacks[:couple_particles] = Callback(Particles.infinitesimal_particle_field_coupling!; callsite = TendencyCallsite()) + +progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) +simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100)) + +filename = "kelp450_11_new" +simulation.output_writers[:profiles] = JLD2OutputWriter(model, merge(model.tracers, model.auxiliary_fields), filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing=true) +simulation.output_writers[:particles] = JLD2OutputWriter(model, (; kelp_particles), filename = "$(filename)_particles.jld2", schedule = TimeInterval(1day), overwrite_existing = true) + + +#simulation.callbacks[:timestep] = Callback(update_timestep!, IterationInterval(1), (c_forcing=0.1, c_adv=0.5, c_diff=0.5, w = 200/day, relaxation=0.95), TimeStepCallsite()) + +scale_negative_tracers = ScaleNegativeTracers(;model, tracers = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON)) #changed +simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) + +# plankton_redfield = model.biogeochemistry.phytoplankton_redfield +# scale_negative_carbon_tracers = ScaleNegativeTracers(tracers = (:P, :Z, :DOC, :sPOC, :bPOC, :DIC), scalefactors = (P = plankton_redfield, Z = plankton_redfield, DOC = 1, sPOC = 1, bPOC = 1, DIC = 1)) +# simulation.callbacks[:neg2] = Callback(scale_negative_carbon_tracers; callsite = UpdateStateCallsite()) + + +#wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 2.0, min_change = 0.5, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) +#simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) + +#simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=SpecifiedTimes([i*years for i=1:12]), prefix=filename) #prefix="kelp_01" +# ## Run! +# Finally we run the simulation +run!(simulation) + + + + + + + + + + + + + + + + + + diff --git a/subpolarnew/kelp450_11_new/kelp_old_example.jl b/subpolarnew/kelp450_11_new/kelp_old_example.jl new file mode 100644 index 000000000..e8e3eea2b --- /dev/null +++ b/subpolarnew/kelp450_11_new/kelp_old_example.jl @@ -0,0 +1,206 @@ +# # Simple active particle example +# In this example we will setup a simple 1D column with the [LOBSTER](@ref LOBSTER) biogeochemical model and active particles modelling the growth of sugar kelp. This demonstrates: +# - How to setup OceanBioME's biogeochemical models +# - How to add biologically active particles which interact with the biodeochemical model +# - How to visualise results + +# This is forced by idealised mixing layer depth and surface photosynthetically available radiation (PAR) which are setup first + +# ## Install dependencies +# First we will check we have the dependencies installed +# ```julia +# using Pkg +# pkg "add OceanBioME, Oceananigans, CairoMakie, JLD2" +# ``` + +# ## Model setup +# First load the required packages +using OceanBioME, Oceananigans, Printf +using Oceananigans.Units +using Oceananigans.Architectures: arch_array + +const year = years = 365days # just for these idealised cases +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 H(t, t₀, t₁) = ifelse(t₀ < t < t₁, 1.0, 0.0) + +@inline fmld1(t) = H(t, 50days, year) * (1 / (1 + exp(-(t - 100days) / 5days))) * (1 / (1 + exp((t - 330days) / 25days))) + +@inline MLD(t) = - (10 + 340 * (1 - fmld1(year - eps(year)) * exp(-mod(t, year) / 25days) - fmld1(mod(t, year)))) + +@inline κₜ(x, y, z, t) = 1e-2 * (1 + tanh((z - MLD(t))/10)) / 2 + 1e-4 + +@inline temp(x, y, z, t) = 2.4 * cos(t * 2π / year + 50day) + 10 +nothing #hide + +architecture = CPU() + +# ## Grid and PAR field +# Define the grid and an extra Oceananigans field for the PAR to be stored in +Lx, Ly = 20meters, 20meters +grid = RectilinearGrid(architecture, size=(1, 1, 50), extent=(Lx, Ly, 200)) + +# Specify the boundary conditions for DIC and O₂ based on the air-sea CO₂ and O₂ flux +CO₂_flux = GasExchange(; gas = :CO₂, temperature = temp, salinity = (args...) -> 35) + +# ## Kelp Particle setup +@info "Setting up kelp particles" + +n = 5 # number of kelp fronds +z₀ = [-21:5:-1;] * 1.0 # depth of kelp fronds + +particles = SLatissima(; architecture, + x = arch_array(architecture, ones(n) * Lx / 2), + y = arch_array(architecture, ones(n) * Ly / 2), + z = arch_array(architecture, z₀), + A = arch_array(architecture, ones(n) * 10.0), + latitude = 57.5, + scalefactor = 500.0, + pescribed_temperature = temp) + +# ## Setup BGC model +biogeochemistry = LOBSTER(; grid, + surface_phytosynthetically_active_radiation = PAR⁰, + carbonates = true, + variable_redfield = true, + particles) + +model = NonhydrostaticModel(; grid, + closure = ScalarDiffusivity(ν = κₜ, κ = κₜ), + biogeochemistry) + +set!(model, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2239.8, Alk = 2409.0) + +# ## Simulation +# Next we setup the simulation along with some callbacks that: +# - Show the progress of the simulation +# - Store the model and particles output +# - Prevent the tracers from going negative from numerical error (see discussion of this in the [positivity preservation](@ref pos-preservation) implementation page) + +simulation = Simulation(model, Δt = 3minutes, stop_time = 100days) + +progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) + +simulation.callbacks[:progress] = Callback(progress_message, TimeInterval(10days)) + +filename = "kelp" +simulation.output_writers[:profiles] = JLD2OutputWriter(model, merge(model.tracers, model.auxiliary_fields), + filename = "$filename.jld2", + schedule = TimeInterval(1day), + overwrite_existing = true) + +simulation.output_writers[:particles] = JLD2OutputWriter(model, (; particles), + filename = "$(filename)_particles.jld2", + schedule = TimeInterval(1day), + overwrite_existing = true) + +scale_negative_tracers = ScaleNegativeTracers(; model, tracers = (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON)) +simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) +nothing #hide + +# ## Run! +# Finally we run the simulation +run!(simulation) + +# Now we can visualise the results with some post processing to diagnose the air-sea CO₂ flux - hopefully this looks different to the example without kelp! + + P = FieldTimeSeries("$filename.jld2", "P") + NO₃ = FieldTimeSeries("$filename.jld2", "NO₃") + Z = FieldTimeSeries("$filename.jld2", "Z") +sPOC = FieldTimeSeries("$filename.jld2", "sPOC") +bPOC = FieldTimeSeries("$filename.jld2", "bPOC") + DIC = FieldTimeSeries("$filename.jld2", "DIC") + Alk = FieldTimeSeries("$filename.jld2", "Alk") + +x, y, z = nodes(P) +times = P.times + +# We compute the air-sea CO₂ flux at the surface (corresponding to vertical index `k = grid.Nz`) and +# the carbon export by computing how much carbon sinks below some arbirtrary depth; here we use depth +# that corresponds to `k = grid.Nz - 20`. +air_sea_CO₂_flux = zeros(length(times)) +carbon_export = zeros(length(times)) + +for (i, t) in enumerate(times) + air_sea_CO₂_flux[i] = CO₂_flux.condition.parameters(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] = sPOC[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.sPOM.w[1, 1, grid.Nz-20] + + bPOC[1, 1, grid.Nz-20, i] * model.biogeochemistry.sinking_velocities.bPOM.w[1, 1, grid.Nz-20] +end + +# Both `air_sea_CO₂_flux` and `carbon_export` are in units `mmol CO₂ / (m² s)`. + +using CairoMakie + +fig = Figure(resolution = (1000, 1500), fontsize = 20) + +axis_kwargs = (xlabel = "Time (days)", ylabel = "z (m)", limits = ((0, times[end] / days), (-85meters, 0))) + +axP = Axis(fig[1, 1]; title = "Phytoplankton concentration (mmol N/m³)", axis_kwargs...) +hmP = heatmap!(times / days, z, interior(P, 1, 1, :, :)', colormap = :batlow) +Colorbar(fig[1, 2], hmP) + +axNO₃ = Axis(fig[2, 1]; title = "Nitrate concentration (mmol N/m³)", axis_kwargs...) +hmNO₃ = heatmap!(times / days, z, interior(NO₃, 1, 1, :, :)', colormap = :batlow) +Colorbar(fig[2, 2], hmNO₃) + +axZ = Axis(fig[3, 1]; title = "Zooplankton concentration (mmol N/m³)", axis_kwargs...) +hmZ = heatmap!(times / days, z, interior(Z, 1, 1, :, :)', colormap = :batlow) +Colorbar(fig[3, 2], hmZ) + +axD = Axis(fig[4, 1]; title = "Detritus concentration (mmol C/m³)", axis_kwargs...) +hmD = heatmap!(times / days, z, interior(sPOC, 1, 1, :, :)' .+ interior(bPOC, 1, 1, :, :)', colormap = :batlow) +Colorbar(fig[4, 2], hmD) + +CO₂_molar_mass = (12 + 2 * 16) * 1e-3 # kg / mol + +axfDIC = Axis(fig[5, 1], xlabel = "Time (days)", ylabel = "Flux (kgCO₂/m²/year)", + title = "Air-sea CO₂ flux and Sinking", limits = ((0, times[end] / days), nothing)) +lines!(axfDIC, times / days, air_sea_CO₂_flux / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Air-sea flux") +lines!(axfDIC, times / days, carbon_export / 1e3 * CO₂_molar_mass * year, linewidth = 3, label = "Sinking export") +Legend(fig[5, 2], axfDIC, framevisible = false) + +fig + +# We can also have a look at how the kelp particles evolve +using JLD2 + +file = jldopen("$(filename)_particles.jld2") + +iterations = keys(file["timeseries/t"]) + +A, N, C = ntuple(n -> zeros(5, length(iterations)), 3) + +times = zeros(length(iterations)) + +for (i, iter) in enumerate(iterations) + particles = file["timeseries/particles/$iter"] + A[:, i] = particles.A + N[:, i] = particles.N + C[:, i] = particles.C + + times[i] = file["timeseries/t/$iter"] +end + +fig = Figure(resolution = (1000, 800), fontsize = 20) + +axis_kwargs = (xlabel = "Time (days)", limits = ((0, times[end] / days), nothing)) + +ax1 = Axis(fig[1, 1]; ylabel = "Frond area (dm²)", axis_kwargs...) +[lines!(ax1, times / day, A[n, :], linewidth = 3) for n in 1:5] + +ax2 = Axis(fig[2, 1]; ylabel = "Total Carbon (gC)", axis_kwargs...) +[lines!(ax2, times / day, (@. A * (N + particles.structural_nitrogen) * particles.structural_dry_weight_per_area)[n, :], linewidth = 3) for n in 1:5] + +ax3 = Axis(fig[3, 1]; ylabel = "Total Nitrogen (gN)", axis_kwargs...) +[lines!(ax3, times / day, (@. A * (C + particles.structural_carbon) * particles.structural_dry_weight_per_area)[n, :], linewidth = 3) for n in 1:5] + +fig diff --git a/test_stretch/readme.rtf b/test_stretch/readme.rtf new file mode 100644 index 000000000..d990989d4 --- /dev/null +++ b/test_stretch/readme.rtf @@ -0,0 +1,17 @@ +{\rtf1\ansi\ansicpg1252\cocoartf2639 +\cocoatextscaling0\cocoaplatform0{\fonttbl\f0\fnil\fcharset134 PingFangSC-Regular;} +{\colortbl;\red255\green255\blue255;} +{\*\expandedcolortbl;;} +\paperw11900\paperh16840\margl1440\margr1440\vieww11520\viewh8400\viewkind0 +\pard\tx720\tx1440\tx2160\tx2880\tx3600\tx4320\tx5040\tx5760\tx6480\tx7200\tx7920\tx8640\pardirnatural\partightenfactor0 + +\f0\fs24 \cf0 1. modify original model in the new OceanBioME, including withoutkelp and with kelp of varying densities, particularly high kelp density. Check whether the results are the same with my old model. Don\'a1\'aft forget to add Carbon limitation modifier. \ +\ +See how the new model deal with negative tracers.\ +\ +\ +2. See if plotting codes need to be modified. \ +\ +3. See if new OceanBioME can run stretched grid.\ +\ +4. Carbon limitation} \ No newline at end of file From 385215a20f1d7b618631293f16998c2c516fd5bd Mon Sep 17 00:00:00 2001 From: sichen Date: Fri, 19 Jan 2024 09:34:14 +0000 Subject: [PATCH 2/5] add package Interpolations --- Manifest.toml | 44 ++++++++++++++++++++++++++++++++++++++++++-- Project.toml | 1 + 2 files changed, 43 insertions(+), 2 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 55fd8bd81..0a740c44c 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.9.3" +julia_version = "1.9.2" manifest_format = "2.0" -project_hash = "9c3a43fcc42b0e8b4f92f9c38f0c935f93e8165c" +project_hash = "557a107566f87409cf6223382af3382ea9ba1b9e" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] @@ -60,6 +60,12 @@ git-tree-sha1 = "c06a868224ecba914baa6942988e2f2aade419be" uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458" version = "0.1.0" +[[deps.AxisAlgorithms]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] +git-tree-sha1 = "01b8ccb13d68535d73d2b0c23e39bd23155fb712" +uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950" +version = "1.1.0" + [[deps.BFloat16s]] deps = ["LinearAlgebra", "Printf", "Random", "Test"] git-tree-sha1 = "dbf84058d0a8cbbadee18d25cf606934b22d7c66" @@ -298,6 +304,18 @@ version = "2023.2.0+0" deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +[[deps.Interpolations]] +deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] +git-tree-sha1 = "88a101217d7cb38a7b481ccd50d21876e1d1b0e0" +uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +version = "0.15.1" + + [deps.Interpolations.extensions] + InterpolationsUnitfulExt = "Unitful" + + [deps.Interpolations.weakdeps] + Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + [[deps.IrrationalConstants]] git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" @@ -641,6 +659,18 @@ git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111" uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" version = "1.5.3" +[[deps.Ratios]] +deps = ["Requires"] +git-tree-sha1 = "1342a47bf3260ee108163042310d26f2be5ec90b" +uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439" +version = "0.4.5" + + [deps.Ratios.extensions] + RatiosFixedPointNumbersExt = "FixedPointNumbers" + + [deps.Ratios.weakdeps] + FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" + [[deps.RealDot]] deps = ["LinearAlgebra"] git-tree-sha1 = "9f0a1b71baaf7650f4fa8a1d168c7fb6ee41f0c9" @@ -712,6 +742,10 @@ git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" version = "1.1.1" +[[deps.SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" @@ -892,6 +926,12 @@ git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868" uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" version = "1.3.0" +[[deps.WoodburyMatrices]] +deps = ["LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "c1a7aa6219628fcd757dede0ca95e245c5cd9511" +uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" +version = "1.0.0" + [[deps.XML2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] git-tree-sha1 = "24b81b59bd35b3c42ab84fa589086e19be919916" diff --git a/Project.toml b/Project.toml index 1d5a043a6..c1c5c1f49 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.9.2" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" From 6820e9fe52fc5c405d71da19e7fc23bb66ac53fe Mon Sep 17 00:00:00 2001 From: sichen Date: Mon, 22 Jan 2024 11:36:18 +0000 Subject: [PATCH 3/5] add one update code --- subpolarnew/kelp450_11_new/kelp450_11_new3.jl | 198 ++++++++++++++++++ 1 file changed, 198 insertions(+) create mode 100644 subpolarnew/kelp450_11_new/kelp450_11_new3.jl diff --git a/subpolarnew/kelp450_11_new/kelp450_11_new3.jl b/subpolarnew/kelp450_11_new/kelp450_11_new3.jl new file mode 100644 index 000000000..5aa8653e1 --- /dev/null +++ b/subpolarnew/kelp450_11_new/kelp450_11_new3.jl @@ -0,0 +1,198 @@ +# # Simple active particle example +# In this example we will setup a simple 1D column with the [LOBSTER](@ref LOBSTER) biogeochemical model and active particles modelling the growth of sugar kelp. This demonstraits: +# - How to setup OceanBioME's biogeochemical models +# - How to setup light attenuation +# - How to add biologically active particles which interact with the biodeochemical model +# - How to include optional tracer sets (carbonate chemistry and oxygen) +# - How to visulise results + +# This is forced by idealised mixing layer depth and surface photosynthetically available radiation (PAR) which are setup first + +# ## Install dependencies +# First we will check we have the dependencies installed +# ```julia +# using Pkg +# pkg"add OceanBioME, Oceananigans, Printf, GLMakie" +# ``` + +# ## Model setup +# We load the packages and choose the default LOBSTER parameter set +using OceanBioME, Oceananigans,Printf +using Oceananigans.Fields: FunctionField, ConstantField +# using Oceananigans.Units: second, minute, minutes, hour, hours, day, days#, year, years +using Oceananigans.Units #####changed +using Interpolations +using JLD2 +using Oceananigans.Architectures: arch_array #####changed +const year = years = 365days #####changed + +# ## 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) + +##PAR⁰(x, y, t) = 60*(1-cos((t+15days)*2π/(365days)))*(1 /(1 +0.2*exp(-((mod(t, 365days)-200days)/50days)^2))) .+ 2 +##H(t, t₀, t₁) = ifelse(t₀ keys(pp["timeseries"]) 15-element Vector{String}: "NO₃" "NH₄" "P" "Z" "sPON" "bPON" "DON" "DIC" "Alk" "O₂" "sPOC" "bPOC" "DOC" "PAR" "t"# + +# ## Simulation +# Next we setup the simulation along with some callbacks that: +# - Couples the particles to the biodeochemical model +# - Update the PAR field from the surface PAR and phytoplankton concentration +# - Show the progress of the simulation +# - Store the model and particles output +# - Prevent the tracers from going negative from numerical error (see discussion of this in the [positivity preservation](@ref pos-preservation) implimentation page) + +simulation = Simulation(model, Δt=3.5minutes, stop_time=duration) + +# simulation.callbacks[:couple_particles] = Callback(Particles.infinitesimal_particle_field_coupling!; callsite = TendencyCallsite()) + +progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) +simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100)) + +filename = "kelp450_11_new2" +simulation.output_writers[:profiles] = JLD2OutputWriter(model, model.tracers, filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing=true) #merge(model.tracers, model.auxiliary_fields), +simulation.output_writers[:particles] = JLD2OutputWriter(model, (; kelp_particles), filename = "$(filename)_particles.jld2", schedule = TimeInterval(1day), overwrite_existing = true) + + +#simulation.callbacks[:timestep] = Callback(update_timestep!, IterationInterval(1), (c_forcing=0.1, c_adv=0.5, c_diff=0.5, w = 200/day, relaxation=0.95), TimeStepCallsite()) + +scale_negative_tracers = ScaleNegativeTracers((:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) #ScaleNegativeTracers(;model, tracers = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) changed (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) +simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) + +# plankton_redfield = model.biogeochemistry.phytoplankton_redfield +# scale_negative_carbon_tracers = ScaleNegativeTracers(tracers = (:P, :Z, :DOC, :sPOC, :bPOC, :DIC), scalefactors = (P = plankton_redfield, Z = plankton_redfield, DOC = 1, sPOC = 1, bPOC = 1, DIC = 1)) +# simulation.callbacks[:neg2] = Callback(scale_negative_carbon_tracers; callsite = UpdateStateCallsite()) + + +#wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 2.0, min_change = 0.5, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) +#simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) + +#simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=SpecifiedTimes([i*years for i=1:12]), prefix=filename) #prefix="kelp_01" +# ## Run! +# Finally we run the simulation +run!(simulation) + + + + + + + + + + + + + + + + + + From 0a897a6c9c4d8969fdc313b36fe8b811364ccdd7 Mon Sep 17 00:00:00 2001 From: sichen Date: Thu, 25 Jan 2024 15:45:49 +0000 Subject: [PATCH 4/5] carbon limitation in the new version of OceanBioME 0.9.2 --- src/Models/Individuals/SLatissima.jl | 24 +- src/Models/Individuals/SLatissima_original.jl | 546 ++++++++++++++++++ 2 files changed, 566 insertions(+), 4 deletions(-) create mode 100644 src/Models/Individuals/SLatissima_original.jl diff --git a/src/Models/Individuals/SLatissima.jl b/src/Models/Individuals/SLatissima.jl index 66594deb5..1b38438f1 100644 --- a/src/Models/Individuals/SLatissima.jl +++ b/src/Models/Individuals/SLatissima.jl @@ -349,6 +349,18 @@ function update_lagrangian_particle_properties!(particles::SLatissima, model, bg particles.custom_dynamics(particles, model, bgc, Δt) end +function C_uptake_modifier(dic) ###change + if dic >= 2200 + return 1.0 + elseif dic >= 1650 && dic <2200 + return 0.5 * dic /550 -1 + elseif dic >=1282 && dic <1650 + return 1/736*(dic -1282) + else + return 0 + end +end + @kernel function _update_lagrangian_particle_properties!(p, light_attenuation, bgc, grid, velocities, tracers, clock, Δt) idx = @index(Global) @@ -365,9 +377,9 @@ end t = clock.time @inbounds if p.A[idx] > 0 - NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(light_attenuation).PAR) + NO₃, NH₄, PAR, u, T, S, DIC = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(light_attenuation).PAR) ###change - photo = photosynthesis(T, PAR, p) + photo = photosynthesis(T, PAR, p) * C_uptake_modifier(DIC) ###change e = exudation(C, p) ν = erosion(A, p) @@ -525,7 +537,11 @@ end else NH₄ = 0.0 end - + if :DIC in bgc_tracers ###change + DIC = _interpolate(tracers.DIC, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) + else + DIC = 2200.0 + end if isnothing(particles.prescribed_velocity) u = sqrt(_interpolate(velocities.u, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) .^ 2 + _interpolate(velocities.v, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) .^ 2 + @@ -540,7 +556,7 @@ end S = _interpolate(tracers.S, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) - return NO₃, NH₄, PAR, u, T, S + return NO₃, NH₄, PAR, u, T, S, DIC ###change NO₃, NH₄, PAR, u, T, S end end #module diff --git a/src/Models/Individuals/SLatissima_original.jl b/src/Models/Individuals/SLatissima_original.jl new file mode 100644 index 000000000..66594deb5 --- /dev/null +++ b/src/Models/Individuals/SLatissima_original.jl @@ -0,0 +1,546 @@ +""" +Sugar kelp model of [Broch2012](@citet) and updated by [Broch2013](@citet), [Fossberg2018](@citet), and [Broch2019](@citet). + +Prognostic properties +===================== +* Area: A (dm²) +* Nitrogen reserve: N (gN/gSW) +* Carbon reserve: C (gC/gSW) + +Tracer dependencies +=================== +* Nitrates: NO₃ (mmol N/m³) +* Photosynthetically available radiation: PAR (einstein/m²/day) + +Optional +======== +* Ammonia: NH₄ (mmol N/m³) +""" +module SLatissimaModel + +using Roots, KernelAbstractions +using OceanBioME.Particles: BiogeochemicalParticles, get_node +using Oceananigans.Units +using Oceananigans: Center, CPU +using Oceananigans.Architectures: arch_array, device, architecture +using Oceananigans.Biogeochemistry: required_biogeochemical_tracers, biogeochemical_auxiliary_fields + +using KernelAbstractions.Extras.LoopInfo: @unroll +using Oceananigans.Operators: volume +using Oceananigans.Grids: AbstractGrid +using Oceananigans.Fields: fractional_indices, _interpolate, datatuple +using Oceananigans.Models: total_velocities + +import Adapt: adapt_structure, adapt +import Oceananigans.Biogeochemistry: update_tendencies! +import Oceananigans.Models.LagrangianParticleTracking: update_lagrangian_particle_properties!, _advect_particles! + +@inline no_dynamics(args...) = nothing + + +""" + SLatissima(; architecture :: AR = CPU(), + growth_rate_adjustment :: FT = 4.5, + photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60), + minimum_carbon_reserve :: FT = 0.01, + structural_carbon :: FT = 0.2, + exudation :: FT = 0.5, + erosion :: FT = 0.22, + saturation_irradiance :: FT = 90 * day/ (10 ^ 6), + structural_dry_weight_per_area :: FT = 0.5, + structural_dry_to_wet_weight :: FT = 0.0785, + carbon_reserve_per_carbon :: FT = 2.1213, + nitrogen_reserve_per_nitrogen :: FT = 2.72, + minimum_nitrogen_reserve :: FT = 0.0126, + maximum_nitrogen_reserve :: FT = 0.0216, + growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)), + growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2, + maximum_specific_growth_rate :: FT = 0.18, + structural_nitrogen :: FT = 0.0146, + photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24, + photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24, + photosynthesis_ref_temp_1 :: FT = 285.0, + photosynthesis_ref_temp_2 :: FT = 288.0, + photoperiod_1 :: FT = 0.85, + photoperiod_2 :: FT = 0.3, + respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24, + respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24, + respiration_ref_temp_1 :: FT = 285.0, + respiration_ref_temp_2 :: FT = 290.0, + photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1), + photosynthesis_low_temp :: FT = 271.0, + photosynthesis_high_temp :: FT = 296.0, + photosynthesis_high_arrhenius_temp :: FT = 1414.87, + photosynthesis_low_arrhenius_temp :: FT = 4547.89, + respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1), + current_speed_for_0p65_uptake :: FT = 0.03, + nitrate_half_saturation :: FT = 4.0, + ammonia_half_saturation :: FT = 1.3, + maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6), + maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6), + current_1 :: FT = 0.72, + current_2 :: FT = 0.28, + current_3 :: FT = 0.045, + respiration_reference_A :: FT = 1.11e-4 * 24, + respiration_reference_B :: FT = 5.57e-5 * 24, + exudation_redfield_ratio :: FT = Inf, + + prescribed_velocity :: U = 0.1, + + #position + x :: P = arch_array(architecture, [0.0]) + y :: P = arch_array(architecture, zeros(Float64, length(x))), + z :: P = arch_array(architecture, zeros(Float64, length(x))), + + #properties + A :: P = arch_array(architecture, ones(Float64, length(x)) * 30), + N :: P = arch_array(architecture, ones(Float64, length(x)) * 0.01), + C :: P = arch_array(architecture, ones(Float64, length(x)) * 0.1), + + #feedback + nitrate_uptake :: P = arch_array(architecture, zeros(Float64, length(x))), + ammonia_uptake :: P = arch_array(architecture, zeros(Float64, length(x))), + primary_production :: P = arch_array(architecture, zeros(Float64, length(x))), + frond_exudation :: P = arch_array(architecture, zeros(Float64, length(x))), + nitrogen_erosion :: P = arch_array(architecture, zeros(Float64, length(x))), + carbon_erosion :: P = arch_array(architecture, zeros(Float64, length(x))), + + custom_dynamics :: F = no_dynamics, + + scalefactor :: FT = 1.0, + latitude :: FT = 57.5) + +Keyword Arguments +================= + +- `architecture`: the architecture to adapt arrays to +- `growth_rate_adjustment`, ..., `exudation_redfield_ratio`: parameter values +- `prescribed_velocity`: functions for the relative velocity `f(x, y, z, t)` +- `x`,`y` and `z`: positions of the particles +- `A`, `N`, and `C`: area, nitrogen, and carbon reserves +- `nitrate_uptake` ... `carbon_erosion`: diagnostic values coupled to tracer fields +- `custom_dynamics`: place to add any function of form `f!(particles, model, bgc, Δt)` +- `scalefactor`: scalar scaling for tracer coupling +- `latitude`: model latitude for seasonal growth modulation +""" +Base.@kwdef struct SLatissima{AR, FT, U, P, F} <: BiogeochemicalParticles + architecture :: AR = CPU() + growth_rate_adjustment :: FT = 4.5 + photosynthetic_efficiency :: FT = 4.15e-5 * 24 * 10^6 / (24 * 60 * 60) + minimum_carbon_reserve :: FT = 0.01 + structural_carbon :: FT = 0.2 + exudation :: FT = 0.5 + erosion :: FT = 0.22 + saturation_irradiance :: FT = 90 * day/ (10 ^ 6) + structural_dry_weight_per_area :: FT = 0.5 + structural_dry_to_wet_weight :: FT = 0.0785 + carbon_reserve_per_carbon :: FT = 2.1213 + nitrogen_reserve_per_nitrogen :: FT = 2.72 + minimum_nitrogen_reserve :: FT = 0.0126 + maximum_nitrogen_reserve :: FT = 0.0216 + growth_adjustment_2 :: FT = 0.039 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) + growth_adjustment_1 :: FT = 0.18 / (2 * (1 - minimum_nitrogen_reserve / maximum_nitrogen_reserve)) - growth_adjustment_2 + maximum_specific_growth_rate :: FT = 0.18 + structural_nitrogen :: FT = 0.0146 + photosynthesis_at_ref_temp_1 :: FT = 1.22e-3 * 24 + photosynthesis_at_ref_temp_2 :: FT = 1.3e-3 * 24 + photosynthesis_ref_temp_1 :: FT = 285.0 + photosynthesis_ref_temp_2 :: FT = 288.0 + photoperiod_1 :: FT = 0.85 + photoperiod_2 :: FT = 0.3 + respiration_at_ref_temp_1 :: FT = 2.785e-4 * 24 + respiration_at_ref_temp_2 :: FT = 5.429e-4 * 24 + respiration_ref_temp_1 :: FT = 285.0 + respiration_ref_temp_2 :: FT = 290.0 + photosynthesis_arrhenius_temp :: FT = (1 / photosynthesis_ref_temp_1 - 1 / photosynthesis_ref_temp_2) ^ -1 * log(photosynthesis_at_ref_temp_2 / photosynthesis_at_ref_temp_1) + photosynthesis_low_temp :: FT = 271.0 + photosynthesis_high_temp :: FT = 296.0 + photosynthesis_high_arrhenius_temp :: FT = 1414.87 + photosynthesis_low_arrhenius_temp :: FT = 4547.89 + respiration_arrhenius_temp :: FT = (1 / respiration_ref_temp_1 - 1 / respiration_ref_temp_2) ^ -1 * log(respiration_at_ref_temp_2 / respiration_at_ref_temp_1) + current_speed_for_0p65_uptake :: FT = 0.03 + nitrate_half_saturation :: FT = 4.0 + ammonia_half_saturation :: FT = 1.3 + maximum_nitrate_uptake :: FT = 10 * structural_dry_weight_per_area * 24 * 14 / (10^6) + maximum_ammonia_uptake :: FT = 12 * structural_dry_weight_per_area * 24 * 14 / (10^6) + current_1 :: FT = 0.72 + current_2 :: FT = 0.28 + current_3 :: FT = 0.045 + respiration_reference_A :: FT = 1.11e-4 * 24 + respiration_reference_B :: FT = 5.57e-5 * 24 + exudation_redfield_ratio :: FT = Inf + + prescribed_velocity :: U = 0.1 + + #position + x :: P = arch_array(architecture, [0.0]) + y :: P = arch_array(architecture, zeros(Float64, length(x))) + z :: P = arch_array(architecture, zeros(Float64, length(x))) + + #properties + A :: P = arch_array(architecture, ones(Float64, length(x)) * 30) + N :: P = arch_array(architecture, ones(Float64, length(x)) * 0.01) + C :: P = arch_array(architecture, ones(Float64, length(x)) * 0.1) + + #feedback + nitrate_uptake :: P = arch_array(architecture, zeros(Float64, length(x))) + ammonia_uptake :: P = arch_array(architecture, zeros(Float64, length(x))) + primary_production :: P = arch_array(architecture, zeros(Float64, length(x))) + frond_exudation :: P = arch_array(architecture, zeros(Float64, length(x))) + nitrogen_erosion :: P = arch_array(architecture, zeros(Float64, length(x))) + carbon_erosion :: P = arch_array(architecture, zeros(Float64, length(x))) + + custom_dynamics :: F = no_dynamics + + scalefactor :: FT = 1.0 + latitude :: FT = 57.5 +end + +adapt_structure(to, kelp::SLatissima) = SLatissima(adapt(to, kelp.architecture), + adapt(to, kelp.growth_rate_adjustment), + adapt(to, kelp.photosynthetic_efficiency), + adapt(to, kelp.minimum_carbon_reserve), + adapt(to, kelp.structural_carbon), + adapt(to, kelp.exudation), + adapt(to, kelp.erosion), + adapt(to, kelp.saturation_irradiance), + adapt(to, kelp.structural_dry_weight_per_area), + adapt(to, kelp.structural_dry_to_wet_weight), + adapt(to, kelp.carbon_reserve_per_carbon), + adapt(to, kelp.nitrogen_reserve_per_nitrogen), + adapt(to, kelp.minimum_nitrogen_reserve), + adapt(to, kelp.maximum_nitrogen_reserve), + adapt(to, kelp.growth_adjustment_2), + adapt(to, kelp.growth_adjustment_1), + adapt(to, kelp.maximum_specific_growth_rate), + adapt(to, kelp.structural_nitrogen), + adapt(to, kelp.photosynthesis_at_ref_temp_1), + adapt(to, kelp.photosynthesis_at_ref_temp_2), + adapt(to, kelp.photosynthesis_ref_temp_1), + adapt(to, kelp.photosynthesis_ref_temp_2), + adapt(to, kelp.photoperiod_1), + adapt(to, kelp.photoperiod_2), + adapt(to, kelp.respiration_at_ref_temp_1), + adapt(to, kelp.respiration_at_ref_temp_2), + adapt(to, kelp.respiration_ref_temp_1), + adapt(to, kelp.respiration_ref_temp_2), + adapt(to, kelp.photosynthesis_arrhenius_temp), + adapt(to, kelp.photosynthesis_low_temp), + adapt(to, kelp.photosynthesis_high_temp), + adapt(to, kelp.photosynthesis_high_arrhenius_temp), + adapt(to, kelp.photosynthesis_low_arrhenius_temp), + adapt(to, kelp.respiration_arrhenius_temp), + adapt(to, kelp.current_speed_for_0p65_uptake), + adapt(to, kelp.nitrate_half_saturation), + adapt(to, kelp.ammonia_half_saturation), + adapt(to, kelp.maximum_nitrate_uptake), + adapt(to, kelp.maximum_ammonia_uptake), + adapt(to, kelp.current_1), + adapt(to, kelp.current_2), + adapt(to, kelp.current_3), + adapt(to, kelp.respiration_reference_A), + adapt(to, kelp.respiration_reference_B), + adapt(to, kelp.exudation_redfield_ratio), + adapt(to, kelp.prescribed_velocity), + adapt(to, kelp.x), + adapt(to, kelp.y), + adapt(to, kelp.z), + adapt(to, kelp.A), + adapt(to, kelp.N), + adapt(to, kelp.C), + adapt(to, kelp.nitrate_uptake), + adapt(to, kelp.ammonia_uptake), + adapt(to, kelp.primary_production), + adapt(to, kelp.frond_exudation), + adapt(to, kelp.nitrogen_erosion), + adapt(to, kelp.carbon_erosion), + adapt(to, kelp.custom_dynamics), + adapt(to, kelp.scalefactor), + adapt(to, kelp.latitude)) + +function update_tendencies!(bgc, particles::SLatissima, model) + num_particles = length(particles) + workgroup = min(num_particles, 256) + worksize = num_particles + + update_tracer_tendencies_kernel! = update_tracer_tendencies!(device(model.architecture), workgroup, worksize) + update_tracer_tendencies_kernel!(bgc, particles, model.timestepper.Gⁿ, model.grid) + + return nothing +end + +@kernel function update_tracer_tendencies!(bgc, p, tendencies, grid::AbstractGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} + idx = @index(Global) + + x = p.x[idx] + y = p.y[idx] + z = p.z[idx] + + bgc_tracers = required_biogeochemical_tracers(bgc) + + i, j, k = fractional_indices(x, y, z, (Center(), Center(), Center()), grid) + + # Convert fractional indices to unit cell coordinates 0 ≤ (ξ, η, ζ) ≤ 1 + # and integer indices (with 0-based indexing). + ξ, i = modf(i) + η, j = modf(j) + ζ, k = modf(k) + + i, j, k = (get_node(TX(), Int(ifelse(ξ < 0.5, i + 1, i + 2)), grid.Nx), + get_node(TY(), Int(ifelse(η < 0.5, j + 1, j + 2)), grid.Ny), + get_node(TZ(), Int(ifelse(ζ < 0.5, k + 1, k + 2)), grid.Nz)) + + node_volume = volume(i, j, k, grid, Center(), Center(), Center()) + + node_scalefactor = p.scalefactor / node_volume #* normfactor / (weight * node_volume) + + @inbounds begin + tendencies.NO₃[i, j, k] -= node_scalefactor * p.nitrate_uptake[idx] + + if :NH₄ in bgc_tracers + tendencies.NH₄[i, j, k] -= node_scalefactor * p.ammonia_uptake[idx] + end + + if :DIC in bgc_tracers + tendencies.DIC[i, j, k] -= node_scalefactor * p.primary_production[idx] + end + + if :O₂ in bgc_tracers + tendencies.O₂[i, j, k] += node_scalefactor * p.primary_production[idx] + end + + if :DOM in bgc_tracers + tendencies.DOM[i, j, k] += node_scalefactor * p.frond_exudation[idx] / p.exudation_redfield_ratio + elseif :DON in bgc_tracers + tendencies.DON[i, j, k] += node_scalefactor * p.frond_exudation[idx] / p.exudation_redfield_ratio + tendencies.DOC[i, j, k] += node_scalefactor * p.frond_exudation[idx] + end + + if :bPOM in bgc_tracers + tendencies.bPOM[i, j, k] += node_scalefactor * p.nitrogen_erosion[idx] + elseif :bPON in bgc_tracers + tendencies.bPON[i, j, k] += node_scalefactor * p.nitrogen_erosion[idx] + tendencies.bPOC[i, j, k] += node_scalefactor * p.carbon_erosion[idx] + end + end +end + +function update_lagrangian_particle_properties!(particles::SLatissima, model, bgc, Δt) + workgroup = min(length(particles), 256) + worksize = length(particles) + + arch = architecture(model) + + # Advect particles + advect_particles_kernel! = _advect_particles!(device(arch), workgroup, worksize) + + advect_particles_kernel!((x = particles.x, y = particles.y, z = particles.z), + 1.0, model.grid, Δt, + total_velocities(model)) + + + update_particle_properties_kernel! = _update_lagrangian_particle_properties!(device(arch), workgroup, worksize) + + tracer_fields = merge(model.tracers, model.auxiliary_fields) + + update_particle_properties_kernel!(particles, bgc.light_attenuation, bgc.underlying_biogeochemistry, model.grid, + total_velocities(model), tracer_fields, model.clock, Δt) + + particles.custom_dynamics(particles, model, bgc, Δt) +end + +@kernel function _update_lagrangian_particle_properties!(p, light_attenuation, bgc, grid, velocities, tracers, clock, Δt) + idx = @index(Global) + + @inbounds begin + x = p.x[idx] + y = p.y[idx] + z = p.z[idx] + + A = p.A[idx] + N = p.N[idx] + C = p.C[idx] + end + + t = clock.time + + @inbounds if p.A[idx] > 0 + NO₃, NH₄, PAR, u, T, S = get_arguments(x, y, z, t, p, bgc, grid, velocities, tracers, biogeochemical_auxiliary_fields(light_attenuation).PAR) + + photo = photosynthesis(T, PAR, p) + e = exudation(C, p) + ν = erosion(A, p) + + fᶜ = f_curr(u, p) + + j_NO₃ = max(0, p.maximum_ammonia_uptake * fᶜ * ((p.maximum_nitrogen_reserve - N) / + (p.maximum_nitrogen_reserve - p.minimum_nitrogen_reserve)) * NO₃ / + (p.nitrate_half_saturation + NO₃)) + + j̃_NH₄ = max(0, p.maximum_ammonia_uptake * fᶜ * NH₄ / (p.ammonia_half_saturation + NH₄)) + + μ_NH₄ = j̃_NH₄ / (p.structural_dry_weight_per_area * (N + p.structural_nitrogen)) + μ_NO₃ = 1 - p.minimum_nitrogen_reserve / N + μ_C = 1 - p.minimum_carbon_reserve / C + + n = floor(Int, mod(t, 364days) / day) + λ = normed_day_length_change(n, p.latitude) + + μ = @inbounds f_area(A, p) * + f_temp(T, p) * + f_photo(λ, p) * + min(μ_C, max(μ_NO₃, μ_NH₄)) + + j_NH₄ = min(j̃_NH₄, μ * p.structural_dry_weight_per_area * (N + p.structural_nitrogen)) + + r = respiration(T, μ, j_NO₃ + j_NH₄, p) + + dA = (μ - ν) * A / day + + dN = ((j_NO₃ + j_NH₄ - photo * e * 14 / (12 * p.exudation_redfield_ratio)) / p.structural_dry_weight_per_area - + μ * (N + p.structural_nitrogen)) / day + + dC = ((photo * (1 - e) - r) / p.structural_dry_weight_per_area - + μ * (C + p.structural_carbon)) / day + + A_new = A + dA * Δt + N_new = N + dN * Δt + C_new = C + dC * Δt + + if C_new < p.minimum_carbon_reserve + A_new *= (1 - (p.minimum_carbon_reserve - C) / p.structural_carbon) + C_new = p.minimum_carbon_reserve + N_new += p.structural_nitrogen * (p.minimum_carbon_reserve - C) / p.structural_carbon + end + + if N_new < p.minimum_nitrogen_reserve + A_new *= (1 - (p.minimum_nitrogen_reserve - N) / p.structural_nitrogen) + N_new = p.minimum_nitrogen_reserve + C_new += p.structural_carbon * (p.minimum_nitrogen_reserve - N) / p.structural_nitrogen + end + + @inbounds begin + p.primary_production[idx] = (photo - r) * A / (day * 12 * 0.001) #gC/dm^2/hr to mmol C/s + + p.frond_exudation[idx] = e * photo * A / (day * 12 * 0.001)#mmol C/s + + p.nitrogen_erosion[idx] = ν * p.structural_dry_weight_per_area * A * (p.structural_nitrogen + N) / (day * 14 * 0.001)#1/day to mmol N/s + + p.carbon_erosion[idx] = ν * p.structural_dry_weight_per_area * A * (p.structural_carbon + C) / (day * 12 * 0.001)#1/hr to mmol C/s + + p.nitrate_uptake[idx] = j_NO₃ * A / (day * 14 * 0.001)#gN/dm^2/hr to mmol N/s + + p.ammonia_uptake[idx] = j_NH₄ * A / (day * 14 * 0.001)#gN/dm^2/hr to mmol N/s + + p.A[idx] = A_new + p.N[idx] = N_new + p.C[idx] = C_new + end + end +end + +@inline function photosynthesis(T, PAR, p) + Tk = T + 273.15 + + pₘₐₓ = p.photosynthesis_at_ref_temp_1 * + exp(p.photosynthesis_arrhenius_temp / p.photosynthesis_ref_temp_1 - + p.photosynthesis_arrhenius_temp / Tk) / + (1 + exp(p.photosynthesis_low_arrhenius_temp / Tk - + p.photosynthesis_low_arrhenius_temp / p.photosynthesis_low_temp) + + exp(p.photosynthesis_high_arrhenius_temp / p.photosynthesis_high_temp - + p.photosynthesis_high_arrhenius_temp / Tk)) + + β = find_zero(β_func, (0, 0.1), Bisection(); p = (; pₘₐₓ, α = p.photosynthetic_efficiency, I_sat = p.saturation_irradiance)) + + pₛ = p.photosynthetic_efficiency * p.saturation_irradiance / log(1 + p.photosynthetic_efficiency / β) + + return pₛ * (1 - exp(- p.photosynthetic_efficiency * PAR / pₛ)) * exp(-β * PAR / pₛ) +end + +@inline β_func(x, p) = p.pₘₐₓ - (p.α * p.I_sat / log(1 + p.α / x)) * + (p.α / (p.α + x)) * (x / (p.α + x))^(x / p.α) + +@inline exudation(C, p) = 1 - exp(p.exudation * (p.minimum_carbon_reserve - C)) + +@inline erosion(A, p) = 1e-6 * exp(p.erosion * A) / (1 + 1e-6 * (exp(p.erosion * A) - 1)) + +@inline respiration(T, μ, j, p) = (p.respiration_reference_A * (μ / p.maximum_specific_growth_rate + + j / (p.maximum_nitrate_uptake + p.maximum_ammonia_uptake)) + + p.respiration_reference_B) * + exp(p.respiration_arrhenius_temp / p.respiration_ref_temp_1 - p.respiration_arrhenius_temp / (T + 273.15)) + +##### +##### Growth limitation +##### + +@inline f_curr(u, p) = p.current_1 * (1 - exp(-u / p.current_3)) + p.current_2 + +@inline f_area(a, p) = p.growth_adjustment_1 * exp(-(a / p.growth_rate_adjustment)^2) + p.growth_adjustment_2 + +@inline function f_temp(T, p) + # should probably parameterise these limits + if -1.8 <= T < 10 + return 0.08 * T + 0.2 + elseif 10 <= T <= 15 + return 1 + elseif 15 < T <= 19 + return 19 / 4 - T / 4 + elseif T > 19 + return 0 + else + return 0 + end +end + +@inline f_photo(λ, p) = p.photoperiod_1 * (1 + sign(λ) * abs(λ) ^ 0.5) + p.photoperiod_2 + +@inline function day_length(n, ϕ) + n -= 171 + M = mod((356.5291 + 0.98560028 * n), 360) + C = 1.9148 * sin(M * π / 180) + 0.02 * sin(2 * M * π / 180) + 0.0003 * sin(3 * M * π / 180) + λ = mod(M + C + 180 + 102.9372, 360) + δ = asin(sin(λ * π / 180) * sin(23.44 * π / 180)) + ω = (sin(-0.83 * π / 180) * sin(ϕ * π / 180) * sin(δ))/(cos(ϕ * π / 180) * cos(δ)) + return ω / 180 +end + +@inline normed_day_length_change(n, ϕ) = (day_length(n, ϕ) - day_length(n - 1, ϕ)) / (day_length(76, ϕ) - day_length(75, ϕ)) + +@inline function get_arguments(x, y, z, t, particles, bgc, grid, velocities, tracers, PAR_field) + bgc_tracers = required_biogeochemical_tracers(bgc) + + i, j, k = fractional_indices(x, y, z, (Center(), Center(), Center()), grid) + + ξ, i = modf(i) + η, j = modf(j) + ζ, k = modf(k) + + # TODO: ADD ALIASING/RENAMING OF TRACERS SO WE CAN USE E.G. N IN STEAD OF NO3 + + NO₃ = _interpolate(tracers.NO₃, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) + PAR = _interpolate(PAR_field, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) * day / (3.99e-10 * 545e12) # W / m² / s to einstein / m² / day + + if :NH₄ in bgc_tracers + NH₄ = _interpolate(tracers.NH₄, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) + else + NH₄ = 0.0 + end + + if isnothing(particles.prescribed_velocity) + u = sqrt(_interpolate(velocities.u, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) .^ 2 + + _interpolate(velocities.v, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) .^ 2 + + _interpolate(velocities.w, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) .^ 2) + elseif isa(particles.prescribed_velocity, Number) + u = particles.prescribed_velocity + else + u = particles.prescribed_velocity(x, y, z, t) + end + + T = _interpolate(tracers.T, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) + + S = _interpolate(tracers.S, ξ, η, ζ, Int(i+1), Int(j+1), Int(k+1)) + + return NO₃, NH₄, PAR, u, T, S +end + +end #module From be2c762e430f05bd8d4bc822f5ab0d985aebd772 Mon Sep 17 00:00:00 2001 From: sichen Date: Thu, 25 Jan 2024 15:46:57 +0000 Subject: [PATCH 5/5] add new kelp model scripts prescribed_velocity = 0.2 --- subpolarnew/kelp450_11_new/kelp450_11_new3.jl | 8 +- subpolarnew/kelp450_11_new/kelp450_11_new4.jl | 201 ++++++++++++++++++ 2 files changed, 206 insertions(+), 3 deletions(-) create mode 100644 subpolarnew/kelp450_11_new/kelp450_11_new4.jl diff --git a/subpolarnew/kelp450_11_new/kelp450_11_new3.jl b/subpolarnew/kelp450_11_new/kelp450_11_new3.jl index 5aa8653e1..fca2309f0 100644 --- a/subpolarnew/kelp450_11_new/kelp450_11_new3.jl +++ b/subpolarnew/kelp450_11_new/kelp450_11_new3.jl @@ -72,7 +72,9 @@ grid = RectilinearGrid(architecture, size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz)) # clock = Clock(; time = 0.0) T = FunctionField{Center, Center, Center}(t_function, grid; clock) -S = ConstantField(35.16) +# S = ConstantField(35.16) # has error: LoadError: type ConstantField has no field grid +S = FunctionField{Center, Center, Center}(s_function, grid; clock) + # ## Kelp Particle setup @info "Setting up kelp particles" n = 100 # number of kelp fronds @@ -156,7 +158,7 @@ progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: prettytime(sim.run_wall_time)) simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100)) -filename = "kelp450_11_new2" +filename = "kelp450_11_new3" simulation.output_writers[:profiles] = JLD2OutputWriter(model, model.tracers, filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing=true) #merge(model.tracers, model.auxiliary_fields), simulation.output_writers[:particles] = JLD2OutputWriter(model, (; kelp_particles), filename = "$(filename)_particles.jld2", schedule = TimeInterval(1day), overwrite_existing = true) @@ -174,7 +176,7 @@ simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateS #wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 2.0, min_change = 0.5, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) #simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) -#simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=SpecifiedTimes([i*years for i=1:12]), prefix=filename) #prefix="kelp_01" +simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=SpecifiedTimes([i*years for i=1:12]), prefix=filename) #prefix="kelp_01" # ## Run! # Finally we run the simulation run!(simulation) diff --git a/subpolarnew/kelp450_11_new/kelp450_11_new4.jl b/subpolarnew/kelp450_11_new/kelp450_11_new4.jl new file mode 100644 index 000000000..8b4a82228 --- /dev/null +++ b/subpolarnew/kelp450_11_new/kelp450_11_new4.jl @@ -0,0 +1,201 @@ +# # Simple active particle example +# In this example we will setup a simple 1D column with the [LOBSTER](@ref LOBSTER) biogeochemical model and active particles modelling the growth of sugar kelp. This demonstraits: +# - How to setup OceanBioME's biogeochemical models +# - How to setup light attenuation +# - How to add biologically active particles which interact with the biodeochemical model +# - How to include optional tracer sets (carbonate chemistry and oxygen) +# - How to visulise results + +# This is forced by idealised mixing layer depth and surface photosynthetically available radiation (PAR) which are setup first + +# ## Install dependencies +# First we will check we have the dependencies installed +# ```julia +# using Pkg +# pkg"add OceanBioME, Oceananigans, Printf, GLMakie" +# ``` + +# ## Model setup +# We load the packages and choose the default LOBSTER parameter set +using OceanBioME, Oceananigans,Printf +using Oceananigans.Fields: FunctionField, ConstantField +# using Oceananigans.Units: second, minute, minutes, hour, hours, day, days#, year, years +using Oceananigans.Units #####changed +using Interpolations +using JLD2 +using Oceananigans.Architectures: arch_array #####changed +const year = years = 365days #####changed + +# ## 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) + +##PAR⁰(x, y, t) = 60*(1-cos((t+15days)*2π/(365days)))*(1 /(1 +0.2*exp(-((mod(t, 365days)-200days)/50days)^2))) .+ 2 +##H(t, t₀, t₁) = ifelse(t₀ keys(pp["timeseries"]) 15-element Vector{String}: "NO₃" "NH₄" "P" "Z" "sPON" "bPON" "DON" "DIC" "Alk" "O₂" "sPOC" "bPOC" "DOC" "PAR" "t"# + +# ## Simulation +# Next we setup the simulation along with some callbacks that: +# - Couples the particles to the biodeochemical model +# - Update the PAR field from the surface PAR and phytoplankton concentration +# - Show the progress of the simulation +# - Store the model and particles output +# - Prevent the tracers from going negative from numerical error (see discussion of this in the [positivity preservation](@ref pos-preservation) implimentation page) + +simulation = Simulation(model, Δt=3.5minutes, stop_time=duration) + +# simulation.callbacks[:couple_particles] = Callback(Particles.infinitesimal_particle_field_coupling!; callsite = TendencyCallsite()) + +progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", + iteration(sim), + prettytime(sim), + prettytime(sim.Δt), + prettytime(sim.run_wall_time)) +simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100)) + +filename = "kelp450_11_new4" +simulation.output_writers[:profiles] = JLD2OutputWriter(model, model.tracers, filename = "$filename.jld2", schedule = TimeInterval(1day), overwrite_existing=true) #merge(model.tracers, model.auxiliary_fields), +simulation.output_writers[:particles] = JLD2OutputWriter(model, (; kelp_particles), filename = "$(filename)_particles.jld2", schedule = TimeInterval(1day), overwrite_existing = true) + + +#simulation.callbacks[:timestep] = Callback(update_timestep!, IterationInterval(1), (c_forcing=0.1, c_adv=0.5, c_diff=0.5, w = 200/day, relaxation=0.95), TimeStepCallsite()) + +scale_negative_tracers = ScaleNegativeTracers((:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) #ScaleNegativeTracers(;model, tracers = (:NO₃, :NH₄, :P, :Z, :sPOM, :bPOM, :DOM)) changed (:NO₃, :NH₄, :P, :Z, :sPON, :bPON, :DON) +simulation.callbacks[:neg] = Callback(scale_negative_tracers; callsite = UpdateStateCallsite()) + +# plankton_redfield = model.biogeochemistry.phytoplankton_redfield +# scale_negative_carbon_tracers = ScaleNegativeTracers(tracers = (:P, :Z, :DOC, :sPOC, :bPOC, :DIC), scalefactors = (P = plankton_redfield, Z = plankton_redfield, DOC = 1, sPOC = 1, bPOC = 1, DIC = 1)) +# simulation.callbacks[:neg2] = Callback(scale_negative_carbon_tracers; callsite = UpdateStateCallsite()) + + +#wizard = TimeStepWizard(cfl = 0.2, diffusive_cfl = 0.2, max_change = 2.0, min_change = 0.5, cell_diffusion_timescale = column_diffusion_timescale, cell_advection_timescale = column_advection_timescale) +#simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10)) + +simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=SpecifiedTimes([i*years for i=1:12]), prefix=filename) #prefix="kelp_01" +# ## Run! +# Finally we run the simulation +run!(simulation) + + + + + + + + + + + + + + + + + +