#include "ShortVegInterceptionCalc.h"

//Note: Consider refactor if snow intercepting, storing, unloading and below canopy simulated differently for grass or other short vegetation
//Note: Consider refactor to call either ShortVegInterceptionCalc or TreeInterceptionCalc, sending the parameters of ShortVeg or Tree 

//Theory: The following Snow equations were developed by Dr. Yang Yang during their PhD work at SUNY ESF w. T. Endreny; see Yang et al. (2011)
//Note: Algorithms originally based on USACE (1998, 1956), updated by Endreny based on work of Lundquist et al. (2021) and Sexstone et al. (2018)
//Note: The PASATH, UnifiedHydro, and HydroPlus previously wrote Eqns using Tair and Tdew as Fahrenheit, and with coefficients errors explained below
//Note: SnowAblationOpenAreaCalc::SnowAblationOnGround function accounts for snow melt and snow sublimation, which combine to form ablation
//Note: Variables with snow depth are liquid and snow water equivalent depths, unless variable name contains SnowCrystal and is a frozen crystalline depth. 

//References:
//USACE, 1998. Engineering and Design: Runoff From Snowmelt. US Army Corps of Engineers Report Engineer Manual 1110-2-1406, Washington, D.C.
//Wang, J., Endreny, T. A., & Nowak, D. J. (2008). Mechanistic simulation of tree effects in an urban water balance model. Journal of the American Water Resources Association, 44(1), 75-85. doi:10.1111/j.1752-1688.2007.00139.x
//Yang, Y., Endreny, T. A., & Nowak, D. J. (2011). i-Tree Hydro: Snow Hydrology Update for the Urban Forest Hydrology Model. Journal of the American Water Resources Association (JAWRA), 47(6), 1211-1218. doi:10.1111/j.1752-1688.2011.00564.x

