Diurnal cycle of radiative convection

This example simulates a dramatic diurnal cycle of moist convection over a tropical land surface. During the day, the sun heats the ground, driving vigorous boundary-layer convection that lofts moisture into towering cumulus clouds. At night, the surface cools rapidly by longwave emission, stabilizing the boundary layer and shutting off convection.

The key ingredient is a time-varying surface temperature that follows the sun: peaking in the early afternoon and dipping well below the air temperature at night. This creates a strong diurnal contrast — afternoon thunderstorms that die at sunset and don't return until the next morning.

Interactive all-sky radiation with the ecCKD gas optics of NumericalRadiation.jl computes spectrally-resolved shortwave and longwave fluxes, with cloud scattering tables mapped onto the same spectral g-points. Saturation-adjustment microphysics diagnoses cloud liquid water that feeds back on the radiation. A stretched vertical grid resolves the cloud layer (100 m spacing below 3 km) while extending to 15 km, and the radiation model continues the column above the grid top through a standard-atmosphere stratosphere to 65 km. A sponge layer above 8 km prevents spurious temperature drift in the coarse upper cells.

using Breezeusing Oceananigansusing Oceananigans.Unitsusing Dates: DateTimeusing Printf, Random, Statisticsusing CairoMakieusing NCDatasets  # Required for the ecCKD lookup tablesusing NumericalRadiation: NumericalRadiationusing CUDARandom.seed!(2025)if CUDA.functional()    CUDA.seed!(2025)end

Grid

We use a 2D vertical slice (x-z) that is periodic in x and bounded in z. The vertical grid is stretched: fine 100 m cells resolve the cloud layer below 3 km, then a smooth transition to 1 km cells carries the column up to 15 km. The radiation model's OverlyingAtmosphereColumn adds the stratosphere above the grid top, so the total cell count stays modest.

Nx = 128Lx = 12800   # 12.8 kmarch = GPU()Oceananigans.defaults.FloatType = Float32z = PiecewiseStretchedDiscretization(    z  = [0, 3000, 8000, 15000],    Δz = [100,  100, 1000,  1000])Nz = length(z) - 1grid = RectilinearGrid(arch;                       size = (Nx, Nz),                       x = (0, Lx),                       z,                       halo = (5, 5),                       topology = (Periodic, Flat, Bounded))

128×1×51 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CUDAGPU with 5×0×5 halo
├── Periodic x ∈ [0.0, 12800.0) regularly spaced with Δx=100.0
├── Flat y                      
└── Bounded  z ∈ [0.0, 15000.0] variably spaced with min(Δz)=100.0, max(Δz)=1000.0

Reference state

p₀ = 101325  # Surface pressure [Pa]θ₀ = 300     # Reference potential temperature [K]constants = ThermodynamicConstants()reference_state = ReferenceState(grid, constants;                                 base_pressure = p₀,                                 potential_temperature = θ₀,                                 vapor_mass_fraction = 0)dynamics = AnelasticDynamics(reference_state)

AnelasticDynamics(p₀=101325.0, θ₀=300.0)
├── pressure_solver_float_type: Float64
└── pressure_anomaly: not materialized

Background atmosphere

The ecCKD gas optics require trace gas concentrations to compute spectral absorption and emission. We specify well-mixed greenhouse gas concentrations and a tropical ozone profile that transitions from low tropospheric values to a stratospheric peak near 25 km. The peak lies above the 15 km grid top, but the radiation model evaluates the ozone profile in its overlying atmosphere too, so the stratospheric ozone still absorbs sunlight before it reaches the grid.

@inline function tropical_ozone(z)    FT = float(typeof(z))   # in the grid's precision, which the radiation kernels evaluate it in    troposphere_O₃ = FT(30e-9) * (1 + FT(0.5) * z / 10_000)    zˢᵗ = FT(25e3)    Hˢᵗ = FT(5e3)    stratosphere_O₃ = FT(8e-6) * exp(-((z - zˢᵗ) / Hˢᵗ)^2)    χˢᵗ = 1 / (1 + exp(-(z - FT(15e3)) / 2))    return troposphere_O₃ * (1 - χˢᵗ) + stratosphere_O₃ * χˢᵗendbackground_atmosphere = BackgroundAtmosphere(    CO₂ = 348e-6,    CH₄ = 1650e-9,    N₂O = 306e-9,    O₃ = tropical_ozone)

