Neutral atmospheric boundary layer (ABL)

This canonical setup is based on the paper by Moeng and Sullivan (1994), which was a demonstration case for the NCAR LES subgrid-scale model development (Sullivan et al., 1994). Sometimes, this model configuration is called "conventionally" neutral (Pedersen et al., 2014) or "conditionally" neutral (Berg et al., 2020), which indicate an idealized dry, shear-driven atmospheric boundary layer, capped by a stable inversion layer, without any surface heating. Forcings come from a specified constant geostrophic wind (i.e., a specified background pressure gradient) and Coriolis forces; the temperature lapse rate in the free atmosphere is maintained with a sponge layer.

In lieu of more sophisticated surface layer modeling in this example, we impose a fixed friction velocity at the bottom boundary.

using Breeze
using Oceananigans: Oceananigans
using Oceananigans.Units
using CUDA
using Printf
using Random
using CairoMakie

Random.seed!(1994)
if CUDA.functional()
    CUDA.seed!(1994)
end

Domain and grid

For faster time to solution, we reduce the numerical precision to Float32.

arch = GPU()
Oceananigans.defaults.FloatType = Float32;

Simulation "S" (shear-driven ABL) domain setup from Moeng and Sullivan (1994):

Nx = Ny = Nz = 96
x = y = (0, 3000)
z = (0, 1000)

grid = RectilinearGrid(arch; x, y, z,
                       size = (Nx, Ny, Nz), halo = (5, 5, 5),
                       topology = (Periodic, Periodic, Bounded))

96×96×96 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── Periodic x ∈ [0.0, 3000.0) regularly spaced with Δx=31.25
├── Periodic y ∈ [0.0, 3000.0) regularly spaced with Δy=31.25
└── Bounded  z ∈ [0.0, 1000.0] regularly spaced with Δz=10.4167

Reference state and formulation

p₀ = 1e5   # Pa
θ₀ = 300   # K

constants = ThermodynamicConstants()

reference_state = ReferenceState(grid, constants,
                                 surface_pressure = p₀,
                                 potential_temperature = θ₀)

dynamics = AnelasticDynamics(reference_state)

AnelasticDynamics(p₀=100000.0, θ₀=300.0)
└── pressure_anomaly: not materialized

Capping inversion for "S" simulation, as in the paper by Moeng and Sullivan (1994): The base of the inversion is at 468 m and has a thickness of 6 grid levels, over which the potential temperature increases by 8 K. Above the cap, the lapse rate is 3 K/km.

Δz  = first(zspacings(grid))
zᵢ₁ = 468        # m
zᵢ₂ = zᵢ₁ + 6Δz  # m
Γᵢ  = 8 / 6Δz    # K/m
Γᵗᵒᵖ = 0.003     # K/m
θᵣ(z) = z < zᵢ₁ ? θ₀ :
        z < zᵢ₂ ? θ₀ + Γᵢ * (z - zᵢ₁) :
        θ₀ + Γᵢ * (zᵢ₂ - zᵢ₁) + Γᵗᵒᵖ * (z - zᵢ₂)

Surface momentum flux (drag)

For testing, we prescribe the surface shear stress. In practice, however, this is not known a priori. A surface layer scheme (i.e., a wall model) will dynamically update $u_★$ based on environmental conditions that include surface roughness and heat fluxes.

u★ = 0.5  # m/s, _result_ from simulation "S" by Moeng and Sullivan (1994)
q₀ = Breeze.Thermodynamics.MoistureMassFractions{eltype(grid)} |> zero
ρ₀ = Breeze.Thermodynamics.density(θ₀, p₀, q₀, constants)

A bulk drag parameterization is applied with friction velocity:

@inline ρu_drag(x, y, t, ρu, ρv, param) = - param.ρ₀ * param.u★^2 * ρu / max(sqrt(ρu^2 + ρv^2), 1e-6)
@inline ρv_drag(x, y, t, ρu, ρv, param) = - param.ρ₀ * param.u★^2 * ρv / max(sqrt(ρu^2 + ρv^2), 1e-6)

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

