Add two-stream radiation model with simplified surface energy balance (force-restore)

CMake/BuildERFExe.cmake

@@ -517,6 +517,8 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/PBL/ERF_ComputeDiffusivityYSU.cpp
        ${SRC_DIR}/PBL/ERF_ComputeDiffusivityMRF.cpp
        ${SRC_DIR}/PBL/ERF_ComputeDiffusivityYSUNew.cpp
+       ${SRC_DIR}/Radiation/ERF_RadiationDiagnostics.cpp
+       ${SRC_DIR}/Radiation/ERF_AdvanceTwoStreamRadiation.cpp
        ${SRC_DIR}/Refinement/ERF_Tagging.cpp
        ${SRC_DIR}/Refinement/ERF_RefineBox.cpp
        ${SRC_DIR}/Refinement/ERF_RefineHurricane.cpp
@@ -627,6 +629,7 @@ function(build_erf_lib erf_lib_name)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Source/IO>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Source/LinearSolvers>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Source/PBL>)
+  target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Source/Radiation>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Source/SourceTerms>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Source/TimeIntegration>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Source/Utils>)

Docs/sphinx_doc/index.rst

@@ -62,6 +62,7 @@ In addition to this documentation, there is API documentation for ERF generated
    theory/DNSvsLES.rst
    theory/PBLschemes.rst
    theory/Forcings.rst
+   theory/Radiation.rst
    Particles.rst
    SuperDroplets.rst
    ProblemInputs.rst

Exec/CanonicalTests/Radiation/LW_Isothermal/README.md

@@ -0,0 +1,53 @@
+# Longwave Isothermal Column Test
+
+## Objective
+
+Validate the gray-gas longwave solver on a physically isothermal column: an atmosphere at a
+uniform absolute temperature over a black surface at the same temperature.
+
+## Test Design
+
+### Configuration
+
+- **Domain**: 3000 m × 3000 m horizontal, 1024 m vertical (64 layers)
+- **Time**: one slow step (0.25 s)
+- **Sounding**: `theta(z) = T0 exp(g z / (c_p T0))` with `T0 = 300 K`, so `T = theta * pi = T0`
+  at every level (the sounding is piecewise linear between five levels, accurate to ~0.05 K)
+- **Surface**: `surface_temp_k = 300`, `surface_emissivity_lw = 1`
+- **Longwave optical depth per layer**: `tau_lw = 1.0` (the 64-layer column is opaque)
+- **Shortwave**: disabled
+
+### Key Physics
+
+Every layer emits `sigma T0^4`. In an opaque isothermal column the upward flux is `sigma T0^4`
+on every interface and the downward flux approaches `sigma T0^4` at the surface, so
+
+```
+LW_up_TOA      = sigma T0^4              ≈ 459.3 W/m²
+LW_net_surface = sigma T0^4 exp(-tau_col) ≈ 0 W/m²
+```
+
+The column still cools to space (the downward flux vanishes at the top), so the heating rate is
+non-zero and strongest in the top layers. This case previously used an `isothermal_test` override
+that forced `F_up = F_down` and zero heating; it now exercises the real solver.
+
+## Files
+
+- `inputs` — control file
+- `input_sounding_lw_isothermal` — isothermal-temperature sounding
+- `check_flux_accuracy.py` — validation script
+- `radiation_lw_diag.dat` — reference diagnostics
+
+## Running the Test
+
+```bash
+cd Exec/CanonicalTests/Radiation/LW_Isothermal
+mpirun -np 1 erf.ex inputs
+python3 check_flux_accuracy.py
+```
+
+## Validation Criteria
+
+1. `LW_up_TOA` within 0.5% of `sigma T0^4`
+2. `|LW_net_surface|` below 0.5% of `sigma T0^4`
+3. `heating_rate_max` finite and non-zero

Exec/CanonicalTests/Radiation/LW_Isothermal/check_flux_accuracy.py

