Baroclinic adjustment

In this example, we simulate the evolution and equilibration of a baroclinically unstable front.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, CairoMakie"

using Oceananigans
using Oceananigans.Units

Grid

We use a three-dimensional channel that is periodic in the x direction:

Lx = 1000kilometers # east-west extent [m]
Ly = 1000kilometers # north-south extent [m]
Lz = 1kilometers    # depth [m]

grid = RectilinearGrid(size = (48, 48, 8),
                       x = (0, Lx),
                       y = (-Ly/2, Ly/2),
                       z = (-Lz, 0),
                       topology = (Periodic, Bounded, Bounded))

48×48×8 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0e6)          regularly spaced with Δx=20833.3
├── Bounded  y ∈ [-500000.0, 500000.0] regularly spaced with Δy=20833.3
└── Bounded  z ∈ [-1000.0, 0.0]        regularly spaced with Δz=125.0

Model

We built a HydrostaticFreeSurfaceModel with an ImplicitFreeSurface solver. Regarding Coriolis, we use a beta-plane centered at 45° South.

model = HydrostaticFreeSurfaceModel(; grid,
                                    coriolis = BetaPlane(latitude = -45),
                                    buoyancy = BuoyancyTracer(),
                                    tracers = :b,
                                    momentum_advection = WENO(),
                                    tracer_advection = WENO())

HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 48×48×8 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
├── free surface: ImplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│   └── solver: FFTImplicitFreeSurfaceSolver
├── advection scheme: 
│   ├── momentum: WENO(order=5)
│   └── b: WENO(order=5)
└── coriolis: BetaPlane{Float64}

We start our simulation from rest with a baroclinically unstable buoyancy distribution. We use ramp(y, Δy), defined below, to specify a front with width Δy and horizontal buoyancy gradient M². We impose the front on top of a vertical buoyancy gradient N² and a bit of noise.

"""
    ramp(y, Δy)

Linear ramp from 0 to 1 between -Δy/2 and +Δy/2.

For example:
```
            y < -Δy/2 => ramp = 0
    -Δy/2 < y < -Δy/2 => ramp = y / Δy
            y >  Δy/2 => ramp = 1
```
"""
ramp(y, Δy) = min(max(0, y/Δy + 1/2), 1)

N² = 1e-5 # [s⁻²] buoyancy frequency / stratification
M² = 1e-7 # [s⁻²] horizontal buoyancy gradient

Δy = 100kilometers # width of the region of the front
Δb = Δy * M²       # buoyancy jump associated with the front
ϵb = 1e-2 * Δb     # noise amplitude

bᵢ(x, y, z) = N² * z + Δb * ramp(y, Δy) + ϵb * randn()

set!(model, b=bᵢ)

Let's visualize the initial buoyancy distribution.

using CairoMakie

# Build coordinates with units of kilometers
x, y, z = 1e-3 .* nodes(grid, (Center(), Center(), Center()))

b = model.tracers.b

fig, ax, hm = heatmap(view(b, 1, :, :),
                      colormap = :deep,
                      axis = (xlabel = "y [km]",
                              ylabel = "z [km]",
                              title = "b(x=0, y, z, t=0)",
                              titlesize = 24))

Colorbar(fig[1, 2], hm, label = "[m s⁻²]")

fig

Simulation

Now let's build a Simulation.

simulation = Simulation(model, Δt=20minutes, stop_time=20days)

Simulation of HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 20 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 20 days
├── Stop iteration : Inf
├── Wall time limit: Inf
├── 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
└── Diagnostics: OrderedDict with no entries

We add a TimeStepWizard callback to adapt the simulation's time-step,

conjure_time_step_wizard!(simulation, IterationInterval(20), cfl=0.2, max_Δt=20minutes)

Also, we add a callback to print a message about how the simulation is going,

using Printf

wall_clock = Ref(time_ns())

