#include "TreeEvaporationCalc.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 TreeEvaporationCalc::calculateTreeEvaporationTemperature(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 TreeEvaporationCalc::calculateTreeEvaporationStatistical to estimate evaporation
		TreeEvaporationCalc::calculateTreeEvaporationStatistical(input, beC, MapPixel_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(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_TreeCanopy_m_p_s (m/s) initialized
	double EvaporationSublimation_TreeCanopy_m_p_s = 0;

	//If Storage_PerviousPondedWater_m or Storage_ImperviousPondedWater_m > Constant_1E_negative6 then presume TreeLEE_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 || beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_ImperviousPondedWater_m") > Constant_1E_negative6) {
		//EvaporationSublimation_SVegCanopy_m_p_s (m/s) derived from TreeLEE_W_p_m2 (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3) ...
		//Note: If depression storage received TreeLEE_W_p_m2 / LAI_BAI_Tree_m2_p_m2, then TreeLEE_W_p_m2 - (TreeLEE_W_p_m2 / LAI_BAI_Tree_m2_p_m2) remains
		//Note: OpenWaterEvaporationCalc::calculateImpDepEvapTemperature computes Evaporation_PerviousPondedWater_m_p_s using ...
		//Note: ... TreeLE_W_p_m2 / LAI_BAI_Tree_m2_p_m2 as energy allowed to evaporate ponded water under the tree canopy
		//Note: ... The TreeLET_W_p_m2 - (TreeLET_W_p_m2 / LAI_BAI_Tree_m2_p_m2) remains available for evapotranspiration 
		//Note: ... The TreeLEE_W_p_m2 - (TreeLEE_W_p_m2 / LAI_BAI_Tree_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_TreeCanopy_m_p_s = (beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLEE_W_p_m2") - input->SafeDivide(beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLEE_W_p_m2"), LAI_BAI_Tree_m2_p_m2)) / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);
	}
	//Else Storage_PerviousPondedWater_m <= Constant_1E_negative6 and use full amount of TreeLEE_W_p_m2 
	else {
		//EvaporationSublimation_TreeCanopy_m_p_s (m/s) derived from TreeLEE_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_TreeCanopy_m_p_s = beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLEE_W_p_m2") / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);
	}
	//Control for negative values
	EvaporationSublimation_TreeCanopy_m_p_s = Inputs::forceZeroIfLessThanOrAlmostEqualZero(EvaporationSublimation_TreeCanopy_m_p_s, Epsilon_Tolerance_1E_negative15);

	//EvaporationSublimation_TreeCanopy_m (m) is product of EvaporationSublimation_TreeCanopy_m_p_s (m/s) and SimulationTimeStep_Duration_sec[timeStep]
	beC->by_key(MapPixel_ID, DataFolder_ID, "EvaporationSublimation_TreeCanopy_m") = EvaporationSublimation_TreeCanopy_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_TreeCanopy_m") > beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m")) {
		beC->by_key(MapPixel_ID, DataFolder_ID, "EvaporationSublimation_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m");
	}

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

void TreeEvaporationCalc::calculateTreeEvaporationStatistical(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	double SublimationSnow_TreeCanopy_m = 0.0;        // Initialize TreeEPSnow (amount of evaporated snow) = 0 
	double EvaporationRain_TreeCanopy_m = 0.0;           // Initialize TreeEPW (amount of evaporated water) = 0
	double EvaporationSublimation_TreeCanopy_m = 0.0;                // Initialize TreeEP (potential evaporation) = 0

	// TreeEPSnow (the amount of evaporated snow) = PotentialSublimation_SnowOnTree_m (the potential evaporation of snow) 
	SublimationSnow_TreeCanopy_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_TreeCanopy_m > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_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_TreeCanopy_m") / SublimationSnow_TreeCanopy_m;
		SublimationSnow_TreeCanopy_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m");

		//If Storage_TreeCanopy_Max_m > 0 then enter for division
		if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m") > 0) {
			//EvaporationRain_TreeCanopy_m based on Eq 6 of Wang et al. (2008), using Deardorff (1978) and Noilhan and Planton(1989)
			//Note: Evaporate all of the snow water, which may be 0 when Storage_Snow_TreeCanopy_m is 0
			//Note: Regulation of evaporation by 2/3 function is used in Liang et al. (1994) VIC paper
			EvaporationRain_TreeCanopy_m = (1.0 - Sublimation_to_EvaporationSublimation_frac) * input->PotentialEvaporation_WaterOnTree_m[timeStep] * pow(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") / beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_TreeCanopy_Max_m"), (2.0 / 3.0));
		}
		//Else remove Storage_TreeCanopy_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_TreeCanopy_m = (1.0 - Sublimation_to_EvaporationSublimation_frac) * input->PotentialEvaporation_WaterOnTree_m[timeStep] * pow(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m"), (2.0 / 3.0));
		}

		//If EvaporationRain_TreeCanopy_m > Storage_Rain_TreeCanopy_m then limit to Storage_Rain_TreeCanopy_m
		if (EvaporationRain_TreeCanopy_m > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m")) {
			EvaporationRain_TreeCanopy_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m");
		}
		//reduce the tree stored water by the tree evaporation of water
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") - EvaporationRain_TreeCanopy_m;
	}

	//reduce the tree stored snow by the tree evaporation of snow
	//if Storage_Snow_TreeCanopy_m was 0, then SublimationSnow_TreeCanopy_m was set to 0 in the above condition
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m") - SublimationSnow_TreeCanopy_m;

	EvaporationSublimation_TreeCanopy_m = SublimationSnow_TreeCanopy_m + EvaporationRain_TreeCanopy_m;

	beC->by_key(DataDrawer_ID, DataFolder_ID, "EvaporationSublimation_TreeCanopy_m") = EvaporationSublimation_TreeCanopy_m;
	beC->by_key(DataDrawer_ID, DataFolder_ID, "EvaporationRain_TreeCanopy_m") = EvaporationRain_TreeCanopy_m;
	beC->by_key(DataDrawer_ID, DataFolder_ID, "SublimationSnow_TreeCanopy_m") = SublimationSnow_TreeCanopy_m;
}



