#include "Flux_LateralGroundWater.h"
#include "../Evaporation/RootZoneEvapoTranspirationCalc.h"
#include "../Drainage/DrainageToEvapotranspirationZone.h"
#include "../Drainage/DrainageToSaturationZone.h"
#include "../Runoff/RunoffSurfaceSaturationExcess.h"

//Note: Catchment average soil moisture deficit (AveSMD) used to obtain available storage in each Topographic Index bin, StorageDeficit_VadoseZone_TI_m
//Note: AveSMD ExponentialDecay (timeStep=0) AveSMD = -mM * log(Q0_mpdt / SZQ) and SZQ = T0_m2pdt * exp(-TI_average); 
//Note: AveSMD PowerDecay (timeStep=0) AveSMD = mM * (1.0 - SZQ) and SZQ = pow(Q0_mpdt / T0_m2pdt, 1 / nN)*(TI_average);
//Note: StorageDeficit_VadoseZone_TI_m ExponentialDecay StorageDeficit_VadoseZone_TI_m = AveSMD + mM * (TI_Value_Average - TI_Value[ia]) 
//Note: StorageDeficit_VadoseZone_TI_m PowerDecay StorageDeficit_VadoseZone_TI_m = AveSMD + mM * pow(Q0_mpdt/T0_m2pdt, 1/nN)*(TI_Value_Average - TI_Value[ia]);

//Note: Theory of updating AveSMD (m) with drainage from vadose zone and runoff from saturated zone based on Wolock (1993).
//Note: Soil volumetric water has 3 zones ranging from wilting point, 0wp to field capacity, 0fc, to saturation, 0sat
//Note: 1) evapotranspiration dynamically contains 0wp to 0sat; 0sat-0fc (gravitational water) releases to lower vadose or saturated zone
//Note: 2) vadose zone dynamically contains 0sat to 0fc; 0sat-0fc drains to saturated zone; it is ground surface to water table
//Note: 3) saturated zone generates subsurface runoff; it can rise into vadose and evapotranspiration zones to fill gravitational voids
//Note: The evapotranspiration and vadose zones are conceptually situated alongside each other, each filling a different volumetric component
//Note: The evapotranspiration zone extends to the limit of evapotranspiration, only processing capillary held water
//Note: The vadose zone is represented by average soil moisture deficit, AveSMD, StorageDeficit_VadoseZone_m, or StorageDeficit_VadoseZone_TI_m
//Note: AveSMD and StorageDeficit_VadoseZone_m are catchment averages, while StorageDeficit_VadoseZone_TI_m varies with topographic index
//Note: For ts=0, AveSMD (m) defined in TerrainProcessor::getAveSoilMoistureDeficit using Q0_mph, T0_m2ph, mM, nN with Exponential or Power Decay
//Note: For ts>0, StorageDeficit_VadoseZone_m & StorageDeficit_VadoseZone_TI_m defined in Flux_LateralGroundWater::calculate as above
//Note: StorageDeficit_FieldCapacity_SoilEvapZone_Max_m (m) = Evapotranspiration_Depth_m * (Soil_FieldCapacity_m3pm3 - Soil_WiltingPoint_m3pm3)
//Note: StorageDeficit_VadoseZone_m (m) = (Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3) * Depth_SoilColumn_Ave_m
//Reference: Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124). 

//Note: Catchment average soil moisture deficit (StorageDeficit_VadoseZone_m, AveSMD) determines water storage in each Topographic Index bin
//Note: Saturation excess runoff is an outcome of AveSMD and local SMD values, based on variable source area and topographic index (TI) theory
//Note: The catchment Bulk Area uses local SMD values to determine saturation excess in each TI bin, StorageDeficit_VadoseZone_TI_m
//Note: The Green Infrastructure (GI) Area uses AveSMD via Depth_to_Groundwater_Table_m = StorageDeficit_VadoseZone_m / Porosity_Drainable_m3pm3
//Note: Green Infrastructure (GI) uses Depth_to_Groundwater_Table_m in Groundwater_StormwaterDevice.cpp to determine saturation by groundwater

