Shallow cumulus convection (BOMEX)
This example simulates shallow cumulus convection following the Barbados Oceanographic and Meteorological Experiment (BOMEX) intercomparison case (Siebesma et al., 2003). BOMEX is a canonical test case for large eddy simulations of shallow cumulus convection over a subtropical ocean.
The case is based on observations from the Barbados Oceanographic and Meteorological Experiment, which documented the structure and organization of trade-wind cumulus clouds. The intercomparison study by Siebesma et al. (2003) brought together results from 10 different large eddy simulation codes to establish benchmark statistics.
Initial and boundary conditions for this case are provided by the wonderfully useful package AtmosphericProfilesLibrary.jl.
using Breeze
using Oceananigans: Oceananigans
using Oceananigans.Units
using AtmosphericProfilesLibrary
using CairoMakie
using CUDA
using Printf
using Random
Random.seed!(938)Random.TaskLocalRNG()Domain and grid
The BOMEX domain is 6.4 km × 6.4 km horizontally with a vertical extent of 3 km (Siebesma et al. (2003); Section 3a). The intercomparison uses 64 × 64 × 75 grid points with 100 m horizontal resolution and 40 m vertical resolution.
For this documentation example, we reduce the numerical precision to Float32. This yields a 10x speed up on an NVidia T4 (which is used to build the docs).
Oceananigans.defaults.FloatType = Float32
Nx = Ny = 64
Nz = 75
x = y = (0, 6400)
z = (0, 3000)
grid = RectilinearGrid(GPU(); x, y, z,
size = (Nx, Ny, Nz), halo = (5, 5, 5),
topology = (Periodic, Periodic, Bounded))64×64×75 RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 5×5×5 halo
├── Periodic x ∈ [0.0, 6400.0) regularly spaced with Δx=100.0
├── Periodic y ∈ [0.0, 6400.0) regularly spaced with Δy=100.0
└── Bounded z ∈ [0.0, 3000.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 = 299.1$ K and surface pressure $p_0 = 1015$ hPa are taken from Siebesma et al. (2003); Appendix B.
constants = ThermodynamicConstants()
reference_state = ReferenceState(grid, constants,
surface_pressure = 101500,
potential_temperature = 299.1)
dynamics = AnelasticDynamics(reference_state)AnelasticDynamics(p₀=101500.0, θ₀=299.1)
└── pressure_anomaly: not materializedSurface fluxes
BOMEX prescribes constant surface sensible and latent heat fluxes (Siebesma et al. (2003), Appendix B, after Eq. B4):
- Sensible heat flux: $\overline{w'\theta'}|_0 = 8 \times 10^{-3}$ K m/s
- Moisture flux: $\overline{w'q_t'}|_0 = 5.2 \times 10^{-5}$ m/s
(Siebesma et al. (2003) refers to the moisture flux as the "latent heat flux". We convert these kinematic fluxes to mass fluxes by multiplying by surface density, which we estimate for a dry state using the pressure and temperature at $z=0$.
w′θ′ = 8e-3 # K m/s (sensible heat flux)
w′qᵗ′ = 5.2e-5 # m/s (moisture flux)
FT = eltype(grid)
p₀ = reference_state.surface_pressure
θ₀ = reference_state.potential_temperature
q₀ = Breeze.Thermodynamics.MoistureMassFractions{FT} |> zero
ρ₀ = Breeze.Thermodynamics.density(p₀, θ₀, q₀, constants)
ρθ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(ρ₀ * w′θ′))
ρqᵗ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(ρ₀ * w′qᵗ′))Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: FluxBoundaryCondition: 6.14847e-5
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)Surface momentum flux (drag)
A bulk drag parameterization is applied with friction velocity $u_* = 0.28$ m/s (Siebesma et al. (2003); Appendix B, after Eq. B4).
u★ = 0.28 # m/s
@inline ρu_drag(x, y, t, ρu, ρv, p) = - p.ρ₀ * p.u★^2 * ρu / sqrt(ρu^2 + ρv^2)
@inline ρv_drag(x, y, t, ρu, ρv, p) = - p.ρ₀ * p.u★^2 * ρv / sqrt(ρu^2 + ρv^2)
ρu_drag_bc = FluxBoundaryCondition(ρu_drag, field_dependencies=(:ρu, :ρv), parameters=(; ρ₀, u★))
ρv_drag_bc = FluxBoundaryCondition(ρv_drag, field_dependencies=(:ρu, :ρv), parameters=(; ρ₀, u★))
ρu_bcs = FieldBoundaryConditions(bottom=ρu_drag_bc)
ρv_bcs = FieldBoundaryConditions(bottom=ρv_drag_bc)Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── east: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── south: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── north: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
├── bottom: FluxBoundaryCondition: ContinuousBoundaryFunction ρv_drag at (Nothing, Nothing, Nothing)
├── top: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)
└── immersed: DefaultBoundaryCondition (FluxBoundaryCondition: Nothing)Large-scale subsidence
The BOMEX case includes large-scale subsidence that advects mean profiles downward. The subsidence velocity profile is prescribed by Siebesma et al. (2003); Appendix B, Eq. B5:
\[w^s(z) = \begin{cases} W^s \frac{z}{z_1} & z \le z_1 \\ W^s \left ( 1 - \frac{z - z_1}{z_2 - z_1} \right ) & z_1 < z \le z_2 \\ 0 & z > z_2 \end{cases}\]
where $W^s = -6.5 \times 10^{-3}$ m/s (note the negative sign for "subisdence"), $z_1 = 1500$ m and $z_2 = 2100$ m.
The subsidence velocity profile is provided by AtmosphericProfilesLibrary,
wˢ = Field{Nothing, Nothing, Face}(grid)
wˢ_profile = AtmosphericProfilesLibrary.Bomex_subsidence(FT)
set!(wˢ, z -> wˢ_profile(z))1×1×76 Field{Nothing, Nothing, Oceananigans.Grids.Face} reduced over dims = (1, 2) on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── grid: 64×64×75 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: Nothing, top: Nothing, immersed: Nothing
└── data: 1×1×86 OffsetArray(::CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}, 1:1, 1:1, -4:81) with eltype Float32 with indices 1:1×1:1×-4:81
└── max=0.0, min=-0.00641333, mean=-0.00224478and looks like:
lines(wˢ; axis = (xlabel = "wˢ (m/s)",))
Subsidence is implemented as an advection of the horizontally-averaged prognostic variables. This implementation –- which requires building Fields to represent horizontal averages and computing it every time step –- is handled by SubsidenceForcing.
subsidence = SubsidenceForcing(wˢ)SubsidenceForcing with wˢ: 1×1×76 Field{Nothing, Nothing, Oceananigans.Grids.Face} reduced over dims = (1, 2) on Oceananigans.Grids.RectilinearGrid on CUDAGPUGeostrophic forcing
The momentum equations include a Coriolis force with prescribed geostrophic wind. The geostrophic wind profiles are given by Siebesma et al. (2003); Appendix B, Eq. B6. Using geostrophic_forcings, we specify the geostrophic velocity profiles as functions of height, and the forcing is automatically materialized with the model's coriolis parameter and reference density.
coriolis = FPlane(f=3.76e-5)
uᵍ = AtmosphericProfilesLibrary.Bomex_geostrophic_u(FT)
vᵍ = AtmosphericProfilesLibrary.Bomex_geostrophic_v(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 (drying)
A prescribed large-scale drying tendency removes moisture above the cloud layer (Siebesma et al. (2003); Appendix B, Eq. B4). This represents the effects of advection by the large-scale circulation.
ρᵣ = reference_state.density
drying = Field{Nothing, Nothing, Center}(grid)
dqdt_profile = AtmosphericProfilesLibrary.Bomex_dqtdt(FT)
set!(drying, z -> dqdt_profile(z))
set!(drying, ρᵣ * drying)
ρqᵗ_drying_forcing = Forcing(drying)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×75 Field{Nothing, Nothing, Oceananigans.Grids.Center} reduced over dims = (1, 2) on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── grid: 64×64×75 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×85 OffsetArray(::CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}, 1:1, 1:1, -4:80) with eltype Float32 with indices 1:1×1:1×-4:80
└── max=0.0, min=-1.41656e-8, mean=-1.86054e-9Radiative cooling
A prescribed radiative cooling profile is applied to the thermodynamic equation (Siebesma et al. (2003); Appendix B, Eq. B3). Below the inversion, radiative cooling of about 2 K/day counteracts the surface heating. We use an energy forcing for radiation to ensure that it is applied to the potential temperature conservation equation consistently (see below for some elaboration about that).
Fρe_field = Field{Nothing, Nothing, Center}(grid)
cᵖᵈ = constants.dry_air.heat_capacity
dTdt_bomex = AtmosphericProfilesLibrary.Bomex_dTdt(FT)
set!(Fρe_field, z -> dTdt_bomex(1, z))
set!(Fρe_field, ρᵣ * cᵖᵈ * Fρe_field)
ρe_radiation_forcing = Forcing(Fρe_field)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×75 Field{Nothing, Nothing, Oceananigans.Grids.Center} reduced over dims = (1, 2) on Oceananigans.Grids.RectilinearGrid on CUDAGPU
├── grid: 64×64×75 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×85 OffsetArray(::CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}, 1:1, 1:1, -4:80) with eltype Float32 with indices 1:1×1:1×-4:80
└── max=-0.000283938, min=-0.0274623, mean=-0.0187429Assembling all the forcings
We build tuples of forcings for all the variables. Note that forcing functions are provided for both ρθ and ρe, which both contribute to the tendency of ρθ in different ways. In particular, the tendency for ρθ is written
\[∂_t (ρ θ) = - \boldsymbol{\nabla \cdot} \, ( ρ \boldsymbol{u} θ ) + F_{ρθ} + \frac{1}{cᵖᵐ Π} F_{ρ e} + \cdots\]
where $F_{ρ e}$ denotes the forcing function provided for ρe (e.g. for "energy density"), $F_{ρθ}$ denotes the forcing function provided for ρθ, and the $\cdots$ denote additional terms.
The geostrophic forcing provides both ρu and ρv components, which we merge with the subsidence forcing.
ρu_forcing = (subsidence, geostrophic.ρu)
ρv_forcing = (subsidence, geostrophic.ρv)
ρqᵗ_forcing = (subsidence, ρqᵗ_drying_forcing)
ρθ_forcing = subsidence
ρe_forcing = ρe_radiation_forcing
forcing = (; ρu=ρu_forcing, ρv=ρv_forcing, ρθ=ρθ_forcing,
ρe=ρe_forcing, ρqᵗ=ρqᵗ_forcing)Model setup
We use warm-phase saturation adjustment microphysics and 9th-order WENO advection.
microphysics = SaturationAdjustment(equilibrium=WarmPhaseEquilibrium())
advection = WENO(order=9)
model = AtmosphereModel(grid; dynamics, coriolis, microphysics, advection, forcing,
boundary_conditions = (ρθ=ρθ_bcs, ρqᵗ=ρqᵗ_bcs, ρu=ρu_bcs, ρv=ρv_bcs))AtmosphereModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 64×64×75 RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on CUDAGPU with 5×5×5 halo
├── dynamics: AnelasticDynamics(p₀=101500.0, θ₀=299.1)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── timestepper: RungeKutta3TimeStepper
├── advection scheme:
│ ├── momentum: WENO{5, Float32, Float32}(order=9)
│ ├── ρθ: WENO{5, Float32, Float32}(order=9)
│ └── ρqᵗ: WENO{5, Float32, Float32}(order=9)
├── tracers: ()
├── coriolis: FPlane{Float32}(f=3.76e-5)
└── microphysics: SaturationAdjustmentInitial conditions
Profiles from AtmosphericProfilesLibrary
Mean profiles are specified as piecewise linear functions by Siebesma et al. (2003), Appendix B, Tables B1 and B2, and include:
- Liquid-ice potential temperature $θ^{\ell i}(z)$ (Table B1)
- Total water specific humidity $q^t(z)$ (Table B1)
- Zonal velocity $u(z)$ (Table B2)
The amazing and convenient AtmosphericProfilesLibrary implements functions that retrieve these profiles.
FT = eltype(grid)
θˡⁱ₀ = AtmosphericProfilesLibrary.Bomex_θ_liq_ice(FT)
qᵗ₀ = AtmosphericProfilesLibrary.Bomex_q_tot(FT)
u₀ = AtmosphericProfilesLibrary.Bomex_u(FT)AtmosphericProfilesLibrary.ZProfile{AtmosphericProfilesLibrary.var"#Bomex_u##0#Bomex_u##1"{Float32}}(AtmosphericProfilesLibrary.var"#Bomex_u##0#Bomex_u##1"{Float32}())The initial profiles are perturbed with random noise below 1600 m to trigger convection. The perturbation amplitudes are specified by Siebesma et al. (2003); Appendix B (third paragraph after Eq. B6):
- Potential temperature perturbation: $δθ = 0.1$ K
- Moisture perturbation: $δqᵗ = 2.5 \times 10^{-5}$ kg/kg
Magnitudes for the random perturbations applied to the initial profiles are given by Siebesma et al. (2003), Appendix B, third paragraph after Eq. B6.
δθ = 0.1 # K
δqᵗ = 2.5e-5 # kg/kg
zδ = 1600 # m
ϵ() = rand() - 1/2
θᵢ(x, y, z) = θˡⁱ₀(z) + δθ * ϵ() * (z < zδ)
qᵢ(x, y, z) = qᵗ₀(z) + δqᵗ * ϵ() * (z < zδ)
uᵢ(x, y, z) = u₀(z)
set!(model, θ=θᵢ, qᵗ=qᵢ, u=uᵢ)Simulation
We run the simulation for 6 hours with adaptive time-stepping.
simulation = Simulation(model; Δt=10, stop_time=6hour)
conjure_time_step_wizard!(simulation, cfl=0.7)Output and progress
We add a progress callback and output the hourly time-averages of the horizontally-averaged profiles for post-processing.
θ = liquid_ice_potential_temperature(model)
qˡ = model.microphysical_fields.qˡ
qᵛ = model.microphysical_fields.qᵛ
function progress(sim)
qˡmax = maximum(qˡ)
qᵗmax = maximum(sim.model.specific_moisture)
wmax = maximum(abs, sim.model.velocities.w)
msg = @sprintf("Iter: %d, t: % 12s, Δt: %s, max|w|: %.2e m/s, max(qᵗ): %.2e, max(qˡ): %.2e",
iteration(sim), prettytime(sim), prettytime(sim.Δt), wmax, qᵗmax, qˡmax)
@info msg
return nothing
end
add_callback!(simulation, progress, IterationInterval(1000))
outputs = merge(model.velocities, model.tracers, (; θ, qˡ, qᵛ))
avg_outputs = NamedTuple(name => Average(outputs[name], dims=(1, 2)) for name in keys(outputs))
filename = "bomex.jld2"
simulation.output_writers[:averages] = JLD2Writer(model, avg_outputs; filename,
schedule = AveragedTimeInterval(1hour),
overwrite_existing = true)JLD2Writer scheduled on TimeInterval(1 hour):
├── filepath: bomex.jld2
├── 6 outputs: (u, v, w, θ, qˡ, qᵛ) averaged on AveragedTimeInterval(window=1 hour, stride=1, interval=1 hour)
├── array_type: Array{Float32}
├── including: [:grid]
├── file_splitting: NoFileSplitting
└── file size: 43.1 KiBOutput horizontal slices at z = 600 m for animation Find the k-index closest to z = 600 m
z = Oceananigans.Grids.znodes(grid, Center())
k = searchsortedfirst(z, 800)
@info "Saving slices at z = $(z[k]) m (k = $k)"
u, v, w = model.velocities
slice_fields = (; w, qˡ)
slice_outputs = (
wxy = view(w, :, :, k),
qˡxy = view(qˡ, :, :, k),
wxz = view(w, :, 1, :),
qˡxz = view(qˡ, :, 1, :),
)
simulation.output_writers[:slices] = JLD2Writer(model, slice_outputs;
filename = "bomex_slices.jld2",
schedule = TimeInterval(30seconds),
overwrite_existing = true)
@info "Running BOMEX simulation..."
run!(simulation)[ Info: Saving slices at z = 820.0 m (k = 21)
[ Info: Running BOMEX simulation...
[ Info: Initializing simulation...
[ Info: Iter: 0, t: 0 seconds, Δt: 8.000 seconds, max|w|: 3.75e-07 m/s, max(qᵗ): 1.70e-02, max(qˡ): 0.00e+00
[ Info: ... simulation initialization complete (31.532 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (6.941 seconds).
[ Info: Iter: 1000, t: 1.052 hours, Δt: 3.381 seconds, max|w|: 4.39e+00 m/s, max(qᵗ): 1.88e-02, max(qˡ): 2.15e-03
[ Info: Iter: 2000, t: 1.920 hours, Δt: 3.438 seconds, max|w|: 4.60e+00 m/s, max(qᵗ): 1.90e-02, max(qˡ): 2.73e-03
[ Info: Iter: 3000, t: 2.683 hours, Δt: 2.686 seconds, max|w|: 6.21e+00 m/s, max(qᵗ): 1.89e-02, max(qˡ): 2.86e-03
[ Info: Iter: 4000, t: 3.472 hours, Δt: 2.556 seconds, max|w|: 7.24e+00 m/s, max(qᵗ): 1.96e-02, max(qˡ): 3.05e-03
[ Info: Iter: 5000, t: 4.177 hours, Δt: 2.934 seconds, max|w|: 5.48e+00 m/s, max(qᵗ): 1.97e-02, max(qˡ): 3.10e-03
[ Info: Iter: 6000, t: 4.911 hours, Δt: 2.778 seconds, max|w|: 6.18e+00 m/s, max(qᵗ): 1.95e-02, max(qˡ): 2.35e-03
[ Info: Iter: 7000, t: 5.690 hours, Δt: 2.843 seconds, max|w|: 5.67e+00 m/s, max(qᵗ): 1.96e-02, max(qˡ): 3.37e-03
[ Info: Simulation is stopping after running for 2.252 minutes.
[ Info: Simulation time 6 hours equals or exceeds stop time 6 hours.
Results: mean profile evolution
We visualize the evolution of horizontally-averaged profiles every hour, similar to Figure 3 in the paper by (Siebesma et al., 2003). The intercomparison study shows that after spin-up, the boundary layer reaches a quasi-steady state with:
- A well-mixed layer below cloud base (~500 m)
- A conditionally unstable cloud layer (~500-1500 m)
- A stable inversion layer (~1500-2000 m)
θt = FieldTimeSeries(filename, "θ")
qᵛt = FieldTimeSeries(filename, "qᵛ")
qˡt = FieldTimeSeries(filename, "qˡ")
ut = FieldTimeSeries(filename, "u")
vt = FieldTimeSeries(filename, "v")1×1×75×7 FieldTimeSeries{Oceananigans.OutputReaders.InMemory} located at (⋅, ⋅, Center) of v at bomex.jld2
├── grid: 64×64×75 RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 5×5×5 halo
├── indices: (:, :, :)
├── time_indexing: Linear()
├── backend: InMemory()
├── path: bomex.jld2
├── name: v
└── data: 1×1×85×7 OffsetArray(::Array{Float32, 4}, 1:1, 1:1, -4:80, 1:7) with eltype Float32 with indices 1:1×1:1×-4:80×1:7
└── max=0.0692929, min=-1.38051, mean=-0.0873838Create a 2×2 panel plot showing the evolution of key variables
fig = Figure(size=(900, 800), fontsize=14)
axθ = Axis(fig[1, 1], xlabel="θ (K)", ylabel="z (m)")
axq = Axis(fig[1, 2], xlabel="qᵛ (kg/kg)", ylabel="z (m)")
axuv = Axis(fig[2, 1], xlabel="u, v (m/s)", ylabel="z (m)")
axqˡ = Axis(fig[2, 2], xlabel="qˡ (kg/kg)", ylabel="z (m)")
times = θt.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"
lines!(axθ, θt[n], color=colors[n], label=label)
lines!(axq, qᵛt[n], color=colors[n])
lines!(axuv, ut[n], color=colors[n], linestyle=:solid)
lines!(axuv, vt[n], color=colors[n], linestyle=:dash)
lines!(axqˡ, qˡt[n], color=colors[n])
endSet axis limits to focus on the boundary layer
for ax in (axθ, axq, axuv, axqˡ)
ylims!(ax, 0, 2500)
end
xlims!(axθ, 298, 310)
xlims!(axq, 3e-3, 18e-3)
xlims!(axuv, -10, 2)Add legends and annotations
axislegend(axθ, position=:rb)
text!(axuv, -8.5, 2200, text="solid: u\ndashed: v", fontsize=12)
fig[0, :] = Label(fig, "BOMEX: Mean profile evolution (Siebesma et al., 2003)", fontsize=18, tellwidth=false)
save("bomex_profiles.png", fig)
fig
The simulation shows the development of a cloudy boundary layer with:
- Warming of the subcloud layer from surface fluxes
- Moistening of the lower troposphere
- Development of cloud water in the conditionally unstable layer
- Westerly flow throughout the domain with weak meridional winds
Animation of horizontal slices
We create an animation showing the evolution of vertical velocity and liquid water at z = 800 m, which is near the cloud base level. We limit the animation to the first two hours, where most of the interesting development occurs.
wxz_ts = FieldTimeSeries("bomex_slices.jld2", "wxz")
qˡxz_ts = FieldTimeSeries("bomex_slices.jld2", "qˡxz")
wxy_ts = FieldTimeSeries("bomex_slices.jld2", "wxy")
qˡxy_ts = FieldTimeSeries("bomex_slices.jld2", "qˡxy")
times = wxz_ts.times
Nt = length(times)
x = xnodes(grid, Center())
z = znodes(grid, Center())20.0f0:40.0f0:2980.0f0Create animation
fig = Figure(size=(900, 700), fontsize=14)
axwxz = Axis(fig[2, 2], aspect=2, xaxisposition=:top, xlabel="x (m)", ylabel="z (m)", title="Vertical velocity w")
axqxz = Axis(fig[2, 3], aspect=2, xaxisposition=:top, xlabel="x (m)", ylabel="z (m)", title="Liquid water qˡ")
axwxy = Axis(fig[3, 2], aspect=1, xlabel="x (m)", ylabel="y (m)", title="@ z = $(z[k]) m")
axqxy = Axis(fig[3, 3], aspect=1, xlabel="x (m)", ylabel="y (m)", title="@ z = $(z[k]) m")Axis with 0 plots:
Determine color limits from the data
wlim = maximum(abs, wxz_ts) / 4
qˡlim = maximum(qˡxz_ts) / 4
n = Observable(1)
wxz_n = @lift wxz_ts[$n]
qˡxz_n = @lift qˡxz_ts[$n]
wxy_n = @lift wxy_ts[$n]
qˡxy_n = @lift qˡxy_ts[$n]
title = @lift "BOMEX slices at t = " * prettytime(times[$n])
hmw = heatmap!(axwxz, wxz_n, colormap=:balance, colorrange=(-wlim, wlim))
hmq = heatmap!(axqxz, qˡxz_n, colormap=Reverse(:Blues_4), colorrange=(0, qˡlim))
hmw = heatmap!(axwxy, wxy_n, colormap=:balance, colorrange=(-wlim, wlim))
hmq = heatmap!(axqxy, qˡxy_n, colormap=Reverse(:Blues_4), colorrange=(0, qˡlim))
for ax in (axwxz, axqxz)
lines!(ax, x, fill(z[k], length(x)), color=:grey, linestyle=:dash)
end
Colorbar(fig[2:3, 1], hmw, label="w (m/s)", tellheight=false, height=Relative(0.7), flipaxis=false)
Colorbar(fig[2:3, 4], hmq, label="qˡ (kg/kg)", tellheight=false, height=Relative(0.7))
fig[1, :] = Label(fig, title, fontsize=18, tellwidth=false)
rowgap!(fig.layout, 1, -50)
rowgap!(fig.layout, 2, -50)Record animation
CairoMakie.record(fig, "bomex_slices.mp4", 1:Nt, framerate=12) do nn
n[] = nn
endJulia 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.