Extension to RRTMGP.jl and implement RadiativeTransferModel

AGENTS.md

@@ -21,10 +21,10 @@ Breeze interfaces with ClimaOcean for coupled atmosphere-ocean simulations.
 1. **Explicit Imports**: Use `ExplicitImports.jl` style - explicitly import all used functions/types
    - Import from Oceananigans explicitly (already done in src/Breeze.jl)
    - Tests automatically check for proper imports
-   
+
 2. **Type Stability**: Prioritize type-stable code for performance
    - All structs must be concretely typed
-   
+
 3. **Kernel Functions**: For GPU compatibility:
    - Use KernelAbstractions.jl syntax for kernels, eg `@kernel`, `@index`
    - Keep kernels type-stable and allocation-free
@@ -33,10 +33,13 @@ Breeze interfaces with ClimaOcean for coupled atmosphere-ocean simulations.
    - Do not put error messages inside kernels.
    - Models _never_ go inside kernels
    - Mark functions inside kernels with `@inline`.
-
+   - **Never use loops outside kernels**: Always replace `for` loops that iterate over grid points
+     with kernels launched via `launch!`. This ensures code works on both CPU and GPU.
+
 4. **Documentation**:
    - Use DocStringExtensions.jl for consistent docstrings
-   - Include `$(SIGNATURES)` for automatic signature documentation
+   - Use `$(TYPEDSIGNATURES)` for automatic typed signature documentation (preferred over `$(SIGNATURES)`)
+   - Never write explicit function signatures in docstrings; always use `$(TYPEDSIGNATURES)`
    - Add examples in docstrings when helpful
 
 5. **Memory leanness**
@@ -75,8 +78,8 @@ Breeze interfaces with ClimaOcean for coupled atmosphere-ocean simulations.
   - `TitleCase` style is reserved for types, type aliases, and constructors.
   - `snake_case` style should be used for functions and variables (instances of types)
   - "Number variables" (`Nx`, `Ny`) should start with capital `N`. For number of time steps use `Nt`.
-    Spatial indices are `i, j, k` and time index is `n`. 
-  
+    Spatial indices are `i, j, k` and time index is `n`.
+
 2. **Import/export style**
   - Write all exported names at the top of a module file, before any other code.
   - For explicit imports, import Oceananigans and Breeze names first. Then write imports for "external" packages.
@@ -125,6 +128,12 @@ Breeze interfaces with ClimaOcean for coupled atmosphere-ocean simulations.
   - The examples and docs have their own `Project.toml` environment. When your run examples you need to use `examples/Project.toml`.
     When you build new examples, please add example-specific packages to `examples/Project.toml`. Do not add example-specific
     packages to the main Breeze Project.toml. You may also need to add relevant packages to AtmosphereProfilesLibrary.
+  - When making plots, do not use `interior(field, i, j, k)` to make a plot. Instead either pass `field`
+    directly, as in `lines!(ax, field)` or `lines!(ax, z, field)` for a 1D plot with either automatic
+    or custom vertical coordinate. You only need to provide the coordinate if it has different units, eg
+    if you have converted z to kilometers. 1D fields work with `lines` and 2D fields work with 2D plots
+    like `heatmap` or `contourf`. For 3D fields in a 2D plane, use `view(field, :, :, k)`
+    (e.g. for a xy-slice).
   - **CRITICAL - Plotting Fields**: NEVER use `interior(field, ...)` for plotting. Makie/CairoMakie
     can plot `Field` objects directly, which automatically handles coordinates and is cleaner.
     Use `view(field, i, j, k)` to window fields if needed. This applies to both static plots and
@@ -133,7 +142,7 @@ Breeze interfaces with ClimaOcean for coupled atmosphere-ocean simulations.
     # WRONG - do not do this:
     data = @lift interior(field_ts[$n], :, 1, :)
     heatmap!(ax, x, z, data, ...)
-    
+
     # CORRECT - pass Field directly:
     field_n = @lift field_ts[$n]
     heatmap!(ax, field_n, ...)
