#include "OpenWaterEvaporationCalc.h"

//OpenWaterEvaporationCalc::calculateImpDepEvapTemperature function is called as calculateImpDepEvap functor to compute Evaporation_ImperviousPondedWater_m (m) for SpatialTemperatureHydro
void OpenWaterEvaporationCalc::calculateImpDepEvapTemperature(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//If timeStep equals 1, then call statistical calculation based on weather input data, given the spatial calculation has not yet been called
	if (timeStep == 0) {
		//Call OpenWaterEvaporationCalc::calculateImpDepEvapStatistical to estimate evaporation
		OpenWaterEvaporationCalc::calculateImpDepEvapStatistical(input, beC, DataDrawer_ID, DataFolder_ID, timeStep);
		//return is C++ function to leave the .cpp
		return;
	}

	//Initialize LAI_BAI_Tree_m2_p_m2 
	double LAI_BAI_Tree_m2_p_m2 = input->LAI_BAI_Tree_m2_p_m2[timeStep];
	//Note: If LAI_BAI_Tree_m2_p_m2 < 1, then set to 1 so that its position as denominator does not enlarge TreeLEE_W_p_m2 when / LAI_BAI_Tree_m2_p_m2
	if (LAI_BAI_Tree_m2_p_m2 < 1) { LAI_BAI_Tree_m2_p_m2 = 1; }

	//Temperature_Air_C (C) is converted from folder Tair_K (K) 
	double Temperature_Air_C = beC->by_key(DataDrawer_ID, DataFolder_ID, "Tair_K") - 273.15;

	//Calculate HeatMetrics with call to HeatMetrics_Calc functions
	//HeatMetrics_Calc HeatMetricsCalc(input) creates pointer to access HeatMetrics_Calc functions with sending input
	HeatMetrics_Calc HeatMetricsCalc(input);

	//DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"))
	double DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(DataDrawer_ID, DataFolder_ID, "Tair_K"));

	//LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(Tair_K)
	double LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(beC->by_key(DataDrawer_ID, DataFolder_ID, "Tair_K"));

	//Evaporation_ImperviousPondedWater_m_p_s (m/s) initialized
	double Evaporation_ImperviousPondedWater_m_p_s = 0;

	//Evaporation_ImperviousPondedWater_m_p_s (m/s) derived as latent energy (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3)
	//Note: Consider that Evaporation_ImperviousPondedWater_m occurs from depression storage on impervious cover with and without tree cover but ...
	//Note:	Calculation of Evaporation_ImperviousPondedWater_m requires use of TreeLE_W_p_m2 for cases where TreeCanopyCover_overImpervious_frac > 0
	//Note: The TreeLE_W_p_m2 / LAI_BAI_Tree_m2_p_m2 is a fraction of energy allowed to evaporate water from ponded water under the tree canopy
	//Note: The TreeLE_W_p_m2 - (TreeLE_W_p_m2 / LAI_BAI_Tree_m2_p_m2) will be available for evaporation of water in the tree canopy and its soil
	//Note: The TreeLEE_W_p_m2 term is not used alone, as that may go to zero while there is still water in Storage_ImperviousPondedWater_m
	//Note: This is exploratively using a single layer energy model to evaporate depression storage from understory areas, maintaining an energy balance
	//Note: RHS numerator is W/m2 = J/s/m2, and RHS denominator is (J/kg) * (kg/m3) = J/m3, therefore RHS quotient is J/s/m2 / J/m3 = m/s.
	Evaporation_ImperviousPondedWater_m_p_s = (beC->by_key(DataDrawer_ID, DataFolder_ID, "ImpLE_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "ImperviousCover_noTreeCanopy_frac") +
		(beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeLE_W_p_m2") / LAI_BAI_Tree_m2_p_m2) * beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeCanopyCover_overImpervious_frac")) /
		(LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);

	//Evaporation_ImperviousPondedWater_m_p_s (m/s) divided by all cover types so it is scaled to cover area
	Evaporation_ImperviousPondedWater_m_p_s = input->SafeDivide(Evaporation_ImperviousPondedWater_m_p_s, (beC->by_key(DataDrawer_ID, DataFolder_ID, "ImperviousCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeCanopyCover_overImpervious_frac")));

	//Control for negative values
	Evaporation_ImperviousPondedWater_m_p_s = Inputs::forceZeroIfLessThanOrAlmostEqualZero(Evaporation_ImperviousPondedWater_m_p_s, Epsilon_Tolerance_1E_negative15);

	//Consider refactor to compute this as two parts, impervious without tree cover and impervious with tree cover
	//Evaporation_ImperviousPondedWater_m (m) is product of Evaporation_ImperviousPondedWater_m_p_s (m/s) and SimulationTimeStep_Duration_sec[timeStep]
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") = Evaporation_ImperviousPondedWater_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
}

