#include "RootZoneEvapoTranspirationCalc.h"

//Root zone evapotranspiration only removes water between field capacity and wilting point, referred to as capillarity held water
//Note: Consider refactor to allow soil areas lower soil moisture limit than wilting point, such as theta residual

//Note: Root zone refers to soil below all pervious areas and any impervious area with tree cover, not only the soil land cover type.
//Note: Soil below any impervious area with tree cover is considered recharged via the lateral redistribution of groundwater. 

//Note: The model StatisticalHydro uses the input data of potential ET to determine atmospheric demand
//Note:	The model SpatialTemperatureHydro uses latent energy to determine atmospheric demand

//Note: EvapoTranspiration_SoilEvapZone_m (m) is reduced from potential ET by ratio of available soil moisture based on theory, including ...
//Note: ... ET_act = ET_pot * (soil storage)/(soil storage max) based on Saxton (1986) Eq in Fisher et al. (2005)
//Note: ... Fisher et al (2005) show potential ET based on Penmam-Monteith is too high relative to observed actual ET. 
//Note: ... The Eqs of Fisher et al. (2005) are not numbered, the same LHS variable is used for multiple Eqs of potential ET
//Note: ... Fisher et al (2005) give the adjustment in text near citation of Saxton (1986), as f(phi); improvement illustrated in Fig 2 vs Fig 1. 
//References:
//Fisher, J. B., DeBiase, T. A., Qi, Y., Xu, M., & Goldstein, A. H. (2005). Evapotranspiration models compared on a Sierra Nevada forest ecosystem. Environmental Modelling & Software, 20(6), 783-796. doi:https://doi.org/10.1016/j.envsoft.2004.04.009

