Commit 51471f70 authored by jouke's avatar jouke
Browse files

Several fixes related to leaf temperature, transpiration 

Several code improvements, removing redundand code,
reducing scope, makine sure all is initiated.
parent 43181757

OpenSimRoot/InputFiles/templates/plantTemplate.C3C4Photosynthesis.xml

+23 −20
Original line number Diff line number Diff line
@@ -33,44 +33,44 @@ You should have received the GNU GENERAL PUBLIC LICENSE v3 with this file in lic 
		<SimulaVariable name="photosynthesis" function="sumOfSunlitAndShaded" unit="g" replacesPreviousDeclaration="true"/>

		<SimulaVariable name="sunlitPhotosynthesis" function="integratePhotosynthesisRate" unit="g"/>

		<SimulaVariable name="shadedPhotosynthesis" function="integratePhotosynthesisRate" unit="g"/>

		<SimulaVariable name="photosynthesisPerLeafArea" function="weighedAverageOfSunlitAndShaded" unit="umol/m2" replacesPreviousDeclaration="true"/>

		<SimulaVariable name="photosynthesisPerLeafArea" function="weightedAverageOfSunlitAndShaded" unit="umol/m2" replacesPreviousDeclaration="true"/>

		<SimulaVariable name="sunlitPhotosynthesisPerLeafArea" function="photosynthesisRateFarquhar" unit="umol/m2" integrationFunction="explicitConvergence"/>

		<SimulaVariable name="shadedPhotosynthesisPerLeafArea" function="photosynthesisRateFarquhar" unit="umol/m2" integrationFunction="explicitConvergence"/>

		<SimulaDerivative name="carbonLimitedPhotosynthesisRate" function="weighedAverageOfSunlitAndShaded" unit="umol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="carbonLimitedPhotosynthesisRate" function="weightedAverageOfSunlitAndShaded" unit="umol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitCarbonLimitedPhotosynthesisRate" function="carbonLimitedPhotosynthesisRate" unit="umol/m2/s"/>

		<SimulaDerivative name="shadedCarbonLimitedPhotosynthesisRate" function="carbonLimitedPhotosynthesisRate" unit="umol/m2/s"/>

		<SimulaDerivative name="lightLimitedPhotosynthesisRate" function="weighedAverageOfSunlitAndShaded" unit="umol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="lightLimitedPhotosynthesisRate" function="weightedAverageOfSunlitAndShaded" unit="umol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitLightLimitedPhotosynthesisRate" function="lightLimitedPhotosynthesisRate" unit="umol/m2/s"/>

		<SimulaDerivative name="shadedLightLimitedPhotosynthesisRate" function="lightLimitedPhotosynthesisRate" unit="umol/m2/s"/>

		<SimulaDerivative name="mesophyllCO2Concentration" function="weighedAverageOfSunlitAndShaded" unit="umol/mol" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="mesophyllCO2Concentration" function="weightedAverageOfSunlitAndShaded" unit="umol/mol" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitMesophyllCO2Concentration" function="mesophyllCO2Concentration" unit="umol/mol"/>

		<SimulaDerivative name="shadedMesophyllCO2Concentration" function="mesophyllCO2Concentration" unit="umol/mol"/>

		<SimulaDerivative name="PEPCarboxylationRate" function="weighedAverageOfSunlitAndShaded" unit="umol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="PEPCarboxylationRate" function="weightedAverageOfSunlitAndShaded" unit="umol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitPEPCarboxylationRate" function="PEPCarboxylationRate" unit="umol/m2/s"/>

		<SimulaDerivative name="shadedPEPCarboxylationRate" function="PEPCarboxylationRate" unit="umol/m2/s"/>

		<SimulaDerivative name="bundleSheathCO2Concentration" function="weighedAverageOfSunlitAndShaded" unit="umol/mol" replacesPreviousDeclaration="true"/><!-- Only for C4 -->

		<SimulaDerivative name="bundleSheathCO2Concentration" function="weightedAverageOfSunlitAndShaded" unit="umol/mol" replacesPreviousDeclaration="true"/><!-- Only for C4 -->

		<SimulaDerivative name="sunlitBundleSheathCO2Concentration" function="bundleSheathCO2Concentration" unit="umol/mol"/><!-- Only for C4 -->

		<SimulaDerivative name="shadedBundleSheathCO2Concentration" function="bundleSheathCO2Concentration" unit="umol/mol"/><!-- Only for C4 -->

		<SimulaDerivative name="CO2Leakage" function="weighedAverageOfSunlitAndShaded" unit="umol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="CO2Leakage" function="weightedAverageOfSunlitAndShaded" unit="umol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitCO2Leakage" function="CO2Leakage" unit="umol/m2/s"/>

		<SimulaDerivative name="shadedCO2Leakage" function="CO2Leakage" unit="umol/m2/s"/>

		<SimulaDerivative name="mesophyllO2Concentration" function="weighedAverageOfSunlitAndShaded" unit="mmol/mol" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="mesophyllO2Concentration" function="weightedAverageOfSunlitAndShaded" unit="mmol/mol" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitMesophyllO2Concentration" function="mesophyllO2Concentration" unit="mmol/mol"/>

		<SimulaDerivative name="shadedMesophyllO2Concentration" function="mesophyllO2Concentration" unit="mmol/mol"/>

		<SimulaDerivative name="bundleSheathO2Concentration" function="weighedAverageOfSunlitAndShaded" unit="mmol/mol" replacesPreviousDeclaration="true"/><!-- Only for C4 -->

		<SimulaDerivative name="bundleSheathO2Concentration" function="weightedAverageOfSunlitAndShaded" unit="mmol/mol" replacesPreviousDeclaration="true"/><!-- Only for C4 -->

		<SimulaDerivative name="sunlitBundleSheathO2Concentration" function="bundleSheathO2Concentration" unit="mmol/mol"/><!-- Only for C4 -->

		<SimulaDerivative name="shadedBundleSheathO2Concentration" function="bundleSheathO2Concentration" unit="mmol/mol"/><!-- Only for C4 -->

		<SimulaDerivative name="O2Leakage" function="weighedAverageOfSunlitAndShaded" unit="mmol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="O2Leakage" function="weightedAverageOfSunlitAndShaded" unit="mmol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitO2Leakage" function="O2Leakage" unit="mmol/m2/s"/>

		<SimulaDerivative name="shadedO2Leakage" function="O2Leakage" unit="mmol/m2/s"/>

		<SimulaDerivative name="stomatalConductance" function="weighedAverageOfSunlitAndShaded" unit="mol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="stomatalConductance" function="weightedAverageOfSunlitAndShaded" unit="mol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitStomatalConductance" function="stomatalConductance" unit="mol/m2/s"/>

		<SimulaDerivative name="shadedStomatalConductance" function="stomatalConductance" unit="mol/m2/s"/>

		<!--SimulaConstant name="leafTemperature" unit="degreesC"> 20 </SimulaConstant>

		<SimulaConstant name="sunlitLeafTemperature" unit="degreesC"> 20 </SimulaConstant>

		<SimulaConstant name="shadedLeafTemperature" unit="degreesC"> 20 </SimulaConstant-->



		<SimulaDerivative name="leafTemperature" function="weighedAverageOfSunlitAndShaded" unit="degreesC" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="leafTemperature" function="weightedAverageOfSunlitAndShaded" unit="degreesC" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitLeafTemperature" function="leafTemperature" unit="degreesC"/>

		<SimulaDerivative name="shadedLeafTemperature" function="leafTemperature" unit="degreesC"/>
