Fix floating invalid exceptions in Process_Library.F90 Bergeron code

[mathomp4]

Important

Status: Proposed Solution — Needs Testing
This pull request contains a candidate fix for floating invalid exceptions encountered during debug runs. It requires validation through testing (e.g. running the 5-day c90 L181 debug simulation with ifort) to confirm stability and verify physical output consistency.

---

Context & Reported Issue

During a c90 L181 run compiled with ifort debug flags, the simulation failed after approximately 5 days of model time with the following traceback:

forrtl: error (65): floating invalid
Image              PC                Routine            Line        Source
libpthread-2.31.s  0000147201754910  Unknown               Unknown  Unknown
GEOSgcm.x          0000000006F67B6F  Unknown               Unknown  Unknown
libGEOSmoist_Grid  000014722B09B3BA  geosmoist_process        3053  Process_Library.F90
libGEOSmoist_Grid  000014722B0999EC  geosmoist_process        2768  Process_Library.F90
libGEOSmoist_Grid  000014722BC71CF1  geos_gfdl_1m_inte         808  GEOS_GFDL_1M_InterfaceMod.F90

---

Detailed Root Cause Analysis

1. Floating Invalid Exception on Real Exponentiation (** 1.94):
   Line 3053 in Process_Library.F90 inside Bergeron_Partition computes:

   diff = (0.211 * 1013.25 / (pl + 0.1)) * (((t_env + 0.1) / MAPL_TICE)**1.94) * 1e-4

   When compiled with Intel Fortran (ifort) debug flags (specifically -fpe0 or floating-point trap options enabled), real-number exponentiation $x^y$ with a real non-integer exponent ($y = 1.94$) is evaluated via $\exp(1.94 \times \ln(x))$. If $x \le 0$, $\ln(x)$ is mathematically undefined for real floating-point arithmetic, immediately trapping with forrtl: error (65): floating invalid.

   If t_env <= -0.1 K (or if t_env becomes negative during solver iterations), (t_env + 0.1) is non-positive, triggering the crash. Note that while (t_env + 0.1) was likely intended as a non-zero guard, + 0.1 does not prevent the base from being negative if t_env <= -0.1.

2. Solver Overshoot in hystpdf:
   In hystpdf (around lines 2808–2818), t_env is updated iteratively inside the do n = 1, nmax loop using a secant method solver:

   t_env = t_old - (f_t_env - t_old) * (t_old - t_env_old_p) / denom

   If denom becomes very small, the secant update step can overshoot dramatically, driving t_env to unphysical negative values (e.g. t_env < 0 K). On the subsequent iteration ($n+1$), hystpdf passes this negative t_env to Bergeron_Partition (line 2768), triggering the exception at line 3053.

3. Secondary Floating Invalid Vulnerabilities:

   • Line 3060 (d_crystal):

     d_crystal = max((q_ice_ls / (n_ice_active * den_ice * MAPL_PI))**(0.333), 20.0e-6)

     Real exponentiation **(0.333) will also trap floating invalid if the term inside parentheses is non-positive.
   • Lines 3023–3025 (f_mass_ice & n_ice_active):
     If liquid mass q_tot_liq becomes slightly negative due to numerical truncation, f_mass_ice = q_tot_ice / q_tot_mass can exceed 1.0, producing negative active ice crystal concentrations (n_ice_active < 0).
   • Line 3054 (den_air):
     den_air = (pl * 100.0) / (MAPL_RGAS * t_env) could divide by zero or yield negative air density if t_env <= 0.

---

Code Changes Summary

• hystpdf (line 2825): Clamped t_env to a physically plausible range [100.0, 400.0] K after secant solver updates to prevent iteration overshoot from propagating negative temperature values into subsequent calls.
• Bergeron_Partition (line 3053 & 3054):
  • Guarded t_env in diffusivity calculation: diff = (0.211 * 1013.25 / max(pl + 0.1, 0.1)) * (((max(t_env, 1.0) + 0.1) / MAPL_TICE)**1.94) * 1e-4
  • Guarded t_env in air density calculation: den_air = (pl * 100.0) / (MAPL_RGAS * max(t_env, 1.0))
• Bergeron_Partition (lines 3024–3025):
  • Clamped f_mass_ice to [0.0, 1.0] and n_ice_active to >= 0.0.
• Bergeron_Partition (line 3060):
  • Ensured den_ice > 0.0 and guarded the base of **(0.333) with max(..., 0.0).

---

Verification & Testing Request

Note

This PR is opened as a Draft for review and testing.

• Recommended Test: Re-run the c90 L181 debug simulation with ifort past the 5-day mark to confirm that the floating invalid error no longer occurs.
• Sanity Checks: Verify that physical output diagnostics remain bit-for-bit identical under normal (non-debug/non-overshoot) conditions.