@@ -180,6 +189,7 @@ Breeze interfaces with ClimaOcean for coupled atmosphere-ocean simulations.
 - **Functions**: snake_case (e.g., `update_atmosphere_model!`, `compute_pressure!`)
 - **Kernels**: "Kernels" (functions prefixed with `@kernel`) may be prefixed with an underscore (e.g., `_kernel_function`)
 - **Variables**: Use _either_ an English long name, or mathematical notation with readable unicode. Variable names should be taken from `docs/src/appendix/notation.md` in the docs. If a new variable is created (or if one doesn't exist), it should be added to the table in notation.md
+- **Avoid abbreviations**: Use full words instead of abbreviations. For example, use `latitude` instead of `lat`, `longitude` instead of `lon`, `temperature` instead of `temp`. This improves code readability and self-documentation.
 
 ### Breeze Module Structure
 ```

Project.toml

@@ -5,27 +5,34 @@ authors = ["NumericalEarth organization (github.com/NumericalEarth) and lovely c
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [weakdeps]
+ClimaComms = "3a4d1b5c-c61d-41fd-a00a-5873ba7a1b0d"
 CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
+RRTMGP = "a01a1ee8-cea4-48fc-987c-fc7878d79da1"
 
 [workspace]
 projects = ["docs", "test"]
 
 [extensions]
 BreezeCloudMicrophysicsExt = "CloudMicrophysics"
+BreezeRRTMGPExt = ["ClimaComms", "RRTMGP"]
 
 [compat]
 Adapt = "4.3.0"
+ClimaComms = "0.6"
 CloudMicrophysics = "0.29.0"
+Dates = "<0.0.1, 1"
 DocStringExtensions = "0.9.5"
 JLD2 = "0.5.13, 0.6"
 KernelAbstractions = "0.9"
 Oceananigans = "0.102.2"
 Printf = "1"
+RRTMGP = "0.21"
 julia = "1.10.10"

docs/Project.toml

@@ -1,24 +1,25 @@
 [deps]
 AtmosphericProfilesLibrary = "86bc3604-9858-485a-bdbe-831ec50de11d"
 Breeze = "660aa2fb-d4c8-4359-a52c-9c057bc511da"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
-CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+RRTMGP = "a01a1ee8-cea4-48fc-987c-fc7878d79da1"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [sources]
 Breeze = {path = ".."}
 
 [compat]
+CUDA = "5"
 CairoMakie = "0.15"
 CloudMicrophysics = "0.29.0"
-CUDA = "5"
 Documenter = "1"
 DocumenterCitations = "1.4"
 InteractiveUtils = "<0.0.1, 1"

docs/make.jl

@@ -42,6 +42,7 @@ example_pages = [
     "Precipitating shallow cumulus (RICO)" => "literated/rico.md",
     "Convection over prescribed sea surface temperature (SST)" => "literated/prescribed_sea_surface_temperature.md",
     "Inertia gravity wave" => "literated/inertia_gravity_wave.md",
+    "Single column gray radiation" => "literated/single_column_radiation.md",
 ]
 
 literate_code(script_path, literated_dir) = """
@@ -94,6 +95,7 @@ makedocs(
                 "Microphysics Interface" => "developer/microphysics_interface.md",
             ],
         ],
+        "Radiative Transfer" => "radiative_transfer.md",
         "Dycore equations and algorithms" => "dycore_equations_algorithms.md",
         "Appendix" => Any[
             "Notation" => "appendix/notation.md",

docs/src/api.md

@@ -63,6 +63,13 @@ Modules = [Breeze.BoundaryConditions]
 Private = false
 ```
 
+### CelestialMechanics
+
+```@autodocs
+Modules = [Breeze.CelestialMechanics]
+Private = false
+```
+
 ## Private API
 
 ```@autodocs
@@ -111,3 +118,10 @@ Public = false
 Modules = [Breeze.BoundaryConditions]
 Public = false
 ```
+
+### CelestialMechanics
+
+```@autodocs
+Modules = [Breeze.CelestialMechanics]
+Public = false
+```

docs/src/appendix/notation.md

@@ -56,6 +56,9 @@ The following table also uses a few conventions that suffuse the source code and
 | ``qᵛ⁺``                             | `qᵛ⁺`  |                                     | Saturation specific humidity over a surface                                    |
 | ``qᵛ⁺ˡ``                            | `qᵛ⁺ˡ` |                                     | Saturation specific humidity over a planar liquid surface                      |
 | ``qᵛ⁺ⁱ``                            | `qᵛ⁺ⁱ` |                                     | Saturation specific humidity over a planar ice surface                         |
+| ``pᵛ``                              | `pᵛ`   |                                     | Vapor pressure (partial pressure of water vapor), ``pᵛ = ρ qᵛ Rᵛ T``           |
+| ``pᵛ⁺``                             | `pᵛ⁺`  |                                     | Saturation vapor pressure                                                      |
+| ``\mathscr{H}``                     | `ℋ`    | `RelativeHumidity(model)`           | Relative humidity, ``ℋ = pᵛ / pᵛ⁺``                                            |
 | ``g``                               | `g`    | `TC.gravitational_acceleration`     | Gravitational acceleration                                                     |
 | ``\mathcal{R}``                     | `ℛ`    | `TC.molar_gas_constant`             | Universal (molar) gas constant                                                 |
 | ``Tᵗʳ``                             | `Tᵗʳ`  | `TC.triple_point_temperature`       | Temperature at the vapor-liquid-ice triple point                               |
@@ -101,3 +104,5 @@ The following table also uses a few conventions that suffuse the source code and
 | ``Cᵛ``                              | `Cᵛ`   |                                     | Surface vapor transfer coefficient (Dalton number)                             |
 | ``T_0``                             | `T₀`   |                                     | Sea surface temperature                                                        |
 | ``qᵛ₀``                             | `qᵛ₀`  |                                     | Saturation specific humidity at sea surface                                    |
+| ``\mathscr{I}``                     | `ℐ`    |                                     | Radiative flux (intensity), W/m²                                               |
+| ``F_{\mathscr{I}}``                 | `Fℐ`   |                                     | Radiative flux divergence (heating rate), K/s                                  |

docs/src/breeze.bib

@@ -154,6 +154,28 @@ @article{Shu09
   pages = {82--126}
 }
 
+@article{OGormanSchneider2008,
+  author = {O'Gorman, Paul A. and Schneider, Tapio},
+  title = {The hydrological cycle over a wide range of climates simulated with an idealized {GCM}},
+  journal = {Journal of Climate},
+  volume = {21},
+  number = {15},
+  pages = {3815--3832},
+  year = {2008},
+  doi = {10.1175/2007JCLI2065.1}
+}
+
+@article{Schneider2004,
+  author = {Schneider, Tapio},
+  title = {The tropopause and the thermal stratification in the extratropics of a dry atmosphere},
+  journal = {Journal of the Atmospheric Sciences},
+  volume = {61},
+  number = {12},
+  pages = {1317--1340},
+  year = {2004},
+  doi = {10.1175/1520-0469(2004)061<1317:TTATTS>2.0.CO;2}
+}
+
 @article{vanZanten2011,
   author = {van Zanten, Margreet C. and Stevens, Bjorn and Nuijens, Louise and Siebesma, A. Pier and Ackerman, Andrew S. and Burnet, Frédéric and Cheng, Anke and Couvreux, Fleur and Jiang, Hongli and Khairoutdinov, Marat and Kogan, Yefim and Lewellen, David C. and Mechem, David and Nakamura, Kozo and Noda, Akira and Shipway, Ben J. and Slawinska, Joanna and Wang, Shuyi and Wyszogrodzki, Andrzej},
   title = {Controls on precipitation and cloudiness in simulations of trade-wind cumulus as observed during {RICO}},

docs/src/radiative_transfer.md

@@ -0,0 +1,124 @@
+# Radiative Transfer
+
+Breeze.jl integrates with [RRTMGP.jl](https://github.com/CliMA/RRTMGP.jl) to provide radiative transfer capabilities for atmospheric simulations. The radiative transfer model computes longwave and shortwave radiative fluxes, which can be incorporated into energy tendency equations.
+
+## Gray Atmosphere Radiation
+
+The simplest radiative transfer option is gray atmosphere radiation, which uses the optical thickness parameterization from [Schneider2004](@citet) and [OGormanSchneider2008](@citet). This approximation treats the atmosphere as having a single effective absorption coefficient rather than computing full spectral radiation.
+
+### Basic Usage
+
+To use gray radiation in a Breeze simulation, create a `GrayRadiativeTransferModel` model and pass it to the [`AtmosphereModel`](@ref) constructor:
+
+```@example
+using Breeze
+using Oceananigans.Units
+using Dates
+using RRTMGP.AtmosphericStates: GrayOpticalThicknessOGorman2008
+
+Nz = 64
+λ, φ = -70.9, 42.5  # longitude, latitude
+grid = RectilinearGrid(size=Nz, x=λ, y=φ, z=(0, 20kilometers),
+                       topology=(Flat, Flat, Bounded))
+
+# Thermodynamic setup
+surface_temperature = 300
+constants = ThermodynamicConstants()
+
+reference_state = ReferenceState(grid, constants;
+                                 surface_pressure = 101325,
+                                 potential_temperature = surface_temperature)
+
+formulation = AnelasticFormulation(reference_state,
+                                   thermodynamics = :LiquidIcePotentialTemperature)
+
+# Create gray radiation model
+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²
+
+# Create atmosphere model with DateTime clock for solar position
+clock = Clock(time=DateTime(2024, 9, 27, 16, 0, 0))
+model = AtmosphereModel(grid; clock, formulation, radiation)
+```
+
+When a `DateTime` clock is used, the solar zenith angle is computed automatically from the time and grid location (longitude and latitude).
+
+### Gray Radiation Model
+
+The [`RadiativeTransferModel`](@ref) model computes:
+
+- **Longwave radiation**: Both upwelling and downwelling thermal radiation using RRTMGP's two-stream solver
+- **Shortwave radiation**: Direct beam solar radiation
+
+The gray atmosphere optical thickness follows the parameterization by [OGormanSchneider2008](@citet):
+
+```math
+τ_{lw} = α \frac{Δp}{p} \left[ f_l \frac{p}{p_0} + 4 (1 - f_l) \left(\frac{p}{p_0}\right)^4 \right] \left[ τ_e + (τ_p - τ_e) \sin^2 φ \right]
+```
+
+where ``φ`` is latitude and ``α``, ``f_l``, ``τ_e``, ``τ_p`` are empirical parameters.
+
+For shortwave:
+```math
+τ_{sw} = 2 τ_0 \frac{p}{p_0} \frac{Δp}{p_0}
+```
+
+where ``τ_0 = 0.22`` is the shortwave optical depth parameter.
+
+### Radiative Fluxes
+
+After running [`set!`](@ref), the radiative fluxes are available from the radiation model:
+
+```julia
+# Longwave fluxes (ZFaceFields)
+ℐ_lw_up = radiation.upwelling_longwave_flux
+ℐ_lw_dn = radiation.downwelling_longwave_flux
+
+# Shortwave flux (direct beam only for non-scattering solver)
+ℐ_sw = radiation.downwelling_shortwave_flux
+```
+
+!!! note "Shortwave Radiation"
+    The gray atmosphere uses a non-scattering shortwave approximation, so only
+    the direct beam flux is computed. There is no diffuse shortwave or upwelling
+    shortwave in this model.
+
+### Solar Zenith Angle
+
+When using a `DateTime` clock, the solar zenith angle is computed from:
+- Grid location (longitude from `x`, latitude from `y` for single-column grids)
+- Date and time from `model.clock.time`
+
+The calculation accounts for:
+- Day of year (for solar declination)
+- Hour angle (based on solar time)
+- Latitude (for observer position)
+
+## Surface Properties
+
+The `RadiativeTransferModel` model requires surface properties:
+
+| Property | Description | Typical Values |
+|----------|-------------|----------------|
+| `surface_temperature` | Temperature at the surface [K] | 280-310 |
+| `surface_emissivity` | Longwave emissivity (0-1) | 0.95-0.99 |
+| `surface_albedo` | Shortwave albedo (0-1) | 0.1-0.3 |
+| `solar_constant` | TOA solar flux [W/m²] | 1361 |
+
+## Integration with dynamics
+
+Radiative fluxes can be used to compute heating rates for the energy equation. The radiative heating rate is computed from flux divergence:
+
+```math
+F_{\mathscr{I}} = -\frac{1}{\rho cᵖᵐ} \frac{\partial \mathscr{I}_{net}}{\partial z}
+```
+
+where ``\mathscr{I}_{net}`` is the net radiative flux (upwelling minus downwelling), ``cᵖᵐ`` is the mixture heat capacity, and ``F_{\mathscr{I}}`` is the radiative flux divergence (heating rate).
+
+## Architecture Support
+
+The radiative transfer implementation supports both CPU and GPU architectures. The column-based RRTMGP solver is called from Oceananigans' field data arrays with appropriate data layout conversions.

examples/Project.toml

@@ -7,6 +7,7 @@ GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+RRTMGP = "a01a1ee8-cea4-48fc-987c-fc7878d79da1"
 
 [sources]
 Breeze = {path = ".."}