From b59996340a5869a79653cee40226e894294f9ca5 Mon Sep 17 00:00:00 2001
From: ftomei <ftomei@arpae.it>
Date: Thu, 19 Dec 2024 15:06:58 +0100
Subject: [PATCH] soil fraction

---
 criteriaModel/criteria1DCase.cpp | 7 ++++---
 crop/root.cpp                    | 3 +--
 crop/root.h                      | 4 +++-
 soil/soil.cpp                    | 2 +-
 soil/soil.h                      | 3 +--
 5 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/criteriaModel/criteria1DCase.cpp b/criteriaModel/criteria1DCase.cpp
index a3407902..e9ef2a70 100644
--- a/criteriaModel/criteria1DCase.cpp
+++ b/criteriaModel/criteria1DCase.cpp
@@ -161,16 +161,17 @@ bool Crit1DCase::initializeNumericalFluxes(std::string &error)
     for (unsigned int horizonIndex = 0; horizonIndex < mySoil.nrHorizons; horizonIndex++)
     {
         soil::Crit3DHorizon horizon = mySoil.horizon[horizonIndex];
-        double soilFraction = (1.0 - horizon.coarseFragments);
+
         result = soilFluxes3D::setSoilProperties(soilIndex, int(horizonIndex),
                             horizon.vanGenuchten.alpha * GRAVITY,
                             horizon.vanGenuchten.n, horizon.vanGenuchten.m,
                             horizon.vanGenuchten.he / GRAVITY,
-                            horizon.vanGenuchten.thetaR * soilFraction,
-                            horizon.vanGenuchten.thetaS * soilFraction,
+                            horizon.vanGenuchten.thetaR * horizon.getSoilFraction(),
+                            horizon.vanGenuchten.thetaS * horizon.getSoilFraction(),
                             (horizon.waterConductivity.kSat * 0.01) / DAY_SECONDS,
                             horizon.waterConductivity.l,
                             horizon.organicMatter, horizon.texture.clay * 0.01);
+
         if (result != CRIT3D_OK)
         {
             error = "Error in setSoilProperties, horizon nr: " + std::to_string(horizonIndex + 1);
diff --git a/crop/root.cpp b/crop/root.cpp
index 3eb4ad0a..17f9f1dd 100644
--- a/crop/root.cpp
+++ b/crop/root.cpp
@@ -588,8 +588,7 @@ namespace root
             int horIndex = currentSoil.getHorizonIndex(layerDepth[l]);
             if (horIndex != int(NODATA))
             {
-                double soilFraction = (1 - currentSoil.horizon[horIndex].coarseFragments);
-                myCrop.roots.rootDensity[l] *= soilFraction;
+                myCrop.roots.rootDensity[l] *= currentSoil.horizon[horIndex].getSoilFraction();
                 rootDensitySumSubset += myCrop.roots.rootDensity[l];
             }
         }
diff --git a/crop/root.h b/crop/root.h
index ddc6d6fd..9846a9d5 100644
--- a/crop/root.h
+++ b/crop/root.h
@@ -47,12 +47,14 @@
 
     namespace root
     {
-        rootDistributionType getRootDistributionType(int rootShape);
         int getRootDistributionNumber(rootDistributionType rootShape);
+
+        rootDistributionType getRootDistributionType(int rootShape);
         rootDistributionType getRootDistributionTypeFromString(const std::string &rootShape);
         std::string getRootDistributionTypeString(rootDistributionType rootType);
 
         double getRootLengthDD(const Crit3DRoot &myRoot, double currentDD, double emergenceDD);
+
         bool computeRootDensity(Crit3DCrop* myCrop, const std::vector<soil::Crit1DLayer> &soilLayers);
 
         bool computeRootDensity3D(Crit3DCrop &myCrop, const soil::Crit3DSoil &currentSoil, unsigned int nrLayers,
diff --git a/soil/soil.cpp b/soil/soil.cpp
index 1244aeb9..1d1e261e 100644
--- a/soil/soil.cpp
+++ b/soil/soil.cpp
@@ -227,7 +227,7 @@ namespace soil
         double waterContentHH = soil::thetaFromSignPsi(hygroscopicHumidity, *horizonPtr);
 
         // [-]
-        soilFraction = (1.0 - horizonPtr->coarseFragments);
+        soilFraction = horizonPtr->getSoilFraction();
 
         // [mm]
         SAT = horizonPtr->vanGenuchten.thetaS * soilFraction * thickness * 1000;
diff --git a/soil/soil.h b/soil/soil.h
index a889c10e..236fff3b 100644
--- a/soil/soil.h
+++ b/soil/soil.h
@@ -154,8 +154,7 @@
 
             Crit3DHorizon();
 
-            double getSoilFraction()
-                { return (1.0 - coarseFragments); }
+            double getSoilFraction() const { return (1.0 - coarseFragments); }
         };