Commit 1add3824 authored by jouke's avatar jouke
Browse files

fixing several errors in the Farquhar code 

-removed caches that stop iterations and convergence
-reduced scope, fixed unintialized use of variables
-fixed wrong use of get vs getRate

note that this runs but photosynthesis is still much off compared to the
lintul version.
parent 26d4d5db

OpenSimRoot/InputFiles/runMaizeFarquhar.xml

+5 −0
Original line number Diff line number Diff line
@@ -44,6 +44,11 @@ 
	</SimulaBase>	



	<SimulaIncludeFile fileName="plantParameters/Maize/Maize/simulationControlParameters.xml"/>



	<SimulaDirective path="/simulationControls/outputParameters/defaults">

		<SimulaConstant	name="timeInterval" replacesPreviousDeclaration="true"> 0.1 </SimulaConstant>

        </SimulaDirective>

    	

	<SimulaIncludeFile fileName="templates/plantTemplateFullModel.xml"/>

	<!-- SimulaIncludeFile fileName="environments/Rocksprings/Rocksprings-maize.xml"/-->

	<SimulaIncludeFile fileName="environments/WageningseBovenBuurt/WageningseBovenBuurt-maize.xml"/>

OpenSimRoot/InputFiles/templates/plantTemplate.C3C4Photosynthesis.xml

+8 −7
Original line number Diff line number Diff line
@@ -26,15 +26,16 @@ You should have received the GNU GENERAL PUBLIC LICENSE v3 with this file in lic 
<!--ShootTemplate-->

	<SimulaDirective path="/shootTemplate">

		<SimulaVariable name="photosynthesisLintul" function="photosynthesisLintulV2" unit="g" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="lightInterception" function="meanLightInterception" unit="W/m2" replacesPreviousDeclaration="true"/>

		<SimulaDerivative name="sunlitLightInterception" function="meanLightInterception" unit="W/m2"/>

		<SimulaDerivative name="shadedLightInterception" function="meanLightInterception" unit="W/m2"/>

		<SimulaDerivative name="lightInterception" function="lightInterception" unit="W/m2" replacesPreviousDeclaration="true"/>

		<!--todo, light interception needs to be updated to deal with shaded and sunlit -->

		<SimulaDerivative name="sunlitLightInterception" function="lightInterception" unit="W/m2"/>

		<SimulaDerivative name="shadedLightInterception" function="lightInterception" unit="W/m2"/>

		<SimulaVariable name="photosynthesis" function="sumOfSunlitAndShaded" unit="g" replacesPreviousDeclaration="true"/>

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

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

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

		<SimulaDerivative name="sunlitPhotosynthesisRate" function="photosynthesisRateFarquhar" unit="umol/m2/s"/>

		<SimulaDerivative name="shadedPhotosynthesisRate" function="photosynthesisRateFarquhar" unit="umol/m2/s"/>

		<SimulaVariable name="photosynthesisPerLeafArea" function="weighedAverageOfSunlitAndShaded" 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="sunlitCarbonLimitedPhotosynthesisRate" function="carbonLimitedPhotosynthesisRate" unit="umol/m2/s"/>

		<SimulaDerivative name="shadedCarbonLimitedPhotosynthesisRate" function="carbonLimitedPhotosynthesisRate" unit="umol/m2/s"/>
@@ -75,7 +76,7 @@ You should have received the GNU GENERAL PUBLIC LICENSE v3 with this file in lic 
		

		<SimulaDerivative name="leafIrradiation" function="leafIrradiation" unit="W/m2"/>

		<SimulaDerivative name="sunlitLeafIrradiation" function="sunlitLeafIrradiation" unit="W/m2"/>

		<SimulaDerivative name="shadedLeafIrradiation" function="leafIrradiation" unit="W/m2"/>

		<SimulaDerivative name="shadedLeafIrradiation" function="shadedLeafIrradiation" unit="W/m2"/>

		

		

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

OpenSimRoot/src/modules/Atmosphere/DiffusionResistance.cpp

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


/********************stomatal resistance****************************/

StomatalResistance::StomatalResistance(SimulaDynamic* const pSV) :

		DerivativeBase(pSV), cachedTime(-10), temperature(25), pressure(101325), invert(false) {

		DerivativeBase(pSV), temperature(25), pressure(101325), invert(false) {



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

	pSD->checkUnit("s/m");

OpenSimRoot/src/modules/Atmosphere/DiffusionResistance.hpp

+1 −1
Original line number Diff line number Diff line
@@ -34,7 +34,7 @@ protected: 
	void calculate(const Time &t, double& r_s);



	SimulaBase *stomatalConductance_, *airTemperature_, *pressure_, *dailyBulkStomatalResistance_, *leafAreaIndex_;

	double cachedTime, temperature, pressure;

	double temperature, pressure;

	bool invert;

};

OpenSimRoot/src/modules/Atmosphere/EvapoEquations.cpp

+2 −20
Original line number Diff line number Diff line
@@ -97,7 +97,7 @@ DerivativeBase * newInstantiationPenmanMonteith( 


/*******************Penman*************************/

 Penman::Penman(SimulaDynamic* const pSV) :

 		ETbaseclass(pSV), cachedTime(-10) {

 		ETbaseclass(pSV) {

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

	airDensity_					= pSD->getSibling("airDensity","kg/m3");

	aerodynamicResistance_		= pSD->getSibling("aerodynamicResistance","s/m");
@@ -116,11 +116,6 @@ DerivativeBase * newInstantiationPenmanMonteith( 
 /***************************************************************************/

	this->checkMode(t);

	//--------------------------------------------------------------------

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

		ETCROP = cachedCropET;

		ETSOIL = cachedSoilET;

		return;

	}

	double temperature;

	dailyTemperature_->get(t, temperature);

	/*********************************************************************/
@@ -169,9 +164,6 @@ DerivativeBase * newInstantiationPenmanMonteith( 
//		interception_->get(t, E_i);

//		ETCROP += E_i;

//	}

	cachedTime = t;

	cachedCropET = ETCROP;

	cachedSoilET = ETSOIL;



}
@@ -179,7 +171,7 @@ DerivativeBase * newInstantiationPenmanMonteith( 
 UseET0fromTable::UseET0fromTable(SimulaDynamic* const pSV) :

		DerivativeBase(pSV), ET0crop_(nullptr), ET0soil_(nullptr),

		interception_(nullptr), leafAreaIndex_(nullptr), splitparam_(nullptr),

		cachedTime(-10), cachedSoilET(0), cachedCropET(0),unitconv_(1.0),mode(0) {

		unitconv_(1.0),mode(0) {

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

	SimulaBase* p = ORIGIN->getChild("environment")->getChild("atmosphere");

 	ET0crop_		= p->existingChild("ET0crop");
@@ -264,11 +256,6 @@ DerivativeBase * newInstantiationPenmanMonteith( 
  /***************************************************************************/

 	//this->checkMode(t);

 	//--------------------------------------------------------------------

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

 		ETCROP = cachedCropET;

 		ETSOIL = cachedSoilET;

 		return;

 	}



 	ET0crop_->get(t,ETCROP);

 	if(ET0soil_){
@@ -279,11 +266,6 @@ DerivativeBase * newInstantiationPenmanMonteith( 
 	ETCROP*=unitconv_;

 	ETSOIL*=unitconv_;





 	cachedTime = t;

 	cachedCropET = ETCROP;

 	cachedSoilET = ETSOIL;



 }