Sponge layer

To enforce an upper-air temperature gradient, we introduce a sponge layer with Gaussian weighting that corresponds to an effective depth of approximately 500 m. At |z - zᵗᵒᵖ| = 500, exp(-0.5 * (500/sponge_width)^2) = 0.04 ~ 0. The sponge_rate (inverse timescale) is an ad hoc value; a higher sponge rate (shorter damping time scale) made no difference in this case, and a weaker sponge rate may be used.

sponge_width = 200  # m
sponge_rate = 0.01  # 1/s
sponge_mask = GaussianMask{:z}(center = last(z), width = sponge_width)

Oceananigans.Forcings.GaussianMask{:z, Int64}(1000, 200)

We relax potential temperature to the initial profile using a discrete forcing.

ρθᵣ = Field{Nothing, Nothing, Center}(grid)
set!(ρθᵣ, z -> θᵣ(z))
set!(ρθᵣ, reference_state.density * ρθᵣ)

ρθᵣ_data = interior(ρθᵣ, 1, 1, :)

@inline function ρθ_sponge_fun(i, j, k, grid, clock, model_fields, p)
    zₖ = znode(k, grid, Center())
    return @inbounds p.rate * p.mask(0, 0, zₖ) * (p.target[k] - model_fields.ρθ[i, j, k])
end

ρθ_sponge = Forcing(
    ρθ_sponge_fun; discrete_form = true,
    parameters = (rate = sponge_rate, mask = sponge_mask, target = ρθᵣ_data)
)

DiscreteForcing{@NamedTuple{rate::Float64, mask::Oceananigans.Forcings.GaussianMask{:z, Int64}, target::SubArray{Float32, 1, CUDA.CuArray{Float32, 3, CUDA.DeviceMemory}, Tuple{Int64, Int64, Base.Slice{Base.OneTo{Int64}}}, true}}}
├── func: ρθ_sponge_fun (generic function with 1 method)
└── parameters: (rate = 0.01, mask = Oceananigans.Forcings.GaussianMask{:z, Int64}(1000, 200), target = Float32[348.28134, 347.9861, 347.69092, 347.39597, 347.10114, 346.8064, 346.51187, 346.21753, 345.92325, 345.62918, 345.33527, 345.04144, 344.74783, 344.45435, 344.16104, 343.86783, 343.5748, 343.28195, 342.98917, 342.69662, 342.40417, 342.11185, 341.81976, 341.52777, 341.23593, 340.94424, 340.65274, 340.36136, 340.0701, 339.77902, 339.4881, 339.19727, 338.90665, 338.61618, 338.32584, 338.03564, 337.7456, 337.45572, 337.166, 336.87637, 336.58698, 336.29767, 336.0085, 335.71954, 335.43066, 335.994, 337.19296, 338.38947, 339.58365, 340.7754, 341.96463, 342.32385, 342.06323, 341.80264, 341.54218, 341.2819, 341.02158, 340.7614, 340.5013, 340.24136, 339.98145, 339.72165, 339.46194, 339.20233, 338.9428, 338.6834, 338.42407, 338.1648, 337.90567, 337.64664, 337.38763, 337.1288, 336.87003, 336.61133, 336.35278, 336.09424, 335.83585, 335.57758, 335.31943, 335.06122, 334.8032, 334.5453, 334.28748, 334.0297, 333.77203, 333.51453, 333.25705, 332.9997, 332.74243, 332.48523, 332.22815, 331.97116, 331.71423, 331.45746, 331.20074, 330.94415])

We also damp out any vertical motions near the top boundary.

ρw_sponge = Relaxation(rate = sponge_rate, mask = sponge_mask) # relaxes to 0 by default

Relaxation{Float64, Oceananigans.Forcings.GaussianMask{:z, Int64}, typeof(Oceananigans.Forcings.zerofunction)}
├── rate: 0.01
├── mask: exp(-(z - 1000)^2 / (2 * 200^2))
└── target: 0

Assembling all the forcings

coriolis = FPlane(f=1e-4)