function print_progress(sim)
    u, v, w = model.velocities
    progress = 100 * (time(sim) / sim.stop_time)
    elapsed = (time_ns() - wall_clock[]) / 1e9

    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            progress, iteration(sim), prettytime(sim), prettytime(elapsed),
            maximum(abs, u), maximum(abs, v), maximum(abs, w), prettytime(sim.Δt))

    wall_clock[] = time_ns()

    return nothing
end

add_callback!(simulation, print_progress, IterationInterval(100))

Diagnostics/Output

Here, we save the buoyancy, $b$, at the edges of our domain as well as the zonal ($x$) average of buoyancy.

u, v, w = model.velocities
ζ = ∂x(v) - ∂y(u)
B = Average(b, dims=1)
U = Average(u, dims=1)
V = Average(v, dims=1)

filename = "baroclinic_adjustment"
save_fields_interval = 0.5day

slicers = (east = (grid.Nx, :, :),
           north = (:, grid.Ny, :),
           bottom = (:, :, 1),
           top = (:, :, grid.Nz))

for side in keys(slicers)
    indices = slicers[side]

    simulation.output_writers[side] = JLD2OutputWriter(model, (; b, ζ);
                                                       filename = filename * "_$(side)_slice",
                                                       schedule = TimeInterval(save_fields_interval),
                                                       overwrite_existing = true,
                                                       indices)
end

simulation.output_writers[:zonal] = JLD2OutputWriter(model, (; b=B, u=U, v=V);
                                                     filename = filename * "_zonal_average",
                                                     schedule = TimeInterval(save_fields_interval),
                                                     overwrite_existing = true)

