soil cracks

ARPA-SIMC · Aug 14, 2024 · e3131f5 · e3131f5
1 parent 8b67014
commit e3131f5
Show file tree

Hide file tree

Showing 3 changed files with 57 additions and 20 deletions.
diff --git a/meteo/meteo.h b/meteo/meteo.h
@@ -108,7 +108,7 @@
     enum criteria3DVariable {volumetricWaterContent, waterTotalPotential, waterMatricPotential,
                               availableWaterContent, degreeOfSaturation, soilTemperature,
                               soilSurfaceMoisture, bottomDrainage, waterDeficit, waterInflow, waterOutflow,
-                              factorOfSafety, minimumFactorOfSafety, surfacePond};
+                              factorOfSafety, minimumFactorOfSafety, surfacePond, maximumVolumetricWaterContent};
 
 
     const std::map<std::string, meteoVariable> MapDailyMeteoVar = {

diff --git a/soilFluxes3D/header/soilFluxes3D.h b/soilFluxes3D/header/soilFluxes3D.h
@@ -59,6 +59,7 @@
     __EXTERN int DLL_EXPORT __STDCALL setWaterSinkSource(long index, double sinkSource);
 
     __EXTERN double DLL_EXPORT __STDCALL getWaterContent(long nodeIndex);
+    __EXTERN double DLL_EXPORT __STDCALL getMaximumWaterContent(long nodeIndex);
     __EXTERN double DLL_EXPORT __STDCALL getAvailableWaterContent(long nodeIndex);
     __EXTERN double DLL_EXPORT __STDCALL getWaterDeficit(long index, double fieldCapacity);
     __EXTERN double DLL_EXPORT __STDCALL getTotalWaterContent();

diff --git a/soilFluxes3D/soilFluxes3D.cpp b/soilFluxes3D/soilFluxes3D.cpp
@@ -640,59 +640,95 @@ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve,
 /*!
  * \brief getWaterContent
  * \param nodeIndex
- * \return water content at the surface: [m] surface water level; and sub-surface: [m^3 m^-3] volumetric water content
+ * \return water content at the surface: [m] surface water level; and sub-surface: [m3 m-3] volumetric water content
  */
  double DLL_EXPORT __STDCALL getWaterContent(long nodeIndex)
  {
-        if (nodeListPtr == nullptr) return(MEMORY_ERROR);
-        if ((nodeIndex < 0) || (nodeIndex >= myStructure.nrNodes)) return(INDEX_ERROR);
+        if (nodeListPtr == nullptr)
+            return MEMORY_ERROR;
+        if ((nodeIndex < 0) || (nodeIndex >= myStructure.nrNodes))
+            return INDEX_ERROR;
 
         if  (nodeListPtr[nodeIndex].isSurface)
+        {
             /*! surface */
             return (nodeListPtr[nodeIndex].H - nodeListPtr[nodeIndex].z);
+        }
         else
+        {
             /*! sub-surface */
-            return (theta_from_Se(nodeIndex));
+            return theta_from_Se(nodeIndex);
+        }
  }
 
 
-/*!
+ /*!
+ * \brief getMaximumWaterContent
+ * \param nodeIndex
+ * \return only sub-surface: [m3 m-3] maximum volumetric water content (theta sat)
+ */
+ double DLL_EXPORT __STDCALL getMaximumWaterContent(long nodeIndex)
+ {
+     if (nodeListPtr == nullptr)
+         return MEMORY_ERROR;
+     if ((nodeIndex < 0) || (nodeIndex >= myStructure.nrNodes))
+         return INDEX_ERROR;
+     if  (nodeListPtr[nodeIndex].isSurface)
+         return INDEX_ERROR;
+
+    return (nodeListPtr[nodeIndex].Soil->Theta_s);
+ }
+
+
+ /*!
  * \brief getAvailableWaterContent
  * \param index
  * \return  available water content (over wilting point)
  * surface: water level [m]; sub-surface: awc [m3 m-3]
  */
  double DLL_EXPORT __STDCALL getAvailableWaterContent(long index)
  {
-        if (nodeListPtr == nullptr) return(MEMORY_ERROR);
-        if ((index < 0) || (index >= myStructure.nrNodes)) return(INDEX_ERROR);
+    if (nodeListPtr == nullptr)
+        return MEMORY_ERROR;
+    if ((index < 0) || (index >= myStructure.nrNodes))
+        return INDEX_ERROR;
 
-        if  (nodeListPtr[index].isSurface)
-            /*! surface */
-            return (nodeListPtr[index].H - nodeListPtr[index].z);
-        else
-            /*! sub-surface */
-            return MAXVALUE(0.0, theta_from_Se(index) - theta_from_sign_Psi(-160, index));
+    if  (nodeListPtr[index].isSurface)
+    {
+        // surface
+        return (nodeListPtr[index].H - nodeListPtr[index].z);
+    }
+    else
+    {
+        // sub-surface
+        return MAXVALUE(0.0, theta_from_Se(index) - theta_from_sign_Psi(-160, index));
+    }
  }
 
 
 /*!
  * \brief getWaterDeficit
  * \param index
  * \param fieldCapacity
- * \return water deficit at surface: 0; sub-surface: [m^3 m^-3]
+ * \return water deficit at surface: 0; sub-surface: [m3 m-3]
  */
 	double DLL_EXPORT __STDCALL getWaterDeficit(long index, double fieldCapacity)
  {
-        if (nodeListPtr == nullptr) return(MEMORY_ERROR);
-        if ((index < 0) || (index >= myStructure.nrNodes)) return(INDEX_ERROR);
+        if (nodeListPtr == nullptr)
+            return MEMORY_ERROR;
+        if ((index < 0) || (index >= myStructure.nrNodes))
+            return INDEX_ERROR;
 
         if  (nodeListPtr[index].isSurface)
-            /*! surface */
-            return (0.0);
+        {
+            // surface
+            return 0;
+        }
         else
-            /*! sub-surface */
+        {
+            // sub-surface
             return (theta_from_sign_Psi(-fieldCapacity, index) - theta_from_Se(index));
+        }
  }