void ShortVegInterceptionCalc::ShortVegInterceptionManager(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//ThroughFall_Rain_SVegCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_SVegCanopy_m") = 0.0;
	//ThroughFall_Snow_SVegCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") = 0.0;
	//Interception_Rain_SVegCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_SVegCanopy_m") = 0.0;
	//Interception_Snow_SVegCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_SVegCanopy_m") = 0.0;
	//Melt_Snow_SVegCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_SVegCanopy_m") = 0.0;
	//Irrigation_SVegInterceptionStorage_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Irrigation_SVegInterceptionStorage_m") = 0;

	//If ShortVegCover_noTreeCanopy_frac > 0 then accumulate ThroughFall and Interception variables
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") > 0) {

		//Storage_SVegCanopy_Max_m (m) = LeafAreaStorage_SVeg_mm * Ratio_m_to_mm * LAI_SVeg_m2_p_m2 + BarkAreaStorage_SVeg_mm * Ratio_m_to_mm * (LAI_BAI_SVeg_m2_p_m2 - LAI_SVeg_m2_p_m2)
		//Note: Storage_SVegCanopy_Max_m is total depth in short veg, not per unit area of short veg; use input->InputXml not beC for access
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m") = input->InputXml["LeafAreaStorage_SVeg_mm"] * Ratio_m_to_mm * input->LAI_SVeg_m2_p_m2[timeStep] + input->InputXml["BarkAreaStorage_SVeg_mm"] * Ratio_m_to_mm * (input->LAI_BAI_SVeg_m2_p_m2[timeStep] - input->LAI_SVeg_m2_p_m2[timeStep]);

		//If Tair_C (C) > 0 C And Precip_mpdt (m/dt) > 0 then precipitation is falling as rain
		if (input->Tair_C[timeStep] > 0.0 && input->Precip_mpdt[timeStep] > 0) {
			//Call Interception_Rain_SVeg function to intercept rain
			Interception_Rain_SVeg(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		}
		//Else If Tair_C (C) <= 0 C And Precip_mpdt (m/dt) > 0 then precipitation is falling as snow
		else if (input->Tair_C[timeStep] <= 0.0 && input->Precip_mpdt[timeStep] > 0) {
			//Call Interception_Snow_SVeg function to intercept snow
			Interception_Snow_SVeg(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		}

		//If input->Precip_mpdt[timeStep] and Storage_Rain_SVegCanopy_m <= 0 then enter to consider irrigation scenario
		if (Inputs::isLessThanOrAlmostEqual(input->Precip_mpdt[timeStep],0) && Inputs::isLessThanOrAlmostEqual(beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m"),0)) {

			//Flag_Scenario_CoolAir_Irrigate_str = to_string(input->SimulationScenarios["Flag_Scenario_CoolAir_Irrigate"])
			string Flag_Scenario_CoolAir_Irrigate_str = to_string(input->SimulationScenarios["Flag_Scenario_CoolAir_Irrigate"]);
			//If input->isReferenceStationPixel(MapPixel_ID) is false and ...
			//... Model_Selection is SpatialTemperatureHydro and ...
			//... Flag_Scenario_CoolAir_Irrigate is for 2 or 5, ShortVeg
			if (!input->isReferenceStationPixel(MapPixel_ID) &&
				input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro" &&
				(Flag_Scenario_CoolAir_Irrigate_str.find('2') != string::npos ||
					Flag_Scenario_CoolAir_Irrigate_str.find('5') != string::npos)) {

				//LandCoverType = "ShortVegCover"
				string LandCoverType = "ShortVegCover";
				//Flag_IsAppropriateLandCover = input->IsAppropriateLandCoverForScenario(MapPixel_ID,Scenario_CoolAir_Irrigate_LandCover_Class)
				//Note: Determine whether to suppress irrigation based on land cover rule
				bool Flag_IsAppropriateLandCover = input->IsAppropriateLandCoverForScenario(MapPixel_ID, input->SimulationScenarios["Scenario_CoolAir_Irrigate_LandCover_Class"]);

				//Flag_IsAppropriateDay_or_Night = input->IsAppropriateDay_or_NightForScenario(timeStep, ElevationAngle_Solar_rad,Scenario_CoolAir_Irrigate_Day_or_Night)
				//Note: Determine whether to allow irrigation for the Day or Night rule; 0=day and night; 1=day; 2=night
				bool Flag_IsAppropriateDay_or_Night = input->IsAppropriateDay_or_NightForScenario(timeStep, beC->by_key(MapPixel_ID, DataFolder_ID, "ElevationAngle_Solar_rad"), input->SimulationScenarios["Scenario_CoolAir_Irrigate_Day_or_Night"]);

				//Flag_IsAppropriateTimesPerDay = input->IsAppropriateTimesPerDay(MapPixel_ID, SimulationDate_YYYYMMDD, LandCoverType, Scenario_CoolAir_Irrigate_MaxPerDay)
				//Note: Determine count of irrigation for each distinct MapPixel, Date, and land cover type
				bool Flag_IsAppropriateTimesPerDay = input->IsAppropriateTimesPerDay(MapPixel_ID, input->SimulationDate_YYYYMMDD[timeStep], LandCoverType, input->SimulationScenarios["Scenario_CoolAir_Irrigate_MaxPerDay"]);

				//If Flag_IsAppropriateLandCover & Flag_IsAppropriateDay_or_Night & Flag_IsAppropriateTimesPerDay = true then appropriate conditions scenario
				if (Flag_IsAppropriateLandCover && Flag_IsAppropriateDay_or_Night && Flag_IsAppropriateTimesPerDay) {

					//Tair_K intialized
					double Tair_K;
					// If !beC->has_var("Tair_K"); checks if beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K") initialized after timeStep 1
					if (!beC->has_var("Tair_K")) {
						//Tair_K = input->Tair_C for reference station and coverted from C to K
						Tair_K = input->Tair_C[timeStep] + 273.15;
					}
					//Else Tair_K initialized
					else {
						//Tair_K = beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K") defined for folder
						Tair_K = beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K");
					}

					//If Scenario_CoolAir_Irrigate_Tair_K < Tair_K then enter for irrigation scenario
					if (Inputs::isLessThanOrAlmostEqual(input->SimulationScenarios["Scenario_CoolAir_Irrigate_Tair_K"], Tair_K)) {

						//DepressionStorageDepth_fraction (frac) is defined by default as 1.0 of total storage
						double DepressionStorageDepth_fraction = 1.0;
						//Irrigated_frac (frac) initialized to 1.0
						double Irrigated_frac = 1.0;
						//If input->map_Irrigation_Area_frac.size() > 0 then redefine 
						if (input->map_Irrigation_Area_frac.size() > 0) {
							//Irrigated_frac = input->map_Irrigation_Area_frac[MapPixel_ID]
							Irrigated_frac = input->map_Irrigation_Area_frac[MapPixel_ID];
						}

						//Irrigation_SVegInterceptionStorage_m is DepressionStorageDepth_fraction * Storage_SVegCanopy_Max_m * Irrigated_frac
						//Note: Scenario activated to provide cooling
						beC->by_key(MapPixel_ID, DataFolder_ID, "Irrigation_SVegInterceptionStorage_m") = DepressionStorageDepth_fraction * beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m") * Irrigated_frac;
						//Storage_Rain_SVegCanopy_m is set to Irrigation_SVegInterceptionStorage_m
						beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Irrigation_SVegInterceptionStorage_m");
						//input->IncreaseTimesPerDay(MapPixel_ID, input->SimulationDate_YYYYMMDD[timeStep], LandCoverType); increases irrigation count
						input->IncreaseTimesPerDay(MapPixel_ID, input->SimulationDate_YYYYMMDD[timeStep], LandCoverType);
					}
				}
			}
		} //End Scenario

		//Note: Consider refactor to have this called only if wind speed above certain threshold
		//Call Unload_Snow_SVeg function to unload snow by wind
		Unload_Snow_SVeg(input, beC, MapPixel_ID, DataFolder_ID, timeStep);

		//If Storage_Snow_SVegCanopy_m (m) is greater than zero, then allow for snowmelt by radiation 
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") > 0.0) {
			//Call SnowMelt_by_Radiation_SVeg to melt snow by net radiation 
			SnowMelt_by_Radiation_SVeg(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		}
	}
	//ThroughFall_RainSnow_SVegCanopy_m (m) is the combined throughfall of rain water throughfall and snow-water equivalent throughfall
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_RainSnow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_SVegCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m");
	//Interception_RainSnow_SVegCanopy_m (m) includes combined interception of rain the water interception and snow-water equivalent interception
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_RainSnow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_SVegCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_SVegCanopy_m");
}

void ShortVegInterceptionCalc::Unload_Snow_SVeg(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//Flux_Snow_SVegCanopy_to_Ground_m (m) from Eq 4 of Yang et al. (2011) gives snow unloaded from canopy and sent to ground
	//Note: Eq 4 gives coefficient 0.0231 * WindSpd_mps at mean canopy height (m/s) * Storage_Snow_SVegCanopy_m (m)
	double Flux_Snow_SVegCanopy_to_Ground_m = 0.0231 * input->WindSpd_mps[timeStep] * beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m");
	
	//If Flux_Snow_SVegCanopy_to_Ground_m (m) greater than Storage_Snow_SVegCanopy_m (m) then set flux to storage
	if (Flux_Snow_SVegCanopy_to_Ground_m > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m")) {
		//Flux_Snow_SVegCanopy_to_Ground_m (m) set equal to Storage_Snow_SVegCanopy_m (m), as flux cannot exceed available storage 
		Flux_Snow_SVegCanopy_to_Ground_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m");
	}
	//If Flux_Snow_SVegCanopy_to_Ground_m < 0 then set to zero
	if (Flux_Snow_SVegCanopy_to_Ground_m < 0) { Flux_Snow_SVegCanopy_to_Ground_m = 0; }

	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") - Flux_Snow_SVegCanopy_to_Ground_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") + Flux_Snow_SVegCanopy_to_Ground_m;

}

void ShortVegInterceptionCalc::SnowMelt_by_Radiation_SVeg(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//Obtain weather data from input functions
	double Tair_C = input->Tair_C[timeStep];
	double Tdew_C = input->Tdew_C[timeStep];
	double Precipitation_m = input->Precip_mpdt[timeStep];
	double Radiation_Shortwave_Direct_Wpm2 = input->Radiation_Shortwave_Direct_Wpm2[timeStep];
	double WindSpd_mps = input->WindSpd_mps[timeStep];
	double AtmPres_kPa = input->AtmPres_kPa[timeStep];
	//Standard pressure at sea level, or standard atmosphere, set to 101.325 kiloPascal; https://en.wikipedia.org/wiki/Atmospheric_pressure
	double AtmPres_SeaLevel_Standard_kPa = 101.325;
	//Precipitation_Rain_m defined and initialized to zero
	double Precipitation_Rain_m = 0.0;

	//Ratio_AtmPres_to_StdPres is the ratio of atmospheric pressure at the location to the atmospheric pressure at the sea level
	double Ratio_AtmPres_to_StdPres = AtmPres_kPa / AtmPres_SeaLevel_Standard_kPa;
	//Albedo_Snow_frac (fraction) defined as 0.85
	//Note: Consider refactor for time varying snow albedo of NREL (1995) p 78; track snow depth and days since snow for each weather station
	double Albedo_Snow_frac = 0.85;
	//SnowMelt_ViaRainfall_m_p_hr (m) is the snow melt by rain, initializing to 0.0
	double SnowMelt_ViaRainfall_m_p_hr = 0.0;

	//SnowMelt_RadiationShortWave_m_p_hr (m/hr) snowmelt by shortwave radiation defined by Eq 6 1st RHS in Yang et al. (2011), derived from Eq 5-2 (USACE, 1998) ...
	//Note: ... Coefficient_1 = 0.00308 mm/day kJ/m2/day * (1/1000)*(1/24) = 0.000000128 m/hr kJ/m2/day
	//Note: ... Coefficient_1 = 0.000000128 m/hr kJ/m2/day * (1/1000)*(24/1)*(3600/1) = 0.00001109 m/hr W/m2
	//Note: Consider refactor to update Albedo_Snow_frac
	double Coeff_01_Mn_m_p_hr_W_p_m2 = 0.00001109;
	double SnowMelt_RadiationShortWave_m_p_hr = Coeff_01_Mn_m_p_hr_W_p_m2 * Radiation_Shortwave_Direct_Wpm2 * (1 - Albedo_Snow_frac);

	//SnowMelt_RadiationLongWave_m_p_hr (m/hr) = the snowmelt by longwave radiation initialized to zero
	double SnowMelt_RadiationLongWave_m_p_hr = 0.0;
	//If Precipitation_m greater than zero, then presume a cloudy day and use appropriate equation for cloudy day snowmelt
	if (Precipitation_m > 0) {
		//SnowMelt_RadiationLongWave_m_p_hr (m) is snowmelt by long wave radiation on cloudy day, replacing 2nd RHS terms of Eq 6 in Yang et al. (2011) ...
		//Note: ... as described in paragraph below Eq 10 (Yang et al., 2011) with 0.000054*Ta, which has rounding error
		//Note: 2nd RHS of Eq 6 in Yang et al, (2011), 0.000054*Ta, was derived from Eq 5-9 (USACE, 1998) for cloudy sky, however ...
		//Note: ... Eq 5-9 (USACE, 1998) Coefficient_1 = 0.029 in/day/F with its derivation illustrated in top line of Figure 5-1
		//Note: ... Coefficient_1 = 0.029 in/day (F)  * (1/12)*(0.3048/1)*(1/24) = 0.000030692 m/hr (F)
		//Note: ... Coefficient_1 = 0.000030692 m/hr (F) * 9/5 degree_F/degree_C = 0.00005525 m/hr (C)
		double Coeff_01_Ml_cloudy_m_p_hr_C = 0.00005525;
		SnowMelt_RadiationLongWave_m_p_hr = Coeff_01_Ml_cloudy_m_p_hr_C * Tair_C;
	}
	//Else If Precipitation zero then presume clear skies and use appropriate equation for clear day snowmelt
	else {
		//SnowMelt_RadiationLongWave_m_p_hr (m/hr) is snowmelt by long wave radiation on clear day, 2nd RHS terms of Eq 6 in Yang et al. (2011) 
		//Note: 2nd RHS of Eq 6 in Yang et al, (2011) was derived from Eq 5-8 (USACE, 1998) for clear sky, however it contains error due to ...
		//Note: ... error in Eq 5-8 (USACE, 1998) Coefficient_2 = 0.84 in/day (F), taken from lower line of Figure 5-1, and should be 0.804 or 0.8
		//Note: ... Hence Yang et al. (2011) should have used Coefficient_2 = 0.0008509 and not Coefficient_2 = 0.00089
		//Note: Figure 5-1 (USACE, 1998) shows air temperature needs to be 70 deg F for melt due to long wave radiation in clear sky conditions
		//Note: Given for clear skies Fig 5-1 shows 0.22 inches/day melts when temperature is 70 F, Eq 5-8 should have 2nd coefficient of 0.8, not 0.84
		//Note: Given so little snow melts under clear skies, we should create overcast conditions based on vapor pressure or cloud data, not just Precipitation_m > 0
		//Note: ... Coefficient_2 = 0.804 in/day (F) * (1/12)*(0.3048/1)*(1/24) = 0.0008509 m/hr (F) = 0.0008509 m/hr (C)  
		//Note: ... Coefficient_1 = 0.0212 in/day (F)  * (1/12)*(0.3048/1)*(1/24) = 0.000022437 m/hr (F)
		//Note: ... Coefficient_1 = 0.000022437 m/hr (F) * 9/5 degree_F/degree_C = 0.000040386 m/hr (C)
		double Coeff_01_Ml_clear_m_p_hr_C = 0.000040386;
		double Coeff_02_Ml_clear_m_p_hr_C = 0.0008509;
		SnowMelt_RadiationLongWave_m_p_hr = Coeff_01_Ml_clear_m_p_hr_C * Tair_C - Coeff_02_Ml_clear_m_p_hr_C;   // Melt by longwave radiation on clear day
	}

	//If melt by long wave radiation < 0, then set it to 0; 
	//Note: Correction needed when clear sky and temperature below 70F
	if (SnowMelt_RadiationLongWave_m_p_hr < 0) {
		SnowMelt_RadiationLongWave_m_p_hr = 0;
	}

	//SnowMelt_Radiation_m_p_hr (m) from Eq 5.6 of Yang et al. (2011) is the total melt by radiation, the shortwave melt plus the long wave melt
	//Note: Eq 5.6 of Yang et al. (2011) derived from Eq 5-1 to 5.9 of USACE (1998), accounting for back radiation as function of air temperature, rain or snow, and cloud cover effects
	double SnowMelt_Radiation_m_p_hr = SnowMelt_RadiationShortWave_m_p_hr + SnowMelt_RadiationLongWave_m_p_hr;

	//If SnowMelt_Radiation_m_p_hr (m) greater than Storage_Snow_SVegCanopy_m (m), then set to Storage_Snow_SVegCanopy_m (m) as available snow to melt
	if (SnowMelt_Radiation_m_p_hr > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m")) {
		SnowMelt_Radiation_m_p_hr = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m");
	}
	//If SnowMelt_Radiation_m_p_hr < 0 then set to zero
	if (SnowMelt_Radiation_m_p_hr < 0) { SnowMelt_Radiation_m_p_hr = 0; }

	//SnowMelt_Radiation_m = SnowMelt_Radiation_m_p_hr * input->SimulationTimeStep_Duration_sec[timeStep] * Ratio_Hour_to_Second
	//Note: Conversion from m/hr to m based on SimulationTimeStep_Duration_sec
	double SnowMelt_Radiation_m = SnowMelt_Radiation_m_p_hr * input->SimulationTimeStep_Duration_sec[timeStep] * Ratio_Hour_to_Second;

	//Melt_Snow_SVegCanopy_m (m) is beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_SVegCanopy_m") + SnowMelt_Radiation_m; increased by SnowMelt_Radiation_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_SVegCanopy_m") + SnowMelt_Radiation_m;
	//Storage_Snow_SVegCanopy_m (m) is reduced by SnowMelt_Radiation_m_p_hr (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") - SnowMelt_Radiation_m;
	//Storage_Rain_SVegCanopy_m (m) is increased by SnowMelt_Radiation_m_p_hr (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") + SnowMelt_Radiation_m;

	//If Storage_Rain_SVegCanopy_m (m) is greater than Storage_SVegCanopy_Max_m (m) then adjust storage and throughfall to fit maximum available
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m")) {
		//ThroughFall_Snow_SVegCanopy_m (m) increased by Storage_Snow_SVegCanopy_m (m) due to overflowing water pushing snow off of canopy
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m");
		//Storage_Snow_SVegCanopy_m (m) set to zero
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") = 0.0;
		//ThroughFall_Rain_SVegCanopy_m (m) is increased by Storage_Rain_SVegCanopy_m (m) and decreased by Storage_SVegCanopy_Max_m (m), that which remains in storage
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_SVegCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m");
		//Storage_Rain_SVegCanopy_m (m) is set to Storage_SVegCanopy_Max_m (m) 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m");
	}

	//Storage_Snow_SVegCanopy_m = max(0.0, beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m"))
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") = max(0.0, beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m"));
	//Storage_Rain_SVegCanopy_m = max(0.0, beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m"))
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") = max(0.0, beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m"));
}

