Commit 43181757 authored by jouke's avatar jouke
Browse files

unit conversions and removal of wrong photoperiod use

parent 1add3824

OpenSimRoot/InputFiles/environments/WageningseBovenBuurt/atmosphere.xml

+0 −6
Original line number Diff line number Diff line
@@ -36,12 +36,6 @@ You should have received the GNU GENERAL PUBLIC LICENSE v3 with this file in lic 
			    		   

              2.8330E-3 Mj/cm2/d sunny summers day in the Netherlands: day 190 year 1999 data KNMI ; NOTE: use PAR/RDD ratio of 0.5-->

			</SimulaConstant>

			<SimulaConstant

				name="PAR/RDD"

				unit="100%">

			<!--note="optional conversion factor from RDD to PAR - normally 0.5; only used if irradiation isn't converted yet otherwise set to 1"-->

				1.

			</SimulaConstant>

			<SimulaTable

				name_column2="precipitation"

				name_column1="time"

OpenSimRoot/src/engine/DataDefinitions/Units.cpp

+46 −0
Original line number Diff line number Diff line
@@ -146,6 +146,52 @@ double Unit::getUnitConversionFactor(const Unit &target)const{ 
		if(factor<0) msg::error("Unit:::getUnitConversiontFactor: factor unknown for units "+this->name+" into "+target.name);

		return factor;

}

double Unit::getUnitConversionFactorIrradiation(const Unit &target)const{

	double r(1.);

	if(*this=="W/m2"){

		r=2.114*60*60*24/10000; //umol/cm2/day

	}else if (*this=="W/cm2"){

		r=2.114*60*60*24;

	}else if (*this=="J/m2/day"){

		r=2.114/10000; //umol/cm2/day

	}else if (*this=="MJ/m2/day"){

		r=2.114*100; //umol/cm2/day

	}else if (*this=="mol/m2/day"){

		r=100;

	}else if (*this=="umol/cm2/day"){

		r=1;

	}else if (*this=="umol/m2/s"){

		r=10000/(60*60*24);

	}else if (*this=="umol/cm2/s"){

		r=1/(60*60*24);

	}else{

		msg::error("Unit: Unknown irradiation unit "+this->name);

	}



	if(target=="W/m2"){

		r/=2.114*60*60*24/10000; //umol/cm2/day

	}else if (target=="W/cm2"){

		r/=2.114*60*60*24;

	}else if (target=="J/m2/day"){

		r/=2.114/10000; //umol/cm2/day

	}else if (target=="MJ/m2/day"){

		r/=2.114*100; //umol/cm2/day

	}else if (target=="mol/m2/day"){

		r/=100.;

	}else if (target=="umol/cm2/day"){

		r/=1.;

	}else if (target=="umol/m2/s"){

		r/=(10000./(60*60*24));

	}else if (target=="umol/cm2/s"){

		r/=(1./(60*60*24));

	}else{

		msg::error("Unit: Unknown target irradiation unit "+target.name);

	}

	return r;

}









