#include "ResistanceCal.h"
//include "HeatMetrics_Calc.h" to access HeatMetrics_Calc.cpp functionality 
#include "HeatMetrics_Calc.h"

//Note: Li and Bou-Zeid (2014) provide alternative resistance formulation through conductance term (Eq 4a), with z0 defined by land use in WRF and temporally variable zoT defined by Eq 5b and 5c (MYJ) and (CZ09)
//Li, D., & Bou-Zeid, E. (2014). Quality and sensitivity of high-resolution numerical simulation of urban heat islands. Environmental Research Letters, 9(5), 055001. doi:10.1088/1748-9326/9/5/055001

//References:
//Allen, R., Pereira, L., Raes, D., Smith, M. (1998). Crop evapotranspiration - Guidelines for computing crop water requirements - FAO Irrigation and drainage paper 56. 
//Choudhury, B. J., & Monteith, J. L. (1988). A four-layer model for the heat budget of homogeneous land surfaces. Quarterly Journal of the Royal Meteorological Society, 114(480), 373-398. doi:https://doi.org/10.1002/qj.49711448006
//Hicks, B. B., Baldocchi, D. D., Meyers, T. P., Hosker, R. P., & Matt, D. R. (1987). A preliminary multiple resistance routine for deriving dry deposition velocities from measured quantities. Water, Air, and Soil Pollution, 36(3), 311-330. doi:10.1007/BF00229675
//Jensen, M. E., & Allen, R. G. (2015). Evaporation, Evapotranspiration, and Irrigation Water Requirements: Manual of Practice 70, Second Edition (2nd ed.). Reston, VA: American Society of Civil Engineers.
//Kjelgren, R., & Montague, T. (1998). Urban tree transpiration over turf and asphalt surfaces. Atmospheric Environment, 32(1), 35-41. doi:https://doi.org/10.1016/S1352-2310(97)00177-5 
//Lee, S. H., & Park, S. U. (2008). A vegetated urban canopy model for meteorological and environmental modelling. Boundary-Layer Meteorology, 126(1), 73-102. doi:10.1007/s10546-007-9221-6
//Lindroth, A. (1993).Aerodynamicand Canopy Resistance of Short - Rotation Forest in Relation to Leaf - Area Index and Climate.Boundary - Layer Meteorology, 66(3), 265 - 279. doi:10.1007 / bf00705478
//Maidment, D. R. (Ed.) (1993). Handbook of Hydrology. New York: McGraw-Hill.
//Shuttleworth, J. W. (1993). Evaporation. In D. R. Maidment (Ed.), Handbook of Hydrology (pp. 4.1-4.5.3). New York: McGraw-Hill.
//Stull, R. B. (2000). Meteorology for Scientists and Engineers (2nd ed.). Pacific Grove, CA: Brooks/Cole.
//Su, Z., Schmugge, T., Kustas, W. P., & Massman, W. J. (2001). An evaluation of two models for estimation of the roughness height for heat transfer between the land surface and the atmosphere. Journal of Applied Meteorology, 40(11), 1933-1951.  
//Yang, Y., Endreny, T.A., Nowak, D.J. (2013).A physically based analytical spatial air temperatureand humidity model.Journal of Geophysical Research - Atmospheres, 118(18), 10449 - 10463. doi:10.1002 / jgrd.50803
//Vandegriend, A. A., & Owe, M. (1994). Bare Soil Surface-Resistance to Evaporation by Vapor Diffusion Under Semiarid Conditions. Water Resources Research, 30(2), 181-188. doi:10.1029/93wr02747

//ResistanceCal constructor of the ResistanceCal class initializes variables shared across void functions, shared globally
ResistanceCal::ResistanceCal(Inputs* input)
{
	//Soil_SaturationPoint_m3pm3 (m3/m3) is soil volumetric saturation from DataFolder
	Soil_SaturationPoint_m3pm3 = input->InputXml["Soil_SaturationPoint_m3pm3"];
	//Soil_WiltingPoint_m3pm3 (m3/m3) is soil volumetric wilting point from DataFolder
	Soil_WiltingPoint_m3pm3 = input->InputXml["Soil_WiltingPoint_m3pm3"];
	//Soil_FieldCapacity_m3pm3 (m3/m3) is soil volumetric field capacity from DataFolder
	Soil_FieldCapacity_m3pm3 = input->InputXml["Soil_FieldCapacity_m3pm3"];
	
	//Coeff_SurfaceResistance_Tree_s_p_m (s/m) is canopy minimum resistance for tree
	Coeff_SurfaceResistance_Tree_s_p_m = input->TemperatureCalculationParams["Resistance_Tree_Min_spm"];
	//Coeff_SurfaceResistance_SVeg_s_p_m (s/m) is canopy minimum resistance for short vegetation
	Coeff_SurfaceResistance_SVeg_s_p_m = input->TemperatureCalculationParams["Resistance_SVeg_Min_spm"];
	//Height_WindMeasurement_m (m) is reference height for wind speed measurements
	Height_WindMeasurement_m = input->SimulationNumericalParams["Height_Sensor_Wind_m"];
	//Leaf_Width_Tree_m (m) is tree leaf width used in resistance calculations
	Leaf_Width_Tree_m = input->TemperatureCalculationParams["Width_Leaf_Tree_m"];
	//Leaf_Width_SVeg_m (m) is short vegetation leaf width used in resistance calculations
	Leaf_Width_SVeg_m = input->TemperatureCalculationParams["Width_Leaf_SVeg_m"];
	
}

