#include "SnowAblationUnderTreeCalc.h"

//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
//Note: Lundquist et al. (2021) updates snow unloading and melt routines. 
//Sanow, J. E., Fassnacht, S. R., Kamin, D. J., Sexstone, G. A., Bauerle, W. L., & Oprea, I. (2018). Geometric Versus Anemometric Surface Roughness for a Shallow Accumulating Snowpack. Geosciences, 8(12), 463. 
//Note: Sanow et al. (2018) demonstrate methods to adjust snow roughness length based on snow depth; median value z0_snow = 0.024 used by UCAR Community Land Model
//Sexstone, G. A., Clow, D. W., Fassnacht, S. R., Liston, G. E., Hiemstra, C. A., Knowles, J. F., & Penn, C. A. (2018). Snow Sublimation in Mountain Environments and Its Sensitivity to Forest Disturbance and Climate Warming. Water Resources Research, 54(2), 1191-1211. doi:https://doi.org/10.1002/2017WR021172
//Note: Sexstone et al. (2018) use SnowModel and demonstrate relative rates of sublimation, melt, and wind on snow ablation for ground and canopy. Canopy loss not always > Ground loss.
//USACE, 1956. Snow Hydrology: Summary Report of the Snow Investigations. US Army Corps of Engineers North Pacific Division Report NTIS PB-151660, Portland, Oregon.
//USACE, 1998. Engineering and Design: Runoff From Snowmelt. US Army Corps of Engineers Report Engineer Manual 1110-2-1406, Washington, D.C.
//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


