Baroclinic wave on the sphere

This example simulates the growth of a baroclinic wave on a near-global LatitudeLongitudeGrid following the DCMIP2016 specification Ullrich et al. (2017), which extends the classic Jablonowski and Williamson (2006) test case. A midlatitude jet in thermal-wind balance with a meridional temperature gradient is seeded with a localized zonal-wind perturbation that triggers baroclinic instability, producing growing Rossby waves over roughly ten days.

This example exercises CompressibleDynamics with SplitExplicitTimeDiscretization (acoustic substepping via AcousticRungeKutta3) and SphericalCoriolis (non-hydrostatic) on a 1° latitude-longitude grid spanning 75° S to 75° N. Acoustic substepping lets the outer time step be set by the advective CFL rather than the much-tighter acoustic CFL — here a time-step wizard floats Δt at advective CFL ≈ 1.4 against the polar Δx_min ≈ 28.8 km, capped at 12 min.

A future moist version (one-moment mixed-phase microphysics + bulk surface fluxes) will be added once the moist substepper supports the larger $Δt$ needed for a tractable runtime. The dry case here runs for 30 days and captures the full BCI life cycle: visible perturbations by day 4, explosive cyclogenesis around day 8, and saturation afterward.

Physical setup

The DCMIP2016 background state is an analytic steady-state solution of the dry, adiabatic, inviscid primitive equations in height coordinates, expressed in virtual temperature $T_v(\varphi, z)$:

\[T_v(φ, z) = \frac{1}{τ_1(z) - τ_2(z)\, F(φ)}\]

where $τ_1$ and $τ_2$ encode the vertical structure and $F(φ) = \cos^K φ - \frac{K}{K+2} \cos^{K+2} φ$ is the meridional shape with jet-width parameter $K = 3$. In the dry case, $T = T_v$.

Balanced zonal jet

The zonal wind is derived analytically from gradient-wind balance, producing a subtropical jet peaking near 30 m/s at 45° latitude in the upper troposphere.

Perturbation

A localized zonal-wind perturbation centered at $(λ_c, φ_c) = (20°\text{E}, 40°\text{N})$ seeds the instability. The perturbation decays exponentially with great-circle distance from the center and is tapered smoothly to zero above 15 km:

