### Install CairoMakie for Visualization

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/quick_start.md

Installs the CairoMakie package, which is used for plotting the simulation results.

```julia
Pkg.add("CairoMakie")
```

--------------------------------

### Install OceanBioME and Oceananigans Dependencies

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/quick_start.md

Installs the necessary Julia packages, OceanBioME and Oceananigans, for running biogeochemical simulations.

```julia
using Pkg
Pkg.add(["OceanBioME", "Oceananigans"])
```

--------------------------------

### Creating a Minimal Working Example

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/contributing.md

This snippet demonstrates how to format code for bug reports using triple backticks, as recommended for creating minimal working examples. This helps in quickly identifying and resolving issues.

```Markdown
``````javascript
```some_code; and_some_more_code;```
``````
```

--------------------------------

### Set Up and Run a Biogeochemical Simulation

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/quick_start.md

Configures a RectilinearGrid, initializes a NonhydrostaticModel with the LOBSTER biogeochemistry model, sets initial tracer conditions, and runs a one-month simulation with constant forcing. Outputs are saved to a JLD2 file.

```julia
using OceanBioME, Oceananigans
using Oceananigans.Units

grid = RectilinearGrid(size = 10, extent = 200meters, topology = (Flat, Flat, Bounded))

model = NonhydrostaticModel(; grid, biogeochemistry = LOBSTER(; grid))

set!(model, P = 0.001, Z = 0.001, NO₃ = 1, NH₄ = 0.01)

simulation = Simulation(model, Δt = 1minute, stop_time = 30days)

simulation.output_writers[:profiles] = JLD2Writer(model, model.tracers,
                                                  filename = "quickstart.jld2",
                                                  schedule = TimeInterval(0.5days),
                                                  overwrite_existing = true)

run!(simulation)
```

--------------------------------

### Instantiate Project Dependencies

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/contributing.md

This command, run within the Julia environment, installs all the necessary dependencies listed in the Project.toml file for OceanBioME.jl.

```Julia
julia --project
] instantiate
```

--------------------------------

### Visualize Simulation Results with CairoMakie

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/quick_start.md

Loads simulation data (Phytoplankton and Nitrate concentrations) from a JLD2 file and visualizes their temporal evolution against depth using heatmaps. Requires the CairoMakie package.

```julia
using CairoMakie

phytoplankton = FieldTimeSeries("quickstart.jld2", "P")
nitrates = FieldTimeSeries("quickstart.jld2", "NO₃")

_, _, z = nodes(nitrates)

fig = Figure()

axis_kwargs = (xlabel = "Day", ylabel = "Depth (m)")
ax1 = Axis(fig[1, 1]; title = "Phytoplankton (mmol N/m³)", axis_kwargs...)
ax2 = Axis(fig[1, 2]; title = "Nitrate (mmol N/m³)", axis_kwargs...)

hm1 = heatmap!(ax1, phytoplankton.times / day, z, interior(phytoplankton , 1, 1, :, :)')
hm2 = heatmap!(ax2,      nitrates.times / day, z, interior(nitrates, 1, 1, :, :)')

fig
```

--------------------------------

### Clone OceanBioME Repository

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/contributing.md

This command clones your fork of the OceanBioME repository to your local machine. Ensure you replace 'your-user-name' with your actual GitHub username.

```Git
git clone https://github.com/your-user-name/OceanBioME.jl.git
```

--------------------------------

### Install CairoMakie for Visualization

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/README.md

Installs the CairoMakie package, which is a dependency for generating visualizations in OceanBioME.jl. This is typically the first step before running visualization code.

```Julia
using Pkg; Pkg.add("CairoMakie")
```

--------------------------------

### Run OceanBioME Tests

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/contributing.md

This command executes all the tests defined for the OceanBioME.jl package. This process may take some time to complete.

```Julia
] test
```

--------------------------------

### Install OceanBioME.jl

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/README.md

This snippet shows how to install the OceanBioME.jl package using the Julia package manager. It requires opening the Julia REPL and executing add commands.

```Julia
using Pkg
julia> Pkg.add("OceanBioME")
```

--------------------------------

### Install Oceananigans

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/index.md

Installs the Oceananigans package using Julia's package manager. This is a prerequisite for running the NPZD model example.

```julia
using Pkg; Pkg.add("Oceananigans")
```

--------------------------------

### Install OceanBioME.jl Package

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/index.md

This snippet demonstrates how to install the OceanBioME.jl package using the Julia package manager. It requires launching Julia and executing add commands.

```Julia
using Pkg
Pkg.add("OceanBioME")
```

--------------------------------

### Build OceanBioME Documentation

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/contributing.md