//ResistanceCal::AerodynamicResistance_CanopyToMesoLayer function computes aerodynamic resistances (s/m) for between canopy layer and mesoscale layer
//Note: Presumes the canopy layer is completely mixed, allowing one average height to begin transmission to mesoscale layer at 50 to 100 m height 
//Note: Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m = ra term in Eq 1 and 2 of Yang et al. (2013), which is estimated w Eq 33 by wind profile at met station
//Note: According to Lee and Park (2008) the overlying urban roughness sublayer (URSL) extends up 50 to 100 m
//Note: Sensitivity of variable: increasing magnitude from order of 10s to 100s introduces large swings in local canopy air temperature
void ResistanceCal::AerodynamicResistance_CanopyToMesoLayer(Inputs* input, CompactRagged* beC, DataOrganizer* organizer, WeatherProcessor* WeatherPro, int MapPixel_ID, int DataFolder_ID, double WindSpeed_reference_m_p_s, double Height_LandCover_Average_m, double Radiation_Shortwave_Direct_W_p_m2, double Radiation_Shortwave_Diffuse_W_p_m2, int timeStep)
{
	/*cte Algorithm Tested but not appropriate for resistance calculation function
	//HeatMetrics_Calc* HeatMetricsCalc to create pointer for accessing variables
	HeatMetrics_Calc HeatMetricsCalc(input);

	//Delta_Height_Tair_m (m) is set to zero (0) by Liljegren et al. (2008) code
	double Delta_Height_Tair_m = 0;
	//Flag_urban is set to 1
	int Flag_urban = 1;
	//Flag_urban is initialized to 1
	int Flag_daytime = 1;
	//Height_WindSpeed_Actual_m (m) is SimulationNumericalParams["Height_Sensor_Wind_m"], reference height for wind speed measurements
	double Height_WindSpeed_Actual_m = input->SimulationNumericalParams["Height_Sensor_Wind_m"];
	//Height_WindSpeed_Target_m = 2 * Height_WindSpeed_Actual_m to obtain stability and speed of air at 2 times station height
	//Note: 2 * Height_WindSpeed_Actual_m is an arbitrary value, if 4 * Height_WindSpeed_Actual_m can lead to doubling wind speed
	double Height_WindSpeed_Target_m = 2 * Height_WindSpeed_Actual_m;
	//ZenithAngle_Solar_rad = beC->by_key(MapPixel_ID, DataFolder_ID, "ZenithAngle_Solar_rad"), will be zero for 1st time step
	//Note: Consider refactor to obtain correct value for 1st time step
	double ZenithAngle_Solar_rad = beC->by_key(MapPixel_ID, DataFolder_ID, "ZenithAngle_Solar_rad");
	//Radiation_Shortwave_Wpm2 is Radiation_Shortwave_Direct_W_p_m2 + Radiation_Shortwave_Diffuse_W_p_m2, combining direct and diffuse
	double Radiation_Shortwave_Wpm2 = Radiation_Shortwave_Direct_W_p_m2 + Radiation_Shortwave_Diffuse_W_p_m2;

	//LandCover_NLCD_Class (Class) defined with input->LandCover_NLCD_Class for MapPixel_ID
	int LandCover_NLCD_Class = input->LandCover_NLCD_Class[MapPixel_ID];
	//if (LandCover_NLCD_Class < 21 || LandCover_NLCD_Class > 24) then LandCover_NLCD_Class not urban
	if (LandCover_NLCD_Class < 21 || LandCover_NLCD_Class > 24) {
		//Flag_urban = 0 and rural wind conditions simulated
		Flag_urban = 0;
	}

	//Cos_ZenithAngle_Solar is cos(ZenithAngle_Solar_rad) for use with direct radiation albedo
	double Cos_ZenithAngle_Solar = cos(ZenithAngle_Solar_rad);
	//Flag_daytime set to True if Cos_ZenithAngle_Solar > 0, or False otherwise
	Flag_daytime = (Cos_ZenithAngle_Solar > 0) ? 1 : 0;
	//StabilityClass_Atmosphere_int is determined from function HeatMetrics_Calc::StabilityClass_Estimator
	int StabilityClass_Atmosphere_int = HeatMetricsCalc.StabilityClass_Estimator(Flag_daytime, WindSpeed_reference_m_p_s, Radiation_Shortwave_Wpm2, Delta_Height_Tair_m);
	//WindSpeed_canopy_to_mesolayer_m_p_s (m/s) computed in HeatMetrics_Calc::WindSpeed_Estimate_mps_Calc 
	//Note: Adjust from measured to resistance height
	double WindSpeed_canopy_to_mesolayer_m_p_s = HeatMetricsCalc.WindSpeed_Estimate_mps_Calc(WindSpeed_reference_m_p_s, Height_WindSpeed_Actual_m, Height_WindSpeed_Target_m, StabilityClass_Atmosphere_int, Flag_urban);
	*/

	//Note: Theory discussed with Dr. Christopher Holst of Kalsruhe Institute of Technology (KIT) Alpine campus in GaPa, developer of PALM-4U model
	//Note: Dr. Holst and post-doc ChangXing Lan (UW PhD) provided insights May-June 2023 at KIT that helped parameterize this energy balance model
	//Note: They assured me that this simple parameterization is as defensible as other approaches given the tremendous uncertainty in the processes
	//Note: Yang et al. (2013) below Eq 33 states no generic formula for the roughness length and zero plane displacement for all landscapes
	//Note: Experimental roughness lengths range 0.4-0.7 m for dense low buildings or 0.7-1.5 m for regularly built towns [Wieringa, 1993; Masson, 2000], 
	//Note: Value approximates one tenth of the building height, or a weighted average landscape roughness element height 
	//Note: Advanced users can specify the values of roughness length and zero plane displacement based on their knowledge of their study area.
	//Note: PASATH Eq 33: Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m = log((Zu- Zd_m) / Z0m_m)*log((Zu- Zd_m) / Z0m_m) / (k*k*WindSpeed_m_p_s);
	//Note: For PASATH model, it hardcoded Eq 33 with parameters: Zu = 10 m, H = 5 m, Zd_m = 4 = 0.8*H, Z0m_m = 0.5 = 0.1*H
	//Note: Height_WindSpeedMeasurement_Standard_m & Height_Landscape_ConstantParameter_m remain constant across all pixels to preserve continuity
	//Note: Height_Landscape_ConstantParameter_m has large influence on value of Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m and Tair_K
	//Note: Height_Landscape_ConstantParameter_m might be adjusted as parameter useful in representing different scenarios and temperature targets

	//Zd_m (m) zero plane displacement height defined as 0.8 * Height_Landscape_ConstantParameter_m, in Yang et al. (2013) for Eq 33
	//Note: Similar to Zd_m = 0.67 * Height_Landscape_ConstantParameter_m used in Eq 4.2.25 of Shuttleworth (1993)
	double Zd_m = 0.8 * Height_Landscape_ConstantParameter_m;
	//Zom_m (m) roughness length for momentum transfer defined as Z0m_m = 0.1 * Height_Landscape_ConstantParameter_m, in Yang et al. (2013) Eq 33
	//Note: Similar to Zom_m = 0.123 * Height_Landscape_ConstantParameter_m used in Eq 4.2.25 of Shuttleworth (1993)
	double Zom_m = 0.1 * Height_Landscape_ConstantParameter_m;

	//WindSpeed_canopy_to_mesolayer_m_p_s (m/s) is presumed equal WindSpeed_reference_m_p_s (m/s) at the reference station, above the urban canopy 
	double WindSpeed_canopy_to_mesolayer_m_p_s = WindSpeed_reference_m_p_s;
	//WindSpeed_minimum_canopy_to_mesolayer_m_p_s (m/s) is the minimum wind speed in this layer which ensures reasonable resistance values
	//double WindSpeed_minimum_canopy_to_mesolayer_m_p_s = 1.75 mps is the minium from 20240530 research with CNR IRET in Italy;
	//Note: Analysis of Verona, Italy for 2003081117 (YYYYMMDDHH) and Torin, Italy for 2003081115 shows 1.75 mps is superior to 1 mps for stability
	double WindSpeed_minimum_canopy_to_mesolayer_m_p_s = 1.75;
	//If WindSpeed_canopy_to_mesolayer_m_p_s < WindSpeed_minimum_canopy_to_mesolayer_m_p_s, then set to modified WindSpeed_minimum_canopy_to_mesolayer_m_p_s
	//Note: WindSpeed_minimum_canopy_to_mesolayer_m_p_s essentially creates a maximum resistance term, and resistance increases as WindSpeed_minimum_canopy_to_mesolayer_m_p_s lowers
	if (WindSpeed_canopy_to_mesolayer_m_p_s < WindSpeed_minimum_canopy_to_mesolayer_m_p_s) { 
		//WindSpeed_canopy_to_mesolayer_m_p_s (m/s) = WindSpeed_minimum_canopy_to_mesolayer_m_p_s + 1/100th of windspeed to obtain unique signature
		WindSpeed_canopy_to_mesolayer_m_p_s = WindSpeed_minimum_canopy_to_mesolayer_m_p_s + WindSpeed_canopy_to_mesolayer_m_p_s / 100.0;
	}

	//Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m (s/m) from Eq 33 of Yang et al. (2013), Equivalently calculated by Eq 1 of Lindroth (1983) 
	//Note: Variables Z_om = Z_ov, Eq 4.2.25 of Maidment (1993) or Shuttleworth (1993)
	//Note: Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m is the same as ra variable used in Eq 1 and Eq 2 of Yang et al. (2013), 
	//Note: Represents resistance between urban canopy layer to the upper mesoscale air layer
	//Note: C++ log function returns natural logarithm
	Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m = (log((Height_WindSpeedMeasurement_Standard_m - Zd_m) / Zom_m) * log((Height_WindSpeedMeasurement_Standard_m - Zd_m) / Zom_m)) / (k_vonKarman * k_vonKarman * WindSpeed_canopy_to_mesolayer_m_p_s);

	//Obtain a constant Tair_K value, across all pixels, to compare against Tair_mesoScale_K
	double Tair_K;
	double Tair_mesoScale_K;
	//If Flag_MultipleStations is true then use Tair_weighted_K and Tair_mesoScale_weighted_K as weighted average of all reference stations
	if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
		Tair_K = WeatherPro->Tair_weighted_K;
		Tair_mesoScale_K = WeatherPro->Tair_mesoScale_weighted_K;
	}
	//Else use Tair_K and Tair_mesoScale_K from reference station
	else {
		//Ref_MapPixel_ID = input->MapPixel_ReferenceStation_IDs[0]; Ref_MapPixel_ID updated to 1st index when Flag_MultipleStations = false
		int Ref_MapPixel_ID = input->MapPixel_ReferenceStation_IDs[0];
		//Ref_DataFolder_ID = 0;
		int Ref_DataFolder_ID = 0;

		Tair_K = beC->by_key(Ref_MapPixel_ID, Ref_DataFolder_ID, "Tair_K");
		Tair_mesoScale_K = beC->by_key(Ref_MapPixel_ID, Ref_DataFolder_ID, "Tair_mesoScale_K");
	}
	double Coeff_ThermalExpansion_p_K = 5;
	double Coeff_LengthCharacteristic_m = Height_WindSpeedMeasurement_Standard_m - Zd_m;
	double StabilityNumber_CM = 0;
	//If timeStep > 0 then compute StabilityNumber_CM, which depends on estimate of mesoScale and local grid temperatures; otherwise it remains zero
	if (timeStep > 0) {
		//StabilityNumber_CM from Eq 19 of Choudhury and Monteith (1988), equals 0 when timeStep=0, Coeff_ThermalExpansion_p_K and Coeff_LengthCharacteristic_m defined above
		//Note: StabilityNumber_CM is only computed for reference station and presumed same for all modeled area
		StabilityNumber_CM = Coeff_ThermalExpansion_p_K * g_GravitationalAcceleration_m_p_s * Coeff_LengthCharacteristic_m * (Tair_K - Tair_mesoScale_K) / (Tair_K * WindSpeed_canopy_to_mesolayer_m_p_s * WindSpeed_canopy_to_mesolayer_m_p_s);
	}

	//StabilityNumber_CM_Max = 10, defined for numerical stability; otherwise resistances become too small
	//Note: Perhaps relax constraint on StabilityNumber_CM range; consider approach of other models
	double StabilityNumber_CM_Max = 10;
	//StabilityNumber_CM_Min = -0.01, defined for numerical stability; otherwise resistances become too large
	//Note: If StabilityNumber_CM is -1.0, resistance becomes infinite in 1st form of Eq 18 of Choudhury and Monteith (1988) 
	double StabilityNumber_CM_Min = -0.01;
	//If StabilityNumber_CM exceeds Max or Min then redefine
	if (StabilityNumber_CM > StabilityNumber_CM_Max) { StabilityNumber_CM = StabilityNumber_CM_Max; }
	if (StabilityNumber_CM < StabilityNumber_CM_Min) { StabilityNumber_CM = StabilityNumber_CM_Min; }

	//If StabilityNumber_CM < 0 then recompute Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m (s/m) 
	//Note: Use 1st form of Eq 18 of Choudhury and Monteith (1988) 
	if (StabilityNumber_CM < 0) {
		//Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m (s/m) from 1st form of Eq 18 of Choudhury and Monteith (1988) 
		Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m = Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m / pow((1 + StabilityNumber_CM), 2);
	}
	//Else If StabilityNumber_CM > 0 then recompute Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m (s/m)
	//Note: Use 2nd form of Eq 18 of Choudhury and Monteith (1988) 
	else if (StabilityNumber_CM >= 0) {
		//Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m (s/m) from 2nd form of Eq 18 of Choudhury and Monteith (1988) 
		Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m = Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m / pow((1 + StabilityNumber_CM), 0.75);
	}

	//Resistance_Aerodynamic_CanopyToMesoLayer_Min_s_p_m (s/m) defined based on numerical stability to a value of 0.5
	double Resistance_Aerodynamic_CanopyToMesoLayer_Min_s_p_m = 0.5;
	double Resistance_Aerodynamic_CanopyToMesoLayer_Max_s_p_m = 5000;
	//Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m is clamped by Resistance_Aerodynamic_CanopyToMesoLayer_Min_s_p_m and _Max_s_p_m
	Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m = clamp(Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m,
			Resistance_Aerodynamic_CanopyToMesoLayer_Min_s_p_m,
			Resistance_Aerodynamic_CanopyToMesoLayer_Max_s_p_m);

	//beC->by_key(MapPixel_ID, DataFolder_ID, "WindSpeed_canopy_to_mesolayer_m_p_s") = WindSpeed_canopy_to_mesolayer_m_p_s; saved for output
	beC->by_key(MapPixel_ID, DataFolder_ID, "WindSpeed_canopy_to_mesolayer_m_p_s") = WindSpeed_canopy_to_mesolayer_m_p_s;
	//beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m") = Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m; saved for internal calculations
	//Note: Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m saved in HeatFluxCal::CollectVariables_CoolAir after updated for actual conditions
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m") = Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m;
}

