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 Breezeusing Oceananigansusing Oceananigans.Unitsusing Printfusing CairoMakieusing 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 = Float32Oceananigans.defaults.gravitational_acceleration = 9.80616Oceananigans.defaults.planet_radius = 6371220Oceananigans.defaults.planet_rotation_rate = 7.29212e-5constants = ThermodynamicConstants(;    gravitational_acceleration = Oceananigans.defaults.gravitational_acceleration,    dry_air_heat_capacity = 1004.5,    dry_air_molar_mass = 8.314462618 / 287)g   = constants.gravitational_accelerationRᵈ  = dry_air_gas_constant(constants)cᵖᵈ = constants.dry_air.heat_capacityκ   = Rᵈ / cᵖᵈp₀  = 1e5    # Pa — surface pressurea   = 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 64 vertical levels, exponentially stretched toward the surface with ExponentialDiscretization: the interfaces are clustered near the ground (bias = :left) so the smallest cells sit at the surface (Δz ≈ 150 m) and coarsen to ≈ 1070 m at the model top. The e-folding scale = H/2 sets how quickly the spacing grows with height.

Nλ = 360Nφ = 150Nz = 64H  = 30kilometersz_faces = ExponentialDiscretization(Nz, 0, H; scale = H/2, bias = :left)grid = LatitudeLongitudeGrid(GPU();                             size = (Nλ, Nφ, Nz),                             halo = (3, 3, 3),                             longitude = (0, 360),                             latitude = (-75, 75),                             z = z_faces)

360×150×64 LatitudeLongitudeGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 halo
├── longitude: Periodic λ ∈ [0.0, 360.0)   regularly spaced with Δλ=1.0
├── latitude:  Bounded  φ ∈ [-75.0, 75.0]  regularly spaced with Δφ=1.0
└── z:         Bounded  z ∈ [0.0, 30000.0] variably spaced with min(Δz)=149.052, max(Δz)=1067.47

Temperature profile parameters

Tᴱ = 310     # K — equatorial surface temperatureTᴾ = 240     # K — polar surface temperatureTᴹ = (Tᴱ + Tᴾ) / 2Γ  = 0.005    # K/m — lapse rateK  = 3        # jet width parameterb  = 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 functionsF(φ)  = 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_perturbationend

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.

Tracer and momentum advection uses fifth-order WENO reconstruction. No explicit closure is applied: WENO's implicit dissipation suffices on this poleward-refined, near-surface-stretched grid.

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(order=5))