This command builds the documentation for OceanBioME.jl locally. The JULIA_DEBUG=Documenter environment variable provides additional information during the build process, which can be helpful for debugging.

```Julia
JULIA_DEBUG=Documenter julia --project=docs/ docs/make.jl
```

--------------------------------

### Setup SimpleMultiGSediment model

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/sediments/simple_multi_g.md

Demonstrates how to initialize the SimpleMultiGSediment model with a RectilinearGrid. It shows the creation of the grid and the sediment model, along with the resulting biogeochemical sediment type and its prognostic, tracked, and coupled fields.

```jldoctest
using OceanBioME, Oceananigans, OceanBioME.Sediments

grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200))

sediment_model = SimpleMultiGSediment(grid)

# output
`BiogeochemicalSediment` with `Single-layer multi-G sediment model (Float64)` biogeochemsitry
    Prognostic fields: (:Ns, :Nf, :Nr)
    Tracked fields: (:NO₃, :NH₄, :O₂, :sPOM, :bPOM)
    Coupled fields: (:NO₃, :NH₄, :O₂)
```

--------------------------------

### Setup CO2 and O2 Gas Exchange Boundary Conditions

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/air-sea-gas.md

Demonstrates how to set up and apply boundary conditions for carbon dioxide and oxygen gas exchange in an OceanBioME.jl model. It initializes the gas exchange conditions and then integrates them into a NonhydrostaticModel with specified biogeochemistry and grid.

```Julia
using OceanBioME
CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition()
O₂_flux  = OxygenGasExchangeBoundaryCondition()
using Oceananigans

grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200));

model = NonhydrostaticModel(; grid,
                              biogeochemistry = LOBSTER(; grid, carbonates = true, oxygen = true),
                              boundary_conditions = (DIC = FieldBoundaryConditions(top = CO₂_flux),
                                                      O₂ = FieldBoundaryConditions(top =  O₂_flux)),
                              tracers = (:T, :S))
```

--------------------------------

### Setup Oceananigans with BGC Model

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/biogeochemical/index.md

Demonstrates how to initialize an Oceananigans NonhydrostaticModel by passing a specified Biogeochemical (BGC) model. This includes the basic structure and optional parameters that can be provided to the BGC model.

```julia
model = NonhydrostaticModel(; grid,
                              ...,
                              biogeochemistry = MODEL_NAME(; grid))
MODEL_NAME(; grid, growth_rate = 10.0)
```

--------------------------------

### Setup PISCES Biogeochemistry in OceanBioME

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/biogeochemical/PISCES/PISCES.md

Demonstrates how to set up the default PISCES biogeochemical model with 24 tracers in OceanBioME using a RectilinearGrid. It also shows how to display the biogeochemistry object and list its required tracers.

```Julia
using OceanBioME, Oceananigans

grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1));

biogeochemistry = PISCES(; grid)

show(biogeochemistry)

Oceananigans.Biogeochemistry.required_biogeochemical_tracers(biogeochemistry)

(:P, :PChl, :PFe, :D, :DChl, :DFe, :DSi, :Z, :M, :DOC, :POC, :GOC, :SFe, :BFe, :PSi, :CaCO₃, :NO₃, :NH₄, :PO₄, :Fe, :Si, :DIC, :Alk, :O₂, :T, :S)
```

--------------------------------

### Near-Global Proof of Concept with OceanBioME.jl

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/paper/paper.md

This example presents a near-global proof of concept for OceanBioME.jl, focusing on light attenuation and the NPZD model. It showcases the package's applicability to larger-scale oceanographic simulations.

```Julia
using OceanBioME

# Example setup for near-global proof of concept
# ... (details of model configuration and parameters) ...

# Light attenuation
# ...

# NPZD model
model = NutrientPhytoplanktonZooplanktonDetritus()

# Further simulation steps...

```

--------------------------------

### Forking and Editing Code

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/contributing.md

This section describes the process of contributing code changes by forking the repository, editing code using git, and creating a pull request. It also mentions the use of the GitHub editor for minor changes.

```Git
git clone https://github.com/OceanBioME/OceanBioME.jl.git
# Make your changes
git commit -m "Your descriptive commit message"
git push origin your-branch-name
```

--------------------------------

### Run NPZD Model with Oceananigans.jl

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/README.md

This code demonstrates setting up and running a Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model using OceanBioME.jl coupled with Oceananigans.jl. It includes grid definition, biogeochemistry setup, model initialization, setting initial conditions, and configuring simulation output.

