trop_cyclone: truncate hol_x to prevent unphysical results of H10 model

climada/hazard/test/test_trop_cyclone.py

@@ -384,21 +384,21 @@ def test_holland_2010_pass(self):
         hol_x = _x_holland_2010(si_track, d_centr, close_centr)
         np.testing.assert_array_almost_equal(hol_x, [
             [0.5, 0.500000, 0.485077, 0.476844, 0.457291],
-            [0.5, 0.500000, 0.410997, 0.289203, 0.000000],
-            [0.5, 0.497620, 0.131072, 0.000000, 0.000000],
+            [0.5, 0.500000, 0.410997, 0.400000, 0.000000],
+            [0.5, 0.497620, 0.400000, 0.400000, 0.400000],
             [0.5, 0.505022, 1.403952, 1.554611, 2.317948],
         ])
 
         v_ang_norm = _stat_holland_2010(si_track, d_centr, close_centr, hol_x)
-        np.testing.assert_array_almost_equal(v_ang_norm, [
+        np.testing.assert_allclose(v_ang_norm, [
             # first column: converge to 0 when approaching storm eye
             # second column: vmax at RMW
             # fourth column: peripheral speed (17 m/s) at peripheral radius (unless x is clipped!)
             [0.0000000, 35.000000, 21.181497, 17.00000, 12.103461],
-            [1.2964800, 34.990037, 21.593755, 17.00000, 0.0000000],
-            [0.3219518, 15.997500, 13.585498, 16.00000, 16.000000],
+            [1.2964800, 34.990037, 21.593755, 12.89131, 0.0000000],
+            [0.3219518, 15.997500, 9.7120060, 8.087240, 6.2289690],
             [24.823469, 89.992938, 24.381965, 17.00000, 1.9292020],
-        ])
+        ], atol=1e-5)
 
     def test_stat_holland_1980(self):
         """Test _stat_holland_1980 function. Compare to MATLAB reference."""

climada/hazard/trop_cyclone.py

@@ -1741,18 +1741,22 @@ def _x_holland_2010(
     r_max_norm = (r_max / r_n)**hol_b
     if vmax_in_brackets:
         x_n = np.log(v_n) / np.log(v_max_s**2 * r_max_norm * np.exp(1 - r_max_norm))
+
+        # Truncate to prevent maximum to be shifted too far away from RMW. The value 1.0 has been
+        # found to be reasonable by manual testing of thresholds. Note that this means that the
+        # peripheral wind speed v_n is not exactly attained in some cases.
+        x_n = np.fmin(x_n, 1.0)
     else:
         x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm))
 
     # linearly interpolate between max exponent and peripheral exponent
     x_max = 0.5
     hol_x[close_centr] = x_max + np.fmax(0, d_centr - r_max) * (x_n - x_max) / (r_n - r_max)
 
-    # Negative hol_x values appear when v_max_s is very close to or even lower than v_n (which
-    # should never happen in theory). In those cases, wind speeds might decrease outside of the eye
-    # wall and increase again towards the peripheral radius (which is actually unphysical).
-    # We clip hol_x to 0, otherwise wind speeds keep increasing indefinitely away from the eye:
-    hol_x[close_centr] = np.fmax(hol_x[close_centr], 0.0)
+    # Truncate to prevent wind speed from increasing again towards the peripheral radius (which is
+    # unphysical). A value of 0.4 has been found to be reasonable by manual testing of thresholds.
+    # Note that this means that the peripheral wind speed v_n is not exactly attained sometimes.
+    hol_x[close_centr] = np.fmax(hol_x[close_centr], 0.4)
 
     return hol_x