GitHub

Biogeochemistry in submesoscale eddies in the Eady model

In this example we setup a 3D model with a constant background buoyancy gradient with associated thermal wind (the Eady model) with the LOBSTER biogeochemical model. This example demonstrates how to use biogeochemistry in a more complicated physical model. The parameters in this example correspond roughly to those used by Taylor (2016) and result to the generation of a single submesoscale eddy.

Install dependencies

First we ensure we have the required dependencies installed

using Pkg
pkg "add OceanBioME, Oceananigans, CairoMakie"

Model setup

We load the required packages. Although not required, we also set the random seed to ensure reproducibility of the results.

using OceanBioME, Oceananigans, Printf, CairoMakie
using Oceananigans.Units

using Random
Random.seed!(11)
Random.TaskLocalRNG()

Construct a grid with uniform grid spacing on CPU

grid = RectilinearGrid(size = (32, 32, 8),
                       extent = (1kilometer, 1kilometer, 100meters))
32×32×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1000.0) regularly spaced with Δx=31.25
├── Periodic y ∈ [0.0, 1000.0) regularly spaced with Δy=31.25
└── Bounded  z ∈ [-100.0, 0.0] regularly spaced with Δz=12.5

Set the Coriolis and buoyancy models.

coriolis = FPlane(f = 1e-4) # [s⁻¹]
buoyancy = SeawaterBuoyancy()
SeawaterBuoyancy{Float64}:
├── gravitational_acceleration: 9.80665
└── equation_of_state: LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078)

Specify parameters that are used to construct the background state.

background_state_parameters = (M = 1e-4,       # s⁻¹, geostrophic shear
                               f = coriolis.f, # s⁻¹, Coriolis parameter
                               N = 1e-4,       # s⁻¹, buoyancy frequency
                               H = grid.Lz,
                               g = buoyancy.gravitational_acceleration,
                               α = buoyancy.equation_of_state.thermal_expansion)
(M = 0.0001, f = 0.0001, N = 0.0001, H = 100.0, g = 9.80665, α = 0.000167)

We assume a background buoyancy $B$ with a constant stratification and also a constant lateral gradient (in the zonal direction). The background velocity components $U$ and $V$ are prescribed so that the thermal wind relationship is satisfied, that is, $f \partial_z U = - \partial_y B$ and $f \partial_z V = \partial_x B$.

T(x, y, z, t, p) = (p.M^2 * x + p.N^2 * (z + p.H)) / (p.g * p.α)
V(x, y, z, t, p) = p.M^2 / p.f * (z + p.H)

V_field = BackgroundField(V, parameters = background_state_parameters)
T_field = BackgroundField(T, parameters = background_state_parameters)
BackgroundField{typeof(Main.var"##277".T), @NamedTuple{M::Float64, f::Float64, N::Float64, H::Float64, g::Float64, α::Float64}}
├── func: T (generic function with 1 method)
└── parameters: (M = 0.0001, f = 0.0001, N = 0.0001, H = 100.0, g = 9.80665, α = 0.000167)

Specify some horizontal and vertical viscosity and diffusivity.

νᵥ = κᵥ = 1e-4 # [m² s⁻¹]
vertical_diffusivity = VerticalScalarDiffusivity(ν = νᵥ, κ = κᵥ)
VerticalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=0.0001)

Setup the biogeochemical model with optional carbonate chemistry turned on.

biogeochemistry = LOBSTER(grid;
                          carbonate_system = CarbonateSystem(),
                          detritus = TwoParticleAndDissolved(grid; open_bottom = true))

DIC_bcs = FieldBoundaryConditions(top = CarbonDioxideGasExchangeBoundaryCondition())
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── top: FluxBoundaryCondition: DiscreteBoundaryFunction with GasExchange
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)

Model instantiation

model = NonhydrostaticModel(grid;
                            biogeochemistry,
                            boundary_conditions = (DIC = DIC_bcs, ),
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            coriolis,
                            tracers = (:T, :S),
                            buoyancy,
                            background_fields = (T = T_field, v = V_field),
                            closure = vertical_diffusivity)
NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 32×32×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
├── tracers: (T, S, NO₃, NH₄, P, Z, sPOM, bPOM, DOM, DIC, Alk)
├── closure: VerticalScalarDiffusivity{ExplicitTimeDiscretization}(ν=0.0001, κ=(T=0.0001, S=0.0001, NO₃=0.0001, NH₄=0.0001, P=0.0001, Z=0.0001, sPOM=0.0001, bPOM=0.0001, DOM=0.0001, DIC=0.0001, Alk=0.0001))
├── buoyancy: SeawaterBuoyancy with g=9.80665 and LinearEquationOfState(thermal_expansion=0.000167, haline_contraction=0.00078) with ĝ = NegativeZDirection()
└── coriolis: FPlane{Float64}(f=0.0001)

