#include "TreeInterceptionCalc.h"


//Note: Consider refactor: Add stemflow algorithm to a modified version of the sparse Rutter model [Valente et al., 1997]. The stemflow will allow for additional infiltration in soils around trees with canopy over impervious cover, and likely require updates to the vadose zone recharge algorithm, StorageDeficit_VadoseZone.cpp. The stemflow can also send excess canopy interception to the tree trunk, i.e., main stem, by way of the tree branches, and thereby allow for a larger atmospheric area to fill of storage in the non-vascular vegetation (NVV) that is vertically stacked along the trunk. While the i-Tree Cool Air model already contains the sparse Rutter model interception and throughfall routines, they do not represent stemflow mechanisms. The new stemflow algorithm will allow funneling of precipitation from the larger tree canopy area to the smaller tree trunk area, and fill storage available in the NVV that is vertically stacked on the trunk. For the fraction of grid in tree canopy (c) the i-Tree Cool Air model computes a maximum canopy storage capacity without NVV, Sc_bare (m). This is based on maximum storage depth on leaves (Sc_L) and bark (Sc_B) of the tree, and the tree leaf area index (LAI) (m2/ m2), or area of upward facing leaf per unit area of ground, and the bark area index (BAI) (m2/ m2) or area of upward facing bark per unit area of ground. The calculation is Sc_bare = Sc_L(LAI) + Sc_L(BAI). The model also tracks the actual canopy storage, Cc (m), filling with precipitation and emptying with evaporation as described by Valente et al. [1997], and at any time step if Cc > Sc_bare, then the excess precipitation becomes canopy drainage, Dc = Cc – Sc. Currently Dc is sent directly to throughfall. The new stemflow algorithm will allow Dc to be partitioned to throughfall or stemflow based on the drainage partitioning coefficient, pd, a fraction between 0 and 1. The amount going to trunk storage, Ct,c (m) is defined as Dc(pd), and if Ct,c exceeds trunk storage capacity, St,c, the excess becomes trunk drainage to the ground. The amount of Dc going to throughfall is Dc(1-pd). To parameterize this drainage partitioning coefficient, we propose pd = k (BAI) / LAI, where k represents the efficiency of water passing along the leaf petiole to the stem. The pd function represents the funneling effect of greater canopy area leading to smaller bark area, where common maximum value of urban BAI is 1.7 and urban LAI is 5.0 [Nowak et al., 2020], and if k = 0.5 then pd = 0.17. This was an idea developed in Endreny and Zhou's 2023 DOE FAIR proposal. 


//Note: Consider refactor: Non-vascular vegetation (NVV) fractional cover and storage will be added to the i-Tree Cool Air model as a new Plant Functional Type (PFTs), using algorithms from the LiBry model. Unfunded collaborator Porada is LiBry developer, and will assist with code transfer and translation (see letter). The algorithms will nearly replicate existing code in i-Tree Cool Air for representing leaf and bark fractional area and storage. The different PFTs cover unique fractions of the total area within the land unit and grid pixel. Community weighted mean values of the NVV fractional coverage on vegetation and bare soil, and that NVV water storage capacity, come from the survival strategies simulated by the LiBry non-vascular dynamic vegetation model [Porada et al., 2018]. NVV storage capacity values in the canopy and on the ground will generally range from 1 to 16 mm [Porada et al., 2018], nearly 10 times greater NVV capacity of depression storage on the ground, or bark and leaf storage in the canopy [Nowak et al., 2020]. The final step of this objective is to tune the i-Tree Cool Air model to work with the NVV updates. The simulations will compare model stemflow and NVV storage to observed values in the supplemental information of Porada et al. [2018], with database guidance provided by unfunded collaborator and NVV field expert Van Stan (see letter). The LiBry estimates were tuned to an observed values from 73 field sites of NVV biomass, storage capacity, interception, and evaporation. The i-Tree Cool Air model will use the LiBry model weather data forcing data to replicate the scenario at each of the seven biomes, and the final estimated fractional cover for NVV, bark, leaves, and ground. Our target in this model tuning exercise is for the i-Tree Cool Air model predictions to be within 10% of the LiBry estimated interception and evaporation rates. References: Nowak, D. J., R. Coville, T. A. Endreny, R. Abdi, and J. T. Van Stan II (2020), Valuing Urban Tree Impacts on Precipitation Partitioning, in Precipitation Partitioning by Vegetation - A Global Synthesis, edited by V. S. I. J.T., pp. 253-267, Springer Nature, Swam, Switzerland. Porada, P., J. T. Van Stan, and A. Kleidon (2018), Significant contribution of non-vascular vegetation to global rainfall interception, Nature Geoscience, 11(8), 563-567, doi: 10.1038/s41561-018-0176-7.


