#include "SnowAblationOpenAreaCalc.h"
#include "../Evaporation/OpenWaterEvaporationCalc.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

//SnowAblationOpenAreaCalc::SnowAblationOnGround function computes ablation from melt first, then sublimation
void SnowAblationOpenAreaCalc::SnowAblationOnGround(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 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;
	//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"];

	//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_RadiationShortWave_m_p_hr = Coeff_01_Mn_m_p_hr_W_p_m2 * Radiation_Shortwave_Direct_Wpm2 * (1 - Albedo_Snow_frac)
	//Note: Snowmelt by shortwave radiation defined by Eq 6 1st RHS in Yang et al. (2011), derived from Eq 5-2 (USACE, 1998) ...
	//Note: ... Coeff_01_Mn_m_p_hr_W_p_m2 = 0.00308 mm/day kJ/m2/day * (1/1000)*(1/24) = 0.000000128 m/hr kJ/m2/day
	//Note: ... Coeff_01_Mn_m_p_hr_W_p_m2 = 0.000000128 m/hr kJ/m2/day * (1/1000)*(24/1)*(3600/1) = 0.00001109 m/hr W/m2
	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
	//Note: Consider refactor to use cloud cover from NOAA NCEI
	if (Precipitation_m > 0) {
		//SnowMelt_RadiationLongWave_m_p_hr = Coeff_01_Ml_cloudy_m_p_hr_C * Tair_C
		//Note: Snowmelt by long wave radiation on cloudy day, replacing 2nd RHS terms of Eq 6 in Yang et al. (2011) ...
		//Note: Consider refactor to use actual longwave radiation
		//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: ... Coeff_01_Ml_cloudy_m_p_hr_C = 0.029 in/day (F)  * (1/12)*(0.3048/1)*(1/24) = 0.000030692 m/hr (F)
		//Note: ... Coeff_01_Ml_cloudy_m_p_hr_C = 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: ... Coeff_02_Ml_clear_m_p_hr_C = 0.804 in/day (F) * (1/12)*(0.3048/1)*(1/24) = 0.0008509 m/hr (F) = 0.0008509 m/hr (C)  
		//Note: ... Coeff_01_Ml_clear_m_p_hr_C = 0.0212 in/day (F)  * (1/12)*(0.3048/1)*(1/24) = 0.000022437 m/hr (F)
		//Note: ... Coeff_01_Ml_clear_m_p_hr_C = 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;

	//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) {
		//Precipitation_Rain_m (m) set to Precipitation_m
		Precipitation_Rain_m = Precipitation_m;
		//SnowMelt_ViaRainfall_m_p_hr = Coeff_01_Mp_m_p_h_C * Precipitation_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 * Precipitation_m * Tair_C;
	}
	
	//If the Tair_C air temperature is below 0C or 32F, then precipitation treated as frozen snow and it accumulated on surface
	if(Tair_C <= 0.0) {
		//Precipitation_Rain_m (m) set to zero
		Precipitation_Rain_m = 0.0;
		//Storage_Snow_noCanopy_m (m) increased by additional snow fall as liquid equivalent
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") + Precipitation_m;
	}

	//Melt_Snow_no_Canopy_m_p_hr = (SnowMelt_Radiation_m_p_hr + SnowMelt_SensibleHeat_m_p_hr + SnowMelt_LatentHeat_m_p_hr + SnowMelt_GroundHeat_m_p_hr + SnowMelt_ViaRainfall_m_p_hr)
	//Note: Equation 5 of Yang et al. (2011) and based on USACE (1998; 1956)
	double Melt_Snow_no_Canopy_m_p_hr = (SnowMelt_Radiation_m_p_hr + SnowMelt_SensibleHeat_m_p_hr + SnowMelt_LatentHeat_m_p_hr + SnowMelt_GroundHeat_m_p_hr + SnowMelt_ViaRainfall_m_p_hr);
	//Melt_Snow_noCanopy_m is Melt_Snow_no_Canopy_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_noCanopy_m") = Melt_Snow_no_Canopy_m_p_hr * input->SimulationTimeStep_Duration_sec[timeStep] * Ratio_Hour_to_Second;

	//If Melt_Snow_noCanopy_m (m) less than 0, then set it to 0
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_noCanopy_m") < 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_noCanopy_m") = 0.0;	}

	//If Melt_Snow_noCanopy_m (m) is > Storage_Snow_noCanopy_m, set it equal to Storage_Snow_noCanopy_m
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_noCanopy_m") > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m")) {
		beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_noCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m");
	}

	//Storage_Snow_noCanopy_m (m) is updated to reflect the loss of the melted snow
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_noCanopy_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 (m) amount of snow that could sublimate is > Storage_Snow_noCanopy_m or snow available then set equal to Storage_Snow_noCanopy_m 
	if (SublimationSnow_Potential_m > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m")) {
		SublimationSnow_Potential_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m");
	}
	
	//Storage_Snow_noCanopy_m (m) is updated to reflect the loss of the sublimated snow, SublimationSnow_Potential_m
	double Storage_Snow_noCanopy_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") - SublimationSnow_Potential_m;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") = Storage_Snow_noCanopy_m;

	//Precipitation_RainSnowMelt_noCanopy_m (m) is updated to equal Melt_Snow_noCanopy_m + Precipitation_Rain_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Precipitation_RainSnowMelt_noCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Melt_Snow_noCanopy_m") + Precipitation_Rain_m;

	//SublimationSnow_noCanopy_m (m) updated with SublimationSnow_Potential_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_noCanopy_m") = SublimationSnow_Potential_m;

	//Control for negative values
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_noCanopy_m") < 0) { beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_noCanopy_m") = 0; }

}