```Julia
using Pkg; Pkg.add("Oceananigans")

using OceanBioME, Oceananigans
using Oceananigans.Units

grid = RectilinearGrid(CPU(), size = (160, 32), extent = (10000meters, 500meters), topology = (Bounded, Flat, Bounded))

biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus(; grid) 

model = NonhydrostaticModel(; grid, biogeochemistry,
                              advection = WENO(),
                              closure = AnisotropicMinimumDissipation(),
			      buoyancy = SeawaterBuoyancy(constant_salinity = true))

@inline front(x, z, μ, δ) = μ + δ * tanh((x - 7000 + 4 * z) / 500)

Pᵢ(x, z) = ifelse(z > -50, 0.03, 0.01)
Nᵢ(x, z) = front(x, z, 2.5, -2)
Tᵢ(x, z) = front(x, z, 9, 0.05)

set!(model, N = Nᵢ, P = Pᵢ, Z = Pᵢ, T = Tᵢ)

simulation = Simulation(model; Δt = 50, stop_time = 4days)

simulation.output_writers[:tracers] = JLD2Writer(model, model.tracers,
                                                 filename = "buoyancy_front.jld2",
                                                 schedule = TimeInterval(24minute),
                                                 overwrite_existing = true)

run!(simulation)
```

--------------------------------

### GPU Grid Construction

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/README.md

Demonstrates how to construct a RectilinearGrid on the GPU for accelerated computations. This involves specifying the GPU device during grid creation, with the rest of the simulation setup remaining unchanged.

```Julia
grid = RectilinearGrid(GPU(), size = (256, 32), extent = (500meters, 100meters), topology = (Bounded, Flat, Bounded))
```

--------------------------------

### Add Upstream Remote

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/contributing.md

This command connects your local repository to the main OceanBioME project repository, allowing you to fetch updates from the upstream. 'oceanbiome' is used as a shorthand for the upstream remote.

```Git
git remote add oceanbiome https://github.com/OceanBioME/OceanBioME.jl.git
```

--------------------------------

### GPU Acceleration with CUDA.jl and KernelAbstractions.jl

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/paper/paper.md

This example illustrates the use of CUDA.jl and KernelAbstractions.jl to enable GPU acceleration for Oceananigans.jl and OceanBioME.jl. It shows how the underlying framework is designed to leverage GPU power for faster computations, a key feature for large-scale ocean modeling.

```julia
using Oceananigans
using OceanBioME
using CUDA
using KernelAbstractions

# Ensure CUDA is available and set the default device
if CUDA.has_cuda()
    CUDA.device!(0) # Select the first available GPU
    @info "Running on CUDA-enabled GPU"
else
    @info "CUDA not available, running on CPU"
end

# Oceananigans and OceanBioME are built using KernelAbstractions.jl,
# which allows them to run on different backends (CPU, GPU).
# When CUDA is available, computations are automatically dispatched to the GPU.

# Example of setting up a model that will run on GPU if CUDA is available:
# The model setup itself doesn't change drastically, but the execution context does.

# Define model parameters
# ...

# Create a model instance. If CUDA is available, this will utilize the GPU backend.
# The underlying KernelAbstractions.jl framework handles the dispatch.
model = HydrostaticModel(
    grid = LatitudeLongitudeGrid(architecture = CPU(), size = (10, 10, 10), lat = (-90, 90), lon = (0, 360), depth = (0, 100)),
    # ... other model parameters ...
)

# If you explicitly want to ensure GPU usage (assuming CUDA is available):
# model = HydrostaticModel(
#     grid = LatitudeLongitudeGrid(architecture = GPU(), size = (10, 10, 10), lat = (-90, 90), lon = (0, 360), depth = (0, 100)),
#     # ... other model parameters ...
# )

# When OceanBioME.jl is used with this model, its kernels will also run on the GPU.
# The integration is seamless due to the shared backend.

# Simulate the model
# simulation = Simulation(model, Δt=1, stop_time=10)
# run!(simulation)

```

--------------------------------

### Initialize CarbonChemistry Model

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/carbon-chemistry.md

Demonstrates how to initialize the CarbonChemistry model from OceanBioME.jl. This is the first step to using the carbon chemistry calculations.

```Julia
using OceanBioME

carbon_chemistry = CarbonChemistry()
```

--------------------------------

### Initialize CarbonChemistry with TEOS-10 Density

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/carbon-chemistry.md

Shows how to initialize the CarbonChemistry model using the TEOS-10 standard for seawater density calculations, which offers improved accuracy at a higher computational cost.

```Julia
using OceanBioME.Models: teos10_density
carbon_chemistry = CarbonChemistry(; density_function = teos10_density)
```

--------------------------------

### Set up and Run Oceananigans Simulation with OceanBioME

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_implementation.md

This snippet demonstrates setting up and running an Oceananigans simulation using the OceanBioME package. It includes defining a Photosynthetically Active Radiation (PAR) profile, temperature function, biogeochemistry model, and simulation parameters, then running the simulation and saving output.

