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.5Set 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: 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 entriesAdapt 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`sNamedTuple 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 CPUand 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-23Periodically 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: 11.458 minutes, CFL: 2.18e-01
[ Info: ... simulation initialization complete (11.291 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.264 seconds).
i: 20, sim time: 4.254 hours, wall time: 17.919 seconds, Δt: 16.776 minutes, CFL: 3.13e-01
i: 40, sim time: 10 hours, wall time: 19.978 seconds, Δt: 24.562 minutes, CFL: 4.74e-01
i: 60, sim time: 18 hours, wall time: 21.509 seconds, Δt: 25.639 minutes, CFL: 5.00e-01
i: 80, sim time: 1.083 days, wall time: 23.121 seconds, Δt: 24.175 minutes, CFL: 5.00e-01
i: 100, sim time: 1.382 days, wall time: 24.902 seconds, Δt: 23.331 minutes, CFL: 5.00e-01
i: 120, sim time: 1.663 days, wall time: 26.146 seconds, Δt: 22.532 minutes, CFL: 5.00e-01
i: 140, sim time: 1.932 days, wall time: 27.365 seconds, Δt: 22.296 minutes, CFL: 5.00e-01
i: 160, sim time: 2.210 days, wall time: 28.561 seconds, Δt: 20.084 minutes, CFL: 5.00e-01
i: 180, sim time: 2.456 days, wall time: 29.842 seconds, Δt: 18.449 minutes, CFL: 5.00e-01
i: 200, sim time: 2.691 days, wall time: 31.066 seconds, Δt: 16.846 minutes, CFL: 5.00e-01
i: 220, sim time: 2.900 days, wall time: 32.248 seconds, Δt: 15.405 minutes, CFL: 5.00e-01
i: 240, sim time: 3.093 days, wall time: 34.249 seconds, Δt: 13.978 minutes, CFL: 5.00e-01
i: 260, sim time: 3.278 days, wall time: 35.518 seconds, Δt: 12.579 minutes, CFL: 5.00e-01
i: 280, sim time: 3.441 days, wall time: 36.870 seconds, Δt: 11.698 minutes, CFL: 5.00e-01
i: 300, sim time: 3.591 days, wall time: 38.080 seconds, Δt: 10.342 minutes, CFL: 5.00e-01
i: 320, sim time: 3.729 days, wall time: 39.366 seconds, Δt: 9.946 minutes, CFL: 5.00e-01
i: 340, sim time: 3.853 days, wall time: 40.554 seconds, Δt: 9.638 minutes, CFL: 5.00e-01
i: 360, sim time: 3.981 days, wall time: 41.732 seconds, Δt: 9.199 minutes, CFL: 5.00e-01
i: 380, sim time: 4.102 days, wall time: 42.918 seconds, Δt: 9.142 minutes, CFL: 5.00e-01
i: 400, sim time: 4.221 days, wall time: 44.094 seconds, Δt: 8.593 minutes, CFL: 5.00e-01
i: 420, sim time: 4.339 days, wall time: 45.284 seconds, Δt: 8.493 minutes, CFL: 5.00e-01
i: 440, sim time: 4.453 days, wall time: 46.471 seconds, Δt: 8.461 minutes, CFL: 5.00e-01
i: 460, sim time: 4.566 days, wall time: 47.732 seconds, Δt: 8.448 minutes, CFL: 5.00e-01
i: 480, sim time: 4.678 days, wall time: 48.923 seconds, Δt: 8.205 minutes, CFL: 5.00e-01
i: 500, sim time: 4.788 days, wall time: 50.132 seconds, Δt: 7.975 minutes, CFL: 5.00e-01
i: 520, sim time: 4.893 days, wall time: 51.372 seconds, Δt: 7.869 minutes, CFL: 5.00e-01
i: 540, sim time: 4.996 days, wall time: 52.559 seconds, Δt: 7.480 minutes, CFL: 5.00e-01
i: 560, sim time: 5.093 days, wall time: 55.544 seconds, Δt: 7.055 minutes, CFL: 5.00e-01
i: 580, sim time: 5.186 days, wall time: 57.273 seconds, Δt: 6.827 minutes, CFL: 5.00e-01
i: 600, sim time: 5.279 days, wall time: 58.962 seconds, Δt: 6.949 minutes, CFL: 5.00e-01
i: 620, sim time: 5.372 days, wall time: 1.007 minutes, Δt: 6.977 minutes, CFL: 5.00e-01
i: 640, sim time: 5.464 days, wall time: 1.042 minutes, Δt: 6.959 minutes, CFL: 5.00e-01
i: 660, sim time: 5.559 days, wall time: 1.067 minutes, Δt: 7.048 minutes, CFL: 5.00e-01
i: 680, sim time: 5.651 days, wall time: 1.090 minutes, Δt: 6.862 minutes, CFL: 5.00e-01
i: 700, sim time: 5.743 days, wall time: 1.113 minutes, Δt: 6.870 minutes, CFL: 5.00e-01
i: 720, sim time: 5.833 days, wall time: 1.149 minutes, Δt: 6.630 minutes, CFL: 5.00e-01
i: 740, sim time: 5.921 days, wall time: 1.192 minutes, Δt: 6.560 minutes, CFL: 5.00e-01
i: 760, sim time: 6.009 days, wall time: 1.223 minutes, Δt: 6.635 minutes, CFL: 5.00e-01
i: 780, sim time: 6.097 days, wall time: 1.243 minutes, Δt: 6.345 minutes, CFL: 5.00e-01
i: 800, sim time: 6.184 days, wall time: 1.263 minutes, Δt: 6.336 minutes, CFL: 5.00e-01
i: 820, sim time: 6.272 days, wall time: 1.284 minutes, Δt: 6.347 minutes, CFL: 5.00e-01
i: 840, sim time: 6.360 days, wall time: 1.306 minutes, Δt: 6.468 minutes, CFL: 5.00e-01
i: 860, sim time: 6.449 days, wall time: 1.327 minutes, Δt: 6.848 minutes, CFL: 5.00e-01
i: 880, sim time: 6.543 days, wall time: 1.348 minutes, Δt: 6.927 minutes, CFL: 5.00e-01
i: 900, sim time: 6.637 days, wall time: 1.368 minutes, Δt: 7.016 minutes, CFL: 5.00e-01
i: 920, sim time: 6.730 days, wall time: 1.393 minutes, Δt: 7.061 minutes, CFL: 5.00e-01
i: 940, sim time: 6.824 days, wall time: 1.415 minutes, Δt: 7.020 minutes, CFL: 5.00e-01
i: 960, sim time: 6.915 days, wall time: 1.435 minutes, Δt: 6.795 minutes, CFL: 5.00e-01
i: 980, sim time: 7.005 days, wall time: 1.459 minutes, Δt: 6.637 minutes, CFL: 5.00e-01
i: 1000, sim time: 7.093 days, wall time: 1.480 minutes, Δt: 6.657 minutes, CFL: 5.00e-01
i: 1020, sim time: 7.185 days, wall time: 1.513 minutes, Δt: 6.713 minutes, CFL: 5.00e-01
i: 1040, sim time: 7.278 days, wall time: 1.534 minutes, Δt: 6.706 minutes, CFL: 5.00e-01
i: 1060, sim time: 7.370 days, wall time: 1.554 minutes, Δt: 6.645 minutes, CFL: 5.00e-01
i: 1080, sim time: 7.458 days, wall time: 1.575 minutes, Δt: 6.458 minutes, CFL: 5.00e-01
i: 1100, sim time: 7.544 days, wall time: 1.595 minutes, Δt: 6.189 minutes, CFL: 5.00e-01
i: 1120, sim time: 7.626 days, wall time: 1.616 minutes, Δt: 5.971 minutes, CFL: 5.00e-01
i: 1140, sim time: 7.707 days, wall time: 1.637 minutes, Δt: 5.739 minutes, CFL: 5.00e-01
i: 1160, sim time: 7.785 days, wall time: 1.658 minutes, Δt: 5.656 minutes, CFL: 5.00e-01
i: 1180, sim time: 7.861 days, wall time: 1.690 minutes, Δt: 5.705 minutes, CFL: 5.00e-01
i: 1200, sim time: 7.941 days, wall time: 1.710 minutes, Δt: 5.851 minutes, CFL: 5.00e-01
i: 1220, sim time: 8.021 days, wall time: 1.731 minutes, Δt: 6.171 minutes, CFL: 5.00e-01
i: 1240, sim time: 8.106 days, wall time: 1.752 minutes, Δt: 6.628 minutes, CFL: 5.00e-01
i: 1260, sim time: 8.201 days, wall time: 1.773 minutes, Δt: 7.359 minutes, CFL: 5.00e-01
i: 1280, sim time: 8.305 days, wall time: 1.792 minutes, Δt: 7.658 minutes, CFL: 5.00e-01
i: 1300, sim time: 8.407 days, wall time: 1.812 minutes, Δt: 7.229 minutes, CFL: 5.00e-01
i: 1320, sim time: 8.500 days, wall time: 1.834 minutes, Δt: 6.869 minutes, CFL: 5.00e-01
i: 1340, sim time: 8.593 days, wall time: 1.854 minutes, Δt: 7.372 minutes, CFL: 5.00e-01
i: 1360, sim time: 8.692 days, wall time: 1.874 minutes, Δt: 7.275 minutes, CFL: 5.00e-01
i: 1380, sim time: 8.791 days, wall time: 1.897 minutes, Δt: 7.414 minutes, CFL: 5.00e-01
i: 1400, sim time: 8.892 days, wall time: 1.918 minutes, Δt: 7.972 minutes, CFL: 5.00e-01
i: 1420, sim time: 9 days, wall time: 1.938 minutes, Δt: 8.279 minutes, CFL: 5.00e-01
i: 1440, sim time: 9.111 days, wall time: 1.958 minutes, Δt: 7.755 minutes, CFL: 5.00e-01
i: 1460, sim time: 9.213 days, wall time: 1.978 minutes, Δt: 7.267 minutes, CFL: 5.00e-01
i: 1480, sim time: 9.311 days, wall time: 1.997 minutes, Δt: 7.265 minutes, CFL: 5.00e-01
i: 1500, sim time: 9.408 days, wall time: 2.017 minutes, Δt: 7.175 minutes, CFL: 5.00e-01
i: 1520, sim time: 9.505 days, wall time: 2.037 minutes, Δt: 7.481 minutes, CFL: 5.00e-01
i: 1540, sim time: 9.604 days, wall time: 2.058 minutes, Δt: 7.280 minutes, CFL: 5.00e-01
i: 1560, sim time: 9.701 days, wall time: 2.078 minutes, Δt: 7.030 minutes, CFL: 5.00e-01
i: 1580, sim time: 9.794 days, wall time: 2.098 minutes, Δt: 6.976 minutes, CFL: 5.00e-01
i: 1600, sim time: 9.886 days, wall time: 2.119 minutes, Δt: 6.763 minutes, CFL: 5.00e-01
i: 1620, sim time: 9.977 days, wall time: 2.138 minutes, Δt: 6.661 minutes, CFL: 5.00e-01
[ Info: Simulation is stopping after running for 2.143 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)Makie.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.