Commit 252aaf14 authored by jouke's avatar jouke
Browse files

Update respiration code 

Merged stem and leaf respiration into one class
Made root respiration part of total base
Added class for growth respiration.
Removed time caching introduced by Ernst
Support for both C and K temperatures
parent f69c0b97

OpenSimRoot/src/modules/CarbonModule/Respiration.cpp

+114 −95
Original line number Diff line number Diff line
@@ -35,7 +35,7 @@ NOTE: The GPL.v3 license requires that all derivative work is distributed under 


//Respiration rate at any point in the rootsystem in g CO2/cm root/day

bool RootSegmentRespirationRate::issueMessage(true);

RootSegmentRespirationRate::RootSegmentRespirationRate(SimulaDynamic* pSD):DerivativeBase(pSD),sizeSimulator(nullptr), relativeRespirationSimulator(nullptr),

RootSegmentRespirationRate::RootSegmentRespirationRate(SimulaDynamic* pSD):TotalBase(pSD),sizeSimulator(nullptr), relativeRespirationSimulator(nullptr),

		factor(nullptr), aerenchymaSimulator(nullptr), aerenchymaCorrectionSimulator(nullptr),mode(0)

{

	//check if unit given in input file agrees with this function
@@ -178,7 +178,7 @@ void RootSegmentRespirationRate::calculate(const Time &t, double &rate) { 
	//current root segment respiration rate in g CO2/day

	rate *= r;



	//multiplication factor

	//multiplication factor associated with stress impact

	if (factor) {

		factor->get(t, r);

		rate *= r;
@@ -211,9 +211,8 @@ DerivativeBase * newInstantiationRootSegmentRespirationRate( 


//shoot respiration

LeafRespirationRate::LeafRespirationRate(SimulaDynamic *pSD) :

		DerivativeBase(pSD),

sizeSimulator(NULL), leafSenescenceSimulator(NULL), relativeRespirationSimulator(NULL), factor(NULL)

{

		DerivativeBase(pSD), sizeSimulator(nullptr), leafSenescenceSimulator(nullptr), relativeRespirationSimulator(nullptr),

		factor(nullptr), growthRespiration(nullptr), rgrowthRespiration(nullptr), refTime(pSD->getStartTime()) {

	//check if unit given in input file agrees with this function

	pSD->checkUnit("g");

	//plant type
@@ -224,6 +223,7 @@ sizeSimulator(NULL), leafSenescenceSimulator(NULL), relativeRespirationSimulator 
	//get the root type parameters

	SimulaBase *parameters(GETSHOOTPARAMETERS(plantType));

	//find simulators in database

	if(pSD->getName().substr(0,4)=="leaf"){

	relativeRespirationSimulator = parameters->existingChild(

			"relativeRespirationRateLeafs", t);

	if (relativeRespirationSimulator) {
@@ -243,12 +243,45 @@ sizeSimulator(NULL), leafSenescenceSimulator(NULL), relativeRespirationSimulator 
				"LeafRespirationRate: relativeRespirationRateLeafs parameter missing, setting respiration to 0.");

	}

	growthRespiration = pSD->existingSibling("leafGrowthRespirationRate","g/day");

	if(!growthRespiration){

		rgrowthRespiration = parameters->existingChild("leafRelativeGrowthRespirationRate", t);

	}

	if(growthRespiration) msg::warning("leafRespiration: adding leafGrowthRespiration, assuming relativeRespirationRateLeafs is maintenance respiration only");



	//check for nutrients effects on respiration

	factor = pSD;

	PLANTTOP(factor);

	factor = factor->existingChild(

			"stressFactor:impactOn:leafRespiration");

	}else{

		relativeRespirationSimulator = parameters->existingChild(

				"relativeRespirationRateStems", t);

		if (relativeRespirationSimulator) {

			if (relativeRespirationSimulator->getUnit() == "g/g/day") {

				sizeSimulator = pSD->getSibling(

						"stemDryWeight");

				leafSenescenceSimulator = pSD->existingSibling("senescedStemDryWeight"); //todo this is not implemented elsewhere

			} else {

				msg::error(

						"StemRespirationRate: relativeRespirationRateStems should be in g/g/day");

			}

		} else {

			msg::warning(

					"stemRespirationRate: relativeRespirationRateStems missing, setting respiration to 0.");

		}

		growthRespiration = pSD->existingSibling("stemGrowthRespirationRate","g/day");

		if(!growthRespiration){

			rgrowthRespiration = parameters->existingChild("stemRelativeGrowthRespirationRate", t);

		}

		if(growthRespiration) msg::warning("stemRespiration: adding stemGrowthRespiration, assuming relativeRespirationRateStems is maintenance respiration only");



		//check for nutrients effects on respiration

		factor = pSD;

		PLANTTOP(factor);

		factor = factor->existingChild(

				"stressFactor:impactOn:stemRespiration");



	}

}

void LeafRespirationRate::calculate(const Time &t, double &rate) {

	//return 0 if repiration parameters are not set
@@ -267,7 +300,7 @@ void LeafRespirationRate::calculate(const Time &t, double &rate) { 


	//current respiration rate in g CO2 /g/day or g CO2 /cm2/day

	double res;

	relativeRespirationSimulator->get(t, res);

	relativeRespirationSimulator->get(t-refTime, res);



	//current shoot respiration rate in g CO2/day

	rate *= res;
@@ -282,73 +315,29 @@ void LeafRespirationRate::calculate(const Time &t, double &rate) { 
		double grr;

		growthRespiration->get(t,grr);

		rate+=grr;

	}else if(rgrowthRespiration){

		double gr;//growth rate

		sizeSimulator->getRate(t, gr);

		double grr;

		rgrowthRespiration->get(t-refTime,grr);

		rate += grr*gr;

	}else{

		//do nothing

	}

}



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

	return "leafRespirationRate";

	return "leafOrStemRespirationRate";

}

DerivativeBase * newInstantiationLeafRespirationRate(SimulaDynamic* const pSD) {

	return new LeafRespirationRate(pSD);

}



//shoot respiration

StemRespirationRate::StemRespirationRate(SimulaDynamic* pSD) :

		DerivativeBase(pSD) {

	//check if unit given in input file agrees with this function

	pSD->checkUnit("g");

	//plant type

	std::string plantType;

	PLANTTYPE(plantType, pSD);

	//get the root type parameters

	SimulaBase *parameters(GETSHOOTPARAMETERS(plantType));

	//find simulators in database

	relativeRespirationSimulator = parameters->existingChild(

			"relativeRespirationRateStems", "g/g/day");

	if (relativeRespirationSimulator) {

		sizeSimulator = pSD->getSibling("stemDryWeight");

	} else {

		msg::warning(

				"StemRespirationRate: relativeRespirationRateStems parameter missing, setting respiration to 0.");

	}

	//check for nutrients effects on respiration

	factor = pSD;

	PLANTTOP(factor);

	factor = factor->existingChild(

			"stressFactor:impactOn:stemRespiration");

}

void StemRespirationRate::calculate(const Time &t, double &rate) {

	//return 0 if repiration parameters are not set

	if (!relativeRespirationSimulator) {

		rate = 0;

		return;

	}



	//shoot size in cm2 leaf area or g dryweight

	sizeSimulator->get(t, rate);



	//current respiration rate in g CO2 /g/day or g CO2 /cm2/day

	double res;

	relativeRespirationSimulator->get(t, res);



	//current shoot respiration rate in g CO2/day

	rate *= res;



	//multiplication factor

	if (factor) {

		factor->get(t, res);

		rate *= res;

	}

}

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

	return "stemRespirationRate";

}



DerivativeBase * newInstantiationStemRespirationRate(SimulaDynamic* const pSD) {

	return new StemRespirationRate(pSD);

	return new LeafRespirationRate(pSD);//these classes have been merged

}



LeafRespirationRateFarquhar::LeafRespirationRateFarquhar(SimulaDynamic* pSD):DerivativeBase(pSD), pLeafDryWeight(nullptr), pStemDryWeight(nullptr), cachedTime(-9999), cachedLeafTemperature(-9999){

LeafRespirationRateFarquhar::LeafRespirationRateFarquhar(SimulaDynamic* pSD):DerivativeBase(pSD), pLeafDryWeight(nullptr), pStemDryWeight(nullptr){

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

	std::string plantType;

	PLANTTYPE(plantType,pSD);
@@ -372,7 +361,6 @@ LeafRespirationRateFarquhar::LeafRespirationRateFarquhar(SimulaDynamic* pSD):Der 
}



void LeafRespirationRateFarquhar::calculate(const Time &t, double &rate) {

	if (std::abs(t - cachedTime) > TIMEERROR){

	pLeafArea->get(t, leafArea);

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

	if(pLeafSenescence){
@@ -396,21 +384,17 @@ void LeafRespirationRateFarquhar::calculate(const Time &t, double &rate) { 
		}

		leafArea = leafArea*stemDryWeight/leafDryWeight;

	}

		pReferenceDayRespiration->get(t, referenceDayRespiration);

		cachedTime = t;

	}

	pReferenceDayRespiration->get(t-pSD->getStartTime(), referenceDayRespiration);



	double leafTemperature; // degrees

	pLeafTemperature->get(t, leafTemperature);

	leafTemperature = std::max(leafTemperature, 0.);

	leafTemperature = std::min(leafTemperature, 100.); // Leaf can not be frozen or boiling, this prevents errors that can occur when default output timestep is larger than 0.1

	leafTemperature = leafTemperature + 273.15; // Need leaf temperature in Kelvin

	if (std::abs(leafTemperature - cachedLeafTemperature) > 0.01){

		cachedLeafTemperature = leafTemperature;

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

		double refTemp = 298.15;

	if(leafTemperature<200) leafTemperature += 273.15; // Need leaf temperature in Kelvin

	leafTemperature = std::max(leafTemperature, 273.15);

	leafTemperature = std::min(leafTemperature, 373.15); // Leaf can not be frozen or boiling, this prevents errors that can occur when default output timestep is larger than 0.1

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

	const double refTemp = 298.15;

	temperatureScalingFactor = exp((leafTemperature-refTemp)*activationEnergyDayRespiration/(refTemp*universalGasConstant*leafTemperature)); // umol/m2/s

	}

	double dayRespiration = referenceDayRespiration*temperatureScalingFactor;

	const double dayRespiration = referenceDayRespiration*temperatureScalingFactor;

	rate = dayRespiration*leafArea*60*60*24*12.0111/1000000; // Convert from umol/m2/s to g/d

}
@@ -421,12 +405,45 @@ DerivativeBase * newInstantiationLeafRespirationRateFarquhar(SimulaDynamic* cons 
	return new LeafRespirationRateFarquhar(pSD);

}



GrowthRespirationRate::GrowthRespirationRate(SimulaDynamic* pSD):DerivativeBase(pSD), pGrowthRate(nullptr), pRelativeRate(nullptr){

	std::string organ=pSD->getName().substr(0,4);

	pGrowthRate=pSD->getSibling(organ+"DryWeight","g");

	SimulaBase *param;

	std::string plantType;

	PLANTTYPE(plantType,pSD);

	if(organ=="root"){

		std::string rootType;

		pSD->getParent(3)->getChild("rootType")->get(rootType);

		param=GETROOTPARAMETERS(plantType, rootType);

	}else{

		param = GETSHOOTPARAMETERS(plantType);

	}

	pRelativeRate=param->getChild(organ+"RelativeGrowthRespirationRate");

}



void GrowthRespirationRate::calculate(const Time &t,double &var){

	pGrowthRate->getRate(t,var);

	double r;

	pRelativeRate->get(t-pSD->getStartTime(), r);

	var*=r;

};



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

	return "growthRespirationRate";

}

DerivativeBase * newInstantiationGrowthRespirationRate(SimulaDynamic* const pSD) {

	return new GrowthRespirationRate(pSD);

}





//Register the module

class AutoRegisterRespirationInstantiationFunctions {

public:

	AutoRegisterRespirationInstantiationFunctions() {

		BaseClassesMap::getDerivativeBaseClasses()["rootSegmentRespirationRate"] =

				newInstantiationRootSegmentRespirationRate;

		BaseClassesMap::getDerivativeBaseClasses()["leafOrStemRespirationRate"] =

				newInstantiationLeafRespirationRate;

		BaseClassesMap::getDerivativeBaseClasses()["leafRespirationRate"] =

				newInstantiationLeafRespirationRate;

		BaseClassesMap::getDerivativeBaseClasses()["stemRespirationRate"] =
@@ -435,6 +452,8 @@ public: 
				newInstantiationLeafRespirationRateFarquhar;

		BaseClassesMap::getDerivativeBaseClasses()["stemRespirationRateFarquhar"] =

				newInstantiationLeafRespirationRateFarquhar;

		BaseClassesMap::getDerivativeBaseClasses()["growthRespirationRate"] =

				newInstantiationGrowthRespirationRate;

	}

};

OpenSimRoot/src/modules/CarbonModule/Respiration.hpp

+14 −3
Original line number Diff line number Diff line
@@ -24,6 +24,7 @@ NOTE: The GPL.v3 license requires that all derivative work is distributed under 
#ifndef RESPIRATION_HPP_

#define RESPIRATION_HPP_

#include "../../engine/BaseClasses.hpp"

#include "../GenericModules/Totals.hpp"



/**

 * Calculation of the respiration of a rootsegment.
@@ -51,7 +52,7 @@ NOTE: The GPL.v3 license requires that all derivative work is distributed under 
 *

 *

 */

class RootSegmentRespirationRate:public DerivativeBase{

class RootSegmentRespirationRate:public TotalBase{

public:

	RootSegmentRespirationRate(SimulaDynamic* pSD);

	std::string getName()const;
@@ -68,7 +69,8 @@ public: 
	std::string getName()const;

protected:

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

	SimulaBase *sizeSimulator, *leafSenescenceSimulator, *relativeRespirationSimulator, *factor, *growthRespiration;

	SimulaBase *sizeSimulator, *leafSenescenceSimulator, *relativeRespirationSimulator, *factor, *growthRespiration, *rgrowthRespiration;

	Time refTime;

};



class StemRespirationRate:public DerivativeBase{
@@ -88,7 +90,16 @@ protected: 
	void calculate(const Time &t,double &var);

	double activationEnergyDayRespiration;

	SimulaBase *pLeafArea, *pLeafSenescence, *pLeafTemperature, *pLeafDryWeight, *pStemDryWeight, *pReferenceDayRespiration;

	double cachedTime, leafArea, cachedLeafTemperature, referenceDayRespiration, temperatureScalingFactor;

	double  leafArea, referenceDayRespiration, temperatureScalingFactor;

};



class GrowthRespirationRate:public DerivativeBase{

public:

	GrowthRespirationRate(SimulaDynamic* pSD);

	std::string getName()const;

protected:

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

	SimulaBase *pGrowthRate, *pRelativeRate;

};



#endif /*RESPIRATION_HPP_*/