BackgroundAtmosphere with 6 active gases:
  N₂ = 0.78084
  O₂ = 0.20946
  CO₂ = 348.0 ppm
  CH₄ = 1.65 ppm
  N₂O = 306.0 ppb
  O₃ = tropical_ozone (generic function with 1 method)

Diurnal surface temperature

Over land, the surface temperature swings dramatically with the sun. We model this as a sinusoidal cycle peaking 2 hours after solar noon: Tₛ(t) = T̄ₛ + ΔTₛ cos(2π(t - t_peak) / 24h), with T̄ₛ = 300 K (mean) and ΔTₛ = 10 K (amplitude). This gives 310 K in early afternoon and 290 K at night — a 20 K diurnal range, typical of tropical semi-arid land.

We start at midnight (t = 0), so the surface starts cold (290 K), warms through the morning, peaks at t = 14 h (2 pm local), and cools at night. A Field stores the surface temperature and a callback updates it each time step, keeping both the bulk fluxes and the radiation in sync.

T̄ₛ = 300   # Mean surface temperature [K]ΔTₛ = 20   # Diurnal amplitude [K]Tₛ = Field{Center, Center, Nothing}(grid)set!(Tₛ, T̄ₛ - ΔTₛ)  # Start at midnight minimum

128×1×1 Field{Center, Center, Nothing} reduced over dims = (3,) on RectilinearGrid on CUDAGPU
├── grid: 128×1×51 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CUDAGPU with 5×0×5 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: Nothing, top: Nothing, immersed: Nothing
└── data: 138×1×1 StaticStridesOffsetArray(::CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}, -4:133, 1:1, 1:1) with eltype Float32 with indices -4:133×1:1×1:1
    └── max=280.0, min=280.0, mean=280.0

Radiation with a diurnal cycle

We place the domain at 15°N latitude on the prime meridian and start at midnight on the spring equinox (March 20). The sun rises at t ≈ 6 h, reaches noon at t ≈ 12 h, and sets at t ≈ 18 h.

latitude = 15 # 15°Noptics = EcCKDOptics(clouds = CloudScatteringTables())radiation = RadiativeTransferModel(grid, optics, constants;                                   surface_temperature = Tₛ,                                   surface_albedo = 0.20,                                   surface_emissivity = 0.95,                                   solar_constant = 1361,                                   background_atmosphere,                                   solar_position = ApparentSolarPosition(coordinate = (0, latitude),                                                                          epoch      = DateTime(2020, 3, 20, 0, 0, 0)),                                   schedule = TimeInterval(5minutes),                                   liquid_effective_radius = ConstantRadiusParticles(10e-6),                                   ice_effective_radius = ConstantRadiusParticles(30e-6))

RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── solar_position: ApparentSolarPosition(coordinate=(0, 15), epoch=2020-03-20T00:00:00)
├── surface_temperature: 128×1×1 Field{Center, Center, Nothing} reduced over dims = (3,) on RectilinearGrid on CUDAGPU
├── surface_emissivity: ConstantField(0.95)
├── direct_surface_albedo: ConstantField(0.2)
├── liquid_effective_radius: Breeze.AtmosphereModels.ConstantRadiusParticles{Float32}(1.0f-5)
├── ice_effective_radius: Breeze.AtmosphereModels.ConstantRadiusParticles{Float32}(3.0f-5)
├── diffuse_surface_albedo: ConstantField(0.2)
├── optics: EcCKDOptics with 32 longwave and 32 shortwave g-points, all sky
└── overlying_atmosphere: 40 layers from 15000.0 m to 65000.0 m

Surface fluxes

Bulk aerodynamic formulae provide surface sensible heat, moisture, and momentum fluxes driven by the time-varying surface temperature. During the day Tₛ > Tair drives strong upward fluxes; at night Tₛ < Tair can produce downward fluxes that cool the boundary layer.