@@ -0,0 +1,184 @@
+#!/usr/bin/env python3
+"""
+Two-Stream Radiation Validation Script
+Longwave Isothermal Column
+
+The sounding prescribes theta(z) = T0 exp(g z / (c_p T0)) so that the
+absolute temperature is T0 = 300 K at every level, and the surface is a
+black body at the same T0. With a gray optical depth of 1 per layer over
+64 layers the column is opaque, so the gray-gas solution gives
+
+    LW_up_TOA      = sigma T0^4                (emission of the top layers)
+    LW_net_surface = sigma T0^4 exp(-tau_col)  ~ 0  (F_down(0) -> sigma T0^4)
+
+and every layer cools (cooling to space, strongest at the top), so
+heating_rate_max is non-zero. This replaces the former isothermal_test
+override, which forced these values instead of computing them.
+
+It reads radiation_lw_diag.dat and checks that:
+1. LW_up_TOA matches sigma*T0^4 within a small tolerance (the sounding is a
+   piecewise-linear approximation of the exponential theta profile)
+2. |LW_net_surface| is negligible compared to sigma*T0^4
+3. heating_rate_max is finite and non-zero (the column radiates to space)
+"""
+
+import sys
+import os
+import numpy as np
+import math
+
+def read_radiation_diag(filename):
+    """Read the radiation diagnostic CSV and return a dict of column lists.
+
+    The file is comma separated with a header line
+    (step,time,call_site,SW_surface,SW_TOA,SW_up_TOA,LW_net_surface,LW_up_TOA,heating_rate_max,...),
+    so columns are looked up by name rather than by position. Non-numeric
+    columns (call_site) are kept as strings; numeric columns are floats.
+    """
+    import csv
+    try:
+        with open(filename, 'r') as f:
+            reader = csv.DictReader(f)
+            if reader.fieldnames is None:
+                print(f"ERROR: No header found in {filename}")
+                return None
+            data = {name.strip(): [] for name in reader.fieldnames}
+            for row in reader:
+                if not any((v or '').strip() for v in row.values()):
+                    continue
+                for name in reader.fieldnames:
+                    key = name.strip()
+                    val = (row.get(name) or '').strip()
+                    if key == 'call_site':
+                        data[key].append(val)
+                    elif key == 'step':
+                        data[key].append(int(float(val)))
+                    else:
+                        try:
+                            data[key].append(float(val))
+                        except ValueError:
+                            data[key].append(float('nan'))
+    except IOError:
+        print(f"ERROR: Could not read {filename}")
+        return None
+
+    if not data.get('step'):
+        print(f"ERROR: No data found in {filename}")
+        return None
+
+    return data
+
+def compute_bb_radiation(T, sigma=5.670374419e-8):
+    """
+    Compute blackbody radiation intensity using Stefan-Boltzmann law.
+
+    Args:
+        T: Temperature [K]
+        sigma: Stefan-Boltzmann constant [W/(m^2·K^4)]
+
+    Returns:
+        Radiative intensity [W/m^2]
+    """
+    if T <= 0:
+        return 0.0
+    return sigma * T**4
+
+def check_lw_isothermal_accuracy():
+    """Check LW isothermal test accuracy."""
+
+    # Read diagnostic file
+    diag_file = "radiation_lw_diag.dat"
+    if not os.path.exists(diag_file):
+        print(f"ERROR: Diagnostic file {diag_file} not found")
+        return False
+
+    data = read_radiation_diag(diag_file)
+    if data is None:
+        return False
+
+    # Test parameters (must match inputs file and sounding)
+    T_iso_K = 300.0  # Isothermal temperature [K] (surface_temp_k and the sounding)
+    sigma = 5.670374419e-8  # Stefan-Boltzmann constant [W/(m^2·K^4)]
+
+    # Expected upwelling/downwelling flux (same in isothermal mode)
+    expected_flux = compute_bb_radiation(T_iso_K, sigma)
+
+    # Tolerance for numerical accuracy
+    # In isothermal mode with exact arithmetic:
+    #   - Upwelling and downwelling should be identical
+    #   - Heating rates should be exactly zero
+    # With finite precision, we allow ~1e-10 relative error
+    # The 5-point sounding approximates the exponential theta profile to
+    # ~0.05 K, i.e. ~0.1% in sigma T^4; allow 0.5%.
+    flux_tolerance = 5e-3
+    net_tolerance = 5e-3      # |LW_net_surface| / sigma T^4
+
+    print(f"\n{'='*70}")
+    print("Two-Stream Radiation: LW Isothermal Column Test")
+    print(f"{'='*70}")
+    print(f"\nTest Parameters:")
+    print(f"  Isothermal temperature T_iso = {T_iso_K:.2f} K")
+    print(f"  Stefan-Boltzmann constant σ = {sigma:.6e} W/(m^2·K^4)")
+    print(f"\nAnalytical Solution:")
+    print(f"  Expected LW_up_TOA = σ*T^4 = {expected_flux:.4f} W/m^2")
+    print(f"  Expected LW_net_surface ≈ 0 W/m^2 (opaque column at the surface temperature)")
+    print(f"  Expected heating_rate_max > 0 (cooling to space from the top layers)")
+
+    # Extract last timestep data
+    last_idx = -1
+    step = data['step'][last_idx]
+    time = data['time'][last_idx]
+    LW_net_surface = data['LW_net_surface'][last_idx]
+    LW_up_TOA = data['LW_up_TOA'][last_idx]
+    heating_rate_max = data['heating_rate_max'][last_idx]
+
+    print(f"\nComputed Values (step {step}, time {time:.4f}s):")
+    print(f"  Computed LW_up_TOA = {LW_up_TOA:.4f} W/m^2")
+    print(f"  Computed LW_net_surface = {LW_net_surface:.4f} W/m^2")
+    print(f"  Maximum heating rate = {heating_rate_max:.4e} K/s")
+
+    # Check results
+    errors = []
+
+    print(f"\nAccuracy Checks:")
+    # Isothermal column: the outgoing LW at the top is sigma*T^4.
+    if expected_flux > 0:
+        up_error = abs(LW_up_TOA - expected_flux) / expected_flux
+        print(f"  LW_up_TOA relative error: {up_error:.4e}", end="")
+        if up_error > flux_tolerance:
+            print(f" [FAIL - exceeds {flux_tolerance:.4e}]")
+            errors.append(f"LW_up_TOA error too large: {up_error:.4e}")
+        else:
+            print(" [PASS]")
+
+    # Isothermal column: up and down fluxes are equal, so the surface net LW is zero.
+    net_error = abs(LW_net_surface) / expected_flux
+    print(f"  |LW_net_surface| / sigma T^4: {net_error:.4e}", end="")
+    if net_error > net_tolerance:
+        print(f" [FAIL - exceeds {net_tolerance:.4e}]")
+        errors.append(f"LW_net_surface not negligible: {LW_net_surface:.4e} W/m^2")
+    else:
+        print(" [PASS]")
+
+    # The column radiates to space, so the heating (cooling) is finite and non-zero
+    print(f"  Heating rate magnitude: {abs(heating_rate_max):.4e} K/s", end="")
+    if not math.isfinite(heating_rate_max) or abs(heating_rate_max) <= 0.0:
+        print(" [FAIL - expected a finite, non-zero cooling rate]")
+        errors.append(f"heating_rate_max not finite/non-zero: {heating_rate_max}")
+    else:
+        print(" [PASS]")
+
+    # Overall result
+    print(f"\n{'='*70}")
+    if errors:
+        print("TEST FAILED")
+        for err in errors:
+            print(f"  - {err}")
+        return False
+    else:
+        print("TEST PASSED - isothermal column verified (LW_up_TOA = σT^4, LW_net_surface ≈ 0, cooling to space)")
+        return True
+
+if __name__ == "__main__":
+    success = check_lw_isothermal_accuracy()
+    sys.exit(0 if success else 1)

