Commit 84009e77 authored by jouke's avatar jouke
Browse files

Verified and improved the airodynamic resistance code now with a minimum

crop height of 10 cm, to account for soil roughness.
parent e8552633
Loading
Loading
Loading
Loading
+20 −16
Original line number Diff line number Diff line
@@ -133,7 +133,9 @@ AerodynamicResistance::AerodynamicResistance(SimulaDynamic* const pSV) :
	pSD->checkUnit("s/m");
	cropHeight_ = pSD->existingPath("/plants/canopyHeight", "cm");
	if (!cropHeight_) msg::warning("AerodynamicResistance: /plants/canopyHeight not found. Assuming 0.5 m");
	cropHeight_->checkUnit("cm");
	windSpeed_ 		 = pSD->existingPath("/environment/atmosphere/windSpeed"); //, "m/s"
	windSpeed_->checkUnit("m/s");
}

std::string AerodynamicResistance::getName() const {
@@ -159,37 +161,39 @@ void AerodynamicResistance::calculate(const Time &t, double& r_a) {
	/// exchange of vapor at the surface when air is calm. This effect occurs when the wind speed is small and
	/// buoyancy of warm air induces air exchange at the surface. Limiting u2 >= 0.5 m/s in the ETo equation improves
	/// the estimation accuracy under the conditions of very low wind speed.
	double windSpeed; // Windspeed at measurement height, which we assume to be 2 meters.
	double windSpeed(2); // Windspeed at measurement height, which we assume to be 2 meters.
	if (windSpeed_) {
		/// --- Wind speed ---  /// --- Weibull distributed ---
		windSpeed_->get(t, windSpeed); ///< U2 being the wind speed at 2-m height in m/s

		// The aerodynamic resistance equation for open air estimations reduces to r_a = 208./U2;
		// An alternative, semi-empirical, equation of ra (Thom and Oliver, 1977)
		// is used for greenhouse conditions (Stanghellini equation), characterized by low wind speed (<1 m s-1)
		// There is always a little bit of wind or convection
		if (windSpeed < 0.5) windSpeed = 0.5;
	} else {
		windSpeed = 2.;
		msg::warning("AerodynamicResistance: No wind speed is set by user. Assuming a wind speed of 2 m/s.");
	}
	// There is always a little bit of wind or convection
	if (windSpeed < 0.5) windSpeed = 0.5;
	// Allen, Richard G., et al. "Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56." Fao, Rome 300.9 (1998): D05109.
	// See also http://www.fao.org/3/X0490E/x0490e06.htm
	double vonKarmanConstant = 0.41;
	const double vonKarmanConstant = 0.41;
	double measurementHeight = 2;
	double zeroPlaneDisplacementHeight = 50;
	double cropHeight = 0.50; //m
	if (cropHeight_) {
		cropHeight_->get(t, zeroPlaneDisplacementHeight);
		if(zeroPlaneDisplacementHeight<1)zeroPlaneDisplacementHeight=1;
		cropHeight_->get(t, cropHeight);
		cropHeight/=100; //cm to m
		if(cropHeight<0.1)cropHeight=0.1; //asume crop height is at least 10 cm, given that the soil will be uneven.
	}
	zeroPlaneDisplacementHeight = zeroPlaneDisplacementHeight/100.; // Convert from cm to m
	double roughnessLengthMomentumTransfer = 0.123*zeroPlaneDisplacementHeight;
	double roughnessLengthHeatAndVapour = 0.1*0.123*zeroPlaneDisplacementHeight;
	zeroPlaneDisplacementHeight *= 2./3.;
	measurementHeight=std::max(measurementHeight,cropHeight+1);// make sure the measurement height is at least 1 m above the crop.

	const double zeroPlaneDisplacementHeight = 2./3. * cropHeight; // Convert from cm to m
	const double roughnessLengthMomentumTransfer = 0.123*cropHeight;
	const double roughnessLengthHeatAndVapour = 0.1*0.123*cropHeight;
	r_a = log((measurementHeight - zeroPlaneDisplacementHeight)/roughnessLengthMomentumTransfer)*log((measurementHeight - zeroPlaneDisplacementHeight)/roughnessLengthHeatAndVapour)/(vonKarmanConstant*vonKarmanConstant*windSpeed);
	if (std::isnan(r_a)) msg::error("AerodynamicResistance: NaN value");
	if (std::isinf(r_a)) msg::error("AerodynamicResistance: inf value");
	if (r_a <= 0) msg::error("AerodynamicResistance: Aerodynamic resistance smaller than or equal to 0. Check if crop height is less than 2, if not, adjust measurement height in the code.");


}