Nonlinear primitive-equation simulations validating and extending Singh, Buckingham & Tandon (2025), Dynamics of Atmospheres and Oceans.
baroclinic_curvature/
├── Project.toml # dependencies
├── src/BaroclinicCurvature.jl # module: grid, base flow, model, output
├── run_straight.jl # Step 1 driver — straight Charney front
├── run_meander.jl # Step 2 driver — sinusoidal meander
├── tinker.jl # parameter playground for REPL workflow
└── plot_results.jl # postprocessing: dispersion, wb, movies
From this directory:
julia --project=. -e 'using Pkg; Pkg.instantiate()'First instantiate will take a while (CUDA, Oceananigans, Makie all compile).
Don't re-launch Julia for each run — startup is slow. Keep one REPL open and
use tinker.jl as the parameter playground.
julia> using Pkg; Pkg.activate(".")
julia> using Revise # picks up edits to src/ without restart
julia> include("tinker.jl") # runs whatever's in the EDIT block
# Edit the params block in tinker.jl in your editor, then:
julia> include("tinker.jl") # re-runs with new parametersOr just construct parameters and run directly at the prompt:
julia> p = ProblemParams(R_c = 80e3, Nx = 64, stop_time_days = 10,
output_prefix = "test1")
julia> sim = build_simulation(p; arch = GPU())
julia> run!(sim)A few convenience functions are defined in tinker.jl:
julia> smoke_test() # ~10 min sanity check, tiny grid
julia> run_experiment(R_c = 100e3, Nx = 96) # any kwargs override defaults
julia> sweep_Rc([300, 150, 100, 80]) # sequence of meander casesAll parameters with defaults are listed in ProblemParams inside
src/BaroclinicCurvature.jl. The most useful knobs:
| Knob | Effect |
|---|---|
Nx, Ny, Nz |
grid resolution (linear in time cost) |
stop_time_days |
how long to integrate |
R_c |
Inf for straight, finite (m) for meander |
L_m |
meander wavelength (m); keep ≫ Charney λ for scale separation |
Λs, N² |
shear and stratification — change the dominant λ and σ |
f₀, β |
latitude (changes Charney depth scale) |
perturbation_amplitude |
bump up if energy doesn't grow cleanly |
adjust_phase_days |
longer if the meander adjustment looks unfinished |
Tip: Always do smoke_test() first after any change to the module file —
catches API errors in ~10 minutes instead of after a 4-hour run.
Reproduces classical Charney results. Target: most-unstable wavelength ~180 km, amplitude growth rate ~0.025 day⁻¹.
julia --project=. run_straight.jl
julia --project=. plot_results.jl straightOutputs:
straight_snapshots.jld2,straight_wbprofile.jld2,straight_midslice.jld2straight_energy.pdf— perturbation KE vs time, log axis. Should be a clean straight line between roughly t=5 and t=25 days.straight_dispersion.pdf— σ(λ). Peak should be near 180 km at 0.025/day.straight_wb.pdf—⟨w'b'⟩(z)profile. Compare to Fig. 8 dashed line.straight_surface_b.mp4— surface buoyancy movie.
Meandering front with R_c set at the crests. Sweep R_c by running multiple cases:
for Rc in 300 150 100 80; do
julia --project=. run_meander.jl $Rc
julia --project=. plot_results.jl meander_Rc${Rc}km
doneThe first 5 days are an adjustment phase with strong Rayleigh damping (tapers quadratically). Treat any analysis as starting at t > 10 days.
| Quantity | Value |
|---|---|
| f₀ | 10⁻⁴ s⁻¹ (45°N) |
| β | 1.62 × 10⁻¹¹ m⁻¹ s⁻¹ |
| N² | 10⁻⁴ s⁻² (N/f = 100) |
| Λs | 10⁻⁴ s⁻¹ (vertical shear) |
| Lz | 2000 m |
Straight front: Lx = 1000 km, Ly = 400 km, 128 × 64 × 64. Meander: Lx = 1200 km (2·L_m), Ly = 500 km, 256 × 96 × 64.
Meander geometry: y_f(x) = A sin(2πx/L_m), L_m = 600 km, A = L_m²/(4π²·R_c). At R_c = 80 km, A ≈ 121 km, crest slope dy/dx ≈ 1.27.
Default architecture is GPU; falls back to CPU if CUDA isn't functional. On an RTX 3060 Ti (8 GB), the meander grid (256 × 96 × 64 ≈ 1.6 M cells) fits comfortably. Expect ~3–6 hours wall time for a 70-day meander run.
- Grid coordinate keys in JLD2 —
plot_results.jlusesxᶜᵃᵃetc; if your Oceananigans version uses different keys, inspect withjldopen(fn, "r") |> keys. The fallback is hard-codedrange(...). - No gradient-wind balance at init. The meander base flow is geostrophic
only. The adjustment phase (5 days, exponential damping with quadratic
taper) is meant to absorb the unbalanced centrifugal term. For a clean
gradient-wind initialization, replace
build_forcingwith an iterative pressure solver before model construction. - WENO advection of buoyancy. Default is
WENO(order=5)on momentum and centered on tracers via the model defaults. For sharper fronts, settracer_advection = WENO(order=5)explicitly. - Output is full 3D snapshots every 6 hours. Storage budget for the
default meander run is ~3–5 GB. Reduce
output_interval_hoursor subset fields if needed.