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 Breeze
using Oceananigans: Oceananigans
using Oceananigans.Units

using AtmosphericProfilesLibrary
using CairoMakie
using CloudMicrophysics
using Printf
using Random
using CUDA

Random.seed!(42)

Random.TaskLocalRNG()

Domain 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 = Float32

Nx = Ny = 128
Nz = 100

x = 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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.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.0

Reference 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 materialized

Surface 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 momentum
Cᵀ = 1.094e-3 # "Temperature" aka sensible heat transfer coefficient
Cᵛ = 1.133e-3 # Moisture flux transfer coefficient
T₀ = 299.8    # Sea surface temperature (K)

299.8

We implement the specified bulk formula with Breeze utilities whose scope currently extends only to constant coefficients (but could expand in the future),

ρθ_flux = BulkSensibleHeatFlux(coefficient=Cᵀ, surface_temperature=T₀)
ρqᵗ_flux = BulkVaporFlux(coefficient=Cᵛ, surface_temperature=T₀)

ρθ_bcs = FieldBoundaryConditions(bottom=ρθ_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, Oceananigans.Forcings.GaussianMask{:z, Int64}, typeof(Oceananigans.Forcings.zerofunction)}
├── rate: 0.125
├── mask: exp(-(z - 3500)^2 / (2 * 500^2))
└── target: 0

Large-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 CUDAGPU

This 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_momentum: #8 (generic function with 1 method)
└── ρv: GeostrophicForcing{YDirection}
    └── geostrophic_momentum: #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).

ρᵣ = reference_state.density
∂t_ρqᵗ_large_scale = Field{Nothing, Nothing, Center}(grid)
dqdt_profile = AtmosphericProfilesLibrary.Rico_dqtdt(FT)
set!(∂t_ρqᵗ_large_scale, z -> dqdt_profile(z))
set!(∂t_ρqᵗ_large_scale, ρᵣ * ∂t_ρqᵗ_large_scale)
∂t_ρ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{CUDA.CUDAKernels.CUDABackend}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 3, CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Nothing, Nothing, Nothing, Nothing, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Nothing, @NamedTuple{bottom_and_top::KernelAbstractions.Kernel{CUDA.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, typeof(Oceananigans.BoundaryConditions.gpu__fill_bottom_and_top_halo!)}, south_and_north::Nothing, west_and_east::Nothing}, @NamedTuple{bottom_and_top::Tuple{Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition}, south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Nothing, Nothing}}}, 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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.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(::CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}, 1:1, 1:1, -4:105) with eltype Float32 with indices 1:1×1:1×-4:105
    └── max=3.67374e-9, min=-1.35993e-8, mean=-2.32697e-9

Radiative 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.

∂t_ρθ_large_scale = Field{Nothing, Nothing, Center}(grid)
∂t_θ_large_scale = - 2.5 / day # K / day
set!(∂t_ρθ_large_scale, ρᵣ * ∂t_θ_large_scale)
ρθ_large_scale_forcing = Forcing(∂t_ρθ_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{CUDA.CUDAKernels.CUDABackend}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 3, CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Nothing, Nothing, Nothing, Nothing, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Nothing, @NamedTuple{bottom_and_top::KernelAbstractions.Kernel{CUDA.CUDAKernels.CUDABackend, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, typeof(Oceananigans.BoundaryConditions.gpu__fill_bottom_and_top_halo!)}, south_and_north::Nothing, west_and_east::Nothing}, @NamedTuple{bottom_and_top::Tuple{Oceananigans.BoundaryConditions.NoFluxBoundaryCondition, Oceananigans.BoundaryConditions.NoFluxBoundaryCondition}, south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Nothing, Nothing}}}, 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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.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(::CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}, 1:1, 1:1, -4:105) with eltype Float32 with indices 1:1×1:1×-4:105
    └── max=-2.42263e-5, min=-3.4308e-5, mean=-2.90939e-5

Assembling forcing and boundary conditions

Fρu = (subsidence, geostrophic.ρu)
Fρv = (subsidence, geostrophic.ρv)
Fρw = sponge
Fρqᵗ = (subsidence, ∂t_ρqᵗ_large_scale_forcing)
Fρθ = (subsidence, ρθ_large_scale_forcing)

forcing = (ρu=Fρu, ρv=Fρv, ρw=Fρw, ρqᵗ=Fρqᵗ, ρθ=Fρθ)
boundary_conditions = (ρθ=ρθ_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. This is a more physically-realistic representation of warm-rain precipitation than the zero-moment scheme.

BreezeCloudMicrophysicsExt = Base.get_extension(Breeze, :BreezeCloudMicrophysicsExt)
using .BreezeCloudMicrophysicsExt: OneMomentCloudMicrophysics

cloud_formation = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())
microphysics = OneMomentCloudMicrophysics(; cloud_formation)

weno = WENO(order=5)
bounds_preserving_weno = WENO(order=5, bounds=(0, 1))

momentum_advection = weno
scalar_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, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 5×5×5 halo
├── dynamics: AnelasticDynamics(p₀=101540.0, θ₀=297.9)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: 
│   ├── momentum: WENO{3, Float32, Float32}(order=5)
│   ├── ρθ: WENO{3, Float32, Float32}(order=5)
│   ├── ρqᵗ: WENO{3, Float32, Float32}(order=5)
│   └── ρqʳ: WENO{3, Float32, Float32}(order=5)
├── tracers: ()
├── coriolis: FPlane{Float32}(f=4.5e-5)
└── microphysics: BulkMicrophysics