//OpenWaterEvaporationCalc::calculateImpDepEvapStatistical function is called as calculateImpDepEvap functor to compute Evaporation_ImperviousPondedWater_m (m) for StatisticalHydro
void OpenWaterEvaporationCalc::calculateImpDepEvapStatistical(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//BulkArea Type will have depth constrained evaporation
	if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "BulkArea") {
		//If Storage_ImperviousPondedWater_m / Storage_ImperviousPondedWater_Max_m < 1 then regulate Evaporation_ImperviousPondedWater_m (m) = fraction of PotentialEvaporation_WaterOnGround_m
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_ImperviousPondedWater_m") < beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_ImperviousPondedWater_Max_m")) {
			//Evaporation_ImperviousPondedWater_m (m) is set to fraction of PotentialEvaporation_WaterOnGround_m based on relative depth of water between fraction of 0 and 1 of maximum depth
			beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep] *
				input->SafeDivide(beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_ImperviousPondedWater_m"), beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_ImperviousPondedWater_Max_m"));
		}
		//If Storage_ImperviousPondedWater_m / Storage_ImperviousPondedWater_Max_m >= 1 then regulate Evaporation_ImperviousPondedWater_m (m) = PotentialEvaporation_WaterOnGround_m
		else
			beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep];
	}
	//GI Type will have unconstrained evaporation
	else {
		beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep];
	}
}

