Merge commit '4c1c69c980601c9420f963c0db8e160ac4835162'

agrolib/hydrall/hydrall.cpp

@@ -698,13 +698,16 @@ void Crit3D_Hydrall::cumulatedResults()
 {
     // taken from Hydrall Model, Magnani UNIBO
     // Cumulate hourly values of gas exchange
-    deltaTime.absorbedPAR = simulationStepInSeconds*(sunlit.absorbedPAR+shaded.absorbedPAR);  //absorbed PAR (mol m-2 yr-1)
+    deltaTime.absorbedPAR = simulationStepInSeconds*(sunlit.absorbedPAR+shaded.absorbedPAR);  //absorbed PAR (mol m-2)
     deltaTime.grossAssimilation = simulationStepInSeconds * treeAssimilationRate ; // canopy gross assimilation (mol m-2)
     deltaTime.respiration = simulationStepInSeconds * Crit3D_Hydrall::plantRespiration() ;
     deltaTime.netAssimilation = deltaTime.grossAssimilation- deltaTime.respiration ;
-    deltaTime.netAssimilation = deltaTime.netAssimilation*12/1000.0/CARBONFACTOR ;
+    deltaTime.netAssimilation = deltaTime.netAssimilation*12/1000.0 ; //CARBONFACTOR DA METTERE dopo convert to kg DM m-2
+
+    statePlant.treeNetPrimaryProduction += deltaTime.netAssimilation ;
+
+    //understorey
 
-    //statePlant.stateGrowth.cumulatedBiomass += deltaTime.netAssimilation ; TODO stand biomass
 
     deltaTime.transpiration = 0.;
 
@@ -812,11 +815,170 @@ bool Crit3D_Hydrall::growthStand()
     // understorey update
     statePlant.understoreycumulatedBiomassFoliage = statePlant.understoreycumulatedBiomass * (1.-understoreyAllocationCoefficientToRoot);    //understorey growth: foliage...
     statePlant.understoreycumulatedBiomassRoot = statePlant.understoreycumulatedBiomass * understoreyAllocationCoefficientToRoot;         //...and roots
+
     // canopy update
     statePlant.treecumulatedBiomassFoliage -= (statePlant.treecumulatedBiomassFoliage/plant.foliageLongevity);
     statePlant.treecumulatedBiomassSapwood -= (statePlant.treecumulatedBiomassSapwood/plant.sapwoodLongevity);
     statePlant.treecumulatedBiomassRoot -= (statePlant.treecumulatedBiomassRoot/plant.fineRootLongevity);
 
 
+    double store;
+    //annual stand growth
+    if (isFirstYearSimulation)
+    {
+        annualGrossStandGrowth = statePlant.treeNetPrimaryProduction / CARBONFACTOR; //kg DM m-2
+        store = 0;
+    }
+    else
+    {
+        annualGrossStandGrowth = (store + statePlant.treeNetPrimaryProduction) / 2 / CARBONFACTOR;
+        store = (store + statePlant.treeNetPrimaryProduction) / 2;
+    }
+
+    if (isFirstYearSimulation)
+    {
+        //optimal
+    }
+    else
+    {
+        double allocationCoeffientFoliageOld = allocationCoefficient.toFoliage;
+        double allocationCoeffientFineRootsOld = allocationCoefficient.toFineRoots;
+        double allocationCoeffientSapwoodOld = allocationCoefficient.toSapwood;
+
+        //optimal
+
+        allocationCoefficient.toFoliage = (allocationCoeffientFoliageOld + allocationCoefficient.toFoliage) / 2;
+        allocationCoefficient.toFineRoots = (allocationCoeffientFineRootsOld + allocationCoefficient.toFineRoots) / 2;
+        allocationCoefficient.toSapwood = (allocationCoeffientSapwoodOld + allocationCoefficient.toSapwood) / 2;
+    }
+
+    if (annualGrossStandGrowth * allocationCoefficient.toFoliage > statePlant.treecumulatedBiomassFoliage/(plant.foliageLongevity - 1))
+
+    isFirstYearSimulation = false;
     return true;
 }