void ShortVegInterceptionCalc::Interception_Rain_SVeg(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//Obtain weather data from input functions
	double Tair_C = input->Tair_C[timeStep];
	//LAI_BAI_SVeg_m2m2 (m2/m2) is a combined LAI and BAI per unit of land
	double LAI_BAI_SVeg_m2m2 = input->LAI_BAI_SVeg_m2_p_m2[timeStep];
	//ThroughFall_Canopy_frac (the short veg throughfall coefficient) =  e^(-0.3 * LAI)
	//Note: Modified form of Eq 2a & Eq 2b from Wang, Endreny, Nowak (2008, JAWRA) based on Dijk and Bruijnzeel (2001) for free throughfall
	//Note: Wang et al. (2008) after Eq 2b assign extinction coefficient k = 0.3 for short veg
	double k_ExtinctionCoefficient = 0.3;
	double ThroughFall_Canopy_frac  =  exp(-k_ExtinctionCoefficient * LAI_BAI_SVeg_m2m2);

	//ThroughFall_Rain_SVegCanopy_m (m) equals ThroughFall_Canopy_frac * Precip_mpdt, from Valente et al., (1997) sparse Rutter model
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_SVegCanopy_m") = ThroughFall_Canopy_frac * input->Precip_mpdt[timeStep];
	//Interception_Rain_SVegCanopy_m (m) equals (1 - ThroughFall_Canopy_frac) * Precip_mpdt, from Valente et al., (1997) sparse Rutter model
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_SVegCanopy_m") = (1 - ThroughFall_Canopy_frac) * input->Precip_mpdt[timeStep];

	//SnowMelt_Rain_m_p_hr (m) is the snow melt by rain, initializing to 0.0
	double SnowMelt_Rain_m_p_hr = 0.0;
	//If Tair_C air temperature is above 0C or 32F, then precipitation treated as liquid rain
	if (Tair_C > 0.0) {
		//SnowMelt_ViaRainfall_m_p_hr = Coeff_01_Mp_m_p_h_C * Interception_Rain_SVegCanopy_m * Tair_C
		//Note: Snowmelt by precipitation, from Eq 10 of Yang et al. (2011), but used erroneous coefficient from USACE (1998)
		//Note: Eq 5-17 of USACE (1998) gives proper Coefficient_1 0.0125 mm/day for rain in mm/day C, which becomes 0.125 m/hr for rain in m/hr C
		//Note: Eq 10 of Yang et al. (2011) was derived from Eq 2-9 (USACE, 1998) gave Qp = Cp_r * rho_w * Precipitation_m * (Tr - Ts)/1000, but ...
		//Note: ... RHS division by 1000 appears erroneous, based on Eq 5-17 coefficient which appears correct
		//Note: Derivation of coefficient with Eq 2-9 and Eq 2-2 (USACE, 1998), showing Eq 2-9 RHS division by 1000 should be removed ...
		//Note: USACE (1998) Eq 5-9 should be Qp = Cp_r * rho_w * P * (Tr - Ts), should not divide by 1000, and then Qp = 840,000 kJ/m2/day when ...
		//Note: ... Cp_r = 4.2 kJ/kg/C specific heat of rain, rho_w = 1000 kg/m3 density of water, P = 20 mm/day of rainfall, ...
		//Note: ... Tr = 10 C temperature of rain, Ts = 0 C, temperature of snow, and if needed, Cp_s = 2.08 kJ/kg/C specific heat of snow
		//Note: ... and Qp = =4.2 * 1000 * 20 * (10 - 0) = 840000 kJ/m2/day when Qp = Cp_r * rho_w * P * (Tr - Ts)
		//Note: Result of Eq 5-9, Qp, goes into Eq 2-2 (USACE, 1988), to determine snowmelt by rainfall, M = Qp / (B * F * rho_w), where
		//Note: ... B = 0.97 thermal quality of snow, F = 334.9 kJ/kg latent heat of fusion, rho_w = 1000 kg/m3 density of water, ...
		//Note: ... and M = 840000 / (0.97 * 334.9 * 1000) = 2.585784955 mm/day = 0.000107741 mm/hr
		double Coeff_01_Mp_m_p_h_C = 0.0125;
		SnowMelt_Rain_m_p_hr = Coeff_01_Mp_m_p_h_C * beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_SVegCanopy_m") * Tair_C;
	}

	//SnowMelt_Rain_m = SnowMelt_Rain_m_p_hr * input->SimulationTimeStep_Duration_sec[timeStep] * Ratio_Hour_to_Second
	//Note: Conversion from m/hr to m based on SimulationTimeStep_Duration_sec
	double SnowMelt_Rain_m = SnowMelt_Rain_m_p_hr * input->SimulationTimeStep_Duration_sec[timeStep] * Ratio_Hour_to_Second;
	
	//If SnowMelt_Rain_m (m) greater than Storage_Snow_SVegCanopy_m (m), set equal as melted snow cannot be larger than the stored snow
	if (SnowMelt_Rain_m > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m")) {
		//SnowMelt_Rain_m (m) set equal to Storage_Snow_SVegCanopy_m (m)
		SnowMelt_Rain_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m");
	}

	//Storage_Snow_SVegCanopy_m (m) is reduced by SnowMelt_Rain_m (m) to reflect the melted snow
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") - SnowMelt_Rain_m;

	//Storage_Rain_SVegCanopy_m (m) is increased by Interception_Rain_SVegCanopy_m (m) and SnowMelt_Rain_m (m)
	double Storage_Rain_SVegCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_SVegCanopy_m") + 
		SnowMelt_Rain_m;

	//If Storage_Rain_SVegCanopy_m (m) greater than Storage_SVegCanopy_Max_m (m) then adjust interception and throughfall 
	if (Storage_Rain_SVegCanopy_m > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m")) {
		//Interception_Rain_SVegCanopy_m (m) is set to Storage_SVegCanopy_Max_m (m) minus Storage_Rain_SVegCanopy_m (m) and SnowMelt_Rain_m (m)
		//Note: Subtraction reduces interception by the amount already in storage
		beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") - SnowMelt_Rain_m;
		//ThroughFall_Rain_SVegCanopy_m (m) is increased by Storage_Rain_SVegCanopy_m (m) and reduced by Storage_SVegCanopy_Max_m (m)
		//Note: Reduction by Storage_SVegCanopy_Max_m accounts for only removing what is not held in maximum storage
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_SVegCanopy_m") + Storage_Rain_SVegCanopy_m - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m");
		//ThroughFall_Snow_SVegCanopy_m (m) is increased by Storage_Snow_SVegCanopy_m (m) that falls as liquid equivalent
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m");
		//Storage_Snow_SVegCanopy_m (m) set to zero after throughfall clears load
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") = 0.0;
		//Storage_Rain_SVegCanopy_m (m) set to Storage_SVegCanopy_Max_m (m)
		Storage_Rain_SVegCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m");
	}

	//Storage_Rain_SVegCanopy_m = Storage_Rain_SVegCanopy_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") = Storage_Rain_SVegCanopy_m;
	//Melt_Snow_SVegCanopy_m (m) is snow melted in short veg canopy due to rain, later added to that melted by radiation
	beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_SVegCanopy_m") = SnowMelt_Rain_m;
	//Control for negative numbers
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_SVegCanopy_m") = max(0.0, beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_SVegCanopy_m"));
}