//ResistanceCal::SurfaceResistance_LandCover function computes surface resistances (s/m) for soil moisture based land cover
void ResistanceCal::SurfaceResistance_LandCover(CompactRagged* beC, int MapPixel_ID, int DataFolder_ID , double Saturation_SoilStorage_frac, double LAI_Tree_m2pm2, double LAI_SVeg_m2pm2)
{
	//Resistance_Surface_Soil_s_p_m (s/m) from Eq 34 of Yang et al. (2013), which cites Sellers et al. (1992) for coefficients 8.206 and 4.255, but this is not evident in the publication
	//Note: Vandegriend and Owe (1994) w Eq 19 and 20 explain origin of exponential form of soil resistance equation, stating it should ...
	//Note: ... have minimum of 0 s/m when saturated and then increase exponentially toward 2000 s/m with drying (see Figure 7).
	Resistance_Surface_Soil_s_p_m = exp(8.206 - 4.255 * (Saturation_SoilStorage_frac * (Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3)) / (Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3)); 

	//Saturation_SoilStorage_zero_frac defined as zero to find maximum likely Resistance_Surface_Soil_s_p_m
	double Saturation_SoilStorage_zero_frac = 0;
	//Resistance_Surface_Soil_Max_s_p_m defined for Saturation_SoilStorage_zero_frac condition
	double Resistance_Surface_Soil_Max_s_p_m = exp(8.206 - 4.255 * (Saturation_SoilStorage_zero_frac * (Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3)) /	(Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3));

	//Resistance_Surface_Impervious_s_p_m (s/m) defined to equal Resistance_Surface_Soil_s_p_m 
	Resistance_Surface_Impervious_s_p_m = Resistance_Surface_Soil_s_p_m;

	//Resistance_Surface_Water_s_p_m (s/m) defined as 0, which is the case for water surface, defined by Chin (2021) after Eq 13.12
	//Note: Resistance_Surface_Water_s_p_m can be used for pervious and impervious depression storage as well as water cover
	Resistance_Surface_Water_s_p_m = 1;

	//Option for calculation of Resistance_Surface_Soil_s_p_m using Eq 20 of Vandegriend and Owe (1994)
	double Saturation_SoilStorage_percent = Saturation_SoilStorage_frac * ratio_percent_to_decimal;
	double Coeff_VO_Eq20_a = 10;
	double Coeff_VO_Eq20_b = 0.3563;
	double Coeff_VO_Eq20_c = 15;
	//Note: Vandegriend and Owe (1994) w Eq 19 and 20 explain origin of exponential form of soil resistance equation, rs = 0 when saturated, then increase exponentially with drying, 
	//Note: Vandegriend and Owe (1994) provide Eq 20 as: rs = 10 * exp [0.3563 * (15 - Soil Moisture Volume_%)], which is for top 1 cm of soil
	//Note: Slight modification of V&O 1994 Eq 20 gives: rs = 10 * exp [ 0.3563 * (15 - (Saturation_SoilStorage_percent * Soil_SaturationPoint_m3pm3))]
	//Resistance_Surface_Soil_s_p_m = Coeff_VO_Eq20_a * exp(Coeff_VO_Eq20_b * (Coeff_VO_Eq20_c - (Saturation_SoilStorage_percent * Soil_SaturationPoint_m3pm3);

	//Resistance_Surface_Tree_Max_s_p_m (s/m) for when soil water content reaches the wilting point, Saturation_SoilStorage_zero_frac = 0 
	//Note: Setting Resistance_Surface_Tree_Max_s_p_m = 95% of Resistance_Surface_Soil_Max_s_p_m
	//Note: Setting Resistance_Surface_SVeg_Max_s_p_m = 98% of Resistance_Surface_Soil_Max_s_p_m
	//Note: From 2020 to circa fall 2024 Resistance_Surface_Tree_Max_s_p_m was 2000 s/m based on Figure 7 of Vandegriend and Owe (1994)
	//Note: From 2013 to circa 2019 was Resistance_Surface_Tree_Max_s_p_m was 5000 s/m, set by Yang et al. (2013) in PASATH, w/o reference
	//Note: Variation could be added between tree and short veg given height to pull water
	double Resistance_Surface_Tree_Max_s_p_m = Resistance_Surface_Soil_Max_s_p_m * 0.95;
	double Resistance_Surface_SVeg_Max_s_p_m = Resistance_Surface_Soil_Max_s_p_m * 0.98;

	//If Saturation_SoilStorage_frac > 0 then water is available for transpiration
	if (Inputs::isGreaterThan(Saturation_SoilStorage_frac, 0.0)) {

		//If LAI_Tree_m2pm2 <= 1 then set to minimum 1 to avoid division by values less than 1
		if (LAI_Tree_m2pm2 <= 1) { LAI_Tree_m2pm2 = 1; }
		//If LAI_SVeg_m2pm2 <= zero then set to minimum 1 to avoid division by values less than 1
		if (LAI_SVeg_m2pm2 <= 1) { LAI_SVeg_m2pm2 = 1; }

		//Resistance_Surface_Tree_s_p_m (s/m) from Eq 35 of Yang et al. (2013), where RHS is g, and g^-1 = (Theta-WP)/(FC-WP)
		//Note: Eq 35 sets Theta to maximum to Soil_SaturationPoint_m3pm3, not FC, to reduce resistances along a drying tail of curve
		//Note: Eq 35 cites Eq 6 of Liang et al. (1994) VIC model, provides flux resistance from tree root to tree canopy
		//Note: Term minimum resistance is misleading for Eq 35, as the function allows resistance to go below that value ...
		//Note: Again, Eq 35 allows Resistance_Surface_Tree_s_p_m to go below Coeff_SurfaceResistance_Tree_s_p_m 
		//Note: e.g., Reports of minimum resistance of 40 s/m are also found in Chpt 11 of Jensen and Allen (2015) 
		//Note: Coeff_SurfaceResistance_Tree_s_p_m = 100 s/m in Table 2 of Liang et al. (1994) VIC model and ...
		//Note: ... Box 5 of Allen et al. (1998) Chapter 2 - FAO Penman-Monteith equation https://www.fao.org/3/x0490e/x0490e06.htm
		//Note: Richard Allen has pointed out Shuttleworth (1993) Eq 4.2.22 is not typically representative of trees ...
		//Note: ... useful for well watered crops that are bred to maximize production
		//Note: ... Trees have evolved prioritize survival over production, and will typically have reduced surface resistance.
		//Note: ... Trees are exposed to stronger winds, have greater roughness, and often less secure water supplies. 
		//Note: Inputs::SafeDivide not needed due to function ensuring Saturation_SoilStorage_frac & LAI_SVeg_m2pm2 > 0
		Resistance_Surface_Tree_s_p_m = Coeff_SurfaceResistance_Tree_s_p_m * (Soil_FieldCapacity_m3pm3 - Soil_WiltingPoint_m3pm3) / (Saturation_SoilStorage_frac * (Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3)) / LAI_Tree_m2pm2; 

		//Resistance_Surface_Tree_s_p_m (s/m) is increased up to two times its value as function of ImperviousCover_frac
		//Note: Theory and equation from T. Endreny in 2022, capturing the phenommena IC reducing soil porosity around roots
		//Note: Effect of IC demonstrated in data collected by V. Shandas for Richmond, VA heatwave of July 13, 2017 
		//Note: ImperviousCover_frac and not TreeCanopyCover_overImpervious_frac is the preferred scaling factor
		//Note: Function used ImperviousCover_frac prior to allowing IC below TC within Cool Air for cases when IC + TC < 100%
		//Resistance_Surface_Tree_s_p_m = Resistance_Surface_Tree_s_p_m * (1 + beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousCover_frac"));
		Resistance_Surface_Tree_s_p_m = Resistance_Surface_Tree_s_p_m * (1 + beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousCover_frac"));

		//Resistance_Surface_Tree_s_p_m (s/m) from Eq 35 of Yang et al. (2013), where RHS is g, and g^-1 = (Theta-WP)/(FC-WP)
		//Note: Eq 35 sets Theta to maximum to Soil_SaturationPoint_m3pm3, not FC, to reduce resistances along a drying tail of curve
		//Note: Eq 35 cites Eq 6 of Liang et al. (1994) VIC model, provides flux resistance from tree root to tree canopy
		//Note: Term minimum resistance is misleading for Eq 35, as the function allows resistance to go below that value ...
		//Note: Again, Eq 35 allows Resistance_Surface_SVeg_s_p_m to go below Coeff_SurfaceResistance_SVeg_s_p_m 
		//Note: Coeff_SurfaceResistance_SVeg_s_p_m = 100 s/m in Table 2 of Liang et al. (1994) VIC model and ...
		//Note: ... Box 5 of Allen et al. (1998) Chapter 2 - FAO Penman-Monteith equation
		//Note: Inputs::SafeDivide not needed due to function ensuring Saturation_SoilStorage_frac & LAI_SVeg_m2pm2 > 0
		Resistance_Surface_SVeg_s_p_m = Coeff_SurfaceResistance_SVeg_s_p_m * (Soil_FieldCapacity_m3pm3 - Soil_WiltingPoint_m3pm3) / (Saturation_SoilStorage_frac * (Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3)) / LAI_SVeg_m2pm2;
	}
	//Else If no soil water is available then set to maximum values
	else {
		//Resistance_Surface_Soil_s_p_m set to maximum of Resistance_Surface_Soil_Max_s_p_m
		Resistance_Surface_Soil_s_p_m = Resistance_Surface_Soil_Max_s_p_m;
		//Resistance_Surface_Tree_s_p_m set to maximum of Resistance_Surface_Tree_Max_s_p_m
		Resistance_Surface_Tree_s_p_m = Resistance_Surface_Tree_Max_s_p_m;
		//Resistance_Surface_SVeg_s_p_m set to maximum of Resistance_Surface_SVeg_Max_s_p_m
		Resistance_Surface_SVeg_s_p_m = Resistance_Surface_SVeg_Max_s_p_m;
	}

	//Resistance_Surface_Tree_s_p_m = clamp(Resistance_Surface_Tree_s_p_m, 0, Resistance_Surface_Tree_Max_s_p_m)
	//Note: Resistance_Surface_Tree_s_p_m cannot go below 0, but can go below Resistance_Surface_Tree_Min_s_p_m
	//Note: Again, Eq 35 allows Resistance_Surface_Tree_s_p_m to go below Coeff_SurfaceResistance_Tree_s_p_m 
	Resistance_Surface_Tree_s_p_m = clamp(Resistance_Surface_Tree_s_p_m, 0.0, Resistance_Surface_Tree_Max_s_p_m);

	//Resistance_Surface_SVeg_s_p_m = clamp(Resistance_Surface_SVeg_s_p_m, 0, Resistance_Surface_SVeg_Max_s_p_m)
	//Note: Resistance_Surface_SVeg_s_p_m cannot go below 0, but can go below Resistance_Surface_SVeg_Min_s_p_m
	//Note: Again, Eq 35 allows Resistance_Surface_SVeg_s_p_m to go below Coeff_SurfaceResistance_SVeg_s_p_m 
	Resistance_Surface_SVeg_s_p_m = clamp(Resistance_Surface_SVeg_s_p_m, 0.0, Resistance_Surface_SVeg_Max_s_p_m);

	//Save resistance values for output
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Surface_Tree_s_p_m") = Resistance_Surface_Tree_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Surface_SVeg_s_p_m") = Resistance_Surface_SVeg_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Surface_Soil_s_p_m") = Resistance_Surface_Soil_s_p_m;
}

