Skip to content

Commit

Permalink
bare soil and update evaporation
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Jun 12, 2024
1 parent 7952ad0 commit 6642632
Show file tree
Hide file tree
Showing 6 changed files with 69 additions and 29 deletions.
Binary file modified DATA/PROJECT/test/data/crop.db
View file
Binary file not shown.
6 changes: 4 additions & 2 deletions agrolib/criteriaModel/carbonNitrogenModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1840,8 +1840,10 @@ void Crit1DCarbonNitrogenProfile::N_harvest(Crit1DCase &myCase) // public functi
// N of leaves is incorporated in litter through the upeer layer with a smoothly rate during the leaf fall

double N_toLitter = 0;
// !!! verificare USR PSR
if (myCase.crop.roots.firstRootLayer == 0 && myCase.crop.roots.lastRootLayer == 0)

if ((myCase.crop.roots.firstRootLayer == 0 && myCase.crop.roots.lastRootLayer == 0)
|| myCase.crop.roots.firstRootLayer == NODATA || myCase.crop.roots.lastRootLayer == NODATA
|| myCase.crop.isBareSoil() )
return;

for (int l = myCase.crop.roots.firstRootLayer; l <= myCase.crop.roots.lastRootLayer; l++) // verificare i cicli for per cambio indici
Expand Down
30 changes: 19 additions & 11 deletions agrolib/criteriaModel/water1D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -373,8 +373,6 @@ double computeCapillaryRise(std::vector<soil::Crit1DLayer> &soilLayers, double w

double computeEvaporation(std::vector<soil::Crit1DLayer> &soilLayers, double maxEvaporation)
{
// TODO extend to geometric soilLayers

// surface evaporation
double surfaceEvaporation = MINVALUE(maxEvaporation, soilLayers[0].waterContent);
soilLayers[0].waterContent -= surfaceEvaporation;
Expand All @@ -390,32 +388,43 @@ double computeEvaporation(std::vector<soil::Crit1DLayer> &soilLayers, double max
// soil evaporation
int nrEvapLayers = int(floor(MAX_EVAPORATION_DEPTH / soilLayers[1].thickness)) +1;
nrEvapLayers = std::min(nrEvapLayers, int(soilLayers.size()-1));
double* coeffEvap = new double[nrEvapLayers];
double layerDepth, coeffDepth;

double sumCoeff = 0;
std::vector<double> coeffEvap;
coeffEvap.resize(nrEvapLayers);

double minDepth = soilLayers[1].depth + soilLayers[1].thickness / 2;
double sumCoeff = 0;
for (int i=1; i <= nrEvapLayers; i++)
{
layerDepth = soilLayers[i].depth + soilLayers[i].thickness / 2.0;
double layerDepth = soilLayers[i].depth + soilLayers[i].thickness / 2.0;

coeffDepth = MAXVALUE((layerDepth - minDepth) / (MAX_EVAPORATION_DEPTH - minDepth), 0);
double coeffDepth = MAXVALUE((layerDepth - minDepth) / (MAX_EVAPORATION_DEPTH - minDepth), 0);
// evaporation coefficient: 1 at depthMin, ~0.1 at MAX_EVAPORATION_DEPTH
coeffEvap[i-1] = exp(-2 * coeffDepth);

sumCoeff += coeffEvap[i-1];
}

// normalize
std::vector<double> coeffThreshold;
coeffThreshold.resize(nrEvapLayers);
for (int i=0; i < nrEvapLayers; i++)
{
coeffThreshold[i] = (1.0 - coeffEvap[i]) * 0.5;
coeffEvap[i] /= sumCoeff;
}

bool isWaterSupply = true;
double sumEvap, evapLayerThreshold, evapLayer;
while ((residualEvaporation > EPSILON) && (isWaterSupply == true))
int nrIteration = 0;
while ((residualEvaporation > EPSILON) && (isWaterSupply == true) && nrIteration < 3)
{
sumEvap = 0.0;

for (int i=1; i <= nrEvapLayers; i++)
{
evapLayer = residualEvaporation * (coeffEvap[i-1] / sumCoeff);
evapLayerThreshold = soilLayers[i].HH + (1 - coeffEvap[i-1]) * (soilLayers[i].FC - soilLayers[i].HH);
evapLayerThreshold = soilLayers[i].HH + (soilLayers[i].FC - soilLayers[i].HH) * coeffThreshold[i-1];

if (soilLayers[i].waterContent > (evapLayerThreshold + evapLayer))
{
Expand All @@ -441,10 +450,9 @@ double computeEvaporation(std::vector<soil::Crit1DLayer> &soilLayers, double max

residualEvaporation -= sumEvap;
actualEvaporation += sumEvap;
nrIteration++;
}

delete[] coeffEvap;

return actualEvaporation;
}

Expand Down
32 changes: 26 additions & 6 deletions agrolib/crop/crop.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,21 +107,36 @@ void Crit3DCrop::initialize(double latitude, unsigned int nrLayers, double total
// initialize root depth
roots.rootDepth = 0;

if (totalSoilDepth == 0 || roots.rootDepthMax < totalSoilDepth)
roots.actualRootDepthMax = roots.rootDepthMax;
if (isBareSoil())
{
roots.rootDepthMax = 0;
roots.actualRootDepthMax = 0;
}
else
roots.actualRootDepthMax = totalSoilDepth;
{
if (totalSoilDepth == 0 || roots.rootDepthMax < totalSoilDepth)
{
roots.actualRootDepthMax = roots.rootDepthMax;
}
else
{
roots.actualRootDepthMax = totalSoilDepth;
}
}

degreeDays = 0;

if (latitude > 0)
{
doyStartSenescence = 305;
}
else
{
doyStartSenescence = 120;
}

LAIstartSenescence = NODATA;
currentSowingDoy = NODATA;

daysSinceIrrigation = NODATA;

// check if the crop is living
Expand All @@ -134,7 +149,7 @@ void Crit3DCrop::initialize(double latitude, unsigned int nrLayers, double total
}
else
{
isLiving = true;
isLiving = !isBareSoil();
}

resetCrop(nrLayers);
Expand Down Expand Up @@ -399,7 +414,7 @@ bool Crit3DCrop::needReset(Crit3DDate myDate, double latitude, double waterTable
void Crit3DCrop::resetCrop(unsigned int nrLayers)
{
// roots
if (! isRootStatic())
if (! isBareSoil() && ! isRootStatic())
{
for (unsigned int i = 0; i < nrLayers; i++)
roots.rootDensity[i] = 0;
Expand Down Expand Up @@ -855,10 +870,13 @@ speciesType getCropType(std::string cropType)
return FALLOW_ANNUAL;
else if (cropType == "tree" || cropType == "fruit_tree")
return TREE;
else if (cropType == "bare" || cropType == "bare_soil")
return BARESOIL;
else
return HERBACEOUS_ANNUAL;
}


std::string getCropTypeString(speciesType cropType)
{
switch (cropType)
Expand All @@ -877,6 +895,8 @@ std::string getCropTypeString(speciesType cropType)
return "fallow_annual";
case TREE:
return "tree";
case BARESOIL:
return "bare_soil";
}

return "No crop type";
Expand Down
5 changes: 4 additions & 1 deletion agrolib/crop/crop.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
#include "root.h"
#endif

enum speciesType {HERBACEOUS_ANNUAL, HERBACEOUS_PERENNIAL, HORTICULTURAL, GRASS, TREE, FALLOW, FALLOW_ANNUAL};
enum speciesType {HERBACEOUS_ANNUAL, HERBACEOUS_PERENNIAL, HORTICULTURAL, GRASS, TREE, FALLOW, FALLOW_ANNUAL, BARESOIL};
#define NR_CROP_SPECIES 7

/*!
Expand Down Expand Up @@ -81,6 +81,9 @@
bool isSowingCrop() const;
bool isRootStatic() const;

bool isBareSoil() const
{ return (type == BARESOIL); }

double getDailyDegreeIncrease(double tmin, double tmax, int doy);

void initialize(double latitude, unsigned int nrLayers, double totalSoilDepth, int currentDoy);
Expand Down
25 changes: 16 additions & 9 deletions agrolib/crop/cropDbTools.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,19 +87,18 @@ bool loadCropParameters(const QSqlDatabase &dbCrop, QString idCrop, Crit3DCrop &
myCrop.kcMax = query.value("kc_max").toDouble();
// [cm]
if (! getValue(query.value("psi_leaf"), &(myCrop.psiLeaf)))
{
// default
myCrop.psiLeaf = 16000;
}

myCrop.stressTolerance = query.value("stress_tolerance").toDouble();

// fraction of Readily Available Water
if (! getValue(query.value("raw_fraction"), &(myCrop.fRAW)))
{
// old version
if (! getValue(query.value("frac_read_avail_water_max"), &(myCrop.fRAW)))
{
errorStr = "Missing RAW_fraction for crop: " + idCropString;
return false;
}
// default
myCrop.fRAW = 0.6;
}

// IRRIGATION
Expand All @@ -109,17 +108,25 @@ bool loadCropParameters(const QSqlDatabase &dbCrop, QString idCrop, Crit3DCrop &
getValue(query.value("doy_start_irrigation"), &(myCrop.doyStartIrrigation));
getValue(query.value("doy_end_irrigation"), &(myCrop.doyEndIrrigation));

// key value for irrigation
// irrigation volume [mm day-1]
if (! getValue(query.value("irrigation_volume"), &(myCrop.irrigationVolume)))
{
// default: no irrigation
myCrop.irrigationVolume = 0;
}

// LAI grass
if (! getValue(query.value("lai_grass"), &(myCrop.LAIgrass)))
{
myCrop.LAIgrass = 0;
}

// max surface puddle
// max surface puddle [mm]
if (! getValue(query.value("max_height_surface_puddle"), &(myCrop.maxSurfacePuddle)))
myCrop.maxSurfacePuddle = 0;
{
// default: 5 mm
myCrop.maxSurfacePuddle = 5;
}

return true;
}
Expand Down

0 comments on commit 6642632

Please sign in to comment.