Skip to content

Commit

Permalink
update transpiration
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Feb 21, 2024
1 parent 4b5468c commit 9bf26cd
Show file tree
Hide file tree
Showing 6 changed files with 121 additions and 104 deletions.
36 changes: 25 additions & 11 deletions agrolib/crop/crop.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -593,18 +593,16 @@ double Crit3DCrop::computeRootLength(double currentDD, double waterTableDepth)
}


/*! \brief updateRootDepth3D
* update current root lenght and root depth
/*! \brief computeRootLength3D
* compute 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 currentDegreeDays [°C]
* \param totalSoilDepth [m]
*/
void Crit3DCrop::updateRootDepth3D(double currentDD, double waterTableDepth, double previousRootDepth, double totalSoilDepth)
void Crit3DCrop::computeRootLength3D(double currentDegreeDays, double totalSoilDepth)
{
// set actualRootDepthMax
if (isEqual(totalSoilDepth, NODATA) || isEqual(totalSoilDepth, 0))
if ( isEqual(totalSoilDepth, NODATA) )
{
roots.actualRootDepthMax = roots.rootDepthMax;
}
Expand All @@ -614,16 +612,32 @@ void Crit3DCrop::updateRootDepth3D(double currentDD, double waterTableDepth, do
}

// set currentRootLength
if (isEqual(previousRootDepth, NODATA))
if (isRootStatic())
{
roots.currentRootLength = 0;
roots.currentRootLength = roots.actualRootDepthMax - roots.rootDepthMin;
}
else
{
roots.currentRootLength = previousRootDepth - roots.rootDepthMin;
if (currentDegreeDays <= 0)
{
roots.currentRootLength = 0.0;
}
else
{
if (currentDegreeDays > roots.degreeDaysRootGrowth)
{
roots.currentRootLength = roots.actualRootDepthMax - roots.rootDepthMin;
}
else
{
// in order to avoid numerical divergences
currentDegreeDays = MAXVALUE(currentDegreeDays, 1);
roots.currentRootLength = root::getRootLengthDD(roots, currentDegreeDays, degreeDaysEmergence);
}
}
}

roots.currentRootLength = computeRootLength(currentDD, waterTableDepth);
// set rootDepth
roots.rootDepth = roots.rootDepthMin + roots.currentRootLength;
}

Expand Down
2 changes: 1 addition & 1 deletion agrolib/crop/crop.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@
void updateRootDepth(double currentDD, double waterTableDepth);
double computeRootLength(double currentDD, double waterTableDepth);

void updateRootDepth3D(double currentDD, double waterTableDepth, double previousRootDepth, double totalSoilDepth);
void computeRootLength3D(double currentDD, double totalSoilDepth);

double computeSimpleLAI(double myDegreeDays, double latitude, int currentDoy);

Expand Down
54 changes: 28 additions & 26 deletions agrolib/crop/root.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -493,37 +493,40 @@ namespace root
}


bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil &currentSoil, unsigned int nrLayers,
bool computeRootDensity3D(Crit3DCrop &myCrop, const soil::Crit3DSoil &currentSoil, unsigned int nrLayers,
const std::vector<double> &layerDepth, const std::vector<double> &layerThickness)
{
// check soil
if (nrLayers == 0)
{
myCrop->roots.firstRootLayer = NODATA;
myCrop->roots.lastRootLayer = NODATA;
myCrop.roots.firstRootLayer = NODATA;
myCrop.roots.lastRootLayer = NODATA;
return false;
}

// Initialize
myCrop->roots.rootDensity.clear();
myCrop->roots.rootDensity.resize(nrLayers);
myCrop.roots.rootDensity.clear();
myCrop.roots.rootDensity.resize(nrLayers);
for (unsigned int i = 0; i < nrLayers; i++)
{
myCrop->roots.rootDensity[i] = 0.0;
myCrop.roots.rootDensity[i] = 0.0;
}

if (myCrop->roots.currentRootLength <= 0 )
if (myCrop.roots.currentRootLength <= 0 )
return true;

if (myCrop->roots.rootShape == GAMMA_DISTRIBUTION)
myCrop->roots.rootShape = CARDIOID_DISTRIBUTION;
// TODO Gamma distribuion
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));
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)
Expand All @@ -543,14 +546,14 @@ namespace root
densityThinLayers[i] = 0.;
}

