Commit c926e151 authored by jouke's avatar jouke
Browse files

Fixed splitting of ET0 into two proportions that are garanteed 100% ET0

parent 2398b52a

OpenSimRoot/src/modules/Atmosphere/ETbase.cpp

+27 −21
Original line number Diff line number Diff line
@@ -56,10 +56,11 @@ NOTE: The GPL.v3 license requires that all derivative work is distributed under 
	 *	EVAPORATION OF A FULLY DEVELOPED PLANT STOCK WITHOUT WATER STRESS

	 */

	/************************************************************************/

ETbaseclass::ETbaseclass(SimulaDynamic* const pSV):DerivativeBase(pSV), splitBySunStatus(false){

ETbaseclass::ETbaseclass(SimulaDynamic* const pSV):DerivativeBase(pSV), splitBySunStatus(0){

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

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

	if (name == "sunlit" || name == "shaded") splitBySunStatus = true;

	if (name == "sunlit") splitBySunStatus = 1;

	if (name == "shaded") splitBySunStatus = 2;

	if(pSD->getType()=="SimulaVariable"){

		pSD->checkUnit("cm");

	}else{
@@ -88,15 +89,13 @@ ETbaseclass::ETbaseclass(SimulaDynamic* const pSV):DerivativeBase(pSV), splitByS 


	}



	extinctionCoef_				= pSD->getPath("/plants/meanExtinctionCoefficient");

	leafAreaIndex_ 				= pSD->getPath("/plants/meanLeafAreaIndex"); //, "m/s"

	if (splitBySunStatus){

		stomatalResistance_ = pSD->existingPath("/plants/" + name + "StomatalResistance", "s/m");

		name.at(0) = std::toupper(name.at(0));

		leafAreaIndex_ = pSD->existingPath("/plants/mean" + name + "LeafAreaIndex");

		netRadiation_ = pSD->getPath("/plants/mean" + name + "LightInterception");

		splitLAI_ = pSD->existingPath("/plants/meanSunLitLeafAreaIndex");//this must be the sunlit, we compute 1-sunlit later

	}

	else{

		leafAreaIndex_ 				= pSD->getPath("/plants/meanLeafAreaIndex"); //, "m/s"

		splitparam_					= pSD->getPath("/plants/meanExtinctionCoefficient");

		stomatalResistance_ 		= pSD->existingPath("/plants/stomatalResistance", "s/m");

		//this radiation corrected for reflection, and since the albedo of the crop is different from that of soil, we have two numbers

		netRadiation_				= pSD->existingSibling("netRadiation");//looks if this is simulated with for example the Radiation class in Radiation.cpp
@@ -175,19 +174,27 @@ ETbaseclass::ETbaseclass(SimulaDynamic* const pSV):DerivativeBase(pSV), splitByS 
//ETbaseclass::MyMap ETbaseclass::mymap;



void ETbaseclass::calculate(const Time &t, double& ET){

	double LAI;

	leafAreaIndex_->get(t,LAI);



	double ETCROP, ETSOIL;

	calculateET(t, ETCROP, ETSOIL );



	if (!splitBySunStatus){

		double KDF;

		splitparam_->get(t,KDF);

		ETCROP = ETCROP*(1.-exp(-KDF *LAI));

	//	ETSOIL = ETSOIL*exp(-KDF *LAI);// exp(..) equals 1-SC, Soil without crop // <- This should not be scaled like this because the soil radiation should be scaled according to leaf area and the component not dependent on radiation should not be scaled

	} else{

		ETCROP = ETCROP*LAI;



	//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);

	ETSOIL *= (split);

	ETCROP *= (1-split);



	if (splitBySunStatus){

		double LAIfrac;

		splitLAI_->get(t,LAIfrac);

		double frac=exp(-KDF *LAIfrac);

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

		ETCROP *= LAIfrac/LAI;

	}



	switch (mode) {

		case 1:

			ET = ETSOIL;
@@ -197,6 +204,7 @@ void ETbaseclass::calculate(const Time &t, double& ET){ 
			break;

		default:

			ET = ETCROP + ETSOIL;

			//jp todo: is this correct? Where does the energy come from?

			if (interception_) {

					double E_i = 0.;

					interception_->get(t, E_i);
@@ -252,16 +260,14 @@ void ETbaseclass::calculateRadiation(const Time &t,double & H_soil, double & H_ 
	} else {

		//msg::warning("ETbaseclass::calculateRadiation: Net Radiation Rn is computed. ");

		//msg::warning("ETbaseclass::calculateRadiation: Soil Heat Flux G is assumed to be zero! ");

		double G = 0.; // for day and ten-day periods

		const double G = 0.; // for day and ten-day periods

		double Rn_soil;

		netRadiationSoil_->get(t, Rn_soil);

		Rn_soil = Rn_soil*convS;

		H_soil = Rn_soil - G;

		H_soil = Rn_soil*convS - G;



		double Rn_crop;

		netRadiation_->get(t, Rn_crop);

		Rn_crop = Rn_crop*conv;

		H_crop = Rn_crop - G;

		H_crop = Rn_crop*conv - G;

	}

}

OpenSimRoot/src/modules/Atmosphere/ETbase.hpp

+2 −2
Original line number Diff line number Diff line
@@ -35,9 +35,9 @@ protected: 
	*slopeVaporPressure_, *airPressure_,*specificHeatCapacityOfAir_,

			*sensibleHeat_, *netRadiation_,*netRadiationSoil_,

			*actualVaporPressure_, *saturatedVaporPressure_,

			*interception_, *rhp_,*leafAreaIndex_, *splitparam_, *stomatalResistance_;

			*interception_, *rhp_,*leafAreaIndex_, *splitLAI_, *extinctionCoef_, *stomatalResistance_;

	double conv, convS,convSH;

	bool splitBySunStatus;

	int splitBySunStatus;





	// associate a certain string with a number for formula selection,