From 8dd89ff04379d758f2d22f6e1f501606f2f8ef13 Mon Sep 17 00:00:00 2001
From: Fausto Tomei <ftomei@arpae.it>
Date: Sat, 2 Dec 2023 19:43:57 +0100
Subject: [PATCH] update 3D

---
 crop/crop.cpp | 49 ++++++++++++++++++++++++++++++++++++++++++++-----
 crop/crop.h   |  2 ++
 2 files changed, 46 insertions(+), 5 deletions(-)

diff --git a/crop/crop.cpp b/crop/crop.cpp
index 9b03f4492..dc95b0658 100644
--- a/crop/crop.cpp
+++ b/crop/crop.cpp
@@ -528,7 +528,7 @@ void Crit3DCrop::updateRootDepth(double currentDD, double waterTableDepth)
 // return current root lenght [m]
 double Crit3DCrop::computeRootLength(double currentDD, double waterTableDepth)
 {
-    double newRootLength = NODATA;
+    double newRootLength;
 
     if (isRootStatic())
     {
@@ -555,19 +555,23 @@ double Crit3DCrop::computeRootLength(double currentDD, double waterTableDepth)
         }
     }
 
+    if (isEqual(waterTableDepth, NODATA))
+    {
+        return newRootLength;
+    }
+
     // WATERTABLE
     // Nel saturo le radici vanno in asfissia
     // per cui si mantengono a distanza dalla falda nella fase di crescita
-    // le radici possono crescere se:
-    // la falda è più bassa o si abbassa (max 2 cm al giorno)
+    // le radici possono crescere (max 2 cm al giorno) se:
+    // la falda è più bassa o si sta abbassando
     // restano invariate se:
     // 1) non sono più in fase di crescita
-    // 2) se sono già dentro la falda
+    // 2) sono già dentro la falda (currentRootDepth > waterTableDepth)
     const double MAX_DAILY_GROWTH = 0.02;             // [m]
     const double MIN_WATERTABLE_DISTANCE = 0.1;       // [m]
 
     if (! isWaterSurplusResistant()
-        && ! isEqual(waterTableDepth, NODATA)
         && ! isEqual(roots.currentRootLength, NODATA)
         && newRootLength > roots.currentRootLength)
     {
@@ -589,6 +593,41 @@ double Crit3DCrop::computeRootLength(double currentDD, double waterTableDepth)
 }
 
 
+/*! \brief updateRootDepth3D
+ *  update current root lenght and root depth
+ *  function for Criteria3D (update key variables)
+ *  \param currentDD:  current degree days sum
+ *  \param waterTableDepth      [m]
+ *  \param previousRootDepth    [m]
+ *  \param totalSoilDepth       [m]
+ */
+void Crit3DCrop::updateRootDepth3D(double currentDD,  double waterTableDepth, double previousRootDepth, double totalSoilDepth)
+{
+    // set actualRootDepthMax
+    if (isEqual(totalSoilDepth, NODATA) || isEqual(totalSoilDepth, 0))
+    {
+        roots.actualRootDepthMax = roots.rootDepthMax;
+    }
+    else
+    {
+        roots.actualRootDepthMax = std::min(roots.rootDepthMax, totalSoilDepth);
+    }
+
+    // set currentRootLength
+    if (isEqual(previousRootDepth, NODATA))
+    {
+        roots.currentRootLength = 0;
+    }
+    else
+    {
+        roots.currentRootLength = previousRootDepth - roots.rootDepthMin;
+    }
+
+    roots.currentRootLength = computeRootLength(currentDD, waterTableDepth);
+    roots.rootDepth = roots.rootDepthMin + roots.currentRootLength;
+}
+
+
 /*! \brief getCoveredSurfaceFraction
  *  \ref Liangxia Zhang, Zhongmin Hu, Jiangwen Fan, Decheng Zhou & Fengpei Tang, 2014
  *  A meta-analysis of the canopy light extinction coefficient in terrestrial ecosystems
diff --git a/crop/crop.h b/crop/crop.h
index 7dbbf1cf9..57b45a4a5 100644
--- a/crop/crop.h
+++ b/crop/crop.h
@@ -91,6 +91,8 @@
         void updateRootDepth(double currentDD, double waterTableDepth);
         double computeRootLength(double currentDD, double waterTableDepth);
 
+        void updateRootDepth3D(double currentDD, double waterTableDepth, double previousRootDepth, double totalSoilDepth);
+
         double computeSimpleLAI(double myDegreeDays, double latitude, int currentDoy);
 
         bool dailyUpdate(const Crit3DDate &myDate, double latitude, const std::vector<soil::Crit3DLayer> &soilLayers,