Precipitating shallow cumulus convection (RICO)
This example simulates precipitating shallow cumulus convection following the Rain in Cumulus over the Ocean (RICO) intercomparison case (van Zanten et al., 2011). RICO is a canonical test case for large eddy simulations of trade-wind cumulus with active warm-rain microphysics.
The case is based on observations from the RICO field campaign conducted in the winter of 2004-2005 near Antigua and Barbuda in the Caribbean. Unlike BOMEX, which is non-precipitating, RICO produces drizzle and light rain from shallow cumulus clouds. The intercomparison study by van Zanten et al. (2011) brought together results from multiple large eddy simulation codes to establish benchmark statistics for precipitating shallow cumulus.
Initial and boundary conditions for this case are provided by the wonderfully useful package AtmosphericProfilesLibrary.jl. For precipitation we use the 1-moment scheme from CloudMicrophysics.jl, which provides prognostic rain mass with autoconversion and accretion processes.
using Breezeusing Oceananigans: Oceananigansusing Oceananigans.Unitsusing AtmosphericProfilesLibraryusing CairoMakieusing CloudMicrophysicsusing Printfusing Randomusing CUDARandom.seed!(42)if CUDA.functional() CUDA.seed!(42)endDomain and grid
The RICO domain is 12.8 km × 12.8 km horizontally with a vertical extent of 4 km (van Zanten et al., 2011). The intercomparison uses 128 × 128 × 100 grid points with 100 m horizontal resolution and 40 m vertical resolution.
Oceananigans.defaults.FloatType = Float32Nx = Ny = 128Nz = 100x = y = (0, 12800)z = (0, 4000)grid = RectilinearGrid(GPU(); x, y, z, size = (Nx, Ny, Nz), halo = (5, 5, 5), topology = (Periodic, Periodic, Bounded))128×128×100 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── Periodic x ∈ [0.0, 12800.0) regularly spaced with Δx=100.0
├── Periodic y ∈ [0.0, 12800.0) regularly spaced with Δy=100.0
└── Bounded z ∈ [0.0, 4000.0] regularly spaced with Δz=40.0Reference state and formulation
We use the anelastic formulation with a dry adiabatic reference state. The surface potential temperature $θ_0 = 297.9$ K and surface pressure $p_0 = 1015.4$ hPa are taken from van Zanten et al. (2011).
constants = ThermodynamicConstants()reference_state = ReferenceState(grid, constants, surface_pressure = 101540, potential_temperature = 297.9)dynamics = AnelasticDynamics(reference_state)AnelasticDynamics(p₀=101540.0, θ₀=297.9)
└── pressure_anomaly: not materializedSurface fluxes
Unlike the BOMEX protocol, which prescribes momentum, moisture, and thermodynamic fluxes, the RICO protocol decrees the computation of fluxes by bulk aerodynamic formulae with constant transfer coefficients (see van Zanten et al. (2011); text surrounding equations 1-4):
Cᴰ = 1.229e-3 # Drag coefficient for momentumCᵀ = 1.094e-3 # Sensible heat transfer coefficientCᵛ = 1.133e-3 # Moisture flux transfer coefficientT₀ = 299.8 # Sea surface temperature (K)299.8We implement the specified bulk formula with Breeze utilities whose scope currently extends only to constant coefficients (but could expand in the future),
ρe_flux = BulkSensibleHeatFlux(coefficient=Cᵀ, surface_temperature=T₀)ρqᵉ_flux = BulkVaporFlux(coefficient=Cᵛ, surface_temperature=T₀)ρe_bcs = FieldBoundaryConditions(bottom=ρe_flux)ρqᵉ_bcs = FieldBoundaryConditions(bottom=ρqᵉ_flux)ρu_bcs = FieldBoundaryConditions(bottom=BulkDrag(coefficient=Cᴰ))ρv_bcs = FieldBoundaryConditions(bottom=BulkDrag(coefficient=Cᴰ))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.001229, gustiness=0)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)Within the canon of Monin-Obukhov similarity theory, these transfer coefficients should be scaled if the vertical grid spacing is changed. Here we can use the values from van Zanten et al. (2011) verbatim because we use the recommended vertical grid spacing of 40 m.
Sponge layer
To prevent spurious wave reflections from the upper boundary, we add a Rayleigh damping sponge layer in the upper 500 m of the domain. The sponge damps vertical velocity toward zero using Oceananigans' Relaxation forcing with a GaussianMask.
sponge_rate = 1/8 # s⁻¹ - relaxation rate (8 s timescale)sponge_mask = GaussianMask{:z}(center=3500, width=500)sponge = Relaxation(rate=sponge_rate, mask=sponge_mask)Relaxation{Float64}
├── rate: 0.125
├── mask: exp(-(z - 3500)^2 / (2 * 500^2))
└── target: 0Large-scale subsidence
The RICO protocol includes large-scale subsidence that advects mean profiles downward. The subsidence velocity profile increases linearly to $-0.005$ m/s at 2260 m and remains constant above (van Zanten et al., 2011),
FT = eltype(grid)wˢ_profile = AtmosphericProfilesLibrary.Rico_subsidence(FT)wˢ = Field{Nothing, Nothing, Face}(grid)set!(wˢ, z -> wˢ_profile(z))subsidence = SubsidenceForcing(wˢ)SubsidenceForcing with wˢ: 1×1×101 Field{Nothing, Nothing, Oceananigans.Grids.Face} reduced over dims = (1, 2) on Oceananigans.Grids.RectilinearGrid on CUDAGPUThis is what it looks like:
lines(wˢ; axis = (xlabel = "wˢ (m/s)",))
Geostrophic forcing
The momentum equations include a Coriolis force with prescribed geostrophic wind. The RICO Coriolis parameter corresponds to latitude around 18°N: $f = 4.5 \times 10^{-5}$ s⁻¹.
coriolis = FPlane(f=4.5e-5)uᵍ = AtmosphericProfilesLibrary.Rico_geostrophic_ug(FT)vᵍ = AtmosphericProfilesLibrary.Rico_geostrophic_vg(FT)geostrophic = geostrophic_forcings(z -> uᵍ(z), z -> vᵍ(z))NamedTuple with 2 GeostrophicForcings:
├── u: GeostrophicForcing{XDirection}
│ └── geostrophic_velocity: #8 (generic function with 1 method)
└── v: GeostrophicForcing{YDirection}
└── geostrophic_velocity: #6 (generic function with 1 method)Moisture tendency
A prescribed large-scale moisture tendency represents the effects of advection by the large-scale circulation (van Zanten et al., 2011). We supply this tendency in specific form; Breeze multiplies by $ρ$ automatically at kernel time when the forcing is keyed under the specific variable name qᵉ.
∂t_qᵉ_large_scale = Field{Nothing, Nothing, Center}(grid)dqdt_profile = AtmosphericProfilesLibrary.Rico_dqtdt(FT)set!(∂t_qᵉ_large_scale, z -> dqdt_profile(z))qᵉ_large_scale_forcing = Forcing(∂t_qᵉ_large_scale)DiscreteForcing{Oceananigans.Fields.Field{Nothing, Nothing, Oceananigans.Grids.Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Architectures.GPU{CUDACore.CUDAKernels.CUDABackend}, Oceananigans.Grids.GridSize{128, 128, 100, 5, 5, 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.NoFluxBoundaryCondition, 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.NoFluxBoundaryCondition}}}, Nothing, Nothing}}
├── func: array_forcing_func (generic function with 1 method)
└── parameters: 1×1×100 Field{Nothing, Nothing, Oceananigans.Grids.Center} reduced over dims = (1, 2) on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── grid: 128×128×100 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 1×1×110 OffsetArray(::CUDACore.CuArray{Float32, 3, CUDACore.DeviceMemory}, 1:1, 1:1, -4:105) with eltype Float32 with indices 1:1×1:1×-4:105
└── max=4.0e-9, min=-1.14696e-8, mean=-1.80108e-9Radiative cooling
A prescribed radiative cooling profile is applied to the thermodynamic equation. The RICO case uses a constant radiative cooling rate of $-2.5$ K/day applied uniformly throughout the domain (van Zanten et al., 2011). This is the key simplification that allows us to avoid interactive radiation. We supply it under the specific potential-temperature key θ.
θ_large_scale_forcing = (x, y, z, t) -> - 2.5 / day#14 (generic function with 1 method)Assembling forcing and boundary conditions
Forcings are keyed under specific prognostic names (u, v, w, θ, qᵉ); Breeze applies the density factor $ρ$ automatically at kernel time. The Oceananigans Relaxation sponge damps the vertical velocity w toward zero.
forcing = (; u = (subsidence, geostrophic.u), v = (subsidence, geostrophic.v), w = sponge, qᵉ = (subsidence, qᵉ_large_scale_forcing), θ = (subsidence, θ_large_scale_forcing))boundary_conditions = (ρe=ρe_bcs, ρqᵉ=ρqᵉ_bcs, ρu=ρu_bcs, ρv=ρv_bcs)Model setup
We use one-moment bulk microphysics from CloudMicrophysics with cloud formatiom modeled with warm-phase saturationa adjustment and 5th-order WENO advection. The one-moment scheme prognoses rain density ρqʳ includes autoconversion (cloud liquid → rain) and accretion (cloud liquid swept up by falling rain) processes.
BreezeCloudMicrophysicsExt = Base.get_extension(Breeze, :BreezeCloudMicrophysicsExt)using .BreezeCloudMicrophysicsExt: OneMomentCloudMicrophysicscloud_formation = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())microphysics = OneMomentCloudMicrophysics(; cloud_formation)weno = WENO(order=5)bounds_preserving_weno = WENO(order=5, bounds=(0, 1))momentum_advection = wenoscalar_advection = (ρθ = weno, ρqᵉ = bounds_preserving_weno, ρqᶜˡ = bounds_preserving_weno, ρqʳ = bounds_preserving_weno)model = AtmosphereModel(grid; dynamics, coriolis, microphysics, momentum_advection, scalar_advection, forcing, boundary_conditions)AtmosphereModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×128×100 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── dynamics: AnelasticDynamics(p₀=101540.0, θ₀=297.9)
├── 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.0f0, 1.0f0))
│ └── ρqʳ: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5, bounds=(0.0f0, 1.0f0))
├── forcing: ρu=>MultipleForcings, ρv=>MultipleForcings, ρw=>SpecificForcing, ρθ=>MultipleForcings, ρqᵉ=>MultipleForcings
├── tracers: ()
├── coriolis: FPlane{Float32}(f=4.5e-5)
└── microphysics: BulkMicrophysicsInitial conditions
Mean profiles are specified as piecewise linear functions by van Zanten et al. (2011):
- Liquid-ice potential temperature $θ^{\ell i}(z)$
- Total water specific humidity $q^t(z)$
- Zonal velocity $u(z)$ and meridional velocity $v(z)$
The profiles are implemented in the wonderfully useful AtmosphericProfilesLibrary package developed by the Climate Modeling Alliance,
θˡⁱ₀ = AtmosphericProfilesLibrary.Rico_θ_liq_ice(FT)qᵗ₀ = AtmosphericProfilesLibrary.Rico_q_tot(FT)u₀ = AtmosphericProfilesLibrary.Rico_u(FT)v₀ = AtmosphericProfilesLibrary.Rico_v(FT)AtmosphericProfilesLibrary.ZProfile{AtmosphericProfilesLibrary.var"#Rico_v##0#Rico_v##1"{Float32}}(AtmosphericProfilesLibrary.var"#Rico_v##0#Rico_v##1"{Float32}())We add a small random perturbation below 1500 m to trigger convection.
zϵ = 1500 # mθᵢ(x, y, z) = θˡⁱ₀(z) + 1e-2 * (rand() - 0.5) * (z < zϵ)qᵢ(x, y, z) = qᵗ₀(z)uᵢ(x, y, z) = u₀(z)vᵢ(x, y, z) = v₀(z)set!(model, θ=θᵢ, qᵗ=qᵢ, u=uᵢ, v=vᵢ)Simulation
We run the simulation for 8 hours with adaptive time-stepping. RICO typically requires longer integration times than BOMEX to develop a quasi-steady precipitating state, and should be run for 24 hours. We choose 8 hours here to save computational costs in building the examples.
simulation = Simulation(model; Δt=2, stop_time=8hour)conjure_time_step_wizard!(simulation, cfl=0.7)Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)Output and progress
We set up a progress callback with hourly messages about interesting quantities,
θ = liquid_ice_potential_temperature(model)qˡ = model.microphysical_fields.qˡ # total liquid (cloud + rain)qᶜˡ = model.microphysical_fields.qᶜˡ # cloud liquid onlyqᵛ = model.microphysical_fields.qᵛqʳ = model.microphysical_fields.qʳ # rain mass fraction (diagnostic)ρqʳ = model.microphysical_fields.ρqʳρqʳ = model.microphysical_fields.ρqʳ # rain mass density (prognostic)# For keeping track of the computational expensewall_clock = Ref(time_ns())function progress(sim) qᶜˡmax = maximum(qᶜˡ) qʳmax = maximum(qʳ) qʳmin = minimum(qʳ) wmax = maximum(abs, model.velocities.w) elapsed = 1e-9 * (time_ns() - wall_clock[]) msg = @sprintf("Iter: %d, t: %s, Δt: %s, wall time: %s, max|w|: %.2e m/s", iteration(sim), prettytime(sim), prettytime(sim.Δt), prettytime(elapsed), wmax) msg *= @sprintf(", max(qᶜˡ): %.2e, extrema(qʳ): (%.2e, %.2e)", qᶜˡmax, qʳmin, qʳmax) @info msg return nothingendadd_callback!(simulation, progress, IterationInterval(1000))In addition to velocities, we output horizontal and time-averages of liquid water mass fraction (cloud and rain separately), specific humidity, and liquid-ice potential temperature,
# Precipitation rate diagnostic from one-moment microphysics# Integrals of precipitation rateP = precipitation_rate(model, :liquid)∫Pdz = Field(Integral(P, dims=3))u, v, w = model.velocitiesoutputs = merge(model.velocities, (; θ, qᶜˡ, qʳ, qᵛ, w² = w^2, uw = u*w, vw = v*w))averaged_outputs = NamedTuple(name => Average(outputs[name], dims=(1, 2)) for name in keys(outputs))filename = "rico.jld2"simulation.output_writers[:averages] = JLD2Writer(model, averaged_outputs; filename, schedule = AveragedTimeInterval(2hour), overwrite_existing = true)JLD2Writer scheduled on TimeInterval(2 hours):
├── filepath: rico.jld2
├── 10 outputs: (u, v, w, θ, qᶜˡ, qʳ, qᵛ, w², uw, vw) averaged on AveragedTimeInterval(window=2 hours, stride=1, interval=2 hours)
├── array_type: Array{Float32}
├── including: [:thermodynamic_constants]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)For an animation, we also output slices,
- xz-slices of qᶜˡ (cloud liquid) and qʳ (rain mass fraction)
- xy-slice of w (vertical velocity) with qˡ contours overlaid
w = model.velocities.wz = Oceananigans.Grids.znodes(grid, Center())k = searchsortedfirst(z, 1500) # cloud layer height for RICO@info "Saving xy slices at z = $(z[k]) m (k = $k)"slice_outputs = ( qᶜˡxz = view(qᶜˡ, :, 1, :), qʳxz = view(qʳ, :, 1, :), wxy = view(w, :, :, k), qˡxy = view(qˡ, :, :, k), qʳxy = view(qʳ, :, :, 1),)filename = "rico_slices.jld2"output_interval = 20secondssimulation.output_writers[:slices] = JLD2Writer(model, slice_outputs; filename, schedule = TimeInterval(output_interval), overwrite_existing = true)JLD2Writer scheduled on TimeInterval(20 seconds):
├── filepath: rico_slices.jld2
├── 5 outputs: (qᶜˡxz, qʳxz, wxy, qˡxy, qʳxy)
├── array_type: Array{Float32}
├── including: [:thermodynamic_constants]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)We're finally ready to run this thing,
run!(simulation)[ Info: Initializing simulation...
[ Info: Iter: 0, t: 0 seconds, Δt: 2.200 seconds, wall time: 1.074 minutes, max|w|: 0.00e+00 m/s, max(qᶜˡ): 0.00e+00, extrema(qʳ): (0.00e+00, 0.00e+00)
[ Info: ... simulation initialization complete (30.579 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (4.965 seconds).
[ Info: Iter: 1000, t: 58.046 minutes, Δt: 2.741 seconds, wall time: 2.098 minutes, max|w|: 5.06e+00 m/s, max(qᶜˡ): 3.50e-03, extrema(qʳ): (-2.35e-06, 8.27e-04)
[ Info: Iter: 2000, t: 1.898 hours, Δt: 3.765 seconds, wall time: 2.599 minutes, max|w|: 2.56e+00 m/s, max(qᶜˡ): 8.19e-04, extrema(qʳ): (-1.44e-07, 7.71e-06)
[ Info: Iter: 3000, t: 2.698 hours, Δt: 3.071 seconds, wall time: 3.059 minutes, max|w|: 4.61e+00 m/s, max(qᶜˡ): 2.17e-03, extrema(qʳ): (-1.96e-06, 1.84e-04)
[ Info: Iter: 4000, t: 3.462 hours, Δt: 3.100 seconds, wall time: 3.791 minutes, max|w|: 4.39e+00 m/s, max(qᶜˡ): 2.29e-03, extrema(qʳ): (-1.80e-06, 3.43e-04)
[ Info: Iter: 5000, t: 4.160 hours, Δt: 3.012 seconds, wall time: 4.507 minutes, max|w|: 4.46e+00 m/s, max(qᶜˡ): 2.30e-03, extrema(qʳ): (-1.90e-06, 3.74e-04)
[ Info: Iter: 6000, t: 4.862 hours, Δt: 2.429 seconds, wall time: 5.304 minutes, max|w|: 7.58e+00 m/s, max(qᶜˡ): 2.82e-03, extrema(qʳ): (-1.86e-06, 5.99e-04)
[ Info: Iter: 7000, t: 5.593 hours, Δt: 2.148 seconds, wall time: 5.876 minutes, max|w|: 8.41e+00 m/s, max(qᶜˡ): 2.93e-03, extrema(qʳ): (-1.82e-06, 3.18e-04)
[ Info: Iter: 8000, t: 6.217 hours, Δt: 1.932 seconds, wall time: 6.316 minutes, max|w|: 1.01e+01 m/s, max(qᶜˡ): 3.42e-03, extrema(qʳ): (-1.90e-06, 5.27e-04)
[ Info: Iter: 9000, t: 6.832 hours, Δt: 2.860 seconds, wall time: 6.714 minutes, max|w|: 5.69e+00 m/s, max(qᶜˡ): 2.94e-03, extrema(qʳ): (-1.78e-06, 6.45e-04)
[ Info: Iter: 10000, t: 7.517 hours, Δt: 2.731 seconds, wall time: 7.112 minutes, max|w|: 5.58e+00 m/s, max(qᶜˡ): 2.41e-03, extrema(qʳ): (-1.80e-06, 4.58e-04)
[ Info: Simulation is stopping after running for 6.690 minutes.
[ Info: Simulation time 8 hours equals or exceeds stop time 8 hours.
Results: mean profile evolution
We visualize the evolution of horizontally-averaged profiles every hour.
averages_filename = "rico.jld2"θts = FieldTimeSeries(averages_filename, "θ")qᵛts = FieldTimeSeries(averages_filename, "qᵛ")qᶜˡts = FieldTimeSeries(averages_filename, "qᶜˡ")qʳts = FieldTimeSeries(averages_filename, "qʳ")uts = FieldTimeSeries(averages_filename, "u")vts = FieldTimeSeries(averages_filename, "v")w²ts = FieldTimeSeries(averages_filename, "w²")uwts = FieldTimeSeries(averages_filename, "uw")vwts = FieldTimeSeries(averages_filename, "vw")fig = Figure(size=(1100, 700), fontsize=14)# Top row: θ, qᵛ, qᶜˡ/qʳaxθ = Axis(fig[1, 1], xlabel="θ (K)", ylabel="z (m)")axqᵛ = Axis(fig[1, 2], xlabel="qᵛ (kg/kg)", ylabel="z (m)")axqˡ = Axis(fig[1, 3], xlabel="qᶜˡ, qʳ (kg/kg)", ylabel="z (m)")# Bottom row: u/v, w², uw/vwaxuv = Axis(fig[2, 1], xlabel="u, v (m/s)", ylabel="z (m)")axw² = Axis(fig[2, 2], xlabel="w² (m²/s²)", ylabel="z (m)")axuw = Axis(fig[2, 3], xlabel="uw, vw (m²/s²)", ylabel="z (m)")times = θts.timesNt = length(times)default_colours = Makie.wong_colors()colors = [default_colours[mod1(i, length(default_colours))] for i in 1:Nt]for n in 1:Nt label = n == 1 ? "initial condition" : "mean over $(Int(times[n-1]/hour))-$(Int(times[n]/hour)) hr" # Top row lines!(axθ, θts[n], color=colors[n], label=label) lines!(axqᵛ, qᵛts[n], color=colors[n]) lines!(axqˡ, qᶜˡts[n], color=colors[n], linestyle=:solid) lines!(axqˡ, qʳts[n], color=colors[n], linestyle=:dash) # Bottom row lines!(axuv, uts[n], color=colors[n], linestyle=:solid) lines!(axuv, vts[n], color=colors[n], linestyle=:dash) lines!(axw², w²ts[n], color=colors[n]) lines!(axuw, uwts[n], color=colors[n], linestyle=:solid) lines!(axuw, vwts[n], color=colors[n], linestyle=:dash)endSet axis limits to focus on the boundary layer
for ax in (axθ, axqᵛ, axqˡ, axuv, axw², axuw) ylims!(ax, -100, 3500)endxlims!(axθ, 296, 318)xlims!(axqᵛ, 0, 1.8e-2)xlims!(axqˡ, -2e-6, 1.2e-5)xlims!(axuv, -12, 2)Add legends and annotations
axislegend(axθ, position=:rb)text!(axuv, -10, 2500, text="solid: u\ndashed: v", fontsize=14)text!(axqˡ, 1e-6, 2500, text="solid: qᶜˡ\ndashed: qʳ", fontsize=14)text!(axuw, 0.01, 2500, text="solid: uw\ndashed: vw", fontsize=14)fig[0, :] = Label(fig, "RICO: Horizontally-averaged profiles", fontsize=18, tellwidth=false)fig
The simulation shows the development of a cloudy, precipitating boundary layer with:
- Deeper cloud layer than BOMEX (tops reaching ~2.5-3 km)
- Higher moisture content supporting warm-rain processes
- Trade-wind flow with stronger westerlies
- Distinct profiles of cloud liquid (qᶜˡ) and rain (qʳ) as in van Zanten et al. (2011)
Animation: cloud structure and dynamics
We create a 4-panel animation showing:
- Top left: xz-slice of cloud liquid water qᶜˡ
- Top right: xz-slice of rain mass fraction qʳ
- Bottom: xy-slice of vertical velocity w with qˡ contours overlaid
wxy_ts = FieldTimeSeries("rico_slices.jld2", "wxy")qᶜˡxz_ts = FieldTimeSeries("rico_slices.jld2", "qᶜˡxz")qʳxz_ts = FieldTimeSeries("rico_slices.jld2", "qʳxz")qˡxy_ts = FieldTimeSeries("rico_slices.jld2", "qˡxy")qʳxy_ts = FieldTimeSeries("rico_slices.jld2", "qʳxy")times = wxy_ts.timesNt = length(times)qᶜˡlim = maximum(qᶜˡxz_ts) / 4qʳlim = maximum(qʳxz_ts) / 4wlim = maximum(abs, wxy_ts) / 24.421843f0Now let's plot the slices and animate them.
fig = Figure(size=(900, 850), fontsize=14)axqᶜˡxz = Axis(fig[2, 1], aspect=2, ylabel="z (m)", xaxisposition=:top)axqʳxz = Axis(fig[2, 2], aspect=2, ylabel="z (m)", yaxisposition=:right, xaxisposition=:top)axwxy = Axis(fig[3, 1], aspect=1, xlabel="x (m)", ylabel="y (m)")axqʳxy = Axis(fig[3, 2], aspect=1, xlabel="x (m)", ylabel="y (m)", yaxisposition=:right)hidexdecorations!(axqᶜˡxz)hidexdecorations!(axqʳxz)n = Observable(1)qᶜˡxz_n = @lift qᶜˡxz_ts[$n]qʳxz_n = @lift qʳxz_ts[$n]wxy_n = @lift wxy_ts[$n]qʳxy_n = @lift qʳxy_ts[$n]qˡxy_n = @lift qˡxy_ts[$n]qˡcontour = @lift maximum(qˡxy_ts[$n]) / 8 # threshold for cloud contourslevels = @lift [$qˡcontour]title = @lift @sprintf("Clouds, rain, and updrafts in RICO at t = %16.3f hours", times[$n] / hour)hmqᶜˡ = heatmap!(axqᶜˡxz, qᶜˡxz_n, colormap=:dense, colorrange=(0, qᶜˡlim))hmqʳ = heatmap!(axqʳxz, qʳxz_n, colormap=:amp, colorrange=(0, qʳlim))hmw = heatmap!(axwxy, wxy_n, colormap=:balance, colorrange=(-wlim, wlim))contour!(axwxy, qˡxy_n; levels, color=(:black, 0.3), linewidth=3)hmqʳ = heatmap!(axqʳxy, qʳxy_n, colormap=:amp, colorrange=(0, qʳlim))contour!(axqʳxy, qˡxy_n; levels, color=(:black, 0.3), linewidth=3)Colorbar(fig[1, 1], hmqᶜˡ, vertical=false, flipaxis=true, label="Cloud liquid qᶜˡ (x, y=0, z)")Colorbar(fig[1, 2], hmqʳ, vertical=false, flipaxis=true, label="Rain mass fraction qʳ (x, y=0, z)")Colorbar(fig[4, 1], hmw, vertical=false, flipaxis=false, label="Vertical velocity w (x, y, z=$(z[k])) with qˡ contours")Colorbar(fig[4, 2], hmqʳ, vertical=false, flipaxis=false, label="Rain mass fraction qʳ (x, y, z=0)")fig[0, :] = Label(fig, title, fontsize=18, tellwidth=false)rowgap!(fig.layout, 2, -60)rowgap!(fig.layout, 3, -80)n₁ = floor(Int, 6hours / output_interval)n₂ = ceil(Int, 8hours / output_interval)CairoMakie.record(fig, "rico_slices.mp4", n₁:n₂; framerate = 12, compression = 23) do nn n[] = nnendJulia 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.8
[660aa2fb] Breeze v0.8.0 `.`
[052768ef] CUDA v6.2.1
[13f3f980] CairoMakie v0.15.13
⌅ [6a9e3e04] CloudMicrophysics v0.37.1
[e30172f5] Documenter v1.17.0
[daee34ce] DocumenterCitations v1.4.1
[b6400b83] DocumenterCodeBlocks v1.1.0
[7da242da] Enzyme v0.13.198
[46192b85] GPUArraysCore v0.2.0
[63c18a36] KernelAbstractions v0.9.42
[98b081ad] Literate v2.21.0
[85f8d34a] NCDatasets v0.14.15
[9e8cae18] Oceananigans v0.110.14
⌅ [a01a1ee8] RRTMGP v0.21.9
[3c362404] Reactant v0.2.278
[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.