```Julia
using OceanBioME, Oceananigans.Units
using Oceananigans.Fields: FunctionField

const year = years = 365days

@inline PAR⁰(t) = 500 * (1 - cos((t + 15days) * 2π / year)) * (1 / (1 + 0.2 * exp(-((mod(t, year) - 200days) / 50days)^2))) + 2

clock = Clock(; time = 0.0)

z = -10 # specify the nominal depth of the box for the PAR profile
@inline PAR_func(t) = PAR⁰(t) * exp(0.2z) # Modify the PAR based on the nominal depth and exponential decay 

PAR = FunctionField{Center, Center, Center}(PAR_func, BoxModelGrid(); clock)

@inline temp(t) = 2.4 * cos(t * 2π / year + 50days) + 26

biogeochemistry = Biogeochemistry(NutrientPhytoplankton(); 
                                  light_attenuation = PrescribedPhotosyntheticallyActiveRadiation(PAR))

model = BoxModel(; biogeochemistry,
                   prescribed_tracers = (; T = temp),
                   clock)

set!(model, N = 15, P = 15)

simulation = Simulation(model; Δt = 5minutes, stop_time = 5years)

simulation.output_writers[:fields] = JLD2Writer(model, model.fields; filename = "box_np.jld2", schedule = TimeInterval(10days), overwrite_existing = true)

# ## Run the model (should only take a few seconds)
@info "Running the model..."
run!(simulation)
```

--------------------------------

### Initialize and Set Particle Properties in Julia

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/individuals/index.md

Demonstrates how to initialize a `BiogeochemicalParticles` object with a specified number of particles and a `GrowingParticles` biogeochemistry. It also shows how to set the initial size and positions of these particles within the simulation grid.

```Julia
using OceanBioME, Oceananigans

Lx, Ly, Lz = 100, 100, 100
grid = RectilinearGrid(; size = (8, 8, 8), extent = (Lx, Ly, Lz))

particles = BiogeochemicalParticles(10; grid, biogeochemistry = GrowingParticles(0.5))

set!(particles, S = 0.1, x = rand(10) * Lx, y = rand(10) * Ly, z = rand(10) * Lz)
```

--------------------------------

### Sub-mesoscale Eddy Simulation with OceanBioME.jl

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/paper/paper.md

This example demonstrates the simulation of a sub-mesoscale eddy using the OceanBioME.jl package. It utilizes the LOBSTER biogeochemical model with an active carbonate model, CO2 exchange with the air, light attenuation, and mass-conserving negativity protection.

```Julia
using OceanBioME

# Example setup for sub-mesoscale eddy simulation
# ... (details of model configuration and parameters) ...

# LOBSTER biogeochemical model
model = LOBSTER()

# Carbonate model active
# ...

# CO2 exchange with the air
# ...

# Light attenuation
# ...

# Mass conserving negativity protection
# ...

# Further simulation steps...

```

--------------------------------

### Run Ocean Model Simulation and Set Output Writers in Julia

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_implementation.md

This snippet sets up and runs the ocean model simulation. It defines the time step, stop time, and configures output writers to save tracer and sediment data to JLD2 files at specified intervals.

```Julia
# run

simulation = Simulation(model, Δt = 9minutes, stop_time = 1years)

simulation.output_writers[:tracers] = JLD2Writer(model, model.tracers,
                                                 filename = "column_np.jld2",
                                                 schedule = TimeInterval(1day),
                                                 overwrite_existing = true)

simulation.output_writers[:sediment] = JLD2Writer(model, model.biogeochemistry.sediment.fields,
                                                  indices = (:, :, 1),
                                                  filename = "column_np_sediment.jld2",
                                                  schedule = TimeInterval(1day),
                                                  overwrite_existing = true)

run!(simulation)
```

--------------------------------

### Load and Visualize Ocean Data with CairoMakie

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/README.md

Loads time-series data for Temperature (T), Nutrients (N), and Phytoplankton (P) from JLD2 files. It then sets up a visualization using CairoMakie to create a GIF animation of these fields over time, displaying temperature, nutrient concentration, and phytoplankton concentration.

