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.70902e-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 (10.408 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (2.866 seconds).
i: 20, sim time: 4.254 hours, wall time: 14.619 seconds, Δt: 16.776 minutes, CFL: 3.13e-01
i: 40, sim time: 10 hours, wall time: 15.817 seconds, Δt: 24.562 minutes, CFL: 4.74e-01
i: 60, sim time: 18 hours, wall time: 17.068 seconds, Δt: 25.639 minutes, CFL: 5.00e-01
i: 80, sim time: 1.083 days, wall time: 18.301 seconds, Δt: 24.175 minutes, CFL: 5.00e-01
i: 100, sim time: 1.382 days, wall time: 19.529 seconds, Δt: 23.331 minutes, CFL: 5.00e-01
i: 120, sim time: 1.663 days, wall time: 20.726 seconds, Δt: 22.532 minutes, CFL: 5.00e-01
i: 140, sim time: 1.932 days, wall time: 21.940 seconds, Δt: 22.296 minutes, CFL: 5.00e-01
i: 160, sim time: 2.210 days, wall time: 23.135 seconds, Δt: 20.084 minutes, CFL: 5.00e-01
i: 180, sim time: 2.456 days, wall time: 24.312 seconds, Δt: 18.449 minutes, CFL: 5.00e-01
i: 200, sim time: 2.691 days, wall time: 25.494 seconds, Δt: 16.846 minutes, CFL: 5.00e-01
i: 220, sim time: 2.900 days, wall time: 26.653 seconds, Δt: 15.405 minutes, CFL: 5.00e-01
i: 240, sim time: 3.093 days, wall time: 27.830 seconds, Δt: 13.978 minutes, CFL: 5.00e-01
i: 260, sim time: 3.278 days, wall time: 28.990 seconds, Δt: 12.579 minutes, CFL: 5.00e-01
i: 280, sim time: 3.441 days, wall time: 30.144 seconds, Δt: 11.698 minutes, CFL: 5.00e-01
i: 300, sim time: 3.591 days, wall time: 31.318 seconds, Δt: 10.342 minutes, CFL: 5.00e-01
i: 320, sim time: 3.729 days, wall time: 32.458 seconds, Δt: 9.946 minutes, CFL: 5.00e-01
i: 340, sim time: 3.853 days, wall time: 33.613 seconds, Δt: 9.638 minutes, CFL: 5.00e-01
i: 360, sim time: 3.981 days, wall time: 34.752 seconds, Δt: 9.199 minutes, CFL: 5.00e-01
i: 380, sim time: 4.102 days, wall time: 35.923 seconds, Δt: 9.142 minutes, CFL: 5.00e-01
i: 400, sim time: 4.221 days, wall time: 37.063 seconds, Δt: 8.593 minutes, CFL: 5.00e-01
i: 420, sim time: 4.339 days, wall time: 38.213 seconds, Δt: 8.493 minutes, CFL: 5.00e-01
i: 440, sim time: 4.453 days, wall time: 39.349 seconds, Δt: 8.461 minutes, CFL: 5.00e-01
i: 460, sim time: 4.566 days, wall time: 40.498 seconds, Δt: 8.448 minutes, CFL: 5.00e-01
i: 480, sim time: 4.678 days, wall time: 41.659 seconds, Δt: 8.205 minutes, CFL: 5.00e-01
i: 500, sim time: 4.788 days, wall time: 42.799 seconds, Δt: 7.975 minutes, CFL: 5.00e-01
i: 520, sim time: 4.893 days, wall time: 43.931 seconds, Δt: 7.869 minutes, CFL: 5.00e-01
i: 540, sim time: 4.996 days, wall time: 45.089 seconds, Δt: 7.480 minutes, CFL: 5.00e-01
i: 560, sim time: 5.093 days, wall time: 46.253 seconds, Δt: 7.055 minutes, CFL: 5.00e-01
i: 580, sim time: 5.186 days, wall time: 47.384 seconds, Δt: 6.827 minutes, CFL: 5.00e-01
i: 600, sim time: 5.279 days, wall time: 48.514 seconds, Δt: 6.949 minutes, CFL: 5.00e-01
i: 620, sim time: 5.372 days, wall time: 49.658 seconds, Δt: 6.977 minutes, CFL: 5.00e-01
i: 640, sim time: 5.464 days, wall time: 50.793 seconds, Δt: 6.959 minutes, CFL: 5.00e-01
i: 660, sim time: 5.559 days, wall time: 51.928 seconds, Δt: 7.048 minutes, CFL: 5.00e-01
i: 680, sim time: 5.651 days, wall time: 53.058 seconds, Δt: 6.862 minutes, CFL: 5.00e-01
i: 700, sim time: 5.743 days, wall time: 54.169 seconds, Δt: 6.870 minutes, CFL: 5.00e-01
i: 720, sim time: 5.833 days, wall time: 55.314 seconds, Δt: 6.630 minutes, CFL: 5.00e-01
i: 740, sim time: 5.921 days, wall time: 56.395 seconds, Δt: 6.560 minutes, CFL: 5.00e-01
i: 760, sim time: 6.009 days, wall time: 57.529 seconds, Δt: 6.635 minutes, CFL: 5.00e-01
i: 780, sim time: 6.097 days, wall time: 58.665 seconds, Δt: 6.345 minutes, CFL: 5.00e-01
i: 800, sim time: 6.184 days, wall time: 59.808 seconds, Δt: 6.336 minutes, CFL: 5.00e-01
i: 820, sim time: 6.272 days, wall time: 1.016 minutes, Δt: 6.347 minutes, CFL: 5.00e-01
i: 840, sim time: 6.360 days, wall time: 1.035 minutes, Δt: 6.468 minutes, CFL: 5.00e-01
i: 860, sim time: 6.449 days, wall time: 1.053 minutes, Δt: 6.848 minutes, CFL: 5.00e-01
i: 880, sim time: 6.543 days, wall time: 1.072 minutes, Δt: 6.927 minutes, CFL: 5.00e-01
i: 900, sim time: 6.637 days, wall time: 1.091 minutes, Δt: 7.016 minutes, CFL: 5.00e-01
i: 920, sim time: 6.730 days, wall time: 1.110 minutes, Δt: 7.061 minutes, CFL: 5.00e-01
i: 940, sim time: 6.824 days, wall time: 1.129 minutes, Δt: 7.020 minutes, CFL: 5.00e-01
i: 960, sim time: 6.915 days, wall time: 1.148 minutes, Δt: 6.795 minutes, CFL: 5.00e-01
i: 980, sim time: 7.005 days, wall time: 1.167 minutes, Δt: 6.637 minutes, CFL: 5.00e-01
i: 1000, sim time: 7.093 days, wall time: 1.186 minutes, Δt: 6.657 minutes, CFL: 5.00e-01
i: 1020, sim time: 7.185 days, wall time: 1.205 minutes, Δt: 6.713 minutes, CFL: 5.00e-01
i: 1040, sim time: 7.278 days, wall time: 1.223 minutes, Δt: 6.706 minutes, CFL: 5.00e-01
i: 1060, sim time: 7.370 days, wall time: 1.242 minutes, Δt: 6.645 minutes, CFL: 5.00e-01
i: 1080, sim time: 7.458 days, wall time: 1.261 minutes, Δt: 6.458 minutes, CFL: 5.00e-01
i: 1100, sim time: 7.544 days, wall time: 1.280 minutes, Δt: 6.189 minutes, CFL: 5.00e-01
i: 1120, sim time: 7.626 days, wall time: 1.299 minutes, Δt: 5.971 minutes, CFL: 5.00e-01
i: 1140, sim time: 7.707 days, wall time: 1.317 minutes, Δt: 5.739 minutes, CFL: 5.00e-01
i: 1160, sim time: 7.785 days, wall time: 1.336 minutes, Δt: 5.656 minutes, CFL: 5.00e-01
i: 1180, sim time: 7.861 days, wall time: 1.355 minutes, Δt: 5.705 minutes, CFL: 5.00e-01
i: 1200, sim time: 7.941 days, wall time: 1.374 minutes, Δt: 5.851 minutes, CFL: 5.00e-01
i: 1220, sim time: 8.021 days, wall time: 1.392 minutes, Δt: 6.171 minutes, CFL: 5.00e-01
i: 1240, sim time: 8.106 days, wall time: 1.411 minutes, Δt: 6.627 minutes, CFL: 5.00e-01
i: 1260, sim time: 8.201 days, wall time: 1.430 minutes, Δt: 7.359 minutes, CFL: 5.00e-01
i: 1280, sim time: 8.305 days, wall time: 1.449 minutes, Δt: 7.658 minutes, CFL: 5.00e-01
i: 1300, sim time: 8.407 days, wall time: 1.468 minutes, Δt: 7.229 minutes, CFL: 5.00e-01
i: 1320, sim time: 8.500 days, wall time: 1.487 minutes, Δt: 6.869 minutes, CFL: 5.00e-01
i: 1340, sim time: 8.593 days, wall time: 1.506 minutes, Δt: 7.372 minutes, CFL: 5.00e-01
i: 1360, sim time: 8.692 days, wall time: 1.525 minutes, Δt: 7.275 minutes, CFL: 5.00e-01
i: 1380, sim time: 8.791 days, wall time: 1.544 minutes, Δt: 7.414 minutes, CFL: 5.00e-01
i: 1400, sim time: 8.892 days, wall time: 1.563 minutes, Δt: 7.972 minutes, CFL: 5.00e-01
i: 1420, sim time: 9 days, wall time: 1.582 minutes, Δt: 8.279 minutes, CFL: 5.00e-01
i: 1440, sim time: 9.111 days, wall time: 1.601 minutes, Δt: 7.755 minutes, CFL: 5.00e-01
i: 1460, sim time: 9.213 days, wall time: 1.620 minutes, Δt: 7.267 minutes, CFL: 5.00e-01
i: 1480, sim time: 9.311 days, wall time: 1.639 minutes, Δt: 7.265 minutes, CFL: 5.00e-01
i: 1500, sim time: 9.408 days, wall time: 1.658 minutes, Δt: 7.175 minutes, CFL: 5.00e-01
i: 1520, sim time: 9.505 days, wall time: 1.677 minutes, Δt: 7.481 minutes, CFL: 5.00e-01
i: 1540, sim time: 9.604 days, wall time: 1.696 minutes, Δt: 7.280 minutes, CFL: 5.00e-01
i: 1560, sim time: 9.701 days, wall time: 1.715 minutes, Δt: 7.030 minutes, CFL: 5.00e-01
i: 1580, sim time: 9.794 days, wall time: 1.734 minutes, Δt: 6.976 minutes, CFL: 5.00e-01
i: 1600, sim time: 9.886 days, wall time: 1.752 minutes, Δt: 6.763 minutes, CFL: 5.00e-01
i: 1620, sim time: 9.977 days, wall time: 1.771 minutes, Δt: 6.661 minutes, CFL: 5.00e-01
[ Info: Simulation is stopping after running for 1.776 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.