//RootZoneEvapoTranspirationCalc::Calculate function handles StatisticalHydro and SpatialTemperatureHydro models
void RootZoneEvapoTranspirationCalc::Calculate(Inputs* input, CompactRagged* beC, int TI_bin, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//ia = TI_bin as current TI bin to get appropriate topographic index properties
	int ia = TI_bin;

	//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; }

	//EvapoTranspiration_SoilEvapZone_m (m) reset to zero at start of each TI value
	double EvapoTranspiration_SoilEvapZone_m = 0;
	double Storage_EvapoTransZone_m = 0;
	//StorageDeficit_FieldCapacity_SoilEvapZone_Max_m (m) is the maximum soil deficit for capillarity held water for ET; = Evapotranspiration_Depth_m * (Soil_FieldCapacity_m3pm3 - Soil_WiltingPoint_m3pm3)
	double StorageDeficit_FieldCapacity_SoilEvapZone_Max_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_Max_m");
	
	//If Model is StatisticalHydro or Model is SpatialTemperatureHydro And timeStep equals 1, then call statistical calculation based on weather input data, given the spatial calculation has not yet been called
	if (input->SimulationStringParams["Model_Selection"] == "StatisticalHydro" || ((input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro" || input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") && timeStep == 0)) {
		//EvapoTranspiration_SoilEvapZone_m (m) = potential ET read in from Evaporation.csv input file
		//Note: PotentialEvapotranspiration_TreeCover_m (m) computed for tree, yet used for tree and pervious areas 
		//Note: ... This is a model simplification for the 1st timestep, while in the next time steps cover types are represented
		//Note: Depths of evapotranspiration are relative to the area of analysis, and do not get scaled by fraction of cover type
		EvapoTranspiration_SoilEvapZone_m = input->PotentialEvapotranspiration_TreeCover_m[timeStep];
	}
	//Else If input->SimulationStringParams["Model_Selection"] is "SpatialTemperatureHydro" or "SpatialBufferwGI" and timeStep > 0
	else {
		//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"));

		//EvapoTranspiration_SoilEvapZone_m_p_s (m/s) initialized
		double EvapoTranspiration_SoilEvapZone_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(DataDrawer_ID, DataFolder_ID, "Storage_PerviousPondedWater_m") > Constant_1E_negative6 || 
			beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_ImperviousPondedWater_m") > Constant_1E_negative6) {
			//EvapoTranspiration_SoilEvapZone_m_p_s (m/s) derived as latent energy (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3)
			//Note: If depression storage received TreeLET_W_p_m2 / LAI_BAI_Tree_m2_p_m2, then TreeLET_W_p_m2 - (TreeLET_W_p_m2 / LAI_BAI_Tree_m2_p_m2) remains
			//Note: EvapoTranspiration_SoilEvapZone_m is weighted by PC_m2 and presumes any trees over impervious are evapotranspiring water from pervious
			//Note: Latent energy is (SoilLE_W_p_m2 * SoilCover_noTreeCanopy_frac + (SVegLET_W_p_m2 - SVegLET_W_p_m2 / LAI_BAI_SVeg_m2_p_m2) * ShortVegCover_noTreeCanopy_frac + (TreeLET_W_p_m2 - TreeLET_W_p_m2 / LAI_BAI_Tree_m2_p_m2) * TreeCover_frac), and ...
			//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: ... 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: TreeCover_frac vs TreeCanopyCover_overPervious_frac is used given tree cover over pervious & impervious areas have roots within soil layer 
			//Note: Soils covered by impervious may receive recharge from groundwater via lateral redistribution.
			//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.
			EvapoTranspiration_SoilEvapZone_m_p_s = (beC->by_key(DataDrawer_ID, DataFolder_ID, "SoilLE_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "SoilCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "ImpLE_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "PermeablePavementCover_noTreeCanopy_frac") + (beC->by_key(DataDrawer_ID, DataFolder_ID, "SVegLET_W_p_m2") - input->SafeDivide(beC->by_key(DataDrawer_ID, DataFolder_ID, "SVegLET_W_p_m2"), LAI_BAI_SVeg_m2_p_m2)) * beC->by_key(DataDrawer_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") + (beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeLET_W_p_m2") - input->SafeDivide(beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeLET_W_p_m2"), LAI_BAI_Tree_m2_p_m2)) * beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeCover_frac")) / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);
		}
		//Else Storage_PerviousPondedWater_m <= Constant_1E_negative6 and use full amount of TreeLEE_W_p_m2 
		else {
			//EvapoTranspiration_SoilEvapZone_m_p_s (m/s) derived as latent energy (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3)
			//Note: If depression storage received TreeLET_W_p_m2 / LAI_BAI_Tree_m2_p_m2, then TreeLET_W_p_m2 - (TreeLET_W_p_m2 / LAI_BAI_Tree_m2_p_m2) remains
			//Note: EvapoTranspiration_SoilEvapZone_m is weighted by PC_m2 and presumes any trees over impervious are evapotranspiring water from pervious
			//Note: Latent energy is (SoilLE_W_p_m2 * SoilCover_noTreeCanopy_frac + SVegLET_W_p_m2 * ShortVegCover_noTreeCanopy_frac + ...
			//Note: ... TreeLET_W_p_m2 * TreeCover_frac), and ...
			//Note: TreeCover_frac is used given tree cover over pervious and impervious areas have roots within the soil layer 
			//Note: Soils covered by impervious may receive recharge from groundwater via lateral redistribution.
			//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.
			EvapoTranspiration_SoilEvapZone_m_p_s = (beC->by_key(DataDrawer_ID, DataFolder_ID, "SoilLE_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "SoilCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "ImpLE_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "PermeablePavementCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "SVegLET_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeLET_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeCover_frac")) / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);
		}

		//EvapoTranspiration_SoilEvapZone_m_p_s (m/s) derived as latent energy (W/m2) divided by (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3)
		//Note: EvapoTranspiration_SoilEvapZone_m is weighted by PC_m2 and presumes any trees over impervious are evapotranspiring water from pervious
		//Note: Latent energy is (SoilLE_W_p_m2 * SoilCover_noTreeCanopy_frac + SVegLET_W_p_m2 * ShortVegCover_noTreeCanopy_frac + ...
		//Note: ... TreeLET_W_p_m2 + TreeCover_frac), and ...
		//Note: TreeCover_frac is used given tree cover over pervious and impervious areas have roots within the soil layer 
		//Note: Soils covered by impervious may receive recharge from groundwater via lateral redistribution.
		//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.
		EvapoTranspiration_SoilEvapZone_m_p_s = (beC->by_key(DataDrawer_ID, DataFolder_ID, "SoilLE_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "SoilCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "ImpLE_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "PermeablePavementCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "SVegLET_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeLET_W_p_m2") * beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeCover_frac")) / (LatentHeat_Vaporization_JpkgK * DensityWater_kgpm3);

		//EvapoTranspiration_SoilEvapZone_m_p_s (m/s) next divided by pervious area + tree over impervious so it is scaled to cover area
		EvapoTranspiration_SoilEvapZone_m_p_s = input->SafeDivide(EvapoTranspiration_SoilEvapZone_m_p_s, (beC->by_key(DataDrawer_ID, DataFolder_ID, "SoilCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "PermeablePavementCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "ShortVegCover_noTreeCanopy_frac") + beC->by_key(DataDrawer_ID, DataFolder_ID, "TreeCover_frac")));

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

		//EvapoTranspiration_SoilEvapZone_m (m) is product of EvapoTranspiration_SoilEvapZone_m_p_s (m/s) and SimulationTimeStep_Duration_sec[timeStep]
		EvapoTranspiration_SoilEvapZone_m = EvapoTranspiration_SoilEvapZone_m_p_s * input->SimulationTimeStep_Duration_sec[timeStep];
	}

	//If Type equals GI do not proceed with updates to topographic index (TI) terms
	//Note: GI or non-BulkArea folder does not use TI_bin value that is passed in function
	//Note: Consider Refactor to implement soil water availability constraints of modified water retention curve, as for BulkArea
	if (beC->by_key_str(DataDrawer_ID, DataFolder_ID, "Type") != "BulkArea") {
		//Storage_EvapoTransZone_m (m) = GI water depth available for ET, defined in Inflow_StormwaterDevice.cpp as beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_EvapoTransZone_m")
		Storage_EvapoTransZone_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m");
		//EvapoTranspiration_SoilEvapZone_m (m) equals EvapoTranspiration_SoilEvapZone_m, potential ET 
		beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") = EvapoTranspiration_SoilEvapZone_m;

		//If the EvapoTranspiration_SoilEvapZone_m (m) is greater than Storage_EvapoTransZone_m (m), then reset to available soil evapotranspiration zone water
		if (beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") > Storage_EvapoTransZone_m) {
			//EvapoTranspiration_SoilEvapZone_TI_m (m) is limited to Storage_EvapoTransZone_m (m), available soil evapotranspiration zone water
			beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") = Storage_EvapoTransZone_m;
		}
	}

	//if Type equals BulkArea then compute Storage_EvapoTransZone_m as function of StorageDeficit_FieldCapacity_SoilEvapZone_TI_m
	if (beC->by_key_str(DataDrawer_ID, DataFolder_ID, "Type") == "BulkArea") {
		//EvapoTranspiration_SoilEvapZone_TI_m (m) reset to zero at start of each TI_Value vector
		beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia] = 0.0;

		//if EvapoTranspiration_SoilEvapZone_m > 0 then remove water from SoilEvapZone
		if (EvapoTranspiration_SoilEvapZone_m > 0.0) {
			//If StorageDeficit_FieldCapacity_SoilEvapZone_Max_m greater than zero, then enter for division by StorageDeficit_FieldCapacity_SoilEvapZone_Max_m
			if (beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_Max_m") > 0) {
				//EvapoTranspiration_SoilEvapZone_TI_m (m) is EvapoTranspiration_SoilEvapZone_m * (1 - (StorageDeficit_FieldCapacity_SoilEvapZone_TI_m / StorageDeficit_FieldCapacity_SoilEvapZone_Max_m))
				//Note: Equal to EvapoTranspiration_SoilEvapZone_m when StorageDeficit_FieldCapacity_SoilEvapZone_TI_m = 0, no capacity remains
				beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia] = EvapoTranspiration_SoilEvapZone_m * (1 - (beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] / beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_Max_m")));
				//Storage_EvapoTransZone_m (m) = StorageDeficit_FieldCapacity_SoilEvapZone_Max_m - StorageDeficit_FieldCapacity_SoilEvapZone_TI_m; 
				//	Example: StorageDeficit_FieldCapacity_SoilEvapZone_Max_m = 0.1 m, StorageDeficit_FieldCapacity_SoilEvapZone_TI_m = 0.08 m, then Storage_EvapoTransZone_m = 0.1 - 0.08 = 0.02 m
				Storage_EvapoTransZone_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_Max_m") - beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia];
			}
			else {
				beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia] = 0;
				//Storage_EvapoTransZone_m (m) equals zero
				Storage_EvapoTransZone_m = 0;
			}

			//If the EvapoTranspiration_SoilEvapZone_TI_m (m) > Storage_EvapoTransZone_m (m), then reset to available soil evapotranspiration zone water
			if (Inputs::isGreaterThan(beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia], Storage_EvapoTransZone_m)) {
				//EvapoTranspiration_SoilEvapZone_TI_m (m) is limited to Storage_EvapoTransZone_m (m), available soil evapotranspiration zone water
				beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia] = Storage_EvapoTransZone_m;
			}
		}

		//StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) is increased by removal of EvapoTranspiration_SoilEvapZone_TI_m
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m decreases due to Infiltration_TI_m in DrainageToEvapotranspirationZone.cpp
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m is a deficit, aka capacity; all TI bins have same thickness
		beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] = beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] + beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia];

		//EvapoTranspiration_SoilEvapZone_m (m) = sum EvapoTranspiration_SoilEvapZone_TI_m for each TI bin for time step, weighted by watershed fraction
		//	HydroPlus weights discrete TI parameter values, i.e, SoilParameter[ia], with TI_binArea_frac, catchment area bounding discrete value
		beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") + beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia] * beC->by_key(DataDrawer_ID, DataFolder_ID, "TI_binArea_frac");

		//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) = sum StorageDeficit_FieldCapacity_SoilEvapZone_TI_m 
		//Note: Summing for each TI bin for time step, weighted by watershed fraction, TI_binArea_frac
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_m ranges frm 0=saturated, 1=dry; available storage up to field capacity
		//HydroPlus weights discrete TI parameter values, i.e, SoilParameter[ia], with TI_binArea_frac, catchment area bounding discrete value
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") + beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] * beC->by_key(DataDrawer_ID, DataFolder_ID, "TI_binArea_frac");
	}
	//Ending of BulkArea calculations
}