```Julia
T = FieldTimeSeries("buoyancy_front.jld2", "T")
N = FieldTimeSeries("buoyancy_front.jld2", "N")
P = FieldTimeSeries("buoyancy_front.jld2", "P")

xc, yc, zc = nodes(T)

times = T.times

using CairoMakie

n = Observable(1)

T_lims = (8.94, 9.06)
N_lims = (0, 4.5)
P_lims = (0.007, 0.02)

Tₙ = @lift interior(T[$n], :, 1, :)
Nₙ = @lift interior(N[$n], :, 1, :)
Pₙ = @lift interior(P[$n], :, 1, :)

fig = Figure(size = (1000, 520), fontsize = 20)

title = @lift "t = $(prettytime(times[$n]))"
Label(fig[0, :], title)

axis_kwargs = (xlabel = "x (m)", ylabel = "z (m)", width = 770, yticks = [-400, -200, 0])
ax1 = Axis(fig[1, 1]; title = "Temperature (°C)", axis_kwargs...)
ax2 = Axis(fig[2, 1]; title = "Nutrients concentration (mmol N / m³)",axis_kwargs...)
ax3 = Axis(fig[3, 1]; title = "Phytoplankton concentration (mmol N / m³)", axis_kwargs...)

hm1 = heatmap!(ax1, xc, zc, Tₙ, colorrange = T_lims, colormap = Reverse(:lajolla), interpolate = true)
hm2 = heatmap!(ax2, xc, zc, Nₙ, colorrange = N_lims, colormap = Reverse(:bamako), interpolate = true)
hm3 = heatmap!(ax3, xc, zc, Pₙ, colorrange = P_lims, colormap = Reverse(:bamako), interpolate = true)

Colorbar(fig[1, 2], hm1, ticks = [8.95, 9.0, 9.05])
Colorbar(fig[2, 2], hm2, ticks = [0, 2, 4])
Colorbar(fig[3, 2], hm3, ticks = [0.01, 0.02, 0.03])

rowgap!(fig.layout, 0)

record(fig, "buoyancy_front.gif", 1:length(times)) do i
    n[] = i
end
```

--------------------------------

### Idealised 1D Kelp Forest Model with OceanBioME.jl

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/paper/paper.md

This example details an idealized 1D model simulating kelp individuals, specifically Saccharina Latissima (sugar kelp), using OceanBioME.jl. It incorporates the LOBSTER biogeochemical model with an active carbonate model, variable Redfield ratio, CO2 exchange, light attenuation, and mass-conserving negativity protection.

```Julia
using OceanBioME

# Example setup for idealized 1D kelp forest model
# ... (details of model configuration and parameters) ...

# LOBSTER biogeochemical model with variable Redfield ratio
model = LOBSTER(redfield = VariableRedfield())

# Carbonate model active
# ...

# CO2 exchange with the air
# ...

# Light attenuation
# ...

# Mass conserving negativity protection
# ...

# Saccharina Latissima (sugar kelp) model
# ...

# Further simulation steps...

```

--------------------------------

### Configure Biogeochemical Model and Closure in Julia

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_implementation.md

This code configures the biogeochemical model, including light attenuation, sediment dynamics, and nutrient-phytoplankton interactions. It also sets up the closure model for diffusion and initializes the model state.

```Julia
# setup the biogeochemical model

light_attenuation = TwoBandPhotosyntheticallyActiveRadiation(; grid, surface_PAR)

sediment = InstantRemineralisationSediment(grid; sinking_tracers = :P)

sinking_velocity = ZFaceField(grid)

w_sink(z) = 2 / day * tanh(z / 5)

set!(sinking_velocity, w_sink)

negative_tracer_scaling = ScaleNegativeTracers((:N, :P))

biogeochemistry = Biogeochemistry(NutrientPhytoplankton(; sinking_velocity);
                                  light_attenuation,
                                  sediment,
                                  modifiers = negative_tracer_scaling) 

κ = CenterField(grid)

set!(κ, κₚ)

# put the model together

model = NonhydrostaticModel(; grid,
                              biogeochemistry,
                              closure = ScalarDiffusivity(ν = κ; κ), 
                              forcing = (; T = ∂ₜT))

set!(model, P = 0.01, N = 15, T = 28)
```

--------------------------------

### Visualize Model Output

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/index.md

Loads simulation output from a JLD2 file and visualizes the Temperature, Nutrients, and Phytoplankton concentrations over time using CairoMakie. It creates animated plots of the buoyancy front.

