From 3d40beb0686a00a844fd49ea78f621be1f2c2585 Mon Sep 17 00:00:00 2001 From: ftomei Date: Thu, 7 Dec 2023 18:45:54 +0100 Subject: [PATCH] root density nel 3D --- agrolib/crop/root.cpp | 127 ++++++++++++++++++++++++++++ agrolib/crop/root.h | 2 + agrolib/soil/soil.cpp | 2 +- agrolib/soil/soil.h | 2 +- bin/CRITERIA3D/shared/project3D.cpp | 7 ++ 5 files changed, 138 insertions(+), 2 deletions(-) diff --git a/agrolib/crop/root.cpp b/agrolib/crop/root.cpp index 3109aef9d..93a0d6707 100644 --- a/agrolib/crop/root.cpp +++ b/agrolib/crop/root.cpp @@ -492,5 +492,132 @@ namespace root return true; } + + bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil ¤tSoil, int nrLayers, + const std::vector &layerDepth, const std::vector &layerThickness) + { + // check soil + if (nrLayers == 0) + { + myCrop->roots.firstRootLayer = NODATA; + myCrop->roots.lastRootLayer = NODATA; + return false; + } + + // Initialize + myCrop->roots.rootDensity.clear(); + myCrop->roots.rootDensity.resize(nrLayers); + for (unsigned int i = 0; i < nrLayers; i++) + { + myCrop->roots.rootDensity[i] = 0.0; + } + + if (myCrop->roots.currentRootLength <= 0 ) + return true; + + if (myCrop->roots.rootShape == GAMMA_DISTRIBUTION) + myCrop->roots.rootShape = CARDIOID_DISTRIBUTION; + + int nrAtoms = int(currentSoil.totalDepth * 100) + 1; + double minimumThickness = 0.01; // [m] + + int numberOfRootedLayers, numberOfTopUnrootedLayers; + numberOfTopUnrootedLayers = int(round(myCrop->roots.rootDepthMin / minimumThickness)); + numberOfRootedLayers = int(round(MINVALUE(myCrop->roots.currentRootLength, currentSoil.totalDepth) / minimumThickness)); + + // roots are still too short + if (numberOfRootedLayers == 0) + return true; + + // check nr of thin layers + if ((numberOfTopUnrootedLayers + numberOfRootedLayers) > nrAtoms) + { + numberOfRootedLayers = nrAtoms - numberOfTopUnrootedLayers; + } + + // initialize thin layers density + std::vector densityThinLayers; + densityThinLayers.resize(nrAtoms); + for (int i=0; i < nrAtoms; i++) + { + densityThinLayers[i] = 0.; + } + + if (myCrop->roots.rootShape == CARDIOID_DISTRIBUTION) + { + cardioidDistribution(myCrop->roots.shapeDeformation, numberOfRootedLayers, + numberOfTopUnrootedLayers, signed(nrAtoms), densityThinLayers); + } + else if (myCrop->roots.rootShape == CYLINDRICAL_DISTRIBUTION) + { + cylindricalDistribution(myCrop->roots.shapeDeformation, numberOfRootedLayers, + numberOfTopUnrootedLayers, signed(nrAtoms), densityThinLayers); + } + + + double maxLayerDepth = layerDepth[nrLayers-1] + layerThickness[nrLayers-1] * 0.5; + int atom = 0; + double currentDepth = atom / 100.; // [m] + double rootDensitySum = 0.; + while (currentDepth <= maxLayerDepth && atom < nrAtoms) + { + for (int l = 0; l < nrLayers; l++) + { + double upperDepth = layerDepth[l] - layerThickness[l] * 0.5; + double lowerDepth = layerDepth[l] + layerThickness[l] * 0.5; + if (currentDepth >= upperDepth && currentDepth <= lowerDepth) + { + myCrop->roots.rootDensity[l] += densityThinLayers[atom]; + rootDensitySum += densityThinLayers[atom]; + break; + } + } + + atom++; + currentDepth = atom / 100.; // [m] + } + + if (rootDensitySum <= EPSILON) + return true; + + double rootDensitySumSubset = 0.; + for (unsigned int l=0 ; l < nrLayers; l++) + { + int horIndex = currentSoil.getHorizonIndex(layerDepth[l]); + if (horIndex != int(NODATA)) + { + double soilFraction = (1 - currentSoil.horizon[horIndex].coarseFragments); + myCrop->roots.rootDensity[l] *= soilFraction; + rootDensitySumSubset += myCrop->roots.rootDensity[l]; + } + } + + if (rootDensitySumSubset != rootDensitySum) + { + double ratio = rootDensitySum / rootDensitySumSubset; + for (unsigned int l=0 ; l < nrLayers; l++) + { + myCrop->roots.rootDensity[l] *= ratio; + } + } + + myCrop->roots.firstRootLayer = 0; + unsigned int layer = 0; + while (layer < nrLayers && myCrop->roots.rootDensity[layer] == 0.0) + { + layer++; + (myCrop->roots.firstRootLayer)++; + } + + myCrop->roots.lastRootLayer = myCrop->roots.firstRootLayer; + while (layer < nrLayers && myCrop->roots.rootDensity[layer] != 0.0) + { + myCrop->roots.lastRootLayer = signed(layer); + layer++; + } + + return true; + } + } diff --git a/agrolib/crop/root.h b/agrolib/crop/root.h index 64330852b..019d9b295 100644 --- a/agrolib/crop/root.h +++ b/agrolib/crop/root.h @@ -54,6 +54,8 @@ double getRootLengthDD(const Crit3DRoot &myRoot, double currentDD, double emergenceDD); bool computeRootDensity(Crit3DCrop* myCrop, const std::vector &soilLayers); + bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil ¤tSoil, int nrLayers, + const std::vector &layerDepth, const std::vector &layerThickness); int getNrAtoms(const std::vector &soilLayers, double &minThickness, std::vector &atoms); } diff --git a/agrolib/soil/soil.cpp b/agrolib/soil/soil.cpp index dbafbad96..a8dea9075 100644 --- a/agrolib/soil/soil.cpp +++ b/agrolib/soil/soil.cpp @@ -188,7 +188,7 @@ namespace soil nrHorizons = nrHorizons - 1; } - int Crit3DSoil::getHorizonIndex(double depth) + int Crit3DSoil::getHorizonIndex(double depth) const { for (unsigned int index = 0; index < nrHorizons; index++) { diff --git a/agrolib/soil/soil.h b/agrolib/soil/soil.h index d2bb33e76..8d8b66fb0 100644 --- a/agrolib/soil/soil.h +++ b/agrolib/soil/soil.h @@ -201,7 +201,7 @@ void cleanSoil(); void addHorizon(int nHorizon, const Crit3DHorizon &newHorizon); void deleteHorizon(int nHorizon); - int getHorizonIndex(double depth); + int getHorizonIndex(double depth) const; bool setSoilLayers(double layerThicknessMin, double geometricFactor, std::vector &soilLayers, std::string &myError); diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp index 4cd808743..86207e474 100644 --- a/bin/CRITERIA3D/shared/project3D.cpp +++ b/bin/CRITERIA3D/shared/project3D.cpp @@ -32,6 +32,7 @@ #include "math.h" #include "utilities.h" +#include "root.h" #include #include @@ -1388,6 +1389,12 @@ double Project3D::assignTranspiration(int row, int col, double lai, double degre if (cropList[cropIndex].roots.currentRootLength <= 0) return actualTranspiration; + // update root density + if (! root::computeRootDensity3D(&(cropList[cropIndex]), soilList[soilIndex], nrLayers, layerDepth, layerThickness)) + { + return actualTranspiration; + } + // TODO assign root density from degreedays /* if (roots.firstRootLayer == NODATA) return 0;