AtmosphereModel{GPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
├── grid: 360×150×64 LatitudeLongitudeGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 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 = 12minutesstop_time = 30dayscfl = 1.4simulation = Simulation(model; Δt, stop_time)conjure_time_step_wizard!(simulation; cfl, 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 nothingendadd_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 (lower troposphere, ~2.9 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)] = owend

Run

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter     0 | t = 0 seconds | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0000 m/s
[ Info:     ... simulation initialization complete (39.522 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (29.511 seconds).
[ Info: Iter    50 | t = 10 hours | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0021 m/s
[ Info: Iter   100 | t = 20 hours | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0013 m/s
[ Info: Iter   150 | t = 1.250 days | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0016 m/s
[ Info: Iter   200 | t = 1.667 days | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0013 m/s
[ Info: Iter   250 | t = 2.083 days | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0011 m/s
[ Info: Iter   300 | t = 2.500 days | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0013 m/s
[ Info: Iter   350 | t = 2.917 days | Δt = 12 minutes | max|u| = 28.0 m/s | max|w| = 0.0011 m/s
[ Info: Iter   400 | t = 3.333 days | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0012 m/s
[ Info: Iter   450 | t = 3.750 days | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0014 m/s
[ Info: Iter   500 | t = 4.167 days | Δt = 12 minutes | max|u| = 28.1 m/s | max|w| = 0.0018 m/s
[ Info: Iter   550 | t = 4.583 days | Δt = 12 minutes | max|u| = 28.2 m/s | max|w| = 0.0024 m/s
[ Info: Iter   600 | t = 5 days | Δt = 12 minutes | max|u| = 28.2 m/s | max|w| = 0.0030 m/s
[ Info: Iter   650 | t = 5.417 days | Δt = 12 minutes | max|u| = 28.3 m/s | max|w| = 0.0040 m/s
[ Info: Iter   700 | t = 5.833 days | Δt = 12 minutes | max|u| = 28.4 m/s | max|w| = 0.0054 m/s
[ Info: Iter   750 | t = 6.250 days | Δt = 12 minutes | max|u| = 28.5 m/s | max|w| = 0.0073 m/s
[ Info: Iter   800 | t = 6.667 days | Δt = 12 minutes | max|u| = 28.6 m/s | max|w| = 0.0102 m/s
[ Info: Iter   850 | t = 7.083 days | Δt = 12 minutes | max|u| = 28.7 m/s | max|w| = 0.0145 m/s
[ Info: Iter   900 | t = 7.500 days | Δt = 12 minutes | max|u| = 28.9 m/s | max|w| = 0.0209 m/s
[ Info: Iter   950 | t = 7.917 days | Δt = 12 minutes | max|u| = 29.2 m/s | max|w| = 0.0309 m/s
[ Info: Iter  1000 | t = 8.333 days | Δt = 12 minutes | max|u| = 29.4 m/s | max|w| = 0.0477 m/s
[ Info: Iter  1050 | t = 8.750 days | Δt = 12 minutes | max|u| = 29.6 m/s | max|w| = 0.0627 m/s
[ Info: Iter  1100 | t = 9.167 days | Δt = 12 minutes | max|u| = 31.4 m/s | max|w| = 0.0767 m/s
[ Info: Iter  1150 | t = 9.583 days | Δt = 12 minutes | max|u| = 38.5 m/s | max|w| = 0.0833 m/s
[ Info: Iter  1200 | t = 10 days | Δt = 12 minutes | max|u| = 42.2 m/s | max|w| = 0.0793 m/s
[ Info: Iter  1250 | t = 10.417 days | Δt = 12 minutes | max|u| = 46.2 m/s | max|w| = 0.0799 m/s
[ Info: Iter  1300 | t = 10.833 days | Δt = 12 minutes | max|u| = 50.6 m/s | max|w| = 0.0763 m/s
[ Info: Iter  1350 | t = 11.250 days | Δt = 12 minutes | max|u| = 58.6 m/s | max|w| = 0.0942 m/s
[ Info: Iter  1400 | t = 11.667 days | Δt = 11.857 minutes | max|u| = 54.7 m/s | max|w| = 0.1945 m/s
[ Info: Iter  1450 | t = 12.057 days | Δt = 12 minutes | max|u| = 60.5 m/s | max|w| = 0.4303 m/s
[ Info: Iter  1500 | t = 12.464 days | Δt = 11.963 minutes | max|u| = 59.0 m/s | max|w| = 0.4036 m/s
[ Info: Iter  1550 | t = 12.867 days | Δt = 12 minutes | max|u| = 56.7 m/s | max|w| = 0.3376 m/s
[ Info: Iter  1600 | t = 13.265 days | Δt = 10.078 minutes | max|u| = 67.5 m/s | max|w| = 0.5707 m/s
[ Info: Iter  1650 | t = 13.576 days | Δt = 7.165 minutes | max|u| = 76.5 m/s | max|w| = 0.5169 m/s
[ Info: Iter  1700 | t = 13.816 days | Δt = 8.209 minutes | max|u| = 76.3 m/s | max|w| = 0.4243 m/s
[ Info: Iter  1750 | t = 14.118 days | Δt = 8.989 minutes | max|u| = 77.5 m/s | max|w| = 0.3549 m/s
[ Info: Iter  1800 | t = 14.414 days | Δt = 8.300 minutes | max|u| = 75.6 m/s | max|w| = 0.5849 m/s
[ Info: Iter  1850 | t = 14.692 days | Δt = 7.859 minutes | max|u| = 68.6 m/s | max|w| = 0.6692 m/s
[ Info: Iter  1900 | t = 14.938 days | Δt = 6.606 minutes | max|u| = 69.0 m/s | max|w| = 0.6253 m/s
[ Info: Iter  1950 | t = 15.162 days | Δt = 6.646 minutes | max|u| = 68.8 m/s | max|w| = 0.7339 m/s
[ Info: Iter  2000 | t = 15.404 days | Δt = 7.298 minutes | max|u| = 67.0 m/s | max|w| = 0.6190 m/s
[ Info: Iter  2050 | t = 15.667 days | Δt = 6.711 minutes | max|u| = 65.0 m/s | max|w| = 0.6883 m/s
[ Info: Iter  2100 | t = 15.892 days | Δt = 7.275 minutes | max|u| = 67.2 m/s | max|w| = 0.6490 m/s
[ Info: Iter  2150 | t = 16.112 days | Δt = 6.248 minutes | max|u| = 68.2 m/s | max|w| = 0.7260 m/s
[ Info: Iter  2200 | t = 16.332 days | Δt = 6.317 minutes | max|u| = 68.6 m/s | max|w| = 0.6041 m/s
[ Info: Iter  2250 | t = 16.599 days | Δt = 10.002 minutes | max|u| = 72.8 m/s | max|w| = 0.3873 m/s
[ Info: Iter  2300 | t = 16.917 days | Δt = 8.226 minutes | max|u| = 79.9 m/s | max|w| = 0.4547 m/s
[ Info: Iter  2350 | t = 17.220 days | Δt = 9.481 minutes | max|u| = 76.0 m/s | max|w| = 0.4502 m/s
[ Info: Iter  2400 | t = 17.512 days | Δt = 8.395 minutes | max|u| = 75.1 m/s | max|w| = 0.7244 m/s
[ Info: Iter  2450 | t = 17.824 days | Δt = 9.938 minutes | max|u| = 69.6 m/s | max|w| = 0.5769 m/s
[ Info: Iter  2500 | t = 18.173 days | Δt = 10.558 minutes | max|u| = 68.8 m/s | max|w| = 0.4572 m/s
[ Info: Iter  2550 | t = 18.512 days | Δt = 9.298 minutes | max|u| = 72.4 m/s | max|w| = 0.3501 m/s
[ Info: Iter  2600 | t = 18.840 days | Δt = 8.172 minutes | max|u| = 70.2 m/s | max|w| = 0.4475 m/s
[ Info: Iter  2650 | t = 19.154 days | Δt = 7.575 minutes | max|u| = 73.9 m/s | max|w| = 0.3454 m/s
[ Info: Iter  2700 | t = 19.430 days | Δt = 7.718 minutes | max|u| = 65.8 m/s | max|w| = 0.2722 m/s
[ Info: Iter  2750 | t = 19.735 days | Δt = 10.247 minutes | max|u| = 71.6 m/s | max|w| = 0.3794 m/s
[ Info: Iter  2800 | t = 20.074 days | Δt = 10.990 minutes | max|u| = 78.6 m/s | max|w| = 0.3345 m/s
[ Info: Iter  2850 | t = 20.473 days | Δt = 12 minutes | max|u| = 81.0 m/s | max|w| = 0.3952 m/s
[ Info: Iter  2900 | t = 20.883 days | Δt = 12 minutes | max|u| = 76.3 m/s | max|w| = 0.6826 m/s
[ Info: Iter  2950 | t = 21.300 days | Δt = 12 minutes | max|u| = 79.6 m/s | max|w| = 0.3814 m/s
[ Info: Iter  3000 | t = 21.717 days | Δt = 12 minutes | max|u| = 72.4 m/s | max|w| = 0.2656 m/s
[ Info: Iter  3050 | t = 22.133 days | Δt = 12 minutes | max|u| = 69.7 m/s | max|w| = 0.2799 m/s
[ Info: Iter  3100 | t = 22.550 days | Δt = 12 minutes | max|u| = 67.9 m/s | max|w| = 0.4474 m/s
[ Info: Iter  3150 | t = 22.967 days | Δt = 12 minutes | max|u| = 68.0 m/s | max|w| = 0.2343 m/s
[ Info: Iter  3200 | t = 23.383 days | Δt = 12 minutes | max|u| = 69.2 m/s | max|w| = 0.4405 m/s
[ Info: Iter  3250 | t = 23.766 days | Δt = 11.554 minutes | max|u| = 71.2 m/s | max|w| = 0.4370 m/s
[ Info: Iter  3300 | t = 24.163 days | Δt = 11.827 minutes | max|u| = 74.2 m/s | max|w| = 0.5726 m/s
[ Info: Iter  3350 | t = 24.567 days | Δt = 12 minutes | max|u| = 76.2 m/s | max|w| = 0.3285 m/s
[ Info: Iter  3400 | t = 24.983 days | Δt = 12 minutes | max|u| = 73.5 m/s | max|w| = 0.3961 m/s
[ Info: Iter  3450 | t = 25.400 days | Δt = 12 minutes | max|u| = 68.9 m/s | max|w| = 0.4526 m/s
[ Info: Iter  3500 | t = 25.786 days | Δt = 10.463 minutes | max|u| = 65.8 m/s | max|w| = 0.4147 m/s
[ Info: Iter  3550 | t = 26.144 days | Δt = 10.674 minutes | max|u| = 66.3 m/s | max|w| = 0.4864 m/s
[ Info: Iter  3600 | t = 26.524 days | Δt = 11.391 minutes | max|u| = 67.6 m/s | max|w| = 0.5367 m/s
[ Info: Iter  3650 | t = 26.919 days | Δt = 11.551 minutes | max|u| = 68.1 m/s | max|w| = 0.3842 m/s
[ Info: Iter  3700 | t = 27.304 days | Δt = 10.525 minutes | max|u| = 71.1 m/s | max|w| = 0.3726 m/s
[ Info: Iter  3750 | t = 27.615 days | Δt = 9.413 minutes | max|u| = 72.8 m/s | max|w| = 0.7580 m/s
[ Info: Iter  3800 | t = 27.898 days | Δt = 8.219 minutes | max|u| = 76.8 m/s | max|w| = 0.6097 m/s
[ Info: Iter  3850 | t = 28.173 days | Δt = 6.867 minutes | max|u| = 78.2 m/s | max|w| = 0.5544 m/s
[ Info: Iter  3900 | t = 28.402 days | Δt = 6.193 minutes | max|u| = 76.2 m/s | max|w| = 0.7781 m/s
[ Info: Iter  3950 | t = 28.627 days | Δt = 7.193 minutes | max|u| = 73.7 m/s | max|w| = 0.6858 m/s
[ Info: Iter  4000 | t = 28.899 days | Δt = 9.408 minutes | max|u| = 72.1 m/s | max|w| = 0.6565 m/s
[ Info: Iter  4050 | t = 29.215 days | Δt = 9.331 minutes | max|u| = 72.1 m/s | max|w| = 0.9330 m/s
[ Info: Iter  4100 | t = 29.539 days | Δt = 7.844 minutes | max|u| = 71.8 m/s | max|w| = 0.4266 m/s
[ Info: Iter  4150 | t = 29.814 days | Δt = 7.066 minutes | max|u| = 72.0 m/s | max|w| = 0.5222 m/s
[ Info: Simulation is stopping after running for 17.511 minutes.
[ Info: Simulation time 30 days equals or exceeds stop time 30 days.

Visualization

We plot three near-surface diagnostics on the sphere: the surface potential temperature $θ_{\rm sfc}$ (the classic diagnostic for the cold/warm sectors), the surface vertical vorticity $ζ$ (which reveals the cyclones and anticyclones), and the lower-tropospheric vertical velocity $w$ at ~2.9 km (the warm conveyor belt).

θ_ts = FieldTimeSeries("baroclinic_wave_k1.jld2",  "θ")ζ_ts = FieldTimeSeries("baroclinic_wave_k1.jld2",  "ζ")w_ts = FieldTimeSeries("baroclinic_wave_k16.jld2", "w")times = θ_ts.timesNt = length(times)k_sfc = 1k_mid = 16

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

sphere_kw = (elevation = π/6, azimuth = π/2, aspect = :data)ζlim = 1e-4wlim = 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 2.9 km", 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)endCairoMakie.record(fig, "baroclinic_wave.mp4", 1:Nt; framerate = 12, compression = 23) do nn    n[] = nnend

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfoversioninfo()

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 PkgPkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
⌃ [86bc3604] AtmosphericProfilesLibrary v0.1.8
  [660aa2fb] Breeze v0.8.0 `.`
  [052768ef] CUDA v6.2.1
  [13f3f980] CairoMakie v0.15.13
⌅ [6a9e3e04] CloudMicrophysics v0.37.1
  [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.4.1
  [b6400b83] DocumenterCodeBlocks v1.1.0
  [7da242da] Enzyme v0.13.198
  [46192b85] GPUArraysCore v0.2.0
  [63c18a36] KernelAbstractions v0.9.42
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.110.14
⌅ [a01a1ee8] RRTMGP v0.21.9
  [3c362404] Reactant v0.2.278
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.1
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.