@@ -81,14 +81,14 @@ You should have received the GNU GENERAL PUBLIC LICENSE v3 with this file in lic 
		

		<SimulaVariable name="cropHeight" function="cropHeight" unit="cm"/>

		<SimulaVariable name="potentialCropHeight" function="cropHeight" unit="cm"/>

		<SimulaVariable name="sunlitLeafArea" unit="cm2" function="sunlitLeafArea"/>

		<SimulaVariable name="shadedLeafArea" unit="cm2" function="sunlitLeafArea"/>

		<SimulaVariable name="sunlitLeafAreaIndex" unit="cm2" function="sunlitLeafAreaIndex"/>

		<SimulaVariable name="shadedLeafAreaIndex" unit="cm2" function="shadedLeafAreaIndex"/>

		<SimulaDerivative name="sunlitLeafArea" unit="cm2" function="sunlitLeafArea"/>

		<SimulaDerivative name="shadedLeafArea" unit="cm2" function="sunlitLeafArea"/>

		<SimulaDerivative name="sunlitLeafAreaIndex" unit="cm2/cm2" function="sunlitLeafAreaIndex"/>

		<SimulaDerivative name="shadedLeafAreaIndex" unit="cm2/cm2" function="shadedLeafAreaIndex"/>

		<SimulaVariable name="leafRespiration" unit="g" function="useDerivative" replacesPreviousDeclaration="true"/>

		<SimulaVariable name="shadedLeafRespiration" unit="g" function="useDerivative"/>

		<SimulaVariable name="sunlitLeafRespiration" unit="g" function="useDerivative"/>

		<SimulaDerivative name="leafRespirationRate" unit="g/day" function="weighedAverageOfSunlitAndShaded" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="leafRespirationRate" unit="g/day" function="weightedAverageOfSunlitAndShaded" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitLeafRespirationRate" unit="g/day" function="leafRespirationRateFarquhar" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="shadedLeafRespirationRate" unit="g/day" function="leafRespirationRateFarquhar"replacesPreviousDeclaration="true"/>

		<SimulaVariable name="potentialTranspiration" unit="cm3" function="useDerivative" replacesPreviousDeclaration="true"/>
