SLUCM surface-temperature solver (TS_SCHEME=1) can diverge to infinity, corrupting output before it crashes

[kylegilchrist-pixel]

Version: v4.8.0 (commit 06d4240) Compiler: GNU Fortran 15.2.1, -O3 -march=znver3 -ftree-vectorize -funroll-loops -fno-trapping-math File: phys/module_sf_urban.F Configuration: sf_urban_physics=1, sf_surface_physics=4 (Noah-MP), bl_pbl_physics=5, sf_sfclay_physics=5, use_wudapt_lcz=1, TS_SCHEME: 1, GROPTION: 0, 3 domains (9/3/1 km), ERA5 forcing

Summary
The Newton-Raphson solver for urban wall and road temperatures applies its correction without any limit on step size or on the resulting iterate. When the Jacobian approaches singularity the iterate runs away; once it reaches infinity the saturation-vapour-pressure expression evaluates inf/inf and injects a NaN into the surface fluxes.

The NaN then leaves SLUCM and, in our configuration, eventually reaches MYNN-SFC where it causes a segmentation fault (filed separately as an unguarded lookup index). The more serious consequence is that the solver writes physically impossible surface temperatures for some time before any crash occurs, so a run that completes is not necessarily correct.

Mechanism
module_sf_urban.F:1534 opens the wall/road iteration:

fortran IF (TS_SCHEME == 1) THEN ! TB, TG Solving Non-Linear Simultaneous Equation by Newton-Rapson DO ITERATION=1,20

and line 1640 computes the correction:

TB = TBP + DTB TG = TGP + DTG ```

(FX*GY-GX*FY) is the Jacobian determinant and can approach zero. Nothing bounds DTB/DTG, and nothing bounds TB/TG. The loop does have a convergence test at line 1665, so a converging cell exits cleanly — the gap is purely on the divergent path.

Once TGP is infinite, line 1541 evaluates

fortran ES=6.11*EXP( (2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )

(inf - 273.15)/(273.15*inf) is inf/inf, i.e. an invalid operation producing NaN. A build with -ffpe-trap=invalid traps precisely here:

3 __module_sf_urban_MOD_urban
at phys/module_sf_urban.f90:1534
4 __noahmpurbandrivermainmod_MOD_noahmp_urban
at noahmp/drivers/wrf/NoahmpUrbanDriverMainMod.F90:541
5 __module_surface_driver_MOD_surface_driver._omp_fn.8

The roof solver at line 1296 (DTR = F/DFDT) and the green-roof loop at line 1349 (DO ITERATION=1,100) share the same unbounded structure.

Evidence that output is corrupted before the crash
Urban-cell skin temperature on the innermost domain, unmodified 4.8.0 versus the same run with a bounded step:

model time	stock	with damped step
00:00	290.09 – 295.12 K	290.09 – 295.12 K (bit-identical)
01:00	290.76 – 295.48 K	290.74 – 295.13 K
02:00	247.06 – 294.41 K	290.81 – 294.59 K
02:33:55	SIGSEGV	run continues
247 K is −26 °C on an urban surface in subtropical Australia in late autumn. The solver had been producing that for at least half an hour of model time before the segfault, so the crash is a late and optional symptom. Configurations that survive may still contain the corruption.

Suggested fix
A damped Newton step with a physical bound on the iterate. Applied to both the wall/road solver and the roof solver:

```fortran IF (.not. (DTB == DTB)) DTB = 0.0 IF (.not. (DTG == DTG)) DTG = 0.0 DTB = MAX(-DT_MAX, MIN(DT_MAX, DTB)) DTG = MAX(-DT_MAX, MIN(DT_MAX, DTG))

TB = TBP + DTB TG = TGP + DTG

