Single column radiation (gray, clear-sky, and all-sky)
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), and compare against clear-sky full-spectrum gas optics, doubled CO₂, and all-sky (cloudy) radiation.
using Breezeusing Oceananigans.Unitsusing CairoMakieusing NCDatasets # For RRTMGP lookup tablesusing RRTMGPGrid and thermodynamics
We create a single column spanning 20 km with 64 layers at a particular place.
Nz = 64λ, φ = -76.13, 39.48grid = RectilinearGrid(size=Nz, x=λ, y=φ, z=(0, 20kilometers), topology=(Flat, Flat, Bounded))1×1×64 RectilinearGrid{Float64, Flat, Flat, Bounded} on 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.5Set up the thermodynamic constants and reference state.
surface_temperature = 300constants = 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 materializedRadiative transfer models
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. We also create clear-sky full-spectrum models with present-day and doubled CO₂ concentrations.
using Datesgray_radiation = RadiativeTransferModel(grid, GrayOptics(), constants; surface_temperature, surface_emissivity = 0.98, surface_albedo = 0.1, solar_constant = 1361) # W/m²RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── solar_position: ApparentSolarPosition(coordinate=(-76.13, 39.48), epoch=<from clock>)
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
└── diffuse_surface_albedo: ConstantField(0.1)Clear-sky with default CO₂ (~420 ppm)
clear_sky_radiation = RadiativeTransferModel(grid, ClearSkyOptics(), constants; surface_temperature, surface_emissivity = 0.98, surface_albedo = 0.1, solar_constant = 1361) # W/m²RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── solar_position: ApparentSolarPosition(coordinate=(-76.13, 39.48), epoch=<from clock>)
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
└── diffuse_surface_albedo: ConstantField(0.1)Clear-sky with doubled CO₂ (~840 ppm) to show the radiative forcing effect
high_co2_atmosphere = BackgroundAtmosphere(CO₂ = 840e-6)high_co2_radiation = RadiativeTransferModel(grid, ClearSkyOptics(), constants; background_atmosphere = high_co2_atmosphere, surface_temperature, surface_emissivity = 0.98, surface_albedo = 0.1, solar_constant = 1361) # W/m²RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── solar_position: ApparentSolarPosition(coordinate=(-76.13, 39.48), epoch=<from clock>)
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
└── diffuse_surface_albedo: ConstantField(0.1)All-sky with cloud scattering optics
all_sky_radiation = RadiativeTransferModel(grid, AllSkyOptics(), constants; surface_temperature, surface_emissivity = 0.98, surface_albedo = 0.1, solar_constant = 1361, liquid_effective_radius = ConstantRadiusParticles(10e-6), ice_effective_radius = ConstantRadiusParticles(30e-6))RadiativeTransferModel
├── solar_constant: 1361.0 W m⁻²
├── solar_position: ApparentSolarPosition(coordinate=(-76.13, 39.48), epoch=<from clock>)
├── surface_temperature: ConstantField(300.0) K
├── surface_emissivity: ConstantField(0.98)
├── direct_surface_albedo: ConstantField(0.1)
├── liquid_effective_radius: Breeze.AtmosphereModels.ConstantRadiusParticles{Float64}(1.0e-5)
├── ice_effective_radius: Breeze.AtmosphereModels.ConstantRadiusParticles{Float64}(3.0e-5)
└── diffuse_surface_albedo: ConstantField(0.1)Atmosphere models
Build the atmosphere models with saturation adjustment microphysics.
clock = Clock(time=DateTime(1950, 11, 1, 12, 0, 0))microphysics = SaturationAdjustment(equilibrium = WarmPhaseEquilibrium())gray_model = AtmosphereModel(grid; clock, dynamics, microphysics, radiation=gray_radiation)clear_sky_model = AtmosphereModel(grid; clock, dynamics, microphysics, radiation=clear_sky_radiation)high_co2_model = AtmosphereModel(grid; clock, dynamics, microphysics, radiation=high_co2_radiation)all_sky_model = AtmosphereModel(grid; clock, dynamics, microphysics, radiation=all_sky_radiation)AtmosphereModel{CPU, RectilinearGrid}(time = 1950-11-01T12:00:00, iteration = 0)
├── grid: 1×1×64 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── dynamics: AnelasticDynamics(p₀=101325.0, θ₀=300.0)
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme:
│ ├── momentum: Centered(order=2)
│ ├── ρθ: Centered(order=2)
│ └── ρqᵉ: Centered(order=2)
├── forcing: @NamedTuple{ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵉ::Returns{Float64}, ρe::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: SaturationAdjustmentInitial condition: idealized tropical profile with a cloud
We prescribe a simple tropical-like temperature profile with a moist boundary layer. To produce clouds for the all-sky comparison, we use high moisture that will saturate in the lower troposphere via saturation adjustment.
θ₀ = reference_state.potential_temperatureq₀ = 0.020 # surface specific humidity (kg/kg) - high enough to saturateHᵗ = 3000 # moisture scale height (m)qᵗᵢ(z) = q₀ * exp(-z / Hᵗ)set!(gray_model; θ=θ₀, qᵗ=qᵗᵢ)set!(clear_sky_model; θ=θ₀, qᵗ=qᵗᵢ)set!(high_co2_model; θ=θ₀, qᵗ=qᵗᵢ)set!(all_sky_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 = gray_model.temperaturepᵣ = reference_state.pressureqᵛ = specific_humidity(gray_model)ℋ = RelativeHumidityField(gray_model)ℐ_lw_up_gray = gray_radiation.upwelling_longwave_fluxℐ_lw_dn_gray = gray_radiation.downwelling_longwave_fluxℐ_sw_gray = gray_radiation.downwelling_shortwave_fluxℐ_net_gray = ℐ_lw_up_gray + ℐ_lw_dn_gray + ℐ_sw_grayℐ_lw_up_clear = clear_sky_radiation.upwelling_longwave_fluxℐ_lw_dn_clear = clear_sky_radiation.downwelling_longwave_fluxℐ_sw_clear = clear_sky_radiation.downwelling_shortwave_fluxℐ_net_clear = ℐ_lw_up_clear + ℐ_lw_dn_clear + ℐ_sw_clearℐ_lw_up_2xco2 = high_co2_radiation.upwelling_longwave_fluxℐ_lw_dn_2xco2 = high_co2_radiation.downwelling_longwave_fluxℐ_sw_2xco2 = high_co2_radiation.downwelling_shortwave_fluxℐ_net_2xco2 = ℐ_lw_up_2xco2 + ℐ_lw_dn_2xco2 + ℐ_sw_2xco2ℐ_lw_up_allsky = all_sky_radiation.upwelling_longwave_fluxℐ_lw_dn_allsky = all_sky_radiation.downwelling_longwave_fluxℐ_sw_allsky = all_sky_radiation.downwelling_shortwave_fluxℐ_net_allsky = ℐ_lw_up_allsky + ℐ_lw_dn_allsky + ℐ_sw_allskyMultiaryOperation at (Center, Center, Face)
├── grid: 1×1×64 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
└── tree:
+ at (Center, Center, Face)
├── 1×1×65 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Face} on Oceananigans.Grids.RectilinearGrid on CPU
├── 1×1×65 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Face} on Oceananigans.Grids.RectilinearGrid on CPU
└── 1×1×65 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Face} on Oceananigans.Grids.RectilinearGrid on CPUGet cloud liquid for visualization
qˡ = all_sky_model.microphysical_fields.qˡset_theme!(fontsize=14, linewidth=2.5)Format altitude ticks in km (but keep internal units in meters).
z_ticks_km = 0:5:20z_ticks_m = ((z_ticks_km .* 1000), string.(z_ticks_km))fig = Figure(size=(1600, 800), fontsize=14)Atmospheric state panels (top row)
ax_T = Axis(fig[1, 1]; xlabel="Temperature (K)", ylabel="Altitude (km)", yticks=z_ticks_m, xticks=200:25:300)ax_q = Axis(fig[1, 2]; xlabel="Specific humidity (kg/kg)", yticks=z_ticks_m)ax_H = Axis(fig[1, 3]; xlabel="Relative humidity (%)", yticks=z_ticks_m)ax_ql = Axis(fig[1, 4]; xlabel="Cloud liquid (g/kg)", yticks=z_ticks_m)Axis with 0 plots:
Radiation panels (bottom row) - one per component
ax_lw_up = Axis(fig[2, 1]; xlabel="LW ↑ (W/m²)", ylabel="Altitude (km)", yticks=z_ticks_m)ax_lw_dn = Axis(fig[2, 2]; xlabel="LW ↓ (W/m²)", yticks=z_ticks_m)ax_sw_dn = Axis(fig[2, 3]; xlabel="SW ↓ (W/m²)", yticks=z_ticks_m)ax_net = Axis(fig[2, 4]; xlabel="Net flux (W/m²)", yticks=z_ticks_m)Axis with 0 plots:
Hide y-axis decorations on inner panels
[hideydecorations!(ax, grid=false) for ax in (ax_q, ax_H, ax_ql, ax_lw_dn, ax_sw_dn, ax_net)]6-element Vector{Bool}:
0
0
0
0
0
0Atmospheric state
lines!(ax_T, T; color=:gray30)lines!(ax_q, qᵛ; color=:gray30)lines!(ax_H, 100ℋ; color=:gray30)lines!(ax_ql, 1000qˡ; color=:lime) # Convert to g/kgLines{Tuple{Vector{Point{2, Float64}}}}Colors for radiation schemes
c_gray = :blackc_clear = :dodgerbluec_2xco2 = :orangeredc_allsky = :lime:limeLW upwelling (positive)
lines!(ax_lw_up, ℐ_lw_up_gray; color=c_gray)lines!(ax_lw_up, ℐ_lw_up_clear; color=c_clear)lines!(ax_lw_up, ℐ_lw_up_2xco2; color=c_2xco2)lines!(ax_lw_up, ℐ_lw_up_allsky; color=c_allsky)Lines{Tuple{Vector{Point{2, Float64}}}}LW downwelling (negative, so we negate for display)
lines!(ax_lw_dn, -ℐ_lw_dn_gray; color=c_gray)lines!(ax_lw_dn, -ℐ_lw_dn_clear; color=c_clear)lines!(ax_lw_dn, -ℐ_lw_dn_2xco2; color=c_2xco2)lines!(ax_lw_dn, -ℐ_lw_dn_allsky; color=c_allsky)Lines{Tuple{Vector{Point{2, Float64}}}}SW downwelling (negative, so we negate for display)
lines!(ax_sw_dn, -ℐ_sw_gray; color=c_gray)lines!(ax_sw_dn, -ℐ_sw_clear; color=c_clear)lines!(ax_sw_dn, -ℐ_sw_2xco2; color=c_2xco2)lines!(ax_sw_dn, -ℐ_sw_allsky; color=c_allsky)Lines{Tuple{Vector{Point{2, Float64}}}}Net flux
lines!(ax_net, ℐ_net_gray; color=c_gray)lines!(ax_net, ℐ_net_clear; color=c_clear)lines!(ax_net, ℐ_net_2xco2; color=c_2xco2)lines!(ax_net, ℐ_net_allsky; color=c_allsky)Lines{Tuple{Vector{Point{2, Float64}}}}Legend
scheme_handles = [ LineElement(color=c_gray, linewidth=3), LineElement(color=c_clear, linewidth=3), LineElement(color=c_2xco2, linewidth=3), LineElement(color=c_allsky, linewidth=3),]scheme_labels = ["Gray", "Clear-sky (420 ppm)", "2×CO₂ (840 ppm)", "All-sky (cloudy)"]Legend(fig[0, :], scheme_handles, scheme_labels; orientation=:horizontal, framevisible=false, tellwidth=false)fig
Heating rates
The RadiativeTransferModel automatically computes the heating tendency Q = -dF_net/dz (W/m³) from the radiative flux divergence. We convert to K/day using dT/dt = Q / (ρ cₚ).
Q_gray = gray_radiation.flux_divergenceQ_clear = clear_sky_radiation.flux_divergenceQ_2xco2 = high_co2_radiation.flux_divergenceQ_allsky = all_sky_radiation.flux_divergence1×1×64 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on CPU
├── grid: 1×1×64 RectilinearGrid{Float64, Flat, Flat, Bounded} on CPU with 0×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│ └── west: Nothing, east: Nothing, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 1×1×70 OffsetArray(::Array{Float64, 3}, 1:1, 1:1, -2:67) with eltype Float64 with indices 1:1×1:1×-2:67
└── max=0.0579947, min=-0.0121296, mean=0.00416759Convert W/m³ → K/day: Q / (ρᵣ cᵖᵈ) × 86400
ρᵣ = reference_state.densitycᵖᵈ = constants.dry_air.heat_capacity / constants.dry_air.molar_mass # J/(kg·K)to_K_per_day = 86400 / cᵖᵈfig2 = Figure(size=(800, 500), fontsize=14)ax_Q = Axis(fig2[1, 1]; xlabel="Heating rate (K/day)", ylabel="Altitude (km)", yticks=z_ticks_m, title="Radiative heating rates")lines!(ax_Q, to_K_per_day * Q_gray / ρᵣ; color=c_gray, label="Gray")lines!(ax_Q, to_K_per_day * Q_clear / ρᵣ; color=c_clear, label="Clear-sky (420 ppm)")lines!(ax_Q, to_K_per_day * Q_2xco2 / ρᵣ; color=c_2xco2, label="2×CO₂ (840 ppm)")lines!(ax_Q, to_K_per_day * Q_allsky / ρᵣ; color=c_allsky, label="All-sky (cloudy)")vlines!(ax_Q, 0; color=:gray50, linestyle=:dash, linewidth=1)axislegend(ax_Q, position=:lt)fig2
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.