+
+
+void Crit3D_Hydrall::resetStandVariables()
+{
+    statePlant.treeNetPrimaryProduction = 0;
+}
+
+void Crit3D_Hydrall::optimal()
+{
+    double allocationCoefficientFoliageOld;
+    double increment;
+    double incrementStart = 5e-2;
+    bool sol = 0;
+    double allocationCoefficientFoliage0;
+    double bisectionMethodIntervalALLF;
+    int jmax = 40;
+    double accuracy = 1e-3;
+
+    for (int j = 0; j < 3; j++)
+    {
+        allocationCoefficientFoliageOld = 1;
+        increment = incrementStart / std::pow(10, j);
+        allocationCoefficient.toFoliage = 1;
+
+        while(! sol && allocationCoefficient.toFoliage > -EPSILON)
+        {
+            rootfind(allocationCoefficient.toFoliage, allocationCoefficient.toFineRoots, allocationCoefficient.toSapwood, sol);
+
+            if (sol)
+                break;
+
+            allocationCoefficientFoliageOld = allocationCoefficient.toFoliage;
+            allocationCoefficient.toFoliage -= increment;
+
+        }
+        if (sol)
+            break;
+    }
+
+    if (sol)
+    {
+        //find optimal allocation coefficients by bisection technique
+
+        double allocationCoefficientFoliageMid,  allocationCoefficientFineRootsMid, allocationCoefficientSapwoodMid;
+        bool solmid = 0;
+
+        //set starting point and range
+        allocationCoefficientFoliage0 = allocationCoefficient.toFoliage;
+        bisectionMethodIntervalALLF = allocationCoefficientFoliageOld - allocationCoefficient.toFoliage;
+
+        //bisection loop
+        int j = 0;
+        while (std::abs(bisectionMethodIntervalALLF) > accuracy && j < jmax)
+        {
+            bisectionMethodIntervalALLF /= 2;
+            allocationCoefficientFoliageMid = allocationCoefficientFoliage0 + bisectionMethodIntervalALLF;
+
+            rootfind(allocationCoefficientFoliageMid, allocationCoefficientFineRootsMid, allocationCoefficientSapwoodMid, solmid);
+
+            if (solmid)
+            {
+                allocationCoefficientFoliage0 = allocationCoefficientFoliageMid;
+                allocationCoefficient.toFoliage = allocationCoefficientFoliageMid;
+                allocationCoefficient.toFineRoots = allocationCoefficientFineRootsMid;
+                allocationCoefficient.toSapwood = allocationCoefficientSapwoodMid;
+            }
+
+            j++;
+        }
+    }
+}
+
+void Crit3D_Hydrall::rootfind(double &allf, double &allr, double &alls, bool &sol)
+{
+    //search for a solution to hydraulic constraint
+    double foliageBiomassNew, heightNew;
+
+    //new foliage biomass of tree after growth
+    if (allf < 0) allf = 0;
+
+    foliageBiomassNew = statePlant.treecumulatedBiomassFoliage + (allf*annualGrossStandGrowth);
+
+    //new tree height after growth
+    if (allf*annualGrossStandGrowth > statePlant.treecumulatedBiomassFoliage/(plant.foliageLongevity-1)) {
+        heightNew = plant.height + (allf*annualGrossStandGrowth-statePlant.treecumulatedBiomassFoliage/(plant.foliageLongevity-1)/plant.foliageDensity);
+    } else {
+        heightNew = plant.height;
+    }
+
+    //soil hydraulic conductivity
+    double ksl;
+
+    //specific hydraulic conductivity of soil+roots
+    double soilRootsSpecificConductivity = 1/(1/KR + 1/ksl);
+    double dum = 0.5151 + 0.0242*soil.temperature;
+    if (dum < 0.5151) dum = 0.5151;
+    soilRootsSpecificConductivity *= dum; //adjust for temp effects on water viscosity
+
+    //new sapwood specific conductivity
+    double sapwoodSpecificConductivity = KSMAX * (1-std::exp(-0.69315*heightNew/H50)); //adjust for height effects
+    dum = 0.5151 + 0.0242*weatherVariable.meanDailyTemp;
+    if (dum < 0.5151) dum = 0.5151;
+    sapwoodSpecificConductivity *= dum;     //adjust for temp effects on water viscosity
+
+    //optimal coefficient of allocation to fine roots and sapwood for set allocation to foliage
+    double quadraticEqCoefficient = std::sqrt(soilRootsSpecificConductivity/sapwoodSpecificConductivity*plant.sapwoodLongevity/plant.fineRootLongevity*RHOS); //TODO: check if rhos or ros
+    allr = (statePlant.treecumulatedBiomassSapwood - quadraticEqCoefficient*heightNew*statePlant.treecumulatedBiomassRoot +
+            annualGrossStandGrowth*(1-allf))/annualGrossStandGrowth/(1+quadraticEqCoefficient*heightNew);
+
+    if (allr < EPSILON) allr = EPSILON; //bracket ALLR between (1-ALLF) and a small value
+    if (allr > (1-allf)) allr = 1-allf;
+
+    alls = 1 - allf - allr;
+    if (alls < EPSILON) alls = EPSILON; //bracket ALLS between 1 and a small value
+    if (alls > 1) alls = 1;
+
+    //resulting fine root and sapwood biomass
+
+
+    //resulting leaf specific resistance (MPa s m2 m-3)
+
+    //resulting minimum leaf water potential
+
+    //check if given value of ALLF satisfies optimality constraint
+}

