From 26cc6be85b7637ce776b3394f9a9e171f501231a Mon Sep 17 00:00:00 2001 From: ftomei Date: Fri, 24 Nov 2023 17:42:05 +0100 Subject: [PATCH] add crop transpiration --- DATA/PROJECT/Montue/DATA/crop_Montue.db | Bin 13312 -> 13312 bytes bin/CRITERIA3D/criteria3DProject.cpp | 81 ++++++++-- bin/CRITERIA3D/criteria3DProject.h | 3 +- bin/CRITERIA3D/mainwindow.cpp | 7 + bin/CRITERIA3D/mainwindow.ui | 6 +- bin/CRITERIA3D/shared/project3D.cpp | 197 ++++++++++++++++++++---- bin/CRITERIA3D/shared/project3D.h | 7 +- 7 files changed, 257 insertions(+), 44 deletions(-) diff --git a/DATA/PROJECT/Montue/DATA/crop_Montue.db b/DATA/PROJECT/Montue/DATA/crop_Montue.db index f44f7f542364d5d678f16f95725d96e183385a69..4e8aff6e49fc025842f65146d0d50f66bd7ae3a0 100644 GIT binary patch delta 17 YcmZq3Xvml#%@{CI#+fl-W5NPs05-b?@&Et; delta 17 YcmZq3Xvml#&FD8##+lJ?W5NPs05+Wk?f?J) diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp index 235bc0e4c..97dd3e8f3 100644 --- a/bin/CRITERIA3D/criteria3DProject.cpp +++ b/bin/CRITERIA3D/criteria3DProject.cpp @@ -50,6 +50,7 @@ void Crit3DProcesses::initialize() computeMeteo = false; computeRadiation = false; computeWater = false; + computeEvaporation = false; computeCrop = false; computeSnow = false; computeSolutes = false; @@ -227,6 +228,45 @@ void Crit3DProject::dailyUpdateCrop() } +/*! + * \brief computeRealET + * assign soil evaporation and crop transpiration for the whole domain + */ +void Crit3DProject::computeRealET() +{ + totalEvaporation = 0; + totalTranspiration = 0; + + double area = DEM.header->cellSize * DEM.header->cellSize; + + for (int row = 0; row < indexMap.at(0).header->nrRows; row++) + { + for (int col = 0; col < indexMap.at(0).header->nrCols; col++) + { + int surfaceIndex = indexMap.at(0).value[row][col]; + if (surfaceIndex != indexMap.at(0).header->flag) + { + float lai = laiMap.value[row][col]; + if (isEqual(lai, NODATA)) + { + lai = 0; + } + + // assign real evaporation + double realEvap = assignEvaporation(row, col, lai); // [mm] + double flow = area * (realEvap / 1000.); // [m3 h-1] + totalEvaporation += flow; + + // assign real transpiration + double realTransp = assignTranspiration(row, col, lai); // [mm] + flow = area * (realTransp / 1000.); // [m3 h-1] + totalTranspiration += flow; + } + } + } +} + + bool Crit3DProject::runModels(QDateTime firstTime, QDateTime lastTime) { // initialize meteo @@ -947,36 +987,55 @@ bool Crit3DProject::modelHourlyCycle(QDateTime myTime, const QString& hourlyOutp qApp->processEvents(); } - if (processes.computeCrop) + if (processes.computeSnow) { - updateDailyTemperatures(); - - if (! hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps)) + // check conflitti con evaporation + // check snowmelt -> surface H0 + if (! computeSnowModel()) { return false; } + qApp->processEvents(); + } + + // initalize sink / source + for (unsigned long i = 0; i < nrNodes; i++) + { + waterSinkSource.at(size_t(i)) = 0.; + } + + + if (processes.computeEvaporation || processes.computeCrop) + { + if (! hourlyMeteoMaps->computeET0PMMap(DEM, radiationMaps)) + return false; if (isSaveOutputRaster()) { saveHourlyMeteoOutput(referenceEvapotranspiration, hourlyOutputPath, myTime); } + } - qApp->processEvents(); - - // TODO compute evap/transp + if (processes.computeEvaporation && (! processes.computeCrop)) + { + computeBareSoilEvaporation(); } - if (processes.computeSnow) + if (processes.computeCrop) { - if (! computeSnowModel()) - return false; + updateDailyTemperatures(); + + computeRealET(); + qApp->processEvents(); } // soil water balance if (processes.computeWater) { - if (! computeWaterSinkSource()) + assignPrecipitation(); + + if (! setSinkSource()) { logError(); return false; diff --git a/bin/CRITERIA3D/criteria3DProject.h b/bin/CRITERIA3D/criteria3DProject.h index 5b4a7ba9c..6bafbf728 100644 --- a/bin/CRITERIA3D/criteria3DProject.h +++ b/bin/CRITERIA3D/criteria3DProject.h @@ -28,7 +28,7 @@ public: bool computeMeteo, computeRadiation, computeWater; - bool computeCrop, computeSnow, computeSolutes; + bool computeEvaporation, computeCrop, computeSnow, computeSolutes; bool computeHeat, computeAdvectiveHeat, computeLatentHeat; Crit3DProcesses(); @@ -71,6 +71,7 @@ bool initializeCriteria3DModel(); void initializeCrop(); void dailyUpdateCrop(); + void computeRealET(); bool runModels(QDateTime firstTime, QDateTime lastTime); diff --git a/bin/CRITERIA3D/mainwindow.cpp b/bin/CRITERIA3D/mainwindow.cpp index 57a781ce8..87d0e4afb 100644 --- a/bin/CRITERIA3D/mainwindow.cpp +++ b/bin/CRITERIA3D/mainwindow.cpp @@ -1922,6 +1922,11 @@ bool MainWindow::startModels(QDateTime firstTime, QDateTime lastTime) if (myProject.processes.computeCrop) { // TODO: check on crop + if (myProject.landUnitList.size() == 0) + { + myProject.logError("load land units map before."); + return false; + } } // Load meteo data @@ -2213,6 +2218,7 @@ void MainWindow::on_actionCriteria3D_compute_current_hour_triggered() myProject.processes.computeMeteo = true; myProject.processes.computeRadiation = true; myProject.processes.computeWater = true; + myProject.processes.computeEvaporation = true; myProject.processes.computeCrop = true; startModels(currentTime, currentTime); @@ -2235,6 +2241,7 @@ void MainWindow::on_actionCriteria3D_run_models_triggered() myProject.processes.computeMeteo = true; myProject.processes.computeRadiation = true; myProject.processes.computeWater = false; + myProject.processes.computeEvaporation = true; myProject.processes.computeCrop = true; startModels(firstTime, lastTime); diff --git a/bin/CRITERIA3D/mainwindow.ui b/bin/CRITERIA3D/mainwindow.ui index 4d884ba40..ff9dd946a 100644 --- a/bin/CRITERIA3D/mainwindow.ui +++ b/bin/CRITERIA3D/mainwindow.ui @@ -1293,7 +1293,7 @@ 0 1752 10 - 2 + 1 @@ -1848,7 +1848,7 @@ 0 0 1280 - 29 + 23 @@ -2221,8 +2221,8 @@ - + diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp index ba94694bd..163ab2380 100644 --- a/bin/CRITERIA3D/shared/project3D.cpp +++ b/bin/CRITERIA3D/shared/project3D.cpp @@ -983,7 +983,8 @@ bool Project3D::aggregateAndSaveDailyMap(meteoVariable myVar, aggregationMethod // compute evaporation [mm] -double Project3D::computeEvaporation(int row, int col, double lai) +// TODO check pond +double Project3D::assignEvaporation(int row, int col, double lai) { double depthCoeff, thickCoeff, layerCoeff; double availableWater; // [mm] @@ -1044,6 +1045,138 @@ double Project3D::computeEvaporation(int row, int col, double lai) } +// assign real crop transpiration +// return sum of crop transpiration over the soil column +double Project3D::assignTranspiration(int row, int col, double lai) +{ + double soilColumnTranspirationSum = 0; + + // check + /*if (idCrop == "" || ! isLiving) return 0; + if (roots.rootDepth <= roots.rootDepthMin) return 0; + if (roots.firstRootLayer == NODATA) return 0; + if (maxTranspiration < EPSILON) return 0; + + double thetaWP; // [m3 m-3] volumetric water content at Wilting Point + double cropWP; // [mm] wilting point specific for crop + double waterSurplusThreshold; // [mm] water surplus stress threshold + double waterScarcityThreshold; // [mm] water scarcity stress threshold + double WSS; // [] water surplus stress + + double TRs=0.0; // [mm] actual transpiration with only water scarsity stress + double TRe=0.0; // [mm] actual transpiration with only water surplus stress + double totRootDensityWithoutStress = 0.0; // [-] + double redistribution = 0.0; // [mm] + + // initialize + unsigned int nrLayers = unsigned(soilLayers.size()); + bool* isLayerStressed = new bool[nrLayers]; + for (unsigned int i = 0; i < nrLayers; i++) + { + isLayerStressed[i] = false; + layerTranspiration[i] = 0; + } + + // water surplus + if (isWaterSurplusResistant()) + WSS = 0.0; + else + WSS = 0.5; + + for (unsigned int i = unsigned(roots.firstRootLayer); i <= unsigned(roots.lastRootLayer); i++) + { + // [mm] + waterSurplusThreshold = soilLayers[i].SAT - (WSS * (soilLayers[i].SAT - soilLayers[i].FC)); + + thetaWP = soil::thetaFromSignPsi(-soil::cmTokPa(psiLeaf), *(soilLayers[i].horizonPtr)); + // [mm] + cropWP = thetaWP * soilLayers[i].thickness * soilLayers[i].soilFraction * 1000.0; + + // [mm] + waterScarcityThreshold = soilLayers[i].FC - fRAW * (soilLayers[i].FC - cropWP); + + if ((soilLayers[i].waterContent - waterSurplusThreshold) > EPSILON) + { + // WATER SURPLUS + layerTranspiration[i] = maxTranspiration * roots.rootDensity[i] * + ((soilLayers[i].SAT - soilLayers[i].waterContent) + / (soilLayers[i].SAT - waterSurplusThreshold)); + + TRe += layerTranspiration[i]; + TRs += maxTranspiration * roots.rootDensity[i]; + isLayerStressed[i] = true; + } + else if (soilLayers[i].waterContent < waterScarcityThreshold) + { + // WATER SCARSITY + if (soilLayers[i].waterContent <= cropWP) + { + layerTranspiration[i] = 0; + } + else + { + layerTranspiration[i] = maxTranspiration * roots.rootDensity[i] * + ((soilLayers[i].waterContent - cropWP) / (waterScarcityThreshold - cropWP)); + } + + TRs += layerTranspiration[i]; + TRe += maxTranspiration * roots.rootDensity[i]; + isLayerStressed[i] = true; + } + else + { + // normal conditions + layerTranspiration[i] = maxTranspiration * roots.rootDensity[i]; + + TRs += layerTranspiration[i]; + TRe += layerTranspiration[i]; + + if ((soilLayers[i].waterContent - layerTranspiration[i]) > waterScarcityThreshold) + { + isLayerStressed[i] = false; + totRootDensityWithoutStress += roots.rootDensity[i]; + } + else + { + isLayerStressed[i] = true; + } + } + } + + // WATER STRESS [-] + double firstWaterStress = 1 - (TRs / maxTranspiration); + + // Hydraulic redistribution + // the movement of water from moist to dry soil through plant roots + // TODO add numerical process + if (firstWaterStress > EPSILON && totRootDensityWithoutStress > EPSILON) + { + // redistribution acts on not stressed roots + redistribution = MINVALUE(firstWaterStress, totRootDensityWithoutStress) * maxTranspiration; + + for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) + { + if (! isLayerStressed[i]) + { + double addLayerTransp = redistribution * (roots.rootDensity[unsigned(i)] / totRootDensityWithoutStress); + layerTranspiration[unsigned(i)] += addLayerTransp; + TRs += addLayerTransp; + } + } + } + + waterStress = 1 - (TRs / maxTranspiration); + + for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) + { + totalTranspiration += layerTranspiration[unsigned(i)]; + } + */ + + return soilColumnTranspirationSum; +} + + // input: timeStep [s] void Project3D::computeWaterBalance3D(double totalTimeStep) { @@ -1155,17 +1288,10 @@ bool Project3D::interpolateAndSaveHourlyMeteo(meteoVariable myVar, const QDateTi } -bool Project3D::computeWaterSinkSource() +void Project3D::assignPrecipitation() { // initialize totalPrecipitation = 0; - totalEvaporation = 0; - totalTranspiration = 0; - - for (unsigned long i = 0; i < nrNodes; i++) - { - waterSinkSource.at(size_t(i)) = 0.; - } double area = DEM.header->cellSize * DEM.header->cellSize; @@ -1191,8 +1317,32 @@ bool Project3D::computeWaterSinkSource() } } } +} + + +bool Project3D::setSinkSource() +{ + for (unsigned long i = 0; i < nrNodes; i++) + { + int myResult = soilFluxes3D::setWaterSinkSource(signed(i), waterSinkSource.at(i)); + if (isCrit3dError(myResult, errorString)) + { + errorString = "waterBalanceSinkSource:" + errorString; + return false; + } + } + + return true; +} + + +void Project3D::computeBareSoilEvaporation() +{ + // initialize + totalEvaporation = 0; + + double area = DEM.header->cellSize * DEM.header->cellSize; - // evaporation for (int row = 0; row < indexMap.at(0).header->nrRows; row++) { for (int col = 0; col < indexMap.at(0).header->nrCols; col++) @@ -1200,19 +1350,21 @@ bool Project3D::computeWaterSinkSource() int surfaceIndex = indexMap.at(0).value[row][col]; if (surfaceIndex != indexMap.at(0).header->flag) { - // TODO crop double lai = 0; - double realEvap = computeEvaporation(row, col, lai); // [mm] + double realEvap = assignEvaporation(row, col, lai); // [mm] double flow = area * (realEvap / 1000.); // [m3 h-1] totalEvaporation += flow; } } } +} - // crop transpiration - // TODO - /* + +// crop transpiration +// TODO +/* + totalTranspiration = 0; for (unsigned int layerIndex=1; layerIndex < nrLayers; layerIndex++) { for (long row = 0; row < indexMap.at(size_t(layerIndex)).header->nrRows; row++) @@ -1235,19 +1387,10 @@ bool Project3D::computeWaterSinkSource() } */ - // assign sink/source - for (unsigned long i = 0; i < nrNodes; i++) - { - int myResult = soilFluxes3D::setWaterSinkSource(signed(i), waterSinkSource.at(i)); - if (isCrit3dError(myResult, errorString)) - { - errorString = "waterBalanceSinkSource:" + errorString; - return false; - } - } - return true; -} + + + bool Project3D::setCriteria3DMap(criteria3DVariable var, int layerIndex) diff --git a/bin/CRITERIA3D/shared/project3D.h b/bin/CRITERIA3D/shared/project3D.h index 2a3427faa..abc594589 100644 --- a/bin/CRITERIA3D/shared/project3D.h +++ b/bin/CRITERIA3D/shared/project3D.h @@ -142,8 +142,11 @@ const QString& dailyPath, const QString& hourlyPath); bool interpolateHourlyMeteoVar(meteoVariable myVar, const QDateTime& myTime); - double computeEvaporation(int row, int col, double lai); - bool computeWaterSinkSource(); + double assignEvaporation(int row, int col, double lai); + double assignTranspiration(int row, int col, double lai); + void computeBareSoilEvaporation(); + void assignPrecipitation(); + bool setSinkSource(); void computeWaterBalance3D(double totalTimeStep); bool updateCrop(QDateTime myTime);