Merge commit '4611b39429af2acffc3b3153814228afd1baba37'

ARPA-SIMC · Dec 19, 2024 · dd97ced · dd97ced
2 parents c9b8782 + 4611b39
commit dd97ced
Show file tree

Hide file tree

Showing 4 changed files with 32 additions and 25 deletions.
diff --git a/agrolib/soil/soil.cpp b/agrolib/soil/soil.cpp
@@ -152,6 +152,7 @@ namespace soil
 
         this->fieldCapacity = NODATA;
         this->wiltingPoint = NODATA;
+        this->waterContentSAT = NODATA;
         this->waterContentFC = NODATA;
         this->waterContentWP = NODATA;
 
@@ -223,19 +224,17 @@ namespace soil
 
         horizonPtr = horizonPointer;
 
-        double hygroscopicHumidity = -2000;     // [kPa]
-        double waterContentHH = soil::thetaFromSignPsi(hygroscopicHumidity, *horizonPtr);
-
-        // [-]
-        soilFraction = horizonPtr->getSoilFraction();
-
         // [mm]
-        SAT = horizonPtr->vanGenuchten.thetaS * soilFraction * thickness * 1000;
-        FC = horizonPtr->waterContentFC * soilFraction * thickness * 1000;
-        WP = horizonPtr->waterContentWP * soilFraction * thickness * 1000;
-        HH = waterContentHH * soilFraction * thickness * 1000;
+        SAT = horizonPtr->waterContentSAT * thickness * 1000;
+        FC = horizonPtr->waterContentFC * thickness * 1000;
+        WP = horizonPtr->waterContentWP * thickness * 1000;
         critical = FC;
 
+        // hygroscopic humidity
+        double hygroscopicHumidity = -2000;     // [kPa]
+        double volWaterContentHH = soil::thetaFromSignPsi(hygroscopicHumidity, *horizonPtr);    // [m3 m-3]
+        HH = volWaterContentHH * horizonPtr->getSoilFraction() * thickness * 1000;
+
         return true;
     }
 
@@ -413,8 +412,10 @@ namespace soil
     // estimate bulk density from total porosity
     double estimateBulkDensity(const Crit3DHorizon &horizon, double totalPorosity, bool increaseWithDepth)
     {
-        if (int(totalPorosity) == int(NODATA))
+        if (isEqual(totalPorosity, NODATA))
+        {
             totalPorosity = (horizon.vanGenuchten.refThetaS);
+        }
 
         double specificDensity = estimateSpecificDensity(horizon.organicMatter);
         double refBulkDensity = (1 - totalPorosity) * specificDensity;
@@ -630,7 +631,7 @@ namespace soil
      * \brief Compute volumetric water content from signed water potential
      * \param signPsi water potential       [kPa]
      * \param horizon
-     * \return volumetric water content     [m^3 m-3]
+     * \return volumetric water content     [m3 m-3]
      */
     double thetaFromSignPsi(double signPsi, const Crit3DHorizon &horizon)
     {     
@@ -710,14 +711,14 @@ namespace soil
 
 
     /*!
-     * \brief return current volumetric water content [m3 m^3]
+     * \brief return current volumetric water content [-]
      */
     double Crit1DLayer::getVolumetricWaterContent()
     {
-        // waterContent [mm]
-        // thickness [m]
-        double theta = waterContent / (thickness * soilFraction * 1000);
-        return theta;
+        // thickness [m] -> mm
+        double soilThickness = thickness * soilFraction * 1000.;
+
+        return waterContent / soilThickness;
     }
 
 
@@ -727,7 +728,7 @@ namespace soil
     double Crit1DLayer::getDegreeOfSaturation()
     {
         double theta = getVolumetricWaterContent();
-        return (theta - horizonPtr->vanGenuchten.thetaR) / (horizonPtr->vanGenuchten.thetaS - horizonPtr->vanGenuchten.thetaR);
+        return SeFromTheta(theta, *horizonPtr);
     }
 
 
@@ -963,8 +964,10 @@ namespace soil
 
         horizon.fieldCapacity = soil::getFieldCapacity(horizon.texture.clay, soil::KPA);
         horizon.wiltingPoint = soil::getWiltingPoint(soil::KPA);
-        horizon.waterContentFC = soil::thetaFromSignPsi(horizon.fieldCapacity, horizon);
-        horizon.waterContentWP = soil::thetaFromSignPsi(horizon.wiltingPoint, horizon);
+
+        horizon.waterContentSAT = horizon.vanGenuchten.thetaS * horizon.getSoilFraction();
+        horizon.waterContentFC = soil::thetaFromSignPsi(horizon.fieldCapacity, horizon) * horizon.getSoilFraction();
+        horizon.waterContentWP = soil::thetaFromSignPsi(horizon.wiltingPoint, horizon) * horizon.getSoilFraction();
 
         return true;
     }
@@ -999,8 +1002,9 @@ 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));
+
+        // add theta sat if minimum observed value is greater than 5 kPa
+        bool addThetaSat = ((thetaMax < horizon.vanGenuchten.thetaS) && (psiMin > 5));
 
         // set values
         unsigned int nrValues = nrObsValues;

diff --git a/agrolib/soil/soil.h b/agrolib/soil/soil.h
@@ -140,8 +140,11 @@
 
             double fieldCapacity;               /*!<  [kPa]     */
             double wiltingPoint;                /*!<  [kPa]     */
+
+            double waterContentSAT;             /*!<  [m3 m-3]*/
             double waterContentFC;              /*!<  [m3 m-3]*/
             double waterContentWP;              /*!<  [m3 m-3]*/
+
             double PH;                          /*!<  [-]       */
             double CEC;                         /*!<  [meq/100g]*/
 

diff --git a/agrolib/soilFluxes3D/header/types.h b/agrolib/soilFluxes3D/header/types.h
@@ -77,7 +77,7 @@
 
         double Roughness;           /*!< [s m-1/3] surface: Manning roughness */
 
-        //for heat
+        // for heat flux
         double organicMatter;       /*!< [-] fraction of organic matter */
         double clay;                /*!< [-] fraction of clay */
     };

diff --git a/agrolib/soilWidget/soilWidget.cpp b/agrolib/soilWidget/soilWidget.cpp
@@ -754,13 +754,13 @@ void Crit3DSoilWidget::setInfoTextural(int nHorizon)
     }
     else
     {
-        if (mySoil.horizon[unsigned(nHorizon)].vanGenuchten.thetaS == NODATA)
+        if (mySoil.horizon[unsigned(nHorizon)].waterContentSAT == NODATA)
         {
             satValue->setText(QString::number(NODATA));
         }
         else
         {
-            satValue->setText(QString::number(mySoil.horizon[unsigned(nHorizon)].vanGenuchten.thetaS, 'f', 3));
+            satValue->setText(QString::number(mySoil.horizon[unsigned(nHorizon)].waterContentSAT, 'f', 3));
         }
 
         if (mySoil.horizon[unsigned(nHorizon)].waterContentFC == NODATA)