//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:
//Lundquist, J. D., Dickerson-Lange, S., Gutmann, E., Jonas, T., Lumbrazo, C., & Reynolds, D. (2021). Snow interception modelling: Isolated observations have led to many land surface models lacking appropriate temperature sensitivities. Hydrological Processes, 35(7), e14274. doi:https://doi.org/10.1002/hyp.14274
//USACE, 1998. Engineering and Design: Runoff From Snowmelt. US Army Corps of Engineers Report Engineer Manual 1110-2-1406, Washington, D.C.
//Valente, F., J. S. David, and J. H. C. Gash (1997), Modeling interception loss for two sparse eucalypt and pine forests in central Portugal using reformulated Rutter and Gash analytical models, Journal of Hydrology, 190, 141-162.
//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 TreeInterceptionCalc::TreeInterceptionManager(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	// initializing these values at the start of each timestep
	//ThroughFall_Rain_TreeCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_TreeCanopy_m") = 0.0;
	//ThroughFall_Snow_TreeCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") = 0.0;
	//Interception_Rain_TreeCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_TreeCanopy_m") = 0.0;
	//Interception_Snow_TreeCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_TreeCanopy_m") = 0.0;
	//Melt_Snow_TreeCanopy_m (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_TreeCanopy_m") = 0;
	//Irrigation_TreeInterceptionStorage_m  (m) reset to zero at each entry to function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Irrigation_TreeInterceptionStorage_m") = 0;

	//If TreeCover_frac > 0 then accumulate ThroughFall and Interception variables
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "TreeCover_frac") > 0) {
		//Storage_TreeCanopy_Max_m (m) = LeafAreaStorage_Tree_mm * Ratio_m_to_mm * LAI_Tree_m2_p_m2 + BarkAreaStorage_Tree_mm * Ratio_m_to_mm * (LAI_BAI_Tree_m2_p_m2 - LAI_Tree_m2_p_m2)
		//Note: Storage_TreeCanopy_Max_m is total depth in tree, not per unit area of tree; use input->InputXml not beC for access
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m") = input->InputXml["LeafAreaStorage_Tree_mm"] * Ratio_m_to_mm * input->LAI_Tree_m2_p_m2[timeStep] + input->InputXml["BarkAreaStorage_Tree_mm"] * Ratio_m_to_mm * (input->LAI_BAI_Tree_m2_p_m2[timeStep] - input->LAI_Tree_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_Tree function to intercept rain
			Interception_Rain_Tree(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_Tree function to intercept snow
			Interception_Snow_Tree(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		}

		//If input->Precip_mpdt[timeStep] and Storage_Rain_TreeCanopy_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_TreeCanopy_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 1 or 5, Tree
			if (!input->isReferenceStationPixel(MapPixel_ID) &&
				input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro" &&
				(Flag_Scenario_CoolAir_Irrigate_str.find('1') != string::npos ||
					Flag_Scenario_CoolAir_Irrigate_str.find('5') != string::npos)) {

				//LandCoverType = "TreeCover"
				string LandCoverType = "TreeCover";

				//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 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_TreeInterceptionStorage_m is DepressionStorageDepth_fraction * Storage_TreeCanopy_Max_m * Irrigated_frac
						//Note: Scenario activated to provide cooling
						beC->by_key(MapPixel_ID, DataFolder_ID, "Irrigation_TreeInterceptionStorage_m") = DepressionStorageDepth_fraction * beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m") * Irrigated_frac;
						//Storage_Rain_TreeCanopy_m is set to Storage_Rain_TreeCanopy_m
						beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Irrigation_TreeInterceptionStorage_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_Tree function to unload snow by wind
		Unload_Snow_Tree(input, beC, MapPixel_ID, DataFolder_ID, timeStep);

		//If Storage_Snow_TreeCanopy_m (m) is greater than zero, then allow for snowmelt by radiation 
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") > 0.0) {
			//Call SnowMelt_by_Radiation_Tree to melt snow by net radiation 
			SnowMelt_by_Radiation_Tree(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		}
	}

	//ThroughFall_RainSnow_TreeCanopy_m (m) includes combined throughfall of rain the water throughfall and snow-water equivalent throughfall
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_RainSnow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m");
	//Interception_RainSnow_TreeCanopy_m (m) includes combined interception of rain the water interception and snow-water equivalent interception
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_RainSnow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_TreeCanopy_m");
}

void TreeInterceptionCalc::Unload_Snow_Tree(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//Flux_Snow_TreeCanopy_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_TreeCanopy_m (m)
	double Flux_Snow_TreeCanopy_to_Ground_m = 0.0231 * input->WindSpd_mps[timeStep] * beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m");

	//If Flux_Snow_TreeCanopy_to_Ground_m (m) greater than Storage_Snow_TreeCanopy_m (m) then set flux to storage
	if (Flux_Snow_TreeCanopy_to_Ground_m > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m")) {
		//Flux_Snow_TreeCanopy_to_Ground_m (m) set equal to Storage_Snow_TreeCanopy_m (m), as flux cannot exceed available storage
		Flux_Snow_TreeCanopy_to_Ground_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m");
	}
	//If Flux_Snow_TreeCanopy_to_Ground_m < 0 then set to zero
	if (Flux_Snow_TreeCanopy_to_Ground_m < 0) { Flux_Snow_TreeCanopy_to_Ground_m = 0; }

	// Update the TreeStoredSWE.  TreeStoredSWE = (TreeStoredSWE - Flux_Snow_TreeCanopy_to_Ground_m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") - Flux_Snow_TreeCanopy_to_Ground_m;
	// Update the TreeThroughSWE.  TreeThroughSWE = (TreeThroughSWE + Flux_Snow_TreeCanopy_to_Ground_m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") + Flux_Snow_TreeCanopy_to_Ground_m;

}

void TreeInterceptionCalc::SnowMelt_by_Radiation_Tree(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
		//Note: Melt by longwave radiation on clear day
		SnowMelt_RadiationLongWave_m_p_hr = Coeff_01_Ml_clear_m_p_hr_C * Tair_C - Coeff_02_Ml_clear_m_p_hr_C;   
	}

	//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_TreeCanopy_m (m), then set to Storage_Snow_TreeCanopy_m (m) as available snow to melt
	if (SnowMelt_Radiation_m_p_hr > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m"))	{
		SnowMelt_Radiation_m_p_hr = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_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_TreeCanopy_m (m) is beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_TreeCanopy_m") + SnowMelt_Radiation_m; increased by SnowMelt_Radiation_m_p_hr (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_TreeCanopy_m") + SnowMelt_Radiation_m;
	//Storage_Snow_TreeCanopy_m (m) is reduced by SnowMelt_Radiation_m_p_hr (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") - SnowMelt_Radiation_m;
	//Storage_Rain_TreeCanopy_m (m) is increased by SnowMelt_Radiation_m_p_hr (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") + SnowMelt_Radiation_m;

	//If Storage_Rain_TreeCanopy_m (m) is greater than Storage_TreeCanopy_Max_m (m) then adjust storage and throughfall to fit maximum available
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m")) {
		//ThroughFall_Snow_TreeCanopy_m (m) increased by Storage_Snow_TreeCanopy_m (m) due to overflowing water pushing snow off of canopy
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m");
		//Storage_Snow_TreeCanopy_m (m) set to zero
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") = 0.0;
		//ThroughFall_Rain_TreeCanopy_m (m) is increased by Storage_Rain_TreeCanopy_m (m) and decreased by Storage_TreeCanopy_Max_m (m), that which remains in storage
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m");
		//Storage_Rain_TreeCanopy_m (m) is set to Storage_TreeCanopy_Max_m (m) 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m");
	}
	//Note: Conditional operator "?" evaluates an expression as true or false, returning 1st value if true, and 2nd value if false
	//Note: Syntax is condition ? result1 : result2
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") = (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") >= 0.0) ? beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") : 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") = (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") >= 0.0) ? beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") : 0.0;
}

