add daily update crop

ARPA-SIMC · Oct 27, 2023 · a5d88ae · a5d88ae
1 parent 9ab8578
commit a5d88ae
Show file tree

Hide file tree

Showing 6 changed files with 61 additions and 5 deletions.
diff --git a/DATA/PROJECT/Montue/DATA/crop_Montue.db b/DATA/PROJECT/Montue/DATA/crop_Montue.db
diff --git a/agrolib/crop/crop.cpp b/agrolib/crop/crop.cpp
@@ -152,11 +152,17 @@ void Crit3DCrop::initialize(double latitude, unsigned int nrLayers, double total
 }
 
 
-double Crit3DCrop::getDailyDegreeIncrease(double tmin, double tmax)
+double Crit3DCrop::getDailyDegreeIncrease(double tmin, double tmax, int doy)
 {
     if (isEqual(tmin, NODATA) || isEqual(tmax, NODATA))
         return NODATA;
 
+    // check crop cycle
+    if (isSowingCrop() && ! isInsideTypicalCycle(doy))
+    {
+        return 0;
+    }
+
     double tmed = (tmin + MINVALUE(tmax, upperThermalThreshold)) * 0.5;
     return MAXVALUE(tmed - thermalThreshold, 0);
 }
@@ -297,7 +303,7 @@ int Crit3DCrop::getDaysFromCurrentSowing(int myDoy) const
 
 bool Crit3DCrop::isInsideTypicalCycle(int myDoy) const
 {
-    return (myDoy >= sowingDoy && getDaysFromTypicalSowing(myDoy) < plantCycle);
+    return (getDaysFromTypicalSowing(myDoy) < plantCycle);
 }
 
 

diff --git a/agrolib/crop/crop.h b/agrolib/crop/crop.h
@@ -81,7 +81,7 @@
         bool isSowingCrop() const;
         bool isRootStatic() const;
 
-        double getDailyDegreeIncrease(double tmin, double tmax);
+        double getDailyDegreeIncrease(double tmin, double tmax, int doy);
 
         void initialize(double latitude, unsigned int nrLayers, double totalSoilDepth, int currentDoy);
         bool needReset(Crit3DDate myDate, double latitude, double waterTableDepth);

diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp
@@ -166,7 +166,7 @@ void Crit3DProject::initializeCrop()
                             float tmax = climateParameters.getClimateVar(dailyAirTemperatureMax, myDate.month,
                                                                          height, quality->getReferenceHeight());
 
-                            double currentDD = cropList[index].getDailyDegreeIncrease(tmin, tmax);
+                            double currentDD = cropList[index].getDailyDegreeIncrease(tmin, tmax, currentDoy);
                             if (! isEqual(currentDD, NODATA))
                             {
                                 degreeDays += currentDD;
@@ -183,6 +183,50 @@ void Crit3DProject::initializeCrop()
 }
 
 
+void Crit3DProject::dailyUpdateCrop()
+{
+    for (int row = 0; row < DEM.header->nrRows; row++)
+    {
+        for (int col = 0; col < DEM.header->nrCols; col++)
+        {
+            float height = DEM.value[row][col];
+            if (! isEqual(height, DEM.header->flag))
+            {
+                // is crop
+                int index = getLandUnitIndexRowCol(row, col);
+                if (index != NODATA && landUnitList[index].idCrop != "")
+                {
+                    float tmin = dailyTminMap.value[row][col];
+                    float tmax = dailyTminMap.value[row][col];
+                    if (! isEqual(tmin, dailyTminMap.header->flag) && ! isEqual(tmax, dailyTmaxMap.header->flag))
+                    {
+                        double dailyDD = cropList[index].getDailyDegreeIncrease(tmin, tmax, currentDate.dayOfYear());
+                        if (! isEqual(dailyDD, NODATA))
+                        {
+                            if (isEqual(degreeDaysMap.value[row][col], degreeDaysMap.header->flag))
+                            {
+                                degreeDaysMap.value[row][col] = float(dailyDD);
+                            }
+                            else
+                            {
+                                degreeDaysMap.value[row][col] += float(dailyDD);
+                            }
+
+                            laiMap.value[row][col] = cropList[index].computeSimpleLAI(degreeDaysMap.value[row][col],
+                                                            gisSettings.startLocation.latitude, currentDate.dayOfYear());
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // clean daily temp maps
+    dailyTminMap.emptyGrid();
+    dailyTmaxMap.emptyGrid();
+}
+
+
 bool Crit3DProject::runModels(QDateTime firstTime, QDateTime lastTime)
 {
     // initialize meteo
@@ -216,6 +260,11 @@ bool Crit3DProject::runModels(QDateTime firstTime, QDateTime lastTime)
     {
         setCurrentDate(myDate);
 
+        if (processes.computeCrop)
+        {
+            dailyUpdateCrop();
+        }
+
         if (isSaveOutputRaster())
         {
             // create directory for hourly raster output

diff --git a/bin/CRITERIA3D/criteria3DProject.h b/bin/CRITERIA3D/criteria3DProject.h
@@ -70,6 +70,7 @@
 
         bool initializeCriteria3DModel();
         void initializeCrop();
+        void dailyUpdateCrop();
 
         bool runModels(QDateTime firstTime, QDateTime lastTime);
 

diff --git a/bin/CRITERIA3D/mainwindow.cpp b/bin/CRITERIA3D/mainwindow.cpp
@@ -2236,7 +2236,7 @@ void MainWindow::on_actionCriteria3D_run_models_triggered()
     myProject.processes.initialize();
     myProject.processes.computeMeteo = true;
     myProject.processes.computeRadiation = true;
-    myProject.processes.computeWater = true;
+    myProject.processes.computeWater = false;
     myProject.processes.computeCrop = true;
 
     startModels(firstTime, lastTime);