Single column radiation

This example sets up a single-column atmospheric model with an idealized temperature and moisture profile. We compute radiative fluxes using RRTMGP's gray atmosphere solver with the optical thickness parameterization by O'Gorman and Schneider (2008).

using Breeze
using Oceananigans.Units
using CairoMakie

using RRTMGP.AtmosphericStates: GrayOpticalThicknessOGorman2008

Grid and thermodynamics

We create a single column spanning 20 km with 64 layers at a particular place.

Nz = 64
λ, φ = -76.13, 39.48

grid = RectilinearGrid(size=Nz, x=λ, y=φ, z=(0, 20kilometers),
                       topology=(Flat, Flat, Bounded))

1×1×64 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 0×0×3 halo
├── Flat x = -76.13             
├── Flat y = 39.48              
└── Bounded  z ∈ [0.0, 20000.0] regularly spaced with Δz=312.5

Set up the thermodynamic constants and reference state.

surface_temperature = 300
constants = ThermodynamicConstants()

reference_state = ReferenceState(grid, constants;
                                 surface_pressure = 101325,
                                 potential_temperature = surface_temperature)

dynamics = AnelasticDynamics(reference_state)

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

Radiative transfer model

We create a gray radiative transfer model using the O'Gorman and Schneider (2008) optical thickness parameterization. The solar zenith angle is computed from the model clock and grid location.

using Dates

optical_thickness = GrayOpticalThicknessOGorman2008(eltype(grid))
radiation = RadiativeTransferModel(grid, constants, optical_thickness;
                                   surface_temperature,
                                   surface_emissivity = 0.98,
                                   surface_albedo = 0.1,
                                   solar_constant = 1361)        # W/m²

RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
└── diffuse_surface_albedo: ConstantField(0.1)

Atmosphere model

Build the atmosphere model with saturation adjustment microphysics.

clock = Clock(time=DateTime(1950, 11, 1, 12, 0, 0))
microphysics = SaturationAdjustment(equilibrium = WarmPhaseEquilibrium())
model = AtmosphereModel(grid; clock, dynamics, microphysics, radiation)

AtmosphereModel{CPU, RectilinearGrid}(time = 1950-11-01T12:00:00, iteration = 0)
├── grid: 1×1×64 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded} on Oceananigans.Architectures.CPU with 0×0×3 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── timestepper: RungeKutta3TimeStepper
├── advection scheme: 
│   ├── momentum: Centered(order=2)
│   ├── ρθ: Centered(order=2)
│   └── ρqᵗ: Centered(order=2)
├── tracers: ()
├── coriolis: Nothing
└── microphysics: SaturationAdjustment

Initial condition: idealized tropical profile with a cloud

We prescribe a simple tropical-like temperature profile with a moist boundary layer and a cloud between 1-2 km altitude.

θ₀ = reference_state.potential_temperature
q₀ = 0.015    # surface specific humidity (kg/kg)
Hᵗ = 2500     # moisture scale height (m)
qᵗᵢ(z) = q₀ * exp(-z / Hᵗ)

set!(model; θ=θ₀, qᵗ=qᵗᵢ)

Visualization

After set!, the radiation has been computed. We build Fields and AbstractOperations to visualize the atmospheric state and radiative fluxes.

T = model.temperature
pᵣ = reference_state.pressure
qᵗ = model.specific_moisture
qˡ = model.microphysical_fields.qˡ
ℋ = RelativeHumidityField(model)

ℐ_lw_up = radiation.upwelling_longwave_flux
ℐ_lw_dn = radiation.downwelling_longwave_flux
ℐ_sw = radiation.downwelling_shortwave_flux
ℐ_net = ℐ_lw_up + ℐ_lw_dn + ℐ_sw

set_theme!(fontsize=14, linewidth=3)
fig = Figure(size=(1200, 400), fontsize=14)

ax_T = Axis(fig[2, 1]; xlabel="Temperature, T (K)", ylabel="Altitude (km)")
ax_p = Axis(fig[2, 2]; xlabel="Pressure, p (hPa)")
ax_q = Axis(fig[2, 3]; xlabel="Specific humidity, q (kg/kg)")
ax_H = Axis(fig[2, 4]; xlabel="Relative humidity, ℋ (%)")
ax_I = Axis(fig[2, 5:6], xlabel="Radiation intensity, ℐ (W/m²)",
            ylabel="Altitude (km)", yaxisposition=:right)

[hideydecorations!(ax, grid=false) for ax in (ax_p, ax_q, ax_H)]
hidespines!(ax_T, :r, :t)
hidespines!(ax_p, :l, :r, :t)
hidespines!(ax_q, :l, :r, :t)
hidespines!(ax_H, :l, :r, :t)
hidespines!(ax_I, :l, :t)


lines!(ax_T, T)
lines!(ax_p, pᵣ / 100)  # Convert Pa to hPa

lines!(ax_q, qᵗ; label="qᵗ (total)")
lines!(ax_q, qˡ; label="qˡ (liquid)")
axislegend(ax_q, position=:rt, framevisible=false)

lines!(ax_H, 100ℋ)  # Convert to %

Lines{Tuple{Vector{Point{2, Float64}}}}

All radiation fluxes in one panel (positive = upward, negative = downward)

lines!(ax_I, ℐ_lw_up; label="LW ↑")
lines!(ax_I, ℐ_lw_dn; label="LW ↓")
lines!(ax_I, ℐ_sw; linestyle=:dash, label="SW ↓")
lines!(ax_I, ℐ_net; linewidth=4, alpha=0.5, color=:black, label="Net")

Legend(fig[1, 6], ax_I, orientation=:horizontal, nbanks=2, framevisible=false)

title = "Single Column Gray Radiation with O'Gorman & Schneider (2008) optical thickness"
fig[1, :] = Label(fig, title, fontsize=18, tellwidth=false)

fig

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD EPYC 7R13 Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_GPG = 3673DF529D9049477F76B37566E3C7DC03D6E495
  JULIA_LOAD_PATH = :@breeze
  JULIA_VERSION = 1.12.2
  JULIA_DEPOT_PATH = /usr/local/share/julia:
  JULIA_PATH = /usr/local/julia
  JULIA_PROJECT = @breeze

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

import Pkg
Pkg.status()

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

This page was generated using Literate.jl.