From 7916330d3f20ec1e04c77b9d3038ee5a78e21d9b Mon Sep 17 00:00:00 2001 From: ftomei Date: Fri, 1 Dec 2023 18:22:08 +0100 Subject: [PATCH] update roots --- agrolib/crop/root.cpp | 42 ++++++++++++++-------------- agrolib/crop/root.h | 3 +- bin/CRITERIA3D/criteria3DProject.cpp | 11 +++++--- bin/CRITERIA3D/shared/project3D.cpp | 21 ++++++++++++-- bin/CRITERIA3D/shared/project3D.h | 2 +- 5 files changed, 49 insertions(+), 30 deletions(-) diff --git a/agrolib/crop/root.cpp b/agrolib/crop/root.cpp index 3cbb53107..816d0a737 100644 --- a/agrolib/crop/root.cpp +++ b/agrolib/crop/root.cpp @@ -33,6 +33,7 @@ #include "commonConstants.h" #include "gammaFunction.h" +#include "basicMath.h" #include "root.h" #include "crop.h" @@ -186,10 +187,10 @@ namespace root const double MAX_DAILY_GROWTH = 0.02; // [m] const double MIN_WATERTABLE_DISTANCE = 0.2; // [m] - if (int(waterTableDepth) != int(NODATA) + if (! isEqual(waterTableDepth, NODATA) && waterTableDepth > 0 - && int(myCrop->roots.rootLength) != int(NODATA) - && !myCrop->isWaterSurplusResistant() + && ! isEqual(myCrop->roots.rootLength, NODATA) + && ! myCrop->isWaterSurplusResistant() && myRootLength > myCrop->roots.rootLength) { // check on growth @@ -288,7 +289,7 @@ namespace root } - int getNrAtoms(const std::vector &soilLayers, double &minThickness, int *atoms) + int getNrAtoms(const std::vector &soilLayers, double &minThickness, std::vector &atoms) { unsigned int nrLayers = unsigned(soilLayers.size()); int multiplicationFactor = 1; @@ -297,7 +298,7 @@ namespace root double tmp = minThickness * 1.001; if (tmp < 1) - multiplicationFactor = int(pow(10.0,-orderOfMagnitude(tmp))); + multiplicationFactor = int(pow(10.0, -orderOfMagnitude(tmp))); if (minThickness < 1) { @@ -312,6 +313,7 @@ namespace root atoms[i] = value; counter += value; } + return counter; } @@ -324,7 +326,7 @@ namespace root */ void cardioidDistribution(double shapeFactor, unsigned int nrLayersWithRoot, unsigned int nrUpperLayersWithoutRoot, unsigned int totalLayers, - double* densityThinLayers) + std::vector &densityThinLayers) { unsigned int i; std::vector lunette, lunetteDensity; @@ -382,15 +384,17 @@ namespace root void cylindricalDistribution(double deformation, unsigned int nrLayersWithRoot, unsigned int nrUpperLayersWithoutRoot, unsigned int totalLayers, - double* densityThinLayers) + std::vector &densityThinLayers) { unsigned int i; - double* cylinderDensity = new double[nrLayersWithRoot*2]; + std::vector cylinderDensity; + cylinderDensity.resize(nrLayersWithRoot*2); + // initialize not deformed cylinder for (i = 0 ; i < (2*nrLayersWithRoot); i++) { cylinderDensity[i]= 1./(2*nrLayersWithRoot); - } // not deformed cylinder + } // linear and ovoidal deformation double deltaDeformation,rootDensitySum; @@ -421,8 +425,6 @@ namespace root { densityThinLayers[nrUpperLayersWithoutRoot+i] = cylinderDensity[2*i] + cylinderDensity[2*i+1]; } - - delete[] cylinderDensity; } @@ -452,26 +454,27 @@ namespace root || (myCrop->roots.rootShape == CYLINDRICAL_DISTRIBUTION)) { double minimumThickness; - int* atoms = new int[nrLayers]; + std::vector atoms; + atoms.resize(nrLayers); int nrAtoms = root::getNrAtoms(soilLayers, minimumThickness, atoms); int numberOfRootedLayers, numberOfTopUnrootedLayers; numberOfTopUnrootedLayers = int(round(myCrop->roots.rootDepthMin / minimumThickness)); numberOfRootedLayers = int(round(MINVALUE(myCrop->roots.rootLength, soilDepth) / minimumThickness)); + // roots are still too short if (numberOfRootedLayers == 0) - { - delete[] atoms; return true; - } - // check nr thin layers + + // check nr of thin layers if ((numberOfTopUnrootedLayers + numberOfRootedLayers) > nrAtoms) { numberOfRootedLayers = nrAtoms - numberOfTopUnrootedLayers; } // initialize thin layers density - double* densityThinLayers = new double[nrAtoms]; + std::vector densityThinLayers; + densityThinLayers.resize(nrAtoms); for (int i=0; i < nrAtoms; i++) { densityThinLayers[i] = 0.; @@ -491,16 +494,13 @@ namespace root int counter = 0; for (unsigned int layer=0; layer < nrLayers; layer++) { - for (unsigned int j = 0; j < unsigned(atoms[layer]); j++) + for (int j = 0; j < atoms[layer]; j++) { if (counter < nrAtoms) myCrop->roots.rootDensity[layer] += densityThinLayers[counter]; counter++; } } - - delete[] atoms; - delete[] densityThinLayers; } else if (myCrop->roots.rootShape == GAMMA_DISTRIBUTION) { diff --git a/agrolib/crop/root.h b/agrolib/crop/root.h index e8aeb30ce..8681a2b26 100644 --- a/agrolib/crop/root.h +++ b/agrolib/crop/root.h @@ -47,7 +47,7 @@ namespace root { - int getNrAtoms(const std::vector &soilLayers, double &minThickness, int *atoms); + int getNrAtoms(const std::vector &soilLayers, double &minThickness, std::vector &atoms); double getRootLengthDD(Crit3DRoot* myRoot, double currentDD, double emergenceDD); rootDistributionType getRootDistributionType(int rootShape); int getRootDistributionNumber(rootDistributionType rootShape); @@ -57,6 +57,7 @@ double computeRootLength(Crit3DCrop* myCrop, double currentDD, double waterTableDepth); double computeRootDepth(Crit3DCrop* myCrop, double currentDD, double waterTableDepth); bool computeRootDensity(Crit3DCrop* myCrop, const std::vector &soilLayers); + } #endif // ROOT_H diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp index 7ee476dc0..60f8a1024 100644 --- a/bin/CRITERIA3D/criteria3DProject.cpp +++ b/bin/CRITERIA3D/criteria3DProject.cpp @@ -241,9 +241,13 @@ void Crit3DProject::assignETreal() totalEvaporation += evapFlow; // assign real transpiration - double realTransp = assignTranspiration(row, col, lai); // [mm] - double traspFlow = area * (realTransp / 1000.); // [m3 h-1] - totalTranspiration += traspFlow; + if (lai > 0) + { + float degreeDays = degreeDaysMap.value[row][col]; + double realTransp = assignTranspiration(row, col, lai, degreeDays); // [mm] + double traspFlow = area * (realTransp / 1000.); // [m3 h-1] + totalTranspiration += traspFlow; + } } } } @@ -1019,7 +1023,6 @@ bool Crit3DProject::modelHourlyCycle(QDateTime myTime, const QString& hourlyOutp waterSinkSource.at(size_t(i)) = 0.; } - if (processes.computeEvaporation || processes.computeCrop) { if (! hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps)) diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp index f20eeb487..9ffa4fe44 100644 --- a/bin/CRITERIA3D/shared/project3D.cpp +++ b/bin/CRITERIA3D/shared/project3D.cpp @@ -24,6 +24,7 @@ */ #include "commonConstants.h" +#include "basicMath.h" #include "cropDbTools.h" #include "project3D.h" #include "soilFluxes3D.h" @@ -1342,13 +1343,27 @@ double Project3D::assignEvaporation(int row, int col, double lai) } -// assign real crop transpiration +// assign actual crop transpiration // return sum of crop transpiration over the soil column -double Project3D::assignTranspiration(int row, int col, double lai) +double Project3D::assignTranspiration(int row, int col, double lai, double degreeDays) { double transpirationSum = 0; - // check + if (lai < EPSILON || isEqual(degreeDays, NODATA)) + return transpirationSum; + + // check land unit + int index = getLandUnitIndexRowCol(row, col); + if (index == NODATA) + return transpirationSum; + + // check crop + if (landUnitList[index].idCrop == "") + return transpirationSum; + + Crit3DCrop myCrop = cropList[index]; + // TODO assign roots from degreedays + /*if (idCrop == "" || ! isLiving) return 0; if (roots.rootDepth <= roots.rootDepthMin) return 0; if (roots.firstRootLayer == NODATA) return 0; diff --git a/bin/CRITERIA3D/shared/project3D.h b/bin/CRITERIA3D/shared/project3D.h index ba07f257a..236f37fb3 100644 --- a/bin/CRITERIA3D/shared/project3D.h +++ b/bin/CRITERIA3D/shared/project3D.h @@ -157,7 +157,7 @@ bool interpolateHourlyMeteoVar(meteoVariable myVar, const QDateTime& myTime); double assignEvaporation(int row, int col, double lai); - double assignTranspiration(int row, int col, double lai); + double assignTranspiration(int row, int col, double lai, double degreeDays); bool setSinkSource(); void computeWaterBalance3D(double totalTimeStep); bool updateCrop(QDateTime myTime);