ERA5 → 12 km land-coupled, convection-permitting hindcast (Breeze + NestedModel + SlabLand)
A limited-area model (LAM) example that downscales ERA5 reanalysis to a ~12 km Breeze compressible atmosphere over the U.S. Southern Great Plains, for the Midlatitude Continental Convective Clouds Experiment (MC3E) 20 May 2011 squall-line case (Fan et al. (2017)) — with an interactive SlabLand surface under the child and all-sky RRTMGP radiation heating the column and the surface, so the diurnal cycle, the clouds, and the land talk to each other.
nested_atmosphere_model(grid, dataset; dates, …) builds the whole nest: an ERA5 "parent" PrescribedAtmosphere on its native 0.25° pressure-level grid, driving a Breeze "child" through open lateral boundary conditions + interior Davies relaxation, both derived on the fly from the parent's raw state. The constructor also initializes the child from the reanalysis and spins it into balance, so a plain Simulation(model) then steps the ready nest.
What this example does
- Downloads ERA5 (pressure + single levels) and regrids it onto a terrain-following
LatitudeLongitudeGrid(ETOPO2022 orography, blended to the parent's at the boundary). - Initializes the prognostic state from ERA5 with a hydrostatic-from-surface pressure, a terrain-consistent
w̃ ≈ 0, and a dynamical-initialization (DFI) pass that balancesρw. - Integrates the compressible equations with split-explicit acoustic substepping, 1-moment mixed-phase microphysics, Coriolis, and Rayleigh damping.
- Couples an interactive
SlabLand(prognostic skin temperature + bucket hydrology) under the child through Monin–Obukhov similarity fluxes, and the child's rain into the land bucket. - Runs all-sky RRTMGP radiative transfer inside the child — interior heating rates on the evolving clouds, observed (CGLS 1 km) surface albedo, surface fluxes into the land energy budget, and the interface skin temperature closing the longwave loop.
- Writes and animates horizontal slices, a vertical section, and the land-surface response.
What it does NOT do (yet)
- Single nest only (ERA5 → 12 km; coarsened 4× from Fan's 3 km Domain 3 for a fast configuration).
- The land is a slab: no soil column, no vegetation, and no boundary-layer or cumulus parameterization — diffusion is numerical, and deep convection is resolved on the grid.
- Snow does not reach the bucket yet (the coupler diagnoses the child's surface rain flux, but no snow analog); immaterial for this warm-season case.
using NumericalEarth
using Oceananigans
using Oceananigans.Units # `minutes` for the output schedule (not re-exported by Oceananigans)
using Breeze
using CopernicusClimateDataStore # activates NumericalEarthCopernicusClimateDataStoreExt (ERA5 + land-albedo downloads)
using CloudMicrophysics # nested_atmosphere_model's default microphysics → 1-moment mixed-phase
using RRTMGP # activates Breeze's radiative-transfer extension (RadiativeTransferModel)
using CairoMakie # loads Makie → NumericalEarthMakieExt (`visualize_nested_domain`)
using NaturalEarth, GeoInterface # together → NumericalEarthNaturalEarthExt (`natural_earth_lines`)
using CUDA
using Printf
using Dates: DateTime, Second[ Info: Oceananigans will use 8 threads
This 12 km LAM (150×136×50 ≈ 1.0M cells, split-explicit) targets a CUDA GPU; switch to CPU() only for a small smoke test.
arch = GPU()Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}(CUDACore.CUDAKernels.CUDABackend(false, true))Set Oceananigans' global default float type, cascading to all grids, fields, and constants.
Oceananigans.defaults.FloatType = Float32Float32Configuration
Domain centered on the DOE Atmospheric Radiation Measurement (ARM) Southern Great Plains (SGP) site in Lamont, OK — the 3 km Domain 3 of the WRF telescoping nest in Fan et al. (2017), coarsened 4× and driven directly by ERA5.
# dates
name = "mc3e"
duration = 3hours
start_date = DateTime(2011, 05, 20, 0)
stop_date = start_date + Second(duration)
# location
φ₀, λ₀ = 36.605, -97.485 # center latitude, longitude (deg)
Lλ, Lφ = 16.7, 15.1
# horizontal resolution: grid spacing in degrees; the cell counts follow from the extent
# (which stays fixed, so the realized spacing Lλ/Nx only approximates Δλ)
Δλ = Δφ = 1/9 # ≈ 12 km in latitude (ERA5's native 0.25° → 2.25× refinement)
Nx = round(Int, Lλ / Δλ)
Ny = round(Int, Lφ / Δφ)
dates = (start_date, stop_date)
λ_west, λ_east = λ₀ .+ [-1, 1] .* Lλ / 2
φ_south, φ_north = φ₀ .+ [-1, 1] .* Lφ / 22-element Vector{Float64}:
29.054999999999996
44.154999999999994Vertical grid matched to Fan et al. (2017)'s WRF nest: Nz = 50 cells, a constant 60 m surface cell, 490 m maximum spacing, and a model top at ~20 km (~50 hPa).
z = ReferenceToStretchedDiscretization(extent = 19525.0,
bias = :left,
bias_edge = 0,
constant_spacing = 60,
constant_spacing_extent = 60,
maximum_spacing = 490,
stretching = LinearStretching(0.15))
Nz = length(z)50LAM grid
A bounded terrain-following LatitudeLongitudeGrid. The terrain itself is handled by nested_atmosphere_model (its terrain keyword): ETOPO2022 is regridded onto the grid, blended toward the parent's ERA5 orography over the outermost relax_width cells (so the terrain at the open boundaries matches the parent state driving them), and materialized onto the grid's terrain-following coordinate (Gal-Chen–Somerville σ): the bottom follows the ground, the top stays flat.
grid = LatitudeLongitudeGrid(arch;
longitude = (λ_west, λ_east),
latitude = (φ_south, φ_north),
z = TerrainFollowingVerticalDiscretization(z),
size = (Nx, Ny, Nz),
halo = (5, 5, 5),
topology = (Bounded, Bounded, Bounded))150×136×50 LatitudeLongitudeGrid{Float32, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on CUDAGPU with 5×5×5 halo
├── longitude: Bounded λ ∈ [-105.835, -89.135] regularly spaced with Δλ=0.111333
├── latitude: Bounded φ ∈ [29.055, 44.155] regularly spaced with Δφ=0.111029
└── z: Bounded z ∈ [0.0, 20003.0] variably spaced with terrain-following min(Δr)=60.0, max(Δr)=490.0Build the nest
nested_atmosphere_model(grid, dataset; dates, …) derives the parent region from the child grid (padded by two native ERA5 cells), loads the parent PrescribedAtmosphere on ERA5's native pressure-level grid, anchors the default split-explicit compressible dynamics at the domain-mean ERA5 surface pressure, derives the parent-driven lateral BCs + Davies relaxation (cosine ramp over relaxation_width cells), materializes the blended terrain, and wraps parent + child in a NestedModel whose time_step! advances the child then ticks the parent clock.
Time window: a (start, end) tuple — the hourly forcing cadence is the dataset's own. The initial outer time step is gentle (the adaptive wizard below ramps it). The Davies relaxation is an explicit nudge dφ/dt = −r(φ − φₚₐᵣₑₙₜ), stable only for r·Δt ≲ 2, so its rate is a fixed 1/300 s⁻¹ — small enough that r·Δt stays well below 1 even at the wizard's max_Δt, rather than tied to Δt (which would climb into the over-relaxation regime as Δt grows). Momentum uses a narrower WENO(5) stencil, which interacts less violently with the open-boundary halo than the default WENO(9).
era5_datadir = "era5" # where ERA5 data is saved locally
dataset = ERA5HourlyPressureLevels()
relax_width = 5
nest = nested_atmosphere_model(grid, dataset;
dates,
dir = era5_datadir,
terrain = ETOPO2022(),
terrain_blend_width = relax_width,
relaxation_rate = 1/300,
relaxation_width = relax_width,
momentum_advection = WENO(order = 5))NestedModel(73×67×37×4 PrescribedAtmosphere{Float32} (λ ∈ [-106.625, -88.375], φ ∈ [28.25, 45.0]) → AtmosphereModel{GPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0))
├── parent: 73×67×37×4 PrescribedAtmosphere{Float32} (λ ∈ [-106.625, -88.375], φ ∈ [28.25, 45.0])
└── child: AtmosphereModel{GPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)No bottom_drag_coefficient: the surface stress (with the heat and moisture fluxes) comes from the land coupling below, so the child's bottom boundary conditions stay the coupler's flux fields.
The realized parent region (child + padding, snapped to the native 0.25° grid) serves the domain map and the ERA5 snapshots below.
parent = nest.parent
era5_region = BoundingBox(parent.grid)NumericalEarth.DataWrangling.BoundingBox{Tuple{Float32, Float32}, Tuple{Float32, Float32}, Nothing}((-106.625f0, -88.375f0), (28.25f0, 45.0f0), nothing)Nested domains
Visualize the nesting before stepping the model — the ERA5 forcing region and the 12 km LAM child, over ETOPO relief with Natural Earth state/country boundaries — so the domain geometry is written even if the run is cut short.
fig = visualize_nested_domain(grid;
parent = era5_region,
padding = 2.5,
title = "ERA5 → 12 km LAM nest (MC3E squall line, ARM SGP)",
label = "12 km LAM (child)",
parent_label = "ERA5 parent",
landmarks = tuple("ARM SGP" => (λ₀, φ₀)))
save("era5_breeze_domains.png", fig)
fig
ERA5 pressure velocity ω (Pa/s), loaded raw onto the parent's own grid — so it shares the parent's time-varying geopotential vertical and composes in the w ≈ −ω/(ρg) operation with the parent state.
ω_metadata = Metadata(:vertical_velocity; dataset, dates, region = era5_region, dir = era5_datadir)
ω_series = FieldTimeSeries(ω_metadata, parent.grid; time_indices_in_memory = length(ω_metadata.dates))73×67×37×4 FieldTimeSeries{NumericalEarth.DataWrangling.DatasetBackend} located at (Center, Center, Center) on Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}
├── grid: 73×67×37 LatitudeLongitudeGrid{Float32, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on CUDAGPU with 3×3×3 halo
├── indices: (:, :, :)
├── time_indexing: Cyclical(period=14400.0)
├── backend: DatasetBackend(1, 4)
└── data: 79×73×43×4 OffsetArray(::CUDACore.CuArray{Float32, 4, CUDACore.DeviceMemory}, -2:76, -2:70, -2:40, 1:4) with eltype Float32 with indices -2:76×-2:70×-2:40×1:4
└── max=3.05137, min=-6.87546, mean=-0.0348586Land surface
The lower boundary is interactive: a SlabLand on the child's horizontal grid — prognostic skin temperature plus bucket hydrology, driven by the coupled turbulent fluxes, the child's rain, and the radiation below. Temperature initializes from ERA5's skin temperature; the bucket starts with ≈ 20 mm of water.
land_grid = LatitudeLongitudeGrid(arch;
longitude = (λ_west, λ_east),
latitude = (φ_south, φ_north),
size = (Nx, Ny),
halo = (5, 5),
topology = (Bounded, Bounded, Flat))
land = SlabLand(land_grid)
skin_temperature = Metadatum(:skin_temperature; dataset = ERA5HourlySingleLevel(),
date = start_date, region = era5_region, dir = era5_datadir)
set!(land.temperature, skin_temperature)
set!(land; M = 20)SlabLand{Float32, GPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
├── grid: 150×136×1 LatitudeLongitudeGrid{Float32, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Flat} on CUDAGPU with 5×5×0 halo
├── energy: ForceRestoreEnergy(dry_heat_capacity=222000.0, liquid_heat_capacity=4186.0, deep_temperature=0.0, deep_time_scale=Inf)
├── hydrology: BucketHydrology(maximum_water_storage=150.0)
├── temperature: 150×136×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.LatitudeLongitudeGrid on CUDAGPU
├── water_storage: 150×136×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.LatitudeLongitudeGrid on CUDAGPU
├── saturation: 150×136×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.LatitudeLongitudeGrid on CUDAGPU
├── fluxes: (:surface_energy_flux, :precipitation, :evaporation)
└── diagnostics: ()RRTMGP radiation
All-sky RRTMGP: interior heating on the evolving clouds, surface fluxes into the land budget. The surface albedo is observed, not constant — passing the dataset materializes the CGLS 1 km blue-sky albedo at the ten-day date nearest epoch, with water pixels (NaN in the land product) defaulting to open water. Everything else is a Breeze default — per-column sun angles from the grid, climatological ozone, emissivity, effective radii. No surface temperature either: the coupled model below binds the atmosphere–land interface skin temperature, re-read every solve. Hourly solves match the ERA5 cadence.
radiation = RadiativeTransferModel(grid, AllSkyOptics(), nest.child.thermodynamic_constants;
solar_position = ApparentSolarPosition(epoch = start_date),
surface_albedo = CopernicusAlbedo(),
schedule = TimeInterval(1hour))RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── solar_position: ApparentSolarPosition(coordinate=<from grid>, epoch=2011-05-20T00:00:00)
├── surface_temperature: Nothing
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: 150×136×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.LatitudeLongitudeGrid on CUDAGPU
├── grid: 150×136×50 LatitudeLongitudeGrid{Float32, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on CUDAGPU with 5×5×5 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: ZeroFlux, east: ZeroFlux, south: ZeroFlux, north: ZeroFlux, bottom: Nothing, top: Nothing, immersed: Nothing
└── data: 160×146×1 OffsetArray(::CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}, -4:155, -4:141, 1:1) with eltype Float32 with indices -4:155×-4:141×1:1
└── max=0.451573, min=0.0377306, mean=0.184505
├── liquid_effective_radius: Breeze.AtmosphereModels.ConstantRadiusParticles{Float32}(1.0f-5)
├── ice_effective_radius: Breeze.AtmosphereModels.ConstantRadiusParticles{Float32}(3.0f-5)
└── diffuse_surface_albedo: 150×136×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.LatitudeLongitudeGrid on CUDAGPU
├── grid: 150×136×50 LatitudeLongitudeGrid{Float32, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded, Oceananigans.Grids.Bounded} on CUDAGPU with 5×5×5 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: ZeroFlux, east: ZeroFlux, south: ZeroFlux, north: ZeroFlux, bottom: Nothing, top: Nothing, immersed: Nothing
└── data: 160×146×1 OffsetArray(::CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}, -4:155, -4:141, 1:1) with eltype Float32 with indices -4:155×-4:141×1:1
└── max=0.451573, min=0.0377306, mean=0.184505Coupled model
AtmosphereLandModel wires nest, land, and radiation together: it binds the interface skin temperature into the RTM, materializes the child's radiation proxy, and each step computes Monin–Obukhov fluxes from the child's lowest-cell state, writing them into the child's bottom flux boundary conditions and the land's accumulators. The adaptive outer Δt is bounded by the advective CFL (acoustic modes are substepped).
Δt = 10
atmosphere = Simulation(nest; Δt) ## the coupled model manages Δt; this sets only the initial value
model = AtmosphereLandModel(atmosphere, land; radiation)
simulation = Simulation(model; Δt, stop_time=duration)
conjure_time_step_wizard!(simulation, IterationInterval(3); cfl=0.5, max_Δt=Δt)┌ Warning: This simulation will run forever as stop iteration = stop time = wall time limit = Inf.
└ @ Oceananigans.Simulations /storage5/github-action-runners/julia-depot/packages/Oceananigans/lCP9j/src/Simulations/simulation.jl:77
┌ Warning: This simulation will run forever as stop iteration = stop time = wall time limit = Inf.
└ @ Oceananigans.Simulations /storage5/github-action-runners/julia-depot/packages/Oceananigans/lCP9j/src/Simulations/simulation.jl:77
Output
Three JLD2Writers save 2-D slices of one shared fields NamedTuple of online diagnostics — never the 3-D fields, which would overflow the disk. Each output is an AbstractOperation (Breeze's diagnostic accessors) that every writer computes and slices at save time via its own indices: horizontal slices at the surface and at k_aloft (the reference level nearest 2 km — on the terrain-following grid a constant reference level ≈ constant height above ground near the surface), and a zonal x-z section at j_section, the latitude row through the ARM SGP site. A fourth writer saves the land surface: skin temperature and bucket saturation (2-D already, no slicing needed — the writer stores each output's own grid).
child = nest.child
k_aloft = searchsortedfirst(Array(znodes(grid, Center())), 2000)
φ_centers = Array(φnodes(grid, Center(), Center(), Center()))
j_section = searchsortedfirst(φ_centers, φ₀)
fields = (θᵛ = VirtualPotentialTemperature(child),
U = sqrt(child.velocities.u^2 + child.velocities.v^2),
w = child.velocities.w,
qᵛ = specific_humidity(child),
qʳ = child.microphysical_fields.qʳ)
surface_filename = name * "_surface.jld2"
aloft_filename = name * "_aloft.jld2"
section_filename = name * "_section.jld2"
schedule = TimeInterval(20minutes)
slice_writer(indices, filename) = JLD2Writer(model, fields; schedule, filename, indices,
overwrite_files = true)
simulation.output_writers[:surface] = slice_writer((:, :, 1), surface_filename)
simulation.output_writers[:aloft] = slice_writer((:, :, k_aloft), aloft_filename)
simulation.output_writers[:section] = slice_writer((:, j_section, :), section_filename)
land_filename = name * "_land.jld2"
land_fields = (Tˡᵃ = land.temperature, 𝒮 = land.saturation)
simulation.output_writers[:land] = JLD2Writer(model, land_fields; schedule,
filename = land_filename,
overwrite_files = true)
function progress(sim)
child = sim.model.atmosphere.model.child
u, v, w = child.velocities
ρ = child.dynamics.total_density
qᵛ = specific_humidity(child)
qʳ = child.microphysical_fields.qʳ
@info @sprintf("iter=%4d, t=%s, Δt=%s, max|u|=(%7.2f, %7.2f, %6.2f), ρ ∈ [%.4f, %.4f], qᵛ ∈ [%.4g, %.4g], qʳ ∈ [%.2g, %.2g]",
sim.model.clock.iteration, prettytime(sim), prettytime(sim.Δt),
maximum(abs, u), maximum(abs, v), maximum(abs, w), minimum(ρ), maximum(ρ),
minimum(qᵛ), maximum(qᵛ), minimum(qʳ), maximum(qʳ))
return nothing
end
add_callback!(simulation, progress, IterationInterval(100))Run
run!(simulation)[ Info: Initializing simulation...
[ Info: iter= 0, t=0 seconds, Δt=10 seconds, max|u|=( 41.81, 49.51, 0.65), ρ ∈ [0.0959, 1.1797], qᵛ ∈ [-3.306e-06, 0.01776], qʳ ∈ [0, 0]
[ Info: ... simulation initialization complete (2.671 minutes)
[ Info: Executing initial time step...
[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (20.638 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.139 minutes).
[ Info: ... initial time step complete (1.882 minutes).
[ Info: iter= 100, t=5.352 minutes, Δt=2.808 seconds, max|u|=( 41.73, 49.30, 15.84), ρ ∈ [0.0950, 1.1871], qᵛ ∈ [-3.518e-06, 0.02312], qʳ ∈ [-9.9e-10, 0.00011]
[ Info: iter= 200, t=9.890 minutes, Δt=2.308 seconds, max|u|=( 41.65, 49.56, 19.12), ρ ∈ [0.0949, 1.1867], qᵛ ∈ [-2.873e-06, 0.02339], qʳ ∈ [-1e-09, 0.00012]
[ Info: iter= 300, t=14.041 minutes, Δt=2.003 seconds, max|u|=( 41.58, 49.66, 19.06), ρ ∈ [0.0949, 1.1870], qᵛ ∈ [-2.167e-06, 0.02344], qʳ ∈ [-1.9e-10, 0.00019]
[ Info: iter= 400, t=17.884 minutes, Δt=1.960 seconds, max|u|=( 43.83, 49.63, 16.50), ρ ∈ [0.0948, 1.1870], qᵛ ∈ [-1.507e-06, 0.02346], qʳ ∈ [-1.1e-10, 0.00023]
[ Info: iter= 500, t=21.359 minutes, Δt=1.888 seconds, max|u|=( 41.46, 49.80, 12.34), ρ ∈ [0.0948, 1.1873], qᵛ ∈ [-8.55e-07, 0.0232], qʳ ∈ [-3e-11, 0.00039]
[ Info: iter= 600, t=25.204 minutes, Δt=2.321 seconds, max|u|=( 41.40, 50.01, 7.49), ρ ∈ [0.0948, 1.1873], qᵛ ∈ [-1.558e-07, 0.02426], qʳ ∈ [-7.3e-11, 0.00064]
[ Info: iter= 700, t=29.073 minutes, Δt=2.744 seconds, max|u|=( 41.33, 50.06, 16.92), ρ ∈ [0.0948, 1.1872], qᵛ ∈ [8.611e-08, 0.02569], qʳ ∈ [-5.9e-11, 0.00098]
[ Info: iter= 800, t=35.954 minutes, Δt=3.318 seconds, max|u|=( 48.21, 50.43, 5.71), ρ ∈ [0.0948, 1.1869], qᵛ ∈ [-3.408e-07, 0.02797], qʳ ∈ [-1.8e-10, 0.002]
[ Info: iter= 900, t=39.610 minutes, Δt=1.966 seconds, max|u|=( 57.89, 50.60, 16.04), ρ ∈ [0.0948, 1.1872], qᵛ ∈ [-4.292e-07, 0.02912], qʳ ∈ [-1.5e-10, 0.0027]
[ Info: iter=1000, t=42.733 minutes, Δt=2.030 seconds, max|u|=( 47.65, 50.89, 21.75), ρ ∈ [0.0948, 1.1874], qᵛ ∈ [-4.956e-07, 0.03011], qʳ ∈ [-4.1e-11, 0.0035]
[ Info: iter=1100, t=47.016 minutes, Δt=2.474 seconds, max|u|=( 43.68, 51.18, 18.62), ρ ∈ [0.0948, 1.1874], qᵛ ∈ [-7.632e-07, 0.03145], qʳ ∈ [-4e-10, 0.005]
[ Info: iter=1200, t=51.040 minutes, Δt=2.035 seconds, max|u|=( 54.64, 51.54, 37.30), ρ ∈ [0.0948, 1.1875], qᵛ ∈ [-9.585e-07, 0.03267], qʳ ∈ [-4.7e-10, 0.006]
[ Info: iter=1300, t=54.873 minutes, Δt=3.105 seconds, max|u|=( 107.20, 77.49, 43.21), ρ ∈ [0.0948, 1.1873], qᵛ ∈ [-1.143e-06, 0.03376], qʳ ∈ [-8.9e-11, 0.0073]
[ Info: iter=1400, t=1.023 hours, Δt=4.471 seconds, max|u|=( 90.64, 81.38, 40.76), ρ ∈ [0.0947, 1.1874], qᵛ ∈ [-1.59e-06, 0.03534], qʳ ∈ [-3.3e-09, 0.0062]
[ Info: iter=1500, t=1.092 hours, Δt=2.001 seconds, max|u|=( 89.55, 80.30, 46.17), ρ ∈ [0.0947, 1.1872], qᵛ ∈ [-1.638e-06, 0.03613], qʳ ∈ [-2.8e-09, 0.0061]
[ Info: iter=1600, t=1.159 hours, Δt=2.234 seconds, max|u|=( 95.51, 92.58, 43.55), ρ ∈ [0.0947, 1.1874], qᵛ ∈ [-1.65e-06, 0.0367], qʳ ∈ [-2e-10, 0.0061]
[ Info: iter=1700, t=1.244 hours, Δt=3.805 seconds, max|u|=( 99.98, 89.66, 42.65), ρ ∈ [0.0945, 1.1874], qᵛ ∈ [-1.722e-06, 0.03706], qʳ ∈ [-2.1e-10, 0.0061]
[ Info: iter=1800, t=1.327 hours, Δt=1.884 seconds, max|u|=( 101.98, 96.67, 43.62), ρ ∈ [0.0938, 1.1875], qᵛ ∈ [-1.777e-06, 0.03722], qʳ ∈ [-1e-10, 0.0062]
[ Info: iter=1900, t=1.379 hours, Δt=1.957 seconds, max|u|=( 102.49, 89.50, 43.90), ρ ∈ [0.0937, 1.1877], qᵛ ∈ [-1.675e-06, 0.03733], qʳ ∈ [-6.5e-11, 0.0062]
[ Info: iter=2000, t=1.447 hours, Δt=2.980 seconds, max|u|=( 106.88, 95.12, 44.02), ρ ∈ [0.0938, 1.1876], qᵛ ∈ [-1.712e-06, 0.03783], qʳ ∈ [-5.4e-10, 0.0061]
[ Info: iter=2100, t=1.538 hours, Δt=4.153 seconds, max|u|=( 103.76, 93.21, 41.85), ρ ∈ [0.0941, 1.1875], qᵛ ∈ [-1.879e-06, 0.03891], qʳ ∈ [-1.7e-08, 0.0055]
[ Info: iter=2200, t=1.638 hours, Δt=3.546 seconds, max|u|=( 108.46, 84.54, 41.63), ρ ∈ [0.0947, 1.1872], qᵛ ∈ [-2.03e-06, 0.04049], qʳ ∈ [-3.3e-10, 0.0056]
[ Info: iter=2300, t=1.698 hours, Δt=2.203 seconds, max|u|=( 93.36, 80.26, 40.30), ρ ∈ [0.0947, 1.1872], qᵛ ∈ [-2.05e-06, 0.04127], qʳ ∈ [-1.9e-10, 0.0055]
[ Info: iter=2400, t=1.761 hours, Δt=2.102 seconds, max|u|=( 90.88, 83.54, 39.47), ρ ∈ [0.0947, 1.1874], qᵛ ∈ [-2.162e-06, 0.04191], qʳ ∈ [-9.2e-11, 0.0055]
[ Info: iter=2500, t=1.821 hours, Δt=2.591 seconds, max|u|=( 96.88, 80.80, 37.56), ρ ∈ [0.0947, 1.1876], qᵛ ∈ [-2.296e-06, 0.04233], qʳ ∈ [-7e-10, 0.0051]
[ Info: iter=2600, t=1.897 hours, Δt=3.939 seconds, max|u|=( 77.37, 84.60, 33.44), ρ ∈ [0.0946, 1.1878], qᵛ ∈ [-2.492e-06, 0.04261], qʳ ∈ [-1e-09, 0.0053]
[ Info: iter=2700, t=2.017 hours, Δt=4.135 seconds, max|u|=( 97.18, 91.44, 30.45), ρ ∈ [0.0945, 1.1872], qᵛ ∈ [-2.789e-06, 0.04249], qʳ ∈ [-9.9e-10, 0.0059]
[ Info: iter=2800, t=2.097 hours, Δt=2.117 seconds, max|u|=( 90.97, 96.45, 32.55), ρ ∈ [0.0946, 1.1873], qᵛ ∈ [-2.767e-06, 0.04206], qʳ ∈ [-5.5e-11, 0.0066]
[ Info: iter=2900, t=2.205 hours, Δt=4.430 seconds, max|u|=( 93.94, 96.31, 40.87), ρ ∈ [0.0947, 1.1877], qᵛ ∈ [-2.884e-06, 0.04128], qʳ ∈ [-6.5e-09, 0.0062]
[ Info: iter=3000, t=2.281 hours, Δt=2.350 seconds, max|u|=( 91.82, 88.66, 41.37), ρ ∈ [0.0947, 1.1877], qᵛ ∈ [-2.959e-06, 0.04066], qʳ ∈ [-1.9e-09, 0.0059]
[ Info: iter=3100, t=2.341 hours, Δt=2.024 seconds, max|u|=( 93.24, 88.18, 40.15), ρ ∈ [0.0947, 1.1877], qᵛ ∈ [-3.001e-06, 0.04014], qʳ ∈ [-1.2e-09, 0.0056]
[ Info: iter=3200, t=2.404 hours, Δt=2.378 seconds, max|u|=( 90.76, 88.59, 36.53), ρ ∈ [0.0947, 1.1878], qᵛ ∈ [-3.05e-06, 0.03961], qʳ ∈ [-3.6e-10, 0.0056]
[ Info: iter=3300, t=2.469 hours, Δt=2.477 seconds, max|u|=( 87.14, 84.98, 35.43), ρ ∈ [0.0947, 1.1878], qᵛ ∈ [-3.239e-06, 0.03908], qʳ ∈ [-6.7e-10, 0.005]
[ Info: iter=3400, t=2.569 hours, Δt=3.212 seconds, max|u|=( 89.90, 89.98, 34.77), ρ ∈ [0.0947, 1.1884], qᵛ ∈ [-4.047e-06, 0.03969], qʳ ∈ [-1.4e-08, 0.0059]
[ Info: iter=3500, t=2.634 hours, Δt=1.756 seconds, max|u|=( 77.43, 81.37, 37.32), ρ ∈ [0.0948, 1.1883], qᵛ ∈ [-4.22e-06, 0.04132], qʳ ∈ [-1.3e-10, 0.0056]
[ Info: iter=3600, t=2.688 hours, Δt=2.077 seconds, max|u|=( 82.04, 83.24, 38.64), ρ ∈ [0.0948, 1.1880], qᵛ ∈ [-4.314e-06, 0.04201], qʳ ∈ [-9.8e-10, 0.0061]
[ Info: iter=3700, t=2.747 hours, Δt=3.422 seconds, max|u|=( 92.45, 88.95, 34.20), ρ ∈ [0.0948, 1.1882], qᵛ ∈ [-4.53e-06, 0.04235], qʳ ∈ [-2.2e-09, 0.0065]
[ Info: iter=3800, t=2.808 hours, Δt=1.814 seconds, max|u|=( 97.63, 92.28, 36.01), ρ ∈ [0.0948, 1.1888], qᵛ ∈ [-4.828e-06, 0.04233], qʳ ∈ [-9.1e-10, 0.0068]
[ Info: iter=3900, t=2.868 hours, Δt=1.926 seconds, max|u|=( 100.36, 88.32, 38.46), ρ ∈ [0.0948, 1.1888], qᵛ ∈ [-5.159e-06, 0.04191], qʳ ∈ [-8e-10, 0.0065]
[ Info: iter=4000, t=2.946 hours, Δt=3.976 seconds, max|u|=( 105.25, 98.17, 41.29), ρ ∈ [0.0947, 1.1888], qᵛ ∈ [-5.796e-06, 0.04078], qʳ ∈ [-1.4e-08, 0.0057]
[ Info: Simulation is stopping after running for 10.408 minutes.
[ Info: Simulation time 3 hours equals or exceeds stop time 3 hours.
Cascade animation
A 2-row × 5-column animation of the downscaling: row 1 the ERA5 parent (dashed rectangle = child extent), row 2 the Breeze child. Columns are near-surface θᵛ′ (referenced pointwise to t=0, so the terrain/stratification background cancels and the cold pool stands out) and wind speed |U|, then w, qᵛ, qʳ at 2 km. The saved slices load straight back as FieldTimeSeries and plot directly — heatmap! reads each 2-D Field's own coordinates and moves it host-side, so there is no manual indexing, and the θᵛ′ anomaly θᵛ[n] - θᵛ[1] is a lazy field operation.
θᵛ_series = FieldTimeSeries(surface_filename, "θᵛ")
U_series = FieldTimeSeries(surface_filename, "U")
w_series = FieldTimeSeries(aloft_filename, "w")
qᵛ_series = FieldTimeSeries(aloft_filename, "qᵛ")
qʳ_series = FieldTimeSeries(aloft_filename, "qʳ")
times = θᵛ_series.times0.0f0:1200.0f0:10800.0f0The parent is prescribed (no output file). Its rows are reconstructed at each frame time from the parent's own raw ERA5 fields — which are full-in-memory FieldTimeSeries — via breeze_prognostic_state, the exported (T, qᵛ, p) → (ρ, θˡⁱ, qᵗ) transform (the same physics the lateral BCs apply). We do NOT read the StateExchanger prognostics here: that FTS is a 2-level window, so post-run per-frame reads would alias to whatever the last exchange_state! left resident. ERA5's w ≈ −ω/(ρg); no rain.
k_parent = searchsortedfirst(Array(znodes(parent.grid, Center())), 2000)
constants = child.thermodynamic_constants
pˢᵗ = child.dynamics.standard_pressure
g = constants.gravitational_acceleration
ε = vapor_gas_constant(constants) / dry_air_gas_constant(constants) - 1 # virtual-temperature coefficient
slice(operation, k) = Field(operation, indices = (:, :, k))
function parent_slices(t)
hydrometeors = parent.microphysical_variables
qˡ = hydrometeors.qᶜˡ[Time(t)] + hydrometeors.qʳ[Time(t)] ## total liquid: cloud + rain
qⁱ = hydrometeors.qᶜⁱ[Time(t)] + hydrometeors.qˢ[Time(t)] ## total ice: cloud ice + snow
(; ρ, θˡⁱ, qᵗ) = breeze_prognostic_state(constants, pˢᵗ,
parent.temperature[Time(t)], parent.specific_humidity[Time(t)],
qˡ, qⁱ, parent.pressure)
u, v = parent.velocities.u[Time(t)], parent.velocities.v[Time(t)]
# `u` and `v` are staggered, so |U| would inherit `u`'s location; force it to centers, where the
# other three frames already live, so every parent panel shares one prototype location.
speed = @at((Center, Center, Center), sqrt(u^2 + v^2))
return (θᵛ = slice(θˡⁱ * (1 + ε * qᵗ), 1),
U = slice(speed, 1),
w = slice(-ω_series[Time(t)] / (ρ * g), k_parent),
qᵛ = slice(qᵗ, k_parent))
end
parent_frames = [parent_slices(t) for t in times]One column each: title, colormap, the child/parent 2-D fields as functions of the frame n, and the parent level k those frames are sliced at. nothing parent ⇒ blank panel (rain). Fixed, illustrative color ranges (this is a demo, not a tuned figure); w alone differs ~10× between parent and child, so it carries a range per row.
g_per_kg(field) = 1f3 * field
columns = [(title = "θᵛ′ₛ (K)", colormap = :balance, child_range = (-6, 6), parent_range = (-6, 6),
k = 1, child_series = θᵛ_series,
child = n -> θᵛ_series[n] - θᵛ_series[1], parent = n -> parent_frames[n].θᵛ - parent_frames[1].θᵛ),
(title = "|U|ₛ (m s⁻¹)", colormap = :speed, child_range = (0, 30), parent_range = (0, 30),
k = 1, child_series = U_series,
child = n -> U_series[n], parent = n -> parent_frames[n].U),
(title = "w₂ₖₘ (m s⁻¹)", colormap = :balance, child_range = (-3, 3), parent_range = (-0.3, 0.3),
k = k_parent, child_series = w_series,
child = n -> w_series[n], parent = n -> parent_frames[n].w),
(title = "qᵛ₂ₖₘ (g kg⁻¹)", colormap = :dense, child_range = (0, 15), parent_range = (0, 15),
k = k_parent, child_series = qᵛ_series,
child = n -> g_per_kg(qᵛ_series[n]), parent = n -> g_per_kg(parent_frames[n].qᵛ)),
(title = "qʳ₂ₖₘ (g kg⁻¹)", colormap = :dense, child_range = (0, 2), parent_range = (0, 2),
k = k_parent, child_series = qʳ_series,
child = n -> g_per_kg(qʳ_series[n]), parent = n -> nothing)]
boxλ = [λ_west, λ_east, λ_east, λ_west, λ_west]
boxφ = [φ_south, φ_south, φ_north, φ_north, φ_south]
fig_cascade = Figure(size = (1500, 700))
cascade_n = Observable(1)
# Each panel plots one plain field (`similar(prototype)`: the prototype's location, grid and indices,
# its own data), refilled per frame and handed back as the *same* object. `heatmap!` reads the field's
# own λ/φ and moves it host-side, so no node arrays and no manual materialization. Reusing one object
# matters twice: Makie's compute graph otherwise compares successive plot inputs elementwise, which on
# the parent's GPU fields is scalar indexing, and `heatmap!` calls `compute!` on what it is handed,
# which on a field carrying an operand would overwrite the frame just copied in.
function panel!(ax, prototype, field_of, colormap, colorrange)
panel = similar(prototype)
frame = Observable(panel)
on(nn -> (set!(panel, compute!(field_of(nn))); frame[] = panel), cascade_n)
set!(panel, compute!(field_of(1)))
return heatmap!(ax, frame; colormap, colorrange)
end
Label(fig_cascade[0, 1:5],
@lift(@sprintf("MC3E 20 May 2011 — ERA5 → 12 km Breeze — t = %.1f h", times[$cascade_n] / 3600)),
fontsize = 20, tellwidth = false)
for (i, column) in enumerate(columns)
parent_ax = Axis(fig_cascade[1, i]; title = column.title, aspect = DataAspect())
child_ax = Axis(fig_cascade[2, i]; aspect = DataAspect())
hmc = panel!(child_ax, column.child_series[1], column.child, column.colormap, column.child_range)
if isnothing(column.parent(1))
text!(parent_ax, 0.5, 0.5; text = "no rain", space = :relative, align = (:center, :center), color = :gray)
Colorbar(fig_cascade[3, i], hmc; vertical = false, flipaxis = false, height = 10)
else
parent_prototype = CenterField(parent.grid, indices = (:, :, column.k))
hmp = panel!(parent_ax, parent_prototype, column.parent, column.colormap, column.parent_range)
lines!(parent_ax, boxλ, boxφ; color = :black, linestyle = :dash, linewidth = 1.5)
if column.parent_range == column.child_range
Colorbar(fig_cascade[3, i], hmc; vertical = false, flipaxis = false, height = 10)
else ## w differs ~10× parent-to-child — a colorbar per row
Colorbar(fig_cascade[3, i], hmp; vertical = false, flipaxis = false, height = 10, label = "ERA5")
Colorbar(fig_cascade[4, i], hmc; vertical = false, flipaxis = false, height = 10, label = "12 km")
end
end
hidedecorations!(parent_ax); hidedecorations!(child_ax)
if i == 1
text!(parent_ax, 0.03, 0.97; text = "ERA5", space = :relative, align = (:left, :top), color = :white, fontsize = 15, font = :bold)
text!(child_ax, 0.03, 0.97; text = "Breeze 12 km", space = :relative, align = (:left, :top), color = :white, fontsize = 15, font = :bold)
end
end
CairoMakie.record(fig_cascade, "era5_cascade_2row.mp4", eachindex(times); framerate = 8) do nn
cascade_n[] = nn ## `record` is also exported by CUDA, so qualify it
endSection animation
The vertical structure along the zonal x-z section through the ARM SGP site: wind speed (the nocturnal Great Plains low-level jet near the surface and the upper-level jet aloft), vertical velocity (the convective updrafts and the squall line's mesoscale circulation), and water vapor (the moist tongue feeding the convection). The section slices load back as FieldTimeSeries of 2-D fields on the terrain-following grid, and heatmap! maps them to physical $(λ, z)$ coordinates directly — one @lift per panel is the whole animation machinery.
U_section = FieldTimeSeries(section_filename, "U")
w_section = FieldTimeSeries(section_filename, "w")
qᵛ_section = FieldTimeSeries(section_filename, "qᵛ")
φ_section = φ_centers[j_section]
# Type-erased observables: the per-frame fields' concrete types can differ in hidden
# parameters, which a `@lift`-typed Observable rejects at the second frame.
section_n = Observable(1)
Uₙ = Observable{Any}(U_section[1])
wₙ = Observable{Any}(w_section[1])
qᵛₙ = Observable{Any}(g_per_kg(qᵛ_section[1]))
on(section_n) do n
Uₙ[] = U_section[n]
wₙ[] = w_section[n]
qᵛₙ[] = g_per_kg(qᵛ_section[n])
end
fig_section = Figure(size = (1000, 900))
Label(fig_section[0, 1:2],
@lift(@sprintf("MC3E 20 May 2011 — vertical structure at %.2f°N — t = %.1f h",
φ_section, times[$section_n] / 3600)),
fontsize = 20, tellwidth = false)
ax_U = Axis(fig_section[1, 1]; ylabel = "z (m)", title = "|U| (m s⁻¹)")
ax_w = Axis(fig_section[2, 1]; ylabel = "z (m)", title = "w (m s⁻¹)")
ax_q = Axis(fig_section[3, 1]; ylabel = "z (m)", xlabel = "longitude (°)", title = "qᵛ (g kg⁻¹)")
hm_U = heatmap!(ax_U, Uₙ; colormap = :speed, colorrange = (0, 40))
hm_w = heatmap!(ax_w, wₙ; colormap = :balance, colorrange = (-3, 3))
hm_q = heatmap!(ax_q, qᵛₙ; colormap = :dense, colorrange = (0, 15))
Colorbar(fig_section[1, 2], hm_U)
Colorbar(fig_section[2, 2], hm_w)
Colorbar(fig_section[3, 2], hm_q)
CairoMakie.record(fig_section, "era5_section_xz.mp4", eachindex(times); framerate = 8) do nn
section_n[] = nn
endLand response
The coupled surface closes the loop: skin-temperature change over the run (daytime heating, plus cooling under the squall line's cold pool and anvil shading) and the final bucket saturation — wetted where the squall line rained, dried by evaporation elsewhere. Both load back as 2-D FieldTimeSeries on the land grid and plot directly.
Tˡᵃ_series = FieldTimeSeries(land_filename, "Tˡᵃ")
𝒮_series = FieldTimeSeries(land_filename, "𝒮")
fig_land = Figure(size = (1100, 420))
ax_T = Axis(fig_land[1, 1]; xlabel = "longitude (°)", ylabel = "latitude (°)",
title = "ΔTˡᵃ over $(prettytime(duration)) (K)")
ax_𝒮 = Axis(fig_land[1, 3]; xlabel = "longitude (°)", title = "final surface saturation 𝒮")
hm_T = heatmap!(ax_T, Tˡᵃ_series[end] - Tˡᵃ_series[1]; colormap = :balance, colorrange = (-8, 8))
hm_𝒮 = heatmap!(ax_𝒮, 𝒮_series[end]; colormap = :dense, colorrange = (0, 1))
Colorbar(fig_land[1, 2], hm_T)
Colorbar(fig_land[1, 4], hm_𝒮)
save("era5_breeze_land_response.png", fig_land)
fig_land
This page was generated using Literate.jl.