uᵍ, vᵍ = 15, 0  # m/s, simulation "S" by Moeng and Sullivan (1994)
geostrophic = geostrophic_forcings(uᵍ, vᵍ)

ρu_forcing = geostrophic.ρu
ρv_forcing = geostrophic.ρv
ρw_forcing = ρw_sponge
ρθ_forcing = ρθ_sponge

forcing = (; ρu=ρu_forcing, ρv=ρv_forcing, ρw=ρw_forcing, ρθ=ρθ_forcing)

Model setup

advection = WENO(order=9)   # WENO(order=5), Centered(order=6) are too dissipative

closure = SmagorinskyLilly()

model = AtmosphereModel(grid; dynamics, coriolis, advection, forcing, closure,
                        boundary_conditions = (ρu=ρu_bcs, ρv=ρv_bcs))

AtmosphereModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 96×96×96 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CUDAGPU with 5×5×5 halo
├── dynamics: AnelasticDynamics(p₀=100000.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float32}
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│   ├── ρθ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│   └── ρqᵛ: WENO{5, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
├── forcing: ρu=>GeostrophicForcing, ρv=>GeostrophicForcing, ρw=>ContinuousForcing, ρθ=>DiscreteForcing
├── tracers: ()
├── coriolis: FPlane{Float32}(f=0.0001)
└── microphysics: Nothing

Initial conditions

The velocity field is initialized to the constant geostrophic wind; the potential temperature field is initialized to the profile defined in the reference state section above. We add velocity and temperature perturbations to help initiate turbulence.

δu = δv = 0.01  # m/s
δθ = 0.1        # K
zδ = 400        # m, < zᵢ₁

ϵ() = rand() - 1/2
uᵢ(x, y, z) =   uᵍ  + δu * ϵ() * (z < zδ)
vᵢ(x, y, z) =   vᵍ  + δv * ϵ() * (z < zδ)
θᵢ(x, y, z) = θᵣ(z) + δθ * ϵ() * (z < zδ)

set!(model, θ=θᵢ, u=uᵢ, v=vᵢ)

Simulation and output

We will run the simulation for 5 hours with adaptive time-stepping.

simulation = Simulation(model; Δt=0.5, stop_time=5hours)
conjure_time_step_wizard!(simulation, cfl=0.7)

Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)

Progress monitor

A progress callback is added to monitor the simulation.

u, v, w = model.velocities
θ = liquid_ice_potential_temperature(model)
νₑ = model.closure_fields.νₑ

# For keeping track of the computational expense
wall_clock = time_ns()

function progress(sim)
    wmax = maximum(abs, sim.model.velocities.w)
    elapsed = 1e-9 * (time_ns() - wall_clock)
    msg = @sprintf("Iter: %d, t: % 12s, Δt: %s, elapsed: %s; max|w|: %.2e m/s",
                   iteration(sim), prettytime(sim), prettytime(sim.Δt), prettytime(elapsed), wmax)
    @info msg
    return nothing
end

add_callback!(simulation, progress, IterationInterval(1000))

Horizontal averaging

Profiles of horizontally averaged quantities are output every 10 minutes for statistical analysis. All outputs are at cell centers.

Note: Higher-order moments are computed at the location of the first field. E.g., u * w results in a BinaryOperation at (Face, Center, Center).

avg_outputs_varlist = (;
    θ, νₑ,
    uu = u^2, vv = v^2, ww = w^2,
    uw = u*w, vw = v*w, θw = θ*w,           # second-order moments for fluxes
    uuw = u^2*w, vvw = v^2*w, www = w^3,    # third-order moments to calculate turbulent transport
    νₑ³ = νₑ^3,                             # SGS dissipation — note: |S̄|² = νₑ² / (Cₛ Δ)⁴) with Smagorinsky model
)
outputs = merge(model.velocities, model.tracers, avg_outputs_varlist)

After computing the output quantity and prior to calculating the slab average, staggered quantities are interpolated to cell centers.

avg_outputs = NamedTuple(name => Average(@at((Center, Center, Center), outputs[name]), dims=(1, 2))
                         for name in keys(outputs))

