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 RRTMGP radiation computes spectrally-resolved shortwave and longwave fluxes. 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 25 km for a realistic atmospheric column. A stratospheric 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 RRTMGP lookup tablesusing RRTMGPusing 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 25 km. This gives RRTMGP a realistic atmospheric column (including the stratosphere) while keeping the total cell count 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_anomaly: not materialized

Background atmosphere

RRTMGP requires 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.

@inline function tropical_ozone(z)    troposphere_O₃ = 30e-9 * (1 + 0.5 * z / 10_000)    zˢᵗ = 25e3    Hˢᵗ = 5e3    stratosphere_O₃ = 8e-6 * exp(-((z - zˢᵗ) / Hˢᵗ)^2)    χˢᵗ = 1 / (1 + exp(-(z - 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 RRTMGP 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{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.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 OffsetArray(::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°Nradiation = RadiativeTransferModel(grid, AllSkyOptics(), 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{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.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)

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 25 km, but the initial stratosphere isn'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 stratosphere 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 / constants.dry_air.molar_mass  # 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∂tendsponge = Forcing(stratospheric_relaxation; discrete_form=true,                 parameters=(; Tᵣ, ρᵣ, cᵖᵈ, τ=τ_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{OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 1, 51, 5, 0, 5}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 3, 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{OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 1, 51, 5, 0, 5}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 3, 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{OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 1, 51, 5, 0, 5}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 3, 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{OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}, OffsetArrays.OffsetVector{Float32, CUDACore.CuArray{Float32, 1, CUDACore.DeviceMemory}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 1, 51, 5, 0, 5}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 3, 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, τ::Float64}}
├── func: stratospheric_relaxation (generic function with 1 method)
└── parameters: (Tᵣ = 1×1×51 Field{Nothing, Nothing, Oceananigans.Grids.Center} reduced over dims = (1, 2) on Oceananigans.Grids.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 OffsetArray(::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, Oceananigans.Grids.Center} reduced over dims = (1, 2) on Oceananigans.Grids.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 OffsetArray(::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ᵖᵈ = 34691.06f0, τ = 21600.0),)

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 × 25 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 × 25 km
[ Info: Initial T range: 238.17134 – 300.12 K
[ Info: Initial qᵗ range: 0.0 – 12.206435 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: 35.666 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 (22.150 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (7.215 seconds).
[ Info: Iter:  1000, t: 1.325 hours, Δt:   3.8s, wall: 20.291 seconds, max|w|:  7.67 m/s, T: [217.9, 302.0] K, max(qˡ): 0.00e+00, Tₛ: 280.3 K, OLR: 350.4 W/m²
[ Info: Iter:  2000, t: 2.528 hours, Δt:   5.1s, wall: 11.561 seconds, max|w|: 12.28 m/s, T: [217.2, 301.9] K, max(qˡ): 0.00e+00, Tₛ: 280.2 K, OLR: 337.8 W/m²
[ Info: Iter:  3000, t: 3.595 hours, Δt:   4.3s, wall: 11.209 seconds, max|w|:  9.04 m/s, T: [217.2, 304.9] K, max(qˡ): 0.00e+00, Tₛ: 281.7 K, OLR: 333.4 W/m²
[ Info: Iter:  4000, t: 4.553 hours, Δt:   3.3s, wall: 10.823 seconds, max|w|: 12.16 m/s, T: [217.3, 305.1] K, max(qˡ): 0.00e+00, Tₛ: 284.3 K, OLR: 334.1 W/m²
[ Info: Iter:  5000, t: 5.603 hours, Δt:   3.7s, wall: 16.770 seconds, max|w|: 10.48 m/s, T: [217.2, 305.8] K, max(qˡ): 0.00e+00, Tₛ: 288.3 K, OLR: 337.0 W/m²
[ Info: Iter:  6000, t: 6.562 hours, Δt:   3.2s, wall: 17.933 seconds, max|w|: 14.17 m/s, T: [217.5, 306.9] K, max(qˡ): 0.00e+00, Tₛ: 292.6 K, OLR: 342.0 W/m²
[ Info: Iter:  7000, t: 7.550 hours, Δt:   3.9s, wall: 18.852 seconds, max|w|: 10.12 m/s, T: [217.6, 303.4] K, max(qˡ): 0.00e+00, Tₛ: 297.6 K, OLR: 349.8 W/m²
[ Info: Iter:  8000, t: 8.586 hours, Δt:   3.8s, wall: 19.153 seconds, max|w|: 11.71 m/s, T: [217.5, 304.3] K, max(qˡ): 0.00e+00, Tₛ: 303.1 K, OLR: 358.2 W/m²
[ Info: Iter:  9000, t: 9.559 hours, Δt:   3.5s, wall: 15.791 seconds, max|w|: 14.01 m/s, T: [217.5, 304.3] K, max(qˡ): 9.50e-04, Tₛ: 307.9 K, OLR: 364.6 W/m²
[ Info: Iter: 10000, t: 10.637 hours, Δt:   3.8s, wall: 11.822 seconds, max|w|: 10.21 m/s, T: [217.2, 304.9] K, max(qˡ): 0.00e+00, Tₛ: 312.7 K, OLR: 372.5 W/m²
[ Info: Iter: 11000, t: 11.646 hours, Δt:   3.2s, wall: 11.460 seconds, max|w|: 16.91 m/s, T: [217.3, 305.1] K, max(qˡ): 5.48e-04, Tₛ: 316.3 K, OLR: 377.1 W/m²
[ Info: Iter: 12000, t: 12.643 hours, Δt:   4.1s, wall: 11.588 seconds, max|w|: 10.72 m/s, T: [217.5, 304.3] K, max(qˡ): 0.00e+00, Tₛ: 318.8 K, OLR: 380.7 W/m²
[ Info: Iter: 13000, t: 13.581 hours, Δt:   2.8s, wall: 11.809 seconds, max|w|: 14.10 m/s, T: [217.4, 303.3] K, max(qˡ): 0.00e+00, Tₛ: 319.9 K, OLR: 380.5 W/m²
[ Info: Iter: 14000, t: 14.549 hours, Δt:   3.4s, wall: 11.027 seconds, max|w|: 12.48 m/s, T: [217.5, 303.9] K, max(qˡ): 1.15e-03, Tₛ: 319.8 K, OLR: 378.6 W/m²
[ Info: Iter: 15000, t: 15.593 hours, Δt:   3.4s, wall: 11.166 seconds, max|w|: 14.21 m/s, T: [217.5, 305.3] K, max(qˡ): 0.00e+00, Tₛ: 318.3 K, OLR: 373.3 W/m²
[ Info: Iter: 16000, t: 16.474 hours, Δt:   3.9s, wall: 10.477 seconds, max|w|:  9.32 m/s, T: [217.6, 305.4] K, max(qˡ): 0.00e+00, Tₛ: 316.0 K, OLR: 367.3 W/m²
[ Info: Iter: 17000, t: 17.479 hours, Δt:   3.8s, wall: 10.962 seconds, max|w|: 11.93 m/s, T: [217.4, 303.0] K, max(qˡ): 7.14e-04, Tₛ: 312.3 K, OLR: 357.6 W/m²
[ Info: Iter: 18000, t: 18.486 hours, Δt:   3.1s, wall: 10.701 seconds, max|w|: 12.00 m/s, T: [217.2, 302.7] K, max(qˡ): 3.94e-04, Tₛ: 307.7 K, OLR: 349.3 W/m²
[ Info: Iter: 19000, t: 19.439 hours, Δt:   2.9s, wall: 11.574 seconds, max|w|: 13.91 m/s, T: [217.6, 303.8] K, max(qˡ): 1.55e-03, Tₛ: 302.9 K, OLR: 339.9 W/m²
[ Info: Iter: 20000, t: 20.399 hours, Δt:   3.8s, wall: 5.602 seconds, max|w|: 11.80 m/s, T: [217.3, 302.7] K, max(qˡ): 5.32e-04, Tₛ: 297.9 K, OLR: 331.9 W/m²
[ Info: Iter: 21000, t: 21.476 hours, Δt:   3.9s, wall: 5.144 seconds, max|w|: 12.51 m/s, T: [217.6, 303.9] K, max(qˡ): 1.20e-03, Tₛ: 292.5 K, OLR: 323.6 W/m²
[ Info: Iter: 22000, t: 22.485 hours, Δt:   4.3s, wall: 5.074 seconds, max|w|: 10.94 m/s, T: [217.8, 304.9] K, max(qˡ): 3.15e-04, Tₛ: 287.9 K, OLR: 315.9 W/m²
[ Info: Iter: 23000, t: 23.364 hours, Δt:   3.6s, wall: 4.961 seconds, max|w|: 11.54 m/s, T: [217.1, 303.0] K, max(qˡ): 6.01e-04, Tₛ: 284.6 K, OLR: 309.8 W/m²
[ Info: Iter: 24000, t: 1.017 days, Δt:   4.7s, wall: 5.109 seconds, max|w|: 10.29 m/s, T: [217.5, 304.1] K, max(qˡ): 0.00e+00, Tₛ: 281.7 K, OLR: 305.4 W/m²
[ Info: Iter: 25000, t: 1.058 days, Δt:   4.4s, wall: 5.179 seconds, max|w|: 13.66 m/s, T: [217.7, 304.8] K, max(qˡ): 0.00e+00, Tₛ: 280.3 K, OLR: 302.5 W/m²
[ Info: Iter: 26000, t: 1.100 days, Δt:   4.6s, wall: 5.499 seconds, max|w|:  8.79 m/s, T: [217.4, 303.3] K, max(qˡ): 0.00e+00, Tₛ: 280.1 K, OLR: 301.2 W/m²
[ Info: Iter: 27000, t: 1.140 days, Δt:   3.2s, wall: 5.090 seconds, max|w|: 11.08 m/s, T: [217.6, 304.5] K, max(qˡ): 0.00e+00, Tₛ: 281.3 K, OLR: 301.6 W/m²
[ Info: Iter: 28000, t: 1.177 days, Δt:   3.3s, wall: 4.987 seconds, max|w|: 12.82 m/s, T: [217.4, 304.1] K, max(qˡ): 0.00e+00, Tₛ: 283.4 K, OLR: 303.7 W/m²
[ Info: Iter: 29000, t: 1.225 days, Δt:   4.4s, wall: 5.361 seconds, max|w|:  8.83 m/s, T: [217.7, 302.7] K, max(qˡ): 0.00e+00, Tₛ: 287.4 K, OLR: 307.9 W/m²
[ Info: Iter: 30000, t: 1.272 days, Δt:   4.0s, wall: 6.496 seconds, max|w|: 10.35 m/s, T: [217.5, 305.4] K, max(qˡ): 3.92e-04, Tₛ: 292.5 K, OLR: 314.2 W/m²
[ Info: Iter: 31000, t: 1.314 days, Δt:   4.5s, wall: 5.379 seconds, max|w|:  8.37 m/s, T: [217.6, 303.4] K, max(qˡ): 0.00e+00, Tₛ: 297.6 K, OLR: 320.6 W/m²
[ Info: Iter: 32000, t: 1.359 days, Δt:   4.8s, wall: 5.435 seconds, max|w|: 10.50 m/s, T: [217.3, 303.7] K, max(qˡ): 0.00e+00, Tₛ: 303.3 K, OLR: 327.9 W/m²
[ Info: Iter: 33000, t: 1.400 days, Δt:   3.4s, wall: 5.404 seconds, max|w|: 14.78 m/s, T: [217.6, 306.8] K, max(qˡ): 1.77e-04, Tₛ: 308.1 K, OLR: 335.4 W/m²
[ Info: Iter: 34000, t: 1.442 days, Δt:   3.0s, wall: 5.320 seconds, max|w|: 13.75 m/s, T: [217.6, 303.9] K, max(qˡ): 0.00e+00, Tₛ: 312.6 K, OLR: 341.8 W/m²
[ Info: Iter: 35000, t: 1.483 days, Δt:   3.9s, wall: 5.351 seconds, max|w|: 11.31 m/s, T: [217.6, 303.9] K, max(qˡ): 6.65e-04, Tₛ: 316.2 K, OLR: 346.5 W/m²
[ Info: Iter: 36000, t: 1.523 days, Δt:   4.1s, wall: 5.181 seconds, max|w|: 11.73 m/s, T: [217.6, 304.9] K, max(qˡ): 1.05e-03, Tₛ: 318.6 K, OLR: 349.6 W/m²
[ Info: Iter: 37000, t: 1.563 days, Δt:   3.6s, wall: 5.294 seconds, max|w|: 13.07 m/s, T: [217.7, 304.2] K, max(qˡ): 2.42e-04, Tₛ: 319.8 K, OLR: 349.8 W/m²
[ Info: Iter: 38000, t: 1.607 days, Δt:   3.4s, wall: 5.332 seconds, max|w|: 13.43 m/s, T: [217.6, 305.0] K, max(qˡ): 9.58e-04, Tₛ: 319.8 K, OLR: 347.8 W/m²
[ Info: Iter: 39000, t: 1.653 days, Δt:   4.5s, wall: 5.602 seconds, max|w|:  8.74 m/s, T: [217.5, 304.1] K, max(qˡ): 2.14e-04, Tₛ: 318.1 K, OLR: 343.7 W/m²
[ Info: Iter: 40000, t: 1.696 days, Δt:   3.3s, wall: 5.339 seconds, max|w|: 15.04 m/s, T: [217.3, 306.0] K, max(qˡ): 2.48e-03, Tₛ: 315.2 K, OLR: 336.8 W/m²
[ Info: Iter: 41000, t: 1.736 days, Δt:   3.2s, wall: 5.230 seconds, max|w|: 12.49 m/s, T: [217.3, 305.1] K, max(qˡ): 1.27e-03, Tₛ: 311.5 K, OLR: 330.8 W/m²
[ Info: Iter: 42000, t: 1.776 days, Δt:   3.2s, wall: 6.352 seconds, max|w|: 11.31 m/s, T: [217.3, 304.9] K, max(qˡ): 1.39e-03, Tₛ: 307.0 K, OLR: 324.1 W/m²
[ Info: Iter: 43000, t: 1.811 days, Δt:   2.7s, wall: 5.121 seconds, max|w|: 16.46 m/s, T: [217.6, 303.9] K, max(qˡ): 2.07e-03, Tₛ: 302.8 K, OLR: 318.7 W/m²
[ Info: Iter: 44000, t: 1.847 days, Δt:   2.8s, wall: 4.972 seconds, max|w|: 16.54 m/s, T: [217.5, 304.9] K, max(qˡ): 3.41e-03, Tₛ: 298.3 K, OLR: 312.6 W/m²
[ Info: Iter: 45000, t: 1.889 days, Δt:   3.5s, wall: 5.110 seconds, max|w|: 12.52 m/s, T: [217.3, 302.5] K, max(qˡ): 2.64e-03, Tₛ: 293.2 K, OLR: 305.9 W/m²
[ Info: Iter: 46000, t: 1.931 days, Δt:   3.8s, wall: 5.281 seconds, max|w|: 14.15 m/s, T: [217.6, 304.8] K, max(qˡ): 1.99e-03, Tₛ: 288.5 K, OLR: 300.1 W/m²
[ Info: Iter: 47000, t: 1.976 days, Δt:   4.2s, wall: 5.192 seconds, max|w|: 12.32 m/s, T: [217.4, 304.9] K, max(qˡ): 0.00e+00, Tₛ: 284.4 K, OLR: 294.9 W/m²
[ Info: Iter: 48000, t: 2.018 days, Δt:   5.0s, wall: 5.144 seconds, max|w|:  8.17 m/s, T: [217.5, 303.3] K, max(qˡ): 1.40e-03, Tₛ: 281.6 K, OLR: 291.4 W/m²
[ Info: Iter: 49000, t: 2.062 days, Δt:   4.5s, wall: 5.139 seconds, max|w|: 11.70 m/s, T: [217.6, 301.7] K, max(qˡ): 1.08e-03, Tₛ: 280.2 K, OLR: 289.6 W/m²
[ Info: Iter: 50000, t: 2.102 days, Δt:   4.8s, wall: 5.106 seconds, max|w|:  8.18 m/s, T: [217.6, 303.7] K, max(qˡ): 7.37e-05, Tₛ: 280.1 K, OLR: 289.0 W/m²
[ Info: Iter: 51000, t: 2.148 days, Δt:   5.0s, wall: 5.197 seconds, max|w|: 10.82 m/s, T: [217.6, 302.2] K, max(qˡ): 9.70e-04, Tₛ: 281.6 K, OLR: 289.5 W/m²
[ Info: Iter: 52000, t: 2.193 days, Δt:   2.9s, wall: 5.148 seconds, max|w|: 15.09 m/s, T: [217.5, 302.2] K, max(qˡ): 1.77e-03, Tₛ: 284.6 K, OLR: 291.7 W/m²
[ Info: Iter: 53000, t: 2.233 days, Δt:   2.7s, wall: 5.100 seconds, max|w|: 16.09 m/s, T: [217.2, 302.1] K, max(qˡ): 3.82e-03, Tₛ: 288.2 K, OLR: 295.2 W/m²
[ Info: Iter: 54000, t: 2.272 days, Δt:   4.0s, wall: 5.105 seconds, max|w|: 10.74 m/s, T: [217.4, 302.1] K, max(qˡ): 2.14e-03, Tₛ: 292.5 K, OLR: 298.9 W/m²
[ Info: Iter: 55000, t: 2.316 days, Δt:   4.3s, wall: 5.462 seconds, max|w|:  9.02 m/s, T: [217.3, 305.1] K, max(qˡ): 1.04e-03, Tₛ: 297.8 K, OLR: 305.1 W/m²
[ Info: Iter: 56000, t: 2.357 days, Δt:   3.2s, wall: 5.203 seconds, max|w|: 11.34 m/s, T: [217.5, 302.0] K, max(qˡ): 2.13e-03, Tₛ: 303.0 K, OLR: 310.9 W/m²
[ Info: Iter: 57000, t: 2.399 days, Δt:   4.6s, wall: 5.347 seconds, max|w|: 10.50 m/s, T: [217.7, 303.7] K, max(qˡ): 6.67e-04, Tₛ: 308.0 K, OLR: 317.3 W/m²
[ Info: Iter: 58000, t: 2.439 days, Δt:   2.9s, wall: 5.334 seconds, max|w|: 17.62 m/s, T: [217.8, 305.0] K, max(qˡ): 1.83e-03, Tₛ: 312.3 K, OLR: 321.6 W/m²
[ Info: Iter: 59000, t: 2.481 days, Δt:   3.1s, wall: 5.327 seconds, max|w|:  9.17 m/s, T: [217.7, 303.9] K, max(qˡ): 9.26e-04, Tₛ: 316.0 K, OLR: 326.5 W/m²
[ Info: Iter: 60000, t: 2.521 days, Δt:   3.6s, wall: 5.373 seconds, max|w|: 12.44 m/s, T: [217.5, 304.4] K, max(qˡ): 3.59e-03, Tₛ: 318.5 K, OLR: 328.8 W/m²
[ Info: Iter: 61000, t: 2.560 days, Δt:   3.9s, wall: 5.215 seconds, max|w|: 12.77 m/s, T: [217.6, 303.6] K, max(qˡ): 2.90e-03, Tₛ: 319.8 K, OLR: 329.2 W/m²
[ Info: Iter: 62000, t: 2.598 days, Δt:   3.0s, wall: 5.278 seconds, max|w|: 15.61 m/s, T: [217.6, 304.9] K, max(qˡ): 3.01e-03, Tₛ: 319.9 K, OLR: 328.3 W/m²
[ Info: Iter: 63000, t: 2.642 days, Δt:   3.8s, wall: 6.522 seconds, max|w|:  9.52 m/s, T: [217.5, 304.6] K, max(qˡ): 8.61e-04, Tₛ: 318.6 K, OLR: 326.0 W/m²
[ Info: Iter: 64000, t: 2.686 days, Δt:   3.0s, wall: 11.111 seconds, max|w|: 12.88 m/s, T: [217.6, 304.9] K, max(qˡ): 6.96e-04, Tₛ: 316.0 K, OLR: 321.5 W/m²
[ Info: Iter: 65000, t: 2.726 days, Δt:   3.8s, wall: 9.186 seconds, max|w|: 10.19 m/s, T: [217.7, 302.6] K, max(qˡ): 5.82e-04, Tₛ: 312.5 K, OLR: 317.4 W/m²
[ Info: Iter: 66000, t: 2.770 days, Δt:   3.1s, wall: 9.181 seconds, max|w|: 12.93 m/s, T: [217.4, 303.5] K, max(qˡ): 3.49e-03, Tₛ: 307.8 K, OLR: 310.5 W/m²
[ Info: Iter: 67000, t: 2.807 days, Δt:   2.9s, wall: 8.621 seconds, max|w|: 14.07 m/s, T: [217.4, 301.6] K, max(qˡ): 2.10e-03, Tₛ: 303.3 K, OLR: 304.9 W/m²
[ Info: Iter: 68000, t: 2.850 days, Δt:   3.9s, wall: 8.268 seconds, max|w|: 10.98 m/s, T: [217.5, 303.0] K, max(qˡ): 1.13e-03, Tₛ: 297.9 K, OLR: 301.0 W/m²
[ Info: Iter: 69000, t: 2.889 days, Δt:   3.2s, wall: 9.011 seconds, max|w|:  9.16 m/s, T: [217.4, 302.8] K, max(qˡ): 2.04e-03, Tₛ: 293.1 K, OLR: 296.0 W/m²
[ Info: Iter: 70000, t: 2.928 days, Δt:   4.1s, wall: 8.840 seconds, max|w|: 11.93 m/s, T: [217.4, 303.8] K, max(qˡ): 4.41e-04, Tₛ: 288.7 K, OLR: 293.0 W/m²
[ Info: Iter: 71000, t: 2.970 days, Δt:   4.0s, wall: 9.062 seconds, max|w|: 11.59 m/s, T: [217.7, 301.5] K, max(qˡ): 6.22e-04, Tₛ: 284.9 K, OLR: 289.6 W/m²
[ Info: Simulation is stopping after running for 9.821 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.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.10
  [660aa2fb] Breeze v0.11.3 `.`
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.15
⌅ [6a9e3e04] CloudMicrophysics v0.40.1
  [e30172f5] Documenter v1.19.0
  [daee34ce] DocumenterCitations v1.5.0
  [b6400b83] DocumenterCodeBlocks v1.5.0
  [1c52b33b] DocumenterLandingPage v0.2.2
  [7da242da] Enzyme v0.13.204
⌅ [46192b85] GPUArraysCore v0.2.0
  [63c18a36] KernelAbstractions v0.9.42
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.113.1
  [a01a1ee8] RRTMGP v1.0.0
  [3c362404] Reactant v0.2.288
  [276daf66] SpecialFunctions v2.9.0
  [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.