add degree of saturation output

Author: ftomei

criteriaModel/criteria1DCase.cpp

@@ -361,7 +361,7 @@ bool Crit1DCase::computeNumericalFluxes(const Crit3DDate &myDate, std::string &e
 }
 
 
-void Crit1DCase::saveWaterContent()
+void Crit1DCase::storeWaterContent()
 {
     prevWaterContent.clear();
     prevWaterContent.resize(soilLayers.size());
@@ -537,7 +537,7 @@ bool Crit1DCase::computeDailyModel(Crit3DDate &myDate, std::string &error)
 
     output.dailyIrrigation = 0;
     // Water fluxes (first computation)
-    saveWaterContent();
+    storeWaterContent();
     if (! computeWaterFluxes(myDate, error)) return false;
     // Irrigation
     double irrigation = checkIrrigationDemand(doy, prec, precTomorrow, output.dailyMaxTranspiration);
@@ -653,11 +653,11 @@ double Crit1DCase::getTotalWaterContent()
 
 
 /*!
- * \brief getWaterContent
+ * \brief getVolumetricWaterContent
  * \param computationDepth = computation soil depth  [cm]
- * \return volumetric water content at specific depth [-]
+ * \return volumetric water content (soil moisture) at specific depth [m3 m-3]
  */
-double Crit1DCase::getWaterContent(double computationDepth)
+double Crit1DCase::getVolumetricWaterContent(double computationDepth)
 {
     computationDepth /= 100;        // [cm] --> [m]
 
@@ -679,6 +679,34 @@ double Crit1DCase::getWaterContent(double computationDepth)
 }
 
 
+/*!
+ * \brief getDegreeOfSaturation
+ * \param computationDepth = computation soil depth  [cm]
+ * \return degree of saturation at specific depth [-]
+ */
+double Crit1DCase::getDegreeOfSaturation(double computationDepth)
+{
+    computationDepth /= 100;        // [cm] --> [m]
+    if (computationDepth <= 0 || computationDepth > mySoil.totalDepth)
+    {
+        return NODATA;
+    }
+
+    double upperDepth, lowerDepth;
+    for (unsigned int i = 1; i < soilLayers.size(); i++)
+    {
+        upperDepth = soilLayers[i].depth - soilLayers[i].thickness * 0.5;
+        lowerDepth = soilLayers[i].depth + soilLayers[i].thickness * 0.5;
+        if (computationDepth >= upperDepth && computationDepth <= lowerDepth)
+        {
+            return soilLayers[i].waterContent / soilLayers[i].SAT;
+        }
+    }
+
+    return NODATA;
+}
+
+
 /*!
  * \brief getWaterPotential
  * \param computationDepth = computation soil depth  [cm]
@@ -778,7 +806,7 @@ double Crit1DCase::getWaterDeficitSum(double computationDepth)
 
 /*!
  * \brief getWaterCapacity
- * \param computationDepth = computation soil depth  [cm]
+ * \param computationDepth = computation soil depth [cm]
  * \return sum of available water capacity (FC-WP) from zero to computationDepth (mm)
  */
 double Crit1DCase::getWaterCapacitySum(double computationDepth)

criteriaModel/criteria1DCase.h

@@ -78,7 +78,8 @@
         bool initializeWaterContent(Crit3DDate myDate);
         bool computeDailyModel(Crit3DDate &myDate, std::string &error);
 
-        double getWaterContent(double computationDepth);
+        double getVolumetricWaterContent(double computationDepth);
+        double getDegreeOfSaturation(double computationDepth);
         double getWaterPotential(double computationDepth);
         double getFractionAW(double computationDepth);
         double getSlopeStability(double computationDepth);
@@ -99,7 +100,7 @@
         bool computeNumericalFluxes(const Crit3DDate &myDate, std::string &error);
         bool computeWaterFluxes(const Crit3DDate &myDate, std::string &error);
         double checkIrrigationDemand(int doy, double currentPrec, double nextPrec, double maxTranspiration);
-        void saveWaterContent();
+        void storeWaterContent();
         void restoreWaterContent();
         double getTotalWaterContent();
 

criteriaModel/criteria1DProject.cpp

@@ -63,9 +63,11 @@ void Crit1DProject::initialize()
     lastSimulationDate = QDate(1800,1,1);
 
     outputString = "";
+
     // specific outputs
     waterDeficitDepth.clear();
     waterContentDepth.clear();
+    degreeOfSaturationDepth.clear();
     waterPotentialDepth.clear();
     availableWaterDepth.clear();
     fractionAvailableWaterDepth.clear();
@@ -257,6 +259,12 @@ bool Crit1DProject::readSettings()
             projectError = "Wrong water content depth in " + configFileName;
             return false;
         }