//Exploring potential refactor with this function; it has never been active for any HydroPlus test cases
//RootZoneEvapoTranspirationCalc::SoilWaterRetentionCurve function handles StatisticalHydro and SpatialTemperatureHydro models
//Theory: Soil water retention curves represent the relationship between the amount of water in the soil and the soil's matric potential
//Theory: Multi-parameter water retention curves were created by Brooks and Corey (1964), van Genuchten (1980), and Saxton and Rawls (2006)
//Note: HydroPlus uses soil water retention curve theory to represent how soil columns dry from top to bottom under evapotranspiration
//Note: Given HydroPlus does not simulate matric potential and uses a uniform vadose zone above the saturated zone, theoretical curves do not fit
//Note: Instead HydroPlus uses a modified water retention curve that is linear to approximate the capillary water above the wilting point
//Note: Consider refactoring to use pedotransfer functions and theoretical curves when resources become available for HydroPlus
//References:
//Brooks, R. H., & Corey, A. T. (1964). Hydraulic Properties of Porous Media Affecting Fluid Flow. Journal of Irrigation and Drainage Engineering, ASCE, 92(IR2), 61-68. 
//Saxton, K. E., & Rawls, W. J. (2006). Soil Water Characteristic Estimates by Texture and Organic Matter for Hydrologic Solutions. Soil Science Society of America Journal, 70(5), 1569-1578. 
//van Genuchten, M. T. (1980). A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal, 44(5), 892-898. 
	/*
	* Illustration of water retention curve (WRC) under two moisture conditions, theta saturation fraction (Theta_Sat_frac) > 0.5 and <= 0.5
	* WRC is anchored to bottom right but moves along top for Theta_sat_frac > 0.5, within a unit range of 1
	* WRC is anchored to top left but moves along bottom for Theta_sat_frac <= 0.5, within a unit range of 1
	* Vertices A (top left), B (top right), C (bottom right), and D (bottom left) bound the WRC within a rectangle in x-z coordinate space
	* The x-axis is relative soil saturation level, Theta_Sat_frac and goes from 0 (left) to 1 
	* The z-axis is effective depth of evapotranspiration zone, Depth_ETZ_m, and goes from 0 (top) to Depth_ETZ_max_m
	* Point E represents the WRC tip as it passes through the x-axis for the TI bin or pixel, at top (Theta_Sat_frac > 0.5) or bottom (Theta_Sat_frac <= 0.5) 
	* Point F represents the depth of the actual Depth_ETZ_m for the TI bin or pixel
	* Point G represents the WRC segment as it passes through the Depth_ETZ_m for the TI bin or pixel
	* 
	* Define Terms:
	* 
	* Theta_Sat_frac = 1 - (beC->by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_Max_m") / beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia]
	* 
	* Depth_to_Groundwater_Table_TI_m = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] / input->InputXml["Porosity_Drainable_m3pm3"];
	* 
	* Depth_ETZ_max_m = maximum(EvapotranspirationDepth_TreeCover_m, EvapotranspirationDepth_SVegCover_m, EvapotranspirationDepth_SoilCover_m, 1 m)
	* Depth_ETZ_max_m = minimum(Depth_to_Groundwater_Table_TI_m, Depth_ETZ_max_m); places Line DE above saturation
	* Depth_ETZ_max_m = maximum(beC->by_key(MapPixel_ID, DataFolder_ID, "Evapotranspiration_Depth_m"), Depth_ETZ_max_m); ensures Point F < Point D
	* 
	* Note: Last Depth_ETZ_max_m adjustment avoids resetting user inputs used to compute Evapotranspiration_Depth_m, and ...
	* Note: 1) Evapotranspiration_Depth_m & StorageDeficit_FieldCapacity_SoilEvapZone_TI_m can extend beneath the groundwater table, ...
	* Note: ... as it represents capillary pore storage, while StorageDeficit_VadoseZone_TI_m represents gravitational pore storage
	* Note: 2) If Evapotranspiration_Depth_m & StorageDeficit_FieldCapacity_SoilEvapZone_TI_m extend beneath the groundwater table, ...
	*
	* Find Point E(x,z) via scaling to length 0 to 1 from 0.5 to 1 (Theta_Sat_frac > 0.5) or from 0 to 0.5 (Theta_Sat_frac <= 0.5): 
	* If Theta_Sat_frac > 0.5, then E(x) = [1-2*(1-Theta_Sat_frac)], E(z) = 0
	* If Theta_Sat_frac <= 0.5, then E(x) = 2*Theta_Sat_frac, E(z) = Depth_ETZ_max_m
	* 
	* Find Point G(x,z) via interpolation between E and C (Theta_Sat_frac > 0.5) or E can A (Theta_Sat_frac <= 0.5):
	* If Theta_Sat_frac > 0.5, then G(x) = C(x) - G(z)*[(C(x) - E(x)) / (C(z) - E(z))], G(z) = Depth_ETZ_m
	* If Theta_Sat_frac <= 0.5, then G(x) = A(x) + G(z)*[(A(x) - E(x)) / (A(z) - E(z))], G(z) = Depth_ETZ_m
	* 
	* Find unit Area of evapotranspiration zone available water for Theta_Sat_frac > 0.5 or Theta_Sat_frac <= 0.5
	* If Theta_Sat_frac > 0.5, then trapezoid AEGF or WRC_ETZ_m = 0.5*[E(x) + G(x)]*F(z), where F(z) = Depth_ETZ_m
	* If Theta_Sat_frac <= 0.5, then triangle AGF or WRC_ETZ_m = 0.5*G(x)*F(z), where F(z) = Depth_ETZ_m
	* 
	* Scale the unit Area of evapotranspiration zone available water, WRC_ETZ_m, based on ...
	* ... unit Area of WRC maximum capillary water capacity within rectangle ABCD, as WRC_max_m = B(x)*C(z)
	* Creating WRC scaling factor for evapotranspiration zone water, Theta_WRC_ETZ_scalingFactor_frac = WRC_ETZ_m / WRC_max_m
	* 
	* Constrain Storage_EvapoTransZone_m water using the WRC scaling factor:
	* beC->vec_by_key(MapPixel_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia] = Theta_WRC_ETZ_scalingFactor_frac * beC->vec_by_key(MapPixel_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m")[ia]
	*
	* 
	* Theta_Sat_frac > 0.5			  Theta_Sat_frac <= 0.5
	A----------\E-----------B		A-----------------------B	Depth = 0
	|		 <--\-->		|		| \						|
	|			 \			|		|  \					|
	|			  \			|		|	\					|
	F==============\G=======|		F====\G=================|	Depth = Depth_ETZ_m
	|				\		|		|	  \					|
	|				 \		|		|	   \				|
	|				  \		|		|		\				|
	|				   \	|		|		 \				|
	|					\	|		|		  \				|
	|					 \	|		|		   \			|
	|					  \	|		|		 <--\-->		|
	D-----------------------C		D------------\E---------C	Depth = Depth_ETZ_max_m
		
	*/