```julia
T = FieldTimeSeries("buoyancy_front.jld2", "T")
N = FieldTimeSeries("buoyancy_front.jld2", "N")
P = FieldTimeSeries("buoyancy_front.jld2", "P")

xc, yc, zc = nodes(T)

times = T.times

using CairoMakie

n = Observable(1)

T_lims = (8.94, 9.06)
N_lims = (0, 4.5)
P_lims = (0.007, 0.02)

Tₙ = @lift interior(T[$n], :, 1, :)
Nₙ = @lift interior(N[$n], :, 1, :)
Pₙ = @lift interior(P[$n], :, 1, :)

fig = Figure(size = (1000, 520), fontsize = 20)

title = @lift "t = $(prettytime(times[$n]))"
Label(fig[0, :], title)

axis_kwargs = (xlabel = "x (m)", ylabel = "z (m)", width = 770, yticks = [-400, -200, 0])
ax1 = Axis(fig[1, 1]; title = "Temperature (°C)", axis_kwargs...)
ax2 = Axis(fig[2, 1]; title = "Nutrients concentration (mmol N / m³)",axis_kwargs...)
ax3 = Axis(fig[3, 1]; title = "Phytoplankton concentration (mmol N / m³)", axis_kwargs...)

hm1 = heatmap!(ax1, xc, zc, Tₙ, colorrange = T_lims, colormap = Reverse(:lajolla), interpolate = true)
hm2 = heatmap!(ax2, xc, zc, Nₙ, colorrange = N_lims, colormap = Reverse(:bamako), interpolate = true)
hm3 = heatmap!(ax3, xc, zc, Pₙ, colorrange = P_lims, colormap = Reverse(:bamako), interpolate = true)

Colorbar(fig[1, 2], hm1, ticks = [8.95, 9.0, 9.05])
Colorbar(fig[2, 2], hm2, ticks = [0, 2, 4])
Colorbar(fig[3, 2], hm3, ticks = [0.01, 0.02, 0.03])

rowgap!(fig.layout, 0)

record(fig, "buoyancy_front.mp4", 1:length(times)) do i
    n[] = i
end

nothing #hide
```

--------------------------------

### Run NPZD Model Simulation

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/index.md

Sets up and runs a two-dimensional NPZD model simulation with a buoyancy front. It configures the grid, biogeochemistry, physics, initial conditions, simulation time, and output writing.

```julia
using OceanBioME, Oceananigans
using Oceananigans.Units

grid = RectilinearGrid(CPU(), size = (160, 32), extent = (10000meters, 500meters), topology = (Bounded, Flat, Bounded))

biogeochemistry = NutrientPhytoplanktonZooplanktonDetritus (; grid) 

model = NonhydrostaticModel(; grid, biogeochemistry,
                              advection = WENO(),
			      closure = AnisotropicMinimumDissipation(),
                              buoyancy = SeawaterBuoyancy(constant_salinity = true))

@inline front(x, z, μ, δ) = μ + δ * tanh((x - 7000 + 4 * z) / 500)

Pᵢ(x, z) = ifelse(z > -50, 0.03, 0.01)
Nᵢ(x, z) = front(x, z, 2.5, -2)
Tᵢ(x, z) = front(x, z, 9, 0.05)

set!(model, N = Nᵢ, P = Pᵢ, Z = Pᵢ, T = Tᵢ)

simulation = Simulation(model; Δt = 50, stop_time = 4days)

simulation.output_writers[:tracers] = JLD2Writer(model, model.tracers,
                                                 filename = "buoyancy_front.jld2",
                                                 schedule = TimeInterval(24minute),
                                                 overwrite_existing = true)

run!(simulation)
```

--------------------------------

### Create Sugar Kelp Particles

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/individuals/slatissima.md

Shows how to initialize SugarKelpParticles for use within the OceanBioME framework. This involves setting up a grid and specifying the number of particles. The particles are initialized with SugarKelp biogeochemistry and track specific fields and coupled tracers.

```jldoctest
using OceanBioME, Oceananigans

grid = RectilinearGrid(size = (1, 1, 1), extent = (1, 1, 1));
particles = SugarKelpParticles(10; grid)

# output
10 BiogeochemicalParticles with SugarKelp{Float64} biogeochemistry:
├── fields: (:A, :N, :C)
└── coupled tracers: (:NO₃, :NH₄, :DIC, :O₂, :DOC, :DON, :bPOC, :bPON)
```

--------------------------------

### Document Light Attenuation Models

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/appendix/library.md

This snippet documents the Light Attenuation Models in OceanBioME.jl. It utilizes the `@autodocs` macro for automatic documentation generation.

```Julia
```@autodocs
Modules = [OceanBioME.Light]
```
```

--------------------------------

### Compute pH using CarbonChemistry

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/carbon-chemistry.md

Demonstrates how to obtain the pH as an output from the CarbonChemistry model by setting the `return_pH` argument to true.

```Julia
carbon_chemistry(; DIC, Alk, T, S, return_pH = true)
```

--------------------------------

### Add GLMakie Package

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/visualization.md

This code snippet demonstrates how to add the GLMakie package to your Julia environment using the Pkg manager. GLMakie is a backend for Makie.jl that provides interactive plotting capabilities.

```julia
using Pkg
pkg"add GLMakie"
```

--------------------------------

### Configure LOBSTER biogeochemistry with sediment model

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/sediments/simple_multi_g.md

Illustrates how to set up the LOBSTER biogeochemical model, enabling carbonate chemistry, oxygen, variable Redfield ratios, and open bottom conditions, while integrating the previously defined SimpleMultiGSediment model.