//SnowAblationUnderTreeCalc::SnowAblationUnderShortVeg 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. 
void SnowAblationUnderTreeCalc::SnowAblationUnderTree(Inputs *input, CompactRagged* beC, WeatherProcessor* WeatherPro, 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 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;
	//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;
	//Height_WindMeasurement_m (m) is reference height for wind speed measurements
	double Height_WeatherStationWindSensor_m = input->SimulationNumericalParams["Height_Sensor_Wind_m"];
	//Height_TairMeasurement_m (m) is reference height for air temperature and relative humidity measurements
	double Height_WeatherStationTemperatureSensor_m = input->SimulationNumericalParams["Height_Sensor_Tair_m"];
	//ThroughFall_Rain_TreeCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_TreeCanopy_m")
	double ThroughFall_Rain_TreeCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Rain_TreeCanopy_m");
	// if "Storage_Rain_TreeCanopy_m" not registered, set it to 0 to prevent access violation
	if (!beC->has_var("Storage_Snow_underTreeCanopy_m")) {
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") = 0.0;
	}
	//Storage_Snow_underTreeCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m")
	double Storage_Snow_underTreeCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "ThroughFall_Snow_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m");
	//beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") = Storage_Snow_underTreeCanopy_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") = Storage_Snow_underTreeCanopy_m;


	//SnowMelt_RadiationShortWave_m_p_hr = 0; initialized to zero
	double SnowMelt_RadiationShortWave_m_p_hr = 0;
	//If Radiation_Shortwave_Diffuse_Wpm2 > 0 then
	if (input->Radiation_Shortwave_Diffuse_Wpm2[timeStep] > 0) {
		//SnowMelt_RadiationShortWave_m_p_hr = 0.0000317; constant rate of melt due to shortwave radiation under the canopy
		//Note: From Eq 11 of Yang et al. (2011) and attributed to USACE (1956)
		SnowMelt_RadiationShortWave_m_p_hr = 0.0000317;
	}
	//SnowMelt_RadiationLongWave_m_p_hr = 0.0000553 * Tair_C; Coefficient was was 0.0000307 for Tair_F, scaled by 9/5
	//Note: From Eq 11 of Yang et al. (2011) and attributed to USACE (1956)
	double SnowMelt_RadiationLongWave_m_p_hr = 0.0000553 * Tair_C;
	//SnowMelt_Radiation_m_p_hr = SnowMelt_RadiationShortWave_m_p_hr + SnowMelt_RadiationLongWave_m_p_hr
	double SnowMelt_Radiation_m_p_hr = SnowMelt_RadiationShortWave_m_p_hr + SnowMelt_RadiationLongWave_m_p_hr;

	//SnowMelt_SensibleHeat_m_p_hr (m/hr) is the snowmelt by sensible heat from Eq 7 of Yang et al. (2011), which had error in coefficient = 0.00000792
	//Note:	Eq 7 of Yang et al. (2011) derived from Eq 5-10 (USACE, 1998) gave LHS term Mc in units in/day, w RHS coefficient 0.00629, product of heights in ft^2, temp in F, winds in miles/hr
	//Note:	Eq 7 (Yang et al., 2011) is Eq 5-10 for units m/hr, and RHS product of heights in m^2, temp in C, winds in m/s
	//Note: ... Coeff_01_Mh_m_p_hr_C_m_p_s = 0.00629 in/day ft, F, mile/hr * (1/12)*0.3048*(1/24) = 0.0000066569 m/hr and ft, F, mile/hr
	//Note: ... Coeff_01_Mh_m_p_hr_C_m_p_s = 0.0000066569 m/hr and ft, F, mile/hr * (1/(5280/1*0.3048/1*1/3600)) = 0.000014891 m/hr and ft, F, m/s
	//Note: ... Coeff_01_Mh_m_p_hr_C_m_p_s = 0.000014891 m/hr and ft, F, m/s * (1/0.3048^2)^(-1/6) = 0.000010021 m/hr and m, F, m/s
	//Note: ... Coeff_01_Mh_m_p_hr_C_m_p_s = 0.000010021 m/hr and m, F, m/s * 9/5 degree_F/degree_C = 0.0000180387 m/hr and m, C, m/s
	double Coeff_01_Mh_m_p_hr_C_m_p_s = 0.0000180387;
	double SnowMelt_SensibleHeat_m_p_hr = Coeff_01_Mh_m_p_hr_C_m_p_s * Ratio_AtmPres_to_StdPres * pow(Height_WeatherStationTemperatureSensor_m * Height_WeatherStationWindSensor_m, -(1.0 / 6.0)) * Tair_C * WindSpd_mps;

	//SnowMelt_LatentHeat_m_p_hr (m/hr) is the snowmelt by latent heat from Eq 8 of Yang et al. (2011), which had error in coefficient = 0.0000223 from USACE (1998)
	//Note:	Eq 8 of Yang et al. (2011) derived from 5-12 to 5-14 (USACE, 1998), with Eq 5-14 using incorrect coefficient = 0.0065 due to constant height
	//Note:	Eq 5-12 (USACE, 1998) gave LHS term SnowMelt_LatentHeat_m_p_hr in units in/day, w RHS coefficient 0.054, ...
	//Note: ... product of heights in ft^2, vapor pressure in mb, winds in miles/hr
	//Note:	Eq 5-13 (USACE, 1998) computes vapor pressure (mb) from air temperature (F) using a coefficient 0.339, changing (ea-es) to Td
	//Note:	Eq 5-14 (USACE, 1998) created by substituting a fixed height into Eq 5-12, and combining Eq 5-12 with Eq 5-13
	//Note: ... Coeff_01_Me_m_p_hr_C_m_p_s = 0.054 in/day and ft, Mb, mile/hr * (1/12)*0.3048*(1/24) = 0.0000571500 m/hr and ft, Mb, mile/hr
	//Note: ... Coeff_01_Me_m_p_hr_C_m_p_s = 0.0000571500 m/hr and ft, Mb, mile/hr * (1/(5280/1*0.3048/1*1/3600)) = 0.000127841 m/hr and ft, Mb, m/s
	//Note: ... Coeff_01_Me_m_p_hr_C_m_p_s = 0.000127841 m/hr and ft, Mb, m/s * (1/0.3048^2)^(-1/6) = 0.000086035 m/hr and m, Mb, m/s
	//Note: ... Coeff_01_Me_m_p_hr_C_m_p_s = 0.000086035 m/hr and m, Mb, m/s * 0.339 degree F / mbar vapor (see Eq 5-13) = 0.000029166 m/hr and m, F, m/s
	//Note: ... Coeff_01_Me_m_p_hr_C_m_p_s = 0.000029166 m/hr and m, F, m/s * 9/5 degree_F/degree_C = 0.0000524985 m/hr and m, C, m/s
	double Coeff_01_Me_m_p_hr_C_m_p_s = 0.0000524985;
	double SnowMelt_LatentHeat_m_p_hr = Coeff_01_Me_m_p_hr_C_m_p_s * pow(Height_WeatherStationTemperatureSensor_m * Height_WeatherStationWindSensor_m, -(1.0 / 6.0)) * Tdew_C * WindSpd_mps;

	//SnowMelt_GroundHeat_m_p_hr = 0.0000277 m/hr
	//Note: Consider refactor for varying value based on ground heat flux
	//Note: SnowMelt_GroundHeat_m_p_hr (m⁄hr) is the snowmelt by the ground heat, from Eq 9 of Yang et al. (2011), 
	//Note: USACE (1998) page 5-6, section 5-2.i. recommends use of a constant value between SnowMelt_GroundHeat_m_p_hr = 0.0254 to 0.0762 cm/day, 
	//Note: ... which is 0.00001058 to 0.00003175 m/hr, an average value of 0.000021 m/hr. 
	//Note: USACE (1998) page 2-5, section 2-2.d. recommends value between Qg = 0 to 5 J/m^2/s, converted in Eq 2-2, M = 0.0000277 m/hr
	//Note: Eq 2-2 of USACE (1998) computed SnowMelt_GroundHeat_mm_p_s = Q / (B + F * rho), where Q = 2.5E-3 kJ/m2/s power to melt by ground heat, ...
	//Note: ... B = 0.97 thermal quality of snow, F = 334.9 kJ/kg latent heat of fusion, rho = 1000 kg/m3 density of water, and ...
	//Note: SnowMelt_GroundHeat_mm_p_s = 7.69579E-09 mm/s => SnowMelt_GroundHeat_m_p_hr = 0.0000277 mm/hr
	double SnowMelt_GroundHeat_m_p_hr = 0.0000277;

	//SnowMelt_ViaRainfall_m_p_hr (m) is the snow melt by rain, initializing to 0.0
	double SnowMelt_ViaRainfall_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 * ThroughFall_Rain_TreeCanopy_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_ViaRainfall_m_p_hr = Coeff_01_Mp_m_p_h_C * ThroughFall_Rain_TreeCanopy_m * Tair_C;
	}

	//Coeff_wind_in_forest = 0.3 + (1.0 - beC->by_key(MapPixel_ID, DataFolder_ID, "TreeCover_frac")) * 0.7
	//Note: Equation by Endreny given Yang et al. (2011) stating Coeff_wind_in_forest ranges from 0.3 in high forest density to 1.0 in low forest density 
	double Coeff_wind_in_forest = 0.3 + (1.0 - beC->by_key(MapPixel_ID, DataFolder_ID, "TreeCover_frac")) * 0.7;
	//SnowMelt_SensibleHeat_and_LatentHeat_m_p_hr = Coeff_wind_in_forest * (SnowMelt_SensibleHeat_m_p_hr + SnowMelt_LatentHeat_m_p_hr)
	//Note: From Eq 12 of Yang et al. (2011) 
	double SnowMelt_SensibleHeat_and_LatentHeat_m_p_hr = Coeff_wind_in_forest * (SnowMelt_SensibleHeat_m_p_hr + SnowMelt_LatentHeat_m_p_hr);

	//Melt_Snow_underTreeCanopy_m_p_hr = (SnowMelt_Radiation_m_p_hr + SnowMelt_SensibleHeat_and_LatentHeat_m_p_hr + SnowMelt_GroundHeat_m_p_hr + SnowMelt_ViaRainfall_m_p_hr)
	//Note: Modified Equation 5 of Yang et al. (2011) for under canopy, combining Sensible and Latent Heat into a single term
	double Melt_Snow_underTreeCanopy_m_p_hr = (SnowMelt_Radiation_m_p_hr + SnowMelt_SensibleHeat_and_LatentHeat_m_p_hr + SnowMelt_GroundHeat_m_p_hr + SnowMelt_ViaRainfall_m_p_hr);
	//Melt_Snow_underTreeCanopy_m is Melt_Snow_underTreeCanopy_m_p_hr * input->SimulationTimeStep_Duration_sec[timeStep] * Ratio_Hour_to_Second
	//Note: Conversion from m/hr to m based on SimulationTimeStep_Duration_sec
	beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m") = Melt_Snow_underTreeCanopy_m_p_hr * input->SimulationTimeStep_Duration_sec[timeStep] * Ratio_Hour_to_Second;

	//controls for zero
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m") < 0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m") = 0.0; }

	//If Melt_Snow_underTreeCanopy_m > Storage_Snow_underTreeCanopy_m then set to Storage_Snow_underTreeCanopy_m
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m") > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m")) {
		beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m");
	}

	//Storage_Snow_underTreeCanopy_m is beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m")
	Storage_Snow_underTreeCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m");
	//Storage_Snow_underTreeCanopy_m is Storage_Snow_underTreeCanopy_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") = Storage_Snow_underTreeCanopy_m;

	//SublimationSnow_Potential_m (m) defined as SublimationSnow_Potential_m, which is generated by functor calculateGroundSublimation in OpenWaterEvaporationCalc::calculateGroundSublimationStatistical or calculateGroundSublimationTemperature
	//Note: From SimulationCoordinator::runCalculations calculateGroundSublimation functor was called to get SublimationSnow_Potential_m (m) 
	double SublimationSnow_Potential_m = beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_Potential_m");

	//If SublimationSnow_Potential_m > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") then set to Storage_Snow_underTreeCanopy_m
	if (SublimationSnow_Potential_m > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m")) {
		SublimationSnow_Potential_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m");
	}

	//Storage_Snow_underTreeCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") - SublimationSnow_Potential_m
	Storage_Snow_underTreeCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") - SublimationSnow_Potential_m;

	//Storage_Snow_underTreeCanopy_m is Storage_Snow_underTreeCanopy_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_underTreeCanopy_m") = Storage_Snow_underTreeCanopy_m;

	//Precipitation_RainSnowMelt_underTreeCanopy_m is beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m") + ThroughFall_Rain_TreeCanopy_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Precipitation_RainSnowMelt_underTreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_underTreeCanopy_m") + ThroughFall_Rain_TreeCanopy_m;

	//SublimationSnow_underTreeCanopy_m (m) is SublimationSnow_Potential_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_underTreeCanopy_m") = SublimationSnow_Potential_m;

	//Control for negative values
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_underTreeCanopy_m") < 0) { beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_underTreeCanopy_m") = 0; }
	
	// if "Storage_Rain_TreeCanopy_m" not registered, set it to 0 to prevent access violation
	if (!beC->has_var("Storage_Rain_TreeCanopy_m")) {
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") = 0.0;
	}

	//Storage_RainSnow_TreeCanopy_m is beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m")
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_RainSnow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m");
}