Exec/CanonicalTests/Radiation/LW_Isothermal/input_sounding_lw_isothermal

@@ -0,0 +1,6 @@
+1000. 300.0 0.0
+   0.0   300.000  0.0  15.0  0.0
+ 256.0   302.511  0.0  15.0  0.0
+ 512.0   305.042  0.0  15.0  0.0
+ 768.0   307.595  0.0  15.0  0.0
+1024.0   310.169  0.0  15.0  0.0

Exec/CanonicalTests/Radiation/LW_Isothermal/inputs

@@ -0,0 +1,123 @@
+# Two-Stream Radiation: Shortwave Clear-Sky Analytical Test
+# Tests Beer-Lambert direct-beam solar radiation without clouds or scattering
+# Verifies against analytical solution: F(z) = S0 * cos(zenith) * exp(-tau_cumulative / cos(zenith))
+erf.prob_name = "ABL"
+
+stop_time = 0.25
+amrex.fpe_trap_invalid = 0
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+geometry.prob_extent = 3000 3000 1024 
+amr.n_cell           =  8   8   64
+amr.max_grid_size_z = 128
+geometry.is_periodic = 1 1 0
+
+
+zlo.type = "surface_layer"
+erf.most.z0   = 0.1      # roughness length [m]
+erf.most.zref = 24.0     # reference height [m]
+erf.most.surf_temp_flux    = 0.0
+zhi.type = "SlipWall"
+zhi.theta_grad = 0.003  # adiabatic upper boundary
+erf.fixed_dt = 0.5  # time step [s]
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v              = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+erf.check_file      = chk_neutral     # root name of checkpoint file
+erf.check_int       = -1             # number of timesteps between checkpoints
+
+# DATA COLLECTION
+erf.data_log    = neutral_hist.dat neutral_profiles.dat
+erf.profile_int = 10
+
+# PLOTFILES
+erf.plot_file_1     = plt_neutral      # prefix of plotfile name
+erf.plot_int_1      = 1000             # number of timesteps between plotfiles
+erf.plot_vars_1     = x_velocity y_velocity z_velocity theta
+
+# SOLVER CHOICE
+erf.use_gravity = true
+
+erf.molec_diff_type = "None"
+# When using a PBL model, LES closure should be disabled
+erf.les_type        = "None"
+erf.Cs = 0.11
+erf.theta_ref = 300.0
+
+# INITIALIZATION - from sounding file
+erf.init_type = "input_sounding"
+erf.sounding_type = Ideal
+
+
+# PBL CONFIGURATION - Baseline MRF (neutral case)
+erf.pbl_type = "MRF"
+
+# MRF PBL PARAMETERS (for MRF model)
+erf.enable_mrf_cloud_adjustment = false
+erf.enable_mrf_countergradient = true
+erf.enable_mrf_unbounded_vpert = false
+erf.mrf_moistvars = false
+erf.pbl_mrf_Ribcr = 0.5
+erf.pbl_mrf_const_b = 7.8
+erf.pbl_mrf_coriolis_freq = 1.0e-4
+erf.pbl_mrf_highres_bounds = false
+erf.pbl_mrf_sf = 0.1
+erf.pbl_mrf_use_zero_ri_extent = false
+
+# ============================================================
+# YSUNew PBL PARAMETERS (uncomment to use YSU model)
+# For YSUNew model:
+#   1. Comment out the MRF parameters above
+#   2. Uncomment the following YSU parameters
+#   3. Change 'erf.pbl_type = "MRF"' to 'erf.pbl_type = "YSUNew"'
+# ============================================================
+# erf.enable_ysu_cloud_pblh = true
+# erf.enable_ysu_countergradient = true
+# erf.enable_ysu_entrainment = true
+# erf.enable_ysu_liquid_theta = true
+# erf.enable_ysu_sat_limiter = false
+# erf.enable_ysu_terrain_pblh_floor = true
+# erf.enable_ysu_topdown = true
+# erf.pbl_ysu_coriolis_freq = 1.0e-4
+# erf.pbl_ysu_force_over_water = false
+# erf.pbl_ysu_land_Ribcr = 0.25
+# erf.pbl_ysu_unst_Ribcr = 0.0
+# erf.pbl_ysu_use_consistent_coriolis = false
+# erf.pbl_ysunew_highres_bounds = false
+# erf.ysu_moistvars = false
+# erf.ysu_qcloud_threshold = 1.0e-4
+
+# Disable countergradient corrections (not relevant for neutral conditions)
+
+# Enable moisture diffusivity for completeness (but no moisture in neutral case)
+
+# CORIOLIS EFFECT
+erf.latitude = 45.0  # mid-latitude (f = 1.0e-4 s^-1)
+
+erf.use_coriolis = true
+
+erf.abl_driver_type = "GeostrophicWind"
+
+erf.abl_geo_wind = 15.0 0.0 0.0
+
+erf.input_sounding_file = "input_sounding_lw_isothermal"
+
+# RADIATION - Two-Stream Longwave (Isothermal mode)
+erf.radiation_type = "TwoStream"
+erf.radiation.sw_enabled = false
+erf.radiation.lw_enabled = true
+erf.radiation.tau_lw_per_layer = 1.0
+# Physically isothermal column: theta(z) = 300 exp(g z / (c_p 300)) gives T = 300 K everywhere,
+# and the surface is black at the same 300 K.
+erf.radiation.surface_temp_k = 300.0
+erf.radiation.surface_emissivity_lw = 1.0
+erf.radiation.v = 1
+erf.radiation.diag_file = "radiation_lw_diag.dat"

Exec/CanonicalTests/Radiation/LW_Isothermal/radiation_lw_diag.dat

@@ -0,0 +1,3 @@
+step,time,call_site,SW_surface,SW_TOA,SW_up_TOA,LW_net_surface,LW_up_TOA,heating_rate_max,SEB_residual_mean,SEB_residual_max,T_s_mean,T_s_max,q_s_mean,q_s_max
+0,0.000000e+00,pre_dycore,0.000000e+00,0.000000e+00,0.000000e+00,-3.986685e-03,4.593046e+02,1.804952e-02,nan,nan,nan,nan,nan,nan
+0,2.500000e-01,post_dycore,0.000000e+00,0.000000e+00,0.000000e+00,-3.986685e-03,4.593046e+02,1.804952e-02,nan,nan,nan,nan,nan,nan