+        depthList = projectSettings->value("degreeOfSaturation").toStringList();
+        if (! setVariableDepth(depthList, degreeOfSaturationDepth))
+        {
+            projectError = "Wrong degree of saturation depth in " + configFileName;
+            return false;
+        }
         depthList = projectSettings->value("waterPotential").toStringList();
         if (! setVariableDepth(depthList, waterPotentialDepth))
         {
@@ -1544,7 +1552,12 @@ bool Crit1DProject::createOutputTable(QString &myError)
     // specific depth variables
     for (unsigned int i = 0; i < waterContentDepth.size(); i++)
     {
-        QString fieldName = "SWC_" + QString::number(waterContentDepth[i]);
+        QString fieldName = "VWC_" + QString::number(waterContentDepth[i]);
+        queryString += ", " + fieldName + " REAL";
+    }
+    for (unsigned int i = 0; i < degreeOfSaturationDepth.size(); i++)
+    {
+        QString fieldName = "DEGSAT_" + QString::number(degreeOfSaturationDepth[i]);
         queryString += ", " + fieldName + " REAL";
     }
     for (unsigned int i = 0; i < waterPotentialDepth.size(); i++)
@@ -1604,7 +1617,12 @@ void Crit1DProject::updateOutput(Crit3DDate myDate, bool isFirst)
         // specific depth variables
         for (unsigned int i = 0; i < waterContentDepth.size(); i++)
         {
-            QString fieldName = "SWC_" + QString::number(waterContentDepth[i]);
+            QString fieldName = "VWC_" + QString::number(waterContentDepth[i]);
+            outputString += ", " + fieldName;
+        }
+        for (unsigned int i = 0; i < degreeOfSaturationDepth.size(); i++)
+        {
+            QString fieldName = "DEGSAT_" + QString::number(degreeOfSaturationDepth[i]);
             outputString += ", " + fieldName;
         }
         for (unsigned int i = 0; i < waterPotentialDepth.size(); i++)
@@ -1669,7 +1687,11 @@ void Crit1DProject::updateOutput(Crit3DDate myDate, bool isFirst)
     // specific depth variables
     for (unsigned int i = 0; i < waterContentDepth.size(); i++)
     {
-        outputString += "," + QString::number(myCase.getWaterContent(waterContentDepth[i]), 'g', 4);
+        outputString += "," + QString::number(myCase.getVolumetricWaterContent(waterContentDepth[i]), 'g', 4);
+    }
+    for (unsigned int i = 0; i < degreeOfSaturationDepth.size(); i++)
+    {
+        outputString += "," + QString::number(myCase.getDegreeOfSaturation(degreeOfSaturationDepth[i]), 'g', 4);
     }
     for (unsigned int i = 0; i < waterPotentialDepth.size(); i++)
     {

criteriaModel/criteria1DProject.h

@@ -91,6 +91,7 @@
 
         // specific output
         std::vector<int> waterContentDepth;
+        std::vector<int> degreeOfSaturationDepth;
         std::vector<int> waterPotentialDepth;
         std::vector<int> waterDeficitDepth;
         std::vector<int> awcDepth;