\[u'(λ, φ, z) = u_p \, \mathcal{T}(z) \, \exp\!\left(-\left(\frac{d}{r_p}\right)^2\right)\]

where $d$ is the great-circle distance, $r_p = 0.1\,a$, $u_p = 1$ m/s, and $\mathcal{T}(z) = 1 - 3(z/z_p)^2 + 2(z/z_p)^3$ for $z < z_p$.

using Breeze
using Oceananigans
using Oceananigans.Units
using Printf
using CairoMakie
using CUDA

DCMIP2016 parameters

All parameters follow the DCMIP2016 test case document Ullrich et al. (2017). We set the Oceananigans defaults and build a custom ThermodynamicConstants matching the DCMIP specification so that the grid, Coriolis, and model thermodynamics are all consistent with the analytic initial conditions.

Oceananigans.defaults.FloatType = Float32
Oceananigans.defaults.gravitational_acceleration = 9.80616
Oceananigans.defaults.planet_radius = 6371220
Oceananigans.defaults.planet_rotation_rate = 7.29212e-5

constants = ThermodynamicConstants(;
    gravitational_acceleration = Oceananigans.defaults.gravitational_acceleration,
    dry_air_heat_capacity = 1004.5,
    dry_air_molar_mass = 8.314462618 / 287)

g   = constants.gravitational_acceleration
Rᵈ  = dry_air_gas_constant(constants)
cᵖᵈ = constants.dry_air.heat_capacity
κ   = Rᵈ / cᵖᵈ
p₀  = 1e5    # Pa — surface pressure
a   = Oceananigans.defaults.planet_radius
Ω   = Oceananigans.defaults.planet_rotation_rate

7.29212e-5

Domain and grid

We use a 1° latitude-longitude grid spanning 75° S to 75° N. Capping the domain at ±75° (rather than the poles) keeps the polar Δx_min manageable: a · cos 75° · 2π/Nλ ≈ 28.8 km. The domain extends from the surface to 30 km with 32 vertical levels (Δz ≈ 940 m).

Nλ = 360
Nφ = 150
Nz = 32
H  = 30kilometers

grid = LatitudeLongitudeGrid(GPU();
                             size = (Nλ, Nφ, Nz),
                             halo = (5, 5, 5),
                             longitude = (0, 360),
                             latitude = (-75, 75),
                             z = (0, H))

# Temperature profile parameters
Tᴱ = 310     # K — equatorial surface temperature
Tᴾ = 240     # K — polar surface temperature
Tᴹ = (Tᴱ + Tᴾ) / 2
Γ  = 0.005    # K/m — lapse rate
K  = 3        # jet width parameter
b  = 2        # vertical half-width parameter

Analytic initial conditions

The DCMIP2016 balanced state is given in virtual temperature $T_v$. We define the IC in terms of $T_v$ so a future moist version (where $T = T_v / (1 + \epsilon q^v)$) can reuse the same functions.

# Vertical structure functions (shallow atmosphere, X = 1)
function τ_and_integrals(z)
    Hˢ = Rᵈ * Tᴹ / g
    η  = z / (b * Hˢ)
    e  = exp(-η^2)

    A = (Tᴹ - Tᴾ) / (Tᴹ * Tᴾ)
    C = (K + 2) * (Tᴱ - Tᴾ) / (2 * Tᴱ * Tᴾ)

    τ₁  = A * (1 - 2η^2) * e + exp(Γ * z / Tᴹ) / Tᴹ
    ∫τ₁ = A * z * e + (exp(Γ * z / Tᴹ) - 1) / Γ

    τ₂  = C * (1 - 2η^2) * e
    ∫τ₂ = C * z * e

    return τ₁, τ₂, ∫τ₁, ∫τ₂
end

# Meridional shape functions
F(φ)  = cosd(φ)^K - K / (K + 2) * cosd(φ)^(K + 2)
dF(φ) = cosd(φ)^(K - 1) - cosd(φ)^(K + 1)

# Virtual temperature: Tᵥ(φ, z) = 1 / (τ₁ - τ₂ F(φ))
function virtual_temperature(λ, φ, z)
    τ₁, τ₂, _, _ = τ_and_integrals(z)
    return 1 / (τ₁ - τ₂ * F(φ))
end

# Pressure: p(φ, z) = p₀ exp(-g/Rᵈ (∫τ₁ - ∫τ₂ F(φ)))
function pressure(λ, φ, z)
    _, _, ∫τ₁, ∫τ₂ = τ_and_integrals(z)
    return p₀ * exp(-g / Rᵈ * (∫τ₁ - ∫τ₂ * F(φ)))
end

# Density (uses Tᵥ in the ideal gas law; in the dry case, T = Tᵥ).
density(λ, φ, z) = pressure(λ, φ, z) / (Rᵈ * virtual_temperature(λ, φ, z))

# Potential temperature: θ = Tᵥ (p₀/p)^κ in the dry case.
potential_temperature(λ, φ, z) = virtual_temperature(λ, φ, z) * (p₀ / pressure(λ, φ, z))^κ

potential_temperature (generic function with 1 method)

Balanced zonal wind

The zonal wind satisfies gradient-wind balance with the temperature field. For the shallow atmosphere ($r = a$):

\[u = -Ω a \cos φ + \sqrt{Ω^2 a^2 \cos^2 φ + a \cos φ \, U(φ, z)}\]

where $U = (g/a) K \int τ_2 \, T_v \, (\cos^{K-1} φ - \cos^{K+1} φ)$.

function zonal_velocity(λ, φ, z)
    _, _, _, ∫τ₂ = τ_and_integrals(z)
    Tᵥ = virtual_temperature(λ, φ, z)

    # Gradient-wind balance
    U = g / a * K * ∫τ₂ * dF(φ) * Tᵥ
    rcosφ  = a * cosd(φ)
    Ωrcosφ = Ω * rcosφ
    u_balanced = -Ωrcosφ + sqrt(Ωrcosφ^2 + rcosφ * U)

    # Localized perturbation (DCMIP2016 §3.3)
    uₚ = 1         # m/s — amplitude
    rₚ = 0.1       # perturbation radius (Earth radii)
    λₚ = π / 9     # 20°E center longitude
    φₚ = 2π / 9    # 40°N center latitude
    zₚ = 15000     # m — height cap

    φʳ = deg2rad(φ)
    λʳ = deg2rad(λ)
    great_circle = acos(sin(φₚ) * sin(φʳ) + cos(φₚ) * cos(φʳ) * cos(λʳ - λₚ)) / rₚ

    taper = ifelse(z < zₚ, 1 - 3 * (z / zₚ)^2 + 2 * (z / zₚ)^3, zero(z))
    u_perturbation = ifelse(great_circle < 1, uₚ * taper * exp(-great_circle^2), zero(z))

    return u_balanced + u_perturbation
end

zonal_velocity (generic function with 1 method)

Model configuration

We use fully compressible dynamics with acoustic substepping via SplitExplicitTimeDiscretization and the AcousticRungeKutta3 (Wicker–Skamarock RK3) outer loop. Acoustic substepping handles the fast acoustic-mode pressure gradient and buoyancy via a vertically-implicit inner loop, so the outer time step is limited only by the advective CFL — about 100× larger than the acoustic-CFL-limited Δt the fully explicit solver requires for the same grid.

We use a hydrostatically-balanced isothermal reference state at T₀ᵣ = 250 K (matching the MPAS convention) so that the substepper's slow tendencies see only perturbations from the background. SphericalCoriolis is the non-hydrostatic spherical form, which retains both the traditional $f = 2Ω \sin φ$ and the non-traditional $2Ω \cos φ$ cross-terms that couple horizontal momentum to $w$. Breeze evolves prognostic $ρw$ so the non-traditional terms are physically required for self-consistent dynamics on the sphere.

coriolis = SphericalCoriolis(rotation_rate=Ω)

T₀ᵣ = 250
θᵣ(z) = T₀ᵣ * exp(g * z / (cᵖᵈ * T₀ᵣ))

dynamics = CompressibleDynamics(SplitExplicitTimeDiscretization();
                                surface_pressure = p₀,
                                reference_potential_temperature = θᵣ)

model = AtmosphereModel(grid; dynamics, coriolis,
                        thermodynamic_constants = constants,
                        advection = WENO())

AtmosphereModel{GPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
├── grid: 360×150×32 LatitudeLongitudeGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on CUDAGPU with 5×5×5 halo
├── dynamics: CompressibleDynamics{SplitExplicitTimeDiscretization}
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float32}
├── timestepper: AcousticRungeKutta3
├── advection scheme: 
│   ├── momentum: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   ├── ρθ: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│   └── ρqᵛ: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
├── forcing: @NamedTuple{ρ::Returns{Float32}, ρu::Returns{Float32}, ρv::Returns{Float32}, ρw::Returns{Float32}, ρθ::Returns{Float32}, ρqᵛ::Returns{Float32}, ρe::Returns{Float32}}
├── tracers: ()
├── coriolis: Oceananigans.Coriolis.SphericalCoriolis{Oceananigans.Advection.EnstrophyConserving{Float32}, Float32, Oceananigans.Coriolis.NonhydrostaticFormulation}
└── microphysics: Nothing

Set initial conditions

set!(model, θ=potential_temperature, u=zonal_velocity, ρ=density)

Time-stepping

Substepping eliminates the acoustic CFL constraint on the outer Δt; only the advective CFL remains. A time-step wizard targets advective CFL ≈ 1.4 against the polar Δx_min ≈ 28.8 km, capped at Δt = 12 min:

\[Δt = \min\!\left(1.4 \cdot Δx_{\min} / U_{\max},\ 720 \text{ s}\right).\]

This is many times larger than the acoustic-CFL-limited Δt a fully explicit solver would require. We run for 30 days to capture the full BCI life cycle.

Δt = 12minutes
stop_time = 30days

simulation = Simulation(model; Δt, stop_time)
conjure_time_step_wizard!(simulation; cfl=1.4, max_Δt=12minutes)
Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)

