Skip to content

Commit

Permalink
Merge commit '34e0fd1c97e2ae10d0424b4935274ec5f2ff39c7'
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Jan 31, 2025
2 parents c7c1b09 + 34e0fd1 commit a413700
Show file tree
Hide file tree
Showing 2 changed files with 139 additions and 37 deletions.
125 changes: 98 additions & 27 deletions hydrall/hydrall.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
/*!
\name hydrall.cpp
\brief
\brief
\authors Antonio Volta, Caterina Toscano
*/


Expand All @@ -14,15 +14,23 @@
#include "furtherMathFunctions.h"
#include "physics.h"

Crit3DHydrallMaps::Crit3DHydrallMaps(const gis::Crit3DRasterGrid& DEM)
Crit3DHydrallMaps::Crit3DHydrallMaps()
{
mapLAI = new gis::Crit3DRasterGrid;
mapLast30DaysTavg = new gis::Crit3DRasterGrid;
}

void Crit3DHydrallMaps::initialize(const gis::Crit3DRasterGrid& DEM)
{
mapLAI->initializeGrid(DEM);
mapLast30DaysTavg->initializeGrid(DEM);
}

Crit3DHydrallMaps::~Crit3DHydrallMaps()
{
mapLAI->clear();
mapLast30DaysTavg->clear();
}