//calculate function computes local soil moisture deficit for each topographic index (TI) bin
void Flux_LateralGroundWater::calculate(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//n is coefficient for power function decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
	//Figure 1 of Wang et al. (2006) shows how Transmissivity varies with n ranging from 0.25 to 16
	double nN = input->InputXml["Parameter_n_KsatPowerDecay"];
	//Q0_mpdt (m/ts) is initial subsurface discharge from catchment, adjusted by Ratio_Hour_to_Second * SimulationTimeStep_Duration_sec[timeStep]
	double Q0_mpdt = input->InputXml["Discharge_Subsurface_Initial_mph"] * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];
	//T0_m2pdt transmissivity (m^2 / timestep) when AveSMD = 0 (watershed fully saturated) adjusted by Ratio_Hour_to_Second * SimulationTimeStep_Duration_sec[timeStep]
	double T0_m2pdt = input->InputXml["VadoseZone_Transmissivity_Max_m2ph"] * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];
	//mM (m) is scaling parameter describing exponential decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
	//Beven, Lamb, et al. (1995) give a physical interpretation of the decay parameter m (m) is that it controls the effective depth, z (m) of the catchment soil profile. This it does interactively in Eq T = T0_m2pdt * exp (-f*z), where f = (Theta_sat - Theta_fc)/m. See Figure 18.1 for values of m ranging from 0.02 to 0.05 creating effective depths z ranging from 1 to 2 m.
	//Beven, K., Lamb, R., Quinn, P., Romanowics, R., & Freer, J. (1995). TOPMODEL. In V. P. Singh (Ed.), Computer models of watershed hydrology (pp. 627-688). Colorado: Water Resources Publications.
	double mM = input->InputXml["Parameter_m_KsatExpDecay"];

	//AveSMD (m) is average soil moisture deficit, used to update local SMD values; = (Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3) * Depth_SoilColumn_Ave_m
	//StorageDeficit_VadoseZone_TS_Prior_m is constant in time, used by all folders while For Loop through DataDrawers iteratively updates AveSMD in StorageDeficit_VadoseZone
	//Note:	StorageDeficit_VadoseZone_TS_Prior_m (m) avoids use of iteratively changing input->RepoDict["StorageDeficit_VadoseZone_m"], updated in StorageDeficit_VadoseZone.cpp for each DataDrawer
	double StorageDeficit_VadoseZone_m = input->InputXml["StorageDeficit_VadoseZone_TS_Prior_m"];
	//Runoff_SatExcess_m(m) is the saturation excess runoff, reset to 0.0
	double Runoff_SatExcess_m = 0.0;
	//Porosity_drainable_m3pm3 is the difference between Soil_SaturationPoint_m3pm3 and Soil_FieldCapacity_m3pm3, HydroPlusConfig.xml inputs
	double Porosity_drainable_m3pm3;

	//Flux_SatExGW_to_SoilFieldCapacity_m (m) reset to zero at start of TI_Value 
	//Note: Flux_SatExGW_to_SoilFieldCapacity_m is groundwater discharge filling soil voids to field capacity in the evapotranspiration zone
	beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_SatExGW_to_SoilFieldCapacity_m") = 0.0;
	//Groundwater_surficial_m (m) reset to zero at start of TI_Value
	//Note: Groundwater_surficial_m (m) used to track surficial expression of variable source area
	beC->by_key(MapPixel_ID, DataFolder_ID, "Groundwater_surficial_m") = 0;
	//Groundwater_surficial_frac (m) reset to zero at start of TI_Value
	//Note: Groundwater_surficial_frac (frac) used to track surficial expression of variable source area
	beC->by_key(MapPixel_ID, DataFolder_ID, "Groundwater_surficial_frac") = 0;
	//Drainage_SoilEvapZone_m (m) reset to zero at start of TI_Value
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_m") = 0;
	//EvapoTranspiration_SoilEvapZone_m (m) reset to zero at start of TI_Value
	beC->by_key(MapPixel_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") = 0.0;
	//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) reset to zero at start of TI_Value
	beC->by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = 0.0;
	//SoilEvapZoneStorageTemp_m (m) reset to zero at start of TI_Value
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SoilEvapZone_m") = 0.0;
	//Drainage_VadoseZone_m (m) initialized reset to zero at start of TI_Value
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_m") = 0.0;
	//Runoff_SatExcess_m (m) is total saturation excess runoff across all ia bins
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_m") = 0.0;
	//Storage_VadoseZone_m (m) reset to zero at start of TI_Value
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_m") = 0.0;

	//Before entering TI loop with potential for removal of Runoff_SatExcess_TI_m removed, record initial values of fluxes entering soil
	//Infiltration_viaSatEx_PreRunoffSatExcess_m (m) set to Infiltration_viaSatEx_m (m) before entering TI loop, where latter term is recomputed based on any Runoff_SatExcess_TI_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_PreRunoffSatExcess_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m");
	//Infiltration_viaInfilEx_PreRunoffSatExcess_m (m) set to Infiltration_viaSatEx_m (m) before entering TI loop, where latter term is recomputed based on any Runoff_SatExcess_TI_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_PreRunoffSatExcess_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m");
	//Drainage_macroPore_PreRunoffSatExcess_m (m) set to Infiltration_viaSatEx_m (m) before entering TI loop, where latter term is recomputed based on any Runoff_SatExcess_TI_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_PreRunoffSatExcess_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_m");
	beC->by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_PreRunoffSatExcess_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_m");
	//Drainage_SoilEvapZone_PreRunoffSatExcess_m (m) set to 0 at start of TI loop, then summed as product of Drainage_SoilEvapZone_TI_m (m) and TI_binArea_frac (fraction),prior to any Runoff_SatExcess_TI_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_PreRunoffSatExcess_m") = 0;
	//Drainage_SoilEvapZone_MacroPore_frac (fraction) initialized to 0.5, will be recomputed during TI loop if needed, when computing Runoff_SatExcess_TI_m in DrainageToSaturationZone 
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_MacroPore_frac") = 0.5;
	//Drainage_SoilEvapZone_Infil_frac (fraction) initialized to 1 - Drainage_SoilEvapZone_MacroPore_frac, will be recomputed in DrainageToSaturationZone as needed
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_Infil_frac") = 1 - beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_MacroPore_frac");
	//Reset to zero the infiltration and drainage to soil scalars, so they can be recomputed as sum of TI bin vector components that account for any Runoff_SatExcess_TI_m
	//Infiltration_viaSatEx_m (m) set to zero before TI bin loop, and during TI bin loop will be recomputed in DrainageToSaturationZone.cpp
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m") = 0;
	//Infiltration_viaInfilEx_m (m) set to zero before TI bin loop, and during TI bin loop will be recomputed in DrainageToSaturationZone.cpp
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m") = 0;
	//Drainage_macroPore_m (m) set to zero before TI bin loop, and during TI bin loop will be recomputed in DrainageToSaturationZone.cpp
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_m") = 0;

	//IndexLoopSize = input->Size_TI_Value_Vector of topographic index (TI) for statistical model, which is 1 less than Size_TI_Area_Vector, which has all TI values in histogram 
	int IndexLoopSize = int(input->Size_TI_Value_Vector);

	//if input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro" then IndexLoopSize = 1 the IndexLoopSize = 1, and input->Size_TI_Value_Vector > input->Size_TI_Area_Vector
	if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") { 
		//IndexLoopSize set for just one topographic index (TI) value in each folder in SpatialTemperatureHydro
		IndexLoopSize = 1; 
		//TI_Area_frac vector set to 1 for first two bins, to ensure beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac") = 1
		if (timeStep == 0) {
			//If TI_Area_frac.size() <= IndexLoopSize then enter and add the IndexLoop size to vector
			if (input->TI_Area_frac.size() <= IndexLoopSize) {
				//TI_Area_frac[1] = 1 for SpatialTemperatureHydro model, in order to compute TI_binArea_frac = 1
				input->TI_Area_frac.push_back(1);
			}
		}
	}

	//Statistically distributed topographic index algorithm taken from TOPMODEL code developed by Keith Beven and used by USGS D.M. Wolock https://pubs.er.usgs.gov/publication/wri934124
	//Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124).
	//Note: The algorithm requires the topographic index (TI) bins (ia) are sorted from highest (wettest) to lowest TI values
	//Note: The algorithm requires when ia = 0, the TI bin TI_Value[0] should be the upper limit and its TI_Area_frac[0] = 0
	//Note: The algorithm presumes for each ia, the area of catchment between discrete TI values TI_Value[ia] and TI_Value[ia+1] is TI_Area_frac[ia]
	//Note: The algorithm requires for each ia, the area of catchment bounding TI_Value[ia] is = to (TI_Area_frac[ia] + TI_Area_frac[ia+1])/2 

	
	//Initialize vectors at start of simulation 
	if (timeStep == 0) {
		//StorageDeficit_VadoseZone_TI_m (m) initialized to 0 for start of simulation for each TI bin
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m", IndexLoopSize, 0.0);
		//StorageDeficit_VadoseZone_TI_temp_m (m) initialized to 0 for start of simulation for each TI bin
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_temp_m", IndexLoopSize, 0.0);
		//Drainage_SoilEvapZone_TI_m (m) initialized to 0 for start of simulation for each TI bin
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m", IndexLoopSize, 0.0);
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m should mbe initialized here with 0 first, then filled with StorageDeficit_FieldCapacity_SoilEvapZone_TS_first_m 
		// ...within the simulation for each drawer at the beninning of DrainageToEvapotranspirationZone::calculate
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m", IndexLoopSize, 0.0);
		//Infiltration_TI_m is infiltration via InfilEx for TI bin, initialized to 0
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "Infiltration_TI_m", IndexLoopSize, 0.0);
		//Drainage_macropore_TI_m is drainage macropore for TI bin, initialized to 0
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "Drainage_macropore_TI_m", IndexLoopSize, 0.0);
		//Discharge_Groundwater_TI_m(m) is groundwater discharge filling the soil voids to field capacity in the evapotranspiration zone; initialized to 0
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "Flux_SatExGW_to_SoilFieldCapacity_TI_m", IndexLoopSize, 0.0);
		//Storage_VadoseZone_TI_m (m) is unsaturated zone storage, initialized to 0; represents only gravitational water
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "Storage_VadoseZone_TI_m", IndexLoopSize, 0.0);
		//Runoff_SatExcess_TI_m (m) is saturation excess runoff, initialized to 0
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_TI_m", IndexLoopSize, 0.0);
		//Drainage_VadoseZone_TI_m (m) is unsaturated zone runoff, initialized to 0; represents release of gravitational water
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "Drainage_VadoseZone_TI_m", IndexLoopSize, 0.0);
		//EvapoTranspiration_SoilEvapZone_TI_m (m) is soil evapotranspiration, initialized to 0
		beC->vec_ensure_length(MapPixel_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_TI_m", IndexLoopSize, 0.0);
	}

	//Loop through all discrete topographic index (TI) values, sorted from largest to smallest, with fractional area for each TI = to average of bounding bin areas
	//Note: For spatially distributed model, the folder area has a single TI value
	for (int ia = 0; ia < IndexLoopSize; ia++) {
	
		//TI_bin = ia as current TI bin to get appropriate topographic index properties
		int TI_bin = ia;
		//ia is the TI bin set one place later, drier, in the sorted list than bin ia; on first pass iaa = 1
		int iaa = ia + 1;

		//SpatialTemperatureHydro the TI_binArea_frac = 1, with 1 assigned to both TI_Area_frac[ia] and TI_Area_frac[iaa]
		//TI_binArea_frac is catchment area (fraction) bounding each discrete TI value
		beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac") = 0.5 * (input->TI_Area_frac[ia] + input->TI_Area_frac[iaa]);

		//Original TOPMODEL weights discrete TI catchment area fractions of histogram with bounding TI parameter values
		//Note: HydroPlus weights discrete TI parameter values, i.e, SoilParameter[ia], with TI_binArea_frac, catchment area bounding discrete value

		//Compute local soil moisture deficit, which will determine which unsaturated zones are saturated
		//Power decay computation of StorageDeficit_VadoseZone_TI_m
		if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") {

			//if Model = "StatisticalHydro" BulkArea, then Size_TI_Area_Vector > 1 and StorageDeficit_VadoseZone_TI_m[ia] (m) function of histogram TI value
			if (input->SimulationStringParams["Model_Selection"] == "StatisticalHydro" || input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {

				//StorageDeficit_VadoseZone_TI_m (m) computed as Power Decay from AveSMD, mM, Q0_mpdt, T0_m2pdt, nN, and TI values based material between Eq 2 and 4 in Wang et al. (2005), using histogram bins
				//Note: the Eq below combines from Wang et al. (2005) S_bar and S_i, and S_i (S = SMD) is using topographic index, not soil TI, and power (1/n), not (1/-n), and (a/tanB)^(1/n) should not be power term
				beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] = StorageDeficit_VadoseZone_m + mM * pow(Q0_mpdt / T0_m2pdt, 1 / nN) * (input->TI_Value_Average - input->TI_Value[ia]);
			}
			//if Model = "SpatialTemperatureHydro", then Size_TI_Area_Vector = 1 and StorageDeficit_VadoseZone_TI_m[ia] (m) function of map TI value
			if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") {
				//StorageDeficit_VadoseZone_TI_m (m) computed as Power Decay from AveSMD, mM, Q0_mpdt, T0_m2pdt, nN, and TI values based material between Eq 2 and 4 in Wang et al. (2005), using map pixels
				//Note:  the Eq below combines from Wang et al. (2005) S_bar and S_i, and S_i (S = SMD) is using topographic index, not soil TI, and power (1/n), not (1/-n), and (a/tanB)^(1/n) should not be power term
				beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] = StorageDeficit_VadoseZone_m + mM * pow(Q0_mpdt / T0_m2pdt, 1 / nN) * (input->TI_Value_Average - input->TI_Value[MapPixel_ID]);
			}

			//Theory of power function decay requires m <= S (AveSMD or StorageDeficit_VadoseZone_TI_m) based on constraint governing decay with depth of transmissivity, Tz = T0_m2pdt*(1-S/m)^n  
			//StorageDeficit_VadoseZone_TI_m (m) is soil moisture deficit between top of soil and groundwater table, defined as Theta_sat-Theta_wp or Theta_sat-Theta_fc (based on Beven and Wood, 1983).
			//Beven and Wood (1983). Catchment Geomorphology and the Dynamics of Runoff Contributing Areas, Journal of Hydrology. 65:139-158. In Special Issue: Scale Problems in Hydrology, Editors: Rodriguez-Iturbe and V.K. Gupta.
			//Theory explained after Eq 5 of Wang, J., Endreny, T. A., & Hassett, J. M. (2006). Power function decay of hydraulic conductivity for a TOPMODEL-based infiltration routine Hydrological Processes, 20(18), 3825-3834. 
			//Wang et al (2006) after Eq 5. explains in power decay conditions the m parameter sets the upper limit for the soil moisture deficit S (m), due to Eq 7, T = T0_m2pdt(1-S/m)^n. If n = 2, then S > m. If n <=1, S>m will generate negative T profiles., if n >=2, S>m will generate non-monotonic T profiles.
			//Figure 1 from Wang et al. (2006) illustrates the decay of transmissivity, Tz, with depth as a function of n, varying from 0.25 to 16
			if (beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] > mM)
				beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] = mM;
		}
		//Exponential decay computation of StorageDeficit_VadoseZone_TI_m
		else {

			//if Model = "StatisticalHydro" BulkArea, then Size_TI_Area_Vector > 1 and StorageDeficit_VadoseZone_TI_m[ia] (m) function of histogram TI value
			if (input->SimulationStringParams["Model_Selection"] == "StatisticalHydro" || input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {

				//StorageDeficit_VadoseZone_TI_m (m) computed as Exponential Decay from AveSMD, mM, and TI values based on Eq 31 from Wolock (1993), using histogram bins
				beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] = StorageDeficit_VadoseZone_m + mM * (input->TI_Value_Average - input->TI_Value[ia]);

			}
			//if Model= "SpatialTemperatureHydro", then Size_TI_Area_Vector = 1 and StorageDeficit_VadoseZone_TI_m[ia] (m) function of map TI value
			if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") {

				//StorageDeficit_VadoseZone_TI_m (m) computed as Exponential Decay from AveSMD, mM, and TI values based on Eq 31 from Wolock (1993), using map pixels
				beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] = StorageDeficit_VadoseZone_m + mM * (input->TI_Value_Average - input->TI_Value[MapPixel_ID]);
			}
		}

		//StorageDeficit_VadoseZone_TI_temp_m set to save temporary record of StorageDeficit_VadoseZone_TI_m for use in RunoffSurfaceSaturationExcess.cpp
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_temp_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia];

		//if StorageDeficit_VadoseZone_TI_m < 0, reset to 0 to avoid erroneous calculation in DrainageToSaturationZone.cpp
		if (beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] < 0.0) {
			//StorageDeficit_VadoseZone_TI_m reset to 0 after generating return flow
			beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_m")[ia] = 0.0;
			//Topmodel code: IF(SD(IA).LT.0.)SD(IA) = 0.
		}
		
		//Compute soil water balance, separating capillarity water for ET and gravitational water for saturation excess or subsurface flow
		//call DrainageToEvapotranspirationZone to add infiltrated water, and release gravitational and macropore water to unsaturated zone in TI bins
		DrainageToEvapotranspirationZone::calculate(input, beC, TI_bin, MapPixel_ID, DataFolder_ID, timeStep);
		//call RootZoneEvapoTranspirationCalc to remove water held by capillarity into evapotranspiration demand in TI bins
		RootZoneEvapoTranspirationCalc::Calculate(input, beC, TI_bin, MapPixel_ID, DataFolder_ID, timeStep);
		//call DrainageToSaturationZone to update gravitational water storage, saturation excess runoff, and drainage to groundwater in TI bins
		DrainageToSaturationZone::calculate(input, beC, TI_bin, MapPixel_ID, DataFolder_ID, timeStep);
		//Note: Consider Refactor to check if Runoff_SatExcess_TI_m * TI_binArea_frac fits within available pervious depression storage
		//call RunoffSurfaceSaturationExcess to sum the saturation excess runoff across all TI bins
		RunoffSurfaceSaturationExcess::calculate(input, beC, TI_bin, MapPixel_ID, DataFolder_ID, timeStep);
	}
	

	//Compute variables at end of TI loop
	//Infiltration_m (m) = Infiltration_viaSatEx_m + Infiltration_viaInfilEx_m, where RHS terms had saturation excess runoff removed	
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m");
	
	//Storage_SoilEvapZone_m (m) is computed as difference of maximum storage deficit StorageDeficit_FieldCapacity_SoilEvapZone_Max_m (m) and actual storage deficit StorageDeficit_FieldCapacity_SoilEvapZone_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_SoilEvapZone_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_Max_m") - beC->by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m");

	//If GI is operational in HydroPlusConfig.xml then implement the transfer of water between BulkArea and GI
	//Note: Determine the saturation excess runoff (m) for the BulkArea, for use as runon in the GI area
	if (input->RepoDict["Flag_GI_Simulated"] == 1 && beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "BulkArea") {

		//Groundwater_surficial_BulkArea_frac (frac) is defined with BulkArea folder Groundwater_surficial_frac (frac)
		//Note: Groundwater_surficial_frac (frac) signifies fraction of topographic index with saturation excess runoff 
		//Note: Groundwater_surficial_frac (frac) can be used in Green Infrastructure routines to represent saturation excess runoff
		//Note: beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_m") might be updated for GI this way
		input->RepoDict["Groundwater_surficial_frac"] = beC->by_key(MapPixel_ID, DataFolder_ID, "Groundwater_surficial_frac");

		//RepoDict is used for transfer between BulkArea and GI area folders
		//Area_BulkArea_m2 (m2) obtained for use to GI
		input->RepoDict["Area_BulkArea_m2"] = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");

		//The subsequent sequence of steps will 1) adjust depths to extract runoff going from BulkArea to GI and 2) adjust depths to move between areas and conserve water
		//1. BulkAreaPerviousDrainingToGI_frac used to extract fraction of pervious area depth going from BulkArea to GI
		double BulkAreaPerviousDrainingToGI_frac = input->RepoDict["BulkAreaPervious_DrainingToGI_frac"];

		//Runoff_SatEx_BAtoGI_m (m) is pervious area depth leaving BulkArea and going to GI via saturation excess runoff
		input->RepoDict["Runoff_SatEx_BAtoGI_m"] = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_m") * BulkAreaPerviousDrainingToGI_frac;
		//Runoff_InfilEx_BAtoGI_m (m) is pervious area depth leaving BulkArea and going to GI via infiltration excess runoff
		input->RepoDict["Runoff_InfilEx_BAtoGI_m"] = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_InfilExcess_m") * BulkAreaPerviousDrainingToGI_frac;

		//Runoff_SatExcess_m (m) is pervious area depth remaining in BulkArea after removing Runoff_SatEx_BAtoGI_m (m) going to GI 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_m") - input->RepoDict["Runoff_SatEx_BAtoGI_m"];
		//Runoff_InfilExcess_m (m) is pervious area depth remaining in BulkArea after removing Runoff_InfilEx_BAtoGI_m (m) going to GI 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_InfilExcess_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_InfilExcess_m") - input->RepoDict["Runoff_InfilEx_BAtoGI_m"];

		//2. Runoff_SatEx_BAtoGI_m (m) is adjusted to folder area depth by multiplication with PerviousCover_frac, with fraction comparable to TotalPerviousCover_m2/TotalArea_m2
		input->RepoDict["Runoff_SatEx_BAtoGI_m"] = input->RepoDict["Runoff_SatEx_BAtoGI_m"] * 
			beC->by_key(MapPixel_ID, DataFolder_ID, "PerviousCover_frac");
		//2. Runoff_InfilEx_BAtoGI_m (m) is adjusted to folder area depth by multiplication with PerviousCover_frac, with fraction comparable to TotalPerviousCover_m2/TotalArea_m2
		input->RepoDict["Runoff_InfilEx_BAtoGI_m"] = input->RepoDict["Runoff_InfilEx_BAtoGI_m"] * 
			beC->by_key(MapPixel_ID, DataFolder_ID, "PerviousCover_frac");

		//RepoDict is used for transfer between BulkArea and GI area folders
		//Runoff_SatExcess_BAtoGI_m3 (m^3) is volume of runon to GI, which is conserved by continuity; derived from folder depth Runoff_SatEx_BAtoGI_m and folder area Area_BulkArea_m2
		input->RepoDict["Runoff_SatExcess_BAtoGI_m3"] = input->RepoDict["Runoff_SatEx_BAtoGI_m"] * 
			input->RepoDict["Area_BulkArea_m2"];
		//Runoff_InfilExcess_BAtoGI_m3 (m^3) is volume of runon to GI, which is conserved by continuity; derived from folder depth Runoff_InfilEx_BAtoGI_m and folder area Area_BulkArea_m2
		input->RepoDict["Runoff_InfilExcess_BAtoGI_m3"] = input->RepoDict["Runoff_InfilEx_BAtoGI_m"] * 
			input->RepoDict["Area_BulkArea_m2"];
	}


	//Drainage_SoilEvapZone_prior_m (m) updates cumulative infiltration in the next time step; it is combined with EvapoTranspiration_SoilEvapZone_prior_m (m). 
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_prior_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_m");
	//EvapoTranspiration_SoilEvapZone_prior_m (m) updates cumulative infiltration in the next time step; it is combined with Drainage_SoilEvapZone_prior_m (m). 
	beC->by_key(MapPixel_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_prior_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m");
}