Initial conditions

Start with a bit of random noise added to the background thermal wind and an arbitary biogeochemical state.

Ξ(z) = randn() * z / grid.Lz * (z / grid.Lz + 1)

Ũ = 1e-3
uᵢ(x, y, z) = Ũ * Ξ(z)
vᵢ(x, y, z) = Ũ * Ξ(z)

set!(model, u=uᵢ, v=vᵢ, P = 0.03, Z = 0.03, NO₃ = 4.0, NH₄ = 0.05, DIC = 2200.0, Alk = 2409.0, S = 35, T = 20)

Setup the simulation

Choose an appropriate initial timestep for this resolution and set up the simulation

Δx = minimum_xspacing(grid, Center(), Center(), Center())
Δy = minimum_yspacing(grid, Center(), Center(), Center())
Δz = minimum_zspacing(grid, Center(), Center(), Center())

Δt₀ = 0.75 * min(Δx, Δy, Δz) / V(0, 0, 0, 0, background_state_parameters)

simulation = Simulation(model, Δt = Δt₀, stop_time = 10days)
Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 15.625 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 10 days
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

Adapt the time step while keeping the CFL number fixed.

wizard = TimeStepWizard(cfl = 0.75, diffusive_cfl = 0.75, max_Δt = 30minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(5))

Create a progress message.

progress(sim) = @printf("i: % 6d, sim time: % 10s, wall time: % 10s, Δt: % 10s, CFL: %.2e\n",
                        sim.model.clock.iteration,
                        prettytime(sim.model.clock.time),
                        prettytime(sim.run_wall_time),
                        prettytime(sim.Δt),
                        AdvectiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(20))
Callback of progress on IterationInterval(20)

Here, we add some diagnostics to calculate and output.

u, v, w = model.velocities # unpack velocity `Field`s
NamedTuple with 3 Fields on 32×32×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo:
├── u: 32×32×8 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
├── v: 32×32×8 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
└── w: 32×32×9 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Face} on Oceananigans.Grids.RectilinearGrid on CPU

and also calculate the vertical vorticity.

ζ = Field(∂x(v) - ∂y(u))
32×32×8 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
├── grid: 32×32×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
├── operand: BinaryOperation at (Face, Face, Center)
├── status: time=0.0
└── data: 38×38×14 OffsetArray(::Array{Float64, 3}, -2:35, -2:35, -2:11) with eltype Float64 with indices -2:35×-2:35×-2:11
    └── max=5.35045e-5, min=-4.94356e-5, mean=2.83309e-23

Periodically save the velocities and vorticity to a file.

simulation.output_writers[:fields] = JLD2Writer(model, merge(model.tracers, (; u, v, w, ζ));
                                                schedule = TimeInterval(2hours),
                                                filename = "eady_turbulence_bgc",
                                                overwrite_existing = true)

Run the simulation

