Better report any convergence errors in harder precip with newton solver

Author: Chrismarsh

src/modules/Harder_precip_phase.cpp

@@ -78,7 +78,7 @@ void Harder_precip_phase::run(mesh_elem& face)
     double Ta = (*face)["t"_s]+273.15; //K
     double T =  (*face)["t"_s];
     double RH = (*face)["rh"_s];
-    double ea = RH/100 * 0.611*exp( (17.3*T) / (237.3+T));
+    double ea = RH/100.0 * 0.611*exp( (17.3*T) / (237.3+T));
 
     // (A.6)
     double D = 2.06 * pow(10,-5) * pow(Ta/273.15,1.75);
@@ -117,7 +117,29 @@ void Harder_precip_phase::run(mesh_elem& face)
     double max = 50;
     double digits = 6;
 
-    double Ti = boost::math::tools::newton_raphson_iterate(fx, guess, min, max, digits);
+    double Ti = 0;
+
+    try
+    {
+       Ti = boost::math::tools::newton_raphson_iterate(fx, guess, min, max, digits);
+    }
+    catch(...)
+    {
+        SPDLOG_ERROR("Ta={}, RH={}, ea={}, guess={}, min={}, max={}",
+            Ta, RH, ea, guess, min, max);
+        try
+        {
+            // try again but instrument it this time to see what went wrong
+#define BOOST_MATH_INSTRUMENT
+            Ti = boost::math::tools::newton_raphson_iterate(fx, guess, min, max, digits);
+#undef BOOST_MATH_INSTRUMENT
+        }
+        catch(...)
+        {
+        }
+        CHM_THROW_EXCEPTION(module_error, "Harder_precip_phase newton_raphson_iterate failed to converge");
+    }
+
 
     double frTi = 1.0 / (1.0+b*pow(c,Ti));