Mv functions in examples to library

README.md

@@ -28,6 +28,7 @@ To activate the project (this takes a while as it installs all the packages):
 ```bash
 julia --project -e 'using Pkg; Pkg.instantiate()'
 ```
+
 You can then use the package interactively, in the terminal:
 
 ```bash

examples/N2P2ZD/box_model.jl

@@ -16,7 +16,9 @@ const year = years = 365day
 # ==================================================
 
 include("tracers.jl")
-bgc_model = Biogeochemistry(N2P2ZD(); light_attenuation=FunctionPAR(; grid=BoxModelGrid()))
+bgc_model = Biogeochemistry(
+    N2P2ZD(); light_attenuation=FunctionFieldPAR(; grid=BoxModelGrid())
+)
 full_model = BoxModel(; biogeochemistry=bgc_model)
 set!(full_model; N=7.0, Z2=0.05, D=0.0, P1=0.1, P2=0.1, Z1=0.05)
 

examples/N2P2ZD/functions.jl

@@ -1,5 +1,5 @@
 # phytoplankton growth
-function menden_limitation(N::Real, nitrogen_half_saturation::Real)
+function monod_limitation(N::Real, nitrogen_half_saturation::Real)
     return N / (nitrogen_half_saturation + N)
 end
 function smith_light_limitation(PAR::Real, alpha::Real, maximum_growth_rate::Real)
@@ -17,7 +17,7 @@ function photosynthetic_growth(
     alpha::Real,
 )
     return maximum_growth_rate *
-           menden_limitation(N, nitrogen_half_saturation) *
+           monod_limitation(N, nitrogen_half_saturation) *
            smith_light_limitation(PAR, alpha, maximum_growth_rate) *
            P
 end

examples/NPZD/functions.jl

@@ -1,21 +1,22 @@
-menden_limitation(R, k) = R / (k + R)
+monod_limitation(R, k) = R / (k + R)
+holling_type_2(R, k) = R / (k + R)
 
 Q₁₀(T) = 1.88^(T / 10) # T in °C
 
 smith_light_limitation(PAR, α, μ₀) = α * PAR / sqrt(μ₀^2 + α^2 * PAR^2)
 
-remineralization(D, r) = r * D
+idealized_remineralization(D, r) = r * D
 
-function photosynthetic_growth(N, P, PAR, μ₀, kₙ, α)
-    return μ₀ * menden_limitation(N, kₙ) * smith_light_limitation(PAR, α, μ₀) * P
+function idealized_photosynthetic_growth(N, P, PAR, μ₀, kₙ, α)
+    return μ₀ * monod_limitation(N, kₙ) * smith_light_limitation(PAR, α, μ₀) * P
 end
 
-predation_loss(P, Z, gₘₐₓ, kₚ) = gₘₐₓ * menden_limitation(P^2, kₚ^2) * Z
+idealized_predation_loss(P, Z, gₘₐₓ, kₚ) = gₘₐₓ * holling_type_2(P^2, kₚ^2) * Z
 
-predation_gain(P, Z, β, gₘₐₓ, kₚ) = β * gₘₐₓ * menden_limitation(P^2, kₚ^2) * Z
+idealized_predation_gain(P, Z, β, gₘₐₓ, kₚ) = β * gₘₐₓ * holling_type_2(P^2, kₚ^2) * Z
 
-function predation_assimilation_loss(P, Z, β, gₘₐₓ, kₚ)
-    return (1 - β) * gₘₐₓ * menden_limitation(P^2, kₚ^2) * Z
+function idealized_predation_assimilation_loss(P, Z, β, gₘₐₓ, kₚ)
+    return (1 - β) * gₘₐₓ * holling_type_2(P^2, kₚ^2) * Z
 end
 
 linear_loss(P, l) = l * P

examples/NPZD/tracers.jl

@@ -19,20 +19,23 @@ parameters = (
 
 tracers = Dict(
     "N" => :(
-        summed_linear_loss([P, Z], lⁿ) + remineralization(D, rᵈⁿ) -
-        photosynthetic_growth(N, P, PAR, μ₀, kₙ, α)
+        summed_linear_loss([P, Z], lⁿ) + idealized_remineralization(D, rᵈⁿ) -
+        idealized_photosynthetic_growth(N, P, PAR, μ₀, kₙ, α)
     ),
     "D" => :(
         linear_loss(P, lᵖᵈ) +
-        predation_assimilation_loss(P, Z, β, gₘₐₓ, kₚ) +
-        quadratic_loss(Z, lᶻᵈ) - remineralization(D, rᵈⁿ)
+        idealized_predation_assimilation_loss(P, Z, β, gₘₐₓ, kₚ) +
+        quadratic_loss(Z, lᶻᵈ) - idealized_remineralization(D, rᵈⁿ)
     ),
     "P" => :(
-        photosynthetic_growth(N, P, PAR, μ₀, kₙ, α) - predation_loss(P, Z, gₘₐₓ, kₚ) -
-        linear_loss(P, lᵖⁿ) - linear_loss(P, lᵖᵈ)
+        idealized_photosynthetic_growth(N, P, PAR, μ₀, kₙ, α) -
+        idealized_predation_loss(P, Z, gₘₐₓ, kₚ) - linear_loss(P, lᵖⁿ) -
+        linear_loss(P, lᵖᵈ)
+    ),
+    "Z" => :(
+        idealized_predation_gain(P, Z, β, gₘₐₓ, kₚ) - linear_loss(Z, lᶻⁿ) -
+        quadratic_loss(Z, lᶻᵈ)
     ),
-    "Z" =>
-        :(predation_gain(P, Z, β, gₘₐₓ, kₚ) - linear_loss(Z, lᶻⁿ) - quadratic_loss(Z, lᶻᵈ)),
 )
 
 NPZD = define_tracer_functions(parameters, tracers; helper_functions="functions.jl")

examples/box_model.jl

@@ -16,7 +16,9 @@ const year = years = 365day
 # ==================================================
 
 include(joinpath("NPZD", "tracers.jl"))
-bgc_model = Biogeochemistry(NPZD(); light_attenuation=FunctionPAR(; grid=BoxModelGrid()))
+bgc_model = Biogeochemistry(
+    NPZD(); light_attenuation=FunctionFieldPAR(; grid=BoxModelGrid())
+)
 full_model = BoxModel(; biogeochemistry=bgc_model)
 set!(full_model; N=7.0, P=0.01, Z=0.05, D=0.0)
 

src/Agate.jl

@@ -3,6 +3,7 @@ module Agate
 include("Library/Library.jl")
 include("Models/Models.jl")
 
+using .Library
 using .Models
 
 export define_tracer_functions

src/Library/growth.jl

@@ -1,25 +1,24 @@
 """
-Modules related to phytoplankton photosynthetic growth
-
+Functions related to phytoplankton photosynthetic growth.
 """
+
 module Growth
 
 "
-    PCⱼ = PCⱼᵐᵃˣ * γⱼⁿᵘᵗ *  γⱼˡⁱᵍʰᵗ * fⱼᵗᵉᵐᵖ *  γⱼᶜᵒ²
-
-Carbon-specific growth rate for plankton j (Default MITgcm-DARWIN formulation).
+    PC = PCᵐᵃˣ * γⁿᵘᵗ *  γˡⁱᵍʰᵗ * fᵗᵉᵐᵖ *  γᶜᵒ²
 
-Where: 
-PCᵐᵃˣⱼ = maximum carbon-specific growth rate for plankton j,
-γⁿᵘᵗⱼ = nutrient limition,
-γˡⁱᵍʰᵗⱼ = light limition,
-fᵗᵉᵐᵖⱼ = temperature limitation,
-γᶜᵒ²ⱼ = carbon dioxide limitation
+Carbon-specific growth rate for plankton (Default MITgcm-DARWIN formulation).
 
+# Arguments
+- `PCᵐᵃˣ`: maximum carbon-specific growth rate
+- `γⁿᵘᵗ`: nutrient limition
+- `γˡⁱᵍʰᵗ`: light limition
+- `fᵗᵉᵐᵖ`: temperature limitation
+- `γᶜᵒ²`: carbon dioxide limitation
 "
-function default_PCⱼ(PCⱼᵐᵃˣ, γⱼⁿᵘᵗ, γⱼˡⁱᵍʰᵗ, fⱼᵗᵉᵐᵖ, γⱼᶜᵒ²)
-    return PCⱼᵐᵃˣ * γⱼⁿᵘᵗ * γⱼˡⁱᵍʰᵗ * fⱼᵗᵉᵐᵖ * γⱼᶜᵒ²
+function default_PC(PCᵐᵃˣ, γⁿᵘᵗ, γˡⁱᵍʰᵗ, fᵗᵉᵐᵖ, γᶜᵒ²)
+    return PCᵐᵃˣ * γⁿᵘᵗ * γˡⁱᵍʰᵗ * fᵗᵉᵐᵖ * γᶜᵒ²
 end
 
-export default_PCⱼ
+export default_PC
 end # module

src/Library/light.jl

@@ -1,5 +1,5 @@
 """
-Modules related to photosynthetically available radiation (PAR)
+Functions related to photosynthetically available radiation (PAR).
 """
 
 module Light
@@ -14,6 +14,8 @@ import Oceananigans.Biogeochemistry:
 
 const year = years = 365day
 
+export cyclical_PAR, FunctionFieldPAR
+
 """
     cyclical_PAR(t, z) -> Float
 
@@ -42,33 +44,31 @@ Light module for PAR defined by simple functions (can be used with box or column
 # Fields
 - `field`: Oceananigans.FunctionField
 """
-struct FunctionPAR
+struct FunctionFieldPAR
     field
 end
 
 """
-    FunctionPAR() -> DataType
+    FunctionFieldPAR() -> DataType
 
 # Keywords
 - `grid`: the geometry to build the model on defined as an Oceananigans grid object
 - `PAR_f`: a PAR function of time (and depth), defaults to `cyclical_PAR`
 """
-function FunctionPAR(; grid, PAR_f=cyclical_PAR(; z=-10))
+function FunctionFieldPAR(; grid, PAR_f=cyclical_PAR(; z=-10))
     clock = Clock(; time=zero(grid))
     PAR_field = FunctionField{Center,Center,Center}(PAR_f, grid; clock)
-    return FunctionPAR(PAR_field)
+    return FunctionFieldPAR(PAR_field)
 end
 
-function update_biogeochemical_state!(model, PAR::FunctionPAR)
+function update_biogeochemical_state!(model, PAR::FunctionFieldPAR)
     PAR.field.clock.time = model.clock.time
     compute!(PAR.field)
     return nothing
 end
 
-function biogeochemical_auxiliary_fields(par::FunctionPAR)
+function biogeochemical_auxiliary_fields(par::FunctionFieldPAR)
     return (PAR=par.field,)
 end
 
-export cyclical_PAR, FunctionPAR
-
 end # module

src/Library/mortality.jl

@@ -1,3 +1,27 @@
 module Mortality
 
+export linear_loss, quadratic_loss
+
+"""
+Linear mortality rate. 
+In this formulation mortality is constant, and can be interpreted as 
+a "closure term" for low density predation and and other death terms.
+
+# Arguments
+- `P`: plankton concentration
+- `l`: mortality rate
+"""
+linear_loss(P, l) = l * P
+
+"""
+Quadratic mortality coefficient.
+In this formulation mortality increases exponentially with plankton biomass
+and is often interpreted to represent viral processes and non-represented density-dependent predation effects.
+
+# Arguments
+- `P`: plankton concentration
+- `l`: mortality rate
+"""
+quadratic_loss(P, l) = l * P^2
+
 end # module

src/Library/nutrients.jl

@@ -1,24 +1,24 @@
 """
-Modules related to plankton nutrient uptake
-
+Functions related to plankton nutrient uptake.
 """
+
 module Nutrients
 
+export monod_limitation
+
 "
     R / (kᵣ + R)
 
 Monod formulation of nutrient limitation, which is based on Michaelis-Menten enzyme kinetics.
 
-Where:
-R = nutrient (e.g. N, P, Si)
-kᵣ = nutrient half saturation constant
-
-Note that sometimes this formulation is also used for Predation.
+# Arguments
+- `R`: nutrient (e.g. N, P, Si)
+- `kᵣ`: nutrient half saturation constant
 
+Note that sometimes this formulation is also used for Predation ('Holling type 2').
 "
 function monod_limitation(R, kᵣ)
     return R / (kᵣ + R)
 end
 
-export monod_limitation
 end # module

src/Library/photosynthesis.jl

@@ -1,40 +1,57 @@
 """
-Modules related to phytoplankton light uptake
-
+Functions related to phytoplankton light uptake.
 """
+
 module Photosynthesis
 
-"
-    γⱼˡⁱᵍʰᵗ = (1 - ℯ^(kⱼˢᵃᵗ*I)) * ℯ^kⱼⁱⁿʰ * nⱼˡⁱᵍʰᵗ
+export γˡⁱᵍʰᵗ, smith_light_limitation, idealized_photosynthetic_growth
 
-Light limitation for plankton j (Default MITgcm-DARWIN formulation).
+"""
+    γˡⁱᵍʰᵗ = (1 - ℯ^(kˢᵃᵗ*I)) * ℯ^kⁱⁿʰ * nˡⁱᵍʰᵗ
 
-Where:
-kⱼˢᵃᵗ = half saturation constant of light saturation of plankton j,
-I = irradiance,
-kⱼⁱⁿʰ = half saturation constant of light inhibition of plankton j,
-nⱼˡⁱᵍʰᵗ = light penalty term of plankton j
+Light limitation for plankton (Default MITgcm-DARWIN formulation).
 
-"
-function γⱼˡⁱᵍʰᵗ(I, kⱼˢᵃᵗ, kⱼⁱⁿʰ, nⱼˡⁱᵍʰᵗ)
-    return (1 - ℯ^(kⱼˢᵃᵗ * I)) * ℯ^kⱼⁱⁿʰ * nⱼˡⁱᵍʰᵗ
+# Arguments
+- `I`: irradiance
+- `kˢᵃᵗ`:  half saturation constant of light saturation
+- `kⁱⁿʰ`: half saturation constant of light inhibition
+- `nˡⁱᵍʰᵗ`: light penalty term
+"""
+function γˡⁱᵍʰᵗ(I, kˢᵃᵗ, kⁱⁿʰ, nˡⁱᵍʰᵗ)
+    return (1 - ℯ^(kˢᵃᵗ * I)) * ℯ^kⁱⁿʰ * nˡⁱᵍʰᵗ
 end
 
 "
     α * PAR / sqrt(μ₀ ^ 2 + α ^ 2 * PAR ^ 2)
 
 Smith 1936 formulation of light limitation (also see Evans and Parslow, 1985).
 
-Where:
-α = Initial photosynthetic slope
-PAR = Photosynthetic Active Radiation
-μ₀ = Maximum growth rate at T = 0 °C (this seems weird?, from Kuhn 2015)
-
+# Arguments
+- `PAR`: photosynthetic active radiation
+- `α`: initial photosynthetic slope
+- `μ₀`: maximum growth rate at T = 0 °C (this seems weird?, from Kuhn 2015)
 "
 function smith_light_limitation(PAR, α, μ₀)
+    # here to avoid division by 0 when α and μ₀ are both 0
+    if alpha == 0
+        return 0.0
+    end
     return α * PAR / sqrt(μ₀^2 + α^2 * PAR^2)
 end
 
-export γⱼˡⁱᵍʰᵗ
-smith_light_limitation
+"""
+Single nutrient monod smith photosynthetic growth (used, for example, in Kuhn 2015).
+
+# Arguments
+- `N`: nutrient concentration
+- `P`: phytoplankton concentration
+- `PAR`: photosynthetic active radiation
+- `α`: initial photosynthetic slope
+- `μ₀`: maximum growth rate at T = 0 °C
+- `kₙ`: nutrient half saturation
+"""
+function idealized_photosynthetic_growth(N, P, PAR, μ₀, kₙ, α)
+    return μ₀ * monod_limitation(N, kₙ) * smith_light_limitation(PAR, α, μ₀) * P
+end
+
 end # module