Cᴰ = Cᵀ = 1e-3Cᵛ = 1.2e-3Uᵍ = 1  # Gustiness [m/s]ρθ_flux = BulkSensibleHeatFlux(coefficient=Cᵀ, gustiness=Uᵍ, surface_temperature=Tₛ)ρqᵗ_flux = BulkVaporFlux(coefficient=Cᵛ, gustiness=Uᵍ, surface_temperature=Tₛ)ρθ_bcs = FieldBoundaryConditions(bottom=ρθ_flux)ρqᵗ_bcs = FieldBoundaryConditions(bottom=ρqᵗ_flux)ρu_bcs = FieldBoundaryConditions(bottom=Breeze.BulkDrag(coefficient=Cᴰ, gustiness=Uᵍ))

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
├── bottom: FluxBoundaryCondition: BulkDragFunction(direction=Nothing, coefficient=0.001, gustiness=1)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: nothing)

Microphysics

Warm-phase saturation adjustment instantly converts supersaturated moisture to cloud liquid water.

microphysics = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())

Breeze.Microphysics.SaturationAdjustment{Breeze.Thermodynamics.WarmPhaseEquilibrium, Breeze.Solvers.SecantSolver{Float32}}(Breeze.Thermodynamics.WarmPhaseEquilibrium(), SecantSolver(reltol=0.0, abstol=0.0001, maxiter=20))

Stratospheric sponge

The domain extends to 15 km, but the initial upper troposphere and lower stratosphere aren't in radiative equilibrium: ozone absorbs shortwave radiation and the coarse upper cells respond strongly. A Newtonian relaxation of temperature toward the initial profile above 8 km keeps the upper cells anchored without affecting the tropospheric dynamics. We apply this as an energy forcing on ρE, which Breeze automatically converts to a ρθ tendency.

Tᵣ = reference_state.temperatureρᵣ = reference_state.densitycᵖᵈ = constants.dry_air.heat_capacity  # J/(kg·K)τ_sponge = 6hours@inline function stratospheric_relaxation(i, j, k, grid, clock, model_fields, p)    @inbounds T = model_fields.T[i, j, k]    @inbounds Tᵣ = p.Tᵣ[i, j, k]    @inbounds ρ = p.ρᵣ[i, j, k]    z = znode(i, j, k, grid, Center(), Center(), Center())    α = clamp((z - 8000) / 4000, 0, 1)    ∂T∂t = -α * (T - Tᵣ) / p.τ    return ρ * p.cᵖᵈ * ∂T∂tend

stratospheric_relaxation (generic function with 1 method)

The parameters take the grid's float type, as the kernels compute in it

FT = eltype(grid)sponge = Forcing(stratospheric_relaxation; discrete_form=true,                 parameters=(; Tᵣ, ρᵣ, cᵖᵈ=convert(FT, cᵖᵈ), τ=convert(FT, τ_sponge)))forcing = (; ρE=sponge)

