From a9b7f68d159167c9b472d36a87ed2565db215f5a Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Sun, 9 Nov 2025 07:53:19 -0700 Subject: [PATCH 1/2] draw up a basic script --- examples/tropical_cyclone_world.jl | 321 +++++++++++++++++++++++++++++ 1 file changed, 321 insertions(+) create mode 100644 examples/tropical_cyclone_world.jl diff --git a/examples/tropical_cyclone_world.jl b/examples/tropical_cyclone_world.jl new file mode 100644 index 00000000..63657709 --- /dev/null +++ b/examples/tropical_cyclone_world.jl @@ -0,0 +1,321 @@ +# # Tropical cyclone world (Cronin & Chavas, 2019) +# +# This example follows the "tropical cyclone world" radiative–convective equilibrium +# experiment described by [Cronin & Chavas (2019)](https://doi.org/10.1175/JAS-D-18-0357.1), +# but implemented with `Breeze.AtmosphereModel`. The goal is to couple a simple +# surface-layer bulk scheme (wind stress, sensible, and latent heat fluxes) to an +# anelastic atmosphere, apply a horizontally-uniform radiative cooling, and allow the +# model to organize convection and tropical cyclones on a doubly-periodic domain. +# +# The script is written in `Literate.jl` style: prose cells introduce each block of +# code, the physics it represents, and how it connects to the Cronin & Chavas setup. +# The numerical values below are close to those used in the paper, but whenever +# possible we pull thermodynamic constants directly from `Breeze.ThermodynamicConstants` +# rather than hard-coding numbers. + +# ## Packages and utilities +# +# We rely on `Breeze` for the model, `Oceananigans` utilities for grids, `Random` +# for small perturbations, and `CairoMakie` for quick-look visualizations. + +using Breeze +using Breeze.Thermodynamics: + ThermodynamicConstants, + ReferenceState, + MoistureMassFractions, + mixture_heat_capacity, + saturation_specific_humidity + +using Breeze.AtmosphereModels: AtmosphereModel +using Oceananigans.BoundaryConditions: fill_halo_regions! +using Printf +using Random +using CairoMakie + +# ## Domain, grid, and reference state +# +# Cronin & Chavas integrate over a square, doubly-periodic domain that spans a few +# thousand kilometers and reaches the lower stratosphere. We adopt a 4 096 × 4 096 km +# horizontal box with a 25 km model top. The "modest" resolution requested corresponds +# to 128 × 128 × 32 cells (≈32 km horizontally, ≈780 m vertically). For CI or +# documentation builds we automatically reduce the resolution to keep runtimes low. + +arch = CPU() # swap to GPU() when accelerators are available +Nx = Ny = 128 +Nz = 32 + +if get(ENV, "CI", "false") == "true" + Nx = Ny = 16 + Nz = 8 +end + +Lx = 4_096_000.0 # meters +Ly = 4_096_000.0 +H = 25_000.0 +halo = (5, 5, 5) + +grid = RectilinearGrid(arch; + size = (Nx, Ny, Nz), + x = (0, Lx), + y = (0, Ly), + z = (0, H), + halo = halo, + topology = (Periodic, Periodic, Bounded)) + +reference_state = ReferenceState(grid; base_pressure=101325, potential_temperature=300) +formulation = AnelasticFormulation(reference_state) + +# ## Surface exchange physics: drag, sensible, and latent heat/moisture fluxes +# +# The Cronin & Chavas case uses a bulk aerodynamic surface scheme. The exchange +# coefficients below (drag `Cd`, sensible `Ch`, and moisture `Ce`) are close to their +# recommended values for a TC-permitting LES. The functions return boundary fluxes for +# the prognostic variables (`ρu`, `ρv`, `ρe`, `ρqᵗ`) and follow Oceananigans' sign +# convention: positive fluxes at the bottom boundary add that quantity to the interior. + +# Compute surface specific humidity and density +# TODO: move this to source code +thermo = reference_state.thermodynamics +T₀ = 302 +p₀ = reference_state.base_pressure +pᵛ⁺₀ = saturation_vapor_pressure(T₀, thermo, thermo.liquid) +Rᵛ = vapor_gas_constant(thermo) +Rᵈ = dry_air_gas_constant(thermo) +ϵᵛᵈ = Rᵛ / Rᵈ - 1 +qᵛ⁺₀ = pᵛ⁺₀ / (p₀ - ϵᵛᵈ * pᵛ⁺₀) +q₀ = MoistureMassFractions(qᵛ⁺₀, zero(qᵛ⁺₀), zero(qᵛ⁺₀)) +ρ₀ = density(p₀, T₀, q₀, thermo) + +surface_exchange = (; ρ₀, T₀, qᵛ⁺₀, + Cd = 1.5e-3, + Ch = 1.1e-3, + Ce = 1.1e-3, + 𝕌₀ = 0.5, + ℒˡᵣ = thermo.liquid.reference_latent_heat, + cᵖᵈ = thermo.dry_air.heat_capacity, + thermo) + +@inline function surface_speed(ρu, ρv, params) + u = ρu / params.ρₛ + v = ρv / params.ρₛ + s = sqrt(u^2 + v^2 + params.𝕌₀^2) + return s, u, v +end + +@inline function drag_flux_x(x, y, t, ρu, ρv, params) + s, u, _ = surface_speed(ρu, ρv, params) + τx = params.ρ₀ * params.Cd * s * u + return -τx +end + +@inline function drag_flux_y(x, y, t, ρu, ρv, params) + s, _, v = surface_speed(ρu, ρv, params) + τy = params.ρ₀ * params.Cd * s * v + return -τy +end + +@inline function surface_moisture_flux(x, y, t, ρqᵗ, ρu, ρv, params) + s, _, _ = surface_speed(ρu, ρv, params) + qᵗ = ρqᵗ / params.ρ₀ + Δq = max(0, params.qᵛ⁺₀ - qᵗ) + return params.ρ₀ * params.Ce * s * (params.qᵛ⁺₀ - qᵗ) +end + +@inline function surface_energy_flux(x, y, t, ρe, ρqᵗ, ρu, ρv, params) + s, _, _ = surface_speed(ρu, ρv, params) + qᵗ = ρqᵗ / params.ρ₀ + q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ)) + cᵖᵐ = mixture_heat_capacity(q, params.thermo) + # ρe = ρᵣ (cᵖᵐ T + qᵗ Lᵥ) at z ≈ 0 ⇒ T = (ρe/ρᵣ - qᵗ Lᵥ) / cᵖᵐ + Tₐ = (ρe / params.ρ₀ - qᵗ * params.ℒˡᵣ) / cᵖᵐ + sensible = params.ρ₀ * params.cᵖᵈ * params.Ch * s * (params.T₀ - Tₐ) + moisture_flux = surface_moisture_flux(x, y, t, ρqᵗ, ρu, ρv, params) + latent = params.ℒˡᵣ * moisture_flux + return sensible + latent +end + +ρu_bcs = FieldBoundaryConditions(bottom = FluxBoundaryCondition(drag_flux_x, + field_dependencies = (:ρu, :ρv), + parameters = surface_exchange)) +ρv_bcs = FieldBoundaryConditions(bottom = FluxBoundaryCondition(drag_flux_y, + field_dependencies = (:ρu, :ρv), + parameters = surface_exchange)) +ρq_bcs = FieldBoundaryConditions(bottom = FluxBoundaryCondition(surface_moisture_flux, + field_dependencies = (:ρqᵗ, :ρu, :ρv), + parameters = surface_exchange)) +ρe_bcs = FieldBoundaryConditions(bottom = FluxBoundaryCondition(surface_energy_flux, + field_dependencies = (:ρe, :ρqᵗ, :ρu, :ρv), + parameters = surface_exchange)) + +boundary_conditions = (; ρu = ρu_bcs, + ρv = ρv_bcs, + ρqᵗ = ρq_bcs, + ρe = ρe_bcs) + +# ## Radiative forcing +# +# Cronin & Chavas impose a vertically varying radiative cooling (≈ -1.5 K day⁻¹ near +# cloud base that tapers with height). We encode this forcing as an energy tendency +# that will be added directly to the `ρe` equation after the model is constructed. + +using Oceananigans.Units: day + +dt_Tˢ = -1.5 / day +λᴿ = 12_000 +ρᵣ = reference_state.density +cᵖᵈ = thermo.dry_air.heat_capacity +dt_Tᴿ = CenterField(grid) +dt_Tᴿ_func = ρᵣ * cᵖᵈ * dt_Tˢ * exp(-z / λᴿ) + +set!(dt_Tᴿ, dt_Tᴿ_func) +radiative_forcing = Forcing(dt_Tᴿ_func) + +# ## Build the AtmosphereModel +# +# We use 9th-order WENO advection, an f-plane with mid-latitude Coriolis parameter, +# and no microphysics for now (warm-rain processes would be the next upgrade). + +coriolis = FPlane(f = 5e-5) +advection = WENO(order=9) + +model = AtmosphereModel(grid; + thermodynamics = thermo, + formulation, + boundary_conditions, + coriolis, + advection) + +model.forcing = merge(model.forcing, (; ρe = radiative_forcing)) + +# ## Initial conditions +# +# We start from a weakly stratified reference potential temperature profile, a moist +# boundary layer that decays exponentially with height, and small random velocity and +# moisture perturbations to seed convection. The helper below keeps the perturbations +# reproducible. + +Random.seed!(42) +N² = 1e-4 +g = thermo.gravitational_acceleration +dθdz = N² * θ_ref / g +q_surface = 0.018 +q_scale_height = 3_000.0 + +θ_profile(z) = θ_ref + dθdz * z +q_profile(z) = q_surface * exp(-z / q_scale_height) + +δθ = 0.5 +δq = 1e-4 +δu = 0.5 + +θᵢ(x, y, z) = θ_profile(z) + δθ * randn() +qᵢ(x, y, z) = max(q_profile(z) + δq * randn(), 0) +Ξᵢ(x, y, z) = δu * randn() + +set!(model, θ = θᵢ, qᵗ = qᵢ, u = Ξᵢ, v = Ξᵢ, w = (x, y, z) -> 0.0) + +# ## Simulation control and diagnostics +# +# We advance the system for ~100 time steps on CPU to validate the setup quickly. +# The CFL wizard keeps the model stable, and a simple callback prints thermodynamic +# and kinematic extrema so long runs can be monitored. + +Δt = 1.0 +stop_iteration = 100 + +if get(ENV, "CI", "false") == "true" + stop_iteration = 20 +end + +simulation = Simulation(model; Δt, stop_iteration) +conjure_time_step_wizard!(simulation, cfl = 0.6) + +ρe = model.energy +u, v, w = model.velocities + +function progress(sim) + ρemin = minimum(ρe) + ρemax = maximum(ρe) + u_max = maximum(abs, u) + v_max = maximum(abs, v) + w_max = maximum(abs, w) + + msg = @sprintf("Iter: %5d, t: %8.2fs, extrema(ρe) = (%.3e, %.3e) J/m³, max|u,v,w| = (%.2f, %.2f, %.2f) m/s", + iteration(sim), sim.model.clock.time, + ρemin, ρemax, u_max, v_max, w_max) + @info msg + return nothing +end + +add_callback!(simulation, progress, IterationInterval(10)) + +# ## Output writers +# +# We capture surface slices and mid-plane sections of vertical velocity, energy, and +# vertical vorticity for documentation plots. The file lives inside `examples/` so it +# can be inspected manually after the script finishes. + +ζ = ∂x(v) - ∂y(u) +outputs = (; u, v, w, ζ, ρe) +output_path = joinpath(@__DIR__, "tropical_cyclone_world_$(Nx)x$(Ny)x$(Nz).jld2") +schedule = IterationInterval(max(1, stop_iteration ÷ 10)) + +simulation.output_writers[:jld2] = JLD2Writer(model, outputs; + filename = output_path, + schedule, + overwrite_existing = true) + +@info "Starting tropical cyclone world spin-up on $(summary(grid))" +run!(simulation) +@info "Finished spin-up, output saved to $(output_path)" + +# ## Quick-look visualization +# +# We make horizontal (surface) and vertical (mid-plane) slices of the vorticity, +# vertical velocity, and moist static energy. The energy slice is a placeholder until +# a dedicated potential-temperature diagnostic is exposed. + +if get(ENV, "CI", "false") == "false" && isfile(output_path) + ζt = FieldTimeSeries(output_path, "ζ") + wt = FieldTimeSeries(output_path, "w") + ρet = FieldTimeSeries(output_path, "ρe") + + n = length(ζt) + ζ_xy = Array(interior(ζt[n], :, :, 1)) + w_xy = Array(interior(wt[n], :, :, 1)) + y_index = size(w_xy, 2) ÷ 2 + w_xz = Array(interior(wt[n], :, y_index, :)) + ρe_xz = Array(interior(ρet[n], :, y_index, :)) + + x = xnodes(ζt) ./ 1000 # km + y = ynodes(ζt) ./ 1000 + z = znodes(ρet) ./ 1000 + + fig = Figure(size = (1100, 900), fontsize = 12) + ax1 = Axis(fig[1, 1], title = "Surface ζ", xlabel = "x (km)", ylabel = "y (km)") + ax2 = Axis(fig[1, 2], title = "Surface w", xlabel = "x (km)", ylabel = "y (km)") + ax3 = Axis(fig[2, 1], title = "Mid-plane w", xlabel = "x (km)", ylabel = "z (km)") + ax4 = Axis(fig[2, 2], title = "Mid-plane ρe", xlabel = "x (km)", ylabel = "z (km)") + + ζ_lim = maximum(abs, ζ_xy) + w_lim_surface = maximum(abs, w_xy) + w_lim_mid = maximum(abs, w_xz) + ζ_map = heatmap!(ax1, x, y, ζ_xy'; + colormap = :balance, colorrange = (-ζ_lim, ζ_lim)) + w_surface = heatmap!(ax2, x, y, w_xy'; + colormap = :balance, colorrange = (-w_lim_surface, w_lim_surface)) + + w_mid_map = heatmap!(ax3, x, z, w_xz'; + colormap = :balance, colorrange = (-w_lim_mid, w_lim_mid)) + ρe_map = heatmap!(ax4, x, z, ρe_xz'; + colormap = :thermal) + + Colorbar(fig[1, 3], ζ_map, label = "ζ (s⁻¹)") + Colorbar(fig[1, 4], w_surface, label = "w (m s⁻¹)") + Colorbar(fig[2, 3], w_mid_map, label = "w (m s⁻¹)") + Colorbar(fig[2, 4], ρe_map, label = "ρe (J m⁻³)") + + png_path = replace(output_path, ".jld2" => "_slices.png") + save(png_path, fig) + @info "Saved slice visualization to $(png_path)" +end From 6659cee50681c1e4715f19e000829e57e33e3f64 Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Sun, 9 Nov 2025 20:55:11 -0700 Subject: [PATCH 2/2] Update examples/tropical_cyclone_world.jl Co-authored-by: Navid C. Constantinou --- examples/tropical_cyclone_world.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tropical_cyclone_world.jl b/examples/tropical_cyclone_world.jl index 63657709..d2c71cb7 100644 --- a/examples/tropical_cyclone_world.jl +++ b/examples/tropical_cyclone_world.jl @@ -1,7 +1,7 @@ # # Tropical cyclone world (Cronin & Chavas, 2019) # # This example follows the "tropical cyclone world" radiative–convective equilibrium -# experiment described by [Cronin & Chavas (2019)](https://doi.org/10.1175/JAS-D-18-0357.1), +# experiment described by [Cronin & Chavas (2019)](@citet Cronin2019), # but implemented with `Breeze.AtmosphereModel`. The goal is to couple a simple # surface-layer bulk scheme (wind stress, sensible, and latent heat fluxes) to an # anelastic atmosphere, apply a horizontally-uniform radiative cooling, and allow the