Refine RSL c2m/c2h with explicit R=beta*hd/elm parameterisation (GH#1055)

[sunt05]

Summary

Implements the revised RSL equations from #1055, removing the assumption $d = \beta^2 L_c$ and using the general displacement height $d = z_H - z_d$ throughout cal_c2 / cal_c2h. This follows the pre-substitution forms from Harman & Finnigan (2008), making the formulation valid for any displacement height definition (including the revised $z_d$ from #302).

Building on @vitorlavor's testing (which confirmed dynamic $c_{2m}$ and $c_{2h}$ are stable), this PR refines the formulas to use explicit intermediate variables: mixing length $l_m = 2\beta^3 L_c$, displacement height $h_d = z_H - z_d$, and the ratio $R = \beta h_d / l_m$.

Revised equations (HF08)

Eqs. 25a,b — RSL correction at canopy top:

$$\hat{\phi}_m(0) = \frac{\kappa , d}{\phi_m , l_m}$$

$$\hat{\phi}_h(0) = \frac{\kappa , S_{cc} , d}{\phi_h , l_m}$$

Eqs. 28a,b — Profile shape parameters (dynamic, replacing fixed 0.5):

$$c_{2m} = \frac{\kappa \left( R + 1 - \frac{d}{\phi_m} \frac{d\phi_m}{dz} \right)}{\beta \phi_m - \kappa R}$$

$$c_{2h} = \frac{\kappa S_{cc} \left( f R + 1 - \frac{d}{\phi_h} \frac{d\phi_h}{dz} \right)}{\beta \phi_h - \kappa S_{cc} R}$$

Eqs. 25c,d — Correction coefficients:

$$c_m = \left(1 - \hat{\phi}_m \right) \exp\left(c_{2m} , R \right)$$

$$c_h = \left(1 - \hat{\phi}_h \right) \exp\left(c_{2h} , R \right)$$

where $R = \beta d / l_m$. Under the old assumption $d = \beta^2 L_c$, $R = 1/2$, recovering the original fixed-coefficient forms (Eqs. 25, 28 of HF08).

Additional changes

• Singularity guard: when $\hat{\phi} \ge 1$ (neutral/stable), RSL correction is set to zero
• Clamping: $c_{2m}$, $c_{2h}$ clamped to [0, 20] to prevent EXP overflow
• Dead code cleanup: ~50 commented-out lines removed from suews_phys_rslprof.f95
• Reference sample_output.csv.gz regenerated to match the refined physics

Test plan

• make test-smoke passes (13/14; Rust CLI test is a pre-existing failure unrelated to this change)
• Physical bounds verified for all 8 key variables
• SEB closes exactly at every timestep ($Q_H = Q^* + Q_F - Q_E - \Delta Q_S$, residual = 0)
• No new NaN values introduced
• Comparison figures and quantitative analysis posted as PR comment
• CI passes on Linux/Windows

Closes #1055

[github-actions]

CI Build Plan

Changed Files

Fortran source (1 file)

• src/suews/src/suews_phys_rslprof.f95

Tests (1 file)

• test/fixtures/data_test/sample_output.csv.gz

Build Configuration

	Configuration
Platforms	Linux x86_64, macOS ARM64, Windows x64
Python	3.9, 3.14
Test tier	standard (all except slow)
UMEP build	Yes (compiled extension may differ)
PR status	Ready (standard matrix)

Rationale

• Fortran source changed -> multiplatform build required
• Test files changed -> validation build
• Compiled extension ABI may differ -> UMEP (NumPy 1.x) build included

---

_{Updated by CI on each push. See path-filters.yml for category definitions.}

[sunt05]

Validation: branch vs master baseline (post-rebase)

This PR refines the RSL c2m/c2h formulas to use an explicit R = beta*hd/elm parameterisation. Since master already includes the dynamic c2m/c2h from #1116, the remaining differences here are purely from the R parameterisation refinement — and they are small.

Summary of differences

• QN: negligible (max |diff| = 0.0002 W/m²)
• QF: mean |diff| = 0.002 W/m², max = 0.06 W/m²
• QS: mean |diff| < 0.001 W/m², max = 0.22 W/m²
• QH: mean |diff| = 0.001 W/m², max = 0.20 W/m²
• QE: mean |diff| < 0.001 W/m², max = 0.10 W/m²
• T2: mean |diff| = 0.001 °C, max = 0.12 °C
• RH2: mean |diff| = 0.004%, max = 0.58%
• U10: mean |diff| = 0.008 m/s, max = 0.55 m/s

No new NaN values. All physical bounds satisfied. SEB closes exactly at every timestep (residual = 0).

How RSL changes propagate to the energy balance

RSLProfile is a post-loop diagnostic — called after the SEB iteration converges, not inside it. So c2m/c2h do not directly alter RA_h. The propagation is cross-timestep, controlled by the RSLLevel parameter:

• RSLLevel = 0: forcing temperature used directly — RSL changes would not propagate to SEB
• RSLLevel = 1 (sample config): RSL-diagnosed 2 m temperature (T2_C) replaces forcing temperature
• RSLLevel = 2: RSL-diagnosed half-building-height temperature (T_half_bldg_C) used instead

With RSLLevel = 1, two cross-timestep pathways are active:

1. SurfaceResistance (driver L3412): Tair = T2_C feeds into stomatal conductance. Changed RS -> changed QE (Penman-Monteith) -> changed QH (residual).
2. AnthropogenicEmissions (driver L1132): same T2_C feeds heating/cooling demand, changing QF.

Time series differences

Scatter comparison (new vs master)

Flux overlay near peak QH difference