update

ARPA-SIMC · Nov 16, 2023 · 18f3164 · 18f3164
1 parent 8352519
commit 18f3164
Showing 1 changed file with 22 additions and 5 deletions.
diff --git a/agrolib/soil/soil.cpp b/agrolib/soil/soil.cpp
@@ -959,7 +959,9 @@ namespace soil
      */
     bool fittingWaterRetentionCurve(Crit3DHorizon &horizon, const Crit3DFittingOptions &fittingOptions)
     {
-        if (! fittingOptions.useWaterRetentionData || horizon.dbData.waterRetention.size() == 0)
+        unsigned int nrObsValues = unsigned(horizon.dbData.waterRetention.size());
+
+        if (! fittingOptions.useWaterRetentionData || nrObsValues == 0)
         {
             // nothing to do
             return true;
@@ -971,8 +973,6 @@ namespace soil
             return false;
         }
 
-        unsigned int nrObsValues = unsigned(horizon.dbData.waterRetention.size());
-
         // search theta max
         double psiMin = 10000;      // [kpa]
         double thetaMax = 0;        // [m3 m-3]
@@ -981,6 +981,7 @@ namespace soil
             psiMin = std::min(psiMin, horizon.dbData.waterRetention[i].water_potential);
             thetaMax = std::max(thetaMax, horizon.dbData.waterRetention[i].water_content);
         }
+        // add theta sat if minimum observed value is greater than 3 kPa
         bool addThetaSat = ((thetaMax < horizon.vanGenuchten.thetaS) && (psiMin > 3));
 
         // set values
@@ -1045,8 +1046,24 @@ namespace soil
         }
         else
         {
-            pmin[2] = 0.01;
-            pmax[2] = 10;
+            // search last value with theta == thetaMax
+            double heMin = 0.01;                            // kPa
+            double heMax = 10;                              // kPa
+            for (unsigned int i = 0; i < nrObsValues; i++)
+            {
+                double delta = (thetaMax - horizon.dbData.waterRetention[i].water_content);
+                double psi = horizon.dbData.waterRetention[i].water_potential;
+                if (delta < 0.001)
+                {
+                   heMin = std::max(heMin, psi);
+                }
+                else
+                {
+                   heMax = std::min(heMax, psi);
+                }
+            }
+            pmin[2] = heMin;
+            pmax[2] = heMax;
         }
 
         // Van Genuchten alpha parameter [kPa^-1]