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
using Oceananigans.Units
using Random
Random.seed!(11)Random.TaskLocalRNG()Construct a grid with uniform grid spacing on GPU (just remove GPU() to run on CPU)
grid = RectilinearGrid(GPU(); size = (32, 32, 8), extent = (1kilometer, 1kilometer, 100meters))32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU 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"##225".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,
carbonates = true,
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 (::OceanBioME.Models.GasExchangeModel.GasExchange{Float64, OceanBioME.Models.GasExchangeModel.ScaledGasTransferVelocity.SchmidtScaledTransferVelocity{OceanBioME.Models.GasExchangeModel.PolynomialParameterisation{2, Tuple{Float64, Float64, Float64}}, OceanBioME.Models.GasExchangeModel.PolynomialParameterisation{4, NTuple{5, Float64}}}, OceanBioME.Models.GasExchangeModel.CarbonDioxideConcentration{OceanBioME.Models.CarbonChemistryModel.CarbonChemistry{OceanBioME.Models.CarbonChemistryModel.K0{Float64}, @NamedTuple{K1::OceanBioME.Models.CarbonChemistryModel.K1{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, K2::OceanBioME.Models.CarbonChemistryModel.K2{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}}, OceanBioME.Models.CarbonChemistryModel.KB{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.KS{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.KF{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, OceanBioME.Models.CarbonChemistryModel.KS{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, @NamedTuple{KP1::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, KP2::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, KP3::OceanBioME.Models.CarbonChemistryModel.KP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}}, OceanBioME.Models.CarbonChemistryModel.KSi{OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, Float64}, OceanBioME.Models.CarbonChemistryModel.KW{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, OceanBioME.Models.CarbonChemistryModel.IonicStrength{Float64}, OceanBioME.Models.CarbonChemistryModel.KSP{Float64, OceanBioME.Models.CarbonChemistryModel.PressureCorrection{Float64}}, typeof(OceanBioME.Models.teos10_polynomial_approximation)}, OceanBioME.Models.GasExchangeModel.PolynomialVirialCoefficientForCarbonDioxide{Float64}, OceanBioME.Models.GasExchangeModel.CrossVirialCoefficientForCarbonDioxide{Float64}, Float64, Nothing}, Float64})
└── 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{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 halo
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: WENO{3, Float64, Float32}(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{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 15.625 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per 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 => 4
│ ├── stop_iteration_exceeded => -
│ ├── wall_time_limit_exceeded => e
│ └── nan_checker => }
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entriesAdapt 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`sNamedTuple with 3 Fields on 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 halo:
├── u: 32×32×8 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── v: 32×32×8 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Face, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CUDAGPU
└── w: 32×32×9 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Face} on Oceananigans.Grids.RectilinearGrid on CUDAGPUand 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 CUDAGPU
├── grid: 32×32×8 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU 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(::CUDA.CuArray{Float64, 3, CUDA.DeviceMemory}, -2:35, -2:35, -2:11) with eltype Float64 with indices -2:35×-2:35×-2:11
└── max=6.41337e-5, min=-5.24045e-5, mean=-1.28213e-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: 17.188 minutes, CFL: 3.25e-01
[ Info: ... simulation initialization complete (2.693 minutes)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (17.084 seconds).
i: 20, sim time: 6 hours, wall time: 2.990 minutes, Δt: 25.164 minutes, CFL: 4.70e-01
i: 40, sim time: 15 hours, wall time: 2.993 minutes, Δt: 30 minutes, CFL: 5.87e-01
i: 60, sim time: 1.042 days, wall time: 2.997 minutes, Δt: 30 minutes, CFL: 6.06e-01
i: 80, sim time: 1.458 days, wall time: 3.000 minutes, Δt: 30 minutes, CFL: 6.34e-01
i: 100, sim time: 1.875 days, wall time: 3.003 minutes, Δt: 30 minutes, CFL: 6.52e-01
i: 120, sim time: 2.292 days, wall time: 3.006 minutes, Δt: 30 minutes, CFL: 7.16e-01
i: 140, sim time: 2.687 days, wall time: 3.009 minutes, Δt: 28.625 minutes, CFL: 7.50e-01
i: 160, sim time: 3 days, wall time: 3.012 minutes, Δt: 25.231 minutes, CFL: 7.50e-01
i: 180, sim time: 3.308 days, wall time: 3.015 minutes, Δt: 19.652 minutes, CFL: 7.50e-01
i: 200, sim time: 3.548 days, wall time: 3.018 minutes, Δt: 16.176 minutes, CFL: 7.50e-01
i: 220, sim time: 3.749 days, wall time: 3.021 minutes, Δt: 14.433 minutes, CFL: 7.50e-01
i: 240, sim time: 3.926 days, wall time: 3.024 minutes, Δt: 13.164 minutes, CFL: 7.50e-01
i: 260, sim time: 4.092 days, wall time: 3.026 minutes, Δt: 13.173 minutes, CFL: 7.50e-01
i: 280, sim time: 4.250 days, wall time: 3.029 minutes, Δt: 12.970 minutes, CFL: 7.50e-01
i: 300, sim time: 4.413 days, wall time: 3.032 minutes, Δt: 12.543 minutes, CFL: 7.50e-01
i: 320, sim time: 4.577 days, wall time: 3.035 minutes, Δt: 12.198 minutes, CFL: 7.50e-01
i: 340, sim time: 4.743 days, wall time: 3.038 minutes, Δt: 11.900 minutes, CFL: 7.50e-01
i: 360, sim time: 4.909 days, wall time: 3.040 minutes, Δt: 11.816 minutes, CFL: 7.50e-01
i: 380, sim time: 5.066 days, wall time: 3.043 minutes, Δt: 11.594 minutes, CFL: 7.50e-01
i: 400, sim time: 5.216 days, wall time: 3.046 minutes, Δt: 11.887 minutes, CFL: 7.50e-01
i: 420, sim time: 5.378 days, wall time: 3.049 minutes, Δt: 12.793 minutes, CFL: 7.50e-01
i: 440, sim time: 5.544 days, wall time: 3.052 minutes, Δt: 12.223 minutes, CFL: 7.50e-01
i: 460, sim time: 5.698 days, wall time: 3.054 minutes, Δt: 11.115 minutes, CFL: 7.50e-01
i: 480, sim time: 5.841 days, wall time: 3.057 minutes, Δt: 10.360 minutes, CFL: 7.50e-01
i: 500, sim time: 5.980 days, wall time: 3.060 minutes, Δt: 9.878 minutes, CFL: 7.50e-01
i: 520, sim time: 6.109 days, wall time: 3.062 minutes, Δt: 8.825 minutes, CFL: 7.50e-01
i: 540, sim time: 6.226 days, wall time: 3.065 minutes, Δt: 8.569 minutes, CFL: 7.50e-01
i: 560, sim time: 6.333 days, wall time: 3.068 minutes, Δt: 8.053 minutes, CFL: 7.50e-01
i: 580, sim time: 6.438 days, wall time: 3.070 minutes, Δt: 7.533 minutes, CFL: 7.50e-01
i: 600, sim time: 6.535 days, wall time: 3.073 minutes, Δt: 7.065 minutes, CFL: 7.50e-01
i: 620, sim time: 6.630 days, wall time: 3.075 minutes, Δt: 6.587 minutes, CFL: 7.50e-01
i: 640, sim time: 6.722 days, wall time: 3.078 minutes, Δt: 6.692 minutes, CFL: 7.50e-01
i: 660, sim time: 6.815 days, wall time: 3.081 minutes, Δt: 6.533 minutes, CFL: 7.50e-01
i: 680, sim time: 6.902 days, wall time: 3.083 minutes, Δt: 6.739 minutes, CFL: 7.50e-01
i: 700, sim time: 6.995 days, wall time: 3.086 minutes, Δt: 7.441 minutes, CFL: 7.50e-01
i: 720, sim time: 7.100 days, wall time: 3.089 minutes, Δt: 8.095 minutes, CFL: 7.50e-01
i: 740, sim time: 7.204 days, wall time: 3.091 minutes, Δt: 7.773 minutes, CFL: 7.50e-01
i: 760, sim time: 7.311 days, wall time: 3.094 minutes, Δt: 7.987 minutes, CFL: 7.50e-01
i: 780, sim time: 7.422 days, wall time: 3.097 minutes, Δt: 8.171 minutes, CFL: 7.50e-01
i: 800, sim time: 7.537 days, wall time: 3.100 minutes, Δt: 9.040 minutes, CFL: 7.50e-01
i: 820, sim time: 7.659 days, wall time: 3.102 minutes, Δt: 8.828 minutes, CFL: 7.50e-01
i: 840, sim time: 7.775 days, wall time: 3.105 minutes, Δt: 8.821 minutes, CFL: 7.50e-01
i: 860, sim time: 7.896 days, wall time: 3.108 minutes, Δt: 9.231 minutes, CFL: 7.50e-01
i: 880, sim time: 8.019 days, wall time: 3.110 minutes, Δt: 8.945 minutes, CFL: 7.50e-01
i: 900, sim time: 8.136 days, wall time: 3.113 minutes, Δt: 8.205 minutes, CFL: 7.50e-01
i: 920, sim time: 8.244 days, wall time: 3.116 minutes, Δt: 7.660 minutes, CFL: 7.50e-01
i: 940, sim time: 8.344 days, wall time: 3.119 minutes, Δt: 7.684 minutes, CFL: 7.50e-01
i: 960, sim time: 8.449 days, wall time: 3.122 minutes, Δt: 7.679 minutes, CFL: 7.50e-01
i: 980, sim time: 8.554 days, wall time: 3.124 minutes, Δt: 7.947 minutes, CFL: 7.50e-01
i: 1000, sim time: 8.660 days, wall time: 3.127 minutes, Δt: 8.957 minutes, CFL: 7.50e-01
i: 1020, sim time: 8.776 days, wall time: 3.131 minutes, Δt: 9.564 minutes, CFL: 7.50e-01
i: 1040, sim time: 8.909 days, wall time: 3.133 minutes, Δt: 9.797 minutes, CFL: 7.50e-01
i: 1060, sim time: 9.035 days, wall time: 3.136 minutes, Δt: 9.899 minutes, CFL: 7.50e-01
i: 1080, sim time: 9.162 days, wall time: 3.139 minutes, Δt: 9.016 minutes, CFL: 7.50e-01
i: 1100, sim time: 9.280 days, wall time: 3.141 minutes, Δt: 8.442 minutes, CFL: 7.50e-01
i: 1120, sim time: 9.391 days, wall time: 3.144 minutes, Δt: 8.132 minutes, CFL: 7.50e-01
i: 1140, sim time: 9.500 days, wall time: 3.146 minutes, Δt: 8.251 minutes, CFL: 7.50e-01
i: 1160, sim time: 9.613 days, wall time: 3.149 minutes, Δt: 8.571 minutes, CFL: 7.50e-01
i: 1180, sim time: 9.733 days, wall time: 3.152 minutes, Δt: 8.967 minutes, CFL: 7.50e-01
i: 1200, sim time: 9.853 days, wall time: 3.155 minutes, Δt: 9.340 minutes, CFL: 7.50e-01
i: 1220, sim time: 9.982 days, wall time: 3.157 minutes, Δt: 9.456 minutes, CFL: 7.50e-01
[ Info: Simulation is stopping after running for 3.158 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 plot.
using CairoMakie
n = Observable(1)
ζₙ = @lift interior( ζ[$n], :, :, grid.Nz)
Nₙ = @lift interior(NO₃[$n], :, :, grid.Nz) .+ interior(NH₄[$n], :, :, grid.Nz)
Pₙ = @lift interior( P[$n], :, :, grid.Nz)
DICₙ = @lift interior(DIC[$n], :, :, grid.Nz)
fig = Figure(size = (1600, 1600), fontsize = 20)
lims = [(minimum(T), maximum(T)) for T in ( ζ[:, :, grid.Nz, :],
NO₃[:, :, grid.Nz, :] .+ NH₄[:, :, grid.Nz, :],
P[:, :, grid.Nz, :],
DIC[:, :, grid.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)
record(fig, "eady.mp4", 1:length(times), framerate = 12) do i
n[] = i
endThis page was generated using Literate.jl.