//OpenWaterEvaporationCalc::calculatePerDepEvapTemperature function is called as calculatePerDepEvap functor to compute Evaporation_PerviousPondedWater_m (m) for SpatialTemperatureHydro
void OpenWaterEvaporationCalc::calculatePerDepEvapTemperature(Inputs* input, CompactRagged* beC,int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//If timeStep equals 1, then call statistical calculation based on weather input data, given the spatial calculation has not yet been called
	if (timeStep == 0) {
		//Call OpenWaterEvaporationCalc::calculatePerDepEvapStatistical to estimate evaporation
		OpenWaterEvaporationCalc::calculatePerDepEvapStatistical(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		//return is C++ function to leave the .cpp
		return;
	}

	//Initialize LAI_BAI_Tree_m2_p_m2 and LAI_BAI_SVeg_m2_p_m2 
	double LAI_BAI_Tree_m2_p_m2 = input->LAI_BAI_Tree_m2_p_m2[timeStep];
	double LAI_BAI_SVeg_m2_p_m2 = input->LAI_BAI_SVeg_m2_p_m2[timeStep];
	//Note: If LAI_BAI_Tree_m2_p_m2 < 1, then set to 1 so that its position as denominator does not enlarge TreeLEE_W_p_m2 when / LAI_BAI_Tree_m2_p_m2
	if (LAI_BAI_Tree_m2_p_m2 < 1) { LAI_BAI_Tree_m2_p_m2 = 1; }
	//Note: If LAI_BAI_SVeg_m2_p_m2  < 1, then set to 1 so that its position as denominator does not enlarge SVegLEE_W_p_m2 when / LAI_BAI_SVeg_m2_p_m2 
	if (LAI_BAI_SVeg_m2_p_m2 < 1) { LAI_BAI_SVeg_m2_p_m2 = 1; }

	//Temperature_Air_C (C) is converted from folder Tair_K (K) 
	double Temperature_Air_C = beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K") - 273.15;

	//Calculate HeatMetrics with call to HeatMetrics_Calc functions
	//HeatMetrics_Calc HeatMetricsCalc(input) creates pointer to access HeatMetrics_Calc functions with sending input
	HeatMetrics_Calc HeatMetricsCalc(input);

	//DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"))
	double DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"));

	//LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(Tair_K)
	double LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"));

	//Evaporation_PerviousPondedWater_m_p_s (m/s) initialized
	double Evaporation_PerviousPondedWater_m_p_s = 0;

	//Evaporation_PerviousPondedWater_m_p_s (m/s) derived as latent energy (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3)
	//Evaporation_PerviousPondedWater_m_p_s (m/s) derived as latent energy (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3)
	//Note: Consider that Evaporation_PerviousPondedWater_m occurs from depression storage on pervious cover with and without tree or short veg cover but ...
	//Note:	Calculation of Evaporation_PerviousPondedWater_m requires use of TreeLE_W_p_m2 and SVegLE_W_p_m2 for cases where Storage_PerviousPondedWater_m > 0
	//Note: The TreeLE_W_p_m2 / LAI_BAI_Tree_m2_p_m2 is a fraction of energy allowed to evaporate water from ponded water under the tree canopy
	//Note: The TreeLE_W_p_m2 - (TreeLE_W_p_m2 / LAI_BAI_Tree_m2_p_m2) will be available for evaporation of water in the tree canopy and its soil
	//Note: The TreeLEE_W_p_m2 term is not used alone, as that may go to zero while there is still water in Storage_PerviousPondedWater_m
	//Note: The SVegLE_W_p_m2 / LAI_BAI_SVeg_m2_p_m2 is a fraction of energy allowed to evaporate water from ponded water under the short veg canopy
	//Note: The SVegLE_W_p_m2 - (SVegLE_W_p_m2 / LAI_BAI_SVeg_m2_p_m2) will be available for evaporation of water in the short veg canopy and its soil
	//Note: The SVegLEE_W_p_m2 term is not used alone, as that may go to zero while there is still water in Storage_PerviousPondedWater_m
	//Note: This is creatively using a single layer energy model to evaporate depression storage from understory areas, maintaining an energy balance
	//Note: RHS numerator is W/m2 = J/s/m2, and RHS denominator is (J/kg) * (kg/m3) = J/m3, therefore RHS quotient is J/s/m2 / J/m3 = m/s.
	Evaporation_PerviousPondedWater_m_p_s = (beC->by_key(MapPixel_ID, DataFolder_ID, "SoilLE_W_p_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "SoilCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "ImpLE_W_p_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "PermeablePavementCover_noTreeCanopy_frac") + (beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLE_W_p_m2") / LAI_BAI_SVeg_m2_p_m2) * beC->by_key(MapPixel_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") + (beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLE_W_p_m2") / LAI_BAI_Tree_m2_p_m2) * beC->by_key(MapPixel_ID, DataFolder_ID, "TreeCanopyCover_overPervious_frac")) / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);

	//Evaporation_PerviousPondedWater_m_p_s (m/s) divided by all cover types so it is scaled to cover area
	Evaporation_PerviousPondedWater_m_p_s = input->SafeDivide(Evaporation_PerviousPondedWater_m_p_s, (beC->by_key(MapPixel_ID, DataFolder_ID, "SoilCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "PermeablePavementCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "TreeCanopyCover_overPervious_frac")));

	//Control for negative values
	Evaporation_PerviousPondedWater_m_p_s = Inputs::forceZeroIfLessThanOrAlmostEqualZero(Evaporation_PerviousPondedWater_m_p_s, Epsilon_Tolerance_1E_negative15);

	//Evaporation_PerviousPondedWater_m (m) is product of Evaporation_PerviousPondedWater_m_p_s (m/s) and SimulationTimeStep_Duration_sec[timeStep]
	beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") = Evaporation_PerviousPondedWater_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
}

//OpenWaterEvaporationCalc::calculatePerDepEvapStatistical function is called as calculatePerDepEvap functor to compute Evaporation_PerviousPondedWater_m (m) for StatisticalHydro
void OpenWaterEvaporationCalc::calculatePerDepEvapStatistical(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//BulkArea Type will have depth constrained evaporation
	if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "BulkArea") {
		//If Storage_PerviousPondedWater_m < Storage_PerviousPondedWater_Max_m then regulate Evaporation_PerviousPondedWater_m (m) = fraction of PotentialEvaporation_WaterOnGround_m
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_PerviousPondedWater_m") < beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_PerviousPondedWater_Max_m")) {
			//Evaporation_PerviousPondedWater_m (m) is set to fraction of PotentialEvaporation_WaterOnGround_m based on relative depth of water between fraction of 0 and 1 of maximum depth
			beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep] * input->SafeDivide(beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_PerviousPondedWater_m"), beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_PerviousPondedWater_Max_m"));
		}
		//If Storage_PerviousPondedWater_m >= Storage_PerviousPondedWater_Max_m then regulate Evaporation_PerviousPondedWater_m (m) = PotentialEvaporation_WaterOnGround_m
		else {
			beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep];
		}
	}
	//GI Type will have unconstrained evaporation
	else {
		beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep];
	}
}