# Calculate derivatives using an AbstractOperation, `∂z`, to facilitate postprocessing later.
∂z_outputs = (; ∂z_u=u, ∂z_v=v, ∂z_θ=θ)
avg_∂z_outputs = NamedTuple(name => Average(@at((Center, Center, Center), ∂z(∂z_outputs[name])), dims=(1, 2))
                            for name in keys(∂z_outputs))

Set up the output writer.

avg_filename = "abl_averages.jld2"
avg_output_interval = 10minutes
simulation.output_writers[:averages] = JLD2Writer(model, merge(avg_outputs, avg_∂z_outputs);
                                                  filename = avg_filename,
                                                  schedule = AveragedTimeInterval(avg_output_interval),
                                                  overwrite_existing = true)

JLD2Writer scheduled on TimeInterval(10 minutes):
├── filepath: abl_averages.jld2
├── 18 outputs: (u, v, w, θ, νₑ, uu, vv, ww, uw, vw, θw, uuw, vvw, www, νₑ³, ∂z_u, ∂z_v, ∂z_θ) averaged on AveragedTimeInterval(window=10 minutes, stride=1, interval=10 minutes)
├── array_type: Array{Float32}
├── including: [:thermodynamic_constants]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

Instantaneous slices for animation

Horizontal (xy) slices: Find the k-index closest to z = 100 m.

z = znodes(grid, Center())
k₁₀₀ = searchsortedfirst(z, 100)
@info "Saving slices at z = $(z[k₁₀₀]) m (k = $k₁₀₀)"