Progress callback

function progress(sim)
    u, v, w = sim.model.velocities
    @info @sprintf("Iter %5d | t = %s | Δt = %s | max|u| = %.1f m/s | max|w| = %.4f m/s",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt),
                   maximum(abs, u), maximum(abs, w))
    return nothing
end

add_callback!(simulation, progress, IterationInterval(50))

Output

We save the velocities, the full potential temperature $θ$ (the classic surface synoptic diagnostic for the cold/warm sectors during cyclogenesis), and the vertical vorticity $ζ$, sliced at two levels: k = 1 (surface) and k = 16 (mid-troposphere, ~5 km).

using Oceananigans.Operators: ζ₃ᶠᶠᶜ
u, v, w = model.velocities
ζ = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, model.grid, u, v)

θ = PotentialTemperature(model)

outputs = merge(model.velocities, (; ζ, θ))

for k in (1, 16)
    filename = "baroclinic_wave_k$k"
    ow = JLD2Writer(model, outputs; filename,
                    indices = (:, :, k),
                    schedule = TimeInterval(6hours),
                    overwrite_existing = true)

    simulation.output_writers[Symbol(filename)] = ow
end

Run

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter     0 | t = 0 seconds | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0000 m/s
[ Info:     ... simulation initialization complete (1.123 minutes)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (35.410 seconds).
[ Info: Iter    50 | t = 10 hours | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0114 m/s
[ Info: Iter   100 | t = 20 hours | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0093 m/s
[ Info: Iter   150 | t = 1.250 days | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0075 m/s
[ Info: Iter   200 | t = 1.667 days | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0063 m/s
[ Info: Iter   250 | t = 2.083 days | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0053 m/s
[ Info: Iter   300 | t = 2.500 days | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0048 m/s
[ Info: Iter   350 | t = 2.917 days | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0040 m/s
[ Info: Iter   400 | t = 3.333 days | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0033 m/s
[ Info: Iter   450 | t = 3.750 days | Δt = 12 minutes | max|u| = 28.2 m/s | max|w| = 0.0028 m/s
[ Info: Iter   500 | t = 4.167 days | Δt = 12 minutes | max|u| = 28.2 m/s | max|w| = 0.0028 m/s
[ Info: Iter   550 | t = 4.583 days | Δt = 12 minutes | max|u| = 28.3 m/s | max|w| = 0.0037 m/s
[ Info: Iter   600 | t = 5 days | Δt = 12 minutes | max|u| = 28.3 m/s | max|w| = 0.0051 m/s
[ Info: Iter   650 | t = 5.417 days | Δt = 12 minutes | max|u| = 28.3 m/s | max|w| = 0.0066 m/s
[ Info: Iter   700 | t = 5.833 days | Δt = 12 minutes | max|u| = 28.5 m/s | max|w| = 0.0086 m/s
[ Info: Iter   750 | t = 6.250 days | Δt = 12 minutes | max|u| = 28.6 m/s | max|w| = 0.0109 m/s
[ Info: Iter   800 | t = 6.667 days | Δt = 12 minutes | max|u| = 28.7 m/s | max|w| = 0.0156 m/s
[ Info: Iter   850 | t = 7.083 days | Δt = 12 minutes | max|u| = 28.8 m/s | max|w| = 0.0237 m/s
[ Info: Iter   900 | t = 7.500 days | Δt = 12 minutes | max|u| = 29.0 m/s | max|w| = 0.0336 m/s
[ Info: Iter   950 | t = 7.917 days | Δt = 12 minutes | max|u| = 29.2 m/s | max|w| = 0.0532 m/s
[ Info: Iter  1000 | t = 8.333 days | Δt = 12 minutes | max|u| = 29.4 m/s | max|w| = 0.0584 m/s
[ Info: Iter  1050 | t = 8.750 days | Δt = 12 minutes | max|u| = 29.7 m/s | max|w| = 0.0846 m/s
[ Info: Iter  1100 | t = 9.167 days | Δt = 12 minutes | max|u| = 30.7 m/s | max|w| = 0.0831 m/s
[ Info: Iter  1150 | t = 9.583 days | Δt = 12 minutes | max|u| = 35.8 m/s | max|w| = 0.0895 m/s
[ Info: Iter  1200 | t = 10 days | Δt = 12 minutes | max|u| = 41.8 m/s | max|w| = 0.0909 m/s
[ Info: Iter  1250 | t = 10.417 days | Δt = 12 minutes | max|u| = 48.8 m/s | max|w| = 0.1009 m/s
[ Info: Iter  1300 | t = 10.833 days | Δt = 12 minutes | max|u| = 57.2 m/s | max|w| = 0.1016 m/s
[ Info: Iter  1350 | t = 11.250 days | Δt = 12 minutes | max|u| = 63.0 m/s | max|w| = 0.1198 m/s
[ Info: Iter  1400 | t = 11.667 days | Δt = 12 minutes | max|u| = 60.6 m/s | max|w| = 0.3247 m/s
[ Info: Iter  1450 | t = 12.083 days | Δt = 12 minutes | max|u| = 57.8 m/s | max|w| = 0.2659 m/s
[ Info: Iter  1500 | t = 12.500 days | Δt = 12 minutes | max|u| = 52.8 m/s | max|w| = 0.4095 m/s
[ Info: Iter  1550 | t = 12.917 days | Δt = 12 minutes | max|u| = 62.4 m/s | max|w| = 0.4937 m/s
[ Info: Iter  1600 | t = 13.326 days | Δt = 10.413 minutes | max|u| = 66.9 m/s | max|w| = 0.4117 m/s
[ Info: Iter  1650 | t = 13.675 days | Δt = 9.860 minutes | max|u| = 70.1 m/s | max|w| = 0.5008 m/s
[ Info: Iter  1700 | t = 14.007 days | Δt = 9.656 minutes | max|u| = 72.1 m/s | max|w| = 0.4166 m/s
[ Info: Iter  1750 | t = 14.333 days | Δt = 9.972 minutes | max|u| = 72.3 m/s | max|w| = 0.6947 m/s
[ Info: Iter  1800 | t = 14.703 days | Δt = 11.046 minutes | max|u| = 68.2 m/s | max|w| = 0.6376 m/s
[ Info: Iter  1850 | t = 15.077 days | Δt = 11.162 minutes | max|u| = 68.0 m/s | max|w| = 0.7555 m/s
[ Info: Iter  1900 | t = 15.457 days | Δt = 10.969 minutes | max|u| = 65.7 m/s | max|w| = 0.8603 m/s
[ Info: Iter  1950 | t = 15.821 days | Δt = 10.204 minutes | max|u| = 67.6 m/s | max|w| = 1.0684 m/s
[ Info: Iter  2000 | t = 16.188 days | Δt = 10.564 minutes | max|u| = 68.2 m/s | max|w| = 0.7622 m/s
[ Info: Iter  2050 | t = 16.566 days | Δt = 12 minutes | max|u| = 68.5 m/s | max|w| = 0.6768 m/s
[ Info: Iter  2100 | t = 16.975 days | Δt = 12 minutes | max|u| = 70.7 m/s | max|w| = 0.6630 m/s
[ Info: Iter  2150 | t = 17.352 days | Δt = 10.417 minutes | max|u| = 66.6 m/s | max|w| = 0.6532 m/s
[ Info: Iter  2200 | t = 17.712 days | Δt = 10.626 minutes | max|u| = 65.0 m/s | max|w| = 0.5332 m/s
[ Info: Iter  2250 | t = 18.077 days | Δt = 11.446 minutes | max|u| = 61.5 m/s | max|w| = 0.6506 m/s
[ Info: Iter  2300 | t = 18.483 days | Δt = 12 minutes | max|u| = 59.5 m/s | max|w| = 0.9059 m/s
[ Info: Iter  2350 | t = 18.900 days | Δt = 12 minutes | max|u| = 60.4 m/s | max|w| = 0.5071 m/s
[ Info: Iter  2400 | t = 19.317 days | Δt = 12 minutes | max|u| = 67.6 m/s | max|w| = 0.3361 m/s
[ Info: Iter  2450 | t = 19.733 days | Δt = 12 minutes | max|u| = 67.1 m/s | max|w| = 0.6113 m/s
[ Info: Iter  2500 | t = 20.150 days | Δt = 12 minutes | max|u| = 69.7 m/s | max|w| = 0.4235 m/s
[ Info: Iter  2550 | t = 20.567 days | Δt = 12 minutes | max|u| = 72.0 m/s | max|w| = 0.3198 m/s
[ Info: Iter  2600 | t = 20.983 days | Δt = 12 minutes | max|u| = 77.0 m/s | max|w| = 0.3053 m/s
[ Info: Iter  2650 | t = 21.400 days | Δt = 12 minutes | max|u| = 76.5 m/s | max|w| = 0.3336 m/s
[ Info: Iter  2700 | t = 21.770 days | Δt = 9.119 minutes | max|u| = 74.9 m/s | max|w| = 0.2524 m/s
[ Info: Iter  2750 | t = 22.094 days | Δt = 11.152 minutes | max|u| = 71.3 m/s | max|w| = 0.3300 m/s
[ Info: Iter  2800 | t = 22.500 days | Δt = 12 minutes | max|u| = 72.4 m/s | max|w| = 0.7530 m/s
[ Info: Iter  2850 | t = 22.917 days | Δt = 10.238 minutes | max|u| = 69.3 m/s | max|w| = 0.8333 m/s
[ Info: Iter  2900 | t = 23.308 days | Δt = 12 minutes | max|u| = 69.6 m/s | max|w| = 0.3375 m/s
[ Info: Iter  2950 | t = 23.725 days | Δt = 12 minutes | max|u| = 74.7 m/s | max|w| = 0.3334 m/s
[ Info: Iter  3000 | t = 24.111 days | Δt = 9.205 minutes | max|u| = 72.7 m/s | max|w| = 0.7853 m/s
[ Info: Iter  3050 | t = 24.394 days | Δt = 7.498 minutes | max|u| = 78.3 m/s | max|w| = 1.0374 m/s
[ Info: Iter  3100 | t = 24.655 days | Δt = 7.371 minutes | max|u| = 84.6 m/s | max|w| = 1.0794 m/s
[ Info: Iter  3150 | t = 24.901 days | Δt = 7.168 minutes | max|u| = 95.4 m/s | max|w| = 0.9790 m/s
[ Info: Iter  3200 | t = 25.140 days | Δt = 6.903 minutes | max|u| = 99.8 m/s | max|w| = 1.1483 m/s
[ Info: Iter  3250 | t = 25.371 days | Δt = 6.313 minutes | max|u| = 107.7 m/s | max|w| = 0.8381 m/s
[ Info: Iter  3300 | t = 25.575 days | Δt = 5.549 minutes | max|u| = 117.9 m/s | max|w| = 1.0251 m/s
[ Info: Iter  3350 | t = 25.771 days | Δt = 6.587 minutes | max|u| = 101.7 m/s | max|w| = 1.5175 m/s
[ Info: Iter  3400 | t = 26.000 days | Δt = 6.694 minutes | max|u| = 94.2 m/s | max|w| = 1.7374 m/s
[ Info: Iter  3450 | t = 26.223 days | Δt = 6.988 minutes | max|u| = 87.0 m/s | max|w| = 1.2804 m/s
[ Info: Iter  3500 | t = 26.462 days | Δt = 8.044 minutes | max|u| = 77.7 m/s | max|w| = 0.6205 m/s
[ Info: Iter  3550 | t = 26.737 days | Δt = 7.155 minutes | max|u| = 87.7 m/s | max|w| = 1.1382 m/s
[ Info: Iter  3600 | t = 26.978 days | Δt = 6.424 minutes | max|u| = 107.9 m/s | max|w| = 0.8559 m/s
[ Info: Iter  3650 | t = 27.199 days | Δt = 6.353 minutes | max|u| = 109.1 m/s | max|w| = 0.9188 m/s
[ Info: Iter  3700 | t = 27.426 days | Δt = 6.952 minutes | max|u| = 101.7 m/s | max|w| = 0.8186 m/s
[ Info: Iter  3750 | t = 27.667 days | Δt = 6.739 minutes | max|u| = 105.4 m/s | max|w| = 1.5611 m/s
[ Info: Iter  3800 | t = 27.860 days | Δt = 4.630 minutes | max|u| = 129.4 m/s | max|w| = 1.0139 m/s
[ Info: Iter  3850 | t = 28.029 days | Δt = 5.609 minutes | max|u| = 109.1 m/s | max|w| = 1.1721 m/s
[ Info: Iter  3900 | t = 28.237 days | Δt = 6.559 minutes | max|u| = 101.4 m/s | max|w| = 2.3608 m/s
[ Info: Iter  3950 | t = 28.450 days | Δt = 5.889 minutes | max|u| = 100.2 m/s | max|w| = 1.9176 m/s
[ Info: Iter  4000 | t = 28.664 days | Δt = 6.465 minutes | max|u| = 104.1 m/s | max|w| = 1.3735 m/s
[ Info: Iter  4050 | t = 28.878 days | Δt = 6.715 minutes | max|u| = 101.1 m/s | max|w| = 1.2962 m/s
[ Info: Iter  4100 | t = 29.091 days | Δt = 5.550 minutes | max|u| = 124.9 m/s | max|w| = 1.3244 m/s
[ Info: Iter  4150 | t = 29.276 days | Δt = 4.628 minutes | max|u| = 149.8 m/s | max|w| = 2.1986 m/s
[ Info: Iter  4200 | t = 29.443 days | Δt = 4.865 minutes | max|u| = 142.5 m/s | max|w| = 2.5609 m/s
[ Info: Iter  4250 | t = 29.609 days | Δt = 4.399 minutes | max|u| = 126.9 m/s | max|w| = 3.2124 m/s
[ Info: Iter  4300 | t = 29.764 days | Δt = 5.495 minutes | max|u| = 120.0 m/s | max|w| = 2.6375 m/s
[ Info: Iter  4350 | t = 29.962 days | Δt = 5.140 minutes | max|u| = 110.5 m/s | max|w| = 1.5379 m/s
[ Info: Simulation is stopping after running for 20.915 minutes.
[ Info: Simulation time 30 days equals or exceeds stop time 30 days.

Visualization

We plot three diagnostics on the sphere: the surface potential temperature $θ_{\rm sfc}$ (the classic synoptic diagnostic for the cold/warm sectors), the surface vertical vorticity $ζ$ (which reveals the cyclones and anticyclones), and the mid-level vertical velocity $w$ (which highlights the warm conveyor belt and the wave's vertical structure).

θ_ts = FieldTimeSeries("baroclinic_wave_k1.jld2",  "θ")
ζ_ts = FieldTimeSeries("baroclinic_wave_k1.jld2",  "ζ")
w_ts = FieldTimeSeries("baroclinic_wave_k16.jld2", "w")
times = θ_ts.times
Nt = length(times)

k_sfc = 1
k_mid = 16

Sphere view: rotate so the developing wave faces the camera.

sphere_kw = (elevation = π/6, azimuth = π/2, aspect = :data)
ζlim = 1e-4
wlim = 0.06

θ_kw = (colormap = :thermal, colorrange = (260, 310))
ζ_kw = (colormap = :balance, colorrange = (-ζlim, ζlim))
w_kw = (colormap = :balance, colorrange = (-wlim, wlim))

(colormap = :balance, colorrange = (-0.06, 0.06))

Animation

n = Observable(1)
θn = @lift view(θ_ts[$n], :, :, k_sfc)
ζn = @lift view(ζ_ts[$n], :, :, k_sfc)
wn = @lift view(w_ts[$n], :, :, k_mid)

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

title = @lift "t = $(prettytime(times[$n]))"
fig[0, 1:6] = Label(fig, title, fontsize=22, tellwidth=false)

ax1 = Axis3(fig[1, 1]; title = "θ at surface", sphere_kw...)
hm1 = surface!(ax1, θn; shading = NoShading, θ_kw...)
Colorbar(fig[1, 2], hm1; label = "θ (K)", height=Relative(0.5))

ax2 = Axis3(fig[1, 3]; title = "ζ at surface", sphere_kw...)
hm2 = surface!(ax2, ζn; shading = NoShading, ζ_kw...)
Colorbar(fig[1, 4], hm2; label = "ζ (1/s)", height=Relative(0.5))

ax3 = Axis3(fig[1, 5]; title = "w at mid-level", sphere_kw...)
hm3 = surface!(ax3, wn; shading = NoShading, w_kw...)
Colorbar(fig[1, 6], hm3; label = "w (m/s)", height=Relative(0.5))

for ax in (ax1, ax2, ax3)
    hidedecorations!(ax)
    hidespines!(ax)
end

CairoMakie.record(fig, "baroclinic_wave.mp4", 1:Nt; framerate = 12) do nn
    n[] = nn
end

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.6
Commit 15346901f00 (2026-04-09 19:20 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD EPYC 7R13 Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_GPG = 3673DF529D9049477F76B37566E3C7DC03D6E495
  JULIA_LOAD_PATH = :@breeze
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_VERSION = 1.12.6
  JULIA_DEPOT_PATH = /usr/local/share/julia:
  JULIA_PATH = /usr/local/julia
  JULIA_PROJECT = @breeze

These were the top-level packages installed in the environment:

import Pkg
Pkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
  [86bc3604] AtmosphericProfilesLibrary v0.1.7
  [660aa2fb] Breeze v0.5.3 `.`
⌅ [052768ef] CUDA v6.0.0
  [13f3f980] CairoMakie v0.15.11
⌅ [6a9e3e04] CloudMicrophysics v0.35.0
  [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.4.1
  [7da242da] Enzyme v0.13.154
  [46192b85] GPUArraysCore v0.2.0
  [63c18a36] KernelAbstractions v0.9.41
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.110.1
  [a01a1ee8] RRTMGP v0.21.9
  [3c362404] Reactant v0.2.264
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.1
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.