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;
                          inorganic_carbon = CarbonateSystem(),
                          scale_negatives = 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(),
                            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:
│   ├── momentum: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── T: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── S: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── NO₃: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── NH₄: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── P: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── Z: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── DOM: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── sPOM: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── bPOM: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── DIC: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   └── Alk: WENO{3, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
├── tracers: (T, S, NO₃, NH₄, P, Z, DOM, sPOM, bPOM, 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, DOM=0.0001, sPOM=0.0001, bPOM=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.5 * 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: 10.417 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.5, diffusive_cfl = 0.5, 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.0773e-5, min=-5.41179e-5, mean=3.14329e-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_files = true)

Run the simulation

run!(simulation)
[ Info: Initializing simulation...
i:      0, sim time:  0 seconds, wall time:  0 seconds, Δt: 11.458 minutes, CFL: 2.17e-01
[ Info:     ... simulation initialization complete (10.645 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (2.627 seconds).
i:     20, sim time: 4.254 hours, wall time: 14.173 seconds, Δt: 16.776 minutes, CFL: 3.13e-01
i:     40, sim time:   10 hours, wall time: 15.124 seconds, Δt: 24.562 minutes, CFL: 4.71e-01
i:     60, sim time:   18 hours, wall time: 16.145 seconds, Δt: 25.233 minutes, CFL: 5.00e-01
i:     80, sim time: 1.083 days, wall time: 17.151 seconds, Δt: 24.535 minutes, CFL: 5.00e-01
i:    100, sim time: 1.399 days, wall time: 18.141 seconds, Δt: 23.404 minutes, CFL: 5.00e-01
i:    120, sim time: 1.667 days, wall time: 19.101 seconds, Δt: 22.918 minutes, CFL: 5.00e-01
i:    140, sim time: 1.946 days, wall time: 20.080 seconds, Δt: 21.037 minutes, CFL: 5.00e-01
i:    160, sim time: 2.223 days, wall time: 21.036 seconds, Δt: 20.209 minutes, CFL: 5.00e-01
i:    180, sim time: 2.484 days, wall time: 21.979 seconds, Δt: 19.010 minutes, CFL: 5.00e-01
i:    200, sim time: 2.715 days, wall time: 22.931 seconds, Δt: 16.973 minutes, CFL: 5.00e-01
i:    220, sim time: 2.928 days, wall time: 23.887 seconds, Δt: 15.309 minutes, CFL: 5.00e-01
i:    240, sim time: 3.123 days, wall time: 24.820 seconds, Δt: 13.833 minutes, CFL: 5.00e-01
i:    260, sim time: 3.294 days, wall time: 25.745 seconds, Δt: 12.165 minutes, CFL: 5.00e-01
i:    280, sim time: 3.448 days, wall time: 26.679 seconds, Δt: 11.360 minutes, CFL: 5.00e-01
i:    300, sim time: 3.591 days, wall time: 27.604 seconds, Δt: 10.569 minutes, CFL: 5.00e-01
i:    320, sim time: 3.729 days, wall time: 28.501 seconds, Δt: 9.827 minutes, CFL: 5.00e-01
i:    340, sim time: 3.854 days, wall time: 29.423 seconds, Δt: 9.643 minutes, CFL: 5.00e-01
i:    360, sim time: 3.980 days, wall time: 30.331 seconds, Δt: 8.872 minutes, CFL: 5.00e-01
i:    380, sim time: 4.095 days, wall time: 31.257 seconds, Δt: 8.674 minutes, CFL: 5.00e-01
i:    400, sim time: 4.208 days, wall time: 32.163 seconds, Δt: 8.609 minutes, CFL: 5.00e-01
i:    420, sim time: 4.321 days, wall time: 33.069 seconds, Δt: 8.286 minutes, CFL: 5.00e-01
i:    440, sim time: 4.428 days, wall time: 33.999 seconds, Δt: 8.332 minutes, CFL: 5.00e-01
i:    460, sim time: 4.541 days, wall time: 34.915 seconds, Δt: 8.124 minutes, CFL: 5.00e-01
i:    480, sim time: 4.652 days, wall time: 35.814 seconds, Δt: 8.099 minutes, CFL: 5.00e-01
i:    500, sim time: 4.756 days, wall time: 36.743 seconds, Δt: 7.761 minutes, CFL: 5.00e-01
i:    520, sim time: 4.860 days, wall time: 37.647 seconds, Δt: 7.770 minutes, CFL: 5.00e-01
i:    540, sim time: 4.966 days, wall time: 38.556 seconds, Δt: 8.063 minutes, CFL: 5.00e-01
i:    560, sim time: 5.073 days, wall time: 39.466 seconds, Δt: 8.078 minutes, CFL: 5.00e-01
i:    580, sim time: 5.184 days, wall time: 40.388 seconds, Δt: 8.240 minutes, CFL: 5.00e-01
i:    600, sim time: 5.296 days, wall time: 41.298 seconds, Δt: 8.196 minutes, CFL: 5.00e-01
i:    620, sim time: 5.408 days, wall time: 42.197 seconds, Δt: 8.391 minutes, CFL: 5.00e-01
i:    640, sim time: 5.518 days, wall time: 43.126 seconds, Δt: 8.461 minutes, CFL: 5.00e-01
i:    660, sim time: 5.627 days, wall time: 44.026 seconds, Δt: 7.774 minutes, CFL: 5.00e-01
i:    680, sim time: 5.729 days, wall time: 44.934 seconds, Δt: 7.343 minutes, CFL: 5.00e-01
i:    700, sim time: 5.826 days, wall time: 45.829 seconds, Δt: 7.449 minutes, CFL: 5.00e-01
i:    720, sim time: 5.927 days, wall time: 46.753 seconds, Δt: 7.576 minutes, CFL: 5.00e-01
i:    740, sim time: 6.032 days, wall time: 47.655 seconds, Δt: 7.829 minutes, CFL: 5.00e-01
i:    760, sim time: 6.139 days, wall time: 48.554 seconds, Δt: 8.203 minutes, CFL: 5.00e-01
i:    780, sim time: 6.250 days, wall time: 49.456 seconds, Δt: 8.095 minutes, CFL: 5.00e-01
i:    800, sim time: 6.354 days, wall time: 50.372 seconds, Δt: 7.490 minutes, CFL: 5.00e-01
i:    820, sim time: 6.459 days, wall time: 51.274 seconds, Δt: 7.616 minutes, CFL: 5.00e-01
i:    840, sim time: 6.563 days, wall time: 52.171 seconds, Δt: 7.627 minutes, CFL: 5.00e-01
i:    860, sim time: 6.665 days, wall time: 53.074 seconds, Δt: 7.144 minutes, CFL: 5.00e-01
i:    880, sim time: 6.760 days, wall time: 53.999 seconds, Δt: 6.607 minutes, CFL: 5.00e-01
i:    900, sim time: 6.846 days, wall time: 54.897 seconds, Δt: 6.014 minutes, CFL: 5.00e-01
i:    920, sim time: 6.925 days, wall time: 55.801 seconds, Δt: 5.518 minutes, CFL: 5.00e-01
i:    940, sim time:     7 days, wall time: 56.681 seconds, Δt: 5.334 minutes, CFL: 5.00e-01
i:    960, sim time: 7.075 days, wall time: 57.577 seconds, Δt: 5.280 minutes, CFL: 5.00e-01
i:    980, sim time: 7.144 days, wall time: 58.477 seconds, Δt: 5.029 minutes, CFL: 5.00e-01
i:   1000, sim time: 7.211 days, wall time: 59.367 seconds, Δt: 4.760 minutes, CFL: 5.00e-01
i:   1020, sim time: 7.272 days, wall time: 1.004 minutes, Δt: 4.550 minutes, CFL: 5.00e-01
i:   1040, sim time: 7.333 days, wall time: 1.019 minutes, Δt: 4.339 minutes, CFL: 5.00e-01
i:   1060, sim time: 7.394 days, wall time: 1.034 minutes, Δt: 4.482 minutes, CFL: 5.00e-01
i:   1080, sim time: 7.455 days, wall time: 1.049 minutes, Δt: 4.697 minutes, CFL: 5.00e-01
i:   1100, sim time: 7.520 days, wall time: 1.064 minutes, Δt: 4.790 minutes, CFL: 5.00e-01
i:   1120, sim time: 7.583 days, wall time: 1.079 minutes, Δt: 5.055 minutes, CFL: 5.00e-01
i:   1140, sim time: 7.655 days, wall time: 1.094 minutes, Δt: 5.270 minutes, CFL: 5.00e-01
i:   1160, sim time: 7.727 days, wall time: 1.109 minutes, Δt: 5.680 minutes, CFL: 5.00e-01
i:   1180, sim time: 7.807 days, wall time: 1.124 minutes, Δt: 5.780 minutes, CFL: 5.00e-01
i:   1200, sim time: 7.887 days, wall time: 1.139 minutes, Δt: 6.219 minutes, CFL: 5.00e-01
i:   1220, sim time: 7.976 days, wall time: 1.154 minutes, Δt: 7.076 minutes, CFL: 5.00e-01
i:   1240, sim time: 8.076 days, wall time: 1.169 minutes, Δt: 7.221 minutes, CFL: 5.00e-01
i:   1260, sim time: 8.172 days, wall time: 1.185 minutes, Δt: 6.830 minutes, CFL: 5.00e-01
i:   1280, sim time: 8.264 days, wall time: 1.200 minutes, Δt: 6.558 minutes, CFL: 5.00e-01
i:   1300, sim time: 8.352 days, wall time: 1.215 minutes, Δt: 6.625 minutes, CFL: 5.00e-01
i:   1320, sim time: 8.446 days, wall time: 1.230 minutes, Δt: 7.241 minutes, CFL: 5.00e-01
i:   1340, sim time: 8.544 days, wall time: 1.245 minutes, Δt: 6.954 minutes, CFL: 5.00e-01
i:   1360, sim time: 8.634 days, wall time: 1.260 minutes, Δt: 6.605 minutes, CFL: 5.00e-01
i:   1380, sim time: 8.721 days, wall time: 1.275 minutes, Δt: 6.351 minutes, CFL: 5.00e-01
i:   1400, sim time: 8.807 days, wall time: 1.290 minutes, Δt: 6.445 minutes, CFL: 5.00e-01
i:   1420, sim time: 8.899 days, wall time: 1.305 minutes, Δt: 6.628 minutes, CFL: 5.00e-01
i:   1440, sim time: 8.988 days, wall time: 1.320 minutes, Δt: 6.231 minutes, CFL: 5.00e-01
i:   1460, sim time: 9.073 days, wall time: 1.335 minutes, Δt: 6.251 minutes, CFL: 5.00e-01
i:   1480, sim time: 9.157 days, wall time: 1.350 minutes, Δt: 6.276 minutes, CFL: 5.00e-01
i:   1500, sim time: 9.239 days, wall time: 1.365 minutes, Δt: 6.030 minutes, CFL: 5.00e-01
i:   1520, sim time: 9.323 days, wall time: 1.380 minutes, Δt: 6.307 minutes, CFL: 5.00e-01
i:   1540, sim time: 9.406 days, wall time: 1.395 minutes, Δt: 5.982 minutes, CFL: 5.00e-01
i:   1560, sim time: 9.486 days, wall time: 1.410 minutes, Δt: 5.691 minutes, CFL: 5.00e-01
i:   1580, sim time: 9.562 days, wall time: 1.425 minutes, Δt: 5.631 minutes, CFL: 5.00e-01
i:   1600, sim time: 9.638 days, wall time: 1.440 minutes, Δt: 5.745 minutes, CFL: 5.00e-01
i:   1620, sim time: 9.714 days, wall time: 1.454 minutes, Δt: 5.801 minutes, CFL: 5.00e-01
i:   1640, sim time: 9.796 days, wall time: 1.470 minutes, Δt: 6.116 minutes, CFL: 5.00e-01
i:   1660, sim time: 9.882 days, wall time: 1.485 minutes, Δt: 6.487 minutes, CFL: 5.00e-01
i:   1680, sim time: 9.972 days, wall time: 1.500 minutes, Δt: 6.543 minutes, CFL: 5.00e-01
[ Info: Simulation is stopping after running for 1.505 minutes.
[ 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)
Label()

and record the movie

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


This page was generated using Literate.jl.