//OpenWaterEvaporationCalc::calculateWaterDepEvapTemperature function is called as calculateWaterDepEvap functor to compute Evaporation_WaterPondedWater_m (m) for SpatialTemperatureHydro
void OpenWaterEvaporationCalc::calculateWaterDepEvapTemperature(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//If timeStep equals 1, then call statistical calculation based on weather input data, given the spatial calculation has not yet been called
	if (timeStep == 0) {
		//Call OpenWaterEvaporationCalc::calculateWaterDepEvapStatistical to estimate evaporation
		OpenWaterEvaporationCalc::calculateWaterDepEvapStatistical(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		//return is C++ function to leave the .cpp
		return;
	}

	//Temperature_Air_C (C) is converted from folder Tair_K (K) 
	double Temperature_Air_C = beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K") - 273.15;

	//Calculate HeatMetrics with call to HeatMetrics_Calc functions
	//HeatMetrics_Calc HeatMetricsCalc(input) creates pointer to access HeatMetrics_Calc functions with sending input
	HeatMetrics_Calc HeatMetricsCalc(input);

	//DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"))
	double DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"));

	//LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(Tair_K)
	double LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"));

	//Evaporation_WaterPondedWater_m_p_s (m/s) initialized
	double Evaporation_WaterPondedWater_m_p_s = 0;

	//Evaporation_WaterPondedWater_m_p_s (m/s) derived from WaterLE_W_p_m2 (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3)
	//Note: RHS numerator is W/m2 = J/s/m2, and RHS denominator is (J/kg) * (kg/m3) = J/m3, therefore RHS quotient is J/s/m2 / J/m3 = m/s.
	Evaporation_WaterPondedWater_m_p_s = beC->by_key(MapPixel_ID, DataFolder_ID, "WaterLE_W_p_m2") / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);

	//Control for negative values
	Evaporation_WaterPondedWater_m_p_s = Inputs::forceZeroIfLessThanOrAlmostEqualZero(Evaporation_WaterPondedWater_m_p_s, Epsilon_Tolerance_1E_negative15);

	//Evaporation_WaterPondedWater_m (m) is product of Evaporation_WaterPondedWater_m_p_s (m/s) and SimulationTimeStep_Duration_sec[timeStep]
	beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") = Evaporation_WaterPondedWater_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
}

//OpenWaterEvaporationCalc::calculateWaterDepEvapStatistical function is called as calculateWaterDepEvap functor to compute Evaporation_WaterPondedWater_m (m) for StatisticalHydro
void OpenWaterEvaporationCalc::calculateWaterDepEvapStatistical(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//BulkArea Type will have depth constrained evaporation
	if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "BulkArea") {
		//If Storage_WaterPondedWater_m / Storage_WaterPondedWaterMax_m < 1 then regulate Evaporation_WaterPondedWater_m (m) = fraction of PotentialEvaporation_WaterOnGround_m
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_WaterPondedWater_m") < beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_WaterPondedWaterMax_m")) {
			//Evaporation_WaterPondedWater_m (m) is set to fraction of PotentialEvaporation_WaterOnGround_m based on relative depth of water between fraction of 0 and 1 of maximum depth
			beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep] * input->SafeDivide(beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_WaterPondedWater_m"), beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_WaterPondedWaterMax_m"));
		}
		//If Storage_ImperviousPondedWater_m / Storage_WaterPondedWaterMax_m >= 1 then regulate Evaporation_WaterPondedWater_m (m) = PotentialEvaporation_WaterOnGround_m
		else
			beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep];
	}
	//GI Type will have unconstrained evaporation
	else {
		beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") = input->PotentialEvaporation_WaterOnGround_m[timeStep];
	}
}