//ResistanceCal::AerodynamicResistance_UrbanCanyon function computes aerodynamic resistances (s/m) for urban canyon area, using aspect ratio
//Note: Inspired by Lee and Park (2008) with Vegetated Urban Canopy Model and Fickian based heat flux with thermal conductivity and temperature gradient, e.g., Eq 2
//Note: Aerodynamic resistancs may need adjustment away from VUCM, as they tend to artificially constrain the minimum temperatures
void ResistanceCal::AerodynamicResistance_UrbanCanyon(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, double WindSpeed_reference_m_p_s, double LAI_Tree_m2pm2, double LAI_SVeg_m2pm2, double AspectRatio_Canyon_frac, double Height_LandCover_Average_m, double SpecificHeat_HumidAir_JpkgK)
{
	//Zd_m (m) zero plane displacement height defined as 0.67*H, see Eq 4.2.25 Shuttleworth (1993)
	double Zd_m = 0.67 * Height_LandCover_Average_m;
	//Zom_m (m) roughness length for momentum transfer defined as 0.123*H, see Eq 4.2.25 Shuttleworth (1993)
	double Zom_m = 0.123 * Height_LandCover_Average_m;
	//WindSpeed_folder_m_p_s (m/s) computed w/ Eq 4.14b of Stull (2000) Wspd_height_b = Wspd_height_a * [log (height_b / roughnesslength_b) / log (height_a / roughnesslength_a)]
	//Note: WindSpeed_folder_m_p_s is very sensitive to value of RoughnessLength_ImpSoil_m, which PASATH set to 0.00137 m
	double WindSpeed_folder_m_p_s = WindSpeed_reference_m_p_s * (log(Height_LandCover_Average_m / RoughnessLength_ImpSoil_m) / log(Height_WindMeasurement_m / RoughnessLength_ImpSoil_m));

	//WindSpeed_minimum_UrbanCanyon_m_p_s (m/s) is the minimum wind speed in this layer which ensures reasonable resistance values
	//Note: Criteria suggested by Lee and Park (2008) in Eq 53, using different vegetation resistance formulation for windspeed < 0.25
	//Note: Yang et al (2013) in PASATH modified Lee and Park (2008) to simply set WindSpeed_minimum_UrbanCanyon_m_p_s rather than use different formulation
	double WindSpeed_minimum_UrbanCanyon_m_p_s = 0.25;
	//If WindSpeed_folder_m_p_s (m/s) < WindSpeed_minimum_UrbanCanyon_m_p_s or 0.25 m/s then set to WindSpeed_minimum_UrbanCanyon_m_p_s or 0.25 m/s
	if (WindSpeed_folder_m_p_s < WindSpeed_minimum_UrbanCanyon_m_p_s) { WindSpeed_folder_m_p_s = WindSpeed_minimum_UrbanCanyon_m_p_s; }

	//WindSpeed_Canyon_m_p_s (m/s) from Eq 47 of Lee and Park (2008) for wind profile in urban canyon was simplified, dropping second RHS term  * exp (0.1 * PerviousArea_frac * VegetationArea_frac * LAI_BAI)
	//Note: Lee and Park (2008) Eq 48 defines first RHS term, WindSpeed_folder_m_p_s, as U_hb, wind speed at building top height, using theory of Hogstrom (1988)
	//WindSpeed_Canyon_m_p_s (m/s) is calculated the windspeed within the canopy based on height to width ratio
	//Note: Consider refactor to improve representation with Lee and Park (2008), fixing WindSpeed_folder_m_p_s to U_hb, and adding last RHS term
	double WindSpeed_Canyon_m_p_s = WindSpeed_folder_m_p_s * exp(-0.386 * AspectRatio_Canyon_frac);
	//Note: AspectRatio_Canyon_frac or height (h) to width (w) variable, which ranges from AspectRatio_Canyon_frac=0.01 for LC=21 to AspectRatio_Canyon_frac=1.2 for LC=24. 
	//Note: Based on Table 2 of Yang et al. (2013), LC=21 had AspectRatio_Canyon_frac=0.01 and h=0.1 and w=10, but the PASATH HeatCal.cpp code for LC=21 used h=1 and AspectRatio_Canyon_frac=0. 
	//Note: NLCD Class 21 AspectRatio_Canyon_frac=0.01 -> WindSpeed_Canyon_m_p_s=99.6% of WindSpeed_folder_m_p_s; Option for deeper canyon WindSpeed_Canyon_m_p_s=96.2% of WindSpeed_folder_m_p_s when AspectRatio_Canyon_frac=0.1
	//Note: NLCD Class 22 AspectRatio_Canyon_frac=0.50 -> WindSpeed_Canyon_m_p_s=82.4% of WindSpeed_folder_m_p_s (shallow urban canyon)
	//Note: NLCD Class 23 AspectRatio_Canyon_frac=0.80 -> WindSpeed_Canyon_m_p_s=73.4% of WindSpeed_folder_m_p_s
	//Note: NLCD Class 23 AspectRatio_Canyon_frac=1.20 -> WindSpeed_Canyon_m_p_s=62.9% of WindSpeed_folder_m_p_s; Option for deeper canyon and WindSpeed_Canyon_m_p_s=42.7% of WindSpeed_folder_m_p_s when AspectRatio_Canyon_frac=2.2 (deep urban canyon)

	//Resistance_Aerodynamic_Impervious_s_p_m (s/m) from Eq 31 Yang et al. (2013) from Eq 52 Lee and Park (2008), which computes canyon heat transfer coefficient, h_C = 11.8 * 4.2 * UC [W/m2/K]
	//Note: Eq 31 Yang et al. (2013) has units s/m, Density_Air_kg_p_m3 * SpecificHeat_DryAir_J_p_kg_K = [kg/m3] * [W/s/kg/K] = [W/s/K/m3], divided by h_C or [W/m2/K] = [s/m] or 1/WindSpeed_Canyon_m_p_s
	//Note: Eq 52 of Lee and Park (2008) for h_C, heat transfer coefficient enters Eq 51 for H_wall, wall sensible heat flux, used with Canyon wall and air temperature gradient
	//Note: Lee and Park (2008) use H_wall in implicit Eq 15 to find Canyon air temperature as function of sensibile heat (H) terms H_wall, H_ground, H_vegetation, H_anthropogenic, ...
	//Note: ... and H_Canyon which depend on Canyon air temperature. Eq 15 also includes Canyon air volume Delta_V_C, and Canyon bottom area, Delta_A_C
	//Resistance_Aerodynamic_Impervious_s_p_m (s/m) is flux resistance from impervious surface to folder canopy level
	Resistance_Aerodynamic_Impervious_s_p_m = Density_Air_kg_p_m3 * SpecificHeat_HumidAir_JpkgK / (11.48 + 4.2 * WindSpeed_Canyon_m_p_s);

	//Eq 32 Yang et al. (2013) from Eq 53 and 54 and LAI* of Lee and Park (2008); Yang used coefficient 12, Lee and Park used coefficient ratio 30/2.5, 30/2.5 = 12
	//Note: Eq 53 and 54 of Lee and Park (2008) combine to give rb = 30 * (1 + 0.55 * LAI_Tree_m2pm2) / sqrt(WindSpeed_Canyon_m_p_s) / LAI*, where ...
	//Note: ... LAI* defined as effective sunlet LAI by Lee and Park (2008) below Eq 17, LAI* = 2.5 * (1 - exp(-0.4 * LAI_Tree_m2pm2)
	//Note: Coefficient 0.4 = K of Eq 1 in Kjelgren and Montague (1998), defined as crown transmissivity to light. Potentially refactor K to vary.
	//Resistance_Aerodynamic_Tree_s_p_m (s/m) flux resistance from tree canopy to folder canopy level
	Resistance_Aerodynamic_Tree_s_p_m = 30 * (1 + 0.55 * LAI_Tree_m2pm2) / sqrt(WindSpeed_Canyon_m_p_s) / (2.5 * (1 - exp(-0.4 * LAI_Tree_m2pm2)));

	//Eq 32 Yang et al. (2013) from Eq 53 and 54 and LAI* of Lee and Park (2008); Yang used coefficient 12, Lee and Park used coefficient ratio 30/2.5, 30/2.5 = 12
	//Note: Eq 53 and 54 of Lee and Park (2008) combine to give rb = 30 * (1 + 0.55 * LAI_Tree_m2pm2) / sqrt(WindSpeed_Canyon_m_p_s) / LAI*, where ...
	//Note: ... LAI* defined as effective sunlet LAI by Lee and Park (2008) below Eq 17, LAI* = 2.5 * (1 - exp(-0.4 * LAI_Tree_m2pm2)
	//Note: Coefficient 0.4 = K of Eq 1 in Kjelgren and Montague (1998), defined as crown transmissivity to light. Potentially refactor K to vary.
	//Resistance_Aerodynamic_SVeg_s_p_m (s/m) flux resistance from short vegetation canopy to folder canopy level
	Resistance_Aerodynamic_SVeg_s_p_m = 30 * (1 + 0.55 * LAI_SVeg_m2pm2) / sqrt(WindSpeed_Canyon_m_p_s) / (2.5 * (1 - exp(-0.4 * LAI_SVeg_m2pm2)));	

	//Resistance_Aerodynamic_Soil_s_p_m (s/m) set equal to Resistance_Aerodynamic_Impervious_s_p_m
	Resistance_Aerodynamic_Soil_s_p_m = Resistance_Aerodynamic_Impervious_s_p_m;
	//Resistance_Aerodynamic_Water_s_p_m (s/m) set to Resistance_Aerodynamic_Impervious_s_p_m
	Resistance_Aerodynamic_Water_s_p_m = Resistance_Aerodynamic_Impervious_s_p_m;

	//Save wind speed and resistance values for output
	beC->by_key(MapPixel_ID, DataFolder_ID, "WindSpeed_folder_m_p_s") = WindSpeed_folder_m_p_s;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Tree_s_p_m") = Resistance_Aerodynamic_Tree_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_SVeg_s_p_m") = Resistance_Aerodynamic_SVeg_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Impervious_s_p_m") = Resistance_Aerodynamic_Impervious_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Soil_s_p_m") = Resistance_Aerodynamic_Soil_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Water_s_p_m") = Resistance_Aerodynamic_Water_s_p_m;
}

