Commit ca14d972 authored by jouke's avatar jouke
Browse files

code simplification, and correct treatment of units

parent a2c9a0b9
Loading
Loading
Loading
Loading
+23 −23
Original line number Diff line number Diff line
@@ -22,26 +22,33 @@ 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) {
		DerivativeBase(pSV), cachedTime(-10), temperature(25), pressure(101325), invert(false) {

	//check if unit given in input file agrees with this function
	pSD->checkUnit("s/m");


	std::string name = pSD->getName().substr(0, 6); // name = sunlit or shaded
	bool splitBySunStatus(false);
	if (name == "sunlit" || name == "shaded"){
		splitBySunStatus = true;
		name.at(0) = std::toupper(name.at(0));
		stomatalConductance_ = pSD->existingSibling("mean" + name + "StomatalConductance", "mol/m2/s");
	}else{
		stomatalConductance_ = pSD->existingSibling("meanStomatalConductance","mol/m2/s");
	}
	//check if unit given in input file agrees with this function
	pSD->checkUnit("s/m");
	if (splitBySunStatus) stomatalConductance_ = pSD->existingSibling("mean" + name + "StomatalConductance", "mol/m2/s");
	else stomatalConductance_ = pSD->existingSibling("meanStomatalConductance","mol/m2/s");
	if (stomatalConductance_){
		airTemperature_ = pSD->existingPath("/environment/atmosphere/averageDailyTemperature", "degreesC");
		pressure_ = pSD->existingPath("/atmosphere/airPressure", "Pa");
		if (!airTemperature_) msg::warning("StomatalResistance: Air temperature not found, defaulting to 25 degrees C");
		if (!pressure_) msg::warning("StomatalResistance: Air pressure not found, defaulting to 101325 Pa");
	}
	}else{
		dailyBulkStomatalResistance_ = pSD->existingSibling("dailyBulkStomatalResistance");
	if(!dailyBulkStomatalResistance_) msg::warning("StomatalResistance: dailyBulkStomatalResistance as sibling of "+pSD->getPath()+"not found, using default value of 100 m/s");
	leafAreaIndex_ 				 = pSD->getSibling("meanLeafAreaIndex"); //, "m/s"
		if(!dailyBulkStomatalResistance_) {
			msg::warning("StomatalResistance: dailyBulkStomatalResistance as sibling of "+pSD->getPath()+"not found, using default value of 100 s/m. Consider defining meanStomatalConductance in mol/m2/s.");
		}else{
			dailyBulkStomatalResistance_->checkUnit("s/m");
		}
	}
	leafAreaIndex_ = pSD->getSibling("meanLeafAreaIndex");
}

std::string StomatalResistance::getName() const {
@@ -74,24 +81,17 @@ void StomatalResistance::calculate(const Time &t, double& r_s) {
//		leaf temperature, and vapor pressure gradient ( Jarvis, 1976, Stewart, 1989 and Price and Black, 1989)
//		and rl increases with environmental stresses such as soil moisture deficit
//		( Stewart, 1988, Stewart and Verma, 1992 and Hatfield and Allen, 1996).
// TODO coupling
		double r_l(100.); // default value, (m/s) is the daily bulk stomatal resistance of the well-illuminated (single) leaf
		double r_l(100.); // default value, (s/m)
		// The bulk stomatal resistance, rl, is the average resistance of an individual leaf,
		// and has a value of about 100 s/m under well-watered conditions.
		// and has a value of about 1/100 s/m under well-watered conditions.
		// TODO could be fitted to a crop average?
		if (dailyBulkStomatalResistance_){
			dailyBulkStomatalResistance_->get(t, r_l);
		}
		double LAI = 3;
		leafAreaIndex_->get(t,LAI);
		double lai_active;
		if(LAI > 0.8/0.3)
			lai_active = 0.5 * LAI;
		else if(LAI >= 1){
			lai_active = LAI / (0.3 * LAI+ 1.2);
		}else{
			lai_active=1./1.5;
		}
		// greater resistance (more conductance) when LAI is larger
		const double frac=std::min(2.0,std::max(1.5,0.3*LAI+1.2));
		// the active (sunlit) leaf area index, which takes into consideration
		/// the fact that generally only the upper half of dense clipped grass
		/// is actively contributing to the surface heat and vapour transfer.
@@ -103,7 +103,7 @@ void StomatalResistance::calculate(const Time &t, double& r_s) {
		/// If the crop is not transpiring at a potential rate, the resistance depends also on the water status of the vegetation.
		/// An acceptable approximation to a much more complex relation of the surface resistance of dense full cover vegetation is:
		//
		r_s = r_l / lai_active;
		r_s = r_l * frac;
		//
		// This procedure also presents difficulties. Firstly, evaporation losses from soil must be
		// negligible. If not, r_s will not be representative of total evaporation losses.
+1 −0
Original line number Diff line number Diff line
@@ -35,6 +35,7 @@ protected:

	SimulaBase *stomatalConductance_, *airTemperature_, *pressure_, *dailyBulkStomatalResistance_, *leafAreaIndex_;
	double cachedTime, temperature, pressure;
	bool invert;
};