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)endDomain and grid
For faster time to solution, we reduce the numerical precision to Float32.
arch = GPU()Oceananigans.defaults.FloatType = Float32Simulation "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.4167Reference state and formulation
p₀ = 1e5 # Paθ₀ = 300 # Kconstants = ThermodynamicConstants()reference_state = ReferenceState(grid, constants, surface_pressure = p₀, potential_temperature = θ₀)dynamics = AnelasticDynamics(reference_state)AnelasticDynamics(p₀=100000.0, θ₀=300.0)
└── pressure_anomaly: not materializedCapping 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 defaultRelaxation{Float64}
├── rate: 0.01
├── mask: exp(-(z - 1000)^2 / (2 * 200^2))
└── target: 0Assembling 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: NothingInitial 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_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: 1.430 minutes; max|w|: 5.08e-03 m/s
[ Info: ... simulation initialization complete (29.629 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (2.858 seconds).
[ Info: Iter: 1000, t: 22.200 minutes, Δt: 1.358 seconds, elapsed: 1.820 minutes; max|w|: 3.62e-01 m/s
[ Info: Iter: 2000, t: 44.649 minutes, Δt: 1.342 seconds, elapsed: 2.008 minutes; max|w|: 3.79e-01 m/s
[ Info: Iter: 3000, t: 1.117 hours, Δt: 1.348 seconds, elapsed: 2.202 minutes; max|w|: 4.64e-01 m/s
[ Info: Iter: 4000, t: 1.482 hours, Δt: 1.283 seconds, elapsed: 2.383 minutes; max|w|: 1.00e+00 m/s
[ Info: Iter: 5000, t: 1.782 hours, Δt: 857.904 ms, elapsed: 2.565 minutes; max|w|: 4.79e+00 m/s
[ Info: Iter: 6000, t: 2.027 hours, Δt: 1.019 seconds, elapsed: 2.746 minutes; max|w|: 2.74e+00 m/s
[ Info: Iter: 7000, t: 2.340 hours, Δt: 1.223 seconds, elapsed: 2.930 minutes; max|w|: 1.88e+00 m/s
[ Info: Iter: 8000, t: 2.685 hours, Δt: 1.246 seconds, elapsed: 3.181 minutes; max|w|: 2.10e+00 m/s
[ Info: Iter: 9000, t: 3.015 hours, Δt: 1.124 seconds, elapsed: 3.445 minutes; max|w|: 2.25e+00 m/s
[ Info: Iter: 10000, t: 3.328 hours, Δt: 1.145 seconds, elapsed: 3.693 minutes; max|w|: 2.35e+00 m/s
[ Info: Iter: 11000, t: 3.649 hours, Δt: 1.101 seconds, elapsed: 3.937 minutes; max|w|: 1.87e+00 m/s
[ Info: Iter: 12000, t: 3.970 hours, Δt: 1.104 seconds, elapsed: 4.181 minutes; max|w|: 2.29e+00 m/s
[ Info: Iter: 13000, t: 4.288 hours, Δt: 1.177 seconds, elapsed: 4.424 minutes; max|w|: 1.92e+00 m/s
[ Info: Iter: 14000, t: 4.604 hours, Δt: 1.145 seconds, elapsed: 4.667 minutes; max|w|: 2.07e+00 m/s
[ Info: Iter: 15000, t: 4.915 hours, Δt: 1.110 seconds, elapsed: 4.926 minutes; max|w|: 2.26e+00 m/s
[ Info: Simulation is stopping after running for 3.901 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.11493Compute 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.0and 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.PrendDefine 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.12.6
Commit 15346901f00 (2026-04-09 19:20 UTC)
Build Info:
Official https://julialang.org release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 8 × AMD EPYC 7R13 Processor
WORD_SIZE: 64
LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
JULIA_GPG = 3673DF529D9049477F76B37566E3C7DC03D6E495
JULIA_LOAD_PATH = :@breeze
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
JULIA_VERSION = 1.12.6
JULIA_DEPOT_PATH = /usr/local/share/julia:
JULIA_PATH = /usr/local/julia
JULIA_PROJECT = @breeze
These were the top-level packages installed in the environment:
import PkgPkg.status()Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
⌃ [86bc3604] AtmosphericProfilesLibrary v0.1.8
[660aa2fb] Breeze v0.8.0 `.`
[052768ef] CUDA v6.2.1
[13f3f980] CairoMakie v0.15.13
⌅ [6a9e3e04] CloudMicrophysics v0.37.1
[e30172f5] Documenter v1.17.0
[daee34ce] DocumenterCitations v1.4.1
[b6400b83] DocumenterCodeBlocks v1.1.0
[7da242da] Enzyme v0.13.198
[46192b85] GPUArraysCore v0.2.0
[63c18a36] KernelAbstractions v0.9.42
[98b081ad] Literate v2.21.0
[85f8d34a] NCDatasets v0.14.15
[9e8cae18] Oceananigans v0.110.14
⌅ [a01a1ee8] RRTMGP v0.21.9
[3c362404] Reactant v0.2.278
[b77e0a4c] InteractiveUtils v1.11.0
[44cfe95a] Pkg v1.12.1
[9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`
This page was generated using Literate.jl.