Initial 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)

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 only
qᵛ = 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 expense
wall_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 nothing
end

add_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 rate
P = precipitation_rate(model, :liquid)
∫Pdz = Field(Integral(P, dims=3))

u, v, w = model.velocities
outputs = 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: [:grid]
├── file_splitting: NoFileSplitting
└── file size: 63.0 KiB

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.w

z = 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 = 20seconds
simulation.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: [:grid]
├── file_splitting: NoFileSplitting
└── file size: 41.7 KiB

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.376 minutes, max|w|: 6.72e-07 m/s, max(qᶜˡ): 0.00e+00, extrema(qʳ): (0.00e+00, 0.00e+00)
[ Info:     ... simulation initialization complete (32.028 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (6.641 seconds).
[ Info: Iter: 1000, t: 53.211 minutes, Δt: 2.584 seconds, wall time: 3.447 minutes, max|w|: 6.08e+00 m/s, max(qᶜˡ): 3.26e-03, extrema(qʳ): (-1.97e-06, 6.02e-04)
[ Info: Iter: 2000, t: 1.810 hours, Δt: 3.534 seconds, wall time: 4.861 minutes, max|w|: 3.39e+00 m/s, max(qᶜˡ): 2.00e-03, extrema(qʳ): (-1.67e-06, 1.12e-04)
[ Info: Iter: 3000, t: 2.591 hours, Δt: 2.998 seconds, wall time: 6.257 minutes, max|w|: 4.37e+00 m/s, max(qᶜˡ): 2.14e-03, extrema(qʳ): (-1.96e-06, 2.96e-04)
[ Info: Iter: 4000, t: 3.308 hours, Δt: 2.464 seconds, wall time: 7.654 minutes, max|w|: 7.09e+00 m/s, max(qᶜˡ): 2.64e-03, extrema(qʳ): (-1.93e-06, 4.28e-04)
[ Info: Iter: 5000, t: 4.008 hours, Δt: 2.963 seconds, wall time: 9.054 minutes, max|w|: 5.23e+00 m/s, max(qᶜˡ): 2.53e-03, extrema(qʳ): (-1.89e-06, 5.32e-04)
[ Info: Iter: 6000, t: 4.686 hours, Δt: 2.712 seconds, wall time: 10.454 minutes, max|w|: 5.75e+00 m/s, max(qᶜˡ): 2.19e-03, extrema(qʳ): (-1.84e-06, 3.42e-04)
[ Info: Iter: 7000, t: 5.351 hours, Δt: 2.535 seconds, wall time: 11.855 minutes, max|w|: 6.52e+00 m/s, max(qᶜˡ): 2.43e-03, extrema(qʳ): (-1.92e-06, 3.03e-04)
[ Info: Iter: 8000, t: 5.970 hours, Δt: 2.172 seconds, wall time: 13.255 minutes, max|w|: 8.73e+00 m/s, max(qᶜˡ): 3.38e-03, extrema(qʳ): (-1.79e-06, 4.93e-04)
[ Info: Iter: 9000, t: 6.570 hours, Δt: 2.162 seconds, wall time: 14.656 minutes, max|w|: 8.19e+00 m/s, max(qᶜˡ): 2.73e-03, extrema(qʳ): (-2.25e-06, 1.20e-03)
[ Info: Iter: 10000, t: 7.235 hours, Δt: 2.828 seconds, wall time: 16.062 minutes, max|w|: 5.52e+00 m/s, max(qᶜˡ): 2.73e-03, extrema(qʳ): (-1.92e-06, 5.03e-04)
[ Info: Iter: 11000, t: 7.880 hours, Δt: 2.785 seconds, wall time: 17.470 minutes, max|w|: 5.38e+00 m/s, max(qᶜˡ): 2.37e-03, extrema(qʳ): (-1.74e-06, 4.48e-04)
[ Info: Simulation is stopping after running for 16.838 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/vw
axuv = 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.times
Nt = 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)
end

Set axis limits to focus on the boundary layer

for ax in (axθ, axqᵛ, axqˡ, axuv, axw², axuw)
    ylims!(ax, -100, 3500)
end

xlims!(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)

save("rico_profiles.png", fig)
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.times
Nt = length(times)

qᶜˡlim = maximum(qᶜˡxz_ts) / 4
qʳlim = maximum(qʳxz_ts) / 4
wlim = maximum(abs, wxy_ts) / 2

4.436622f0

Now 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 contours
levels = @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) do nn
    n[] = nn
end

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 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_VERSION = 1.12.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 Pkg
Pkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
  [86bc3604] AtmosphericProfilesLibrary v0.1.7
  [660aa2fb] Breeze v0.2.1 `.`
  [052768ef] CUDA v5.9.5
  [13f3f980] CairoMakie v0.15.8
  [6a9e3e04] CloudMicrophysics v0.29.1
  [e30172f5] Documenter v1.16.1
  [daee34ce] DocumenterCitations v1.4.1
  [98b081ad] Literate v2.21.0
⌅ [9e8cae18] Oceananigans v0.102.5
  [a01a1ee8] RRTMGP v0.21.6
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.0
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.