[ Info: Saving slices at z = 109.375 m (k = 11)

Vertical (xz) slices: Find the j-index closest to the domain center.

y = ynodes(grid, Center())
jmid = Ny ÷ 2
@info "Saving slices at y = $(y[jmid]) m (j = $jmid)"

[ Info: Saving slices at y = 1484.375 m (j = 48)

Set up another output writer.

slice_fields = (; u, v, w, θ)
slice_outputs = (
    u_xy = view(u, :, :, k₁₀₀),
    v_xy = view(v, :, :, k₁₀₀),
    w_xy = view(w, :, :, k₁₀₀),
    u_xz = view(u, :, jmid, :),
    w_xz = view(w, :, jmid, :),
    θ_xz = view(θ, :, jmid, :),
)

simulation.output_writers[:slices] = JLD2Writer(model, slice_outputs;
                                                filename = "abl_slices.jld2",
                                                schedule = TimeInterval(5minutes),
                                                overwrite_existing = true)

JLD2Writer scheduled on TimeInterval(5 minutes):
├── filepath: abl_slices.jld2
├── 6 outputs: (u_xy, v_xy, w_xy, u_xz, w_xz, θ_xz)
├── array_type: Array{Float32}
├── including: [:thermodynamic_constants]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

Go time

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter: 0, t:    0 seconds, Δt: 550.000 ms, elapsed: 2.739 minutes; max|w|: 4.56e-03 m/s
[ Info:     ... simulation initialization complete (45.364 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (4.996 seconds).
[ Info: Iter: 1000, t: 22.254 minutes, Δt: 1.364 seconds, elapsed: 4.456 minutes; max|w|: 3.26e-01 m/s
[ Info: Iter: 2000, t: 44.744 minutes, Δt: 1.346 seconds, elapsed: 5.812 minutes; max|w|: 3.76e-01 m/s
[ Info: Iter: 3000, t:  1.118 hours, Δt: 1.330 seconds, elapsed: 7.288 minutes; max|w|: 5.32e-01 m/s
[ Info: Iter: 4000, t:  1.470 hours, Δt: 1.100 seconds, elapsed: 8.751 minutes; max|w|: 1.57e+00 m/s
[ Info: Iter: 5000, t:  1.729 hours, Δt: 826.240 ms, elapsed: 10.945 minutes; max|w|: 3.97e+00 m/s
[ Info: Iter: 6000, t:  1.992 hours, Δt: 1.069 seconds, elapsed: 13.717 minutes; max|w|: 2.43e+00 m/s
[ Info: Iter: 7000, t:  2.320 hours, Δt: 1.277 seconds, elapsed: 16.436 minutes; max|w|: 1.70e+00 m/s
[ Info: Iter: 8000, t:  2.662 hours, Δt: 1.216 seconds, elapsed: 18.637 minutes; max|w|: 1.77e+00 m/s
[ Info: Iter: 9000, t:  2.988 hours, Δt: 1.198 seconds, elapsed: 20.855 minutes; max|w|: 1.92e+00 m/s
[ Info: Iter: 10000, t:  3.312 hours, Δt: 1.160 seconds, elapsed: 23.574 minutes; max|w|: 2.00e+00 m/s
[ Info: Iter: 11000, t:  3.635 hours, Δt: 1.180 seconds, elapsed: 25.802 minutes; max|w|: 1.81e+00 m/s
[ Info: Iter: 12000, t:  3.949 hours, Δt: 1.072 seconds, elapsed: 28.005 minutes; max|w|: 2.19e+00 m/s
[ Info: Iter: 13000, t:  4.257 hours, Δt: 1.154 seconds, elapsed: 30.372 minutes; max|w|: 2.06e+00 m/s
[ Info: Iter: 14000, t:  4.571 hours, Δt: 1.126 seconds, elapsed: 32.805 minutes; max|w|: 2.26e+00 m/s
[ Info: Iter: 15000, t:  4.878 hours, Δt: 1.069 seconds, elapsed: 34.789 minutes; max|w|: 2.10e+00 m/s
[ Info: Simulation is stopping after running for 33.308 minutes.
[ Info: Simulation time 5 hours equals or exceeds stop time 5 hours.

Load output and visualize

Let's load the saved output.

u_ts  = FieldTimeSeries(avg_filename, "u")
v_ts  = FieldTimeSeries(avg_filename, "v")
w_ts  = FieldTimeSeries(avg_filename, "w")
θ_ts  = FieldTimeSeries(avg_filename, "θ")
uu_ts = FieldTimeSeries(avg_filename, "uu")
vv_ts = FieldTimeSeries(avg_filename, "vv")
ww_ts = FieldTimeSeries(avg_filename, "ww")
uw_ts = FieldTimeSeries(avg_filename, "uw")
vw_ts = FieldTimeSeries(avg_filename, "vw")
θw_ts = FieldTimeSeries(avg_filename, "θw")
νₑ_ts = FieldTimeSeries(avg_filename, "νₑ")
∂z_u_ts = FieldTimeSeries(avg_filename, "∂z_u")
∂z_v_ts = FieldTimeSeries(avg_filename, "∂z_v")
∂z_θ_ts = FieldTimeSeries(avg_filename, "∂z_θ")

grid = u_ts.grid
times = u_ts.times
Nt = length(times)

ρᵣ = Oceananigans.on_architecture(CPU(), reference_state.density)

1×1×96 Field{Nothing, Nothing, Oceananigans.Grids.Center} reduced over dims = (1, 2) on Oceananigans.Grids.RectilinearGrid on CPU
├── grid: 96×96×96 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CPU with 5×5×5 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: Value, top: ZeroFlux, immersed: Nothing
└── data: 1×1×106 OffsetArray(::Array{Float32, 3}, 1:1, 1:1, -4:101) with eltype Float32 with indices 1:1×1:1×-4:101
    └── max=1.16094, min=1.06966, mean=1.11493

Compute diagnostics at each saved time. First, we create a few empty timeseries to save the computed diagnostics.

loc = (nothing, nothing, Center())

WS_mean_ts = FieldTimeSeries(loc, grid, times)
WD_mean_ts = FieldTimeSeries(loc, grid, times)
 uu_var_ts = FieldTimeSeries(loc, grid, times)
 vv_var_ts = FieldTimeSeries(loc, grid, times)
 ww_var_ts = FieldTimeSeries(loc, grid, times)
 uw_res_ts = FieldTimeSeries(loc, grid, times)
 vw_res_ts = FieldTimeSeries(loc, grid, times)
 θw_res_ts = FieldTimeSeries(loc, grid, times)
 uw_sgs_ts = FieldTimeSeries(loc, grid, times)
 vw_sgs_ts = FieldTimeSeries(loc, grid, times)
 θw_sgs_ts = FieldTimeSeries(loc, grid, times)

1×1×96×31 FieldTimeSeries{Oceananigans.OutputReaders.InMemory} located at (⋅, ⋅, Center) on Oceananigans.Architectures.CPU
├── grid: 96×96×96 RectilinearGrid{Float32, Periodic, Periodic, Bounded} on CPU with 5×5×5 halo
├── indices: (:, :, :)
├── time_indexing: Clamp()
├── backend: InMemory()
└── data: 1×1×106×31 OffsetArray(::Array{Float32, 4}, 1:1, 1:1, -4:101, 1:31) with eltype Float32 with indices 1:1×1:1×-4:101×1:31
    └── max=0.0, min=0.0, mean=0.0

and then we loop over all saved fields and compute what we want.

for n in 1:Nt
    u_n = u_ts[n]
    v_n = v_ts[n]
    w_n = w_ts[n]
    θ_n = θ_ts[n]
    νₑ_n = νₑ_ts[n]
    uu_n = uu_ts[n]
    vv_n = vv_ts[n]
    ww_n = ww_ts[n]
    uw_n = uw_ts[n]
    vw_n = vw_ts[n]
    θw_n = θw_ts[n]
    ∂z_u_n = ∂z_u_ts[n]
    ∂z_v_n = ∂z_v_ts[n]
    ∂z_θ_n = ∂z_θ_ts[n]

    WS_mean_ts[n] .= sqrt(u_n^2 + v_n^2)
    WD_mean_ts[n] .= mod(270 - atand(v_n, u_n), 360)

    uu_var_ts[n] .= (uu_n - u_n^2) / u★^2
    vv_var_ts[n] .= (vv_n - v_n^2) / u★^2
    ww_var_ts[n] .= (ww_n - w_n^2) / u★^2
    uw_res_ts[n] .= ρᵣ * (uw_n - u_n * w_n) / (ρ₀ * u★^2)
    vw_res_ts[n] .= ρᵣ * (vw_n - v_n * w_n) / (ρ₀ * u★^2)
    θw_res_ts[n] .= θw_n - θ_n * w_n

    uw_sgs_ts[n] .= -ρᵣ * νₑ_n * ∂z_u_n / (ρ₀ * u★^2)
    vw_sgs_ts[n] .= -ρᵣ * νₑ_n * ∂z_v_n / (ρ₀ * u★^2)
    θw_sgs_ts[n] .= -νₑ_n * ∂z_θ_n / closure.Pr  # not normalized
end

Define a colormap for each time.

cmap = cgrad(:viridis)
colors = [cmap[(n-1)/max(Nt-1, 1)] for n in 1:Nt]

Note that the AveragedTimeInterval schedule outputs the initial snapshot followed by averaged snapshots at the requested output interval.

smart_label(n) = prettytime(n * avg_output_interval)
labels = [n == 1 ? "initial condition" : smart_label(n-1) for n in 1:Nt]

We are now ready to plot.

plot_interval = 1hour  # should be a multiple of avg_output_interval
plot_skip = Int(plot_interval / avg_output_interval)

First, we plot the mean profiles (wind speed, wind direction, potential temperature).

fig1 = Figure(size=(1000, 500), fontsize=14)

ax1a = Axis(fig1[1, 1], xlabel="√(U² + V²) (m/s)", ylabel="z (m)", title="Horizontal wind speed")
ax1b = Axis(fig1[1, 2], xlabel="WD (° from N)", ylabel="z (m)", title="Wind direction")
ax1c = Axis(fig1[1, 3], xlabel="θ (K)", ylabel="z (m)", title="Potential temperature")

for n in 1:plot_skip:Nt
    lines!(ax1a, WS_mean_ts[n], color=colors[n], label=labels[n])
    lines!(ax1b, WD_mean_ts[n], color=colors[n])
    lines!(ax1c,       θ_ts[n], color=colors[n])
end

linkyaxes!(ax1a, ax1b, ax1c)
hideydecorations!(ax1b, grid=false)
hideydecorations!(ax1c, grid=false)

Legend(fig1[1, 4], ax1a, "Time", framevisible=false)
fig1

Next, we plot the velocity variances normalized by $u_★^2$

fig2 = Figure(size=(1000, 500), fontsize=14)

ax2a = Axis(fig2[1, 1], xlabel="⟨u′u′⟩ / u_★²", ylabel="z (m)", title="u variance")
ax2b = Axis(fig2[1, 2], xlabel="⟨v′v′⟩ / u_★²", ylabel="z (m)", title="v variance")
ax2c = Axis(fig2[1, 3], xlabel="⟨w′w′⟩ / u_★²", ylabel="z (m)", title="w variance")

for n in 1:plot_skip:Nt
    lines!(ax2a, uu_var_ts[n], color=colors[n], label=labels[n])
    lines!(ax2b, vv_var_ts[n], color=colors[n])
    lines!(ax2c, ww_var_ts[n], color=colors[n])
end

linkyaxes!(ax2a, ax2b, ax2c)
hideydecorations!(ax2b, grid=false)
hideydecorations!(ax2c, grid=false)

Legend(fig2[1, 4], ax2a, "Time", framevisible=false)
fig2

Last, we plot the resolved and the SGS fluxes.

fig3 = Figure(size=(1000, 500), fontsize=14)

ax3a = Axis(fig3[1, 1], xlabel="τˣ / ρ₀u_★²", ylabel="z (m)", title="x-momentum flux")
ax3b = Axis(fig3[1, 2], xlabel="τʸ / ρ₀u_★²", ylabel="z (m)", title="y-momentum flux")
ax3c = Axis(fig3[1, 3], xlabel="Jᶿ (K m/s)", ylabel="z (m)", title="Potential temperature flux")

for n in 1:plot_skip:Nt
    lines!(ax3a, uw_res_ts[n] + uw_sgs_ts[n], color=colors[n], label=labels[n])
    lines!(ax3a, uw_res_ts[n], color=colors[n], linestyle=:dash)
    lines!(ax3a, uw_sgs_ts[n], color=colors[n], linestyle=:dot)

    lines!(ax3b, vw_res_ts[n] + vw_sgs_ts[n], color=colors[n])
    lines!(ax3b, vw_res_ts[n], color=colors[n], linestyle=:dash)
    lines!(ax3b, vw_sgs_ts[n], color=colors[n], linestyle=:dot)

    lines!(ax3c, θw_res_ts[n] + θw_sgs_ts[n], color=colors[n])
    lines!(ax3c, θw_res_ts[n], color=colors[n], linestyle=:dash)
    lines!(ax3c, θw_sgs_ts[n], color=colors[n], linestyle=:dot)
end

for ax in (ax3a, ax3b, ax3c)
    vlines!(ax, 0, color=:grey, linewidth=0.5)
end

# Legends: line style (inside panel a) and time (right)
style_entries = [LineElement(color=:black, linestyle=:solid),
                 LineElement(color=:black, linestyle=:dash),
                 LineElement(color=:black, linestyle=:dot)]
axislegend(ax3a, style_entries, ["total", "resolved", "SGS"],
           position=:lt, framevisible=false)

linkyaxes!(ax3a, ax3b, ax3c)
hideydecorations!(ax3b, grid=false)
hideydecorations!(ax3c, grid=false)

Legend(fig3[1, 4], ax3a, "Time", framevisible=false)
fig3

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

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 Pkg
Pkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
  [86bc3604] AtmosphericProfilesLibrary v0.1.7
  [660aa2fb] Breeze v0.4.7 `.`
⌃ [052768ef] CUDA v5.11.0
  [13f3f980] CairoMakie v0.15.9
⌅ [6a9e3e04] CloudMicrophysics v0.32.0
  [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.4.1
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.107.3
  [a01a1ee8] RRTMGP v0.21.7
  [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.