void TreeInterceptionCalc::Interception_Rain_Tree(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_Tree_m2pm2 (m2/m2) is a combined LAI and BAI per unit of land
	double LAI_BAI_Tree_m2pm2 = input->LAI_BAI_Tree_m2_p_m2[timeStep];
	//ThroughFall_Canopy_frac (the tree throughfall coefficient) =  e^(-0.7 * LAI) becomes 1 with no LAI or tree cover
	//Note: 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.7 to trees
	double expTerm = exp(-0.7 * LAI_BAI_Tree_m2pm2);
	double ThroughFall_Canopy_frac = (1 - (1 - expTerm));

	// Compute initial values before updating beC->by_key
	//ThroughFall_Rain_TreeCanopy_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_TreeCanopy_m") = ThroughFall_Canopy_frac * input->Precip_mpdt[timeStep];
	//Interception_Rain_TreeCanopy_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_TreeCanopy_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_TreeCanopy_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_TreeCanopy_m")) {
		//SnowMelt_Rain_m (m) set equal to Storage_Snow_TreeCanopy_m (m)
		SnowMelt_Rain_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m");
	}

	//Storage_Snow_TreeCanopy_m (m) is reduced by SnowMelt_Rain_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") - SnowMelt_Rain_m;

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

	// Retrieve maximum canopy storage once
	//If Storage_Rain_TreeCanopy_m (m) greater than Storage_TreeCanopy_Max_m (m) then adjust interception and throughfall 
	if (Storage_Rain_TreeCanopy_m > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m")) {
		//Interception_Rain_TreeCanopy_m (m) is set to Storage_TreeCanopy_Max_m (m) minus Storage_Rain_TreeCanopy_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_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") - SnowMelt_Rain_m;
		//ThroughFall_Rain_TreeCanopy_m (m) is increased by Storage_Rain_TreeCanopy_m (m) and reduced by Storage_TreeCanopy_Max_m (m)
		//Note: Reduction by Storage_TreeCanopy_Max_m accounts for only removing what is not held in maximum storage
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_TreeCanopy_m") + Storage_Rain_TreeCanopy_m - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m");
		//ThroughFall_Snow_TreeCanopy_m (m) is increased by Storage_Snow_TreeCanopy_m (m) that falls as liquid equivalent
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m");
		//Storage_Snow_TreeCanopy_m (m) set to zero after throughfall clears load
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") = 0.0;
		//Storage_Rain_TreeCanopy_m (m) set to Storage_TreeCanopy_Max_m (m)
		Storage_Rain_TreeCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m");
	}

	//Storage_Rain_TreeCanopy_m (m) updated as <= Storage_TreeCanopy_Max_m
	//Note: Storage_Rain_TreeCanopy_m next sent to TreeEvaporationCalc.cpp
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") = Storage_Rain_TreeCanopy_m;
	//Melt_Snow_TreeCanopy_m (m) is snow melted in tree canopy due to rain, later added to that melted by radiation
	beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_TreeCanopy_m") = SnowMelt_Rain_m;
	//Control against negative values
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_TreeCanopy_m") = max(0.0, beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Rain_TreeCanopy_m"));
}

