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

[hgopalan]

Summary

Adds a two-stream atmospheric radiation model as an alternative to RRTMGP, together with an optional simplified surface energy balance (SEB) module that evolves surface temperature and moisture by a force-restore method where Noah-MP is not driving them. The model is selected with erf.radiation_type = "TwoStream" and is independent of the RRTMGP path, which is still selected with erf.radiation_model.

Everything new is off by default, so a run that does not ask for two-stream radiation is unaffected.

Shortwave

• Direct beam by Beer-Lambert attenuation, F_dir(z) = S0 cos(zenith) exp(-tau_cum(z) / cos(zenith)).
• Diffuse field from the Meador-Weaver (1980) two-stream solution with Zdunkowski PIFM coefficients, parameterized by single-scattering albedo and asymmetry factor, with independent clear-sky and cloud values.
• The layer solutions are combined by the adding method, so the surface albedo reflects the beam back into the column rather than being applied only as an absorbed fraction. Direct and diffuse albedo can differ (surface_albedo_sw_diffuse).
• Optional Earth-Sun distance factor (Spencer 1971) scaling S0 by day of year.

Longwave

• Gray-gas two-stream sweep on layer interfaces (Toon et al. 1989), with Stefan-Boltzmann surface emission at a configurable emissivity and reflection of the downwelling flux at the surface.
• Cooling to space is set by the flux divergence, with no isothermal override.

Optical depth

Two models, selected by erf.radiation.tau_model:

• per_layer (default): a fixed optical depth per layer, so the column value scales with the number of vertical cells.
• mass: built from each layer's mass path rho dz, so the column optical depth is set by the atmosphere rather than by the grid. Shortwave mixes dry absorption, Rayleigh scattering, water vapour and cloud water with extinction-weighted single-scattering albedo and asymmetry factor. Longwave uses the corresponding mass absorption coefficients, with the cloud term expressed as the Stephens (1978) emissivity.

On top of either model: a prescribed cloud band, moisture- and cloud-aware dynamic optical depth diagnosed from qv and qc, prognostic cloud fraction diagnosed from relative humidity and cloud water, and a prescribed bulk aerosol optical depth that is constant or exponential with height.

Coupling

• Heating rates are stored as potential-temperature tendencies, dtheta/dt = (dT/dt) / pi, matching the convention the RhoTheta source term and the RRTMGP path already use. They are written to qheating_rates and exposed as the qsrc_sw and qsrc_lw plot variables.
• Applied once per slow step, during slow-RHS construction only, so the tendency is consistent across all fast substeps.
• Time-varying solar zenith angle from simulation time, latitude, longitude and day of year, with a static fallback.
• Per-column surface albedo, emissivity and temperature from an LSM when present, falling back to scalars.
• Non-uniform vertical grids through the physical cell-centre heights.
• Optional smoothing and limiting of the radiative tendency handed to YSUNew.

Simplified surface energy balance

• Allocates the surface flux fields, passing Noah-MP values through when it is active and using scalar fallbacks otherwise.
• Diagnostic mode reports the residual R_net - H - LE - G.
• Prognostic mode evolves surface temperature and moisture by force-restore with bounded explicit Euler updates, and switches itself off at any level where Noah-MP is already driving the surface state.

Grid requirement

The column sweep integrates a whole atmospheric column in one pass, so each grid must span the domain in the vertical. ERF only decomposes in z when amr.max_grid_size_z is smaller than the number of cells in z; setting it at least that large is sufficient. A vertically decomposed grid aborts with a message naming the input rather than returning heating rates that look plausible and are wrong. Results do not depend on the horizontal decomposition or on the MFIter tile size.

Input validation

Several misconfigurations used to be accepted in silence. None of them aborts, so existing inputs keep running, but each now prints a warning naming the input, the valid values, and what the run will actually do instead:

• an unrecognized tau_model, tau_profile_type, aerosol_profile_type or diag_callsite_mode falls back to its default, which is stated in the warning;
• aerosol_profile_type = table is not implemented and contributes no aerosol optical depth;
• seb_prognostic_enable reports that it also switches on seb_enable and seb_diagnostic_enable, which it depends on;
• erf.radiation_type = RRTMGP is not wired to a solver, so no radiation is computed. RRTMGP is selected with erf.radiation_model.

Cost and scaling

Exec/CanonicalTests/Radiation/TwoStream_vs_RRTMGP_Timing/ measures what each solver costs per radiation call and how that cost scales, with everything except the solver held fixed.

Measured on one CPU core, marginal cost per call:

Grid	Cells	Two-stream	RRTMGP
42 x 4 x 42	7,056	1.08 ms	433 ms
84 x 4 x 42	14,112	2.07 ms	881 ms
84 x 4 x 84	28,224	4.27 ms	1,812 ms
168 x 4 x 84	56,448	8.55 ms	3,581 ms

The two solvers compute different things: RRTMGP resolves 224 shortwave and 256 longwave g-points with gas optics from lookup tables, while the two-stream solver does a single gray sweep per column. The difference is the cost of spectral detail, so this is a guide to choosing a solver for a given run rather than a ranking. Cost per cell is flat for both across a factor of eight in problem size, which means either one extrapolates predictably to a production grid.

Both radiation paths previously had no BL_PROFILE region, so radiation cost could not be separated from the dycore for either solver. One is now added at each entry point, and the measurement reads exclusive time from AMReX's TinyProfiler.

Keeping it fair: both inputs pull in one shared settings block through AMReX's FILE include, pinning the grid, timestep, sounding, moisture model and surface temperature; radiation is called every slow step on both sides; plotfiles and both diagnostic logs are off so file I/O stays outside the measured region. The reported figure is the marginal cost across two step counts, which cancels one-time work, notably the roughly 150 ms RRTMGP spends reading its lookup tables on the first call. check_timing_fairness.py enforces those invariants and needs no RRTMGP build.

Caveats: a single CPU core, RRTMGP targets GPUs, and cost per call is not cost per unit of accuracy.

Testing

• 34 GoogleTest cases covering the shortwave and longwave kernels and the column sweep, including energy-budget closure, the Exner convention, surface reflection, and resolution independence of the mass model.
• Two ctest regressions: a column-heating check that verifies the vertical structure and orientation of the heating rates, and a plotfile check that qsrc_sw and qsrc_lw are written.
• 17 canonical cases under Exec/CanonicalTests/Radiation/, each with a Python checker, covering clear-sky and cloudy shortwave, scattering, longwave, non-uniform spacing, dynamic optical depth, prognostic cloud fraction, aerosol, diurnal solar geometry, surface heterogeneity, the mass optical-depth model, RhoTheta and YSUNew coupling, and all three SEB modes.

Documentation

• Docs/sphinx_doc/theory/Radiation.rst covers the shortwave and longwave formulations, the optical-depth models, the potential-temperature convention, the SEB equations and the grid requirement.
• Docs/sphinx_doc/Inputs.rst documents every input.
• Exec/CanonicalTests/Radiation/inputs_twostream_master is the reference input file and lists all 79 erf.radiation.* parameters with defaults, units and descriptions.