//OpenWaterEvaporationCalc::calculateGroundSublimationTemperature function is called as calculateGroundSublimation functor to compute SublimationSnow_Potential_m (m) for SpatialTemperatureHydro
void OpenWaterEvaporationCalc::calculateGroundSublimationTemperature(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//If timeStep equals 1, then call statistical calculation based on weather input data, given the spatial calculation has not yet been called
	if (timeStep == 0) {
		//Call OpenWaterEvaporationCalc::calculateGroundSublimationStatistical to estimate evaporation
		OpenWaterEvaporationCalc::calculateGroundSublimationStatistical(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		//return is C++ function to leave the .cpp
		return;
	}

	//If Storage_Snow_noCanopy_m > 0 then compute SublimationSnow_Potential_m
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") > 0) {
		double Ground_LE_W_p_m2 = 0.0;

		//Ground_LE_W_p_m2 (W/m2) is weighted average of all ground areas not under canopy, or zero
		Ground_LE_W_p_m2 = beC->by_key(MapPixel_ID, DataFolder_ID, "SoilLE_W_p_m2") * (beC->by_key(MapPixel_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "SoilCover_noTreeCanopy_frac")) +	beC->by_key(MapPixel_ID, DataFolder_ID, "ImpLE_W_p_m2") * (beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "PermeablePavementCover_noTreeCanopy_frac")) + beC->by_key(MapPixel_ID, DataFolder_ID, "WaterLE_W_p_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "WaterCover_noTreeCanopy_frac");

		//Ground_LE_W_p_m2 (W/m2) next divided by all cover types with sky view so it is scaled to cover area
		//Note: Scaling is used to get estimate of ground sublimation rate, which is later used for both open and vegegtation covered areas
		Ground_LE_W_p_m2 = input->SafeDivide(Ground_LE_W_p_m2, (beC->by_key(MapPixel_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "SoilCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "PermeablePavementCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousCover_noTreeCanopy_frac") + beC->by_key(MapPixel_ID, DataFolder_ID, "WaterCover_noTreeCanopy_frac")));

		//Temperature_Air_C (C) is converted from folder Tair_K (K) 
		double Temperature_Air_C = beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K") - 273.15;

		//Calculate HeatMetrics with call to HeatMetrics_Calc functions
		//HeatMetrics_Calc HeatMetricsCalc(input) creates pointer to access HeatMetrics_Calc functions with sending input
		HeatMetrics_Calc HeatMetricsCalc(input);

		//DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"))
		double DensityWater_kgpm3 = HeatMetricsCalc.Density_Water_kgpm3_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"));

		//LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(Tair_K)
		double LatentHeat_Vaporization_JpkgK = HeatMetricsCalc.LatentHeat_Vaporization_JpkgK_Calc(beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K"));

		//Evaporation_WaterPondedWater_m_p_s (m/s) initialized
		double SublimationSnow_Potential_m_p_s = 0;
		//SublimationSnow_Potential_m_p_s (m/s) derived from Ground_LE_W_p_m2 (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3)
		//Note: RHS numerator is W/m2 = J/s/m2, and RHS denominator is (J/kg) * (kg/m3) = J/m3, therefore RHS quotient is J/s/m2 / J/m3 = m/s.
		SublimationSnow_Potential_m_p_s = Ground_LE_W_p_m2 / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);

		//Control for negative values
		SublimationSnow_Potential_m_p_s = Inputs::forceZeroIfLessThanOrAlmostEqualZero(SublimationSnow_Potential_m_p_s, Epsilon_Tolerance_1E_negative15);

		//SublimationSnow_Potential_m (m) is product of SublimationSnow_Potential_m_p_s (m/s) and SimulationTimeStep_Duration_sec[timeStep]
		beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_Potential_m") = SublimationSnow_Potential_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
	}
}

//OpenWaterEvaporationCalc::calculateGroundSublimationStatistical function is called as calculateGroundSublimation functor to compute SublimationSnow_Potential_m (m) for StatisticalHydro
void OpenWaterEvaporationCalc::calculateGroundSublimationStatistical(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//SublimationSnow_Potential_m (m) is input->PotentialSublimation_SnowOnGround_m[timeStep], potential sublimation depth (m) from input time series
	beC->by_key(MapPixel_ID, DataFolder_ID, "SublimationSnow_Potential_m") = input->PotentialSublimation_SnowOnGround_m[timeStep];
}
