#include "ShortVegEvaporationCalc.h"

//References:
//McCutcheon, S.C., Martin, J. and T.O. Barnwell Jr. (1993). Water Quality. In D. R. Maidment (Ed.), Handbook of Hydrology (pp. 11.1-11.73). New York: McGraw-Hill.

void ShortVegEvaporationCalc::calculateShortVegEvaporationTemperature(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 ShortVegEvaporationCalc::calculateShortVegEvaporationStatistical to estimate evaporation
		ShortVegEvaporationCalc::calculateShortVegEvaporationStatistical(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		//return is C++ function to leave the .cpp
		return;
	}

	//Initialize LAI_BAI_SVeg_m2_p_m2 
	double LAI_BAI_SVeg_m2_p_m2 = input->LAI_BAI_SVeg_m2_p_m2[timeStep];
	//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"));

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

	//If Storage_PerviousPondedWater_m > Constant_1E_negative6 then presume SVegLEE_W_p_m2 needs to be reduced due to some used to evaporate depression storage
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_PerviousPondedWater_m") > Constant_1E_negative6) {
		//EvaporationSublimation_SVegCanopy_m_p_s (m/s) derived from SVegLEE_W_p_m2 (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3) ...
		//Note: If depression storage received SVegLEE_W_p_m2 / LAI_BAI_SVeg_m2_p_m2, then SVegLEE_W_p_m2 - (SVegLEE_W_p_m2 / LAI_BAI_SVeg_m2_p_m2) remains
		//Note: OpenWaterEvaporationCalc::calculateImpDepEvapTemperature computes Evaporation_PerviousPondedWater_m_p_s using ...
		//Note: ... SVegLE_W_p_m2 / LAI_BAI_SVeg_m2_p_m2 as energy allowed to evaporate ponded water under the short veg canopy
		//Note: ... The SVegLET_W_p_m2 - (SVegLET_W_p_m2 / LAI_BAI_SVeg_m2_p_m2) remains available for evapotranspiration
		//Note: ... The SVegLEE_W_p_m2 - (SVegLEE_W_p_m2 / LAI_BAI_SVeg_m2_p_m2) remains available for canopy evaporation
		//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.
		EvaporationSublimation_SVegCanopy_m_p_s = (beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLEE_W_p_m2") - (beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLEE_W_p_m2") / LAI_BAI_SVeg_m2_p_m2)) / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);
	}
	//Else Storage_PerviousPondedWater_m <= Constant_1E_negative6 and use full amount of SVegLEE_W_p_m2 
	else {
		//EvaporationSublimation_SVegCanopy_m_p_s (m/s) derived from SVegLEE_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.
		EvaporationSublimation_SVegCanopy_m_p_s = beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLEE_W_p_m2") / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);
	}
	//Control for negative values
	if (EvaporationSublimation_SVegCanopy_m_p_s <= 0) { EvaporationSublimation_SVegCanopy_m_p_s = 0; }

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

	//If EvaporationSublimation_SVegCanopy_m (m) is greater than Storage_Rain_SVegCanopy_m, then set to Storage_Rain_SVegCanopy_m
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "EvaporationSublimation_SVegCanopy_m") > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m")) {
		beC->by_key(MapPixel_ID, DataFolder_ID, "EvaporationSublimation_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m");
	}

	//Storage_Rain_SVegCanopy_m (m) is updated, defined as Storage_Rain_SVegCanopy_m - EvaporationSublimation_SVegCanopy_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "EvaporationSublimation_SVegCanopy_m");
}

void ShortVegEvaporationCalc::calculateShortVegEvaporationStatistical(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	double SublimationSnow_SVegCanopy_m = 0.0;        // Initialize SublimationSnow_SVegCanopy_m (amount of evaporated snow) = 0 
	double EvaporationRain_SVegCanopy_m = 0.0;           // Initialize EvaporationRain_SVegCanopy_m (amount of evaporated water) = 0
	double EvaporationSublimation_SVegCanopy_m = 0.0;                // Initialize ShortVegEP (potential evaporation) = 0

	// SublimationSnow_SVegCanopy_m (the amount of evaporated snow) = PEShortVegSnow (the potential evaporation of snow) 
	SublimationSnow_SVegCanopy_m = input->PotentialSublimation_SnowOnTree_m[timeStep];

	//routine to evaporate liquid water after evaporating any snow water equivalent
	//enter conditional if potential evaporation for snow water exceeds stored snow water
	if (SublimationSnow_SVegCanopy_m > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m")) {

		//calcluate fraction of potential evaporation utilzed by snow water evaporation
		double Sublimation_to_EvaporationSublimation_frac = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") / SublimationSnow_SVegCanopy_m;

		//evaporate all of the snow water, which may be 0 when ShortVegStoredSWE is 0
		SublimationSnow_SVegCanopy_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m");

		//If Storage_SVegCanopy_Max_m > 0 then enter for division
		if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m") > 0) {
			//EvaporationRain_TreeCanopy_m based on Eq 6 of Wang et al. (2008), using Deardorff (1978) and Noilhan and Planton(1989)
			//Note: Regulation of evaporation by 2/3 function is used in Liang et al. (1994) VIC paper
			EvaporationRain_SVegCanopy_m = (1 - Sublimation_to_EvaporationSublimation_frac) * input->PotentialEvaporation_WaterOnTree_m[timeStep] *
				pow(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") / beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_SVegCanopy_Max_m"), (2.0 / 3.0));
		}
		//Else remove Storage_SVegCanopy_Max_m from denominator
		else {
			//EvaporationRain_TreeCanopy_m based on Eq 6 of Wang et al. (2008), using Deardorff (1978) and Noilhan and Planton(1989)
			EvaporationRain_SVegCanopy_m = (1 - Sublimation_to_EvaporationSublimation_frac) * input->PotentialEvaporation_WaterOnTree_m[timeStep] *
				pow(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m"), (2.0 / 3.0));
		}
		//If EvaporationRain_SVegCanopy_m > Storage_Rain_SVegCanopy_m then limit to Storage_Rain_SVegCanopy_m
		if (EvaporationRain_SVegCanopy_m > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m")) {
			EvaporationRain_SVegCanopy_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m");
		}
		//reduce the shortVeg stored water by the tree evaporation of water
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_SVegCanopy_m") - EvaporationRain_SVegCanopy_m;
	}
	//reduce the shortVeg stored snow by the shortVeg evaporation of snow
	//if Storage_Snow_SVegCanopy_m was 0, then SublimationSnow_SVegCanopy_m was set to 0 in the above condition
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_SVegCanopy_m") - SublimationSnow_SVegCanopy_m;

	EvaporationSublimation_SVegCanopy_m = SublimationSnow_SVegCanopy_m + EvaporationRain_SVegCanopy_m;

	beC->by_key(DataDrawer_ID, DataFolder_ID, "EvaporationSublimation_SVegCanopy_m") = EvaporationSublimation_SVegCanopy_m;
	beC->by_key(DataDrawer_ID, DataFolder_ID, "EvaporationRain_SVegCanopy_m") = EvaporationRain_SVegCanopy_m;
	beC->by_key(DataDrawer_ID, DataFolder_ID, "SublimationSnow_SVegCanopy_m") = SublimationSnow_SVegCanopy_m;

}