run!(simulation)
[ Info: Initializing simulation...
i:      0, sim time:  0 seconds, wall time:  0 seconds, Δt: 17.188 minutes, CFL: 3.27e-01
[ Info:     ... simulation initialization complete (10.379 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.054 seconds).
i:     20, sim time:    6 hours, wall time: 14.550 seconds, Δt: 25.164 minutes, CFL: 4.77e-01
i:     40, sim time:   15 hours, wall time: 15.482 seconds, Δt: 30 minutes, CFL: 5.76e-01
i:     60, sim time: 1.042 days, wall time: 16.401 seconds, Δt: 30 minutes, CFL: 6.17e-01
i:     80, sim time: 1.458 days, wall time: 17.323 seconds, Δt: 30 minutes, CFL: 6.40e-01
i:    100, sim time: 1.875 days, wall time: 18.243 seconds, Δt: 30 minutes, CFL: 6.66e-01
i:    120, sim time: 2.292 days, wall time: 19.160 seconds, Δt: 29.178 minutes, CFL: 7.50e-01
i:    140, sim time: 2.621 days, wall time: 20.068 seconds, Δt: 25.958 minutes, CFL: 7.50e-01
i:    160, sim time: 2.933 days, wall time: 20.969 seconds, Δt: 22.873 minutes, CFL: 7.50e-01
i:    180, sim time: 3.210 days, wall time: 21.862 seconds, Δt: 20.092 minutes, CFL: 7.50e-01
i:    200, sim time: 3.453 days, wall time: 22.753 seconds, Δt: 17.223 minutes, CFL: 7.50e-01
i:    220, sim time: 3.667 days, wall time: 23.665 seconds, Δt: 15.425 minutes, CFL: 7.50e-01
i:    240, sim time: 3.863 days, wall time: 24.580 seconds, Δt: 14.032 minutes, CFL: 7.50e-01
i:    260, sim time: 4.046 days, wall time: 25.417 seconds, Δt: 13.577 minutes, CFL: 7.50e-01
i:    280, sim time: 4.230 days, wall time: 26.106 seconds, Δt: 12.835 minutes, CFL: 7.50e-01
i:    300, sim time: 4.396 days, wall time: 26.812 seconds, Δt: 12.462 minutes, CFL: 7.50e-01
i:    320, sim time: 4.563 days, wall time: 27.506 seconds, Δt: 12.679 minutes, CFL: 7.50e-01
i:    340, sim time: 4.726 days, wall time: 28.198 seconds, Δt: 12.097 minutes, CFL: 7.50e-01
i:    360, sim time: 4.883 days, wall time: 28.890 seconds, Δt: 11.536 minutes, CFL: 7.50e-01
i:    380, sim time: 5.031 days, wall time: 29.581 seconds, Δt: 10.889 minutes, CFL: 7.50e-01
i:    400, sim time: 5.174 days, wall time: 30.270 seconds, Δt: 10.302 minutes, CFL: 7.50e-01
i:    420, sim time: 5.315 days, wall time: 30.958 seconds, Δt: 10.453 minutes, CFL: 7.50e-01
i:    440, sim time: 5.453 days, wall time: 31.652 seconds, Δt: 10.508 minutes, CFL: 7.50e-01
i:    460, sim time: 5.591 days, wall time: 32.341 seconds, Δt: 10.420 minutes, CFL: 7.50e-01
i:    480, sim time: 5.732 days, wall time: 33.026 seconds, Δt: 10.228 minutes, CFL: 7.50e-01
i:    500, sim time: 5.868 days, wall time: 33.724 seconds, Δt: 9.901 minutes, CFL: 7.50e-01
i:    520, sim time: 5.999 days, wall time: 34.429 seconds, Δt: 10.011 minutes, CFL: 7.50e-01
i:    540, sim time: 6.123 days, wall time: 35.134 seconds, Δt: 9.489 minutes, CFL: 7.50e-01
i:    560, sim time: 6.250 days, wall time: 35.823 seconds, Δt: 9.596 minutes, CFL: 7.50e-01
i:    580, sim time: 6.380 days, wall time: 36.516 seconds, Δt: 9.800 minutes, CFL: 7.50e-01
i:    600, sim time: 6.515 days, wall time: 37.211 seconds, Δt: 10.390 minutes, CFL: 7.50e-01
i:    620, sim time: 6.656 days, wall time: 37.901 seconds, Δt: 10.550 minutes, CFL: 7.50e-01
i:    640, sim time: 6.794 days, wall time: 38.593 seconds, Δt: 10.560 minutes, CFL: 7.50e-01
i:    660, sim time: 6.931 days, wall time: 39.284 seconds, Δt: 10.087 minutes, CFL: 7.50e-01
i:    680, sim time: 7.069 days, wall time: 39.971 seconds, Δt: 9.995 minutes, CFL: 7.50e-01
i:    700, sim time: 7.202 days, wall time: 40.665 seconds, Δt: 10.060 minutes, CFL: 7.50e-01
i:    720, sim time: 7.340 days, wall time: 41.361 seconds, Δt: 10.001 minutes, CFL: 7.50e-01
i:    740, sim time: 7.471 days, wall time: 42.049 seconds, Δt: 9.678 minutes, CFL: 7.50e-01
i:    760, sim time: 7.596 days, wall time: 42.742 seconds, Δt: 9.136 minutes, CFL: 7.50e-01
i:    780, sim time: 7.716 days, wall time: 43.433 seconds, Δt: 8.619 minutes, CFL: 7.50e-01
i:    800, sim time: 7.833 days, wall time: 44.123 seconds, Δt: 8.463 minutes, CFL: 7.50e-01
i:    820, sim time: 7.947 days, wall time: 44.822 seconds, Δt: 8.806 minutes, CFL: 7.50e-01
i:    840, sim time: 8.071 days, wall time: 45.516 seconds, Δt: 9.631 minutes, CFL: 7.50e-01
i:    860, sim time: 8.203 days, wall time: 46.216 seconds, Δt: 11.073 minutes, CFL: 7.50e-01
i:    880, sim time: 8.357 days, wall time: 46.919 seconds, Δt: 11.282 minutes, CFL: 7.50e-01
i:    900, sim time: 8.500 days, wall time: 47.616 seconds, Δt: 10.309 minutes, CFL: 7.50e-01
i:    920, sim time: 8.644 days, wall time: 48.322 seconds, Δt: 10.955 minutes, CFL: 7.50e-01
i:    940, sim time: 8.796 days, wall time: 49.026 seconds, Δt: 11.128 minutes, CFL: 7.50e-01
i:    960, sim time: 8.950 days, wall time: 49.737 seconds, Δt: 12.503 minutes, CFL: 7.50e-01
i:    980, sim time: 9.116 days, wall time: 50.442 seconds, Δt: 11.607 minutes, CFL: 7.50e-01
i:   1000, sim time: 9.265 days, wall time: 51.146 seconds, Δt: 10.902 minutes, CFL: 7.50e-01
i:   1020, sim time: 9.409 days, wall time: 51.841 seconds, Δt: 10.772 minutes, CFL: 7.50e-01
i:   1040, sim time: 9.554 days, wall time: 52.540 seconds, Δt: 11.117 minutes, CFL: 7.50e-01
i:   1060, sim time: 9.696 days, wall time: 53.238 seconds, Δt: 10.525 minutes, CFL: 7.50e-01
i:   1080, sim time: 9.833 days, wall time: 53.931 seconds, Δt: 10.327 minutes, CFL: 7.50e-01
i:   1100, sim time: 9.972 days, wall time: 54.627 seconds, Δt: 10.005 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 54.762 seconds.
[ Info: Simulation time 10 days equals or exceeds stop time 10 days.

Now load the saved output,

  ζ = FieldTimeSeries("eady_turbulence_bgc.jld2", "ζ")
  P = FieldTimeSeries("eady_turbulence_bgc.jld2", "P")
NO₃ = FieldTimeSeries("eady_turbulence_bgc.jld2", "NO₃")
NH₄ = FieldTimeSeries("eady_turbulence_bgc.jld2", "NH₄")
DIC = FieldTimeSeries("eady_turbulence_bgc.jld2", "DIC")

times = ζ.times

xζ, yζ, zζ = nodes(ζ)
xc, yc, zc = nodes(P)

and build the frames,

n = Observable(1)

Nz = grid.Nz

  ζₙ = @lift interior(  ζ[$n], :, :, Nz)
  Nₙ = @lift interior(NO₃[$n], :, :, Nz) .+ interior(NH₄[$n], :, :, Nz)
  Pₙ = @lift interior(  P[$n], :, :, Nz)
DICₙ = @lift interior(DIC[$n], :, :, Nz)
Observable([2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0; 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0 2200.0])

now plot

fig = Figure(size = (1600, 1600), fontsize = 20)

lims = [(minimum(T), maximum(T)) for T in (  ζ[:, :, Nz, :],
                                           NO₃[:, :, Nz, :] .+ NH₄[:, :, Nz, :],
                                             P[:, :, Nz, :],
                                           DIC[:, :, Nz, :])]

axis_kwargs = (xlabel = "x (m)", ylabel = "y (m)", aspect = DataAspect())

ax1 = Axis(fig[1, 1]; title = "Vertical vorticity (1 / s)", axis_kwargs...)
hm1 = heatmap!(ax1, xζ, yζ, ζₙ, colormap = :balance, colorrange = lims[1])
Colorbar(fig[1, 2], hm1)

ax2 = Axis(fig[1, 3]; title = "Nutrient (NO₃ + NH₄) concentration (mmol N / m³)", axis_kwargs...)
hm2 = heatmap!(ax2, xc, yc, Nₙ, colormap = Reverse(:bamako), colorrange = lims[2])
Colorbar(fig[1, 4], hm2)

ax3 = Axis(fig[2, 1]; title = "Phytoplankton concentration (mmol N / m³)", axis_kwargs...)
hm3 = heatmap!(ax3, xc, yc, Pₙ, colormap = Reverse(:batlow), colorrange = lims[3])
Colorbar(fig[2, 2], hm3)

ax4 = Axis(fig[2, 3]; title = "Dissolved inorganic carbon (mmol C / m³)", axis_kwargs...)
hm4 = heatmap!(ax4, xc, yc, DICₙ, colormap = Reverse(:devon), colorrange = lims[4])
Colorbar(fig[2, 4], hm4)

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

and record the movie record(fig, "eady.mp4", 1:length(times), framerate = 12) do i n[] = i end=# # temporatily disable n[] = length(times) fig #= nothing #hide

This page was generated using Literate.jl.