//Interception_Snow_SVeg function will handle snow interception in the SVeg canopy
//Note: Consider refactor to improve equations for SnowDensity_kg_p_m3 and Storage_SnowCrystal_Leaf_Max_m
void ShortVegInterceptionCalc::Interception_Snow_SVeg(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//LAI_BAI_SVeg_m2m2 (m2/m2) is a combined LAI and BAI per unit of land
	double LAI_BAI_SVeg_m2m2 = input->LAI_BAI_SVeg_m2_p_m2[timeStep];
	//Tair_C (deg C) is Tair_C 
	double Tair_C = input->Tair_C[timeStep]; 
	//SnowDensity_kg_p_m3 (kg/m3) equals 67.92 + 51.25 * exp(Tair_C/2.59), source of equation not provided by Yang
	//Note: SnowDensity simplified for model, and depends on more than Tair_C: https://www.vcalc.com/wiki/KurtHeckman/Snow+Water+Density+%28SWD%29
	double SnowDensity_kg_p_m3 = 67.92 + 51.25 * exp(Tair_C/2.59);

	//Note: Variables with snow depth are liquid and snow water equivalent depths, unless variable name contains SnowCrystal and is a frozen crystalline depth. 
	//Note: Coefficient 6, 0.27, and 46.0 need documentation. Evidently 6 is species specific and can be changed to fit data
	double Storage_SnowCrystal_Leaf_Max_m = 6.0 * (0.27 + 46.0/SnowDensity_kg_p_m3)/SnowDensity_kg_p_m3;
	double Storage_SnowCrystal_SVegCanopy_Max_m = Storage_SnowCrystal_Leaf_Max_m * LAI_BAI_SVeg_m2m2;
	//Note: Consider refactor to explore logic for division by 1000; was Storage_SnowCrystal_SVegCanopy_Max_m considered a volume?
	double Storage_Snow_SVegCanopy_Max_m = Storage_SnowCrystal_SVegCanopy_Max_m * (SnowDensity_kg_p_m3 / 1000);

	//Eq 2a & Eq 2b from Wang, Endreny, Nowak (2008, JAWRA) based on Dijk and Bruijnzeel (2001) for free throughfall
	//Wang et al. (2008) after Eq 2b assign extinction coefficient k = 0.3 to short vegetation
	double ThroughFall_Canopy_frac = (1-(1-exp(-0.3 * LAI_BAI_SVeg_m2m2)));
	double Precipitation_Snow_m = input->Precip_mpdt[timeStep];
	
	//Note: Theory from Valente et al., (1997)
	//ThroughFall_Snow_SVegCanopy_m (m) equals ThroughFall_Canopy_frac * Precipitation_Snow_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") = ThroughFall_Canopy_frac * Precipitation_Snow_m;
	//Interception_Snow_SVegCanopy_m (m) equals (1 - ThroughFall_Canopy_frac) * Precipitation_Snow_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_SVegCanopy_m") = (1-ThroughFall_Canopy_frac) * Precipitation_Snow_m;

	//Storage_Snow_SVegCanopy_m (m) is defined as sum of Storage_Snow_SVegCanopy_m (m) and Interception_Snow_SVegCanopy_m (m)
	double Storage_Snow_SVegCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_SVegCanopy_m");

	//If Storage_Snow_SVegCanopy_m (m) greater than Storage_Snow_SVegCanopy_Max_m (m) then reduce interception, increase throughfall
	if (Storage_Snow_SVegCanopy_m > Storage_Snow_SVegCanopy_Max_m)	{
		//Interception_Snow_SVegCanopy_m is Storage_Snow_SVegCanopy_Max_m - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m")
		beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_SVegCanopy_m") = Storage_Snow_SVegCanopy_Max_m - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m");
		//ThroughFall_Snow_SVegCanopy_m is beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") + (Storage_Snow_SVegCanopy_m - Storage_Snow_SVegCanopy_Max_m)
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_SVegCanopy_m") + (Storage_Snow_SVegCanopy_m - Storage_Snow_SVegCanopy_Max_m);
		//Storage_Snow_SVegCanopy_m = Storage_Snow_SVegCanopy_Max_m; set to maximum
		Storage_Snow_SVegCanopy_m = Storage_Snow_SVegCanopy_Max_m;
	} 

	//Storage_Snow_SVegCanopy_m (m) set to Storage_Snow_SVegCanopy_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") = Storage_Snow_SVegCanopy_m;

	//If Interception_Snow_SVegCanopy_m (m) less than zero, set to zero 
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_SVegCanopy_m") = max(0.0, beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_SVegCanopy_m"));
}