```julia
biogeochemistry = LOBSTER(; grid,
                            carbonates = true, oxygen = true, variable_redfield = true,
                            open_bottom = true,
                            sediment_model)
```

--------------------------------

### Visualize Simulation Results with CairoMakie

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_implementation.md

This code snippet visualizes the results of an OceanBioME simulation using CairoMakie. It loads time-series data for nutrients (N) and phytoplankton (P) from a JLD2 file and plots them against time, along with the PAR and temperature profiles.

```Julia
using CairoMakie

P = FieldTimeSeries("box_np.jld2", "P")
N = FieldTimeSeries("box_np.jld2", "N")

times = P.times

# ## And plot

fig = Figure(size = (1200, 480), fontsize = 20)

axN= Axis(fig[1, 1], ylabel = "Nutrient \n(mmol N / m³)")
lines!(axN, times / year, N[1, 1, 1, :], linewidth = 3)

axP = Axis(fig[1, 2], ylabel = "Phytoplankton \n(mmol N / m³)")
lines!(axP, times / year, P[1, 1, 1, :], linewidth = 3)

axPAR= Axis(fig[2, 1], ylabel = "PAR (einstein / m² / s)", xlabel = "Time (years)")
lines!(axPAR, times / year, PAR_func.(times), linewidth = 3)

axT = Axis(fig[2, 2], ylabel = "Temperature (°C)", xlabel = "Time (years)")
lines!(axT, times / year, temp.(times), linewidth = 3)

fig
```

--------------------------------

### Initialize LOBSTER Model with Carbonates

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/biogeochemical/LOBSTER.md

Demonstrates how to initialize the LOBSTER biogeochemical model with carbonate chemistry enabled using OceanBioME and Oceananigans libraries in Julia. It sets up a rectilinear grid and configures the LOBSTER model with specific parameters.

```julia
using OceanBioME, Oceananigans

grid = RectilinearGrid(size=(3, 3, 30), extent=(10, 10, 200))

bgc_model = LOBSTER(; grid, carbonates = true)

```

--------------------------------

### Visualize Ocean Model Simulation Results in Julia

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_implementation.md

This code visualizes the simulation results by loading tracer and sediment data from JLD2 files and plotting them using heatmaps and line plots. It sets up axes, plots the data, and adds colorbars for clarity.

```Julia
N = FieldTimeSeries("column_np.jld2", "N")
P = FieldTimeSeries("column_np.jld2", "P")

sed = FieldTimeSeries("column_np_sediment.jld2", "storage")

fig = Figure()

axN = Axis(fig[1, 1], ylabel = "z (m)")
axP = Axis(fig[2, 1], ylabel = "z (m)")
axSed = Axis(fig[3, 1:2], ylabel = "Sediment (mmol N / m²)", xlabel = "Time (years)")

_, _, zc = nodes(grid, Center(), Center(), Center())
times = N.times

hmN = heatmap!(axN, times ./ year, zc, N[1, 1, 1:grid.Nz, 1:end]',
               interpolate = true, colormap = Reverse(:batlow))

hmP = heatmap!(axP, times ./ year, zc, P[1, 1, 1:grid.Nz, 1:end]',
               interpolate = true, colormap = Reverse(:batlow))

lines!(axSed, times ./ year, sed[1, 1, 1, :])

Colorbar(fig[1, 2], hmN, label = "Nutrient (mmol N / m³)")
Colorbar(fig[2, 2], hmP, label = "Phytoplankton (mmol N / m³)")

fig
```

--------------------------------

### Define NutrientPhytoplankton Model Structure

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_implementation.md

Defines the `NutrientPhytoplankton` struct, inheriting from `AbstractContinuousFormBiogeochemistry`. It includes parameters for growth rate, nutrient and light saturation, temperature effects, mortality, and sinking velocity, with default values provided using `@kwdef`.

```Julia
@kwdef struct NutrientPhytoplankton{FT, W} <: AbstractContinuousFormBiogeochemistry
            base_growth_rate :: FT = 1.27 / day              # 1 / seconds
    nutrient_half_saturation :: FT = 0.025 * 1000 / 14       # mmol N / m³
       light_half_saturation :: FT = 300.0                   # micro einstein / m² / s
        temperature_exponent :: FT = 0.24                    # 1
     temperature_coefficient :: FT = 1.57                    # 1
         optimal_temperature :: FT = 28.0                    # °C
              mortality_rate :: FT = 0.15 / day              # 1 / seconds
     crowding_mortality_rate :: FT = 0.004 / day / 1000 * 14 # 1 / seconds / mmol N / m³
            sinking_velocity :: W  = 2 / day
end
```

--------------------------------