(ρE = DiscreteForcing{@NamedTuple{Tᵣ::Oceananigans.Fields.Field{Nothing, Nothing, Oceananigans.Grids.Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (62,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (61,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (62,), (-6,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (61,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (138,), (-5,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (138,), (-5,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 1, 51, 5, 0, 5}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (1, 1, 61), (0, 0, -5), CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Nothing, Nothing, Nothing, Nothing, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition{Oceananigans.BoundaryConditions.Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition{Oceananigans.BoundaryConditions.Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Nothing, @NamedTuple{south_and_north::Nothing, west_and_east::Nothing, bottom_and_top::KernelAbstractions.Kernel{CUDACore.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.gpu__fill_bottom_and_top_halo!)}}, @NamedTuple{south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Nothing, Nothing}, bottom_and_top::Tuple{Oceananigans.BoundaryConditions.NoFluxBoundaryCondition{Oceananigans.BoundaryConditions.Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition{Oceananigans.BoundaryConditions.Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}}}}, Nothing, Nothing}, ρᵣ::Oceananigans.Fields.Field{Nothing, Nothing, Oceananigans.Grids.Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (62,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (61,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (62,), (-6,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (61,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (138,), (-5,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (138,), (-5,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 1, 51, 5, 0, 5}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (1, 1, 61), (0, 0, -5), CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Nothing, Nothing, Nothing, Nothing, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Value{Nothing}, Oceananigans.Fields.Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (62,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (61,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (62,), (-6,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (61,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (138,), (-5,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (138,), (-5,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 1, 51, 5, 0, 5}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (138, 1, 1), (-5, 0, 0), CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, Nothing, @NamedTuple{bottom_and_top::Nothing, south_and_north::Nothing, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{CUDACore.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.gpu__fill_periodic_west_and_east_halo!)}, 128, 5}}, @NamedTuple{bottom_and_top::Tuple{Nothing, Nothing}, south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}}, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition{Oceananigans.BoundaryConditions.Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Nothing, @NamedTuple{south_and_north::Nothing, west_and_east::Nothing, bottom_and_top::KernelAbstractions.Kernel{CUDACore.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.gpu__fill_bottom_and_top_halo!)}}, @NamedTuple{south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Nothing, Nothing}, bottom_and_top::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Value{Nothing}, Oceananigans.Fields.Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (62,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (61,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (62,), (-6,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (61,), (-5,), CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (138,), (-5,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (138,), (-5,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 1, 51, 5, 0, 5}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (138, 1, 1), (-5, 0, 0), CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, Nothing, @NamedTuple{bottom_and_top::Nothing, south_and_north::Nothing, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{CUDACore.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.gpu__fill_periodic_west_and_east_halo!)}, 128, 5}}, @NamedTuple{bottom_and_top::Tuple{Nothing, Nothing}, south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}}, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition{Oceananigans.BoundaryConditions.Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}}}}, Nothing, Nothing}, cᵖᵈ::Float32, τ::Float32}}
├── func: stratospheric_relaxation (generic function with 1 method)
└── parameters: (Tᵣ = 1×1×51 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on RectilinearGrid on CUDAGPU
├── grid: 128×1×51 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CUDAGPU with 5×0×5 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 1×1×61 StaticStridesOffsetArray(::CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}, 1:1, 1:1, -4:56) with eltype Float32 with indices 1:1×1:1×-4:56
    └── max=300.642, min=159.593, mean=264.678, ρᵣ = 1×1×51 Field{Nothing, Nothing, Center} reduced over dims = (1, 2) on RectilinearGrid on CUDAGPU
├── grid: 128×1×51 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CUDAGPU with 5×0×5 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Value, top: ZeroFlux, immersed: Nothing
└── data: 1×1×61 StaticStridesOffsetArray(::CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}, 1:1, 1:1, -4:56) with eltype Float32 with indices 1:1×1:1×-4:56
    └── max=1.16766, min=0.239471, mean=0.876104, cᵖᵈ = 1005.0f0, τ = 21600.0f0),)

Model assembly

coriolis = FPlane(; latitude)boundary_conditions = (ρθ=ρθ_bcs, ρqᵗ=ρqᵗ_bcs, ρu=ρu_bcs)weno_order = 5momentum_advection = WENO(order=weno_order)scalar_advection = (ρθ  = WENO(order=weno_order),                    ρqᵉ = WENO(order=weno_order, bounds=(0, 1)))model = AtmosphereModel(grid; dynamics, microphysics, radiation, forcing,                        momentum_advection, scalar_advection,                        boundary_conditions, coriolis)

AtmosphereModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×51 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CUDAGPU with 5×0×5 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float32}
├── timestepper: SSPRungeKutta3
├── 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, bounds=(0.0, 1.0))
├── forcing: ρE=>DiscreteForcing
├── tracers: ()
├── coriolis: FPlane{Float32}(f=3.77468e-5)
└── microphysics: SaturationAdjustment

Initial conditions

The sounding has a dry-adiabatic sub-cloud layer (0–1 km) capped by a conditionally unstable troposphere (5 K/km lapse rate) that transitions to an isothermal stratosphere at 210 K. Moisture is 20 g/kg at the surface with a 2.5 km scale height, typical of the tropical maritime boundary layer.

function Tᵇᵍ(z)    T₀ = 300    Tˢᵗ = 210 # stratosphere temperature    T = T₀ - 1e-3 * max(z, 1000) - 5e-3 * max(0, z - 1000)    return max(T, Tˢᵗ)enduᵢ(x, z) = -5 * max(1 - z / 3000, 0)

uᵢ (generic function with 1 method)

Random perturbations in the lowest 1 km trigger convection.

δT = 2δℋ = 1e-2zδ = 1000ϵ() = rand() - 0.5Tᵢ(x, z) = Tᵇᵍ(z) + δT * ϵ() * (z < zδ)ℋᵢ(x, z) = (0.5 + δℋ * ϵ()) * (z < zδ)

ℋᵢ (generic function with 1 method)

After setting initial conditions, we recompute the reference state from the horizontally-averaged model state. This is important for Float32 accuracy in tall domains: the default dry-adiabat reference state diverges from the actual stratospheric profile, causing large density errors that overwhelm Float32 precision. set_to_mean! adjusts ρᵣ to match the current state while rescaling density-weighted prognostic fields to preserve specific quantities.

set!(model; T=Tᵢ, ℋ=ℋᵢ, u=uᵢ)reference_state = model.dynamics.reference_stateset_to_mean!(reference_state, model, rescale_densities=true)T = model.temperatureqᵗ = specific_prognostic_moisture(model)u, w = model.velocities.u, model.velocities.wqˡ = model.microphysical_fields.qˡ@info "Diurnal Radiative Convection (2D)"@info "Grid: $(Nx) × $(Nz) (stretched), domain: $(Lx/1000) km × 15 km"@info "Initial T range: $(minimum(T)) – $(maximum(T)) K"@info "Initial qᵗ range: $(minimum(qᵗ)*1000) – $(maximum(qᵗ)*1000) g/kg"

[ Info: Diurnal Radiative Convection (2D)
[ Info: Grid: 128 × 51 (stretched), domain: 12.8 km × 15 km
[ Info: Initial T range: 238.17137 – 300.11438 K
[ Info: Initial qᵗ range: 0.0 – 12.133494 g/kg

Simulation

We run for two full diurnal cycles (48 hours) starting at midnight, so the on/off pattern of convection repeats convincingly.

simulation = Simulation(model; Δt=1, stop_time=3days)conjure_time_step_wizard!(simulation, cfl=0.7)Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)

Surface temperature callback

At each time step we update the surface temperature field following a cosine curve that peaks 2 hours after solar noon (t = 14 h local). The period is 24 hours and the simulation starts at midnight.

function update_surface_temperature!(sim)    t = time(sim)    t_peak = 14hour  # peak at 14:00 local (2 pm)    Tₛ_now = T̄ₛ + ΔTₛ * cos(2π * (t - t_peak) / day)    set!(Tₛ, Tₛ_now)    return nothingendadd_callback!(simulation, update_surface_temperature!, TimeInterval(1minute))wall_clock = Ref(time_ns())function progress(sim)    elapsed = 1e-9 * (time_ns() - wall_clock[])    wmax = maximum(abs, w)    Tmin, Tmax = extrema(T)    qˡmax = maximum(qˡ)    Tₛ_now = T̄ₛ + ΔTₛ * cos(2π * (time(sim) - 14hours) / 24hours)    OLR = mean(view(radiation.upwelling_longwave_flux, :, 1, Nz+1))    msg = @sprintf("Iter: %5d, t: %8s, Δt: %5.1fs, wall: %8s",                   iteration(sim), prettytime(sim), sim.Δt, prettytime(elapsed))    msg *= @sprintf(", max|w|: %5.2f m/s, T: [%5.1f, %5.1f] K, max(qˡ): %.2e",                   wmax, Tmin, Tmax, qˡmax)    msg *= @sprintf(", Tₛ: %.1f K, OLR: %.1f W/m²", Tₛ_now, OLR)    @info msg    wall_clock[] = time_ns()    return nothingendadd_callback!(simulation, progress, IterationInterval(1000))

Output

Horizontally-averaged profiles are saved every hour (time-averaged) and 2D slices every 10 minutes for animation.

qᵛ = model.microphysical_fields.qᵛFᴿ = radiation.flux_divergenceoutputs = (; u, w, T, qˡ, qᵛ, Fᴿ)avg_outputs = NamedTuple(name => Average(outputs[name], dims=1) for name in keys(outputs))filename = "radiative_convection"averages_filename = filename * "_averages.jld2"slices_filename = filename * "_slices.jld2"simulation.output_writers[:averages] = JLD2Writer(model, avg_outputs;                                                  filename = averages_filename,                                                  schedule = AveragedTimeInterval(1hour),                                                  overwrite_files = true)slice_outputs = (; w, qᵛ, T)simulation.output_writers[:slices] = JLD2Writer(model, slice_outputs;                                                filename = slices_filename,                                                schedule = TimeInterval(10minutes),                                                overwrite_files = true)@info "Starting simulation..."run!(simulation)@info "Simulation completed!"

[ Info: Starting simulation...
[ Info: Initializing simulation...
[ Info: Iter:     0, t: 0 seconds, Δt:   1.1s, wall: 52.542 seconds, max|w|:  0.00 m/s, T: [238.2, 300.1] K, max(qˡ): 0.00e+00, Tₛ: 282.7 K, OLR: 373.9 W/m²
[ Info:     ... simulation initialization complete (40.428 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (10.398 seconds).
[ Info: Iter:  1000, t: 2.762 hours, Δt:  14.1s, wall: 20.569 seconds, max|w|:  1.76 m/s, T: [230.5, 299.6] K, max(qˡ): 0.00e+00, Tₛ: 280.4 K, OLR: 359.9 W/m²
[ Info: Iter:  2000, t: 6.095 hours, Δt:  14.3s, wall: 4.368 seconds, max|w|:  2.27 m/s, T: [225.1, 299.3] K, max(qˡ): 0.00e+00, Tₛ: 290.4 K, OLR: 365.1 W/m²
[ Info: Iter:  3000, t: 10.490 hours, Δt:  26.4s, wall: 5.168 seconds, max|w|:  1.47 m/s, T: [221.3, 298.9] K, max(qˡ): 1.02e-03, Tₛ: 312.1 K, OLR: 394.3 W/m²
[ Info: Iter:  4000, t: 14.000 hours, Δt:   8.4s, wall: 4.423 seconds, max|w|:  5.23 m/s, T: [219.9, 304.2] K, max(qˡ): 0.00e+00, Tₛ: 320.0 K, OLR: 406.1 W/m²
[ Info: Iter:  5000, t: 15.856 hours, Δt:   5.9s, wall: 3.574 seconds, max|w|:  7.05 m/s, T: [219.1, 305.9] K, max(qˡ): 0.00e+00, Tₛ: 317.7 K, OLR: 396.9 W/m²
[ Info: Iter:  6000, t: 17.481 hours, Δt:   5.1s, wall: 3.653 seconds, max|w|:  6.59 m/s, T: [218.7, 305.7] K, max(qˡ): 0.00e+00, Tₛ: 312.3 K, OLR: 385.3 W/m²
[ Info: Iter:  7000, t: 19.262 hours, Δt:   7.1s, wall: 4.121 seconds, max|w|:  5.10 m/s, T: [218.4, 305.2] K, max(qˡ): 0.00e+00, Tₛ: 303.8 K, OLR: 368.3 W/m²
[ Info: Iter:  8000, t: 21.555 hours, Δt:   9.3s, wall: 6.079 seconds, max|w|:  3.03 m/s, T: [218.3, 304.8] K, max(qˡ): 0.00e+00, Tₛ: 292.1 K, OLR: 348.6 W/m²
[ Info: Iter:  9000, t: 1.003 days, Δt:  10.0s, wall: 3.636 seconds, max|w|:  2.73 m/s, T: [218.2, 304.7] K, max(qˡ): 0.00e+00, Tₛ: 282.5 K, OLR: 332.9 W/m²
[ Info: Iter: 10000, t: 1.108 days, Δt:   9.7s, wall: 3.928 seconds, max|w|:  3.90 m/s, T: [218.1, 303.7] K, max(qˡ): 0.00e+00, Tₛ: 280.2 K, OLR: 326.8 W/m²
[ Info: Iter: 11000, t: 1.219 days, Δt:  10.9s, wall: 3.862 seconds, max|w|:  4.01 m/s, T: [217.9, 302.8] K, max(qˡ): 0.00e+00, Tₛ: 286.8 K, OLR: 331.2 W/m²
[ Info: Iter: 12000, t: 1.344 days, Δt:  14.0s, wall: 4.223 seconds, max|w|:  2.66 m/s, T: [217.9, 302.6] K, max(qˡ): 6.82e-05, Tₛ: 301.3 K, OLR: 347.4 W/m²
[ Info: Iter: 13000, t: 1.548 days, Δt:   9.4s, wall: 5.167 seconds, max|w|:  5.76 m/s, T: [218.0, 304.2] K, max(qˡ): 3.07e-04, Tₛ: 319.5 K, OLR: 355.7 W/m²
[ Info: Iter: 14000, t: 1.626 days, Δt:   6.3s, wall: 3.752 seconds, max|w|:  8.88 m/s, T: [217.7, 305.7] K, max(qˡ): 9.28e-04, Tₛ: 319.3 K, OLR: 343.5 W/m²
[ Info: Iter: 15000, t: 1.695 days, Δt:   6.8s, wall: 3.929 seconds, max|w|:  6.04 m/s, T: [217.7, 305.5] K, max(qˡ): 8.58e-04, Tₛ: 315.3 K, OLR: 348.1 W/m²
[ Info: Iter: 16000, t: 1.773 days, Δt:   8.6s, wall: 3.963 seconds, max|w|:  4.48 m/s, T: [217.6, 305.2] K, max(qˡ): 5.20e-04, Tₛ: 307.4 K, OLR: 338.6 W/m²
[ Info: Iter: 17000, t: 1.895 days, Δt:  16.3s, wall: 4.728 seconds, max|w|:  1.46 m/s, T: [217.7, 304.1] K, max(qˡ): 0.00e+00, Tₛ: 292.4 K, OLR: 329.5 W/m²
[ Info: Iter: 18000, t: 2.079 days, Δt:  20.9s, wall: 5.884 seconds, max|w|:  0.85 m/s, T: [217.8, 303.2] K, max(qˡ): 1.02e-04, Tₛ: 280.0 K, OLR: 312.4 W/m²
[ Info: Iter: 19000, t: 2.306 days, Δt:  46.0s, wall: 8.105 seconds, max|w|:  0.61 m/s, T: [217.9, 301.4] K, max(qˡ): 0.00e+00, Tₛ: 296.6 K, OLR: 326.0 W/m²
[ Info: Iter: 20000, t: 2.597 days, Δt:   7.3s, wall: 9.010 seconds, max|w|:  4.73 m/s, T: [218.0, 304.2] K, max(qˡ): 0.00e+00, Tₛ: 319.9 K, OLR: 356.4 W/m²
[ Info: Iter: 21000, t: 2.670 days, Δt:   8.0s, wall: 7.158 seconds, max|w|:  6.82 m/s, T: [217.0, 305.3] K, max(qˡ): 1.09e-03, Tₛ: 317.1 K, OLR: 311.3 W/m²
[ Info: Iter: 22000, t: 2.757 days, Δt:   7.9s, wall: 3.653 seconds, max|w|:  5.44 m/s, T: [217.5, 305.2] K, max(qˡ): 2.37e-04, Tₛ: 309.2 K, OLR: 331.3 W/m²
[ Info: Iter: 23000, t: 2.861 days, Δt:  10.8s, wall: 3.735 seconds, max|w|:  3.18 m/s, T: [217.7, 304.6] K, max(qˡ): 0.00e+00, Tₛ: 296.6 K, OLR: 318.3 W/m²
[ Info: Simulation is stopping after running for 2.761 minutes.
[ Info: Simulation time 3 days equals or exceeds stop time 3 days.
[ Info: Simulation completed!

Hovmöller diagrams of mean profile evolution

Time–height Hovmöller diagrams show how the horizontally-averaged temperature, cloud liquid water, and radiative flux divergence evolve over two diurnal cycles. The vertical axis is zoomed to the lowest 6 km where the convective dynamics live.

Tts  = FieldTimeSeries(averages_filename, "T")qᵛts = FieldTimeSeries(averages_filename, "qᵛ")Fᴿts = FieldTimeSeries(averages_filename, "Fᴿ")times = Tts.timesNt = length(times)z = Oceananigans.Grids.znodes(Tts.grid, Center())t_hours = times ./ 3600

0.0:1.0:72.0

Build time–height matrices from the averaged profiles.

T_zt  = permutedims(interior(Tts,  1, 1, :, :))T_zt .-= mean(T_zt, dims=1)qᵛ_zt = permutedims(interior(qᵛts, 1, 1, :, :))Fᴿ_zt = permutedims(interior(Fᴿts, 1, 1, :, :))fig = Figure(size=(900, 800), fontsize=14)zmax = 12000axT  = Axis(fig[1, 1]; ylabel="z (m)", title="Temperature anomaly (K)",            limits=((t_hours[1], t_hours[end]), (0, zmax)))axqᵛ = Axis(fig[2, 1]; ylabel="z (m)", title="Specific humidity (kg/kg)",            limits=((t_hours[1], t_hours[end]), (0, zmax)))axF  = Axis(fig[3, 1]; ylabel="z (m)", xlabel="time (hours)", title="Radiative flux divergence (W/m³)",            limits=((t_hours[1], t_hours[end]), (0, zmax)))hidexdecorations!(axT;  grid=false)hidexdecorations!(axqᵛ; grid=false)hmT  = heatmap!(axT,  t_hours, z, T_zt;  colormap=:balance, colorrange=(-2, 2))hmqᵛ = heatmap!(axqᵛ, t_hours, z, qᵛ_zt; colormap=:dense, colorrange=(0, 1e-2))hmF  = heatmap!(axF,  t_hours, z, Fᴿ_zt; colormap=:balance, colorrange=(-0.1, 0.1))Colorbar(fig[1, 2], hmT;  label="T′ (K)")Colorbar(fig[2, 2], hmqᵛ; label="qᵛ (kg/kg)")Colorbar(fig[3, 2], hmF;  label="Fᴿ (W/m³)")fig

Animation of cloud structure

We animate xz slices of vertical velocity and cloud liquid water, zoomed to the lowest 5 km where the convective dynamics and clouds live.

wts  = FieldTimeSeries(slices_filename, "w")qᵛts = FieldTimeSeries(slices_filename, "qᵛ")times = wts.timesNt = length(times)wlim  = maximum(abs, wts) / 4qᵛlim = maximum(qᵛts) / 2fig = Figure(size=(1000, 600), fontsize=14)n = Observable(Nt)axw  = Axis(fig[1, 1]; xlabel="x (km)", ylabel="z (km)", title="w (m/s)", limits=(nothing, (0, 5e3)))axqᵛ = Axis(fig[1, 2]; xlabel="x (km)", ylabel="z (km)", title="qᵛ (kg/kg)", limits=(nothing, (0, 5e3)))title = @lift "Diurnal Radiative Convection at t = " * prettytime(times[$n])fig[0, :] = Label(fig, title, fontsize=16, tellwidth=false)wn  = @lift wts[$n]qᵛn = @lift qᵛts[$n]hmw  = heatmap!(axw,  wn;  colormap=:balance, colorrange=(-wlim, wlim))hmqᵛ = heatmap!(axqᵛ, qᵛn; colormap=:dense,   colorrange=(0, qᵛlim))Colorbar(fig[2, 1], hmw;  vertical=false, label="w (m/s)")Colorbar(fig[2, 2], hmqᵛ; vertical=false, label="qᵛ (g/kg)")hideydecorations!(axqᵛ; grid=false)CairoMakie.record(fig, "radiative_convection.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.13.2
Commit 45b6d15e34b (2026-10-14 20:33 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-20.1.8 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 5 GC (on 8 virtual cores)
Environment:
  JULIA_NUM_GC_THREADS = 4,1
  JULIA_GPG = 64B779A570972FFF7BFC2B54EAD471E1A1F2C10A
  JULIA_LOAD_PATH = :@breeze
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_VERSION = 1.13.2
  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.15.1 `.`
⌃ [052768ef] CUDA v6.4.1
  [13f3f980] CairoMakie v0.15.17
  [6a9e3e04] CloudMicrophysics v0.49.1
  [e30172f5] Documenter v1.20.0
  [daee34ce] DocumenterCitations v1.5.0
  [b6400b83] DocumenterCodeBlocks v1.5.1
  [1c52b33b] DocumenterLandingPage v0.2.2
  [7da242da] Enzyme v0.13.219
⌅ [46192b85] GPUArraysCore v0.2.1
⌅ [63c18a36] KernelAbstractions v0.9.47
  [98b081ad] Literate v2.21.1
  [85f8d34a] NCDatasets v0.14.15
  [cd8119b0] NumericalRadiation v0.2.1
  [9e8cae18] Oceananigans v0.118.1
  [a01a1ee8] RRTMGP v1.1.1
  [3c362404] Reactant v0.2.308
  [276daf66] SpecialFunctions v2.9.2
  [10745b16] Statistics v1.11.5
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.13.0
  [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.