if (myCrop->roots.rootShape == CARDIOID_DISTRIBUTION)
if (myCrop.roots.rootShape == CARDIOID_DISTRIBUTION)
{
cardioidDistribution(myCrop->roots.shapeDeformation, numberOfRootedLayers,
cardioidDistribution(myCrop.roots.shapeDeformation, numberOfRootedLayers,
numberOfTopUnrootedLayers, signed(nrAtoms), densityThinLayers);
}
else if (myCrop->roots.rootShape == CYLINDRICAL_DISTRIBUTION)
else if (myCrop.roots.rootShape == CYLINDRICAL_DISTRIBUTION)
{
cylindricalDistribution(myCrop->roots.shapeDeformation, numberOfRootedLayers,
cylindricalDistribution(myCrop.roots.shapeDeformation, numberOfRootedLayers,
numberOfTopUnrootedLayers, signed(nrAtoms), densityThinLayers);
}

Expand All @@ -566,7 +569,7 @@ namespace root
double lowerDepth = layerDepth[l] + layerThickness[l] * 0.5;
if (currentDepth >= upperDepth && currentDepth <= lowerDepth)
{
myCrop->roots.rootDensity[l] += densityThinLayers[atom];
myCrop.roots.rootDensity[l] += densityThinLayers[atom];
rootDensitySum += densityThinLayers[atom];
break;
}
Expand All @@ -586,8 +589,8 @@ namespace root
if (horIndex != int(NODATA))
{
double soilFraction = (1 - currentSoil.horizon[horIndex].coarseFragments);
myCrop->roots.rootDensity[l] *= soilFraction;
rootDensitySumSubset += myCrop->roots.rootDensity[l];
myCrop.roots.rootDensity[l] *= soilFraction;
rootDensitySumSubset += myCrop.roots.rootDensity[l];
}
}

Expand All @@ -597,28 +600,27 @@ namespace root
double ratio = rootDensitySum / rootDensitySumSubset;
for (unsigned int l=0 ; l < nrLayers; l++)
{
myCrop->roots.rootDensity[l] *= ratio;
myCrop.roots.rootDensity[l] *= ratio;
}
}

// find first and last root layers
myCrop->roots.firstRootLayer = 0;
myCrop.roots.firstRootLayer = 0;
unsigned int layer = 0;
while (layer < nrLayers && myCrop->roots.rootDensity[layer] == 0.0)
while (layer < nrLayers && myCrop.roots.rootDensity[layer] == 0.0)
{
layer++;
(myCrop->roots.firstRootLayer)++;
(myCrop.roots.firstRootLayer)++;
}

myCrop->roots.lastRootLayer = myCrop->roots.firstRootLayer;
while (layer < nrLayers && myCrop->roots.rootDensity[layer] != 0.0)
myCrop.roots.lastRootLayer = myCrop.roots.firstRootLayer;
while (layer < nrLayers && myCrop.roots.rootDensity[layer] != 0.0)
{
myCrop->roots.lastRootLayer = signed(layer);
myCrop.roots.lastRootLayer = signed(layer);
layer++;
}

return true;
}

}

3 changes: 2 additions & 1 deletion agrolib/crop/root.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,8 @@

double getRootLengthDD(const Crit3DRoot &myRoot, double currentDD, double emergenceDD);
bool computeRootDensity(Crit3DCrop* myCrop, const std::vector<soil::Crit3DLayer> &soilLayers);
bool computeRootDensity3D(Crit3DCrop* myCrop, const soil::Crit3DSoil &currentSoil, unsigned int nrLayers,

bool computeRootDensity3D(Crit3DCrop &myCrop, const soil::Crit3DSoil &currentSoil, unsigned int nrLayers,
const std::vector<double> &layerDepth, const std::vector<double> &layerThickness);

int getNrAtoms(const std::vector<soil::Crit3DLayer> &soilLayers, double &minThickness, std::vector<int> &atoms);
Expand Down
Loading

0 comments on commit 9bf26cd

Please sign in to comment.