@@ -105,10 +105,10 @@ You should have received the GNU GENERAL PUBLIC LICENSE v3 with this file in lic 
		<SimulaDerivative name="meanShadedLeafAreaIndex" function="meanLeafAreaIndex" unit="cm2/cm2"/>

		<SimulaDerivative name="sunlitLeafArea" function="sumOverAllPlantShoots" unit="cm2"/>

		<SimulaDerivative name="shadedLeafArea" function="sumOverAllPlantShoots" unit="cm2"/>

		<SimulaDerivative name="stomatalResistance" function="weighedAverageOfSunlitAndShaded" unit="s/m" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="stomatalResistance" function="weightedAverageOfSunlitAndShaded" unit="s/m" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitStomatalResistance" function="stomatalResistance" unit="s/m"/>

		<SimulaDerivative name="shadedStomatalResistance" function="stomatalResistance" unit="s/m"/>

		<SimulaDerivative name="meanStomatalConductance" function="weighedAverageOfSunlitAndShaded" unit="mol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="meanStomatalConductance" function="weightedAverageOfSunlitAndShaded" unit="mol/m2/s" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="meanSunlitStomatalConductance" function="meanStomatalConductance" unit="mol/m2/s"/>

		<SimulaDerivative name="meanShadedStomatalConductance" function="meanStomatalConductance" unit="mol/m2/s"/>

		<SimulaDerivative name="canopyHeight" function="maximumCanopyHeight" unit="cm" replacesPreviousDeclaration="true"/>
@@ -117,7 +117,10 @@ You should have received the GNU GENERAL PUBLIC LICENSE v3 with this file in lic 
	<SimulaDirective path="/environment/atmosphere">

		<SimulaDerivative name="sineSolarElevationAngle" function="sineSolarElevationAngle" unit="cm" replacesPreviousDeclaration="true"/>

		<!--SimulaDerivative name="solarElevationAngle" function="solarElevationAngle"/-->

			<SimulaDerivative name="irradiation" function="diurnalRadiationSimulator" unit="W/m2"/>

			<!--SimulaDerivative name="shortWaveRadiation" function="diurnalRadiationSimulator" unit="W/m2" replacesPreviousDeclaration="true"/-->

			<!-- just to make sure these are not used -->

			<SimulaConstant name="shortwaveRadiation" unit="W/m2" replacesPreviousDeclaration="true">-999</SimulaConstant>

			<SimulaConstant name="irradiation" unit="W/m2"replacesPreviousDeclaration="true">-999</SimulaConstant>

			<SimulaDerivative name="beamIrradiation" function="beamRadiationSimulator" unit="W/m2"/>

			<SimulaDerivative name="diffuseIrradiation" function="diffuseRadiationSimulator" unit="W/m2"/>

			<SimulaDerivative name="dayLightHours" function="photoPeriod" unit="hour"/>

OpenSimRoot/src/modules/Atmosphere/DiffusionResistance.cpp

