From e3131f5ceb5db51ab3243daf8ec07312625b65eb Mon Sep 17 00:00:00 2001
From: ftomei <ftomei@arpae.it>
Date: Wed, 14 Aug 2024 13:02:47 +0200
Subject: [PATCH] soil cracks

---
 meteo/meteo.h                      |  2 +-
 soilFluxes3D/header/soilFluxes3D.h |  1 +
 soilFluxes3D/soilFluxes3D.cpp      | 74 ++++++++++++++++++++++--------
 3 files changed, 57 insertions(+), 20 deletions(-)

diff --git a/meteo/meteo.h b/meteo/meteo.h
index 312795546..5c6e4a78d 100644
--- 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
index d9724bc94..cd3550b56 100644
--- 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
index a46aa6537..f601d69ec 100644
--- a/soilFluxes3D/soilFluxes3D.cpp
+++ b/soilFluxes3D/soilFluxes3D.cpp
@@ -640,23 +640,47 @@ 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)
@@ -664,15 +688,21 @@ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve,
  */
  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));
+    }
  }
 
 
@@ -680,19 +710,25 @@ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve,
  * \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));
+        }
  }