//Interception_Snow_Tree function will handle snow interception in the Tree canopy
//Note: Consider refactor to improve equations for SnowDensity_kg_p_m3 and Storage_SnowCrystal_Leaf_Max_m
void TreeInterceptionCalc::Interception_Snow_Tree(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//LAI_BAI_Tree_m2pm2 (m2/m2) is a combined LAI and BAI per unit of land
	double LAI_BAI_Tree_m2pm2 = input->LAI_BAI_Tree_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;
	// calculates the current snow storage max - Storage_SnowCrystal_TreeCanopy_Max_m
	double Storage_SnowCrystal_TreeCanopy_Max_m = Storage_SnowCrystal_Leaf_Max_m * LAI_BAI_Tree_m2pm2;
	//Note: Consider refactor to explore logic for division by 1000; was Storage_SnowCrystal_TreeCanopy_Max_m considered a volume?
	double Storage_Snow_TreeCanopy_Max_m = Storage_SnowCrystal_TreeCanopy_Max_m * (SnowDensity_kg_p_m3/1000);

	//ThroughFall_Canopy_frac (the tree throughfall coefficient) =  e^(-0.7 * LAI) becomes 1 with no LAI or tree cover
	//Note: 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.7 to trees
	double ThroughFall_Canopy_frac = exp(-0.7 * LAI_BAI_Tree_m2pm2);
	// transfers precipitation input (snow since temp < 32F)
	double Precipitation_Snow_m = input->Precip_mpdt[timeStep];
	
	//Note: Theory from Valente et al., (1997)
	//ThroughFall_Snow_TreeCanopy_m (m) equals ThroughFall_Canopy_frac * Precipitation_Snow_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") = ThroughFall_Canopy_frac * Precipitation_Snow_m;
	//Interception_Snow_TreeCanopy_m (m) equals (1 - ThroughFall_Canopy_frac) * Precipitation_Snow_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_TreeCanopy_m") = (1 - ThroughFall_Canopy_frac) * Precipitation_Snow_m;

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

	//If Storage_Snow_TreeCanopy_m (m) greater than Storage_Snow_TreeCanopy_Max_m (m) then reduce interception, increase throughfall
	if (Storage_Snow_TreeCanopy_m > Storage_Snow_TreeCanopy_Max_m)	{
		// CurrentTreeInterceptedSWE = Storage_Snow_TreeCanopy_Max_m - incoming TreeStoredSWE.  So the CurrentTreeInterceptedSWE = available SWE space
		beC->by_key(MapPixel_ID, DataFolder_ID, "Interception_Snow_TreeCanopy_m") = Storage_Snow_TreeCanopy_Max_m - beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m");	
		// SWE Throughfall  = calculated SWE throughfall + (temp stored SWE - Storage_TreeCanopy_Max_m) (the amount of water over Storage_TreeCanopy_Max_m)
		beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") + (Storage_Snow_TreeCanopy_m - Storage_Snow_TreeCanopy_Max_m);
		// Reset the tempTreeStoredSWE to be Storage_Snow_TreeCanopy_Max_m
		Storage_Snow_TreeCanopy_m = Storage_Snow_TreeCanopy_Max_m;
	}  
	
	//Storage_Snow_TreeCanopy_m (m) set to Storage_Snow_TreeCanopy_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") = Storage_Snow_TreeCanopy_m;

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