agrolib/hydrall/hydrall.h

@@ -83,7 +83,7 @@
 
     struct TstatePlant
     {
-        double treecumulatedBiomass;
+        double treeNetPrimaryProduction;
         double treecumulatedBiomassFoliage;
         double treecumulatedBiomassRoot;
         double treecumulatedBiomassSapwood;
@@ -117,6 +117,7 @@
         //double meanDailyTemperature;
         double vaporPressureDeficit;
         double last30DaysTAvg;
+        double meanDailyTemp;
 
 
     };
@@ -136,6 +137,7 @@
         double foliageLongevity;
         double sapwoodLongevity;
         double fineRootLongevity;
+        double foliageDensity;
 
 
     };
@@ -218,6 +220,13 @@
         double fineRoot ;
     };
 
+    struct TallocationCoefficient {
+
+        double toFoliage;
+        double toFineRoots;
+        double toSapwood;
+    };
+
 
     class Crit3DHydrallMaps
     {
@@ -238,11 +247,13 @@
     class Crit3D_Hydrall{
     public:
 
-        // Crit3D_Hydrall();
-        // ~Crit3D_Hydrall();
+        //Crit3D_Hydrall();
+        //~Crit3D_Hydrall();
 
         void initialize();
-        bool writeHydrallMaps;
+        bool firstDayOfMonth;
+        int firstMonthVegetativeSeason;
+        bool isFirstYearSimulation;
 
         TbigLeaf sunlit,shaded, understorey;
         TweatherVariable weatherVariable;
@@ -256,12 +267,15 @@
         ThydrallNitrogen nitrogenContent;
         ThydrallBiomass treeBiomass, understoreyBiomass;
         TstatePlant statePlant;
+        TallocationCoefficient allocationCoefficient;
 
 
         double elevation;
         int simulationStepInSeconds;
         double leafAreaIndex;
 
+        double annualGrossStandGrowth;
+
         //gasflux results
         std::vector<double> treeTranspirationRate;          //molH2O m^-2 s^-1
         double treeAssimilationRate;
@@ -294,6 +308,9 @@
         double soilTemperatureModel();
         double temperatureMoistureFunction(double temperature);
         bool growthStand();
+        void resetStandVariables();
+        void optimal();
+        void rootfind(double &allf, double &allr, double &alls, bool &sol);
 
     };
 

agrolib/proxyWidget/localProxyWidget.cpp

@@ -20,7 +20,12 @@ Crit3DLocalProxyWidget::Crit3DLocalProxyWidget(double x, double y, double zDEM,
     localUtmPoint.y = y;
     gis::getLatLonFromUtm(gisSettings, localUtmPoint, localGeoPoint);
 
-    this->setWindowTitle("Local proxy analysis for point of coordinates (" + QString::number(localGeoPoint.latitude) + ", " + QString::number(localGeoPoint.longitude) + ")." + " z value: " + QString::number(zDEM) + " (DEM), macro area nr. " + QString::number(gis::getValueFromXY(*interpolationSettings->getMacroAreasMap(), x, y))); // + QString::number(zGrid) + " (Grid)");
+    QString windowTitle = "Local proxy analysis for point of coordinates (" + QString::number(localGeoPoint.latitude) + ", " + QString::number(localGeoPoint.longitude) + ")." + " z value: " + QString::number(zDEM) + " (DEM)";
+
+    if (interpolationSettings->getUseGlocalDetrending())
+        windowTitle += "macro area nr. " + QString::number(gis::getValueFromXY(*interpolationSettings->getMacroAreasMap(), x, y));
+
+    this->setWindowTitle(windowTitle); // + QString::number(zGrid) + " (Grid)");
     this->resize(1024, 700);
     this->setSizePolicy(QSizePolicy::Preferred, QSizePolicy::Preferred);
     this->setAttribute(Qt::WA_DeleteOnClose);