std::string Unit::formattedName()const{

	std::string r;

OpenSimRoot/src/engine/DataDefinitions/Units.hpp

+1 −0
Original line number Diff line number Diff line
@@ -39,6 +39,7 @@ public: 
	Unit(const Unit &copyThis);

	//unit conversion

	double getUnitConversionFactor(const Unit &target)const;

	double getUnitConversionFactorIrradiation(const Unit &target)const;

	void check(const Unit &check,const std::string &parameterName,const std::string &moduleName="")const;

	std::string name;

	double factor;

OpenSimRoot/src/modules/ShootModule/LeafGasExchange.cpp

+17 −19
Original line number Diff line number Diff line
@@ -172,36 +172,35 @@ DerivativeBase * newInstantiationPEPCarboxylationRate(SimulaDynamic* const pSD){ 


BundleSheathCO2Concentration::BundleSheathCO2Concentration(SimulaDynamic* pSD):DerivativeBase(pSD){

	std::string name = pSD->getName().substr(0, 6); // name = sunlit or shaded

	pIrradiation = ORIGIN->getPath("/environment/atmosphere/irradiation");

	SimulaBase *probe = ORIGIN->getPath("/environment/atmosphere/CO2Concentration", "umol/mol");

	probe->get(atmosphericCO2Concentration);

	pPhotosynthesisRate = pSD->existingSibling(name + "PhotosynthesisPerLeafArea", "umol/m2");

	if (!pPhotosynthesisRate) pPhotosynthesis = pSD->getSibling(name + "Photosynthesis", "g");



	pLeafArea = pSD->getSibling(name + "LeafArea", "cm2");



	pLeafRespirationRate = pSD->existingSibling(name + "LeafRespirationRate", "g/day");

	if (!pLeafRespirationRate) pLeafRespiration = pSD->getSibling(name + "LeafRespiration", "g");



	pLeafTemperature = pSD->getSibling(name + "LeafTemperature", "degreesC");

	pMesophyllC = pSD->getSibling(name + "MesophyllCO2Concentration", "umol/mol");

	pPEPCarboxylation = pSD->getSibling(name + "PEPCarboxylationRate", "umol/m2/s");



	std::string plantType;

	PLANTTYPE(plantType,pSD);

	SimulaBase *shootParameters = GETSHOOTPARAMETERS(plantType);

	probe = shootParameters->getChild("bundleSheathConductanceAt25C", "mol/m2/s");

	probe->get(sheathConductanceAt25C);

	probe = shootParameters->getChild("bundleSheathConductanceActivationEnergy", "J/mol");

	probe->get(sheathConductanceActivationEnergy);

	probe = shootParameters->getChild("bundleSheathConductanceDeactivationEnergy", "J/mol");

	probe->get(sheathConductanceDeactivationEnergy);

	probe = shootParameters->getChild("bundleSheathConductanceEntropyTerm", "J/K/mol");

	probe->get(sheathConductanceEntropyTerm);

	probe = shootParameters->getChild("dayRespirationMesophyllFraction");

	probe->get(dayRespirationMesophyllFraction);

	shootParameters->getChild("bundleSheathConductanceAt25C", "mol/m2/s")->get(sheathConductanceAt25C);

	shootParameters->getChild("bundleSheathConductanceActivationEnergy", "J/mol")->get(sheathConductanceActivationEnergy);

	shootParameters->getChild("bundleSheathConductanceDeactivationEnergy", "J/mol")->get(sheathConductanceDeactivationEnergy);

	shootParameters->getChild("bundleSheathConductanceEntropyTerm", "J/K/mol")->get(sheathConductanceEntropyTerm);

	shootParameters->getChild("dayRespirationMesophyllFraction")->get(dayRespirationMesophyllFraction);

}



void BundleSheathCO2Concentration::calculate(const Time &t, double &sheathC){

	double leafArea;

	pLeafArea->get(t, leafArea);

	leafArea = leafArea/10000; // convert from cm2 to m2

	if (leafArea < 1e-4){

	if (leafArea < 1e-8){

		sheathC = atmosphericCO2Concentration;

		return;

	}
@@ -220,12 +219,11 @@ void BundleSheathCO2Concentration::calculate(const Time &t, double &sheathC){ 
	leafRespiration = leafRespiration*1000000/(12.0111*60*60*24*leafArea); // convert from g/day to umol/(m2s)

	pMesophyllC->get(t, mesophyllC);

	pPEPCarboxylation->get(t, pepCarboxylation);

	double sheathRespiration = (1 - dayRespirationMesophyllFraction)*leafRespiration;

	const double sheathRespiration = (1 - dayRespirationMesophyllFraction)*leafRespiration;



		double refTemp = 298.15;

		double universalGasConstant = 8.3144598; // J/(mol*K)

		sheathConductance = sheathConductanceAt25C*exp((leafTemperature - refTemp)*sheathConductanceActivationEnergy/(refTemp*universalGasConstant*leafTemperature))*(1 + exp((sheathConductanceEntropyTerm - sheathConductanceDeactivationEnergy/refTemp)/universalGasConstant))/(1 + exp((sheathConductanceEntropyTerm - sheathConductanceDeactivationEnergy/leafTemperature)/universalGasConstant));

		if (sheathConductance == 0.) sheathConductance = 1e-10;

	const double refTemp = 298.15;

	const double universalGasConstant = 8.3144598; // J/(mol*K)

	const double sheathConductance = std::max(1e-10,sheathConductanceAt25C*exp((leafTemperature - refTemp)*sheathConductanceActivationEnergy/(refTemp*universalGasConstant*leafTemperature))*(1 + exp((sheathConductanceEntropyTerm - sheathConductanceDeactivationEnergy/refTemp)/universalGasConstant))/(1 + exp((sheathConductanceEntropyTerm - sheathConductanceDeactivationEnergy/leafTemperature)/universalGasConstant)));



	sheathC = (pepCarboxylation - photoRate + sheathRespiration)/sheathConductance + mesophyllC;

	if(std::isnan(sheathC)) msg::error("sheathC: sheathC is NaN ");

OpenSimRoot/src/modules/ShootModule/LeafGasExchange.hpp

+1 −2
Original line number Diff line number Diff line
@@ -57,9 +57,8 @@ public: 
protected:

	void getDefaultValue(const Time &t, double &var);

	void calculate(const Time &t, double &concentration);

	SimulaBase *pIrradiation, *pPhotosynthesis, *pPhotosynthesisRate, *pStomatalConductance, *pLeafArea, *pLeafRespiration, *pLeafRespirationRate, *pLeafTemperature, *pMesophyllC, *pPEPCarboxylation;

	SimulaBase *pPhotosynthesis, *pPhotosynthesisRate, *pLeafArea, *pLeafRespiration, *pLeafRespirationRate, *pLeafTemperature, *pMesophyllC, *pPEPCarboxylation;

	double atmosphericCO2Concentration, sheathConductanceAt25C, sheathConductanceActivationEnergy, sheathConductanceDeactivationEnergy, sheathConductanceEntropyTerm, dayRespirationMesophyllFraction;

	double sheathConductance, leafArea;

};



class CO2Leakage:public DerivativeBase {