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 Breezeusing Oceananigans: Oceananigansusing Oceananigans.Unitsusing CUDAusing Printfusing Randomusing CairoMakieRandom.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 = 96x = 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   # Kconstants = ThermodynamicConstants()reference_state = ReferenceState(grid, constants,                                 base_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        # mzᵢ₂ = 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  # msponge_rate = 0.01  # 1/ssponge_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, CUDACore.CuArray{Float32, 3, CUDACore.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}
├──   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ᵍ)forcing = (; u=geostrophic.u, v=geostrophic.v, ρw=ρw_sponge, ρθ=ρθ_sponge)

Model setup

advection = WENO(order=9)   # WENO(order=5), Centered(order=6) are too dissipativeclosure = 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=>SpecificForcing, ρv=>SpecificForcing, ρw=>Relaxation, ρθ=>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        # Kzδ = 400        # m, < zᵢ₁ϵ() = rand() - 1/2uᵢ(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 expensewall_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 nothingendadd_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 = 10minutessimulation.output_writers[:averages] = JLD2Writer(model, merge(avg_outputs, avg_∂z_outputs);                                                  filename = avg_filename,                                                  schedule = AveragedTimeInterval(avg_output_interval),                                                  overwrite_files = 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_files = 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: 1.425 minutes; max|w|: 4.92e-03 m/s
[ Info:     ... simulation initialization complete (44.086 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.445 seconds).
[ Info: Iter: 1000, t: 22.234 minutes, Δt: 1.366 seconds, elapsed: 2.104 minutes; max|w|: 3.30e-01 m/s
[ Info: Iter: 2000, t: 44.718 minutes, Δt: 1.344 seconds, elapsed: 2.613 minutes; max|w|: 4.04e-01 m/s
[ Info: Iter: 3000, t:  1.121 hours, Δt: 1.356 seconds, elapsed: 3.126 minutes; max|w|: 4.59e-01 m/s
[ Info: Iter: 4000, t:  1.463 hours, Δt: 1.071 seconds, elapsed: 3.643 minutes; max|w|: 2.28e+00 m/s
[ Info: Iter: 5000, t:  1.719 hours, Δt: 852.233 ms, elapsed: 4.153 minutes; max|w|: 3.99e+00 m/s
[ Info: Iter: 6000, t:  1.977 hours, Δt: 962.837 ms, elapsed: 4.830 minutes; max|w|: 3.05e+00 m/s
[ Info: Iter: 7000, t:  2.292 hours, Δt: 1.204 seconds, elapsed: 5.371 minutes; max|w|: 2.05e+00 m/s
[ Info: Iter: 8000, t:  2.624 hours, Δt: 1.194 seconds, elapsed: 5.883 minutes; max|w|: 2.18e+00 m/s
[ Info: Iter: 9000, t:  2.948 hours, Δt: 1.142 seconds, elapsed: 6.390 minutes; max|w|: 2.09e+00 m/s
[ Info: Iter: 10000, t:  3.259 hours, Δt: 1.149 seconds, elapsed: 6.893 minutes; max|w|: 1.92e+00 m/s
[ Info: Iter: 11000, t:  3.570 hours, Δt: 1.185 seconds, elapsed: 7.403 minutes; max|w|: 2.43e+00 m/s
[ Info: Iter: 12000, t:  3.886 hours, Δt: 1.133 seconds, elapsed: 7.898 minutes; max|w|: 1.95e+00 m/s
[ Info: Iter: 13000, t:  4.195 hours, Δt: 1.064 seconds, elapsed: 8.397 minutes; max|w|: 2.87e+00 m/s
[ Info: Iter: 14000, t:  4.505 hours, Δt: 1.113 seconds, elapsed: 8.896 minutes; max|w|: 1.99e+00 m/s
[ Info: Iter: 15000, t:  4.812 hours, Δt: 1.071 seconds, elapsed: 9.396 minutes; max|w|: 2.27e+00 m/s
[ Info: Simulation is stopping after running for 8.912 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.gridtimes = u_ts.timesNt = 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())wind_speed_ts = FieldTimeSeries(loc, grid, times)wind_direction_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    # Load saved output    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]    # Compute other quantities    wind_speed_ts[n] .= sqrt(u_n^2 + v_n^2)    wind_direction_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.Prend

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  # a multiple of avg_output_intervalplot_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, wind_speed_ts[n], color=colors[n], label=labels[n])    lines!(ax1b, wind_direction_ts[n], color=colors[n])    lines!(ax1c,       θ_ts[n], color=colors[n])endlinkyaxes!(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])endlinkyaxes!(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)endfor 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: versioninfoversioninfo()

Julia Version 1.13.0
Commit d1c37793dd2 (2026-09-09 19:00 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-20.1.8 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 5 GC (on 8 virtual cores)
Environment:
  JULIA_NUM_GC_THREADS = 4,1
  JULIA_GPG = 64B779A570972FFF7BFC2B54EAD471E1A1F2C10A
  JULIA_LOAD_PATH = :@breeze
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_VERSION = 1.13.0
  JULIA_DEPOT_PATH = /usr/local/share/julia:
  JULIA_PATH = /usr/local/julia
  JULIA_PROJECT = @breeze

These were the top-level packages installed in the environment:

import PkgPkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
  [86bc3604] AtmosphericProfilesLibrary v0.1.10
  [660aa2fb] Breeze v0.11.3 `.`
⌃ [052768ef] CUDA v6.2.2
  [13f3f980] CairoMakie v0.15.15
  [6a9e3e04] CloudMicrophysics v0.41.0
  [e30172f5] Documenter v1.19.0
  [daee34ce] DocumenterCitations v1.5.0
  [b6400b83] DocumenterCodeBlocks v1.5.1
  [1c52b33b] DocumenterLandingPage v0.2.2
  [7da242da] Enzyme v0.13.206
⌅ [46192b85] GPUArraysCore v0.2.0
  [63c18a36] KernelAbstractions v0.9.43
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.113.4
  [a01a1ee8] RRTMGP v1.0.1
  [3c362404] Reactant v0.2.289
  [276daf66] SpecialFunctions v2.9.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.13.0
  [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.