//ResistanceCal::AerodynamicResistance_OpenArea function computes aerodynamic resistances (s/m) for open area without urban canyons
//Note: Algorithm components from Choudhury and Monteith (1988) as in https://rdrr.io/github/lhmet-forks/bigleaf/src/R/boundary_layer_conductance.r	
void ResistanceCal::AerodynamicResistance_OpenArea(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, double WindSpeed_reference_m_p_s, double LAI_Tree_m2pm2, double LAI_SVeg_m2pm2, double Height_LandCover_Average_m, int timeStep)
{
	//TreeCover_frac set to folder value of TreeCover_frac
	TreeCover_frac = beC->by_key(MapPixel_ID, DataFolder_ID, "TreeCover_frac");
	//ShortVegCover_noTreeCanopy_frac set to folder value of ShortVegCover_noTreeCanopy_frac
	ShortVegCover_noTreeCanopy_frac = beC->by_key(MapPixel_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac");
	//If TreeCover_frac <= zero then set to minimum 0.01 to avoid division by zero; note smaller fraction inadvertently generates larger aerodynamic resistance
	if (TreeCover_frac <= 0) { TreeCover_frac = 0.01; }
	//If ShortVegCover_noTreeCanopy_frac <= zero then set to minimum 0.01 to avoid division by zero; note smaller fraction inadvertently generates larger aerodynamic resistance 
	if (ShortVegCover_noTreeCanopy_frac <= 0) { ShortVegCover_noTreeCanopy_frac = 0.01; }
	//If LAI_Tree_m2pm2 <= zero then set to minimum 0.1 to avoid division by zero
	if (LAI_Tree_m2pm2 <= 0) { LAI_Tree_m2pm2 = 0.1; }
	//If LAI_SVeg_m2pm2 <= zero then set to minimum 0.1 to avoid division by zero 
	if (LAI_SVeg_m2pm2 <= 0) { LAI_SVeg_m2pm2 = 0.1; }
	//If Leaf_Width_Tree_m <= zero then set to minimum 0.05 to avoid division by zero
	if (Leaf_Width_Tree_m <= 0) { Leaf_Width_Tree_m = 0.05; }
	//If Leaf_Width_SVeg_m <= zero then set to minimum 0.01 to avoid division by zero
	if (Leaf_Width_SVeg_m <= 0) { Leaf_Width_SVeg_m = 0.01; }

	//Note: Consider refactor use Zom_m formula of Choudhury and Monteith (1983) Eq 20 and 21, noting Zom_soil_m and Zom_veg_m are distinct
	//Note: Consider refactor to follow Su et al. (2001) and github bigleaf to compute Zoh_m for soil and vegetation as function of kB^-1
	//Zd_m (m) zero plane displacement height defined as 0.67*H, see Eq 4.2.25 Shuttleworth (1993)
	double Zd_m = 0.67 * Height_LandCover_Average_m;
	//Zom_m (m) roughness length for momentum transfer defined as 0.123*H, see Eq 4.2.25 Shuttleworth (1993)
	double Zom_m = 0.123 * Height_LandCover_Average_m;
	//Zom_m (m) roughness length for heat transfer defined as 0.0123*H, see Eq 4.2.25 Shuttleworth (1993)
	double Zoh_m = 0.0123 * Height_LandCover_Average_m;

	//WindSpeed_folder_m_p_s (m/s) computed w/ Eq 4.14b of Stull (2000) Wspd_height_b = Wspd_height_a * [log (height_b / roughnesslength_b) / log (height_a / roughnesslength_a)]
	//Note: Choudbury and Monteith (1983) expect windspeed at height of vegetation or soil or other ground cover, e.g. building
	//Note: RoughnessLength_ImpSoil_m is set as effective roughness length for all land cover combinations to transfer momentum immediately above surface, which PASATH set to 0.00137 m
	//Note: Vortices or turbulence within urban canyons with large aspect ratio and greater actual roughness presumed to completely mix air temperature
	double WindSpeed_folder_m_p_s = WindSpeed_reference_m_p_s * (log(Height_LandCover_Average_m / RoughnessLength_ImpSoil_m) /
		log(Height_WindMeasurement_m / RoughnessLength_Airport_m));

	//WindSpeed_minimum_OpenArea_m_p_s (m/s) is 0.25 m/s, the minimum wind speed in this layer which ensures reasonable resistance values
	double WindSpeed_minimum_OpenArea_m_p_s = 0.25;
	//If WindSpeed_folder_m_p_s (m/s) < WindSpeed_minimum_OpenArea_m_p_s or 0.25 m/s then set to WindSpeed_minimum_OpenArea_m_p_s or 0.25 m/s
	if (WindSpeed_folder_m_p_s < WindSpeed_minimum_OpenArea_m_p_s) { WindSpeed_folder_m_p_s = WindSpeed_minimum_OpenArea_m_p_s; }

	//Note: Algorithm developed using observed sensible and heat flux data from Capodimonte Park of Naples, Italy from Dr. Carlo Calfapietra's team
	//Resistance_Aerodynamic_Soil_s_p_m algorithm from Eq 25 of Choudhury and Monteith (1988), where momentum, heat and water vapor transfer treated equally (see below Eq 25)
	//K_H_soil from Eq 24 of Choudhury and Monteith (1988) as function relating turbulent transfer coefficient, friction velocity, and height; used in Eq 25
	//Note: log is ln in C++, and Eq 24 z = H, so Height_SVeg_m used twice in Eq
	K_H_soil = k_vonKarman * k_vonKarman * (Height_LandCover_Average_m - Zd_m) * WindSpeed_folder_m_p_s / log((Height_LandCover_Average_m - Zd_m) / Zom_m);
	//Coefficient_alpha from Eq 23 of Choudhury and Monteith (1988) as relation between the turbulent transfer coefficient, friction velocity and heigh
	//Note: Choudhury and Monteith (1988) defined as alpha=2 in Table 2; Yang in PASARH defined coefficient alpha as ak_attenuation = 2.25
	Coefficient_alpha = 2.25;
	//Resistance_Aerodynamic_Soil_s_p_m (s/m) from Eq 25 of Choudhury and Monteith (1988), where momentum, heat and water vapor transfer treated equally (see below Eq 25)
	Resistance_Aerodynamic_Soil_s_p_m = Height_LandCover_Average_m * exp(Coefficient_alpha) / (Coefficient_alpha * K_H_soil) *
		(exp(-Coefficient_alpha * RoughnessLength_ImpSoil_m / Height_LandCover_Average_m) - exp(-Coefficient_alpha * (RoughnessLength_ImpSoil_m + Zom_m) / Height_LandCover_Average_m));

	//Resistance_Aerodynamic_Tree_s_p_m algorithm from Eq 29 and 30 of Choudhury and Monteith (1988), where momentum, heat and water vapor transfer treated equally (see below Eq 21)
	//If LAI_Tree_m2pm2 * TreeCover_frac < 1 then ensure Coefficient_a_spm and Coefficient_alphaPrime take minimum values
	if (LAI_Tree_m2pm2 * TreeCover_frac < 1) {
		//Coefficient_a_spm (s/m) defined as 0.01 s/m below Eq 26 in Choudhury and Monteith (1988), now function of LAI and cover fraction ... 
		//Note: Yang in PASATH set as a_attenuation = 0.00662, which is rounded to 0.007 and then increased as function of LAI and cover fraction
		Coefficient_a_spm = 0.007;
		//Coefficient_alphaPrime defined as 3.0 in Table 2 Choudhury and Monteith (1988), now function of LAI and cover fraction ...
		//Note: Github set alphaPrime (denoted alpha) as function of LAI in https://rdrr.io/github/lhmet-forks/bigleaf/man/Gb.Choudhury.html: αlpha = 4.39 - 3.97*exp(-0.258*LAI)
		Coefficient_alphaPrime = 4.39 - 3.97 * exp(-0.258);
	}
	//Else If LAI_Tree_m2pm2 * TreeCover_frac >= 1 then compute Coefficient_a_spm and Coefficient_alphaPrime as function of LAI and cover fraction
	else {
		//Coefficient_a_spm (s/m) defined as 0.01 s/m below Eq 26 in Choudhury and Monteith (1988), now function of LAI and cover fraction ... 
		//Note: Yang in PASATH set as a_attenuation = 0.00662, which is rounded to 0.007 and then increased as function of LAI and cover fraction
		Coefficient_a_spm = 0.007 * LAI_Tree_m2pm2 * TreeCover_frac;
		//Coefficient_alphaPrime defined as 3.0 in Table 2 Choudhury and Monteith (1988), now function of LAI and cover fraction ...
		//Note: Github set alphaPrime (denoted alpha) as function of LAI in https://rdrr.io/github/lhmet-forks/bigleaf/man/Gb.Choudhury.html: αlpha = 4.39 - 3.97*exp(-0.258*LAI)
		Coefficient_alphaPrime = 4.39 - 3.97 * exp(-0.258 * LAI_Tree_m2pm2 * TreeCover_frac);
	}
	//Resistance_Aerodynamic_Tree_s_p_m algorithm from Eq 29 and 30 of Choudhury and Monteith (1988), where momentum, heat and water vapor transfer treated equally (see below Eq 21)
	//Conductance_Aerodynamic_Tree_m_p_s (m/s) from Eq 29 of Choudhury and Monteith (1988) for entire LAI, converting with LAI_Tree_m2pm2 from single leaf
	//Note: WindSpeed_folder_m_p_s (m/s) should be wind speed at canopy surface, where momentum, heat and water vapor transfer treated equally (see below Eq 21)
	//Note: Eq 29 conductance of Choudhury and Monteith (1988) in lhmet-forks/bigleaf is extended using Hicks et al. (1987) to obtain conductance rate for dry deposition 
	Conductance_Aerodynamic_Tree_m_p_s = LAI_Tree_m2pm2 * (2 * Coefficient_a_spm / Coefficient_alphaPrime) * pow((WindSpeed_folder_m_p_s / Leaf_Width_Tree_m), 0.5) *
		(1 - exp(-Coefficient_alphaPrime / 2));

	//Resistance_Aerodynamic_Tree_s_p_m (s/m) from Eq 30 of Choudhury and Monteith (1988) as 1 / Conductance_Aerodynamic_Tree_m_p_s
	Resistance_Aerodynamic_Tree_s_p_m = 1 / Conductance_Aerodynamic_Tree_m_p_s;

	//Resistance_Aerodynamic_SVeg_s_p_m algorithm from Eq 29 and 30 of Choudhury and Monteith (1988), where momentum, heat and water vapor transfer treated equally (see below Eq 21)
	//If LAI_SVeg_m2pm2 * ShortVegCover_noTreeCanopy_frac < 1 then ensure Coefficient_a_spm and Coefficient_alphaPrime take minimum values
	if (LAI_SVeg_m2pm2 * ShortVegCover_noTreeCanopy_frac < 1) {
		//Coefficient_a_spm (s/m) defined as 0.01 s/m below Eq 26 in Choudhury and Monteith (1988), now function of LAI and cover fraction ... 
		//Note: Yang in PASATH set as a_attenuation = 0.00662, which is rounded to 0.007 and then increased as function of LAI and cover fraction
		Coefficient_a_spm = 0.007;
		//Coefficient_alphaPrime defined as 3.0 in Table 2 Choudhury and Monteith (1988), now function of LAI and cover fraction ...
		//Note: Github set alphaPrime (denoted alpha) as function of LAI in https://rdrr.io/github/lhmet-forks/bigleaf/man/Gb.Choudhury.html: αlpha = 4.39 - 3.97*exp(-0.258*LAI)
		Coefficient_alphaPrime = 4.39 - 3.97 * exp(-0.258);
	}
	//Else If LAI_SVeg_m2pm2 * ShortVegCover_noTreeCanopy_frac >= 1 then compute Coefficient_a_spm and Coefficient_alphaPrime as function of LAI and cover fraction
	else {
		//Coefficient_a_spm (s/m) defined as 0.01 s/m below Eq 26 in Choudhury and Monteith (1988), now function of LAI and cover fraction ... 
		//Note: Yang in PASATH set as a_attenuation = 0.00662, which is rounded to 0.007 and then increased as function of LAI and cover fraction
		Coefficient_a_spm = 0.007 * LAI_SVeg_m2pm2 * ShortVegCover_noTreeCanopy_frac;
		//Coefficient_alphaPrime defined as 3.0 in Table 2 Choudhury and Monteith (1988), now function of LAI and cover fraction ...
		//Note: Github set alphaPrime (denoted alpha) as function of LAI in https://rdrr.io/github/lhmet-forks/bigleaf/man/Gb.Choudhury.html: αlpha = 4.39 - 3.97*exp(-0.258*LAI)
		Coefficient_alphaPrime = 4.39 - 3.97 * exp(-0.258 * LAI_SVeg_m2pm2 * ShortVegCover_noTreeCanopy_frac);
	}

	//Conductance_Aerodynamic_SVeg_m_p_s (m/s) from Eq 29 of Choudhury and Monteith (1988) for entire LAI, converting with LAI_SVeg_m2pm2 from single leaf
	//Note: WindSpeed_folder_m_p_s (m/s) should be wind speed at canopy surface, where momentum, heat and water vapor transfer treated equally (see below Eq 21)
	//Note: Eq 29 conductance of Choudhury and Monteith (1988) in lhmet-forks/bigleaf is extended using Hicks et al. (1987) to obtain conductance rate for dry deposition 
	Conductance_Aerodynamic_SVeg_m_p_s = LAI_SVeg_m2pm2 * (2 * Coefficient_a_spm / Coefficient_alphaPrime) * pow((WindSpeed_folder_m_p_s / Leaf_Width_SVeg_m), 0.5) *
		(1 - exp(-Coefficient_alphaPrime / 2));
	//Resistance_Aerodynamic_SVeg_s_p_m (s/m) from Eq 30 of Choudhury and Monteith (1988) as 1 / Conductance_Aerodynamic_SVeg_m_p_s
	Resistance_Aerodynamic_SVeg_s_p_m = 1 / Conductance_Aerodynamic_SVeg_m_p_s;

	//Resistance_Surface_Soil_s_p_m (s/m) from Eq 34 of Yang et al. (2013), which cites Sellers et al. (1992) for coefficients 8.206 and 4.255, but this is not evident in the publication
	//Note: Vandegriend and Owe (1994) w Eq 19 and 20 explain origin of exponential form of soil resistance equation, stating it should ...
	//Note: ... have minimum of 0 s/m when saturated and then increase exponentially toward 2000 s/m with drying (see Figure 7).
	double Saturation_SoilStorage_frac = 0.1;
	Resistance_Surface_Soil_s_p_m = exp(8.206 - 4.255 * (Saturation_SoilStorage_frac * (Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3)) /	(Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3));
	//Saturation_SoilStorage_zero_frac defined as zero to find maximum likely Resistance_Surface_Soil_s_p_m
	double Saturation_SoilStorage_zero_frac = 0;
	//Resistance_Surface_Soil_Max_s_p_m defined for Saturation_SoilStorage_zero_frac condition
	double Resistance_Surface_Soil_Max_s_p_m = exp(8.206 - 4.255 * (Saturation_SoilStorage_zero_frac * (Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3)) /	(Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3));

	//Resistance_Aerodynamic_Impervious_s_p_m (s/m) flux resistance from impervious to canopy level equal to Resistance_Aerodynamic_Soil_s_p_m, as in Yang et al. (2013)
	//Note: Consider refactor to explore impervious and water has greater or lesser resistance than soil
	Resistance_Aerodynamic_Impervious_s_p_m = Resistance_Aerodynamic_Soil_s_p_m;
	Resistance_Aerodynamic_Water_s_p_m = Resistance_Aerodynamic_Soil_s_p_m;

	//Consider putting Resistance_Aerodynamic_base_s_p_m into HydroPlusConfig.xml to ensure a minimum resistance value
	//Resistance_Aerodynamic_base_s_p_m (s/m) is the minimum aerodynamic resistance to allow values comparable to Choudhury and Monteith (1988) 
	//Note: Algorithm is modified from Choudhury and Monteith (1988), who combined soil and tree resistances
	double Resistance_Aerodynamic_base_s_p_m = 0;

	//Resistance_Aerodynamic_Tree_s_p_m (s/m) combined with Resistance_Aerodynamic_base_s_p_m to ensure values never below base, keeping LE reasonable	
	Resistance_Aerodynamic_Tree_s_p_m = Resistance_Aerodynamic_Tree_s_p_m + Resistance_Aerodynamic_base_s_p_m;
	Resistance_Aerodynamic_SVeg_s_p_m = Resistance_Aerodynamic_SVeg_s_p_m + Resistance_Aerodynamic_base_s_p_m;
	Resistance_Aerodynamic_Impervious_s_p_m = Resistance_Aerodynamic_Impervious_s_p_m + Resistance_Aerodynamic_base_s_p_m;
	Resistance_Aerodynamic_Soil_s_p_m = Resistance_Aerodynamic_Soil_s_p_m + Resistance_Aerodynamic_base_s_p_m;
	Resistance_Aerodynamic_Water_s_p_m = Resistance_Aerodynamic_Water_s_p_m + Resistance_Aerodynamic_base_s_p_m;
	
	//Save wind speed and resistance values for output
	beC->by_key(MapPixel_ID, DataFolder_ID, "WindSpeed_folder_m_p_s") = WindSpeed_folder_m_p_s;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Tree_s_p_m") = Resistance_Aerodynamic_Tree_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_SVeg_s_p_m") = Resistance_Aerodynamic_SVeg_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Impervious_s_p_m") = Resistance_Aerodynamic_Impervious_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Soil_s_p_m") = Resistance_Aerodynamic_Soil_s_p_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Water_s_p_m") = Resistance_Aerodynamic_Water_s_p_m;
}