### Julia: P-I Slope Adaptation

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/biogeochemical/PISCES/notable_differences.md

Explains how the P-I slope (α) can be adapted to depth in the Oceanbiome Julia project by setting 'low_light_adaptation' in growth rate parameterizations.

```Julia
low_light_adaptation
```

--------------------------------

### Extending OceanBioME.jl with New Formulations

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/paper/paper.md

OceanBioME.jl facilitates rapid prototyping by allowing models to be implemented and swapped easily through the extension of key functions. This flexibility enables users to explore different biogeochemical model formulations without modifying other model configurations.

```Julia
# Extend existing model functions or define new ones
# For example, to add a new tracer or modify a process rate:
# function update_process_rates!(model, p, t)
#     # ... custom logic ...
# end
```

--------------------------------

### Import Oceananigans and OceanBioME

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_implementation.md

This is a basic import statement for the Oceananigans and OceanBioME packages, typically used at the beginning of a script or module to make their functionalities available.

```Julia
using Oceananigans, OceanBioME
```

--------------------------------

### Julia: Iron Equilibrium Constant Formula

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/biogeochemical/PISCES/notable_differences.md

Provides the formula for K_Fe_eq in equation 65, as taken from the NEMO source code and used in the Oceanbiome Julia project.

```Julia
K_Fe_eq = exp(16.27 - 1565.7 / max(T + 273.15, 5))
```

--------------------------------

### PISCES Model Unit Conversions and Parameters

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_components/biogeochemical/PISCES/notable_differences.md

This section details unit conversions and parameter value discrepancies between the OceanBioME PISCES implementation and other sources like Aumont2015 and NEMO. It covers aggregation factors, flux feeding rates, bacterial iron uptake efficiency, and maximum iron ratios for different biological components.

```Julia
aggregation_parameters: from the NEMO source code, all of these parameters need a factor of ``10^{-6}`` to convert them from units of 1 / (mol / L) to 1 / (μmol / L). Additionally, all the parameters require a factor of ``1 / 86400``, for the parameters not multiplied by shear this straight forwardly is because they have time units of 1 / day in the NEMO code, but for those multiplied by shear this is because they implicitly have a factor of seconds per day in them to result in an aggregation in mmol C / day
```

```Julia
flux feeding rate for zooplankton ``g_{FF}^M`` is in 1 / (m mol / L) where that `m` in `m`eters, so we need to multiply by a factor of ``10^{-6}`` to get 1 / (m μmol / L)
```

```Julia
fraction of bacterially consumed iron going to small and big particles, ``kappa^{BFe}_{Bact}`` and ``kappa^{SFe}_{Bact}``, in equations 48 and 49 are not recorded but from the NEMO source code we can see that they are `0.04` and `0.12` respectively. Additionally, we need to multiply by a factor of `0.16` (``bacterial_iron_uptake_efficiency``) to the drawdown from the iron pool due to the instantaneous return (as per the NEMO version)
```

```Julia
``	heta_{max}^{Fe, Bact}`` is not recorded so the value `0.06` μmol Fe / mmol C is taken from the NEMO source code
```

```Julia
``	heta^{Fe, Z}`` and ``	heta^{Fe, M}`` are taken from the NEMO source code to be 0.01 and 0.015 μmol Fe / mmol C
```

```Julia
``	heta_{max}^{Fe, P}`` is taken from the NEMO source code to be `0.06` μmol Fe / mmol C, we note that this isn't actually the maximum in the sense that the ratio could (probably) go above this value
```

```Julia
``K^{B, 1}_{Fe}`` is not recorded so the value `0.3` μmol Fe / m³ is taken from the NEMO source code
```

--------------------------------

### Specify Required Biogeochemical Tracers and Fields

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/model_implementation.md

Declares the required tracers (`:N`, `:P`, `:T`) and auxiliary fields (`:PAR`) for the `NutrientPhytoplankton` model. This informs OceanBioME and Oceananigans about the state variables and external inputs the model uses.

```Julia
required_biogeochemical_tracers(::NutrientPhytoplankton) = (:N, :P, :T)

required_biogeochemical_auxiliary_fields(::NutrientPhytoplankton) = (:PAR, )
```

--------------------------------

### Document Gas Exchange Models

Source: https://github.com/oceanbiome/oceanbiome.jl/blob/main/docs/src/appendix/library.md

This snippet documents the Gas exchange boundary conditions in OceanBioME.jl, including the main GasExchangeModel and ScaledGasTransferVelocity. It uses the `@autodocs` macro for documentation.

```Julia
```@autodocs
Modules = [OceanBioME.Models.GasExchangeModel, OceanBioME.Models.GasExchangeModel.ScaledGasTransferVelocity]
```
```