TB = MAX(T_MIN, MIN(T_MAX, TB)) TG = MAX(T_MIN, MIN(T_MAX, TG)) ```

We used DT_MAX = 5.0 K, T_MIN = 200 K, T_MAX = 350 K.

Two observations from testing this locally:

The first output frame was bit-identical across all 5184 cells to the unmodified build, and frame 1 differed by only 0.02–0.35 K. A converging cell takes steps far below 5 K, so the limit is inert on healthy cells.
Zero cells reached the 200 K or 350 K bounds. The hard clamp never fired; the step damping alone was sufficient. That suggests damping is the substantive part and the temperature bounds are only a backstop.
With the guard applied, the previously crashing case completed a 4-hour run and a 24-hour run with no NaN in any output field.

Caveat on reproducibility
This is a marginal instability and we want to be straightforward about that. Whether a given cell diverges on a given timestep was sensitive to optimization flags, OpenMP thread count, and even to adding a few lines of diagnostic printing — one instrumented build ran clean through the same period where the uninstrumented build crashed.

We therefore cannot offer a small portable reproducer, and we would not expect the exact failure time to reproduce on different hardware or compiler settings. What we can state confidently is the structural defect: the correction is unbounded, the iterate is unbounded, and there is a code path from that to inf/inf at line 1541. The fix is defensible on general numerical grounds regardless of whether any particular case triggers it.

Happy to supply the FPE-trap backtraces, the full run matrix, or the input configuration if useful.

[weiwangncar]

@kylegilchrist-pixel Again thank you very much for your report. We will have the experts take a look.
@cenlinhe If you know the contributors to this code, please invite them to take a look at this report. Thanks!

[cenlinhe]

This issue and fix look reasonable to me. The question would be how we should set a reasonable bound.
@tslin2 @chenghaow Could you please also take a look and let us know what you think? Thanks!

[tslin2]

@kylegilchrist-pixel Thanks for sharing this issue, it seems that you also looked at the green-roof loop, does that also have the same issue?
I had a test, which seem working by setting DTG/6 and DTB/6 (default used for green roof)
https://forum.mmm.ucar.edu/threads/resolved-wrf-single-layer-ucm-crashed.13453/
agree with @cenlinhe on how to set a reasonable bound is not easy as these are building wall and road temperatures

Also invite @DanLi-BU, @epn09 and @jyang104 to take a look as they were partially involved in developing this model

[tslin2]

juse see that new urban-nbs in WRF4.8 include some bounds, if this part is added for urban tree, then it should also be included for no tree cases (original codes):

WRF/phys/module_sf_urban.F

Lines 2074 to 2146 in 06d4240

	! Near-singular check: regularize diagonals and relax the update
	IF (ABS(det) < det_tol) THEN
	lambda = 1.0D-6 * det_scale
	A_use = A + lambda
	FF_use = FF + lambda
	KK_use = KK + lambda
	det = A_use*(FF_use*KK_use-GG*J)-B*(E*KK_use-GG*I)+C*(E*J-FF_use*I)
	relax = 0.5D0
	IF (ABS(det) < det_tol) THEN
	det = SIGN(det_tol, det)
	END IF
	WRITE_MESSAGE("WARNING: Near-singular matrix in Urban Tree Energy Budget; regularizing solve and relaxing update")
	END IF
	CHARA = 1.0D0 / det
	DTB = -CHARA* (WF*(FF_use*KK_use-GG*J)+GF*(C*J-B*KK_use)+VGF*(B*GG-C*FF_use))
	DTG = -CHARA* (WF*(GG*I-E*KK_use)+GF*(A_use*KK_use-C*I)+VGF*(C*E-A_use*GG))
	DTVG= -CHARA* (WF*(E*J-FF_use*I)+GF*(B*I-A_use*J)+VGF*(A_use*FF_use-B*E))
	DTB = DTB * relax
	DTG = DTG * relax
	DTVG = DTVG * relax
	DTB = SIGN(MIN(ABS(DTB), max_step), DTB)
	DTG = SIGN(MIN(ABS(DTG), max_step), DTG)
	DTVG = SIGN(MIN(ABS(DTVG), max_step), DTVG)
	TB = TBP + DTB
	TG = TGP + DTG
	TVG= TVGP+DTVG
	TGE = (1 - FVG) * TG + FVG * TVG
	TBP = TB
	TGP = TG
	TVGP= TVG
	TGEP=TGE
	! After updating the temperatures (e.g., TBP = TB) (urban-nbs)
	TBP = MAX(223.15, MIN(353.15, TBP)) !-50C to +80C
	TGP = MAX(223.15, MIN(353.15, TGP))
	TVGP = MAX(223.15, MIN(353.15, TVGP))
	!IF (distributed_aerodynamics_option) THEN
	!DTC = 0.0
	!ELSE
	TC1 = RW * ALPHAC + RW * FVG * ALPHAVG + RW * (1 - FVG) * ALPHAG + W * ALPHAB
	TC2 = RW * ALPHAC * TA + RW * FVG * ALPHAVG * TVGP + RW * (1 - FVG) * ALPHAG * TGP + W * ALPHAB * TBP + 2 * RTREE * LAI_TREE * (HT/100)
	TC = TC2 / TC1
	TC = MAX(223.15, MIN(353.15, TC))
	TCP = TC ! update TCP (urban-nbs)
	! Chenghao: HT unit is cgs, keep it as is [RHO - g/cm3, CP - cal/g/K]
	QC1 = W * ALPHAB * BETB + RW * (1 - FVG) * ALPHAG * BETG + RW * ALPHAC
	QC2 = W * ALPHAB * BETB * QS0B + RW * (1 - FVG) * ALPHAG * BETG * QS0G + RW * FVG * ETAG *CP / 1000 + RW * ALPHAC * QA + (2 * RTREE * LAI_TREE * CP * (ELET/697.7/60.)/100/EL)
	QC=QC2/QC1
	! Chenghao: ELET - [W/m2], convert it to cgs to be consistent with all other terms (e.g., ALPHAB)
	ES=6.11*EXP( (2.5*10**6/461.51)*(TC-273.15)/(273.15*TC) )
	DESDT=(2.5*10**6/461.51)*ES/(TC**2)
	QS0C=0.622*ES/(PS-0.378*ES)
	QC = MIN(QC, QS0C)
	QC = MAX(QC, 0.0)
	IF (QC <= 1.E-9 .AND. QCP > 1.E-9) THEN
	QC = QCP
	ELSE IF (QC <= 1.E-9 .AND. QCP <= 1.E-9) THEN
	WARN_MSG = 'WARNING: Canyon air has become completely dry.'
	WRITE_MESSAGE(WARN_MSG)
	END IF

[kylegilchrist-pixel]

@tslin2 Thanks — I checked the green-roof loop directly.

It's not fully unbounded like the wall/road case: DTGR = FV/DFDVT/6 (line ~1349) already has the /6 damping — it looks like this is the same factor from the forum thread you linked, and it's already shipped in stock WRF, not something new to add. So the green-roof loop is partially protected.

What it's still missing is a hard temperature clamp after the update — the tree branch (lines 2096–2146) clamps TBP/TGP/TVGP to 223.15–353.15 K in addition to damping the step, but the green-roof loop has no equivalent clamp on TGRP. So a persistently near-singular case there could still drift toward an unphysical value, just ~6× slower, with nothing to backstop it if it does.

The plain roof solver (DTR = F/DFDT, line ~1296) has neither damping nor a clamp — no /6, no bound. Of the three, it's the closest structural match to the original wall/road bug and the most exposed.

Loop Step damping Temp clamp

Roof (1296) None None
Green roof (1349) Implicit /6 (existing) None
Wall/road, no-tree (1640, this issue) None None
Wall/road, tree (2096) SIGN(MIN(ABS(DTB), max_step), DTB) 223.15-353.15 K

Given the tree branch is the only one of the four with full protection, it might be cleanest to treat this as one fix — backport its damping + clamp pattern to the other three — rather than three separate patches with different bounds. Happy to help test if that's useful.

[tslin2]

@kylegilchrist-pixel Thanks for providing the insightful thoughts; I agree with these ideas.

[cenlinhe]

The proposed consistent fix sounds good to me.