Commit 38da4143 authored by jouke's avatar jouke
Browse files

Small fix for secondary growth in the carbon scaling factor

parent 5a750957

OpenSimRoot/src/modules/CarbonModule/CarbonAllocation.cpp

+16 −22
Original line number Diff line number Diff line
@@ -633,10 +633,11 @@ void RootGrowthScalingFactor::calculate(const Time &t,double &scale){ 
	}



	//carbon needed for secondary growth

	double s(0),sp(0);//actual and potential carbon for secondary growth

	double s(0);//actual and potential carbon for secondary growth

	if(c2SecondaryGrowth && a>0){

		c2SecondaryGrowth->getRate(t,sp);

		s=sp;//s is potential

		c2SecondaryGrowth->getRate(t,s);

	}

	const double sp(s);

	//make sure secondary growth is not larger than threshold

	double thrh(1.0);//don't change this a backward compatibility default

	if(threshold)threshold->get(t-plantingTime,thrh);
@@ -645,7 +646,6 @@ void RootGrowthScalingFactor::calculate(const Time &t,double &scale){ 
		s=thrh;

		msg::warning("RootGrowthScaling Factor: parameter maxCarbonAllocation2SecondaryGrowth limiting secondary growth.");

	}

	}



	//secondary axis.

	//preferences secondary axis: see borch et al., 1999 (bean) and french paper on maize
@@ -653,30 +653,24 @@ void RootGrowthScalingFactor::calculate(const Time &t,double &scale){ 
	if(mode!=2){

		//Subtract allocation for secondary growth

		a-=s;

		p-=s;

		p-=sp;

		//subtract carbon needed for major axis

		if(potentialRootGrowthSimulatorMajorAxis && p>0 ){

			potentialRootGrowthSimulatorMajorAxis->getRate(t,mp);

			if(mp>p) mp=p;//should not happen but just in case there are some round of errors.

			//see if available carbon does not match sink

			if(a<p*0.9995){

			//shortage of carbon

				//sink fine roots

				double fp((p-mp)/p);

			if(a<p){

				//carbon shortage, decrease ma

				double thrsh(1-(fp*0.5));

				double thrsh(0.5+0.5*mp/p);

				if(mp<thrsh*a){

					//give preference to major axis

					ma=mp;

				}else{

					ma=thrsh*a;

				}

			}else if(a>p*1.0005){

			}else{

				//carbon extra increase ma

			//jp should not happen any longer given that a>p is not really allowed above.	

				ma=mp*a/p;

			}else{

				ma=mp;

			}

		}

	}