bool computeHydrallPoint(Crit3DDate myDate, double myTemperature, double myElevation, int secondPerStep)
{
Expand All @@ -47,11 +55,11 @@ bool computeHydrallPoint(Crit3DDate myDate, double myTemperature, double myEleva

double getCO2(Crit3DDate myDate, double myTemperature, double myElevation)
{
double atmCO2 ; //https://www.eea.europa.eu/data-and-maps/daviz/atmospheric-concentration-of-carbon-dioxide-5/download.table
double atmCO2 = 400 ; //https://www.eea.europa.eu/data-and-maps/daviz/atmospheric-concentration-of-carbon-dioxide-5/download.table
double year[24] = {1750,1800,1850,1900,1910,1920,1930,1940,1950,1960,1970,1980,1990,2000,2010,2020,2030,2040,2050,2060,2070,2080,2090,2100};
double valueCO2[24] = {278,283,285,296,300,303,307,310,311,317,325,339,354,369,389,413,443,473,503,530,550,565,570,575};

interpolation::linearInterpolation(double(myDate.year), year, valueCO2, 24);
atmCO2 = interpolation::linearInterpolation(double(myDate.year), year, valueCO2, 24);

// exponential fitting Mauna Loa
/*if (myDate.year < 1990)
Expand Down Expand Up @@ -85,16 +93,87 @@ double photosynthesisAndTranspiration()
return 0;
}

void weatherVariables(double myInstantTemp,double myRelativeHumidity,double myCloudiness,TweatherDerivedVariable weatherDerivedVariable)
void Crit3D_Hydrall::initialize()
{
// taken from Hydrall Model, Magnani UNIBO
weatherDerivedVariable.airVapourPressure = saturationVaporPressure(myInstantTemp)*myRelativeHumidity/100.;
weatherDerivedVariable.emissivitySky = 1.24 * pow((weatherDerivedVariable.airVapourPressure/100.0) / (myInstantTemp+ZEROCELSIUS),(1.0/7.0))*(1 - 0.84*myCloudiness)+ 0.84*myCloudiness;
weatherDerivedVariable.longWaveIrradiance = pow(myInstantTemp+ZEROCELSIUS,4) * weatherDerivedVariable.emissivitySky * STEFAN_BOLTZMANN ;
weatherDerivedVariable.slopeSatVapPressureVSTemp = 2588464.2 / pow(240.97 + myInstantTemp, 2) * exp(17.502 * myInstantTemp / (240.97 + myInstantTemp)) ;
//initializeLeaf(sunlit);
//initializeLeaf(shaded);
//weatherVariable.initialize();
myChlorophyllContent = NODATA;
elevation = NODATA;
}

void hydrall::radiationAbsorption(double mySunElevation, double leafAreaIndex)
void initializeLeaf(TbigLeaf myLeaf)
{
myLeaf.leafAreaIndex = NODATA;
myLeaf.leafTemperature = NODATA;
myLeaf.absorbedPAR = NODATA;
myLeaf.aerodynamicConductanceCO2Exchange = NODATA;
myLeaf.aerodynamicConductanceHeatExchange = NODATA;
myLeaf.totalConductanceHeatExchange = NODATA;
myLeaf.minimalStomatalConductance = NODATA;
myLeaf.assimilation = NODATA;
myLeaf.isothermalNetRadiation = NODATA;
myLeaf.darkRespiration = NODATA;
myLeaf.maximalCarboxylationRate = NODATA;
myLeaf.maximalElectronTrasportRate = NODATA;
myLeaf.carbonMichaelisMentenConstant = NODATA;
myLeaf.oxygenMichaelisMentenConstant = NODATA;
myLeaf.compensationPoint = NODATA;
myLeaf.convexityFactorNonRectangularHyperbola = NODATA;
myLeaf.quantumYieldPS2 = NODATA;
myLeaf.assimilation = NODATA;
myLeaf.transpiration = NODATA;
myLeaf.stomatalConductance = NODATA;

}

void Crit3D_Hydrall::setHourlyVariables(double temp, double irradiance , double prec , double relativeHumidity , double windSpeed, double directIrradiance, double diffuseIrradiance, double cloudIndex)
{
setWeatherVariables(temp, irradiance, prec, relativeHumidity, windSpeed, directIrradiance, diffuseIrradiance, cloudIndex);
}

bool Crit3D_Hydrall::setWeatherVariables(double temp, double irradiance , double prec , double relativeHumidity , double windSpeed, double directIrradiance, double diffuseIrradiance, double cloudIndex)
{

bool isReadingOK = false ;
weatherVariable.irradiance = irradiance ;
weatherVariable.myInstantTemp = temp ;
weatherVariable.prec = prec ;
weatherVariable.relativeHumidity = relativeHumidity ;
weatherVariable.windSpeed = windSpeed ;
weatherVariable.atmosphericPressure = getPressureFromElevation(weatherVariable.myInstantTemp, elevation) ;
//weatherVariable.meanDailyTemperature = meanDailyTemp;
double deltaRelHum = MAXVALUE(100.0 - weatherVariable.relativeHumidity, 0.01);
weatherVariable.vaporPressureDeficit = 0.01 * deltaRelHum * 613.75 * exp(17.502 * weatherVariable.myInstantTemp / (240.97 + weatherVariable.myInstantTemp));

setDerivedWeatherVariables(directIrradiance, diffuseIrradiance, cloudIndex);

if ((int(prec) != NODATA) && (int(temp) != NODATA) && (int(windSpeed) != NODATA)
&& (int(irradiance) != NODATA) && (int(relativeHumidity) != NODATA) && (int(weatherVariable.atmosphericPressure) != NODATA))
isReadingOK = true;

return isReadingOK;
}

void Crit3D_Hydrall::setDerivedWeatherVariables(double directIrradiance, double diffuseIrradiance, double cloudIndex)
{
weatherVariable.derived.airVapourPressure = saturationVaporPressure(weatherVariable.myInstantTemp)*weatherVariable.relativeHumidity/100.;
weatherVariable.derived.slopeSatVapPressureVSTemp = 2588464.2 / pow(240.97 + weatherVariable.myInstantTemp, 2) * exp(17.502 * weatherVariable.myInstantTemp / (240.97 + weatherVariable.myInstantTemp)) ;
weatherVariable.derived.myDirectIrradiance = directIrradiance;
weatherVariable.derived.myDiffuseIrradiance = diffuseIrradiance;
double myCloudiness = MINVALUE(1,MAXVALUE(0,cloudIndex));
weatherVariable.derived.myEmissivitySky = 1.24 * pow((weatherVariable.derived.airVapourPressure/100.0) / (weatherVariable.myInstantTemp+ZEROCELSIUS),(1.0/7.0))*(1 - 0.84*myCloudiness)+ 0.84*myCloudiness;
weatherVariable.derived.myLongWaveIrradiance = pow(weatherVariable.myInstantTemp+ZEROCELSIUS,4) * weatherVariable.derived.myEmissivitySky * STEFAN_BOLTZMANN ;

return;
}

void Crit3D_Hydrall::setPlantVariables(double chlorophyllContent)
{
myChlorophyllContent = chlorophyllContent;
}

void Crit3D_Hydrall::radiationAbsorption(double mySunElevation, double leafAreaIndex)
{
// taken from Hydrall Model, Magnani UNIBO
double sineSolarElevation;
Expand Down Expand Up @@ -133,7 +212,7 @@ void hydrall::radiationAbsorption(double mySunElevation, double leafAreaIndex)
shaded.leafAreaIndex = leafAreaIndex - sunlit.leafAreaIndex ;
//Extinction coefficients for direct and diffuse PAR and NIR radiation, scattering leaves
//Based on approximation by Goudriaan 1977 (in Goudriaan & van Laar 1994)
double exponent= -pow(10,0.28 + 0.63*log10(chlorophyllContent*0.85/1000));
double exponent= -pow(10,0.28 + 0.63*log10(myChlorophyllContent*0.85/1000));
leafAbsorbancePAR = 1 - pow(10,exponent);//from Agusti et al (1994), Eq. 1, assuming Chl a = 0.85 Chl (a+b)
scatteringCoefPAR = 1.0 - leafAbsorbancePAR ; //scattering coefficient for PAR
scatteringCoefNIR = 1.0 - leafAbsorbanceNIR ; //scattering coefficient for NIR
Expand All @@ -150,9 +229,9 @@ void hydrall::radiationAbsorption(double mySunElevation, double leafAreaIndex)
directReflectionCoefficientNIR = directReflectionCoefficientPAR = dum[4] * dum[2];
diffuseReflectionCoefficientNIR = diffuseReflectionCoefficientPAR = dum[4] * dum[3];
//Incoming direct PAR and NIR (W m-2)
directIncomingNIR = directIncomingPAR = myDirectIrradiance * 0.5 ;
directIncomingNIR = directIncomingPAR = weatherVariable.derived.myDirectIrradiance * 0.5 ;
//Incoming diffuse PAR and NIR (W m-2)
diffuseIncomingNIR = diffuseIncomingPAR = myDiffuseIrradiance * 0.5 ;
diffuseIncomingNIR = diffuseIncomingPAR = weatherVariable.derived.myDiffuseIrradiance * 0.5 ;
//Preliminary computations
dum[5]= diffuseIncomingPAR * (1.0-diffuseReflectionCoefficientPAR) * diffuseLightKPAR ;
dum[6]= directIncomingPAR * (1.0-directReflectionCoefficientPAR) * directLightKPAR ;
Expand All @@ -172,12 +251,12 @@ void hydrall::radiationAbsorption(double mySunElevation, double leafAreaIndex)
sunlitAbsorbedNIR = dum[8]*dum[13]+dum[9]*dum[14]+dum[10]*dum[15];
shadedAbsorbedNIR = dum[8]*(UPSCALINGFUNC(diffuseLightKNIR,leafAreaIndex)-dum[13])+dum[9]*(UPSCALINGFUNC(directLightKNIR,leafAreaIndex)- dum[14]) - dum[10] * dum[15];
// Long-wave radiation balance by sunlit (1) and shaded (2) big-leaf (W m-2) from Wang & Leuning 1998
dum[16]= myLongWaveIrradiance -STEFAN_BOLTZMANN*pow(myInstantTemp+ZEROCELSIUS,4); //negativo
dum[16]= weatherVariable.derived.myLongWaveIrradiance -STEFAN_BOLTZMANN*pow(weatherVariable.myInstantTemp+ZEROCELSIUS,4); //negativo
dum[16] *= diffuseLightK ;
double emissivityLeaf, emissivitySoil;
emissivityLeaf = 0.96 ; // supposed constant because variation is very small
emissivitySoil= 0.94 ; // supposed constant because variation is very small
sunlitAbsorbedLW = (dum[16] * UPSCALINGFUNC((directLightK+diffuseLightK),leafAreaIndex))*emissivityLeaf+(1.0-emissivitySoil)*(emissivityLeaf-myEmissivitySky)* UPSCALINGFUNC((2*diffuseLightK),leafAreaIndex)* UPSCALINGFUNC((directLightK-diffuseLightK),leafAreaIndex);
sunlitAbsorbedLW = (dum[16] * UPSCALINGFUNC((directLightK+diffuseLightK),leafAreaIndex))*emissivityLeaf+(1.0-emissivitySoil)*(emissivityLeaf-weatherVariable.derived.myEmissivitySky)* UPSCALINGFUNC((2*diffuseLightK),leafAreaIndex)* UPSCALINGFUNC((directLightK-diffuseLightK),leafAreaIndex);
shadedAbsorbedLW = dum[16] * UPSCALINGFUNC(diffuseLightK,leafAreaIndex) - sunlitAbsorbedLW ;
// Isothermal net radiation for sunlit (1) and shaded (2) big-leaf
sunlit.isothermalNetRadiation= sunlit.absorbedPAR + sunlitAbsorbedNIR + sunlitAbsorbedLW ;
Expand All @@ -196,7 +275,7 @@ void hydrall::radiationAbsorption(double mySunElevation, double leafAreaIndex)
shaded.leafAreaIndex = leafAreaIndex;
shaded.absorbedPAR = 0.0 ;
shadedAbsorbedNIR = 0.0 ;
dum[16]= myLongWaveIrradiance -STEFAN_BOLTZMANN*pow(myInstantTemp + ZEROCELSIUS,4) ;
dum[16]= weatherVariable.derived.myLongWaveIrradiance -STEFAN_BOLTZMANN*pow(weatherVariable.myInstantTemp + ZEROCELSIUS,4) ;
dum[16] *= diffuseLightK ;
shadedAbsorbedLW= dum[16] * (UPSCALINGFUNC(diffuseLightK,leafAreaIndex) - UPSCALINGFUNC(directLightK+diffuseLightK,leafAreaIndex)) ;
shaded.isothermalNetRadiation = shaded.absorbedPAR + shadedAbsorbedNIR + shadedAbsorbedLW ;
Expand All @@ -206,12 +285,4 @@ void hydrall::radiationAbsorption(double mySunElevation, double leafAreaIndex)
sunlit.absorbedPAR *= 4.57E-6 ;
shaded.absorbedPAR *= 4.57E-6 ;
}
/*
double meanLastMonthTemperature(double previousLastMonthTemp, double simulationStepInSeconds, double myInstantTemp)
{
double newTemperature;
double monthFraction;
monthFraction = simulationStepInSeconds/(2592000.0); // seconds of 30 days
newTemperature = previousLastMonthTemp * (1 - monthFraction) + myInstantTemp * monthFraction ;
return newTemperature;
}*/

51 changes: 41 additions & 10 deletions hydrall/hydrall.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,11 +27,31 @@

#define NOT_INITIALIZED_VINE -1


struct TweatherDerivedVariable {
double airVapourPressure;
double emissivitySky;
double longWaveIrradiance;
double slopeSatVapPressureVSTemp;
double myDirectIrradiance;
double myDiffuseIrradiance;
double myEmissivitySky;
double myLongWaveIrradiance;

};

struct TweatherVariable {
TweatherDerivedVariable derived;

double myInstantTemp;
double prec;
double irradiance;
double relativeHumidity;
double windSpeed;
double atmosphericPressure;
//double meanDailyTemperature;
double vaporPressureDeficit;


};

Expand All @@ -51,6 +71,8 @@
double compensationPoint, convexityFactorNonRectangularHyperbola ;
double quantumYieldPS2 ;
double assimilation,transpiration,stomatalConductance ;


};

class Crit3DHydrallMaps
Expand All @@ -61,22 +83,31 @@
gis::Crit3DRasterGrid* mapLAI;
gis::Crit3DRasterGrid* mapLast30DaysTavg;

Crit3DHydrallMaps(const gis::Crit3DRasterGrid& DEM);
Crit3DHydrallMaps();
~Crit3DHydrallMaps();

void clear();
void initialize();
void initialize(const gis::Crit3DRasterGrid& DEM);
};
class hydrall{
class Crit3D_Hydrall{
public:

// Crit3D_Hydrall();
// ~Crit3D_Hydrall();

void initialize();
void initializeLeaf(TbigLeaf myLeaf);
//gis::Crit3DRasterGrid* stateMaps;

TbigLeaf sunlit,shaded;
double myLongWaveIrradiance;
double myInstantTemp;
double myDirectIrradiance;
double myDiffuseIrradiance;
double myEmissivitySky;
double chlorophyllContent;
TweatherVariable weatherVariable;
double myChlorophyllContent;
double elevation;

void radiationAbsorption(double mySunElevation, double leafAreaIndex);
void setHourlyVariables(double temp, double irradiance , double prec , double relativeHumidity , double windSpeed, double directIrradiance, double diffuseIrradiance, double cloudIndex);
bool setWeatherVariables(double temp, double irradiance , double prec , double relativeHumidity , double windSpeed, double directIrradiance, double diffuseIrradiance, double cloudIndex);
void setDerivedWeatherVariables(double directIrradiance, double diffuseIrradiance, double cloudIndex);
void setPlantVariables(double chlorophyllContent);

};

Expand Down

0 comments on commit a413700

Please sign in to comment.