JLD2OutputWriter scheduled on TimeInterval(12 hours):
├── filepath: baroclinic_adjustment_zonal_average.jld2
├── 3 outputs: (b, u, v)
├── array type: Array{Float64}
├── including: [:grid, :coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 31.6 KiB

Now we're ready to run.

@info "Running the simulation..."

run!(simulation)

@info "Simulation completed in " * prettytime(simulation.run_wall_time)

[ Info: Running the simulation...
[ Info: Initializing simulation...
[00.00%] i: 0, t: 0 seconds, wall time: 25.845 seconds, max(u): (0.000e+00, 0.000e+00, 0.000e+00) m/s, next Δt: 20 minutes
[ Info:     ... simulation initialization complete (24.328 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (18.478 seconds).
[06.94%] i: 100, t: 1.389 days, wall time: 35.541 seconds, max(u): (1.312e-01, 1.245e-01, 1.573e-03) m/s, next Δt: 20 minutes
[13.89%] i: 200, t: 2.778 days, wall time: 804.097 ms, max(u): (2.262e-01, 1.711e-01, 1.708e-03) m/s, next Δt: 20 minutes
[20.83%] i: 300, t: 4.167 days, wall time: 787.438 ms, max(u): (3.017e-01, 2.461e-01, 1.839e-03) m/s, next Δt: 20 minutes
[27.78%] i: 400, t: 5.556 days, wall time: 864.149 ms, max(u): (3.681e-01, 3.679e-01, 1.861e-03) m/s, next Δt: 20 minutes
[34.72%] i: 500, t: 6.944 days, wall time: 783.264 ms, max(u): (4.305e-01, 4.814e-01, 2.035e-03) m/s, next Δt: 20 minutes
[41.67%] i: 600, t: 8.333 days, wall time: 826.313 ms, max(u): (5.492e-01, 6.507e-01, 2.148e-03) m/s, next Δt: 20 minutes
[48.61%] i: 700, t: 9.722 days, wall time: 801.070 ms, max(u): (7.357e-01, 9.356e-01, 2.846e-03) m/s, next Δt: 20 minutes
[55.56%] i: 800, t: 11.111 days, wall time: 717.105 ms, max(u): (1.142e+00, 1.199e+00, 4.107e-03) m/s, next Δt: 20 minutes
[62.50%] i: 900, t: 12.500 days, wall time: 721.240 ms, max(u): (1.380e+00, 1.097e+00, 5.251e-03) m/s, next Δt: 20 minutes
[69.44%] i: 1000, t: 13.889 days, wall time: 749.771 ms, max(u): (1.468e+00, 1.115e+00, 4.187e-03) m/s, next Δt: 20 minutes
[76.39%] i: 1100, t: 15.278 days, wall time: 709.109 ms, max(u): (1.308e+00, 1.193e+00, 3.267e-03) m/s, next Δt: 20 minutes
[83.33%] i: 1200, t: 16.667 days, wall time: 676.973 ms, max(u): (1.376e+00, 1.100e+00, 2.512e-03) m/s, next Δt: 20 minutes
[90.28%] i: 1300, t: 18.056 days, wall time: 769.307 ms, max(u): (1.357e+00, 1.174e+00, 2.859e-03) m/s, next Δt: 20 minutes
[97.22%] i: 1400, t: 19.444 days, wall time: 712.918 ms, max(u): (1.357e+00, 1.388e+00, 3.605e-03) m/s, next Δt: 20 minutes
[ Info: Simulation is stopping after running for 56.551 seconds.
[ Info: Simulation time 20 days equals or exceeds stop time 20 days.
[ Info: Simulation completed in 56.583 seconds

Visualization

All that's left is to make a pretty movie. Actually, we make two visualizations here. First, we illustrate how to make a 3D visualization with Makie's Axis3 and Makie.surface. Then we make a movie in 2D. We use CairoMakie in this example, but note that using GLMakie is more convenient on a system with OpenGL, as figures will be displayed on the screen.

using CairoMakie

Three-dimensional visualization

We load the saved buoyancy output on the top, north, and east surface as FieldTimeSerieses.

filename = "baroclinic_adjustment"

sides = keys(slicers)

slice_filenames = NamedTuple(side => filename * "_$(side)_slice.jld2" for side in sides)

b_timeserieses = (east   = FieldTimeSeries(slice_filenames.east, "b"),
                  north  = FieldTimeSeries(slice_filenames.north, "b"),
                  top    = FieldTimeSeries(slice_filenames.top, "b"))

B_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "b")

times = B_timeseries.times
grid = B_timeseries.grid

48×48×8 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×3 halo
├── Periodic x ∈ [0.0, 1.0e6)          regularly spaced with Δx=20833.3
├── Bounded  y ∈ [-500000.0, 500000.0] regularly spaced with Δy=20833.3
└── Bounded  z ∈ [-1000.0, 0.0]        regularly spaced with Δz=125.0

We build the coordinates. We rescale horizontal coordinates to kilometers.

xb, yb, zb = nodes(b_timeserieses.east)

xb = xb ./ 1e3 # convert m -> km
yb = yb ./ 1e3 # convert m -> km

Nx, Ny, Nz = size(grid)

x_xz = repeat(x, 1, Nz)
y_xz_north = y[end] * ones(Nx, Nz)
z_xz = repeat(reshape(z, 1, Nz), Nx, 1)

x_yz_east = x[end] * ones(Ny, Nz)
y_yz = repeat(y, 1, Nz)
z_yz = repeat(reshape(z, 1, Nz), grid.Ny, 1)

x_xy = x
y_xy = y
z_xy_top = z[end] * ones(grid.Nx, grid.Ny)

Then we create a 3D axis. We use zonal_slice_displacement to control where the plot of the instantaneous zonal average flow is located.

fig = Figure(size = (1600, 800))

zonal_slice_displacement = 1.2

ax = Axis3(fig[2, 1],
           aspect=(1, 1, 1/5),
           xlabel = "x (km)",
           ylabel = "y (km)",
           zlabel = "z (m)",
           xlabeloffset = 100,
           ylabeloffset = 100,
           zlabeloffset = 100,
           limits = ((x[1], zonal_slice_displacement * x[end]), (y[1], y[end]), (z[1], z[end])),
           elevation = 0.45,
           azimuth = 6.8,
           xspinesvisible = false,
           zgridvisible = false,
           protrusions = 40,
           perspectiveness = 0.7)

Axis3()

We use data from the final savepoint for the 3D plot. Note that this plot can easily be animated by using Makie's Observable. To dive into Observables, check out Makie.jl's Documentation.

n = length(times)

41

Now let's make a 3D plot of the buoyancy and in front of it we'll use the zonally-averaged output to plot the instantaneous zonal-average of the buoyancy.

b_slices = (east   = interior(b_timeserieses.east[n], 1, :, :),
            north  = interior(b_timeserieses.north[n], :, 1, :),
            top    = interior(b_timeserieses.top[n], :, :, 1))

# Zonally-averaged buoyancy
B = interior(B_timeseries[n], 1, :, :)

clims = 1.1 .* extrema(b_timeserieses.top[n][:])

kwargs = (colorrange=clims, colormap=:deep, shading=NoShading)

surface!(ax, x_yz_east, y_yz, z_yz;  color = b_slices.east, kwargs...)
surface!(ax, x_xz, y_xz_north, z_xz; color = b_slices.north, kwargs...)
surface!(ax, x_xy, y_xy, z_xy_top;   color = b_slices.top, kwargs...)

sf = surface!(ax, zonal_slice_displacement .* x_yz_east, y_yz, z_yz; color = B, kwargs...)

contour!(ax, y, z, B; transformation = (:yz, zonal_slice_displacement * x[end]),
         levels = 15, linewidth = 2, color = :black)

Colorbar(fig[2, 2], sf, label = "m s⁻²", height = Relative(0.4), tellheight=false)

title = "Buoyancy at t = " * string(round(times[n] / day, digits=1)) * " days"
fig[1, 1:2] = Label(fig, title; fontsize = 24, tellwidth = false, padding = (0, 0, -120, 0))

rowgap!(fig.layout, 1, Relative(-0.2))
colgap!(fig.layout, 1, Relative(-0.1))

save("baroclinic_adjustment_3d.png", fig)

Two-dimensional movie

We make a 2D movie that shows buoyancy $b$ and vertical vorticity $ζ$ at the surface, as well as the zonally-averaged zonal and meridional velocities $U$ and $V$ in the $(y, z)$ plane. First we load the FieldTimeSeries and extract the additional coordinates we'll need for plotting

ζ_timeseries = FieldTimeSeries(slice_filenames.top, "ζ")
U_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "u")
B_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "b")
V_timeseries = FieldTimeSeries(filename * "_zonal_average.jld2", "v")

xζ, yζ, zζ = nodes(ζ_timeseries)
yv = ynodes(V_timeseries)

xζ = xζ ./ 1e3 # convert m -> km
yζ = yζ ./ 1e3 # convert m -> km
yv = yv ./ 1e3 # convert m -> km

49-element Vector{Float64}:
 -500.0
 -479.1666666666667
 -458.3333333333333
 -437.5
 -416.6666666666667
 -395.8333333333333
 -375.0
 -354.1666666666667
 -333.3333333333333
 -312.5
 -291.6666666666667
 -270.8333333333333
 -250.0
 -229.16666666666666
 -208.33333333333334
 -187.5
 -166.66666666666666
 -145.83333333333334
 -125.0
 -104.16666666666667
  -83.33333333333333
  -62.5
  -41.666666666666664
  -20.833333333333332
    0.0
   20.833333333333332
   41.666666666666664
   62.5
   83.33333333333333
  104.16666666666667
  125.0
  145.83333333333334
  166.66666666666666
  187.5
  208.33333333333334
  229.16666666666666
  250.0
  270.8333333333333
  291.6666666666667
  312.5
  333.3333333333333
  354.1666666666667
  375.0
  395.8333333333333
  416.6666666666667
  437.5
  458.3333333333333
  479.1666666666667
  500.0

Next, we set up a plot with 4 panels. The top panels are large and square, while the bottom panels get a reduced aspect ratio through rowsize!.

set_theme!(Theme(fontsize=24))

fig = Figure(size=(1800, 1000))

axb = Axis(fig[1, 2], xlabel="x (km)", ylabel="y (km)", aspect=1)
axζ = Axis(fig[1, 3], xlabel="x (km)", ylabel="y (km)", aspect=1, yaxisposition=:right)

axu = Axis(fig[2, 2], xlabel="y (km)", ylabel="z (m)")
axv = Axis(fig[2, 3], xlabel="y (km)", ylabel="z (m)", yaxisposition=:right)

rowsize!(fig.layout, 2, Relative(0.3))

To prepare a plot for animation, we index the timeseries with an Observable,

n = Observable(1)

b_top = @lift interior(b_timeserieses.top[$n], :, :, 1)
ζ_top = @lift interior(ζ_timeseries[$n], :, :, 1)
U = @lift interior(U_timeseries[$n], 1, :, :)
V = @lift interior(V_timeseries[$n], 1, :, :)
B = @lift interior(B_timeseries[$n], 1, :, :)

Observable([-0.009370714286853548 -0.008107382223240448 -0.0068704761641474365 -0.005614815081906362 -0.004403361139284586 -0.003142470612575701 -0.0018820627193703269 -0.0006221412687009454; -0.009384164507486184 -0.008125555547621887 -0.006891030195473702 -0.005620469853653107 -0.004359698070191914 -0.003099659420454004 -0.0018768196059675932 -0.0006180153622469953; -0.009379497143223994 -0.00812759008515642 -0.006866638475558347 -0.00561950742724678 -0.0043742750805188324 -0.00312270746315174 -0.0018931996159332765 -0.0006038467768436517; -0.009408363877129905 -0.00810323309305136 -0.0068921285238623745 -0.0056287590452781325 -0.004377623352613535 -0.0031139134090456175 -0.0018926170720303442 -0.000631544584983182; -0.009363894431592315 -0.008139348768753938 -0.006874324232811889 -0.005619420066172835 -0.004367781601045148 -0.003100450316576228 -0.0018900668315826796 -0.0006392176461966326; -0.009356679510805167 -0.008139665051477318 -0.006885805452892242 -0.00560404250570115 -0.004361093867736857 -0.003124051879568497 -0.0018762311950616406 -0.0006168432153860966; -0.009374136614533539 -0.008103147981759663 -0.006879428485977453 -0.005624622136479479 -0.004360683887004083 -0.0031189544079585564 -0.0018671481449340447 -0.0006385620903830123; -0.009358787198175489 -0.008133492829870052 -0.006861865172564604 -0.005602906848233415 -0.0043779921089819205 -0.003133211374986223 -0.001870621839667994 -0.0006397654138934715; -0.009358590990180264 -0.008108248610126078 -0.0068969040343947905 -0.005639310508611002 -0.004384863945731616 -0.003139117568485441 -0.0018842625738406896 -0.0006476255295283973; -0.009388075573280705 -0.00809951727794605 -0.006876424150390016 -0.0055993616335434536 -0.0043808557596287465 -0.0031303812820619274 -0.0018795263839168873 -0.0006200511929716864; -0.009379197610874655 -0.008136661386451116 -0.006876598542693932 -0.005619614084175281 -0.004397385820365898 -0.0031340334611436906 -0.0018942005335665663 -0.0005956209992063905; -0.009360517477648406 -0.00814295802916711 -0.0068753813871917065 -0.005623640092096732 -0.00436315370665388 -0.0031237735114468903 -0.0018661566412034028 -0.0006315319520222026; -0.009362301809061642 -0.008111295042002423 -0.006871246933301409 -0.0056017823159214275 -0.00438632917185995 -0.0031204213144566673 -0.0018836222442997184 -0.0006151472542203412; -0.009357078983182618 -0.00811844395674663 -0.006869805365155729 -0.00563886636623302 -0.004372604309645026 -0.0031452329377004234 -0.0018566081045587985 -0.000610397867178192; -0.009364310120721099 -0.008142322520278631 -0.006869354709951419 -0.005630396632080329 -0.004376731108171713 -0.0031332591245694986 -0.0018736967145424489 -0.0006206505705951438; -0.009364410571656623 -0.00811344818768709 -0.0068755413342187445 -0.005624890806021189 -0.004381378674495677 -0.0030933192081986313 -0.001904002548716899 -0.0006339507357221069; -0.00936474136501162 -0.00814922135098148 -0.006868893475658745 -0.005647109675659276 -0.004390603707089291 -0.003132255782145274 -0.0018705872753701739 -0.0006102497766675899; -0.009376547173054847 -0.008115238183147886 -0.006873612456093352 -0.00563752786320369 -0.00436301510236017 -0.003147828983452528 -0.0018613758947917105 -0.0006121230541251036; -0.009400130170396944 -0.008107601992181745 -0.006852763209792379 -0.0056304175822070104 -0.004393521200570923 -0.0031158808352383057 -0.0018787970383254478 -0.0006312414296736813; -0.009383017710692165 -0.008126951920425852 -0.006889385251736627 -0.005594217038445523 -0.0043871122513992 -0.0031391202991713344 -0.00185047108108823 -0.0006306457788392921; -0.009392394969404243 -0.00811137826070491 -0.006884049638465336 -0.005615202044961193 -0.004376597322175951 -0.0031415957939643392 -0.0018651965573836232 -0.0006255204280365875; -0.009389003996588994 -0.00811588683707575 -0.006862963884370289 -0.005606531276274967 -0.00436816255197352 -0.003118393487911411 -0.0018434670695093348 -0.0006032441844612287; -0.00749471721972772 -0.006252579890931612 -0.004999308685338589 -0.0037578502221680657 -0.002496166122672529 -0.00125583785293072 8.249004967683996e-6 0.0012324864976722182; -0.005421136701875442 -0.004167346294524058 -0.002905578926001804 -0.0016855532466705774 -0.00042346849893333587 0.0008758622364582221 0.002093601641010413 0.0033240062290908305; -0.00333934915697208 -0.0020885435996587723 -0.0008239733827774864 0.00040528708842631255 0.001651450929615535 0.0029436696241923535 0.004168951439122833 0.0054174357677177695; -0.0012236388871952823 8.595838768805774e-6 0.0012503575982497305 0.0024986003708883206 0.003746321401529433 0.0049971419948688514 0.0062857673548502606 0.007492791161334657; 0.000610530875867725 0.0018609945794045718 0.003138944300242238 0.004369732842708196 0.0056200314292535886 0.006876264765824124 0.008126331046674666 0.009366884167057868; 0.0006226208361089939 0.0018919606141538931 0.0031210581864724255 0.0043829772169775355 0.005629960788760227 0.006884737010923244 0.008100366976019305 0.009373817127773054; 0.000628125148776921 0.001879683130268893 0.0031136658795589793 0.004387820621207462 0.005618549075588773 0.006875001784686027 0.008140604771354443 0.009372849205687291; 0.0006128181510139632 0.0018479073070187882 0.0031017344335975646 0.004397397254087351 0.005653419664252368 0.006893313840472922 0.008118751927846826 0.009377118337631074; 0.000598237640080428 0.0018851143839445778 0.003110588154518503 0.004372125703107984 0.005619850264171953 0.006872607013840292 0.008118304614962322 0.00937776126674965; 0.0006096441203784907 0.0018841183751910618 0.0031377485687941278 0.004380445545221888 0.005608631137690434 0.0068635364468225 0.008122645823674576 0.009388337838312881; 0.0006071373655159732 0.00188037949026733 0.0031178240823313747 0.004361408963764854 0.005640267859091549 0.006900291729586089 0.008126983548996313 0.00939209039320754; 0.0006223440134425433 0.0018910439369302899 0.0031040919325891463 0.004372421864898174 0.005644394185815377 0.006863949885862734 0.008113843716107122 0.009395085168458125; 0.0006175981849776243 0.001878330031720889 0.0031316318046159986 0.004379060886389844 0.005630091054730617 0.006872698487101456 0.008116423773999461 0.009371912114634113; 0.0006474678653435505 0.0018748999509249481 0.003117515683234486 0.004353375533533474 0.00561591751115443 0.006877426677499745 0.008109870412569933 0.009368331619452197; 0.00063567100254849 0.0018776406055705172 0.0031246258558528446 0.004382886898472364 0.005644645724265991 0.006868988691731653 0.008120596636253202 0.009378403536521267; 0.0006392337826787543 0.0018465983102374701 0.0031123868173563492 0.004371683878544486 0.0056218339448897425 0.006881268466095258 0.00810000949881602 0.009364542789053862; 0.0006191539600330812 0.0018671016345360612 0.0031385616382176417 0.004391809504567911 0.005628015998309835 0.006871364498763355 0.008110543081735784 0.00937010633054611; 0.0006162269389263083 0.0018573023323233313 0.00311989449096209 0.004365534064697949 0.005623335346883887 0.006878528033071767 0.008114963411414744 0.009390667537255834; 0.0006166364632865963 0.0018587179049916471 0.0031360932303661974 0.004378056190575513 0.005631244817878178 0.0068531168448409605 0.008117680559395742 0.009398931856707154; 0.0006352655540465111 0.0018591055809809385 0.0031206427049220407 0.004363910782959294 0.005625948839821778 0.006855104092244192 0.008128115139302962 0.009397045184145676; 0.0006276047206541083 0.0018837935426281325 0.0031075623193142377 0.00435925677360753 0.00561264885266842 0.006893754464530062 0.00812597223542596 0.009384726609900418; 0.000632449973175823 0.0018709561693896273 0.003120199581539701 0.0043664775768926965 0.005643960519478293 0.006883790219901151 0.008087648996718773 0.009380117821996263; 0.0006150915563872117 0.0018514633411515682 0.003148005667812906 0.004380746573321521 0.005633582288047816 0.006881337541627802 0.00812812811432257 0.009410452186848675; 0.0006489576115477866 0.0019001961515279064 0.003135236588143485 0.004372834729926986 0.005642054948034815 0.006891182273244482 0.008133326664983889 0.009398263851125601; 0.0006176645306695701 0.001878734268714765 0.0031235767689220713 0.004374740513178681 0.0056054721788179876 0.006861201406935177 0.008105148568509546 0.00936465943736835; 0.0006273099030229836 0.0018692910388250807 0.003141798019693491 0.004406815954882758 0.0056244670197747315 0.006879924237557451 0.008119987385647178 0.009389131927712313])

and then build our plot:

hm = heatmap!(axb, xb, yb, b_top, colorrange=(0, Δb), colormap=:thermal)
Colorbar(fig[1, 1], hm, flipaxis=false, label="Surface b(x, y) (m s⁻²)")

hm = heatmap!(axζ, xζ, yζ, ζ_top, colorrange=(-5e-5, 5e-5), colormap=:balance)
Colorbar(fig[1, 4], hm, label="Surface ζ(x, y) (s⁻¹)")

hm = heatmap!(axu, yb, zb, U; colorrange=(-5e-1, 5e-1), colormap=:balance)
Colorbar(fig[2, 1], hm, flipaxis=false, label="Zonally-averaged U(y, z) (m s⁻¹)")
contour!(axu, yb, zb, B; levels=15, color=:black)

hm = heatmap!(axv, yv, zb, V; colorrange=(-1e-1, 1e-1), colormap=:balance)
Colorbar(fig[2, 4], hm, label="Zonally-averaged V(y, z) (m s⁻¹)")
contour!(axv, yb, zb, B; levels=15, color=:black)

Finally, we're ready to record the movie.

frames = 1:length(times)

record(fig, filename * ".mp4", frames, framerate=8) do i
    n[] = i
end

This page was generated using Literate.jl.