+6 −1
Original line number Diff line number Diff line
@@ -73,6 +73,9 @@ void StomatalResistance::calculate(const Time &t, double& r_s) { 
		if (airTemperature_) airTemperature_->get(t, temperature);

		if (pressure_) pressure_->get(t, pressure);

		r_s = pressure/(g_sV*8.314*(temperature + 273.15)); // Converting to s/m using the ideal gas law

		if(std::isnan(r_s)){

			msg::error("StomatalResistance: value is NaN");

		}

	} else {

//		Parameter rl is the inverse of the stomatal conductance per unit leaf area.

//		Several earlier studies have fixed the value of rl at 100 s m−1 for well-watered
@@ -121,7 +124,9 @@ void StomatalResistance::calculate(const Time &t, double& r_s) { 
		// It approximates to 45 s/m for a 0.50-m crop height.

		// The hourly r_c (=r_s) values of 50 (clipped grass) or 30 (alfalfa) and 200 s/m (both crops), for day and nighttime respectively,

		// were concluded to be fairly accurate in matching reference evapotranspiration calculated with daily data (Walter et al., 2002).



		if(std::isnan(r_s)){

			msg::error("StomatalResistance: value is NaN");

		}

	}

}

OpenSimRoot/src/modules/Atmosphere/ETbase.cpp

+27 −12
Original line number Diff line number Diff line
@@ -180,24 +180,28 @@ void ETbaseclass::calculate(const Time &t, double& ET){ 
	calculateET(t, ETCROP, ETSOIL );







	//note that this split is afterwards. Theoretically we could split radiation energy, but the sum would not be 1 and the equation are not made for that.

	double KDF,LAI;

	extinctionCoef_->get(t,KDF);

	leafAreaIndex_->get(t,LAI);

	const double split=exp(-KDF *LAI);

	double split=exp(-KDF *LAI);

	ETSOIL *= (split);

	ETCROP *= (1-split);



	if (splitBySunStatus){

		if(LAI>0.5){

			double LAIfrac;

			splitLAI_->get(t,LAIfrac);

			double frac=exp(-KDF *LAIfrac);

			if(splitBySunStatus==1) frac=1-frac;

			ETCROP *= LAIfrac/LAI;

		}else{

			if(splitBySunStatus==2) ETCROP=0.;

	split = 1-split; //crop fraction



	if (splitBySunStatus && LAI>0.){

		double LAIfrac(0);

		splitLAI_->get(t,LAIfrac);//sunlit LAI

		if(LAIfrac<0)LAIfrac=0;

		if(LAIfrac>LAI)LAIfrac=LAI;

		double frac=1-exp(-KDF *LAIfrac);//frac of light captured by sunlit LAI

		if(splitBySunStatus==2) frac=split-frac;//fraction of light captured by shaded LAI

		ETCROP *= frac;

		if(ETCROP<0){

			ETCROP=0;

		}

	}else{

		ETCROP *= split;

	}



	switch (mode) {
@@ -217,6 +221,17 @@ void ETbaseclass::calculate(const Time &t, double& ET){ 
			}

			break;

	}

	if(std::isnan(ET))

			msg::error("ETbaseclass::calculate: values are NaN ");



	if(ETCROP<0){

		ETCROP=0;

	}

	if(ETSOIL<0){

		ETSOIL=0;

	}





}



double ETbaseclass::calculateLatentHeat(double &temperature) { //j/kg

OpenSimRoot/src/modules/Atmosphere/EvapoEquations.cpp

+12 −0
Original line number Diff line number Diff line
@@ -92,6 +92,18 @@ DerivativeBase * newInstantiationPenmanMonteith( 
    /// Monteith, J.L., Evaporation and Environment, Proc. Symp. Soc. Exp. Biol. 19 (1965), 205-234

	ETCROP = (delta*H_crop + conversion_day*airDensity*airHeatCapacity*VPD/r_a)/(lambda*(delta + gamma*(1 + r_s/r_a))); // kg/m2/d

	ETCROP = ETCROP*0.1; // convert from kg/m2/d to cm3/cm2/d



	if(ETCROP<0){

		ETCROP=0;

	}

	if(ETSOIL<0){

		ETSOIL=0;

	}



	if(std::isnan(ETCROP) || std::isnan(ETSOIL)){

		stomatalResistance_->get(t, r_s);

		msg::error("PenmanMonteith::calculateET: values are NaN ");

	}

}

OpenSimRoot/src/modules/Atmosphere/Transpiration.cpp

+0 −3
Original line number Diff line number Diff line
@@ -338,9 +338,6 @@ void StomatalConductance::calculate(const Time &t,double &sc){ 
	// Unit of sc = mol/(m^2*s)

}



void StomatalConductance::getDefaultValue(const Time &t, double &var){

	var = residualConductance;

}

std::string StomatalConductance::getName()const{

	return "stomatalConductance";

}