 #include "InfiltrationExponentialDecay.h"
#include <iostream>

//Partitions Flux_to_Infiltration_Drainage_m (m) into Infiltration_viaInfilEx_m (m) and Runoff_InfilExcess_m (m)
//Note: Infiltration_viaInfilEx_m is potential and may be reduced if saturated areas exists
//Note: HydroPlusConfig.xml DataFolder values of K0_mph are combined with input->InputXml["Parameter_n_KsatPowerDecay"] & ["Parameter_m_KsatExpDecay"]
//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: Green Ampt infiltration routines are based on USGS TOPMODEL Fortran Code from D. Wolock (USGS) and K. Bevin, "C  TOPMODEL DEMONSTRATION PROGRAM VERSION 95.02"
//Note: Theory from Wang et al. (2005). and from Eq. 2 to 8 in Beven (1984)
//Note:	Before ponding, potential infiltration rate > rainfall rate
//Note: Potential infiltration rate declines from maximum when cumulative infiltration = 0 toward minimum capacity, defined as K0
//Note: Rainfall rate has to exceed potential infiltration rate for ponding to occur

//Note: Cumulative infiltration determines performance of Green Ampt function.
//Note: Wang et al. (2008) code used temp = total infiltrated so far (CumInfilQ) - 0.1 * (BaseFlow + SumRET), SumRET = SoilEvapZone Evapotranspiration
//Note: Yang et al. (2011) code used cumulative infiltration = initial SoilEvapZone deficit  - current (previous ts) soil evapotranspiration zone deficit

//References: 
//Beven (1984). Infiltration into a class of vertically non-uniform soils. Hydrological Sciences Journal 29: 425–434.
//Beven and Wood (1983). Catchment Geomorphology and the Dynamics of Runoff Contributing Areas, Journal of Hydrology. 65:139-158. 
//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, J., Endreny, T. A., & Hassett, J. M. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
//Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124). 

void InfiltrationExponentialDecay::calculate(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//DWSMI (m3/m3) is Green Ampt infiltration parameter for delta soil moisture across wetting front, safely set to (Soil_SaturationPoint_m3pm3 - Soil_WiltingPoint_m3pm3)
	beC->by_key(MapPixel_ID, DataFolder_ID, "DWSMI") = beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_SaturationPoint_m3pm3") - beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_FieldCapacity_m3pm3");

	//if Type is not PermeablePavement then HydraulicConductivity_Infiltration_mpdt is function of K0_mph
	if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") != "PermeablePavement") {
		//HydraulicConductivity_Infiltration_mpdt (m/timeStep) is HydroPlusConfig input of hydraulic conductivity K0_mph (m/hr) times Ratio_Hour_to_Second and SimulationTimeStep_Duration_sec[timeStep]
		beC->by_key(MapPixel_ID, DataFolder_ID, "HydraulicConductivity_Infiltration_mpdt") = beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Ksat_mph") * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];
	}
	//else if Type is PermeablePavement then HydraulicConductivity_Infiltration_mpdt is function of Pavement_HydraulicConductivity_mph
	//Note: GI PermeablePavement sits above soil layer, and hence is the layer for this infiltration routine
	else {
		//HydraulicConductivity_Infiltration_mpdt (m/timeStep) is Pavement_HydraulicConductivity_Decayed_mph times Ratio_Hour_to_Second and SimulationTimeStep_Duration_sec[timeStep]
		//Note: Pavement_HydraulicConductivity_Decayed_mph is reduced in Inflow_StormwaterDevice::GI_RateOfPavementPermeability from HydroPlusConfig.xml element Pavement_HydraulicConductivity_mph
		beC->by_key(MapPixel_ID, DataFolder_ID, "HydraulicConductivity_Infiltration_mpdt") = beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_HydraulicConductivity_Decayed_mph") * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];
	}

	//Flux_to_Infiltration_Drainage_m (m) is the sum of 3 terms, all adjusted previously to pervious area depth from folder area depth
	//Note: 3 terms: Flux_to_PerviousArea_Rain_SnowMelt_Irrigation_ImpRunon_m (m), Flux_to_GI_PerviousArea_BARunon_m (m), and Ponding_on_PerviousArea_PriorTS_m (m)
	beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_Infiltration_Drainage_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_PerviousArea_Rain_SnowMelt_Irrigation_ImpRunon_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_GI_PerviousArea_BARunon_m") +
		beC->by_key(MapPixel_ID, DataFolder_ID, "Ponding_on_PerviousArea_PriorTS_m");
	
	//macroPore_frac (fraction) is macropore fraction (0 to 1) within unsaturated zone soils 
	//Note: macroPore_frac > 0 allows a fraction of surface water to bypass infiltration and SoilEvapZone evapotranspiration, and directly enter groundwater
	double macroPore_frac = beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Macropore_frac");
	//Drainage_macroPore_m (m) is the water from Flux_to_Infiltration_Drainage_m that bypasses infiltration and soil evapotranspiration zone, going directly to subsurface flow
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_Infiltration_Drainage_m") * macroPore_frac;
	//Flux_to_Infiltration_m (m) is the water from Flux_to_Infiltration_Drainage_m that bypasses macropores, going directly to soil evapotranspiration zone
	beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_Infiltration_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_Infiltration_Drainage_m") * (1 - macroPore_frac);
	//Infiltration_viaSatEx_m (m) is Flux_to_Infiltration_m reduced by (1-InfiltExcessGovernedArea_frac), regulated only later by saturation excess and not by Green Ampt infiltration
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_Infiltration_m") * (1 - beC->by_key(MapPixel_ID, DataFolder_ID, "InfiltExcessGovernedArea_frac"));
	//Flux_to_Infiltration_viaInfilEx_m (m) is Flux_to_Infiltration_m reduced by InfiltExcessGovernedArea_frac, regulated by Green Ampt infiltration and later by saturation excess
	double Flux_to_Infiltration_viaInfilEx_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_Infiltration_m") * beC->by_key(MapPixel_ID, DataFolder_ID, "InfiltExcessGovernedArea_frac");
	//Delta_timeStep is simulation time interval, such as an 1 hour timestep or a 5 min timestep; determined by timestep of Weather.csv inputs
	double Delta_timeStep = 1.0;

	//If InfiltExcessGovernedArea_frac equals 1 then define remaining infiltration terms and return
	//Note: This allows HydroPlusConfig.xml to have K0_mph = 0 for some GI Types, yet utilize Infiltration_viaSatEx_m to convey water into storage
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "InfiltExcessGovernedArea_frac") == 0) {
		beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m") = 0;
		beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m");
		return;
	}

	//consider resetting cumulative infiltration if not first timestep of simulation
	if (timeStep > 0) {
		//Outflow_SoilEvapZone_prior_m (m) represents water leaving soil evapotranspiration zone by drainage and evapotranspiration
		double Outflow_SoilEvapZone_prior_m = (beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_prior_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_prior_m"));

		//Infiltration_viaInfilEx_cumulative_m (m) is reduced by the amount of water that left soil evapotranspiration zone
		beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") - Outflow_SoilEvapZone_prior_m;

		//if Infiltration_viaInfilEx_cumulative_m (m) < 0 then set to 0
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") < 0.0) {
			beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") = 0.0;
		}
	}

	//Infiltration decision tree with multiple if-else options
	//Note: The inPonding function is never called from the calculate function, and is always called by the prepPonding function, which is called by two functions: 
	//Note: a) the estimateTP function calls the prepPonding function when entering wet weather; 
	//Note: b) the beforePonding function calls the prepPonding function in each subsequent time step of wet weather.
	//Option 1: Soil treated as impervious area when hydraulic conductivity = 0; HydraulicConductivity_Infiltration_mpdt (m/timeStep)
	if (beC->by_key(MapPixel_ID, DataFolder_ID, "HydraulicConductivity_Infiltration_mpdt") == 0.0) {
		//infilQ (m) defined as zero when HydraulicConductivity_Infiltration_mpdt equals 0
		beC->by_key(MapPixel_ID, DataFolder_ID, "infilQ") = 0;
		//Flux_InfilEx_to_Ponding_m (m) defined as Flux_to_Infiltration_viaInfilEx_m (m), removing zero infilQ
		beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_InfilEx_to_Ponding_m") = Flux_to_Infiltration_viaInfilEx_m - beC->by_key(MapPixel_ID, DataFolder_ID, "infilQ");

		//Infiltration excess runoff (m) is runoff to infiltration that arrives at a rate greater than the infiltration rate
		//Infiltration_viaInfilEx_cumulative_m (m) is cumulative infiltration, set to 0 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") = 0.0;
	}
	//Option 2: No infiltration occurring; this option for no ponding is also used later when Flux_to_Infiltration_viaInfilEx_m > 0
	//Note: this may need to be modified to allow for intentional delayed infiltration of ponded water 
	else if (Flux_to_Infiltration_viaInfilEx_m == 0.0) {
		//noPonding routine uses Flux_to_Infiltration_viaInfilEx_m = Flux_to_Infiltration_m, and will record 0 infiltration in time series
		noPonding(Flux_to_Infiltration_viaInfilEx_m, input, beC, MapPixel_ID, DataFolder_ID, timeStep);
	}
	//Option 3: Infiltration starting anew, with cumulative infiltration = 0 and the need to estimate time to ponding
	else if (beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") == 0.0) {
		//estimateTP routine will estimate time to ponding based on Flux_to_Infiltration_m and timestep and soil properties in routine
		estimateTP(timeStep, Flux_to_Infiltration_viaInfilEx_m, Delta_timeStep, input, beC, MapPixel_ID, DataFolder_ID, timeStep);
	}
	//Option 4: Infiltration continues without ponding condition
	else {
		//beforePonding routine	uses Flux_to_Infiltration_m and timestep to determine infiltration and check if ponding will occur
		beforePonding(timeStep, Flux_to_Infiltration_viaInfilEx_m, Delta_timeStep, input, beC, MapPixel_ID, DataFolder_ID, timeStep);
	}

	//Note: Consider refactor to remove ponding (flag), set to false (0) to avoid chance for anomalous behavior
	//beC->by_key(DataDrawer_ID, DataFolder_ID, "ponding") = 0;

	//Note: Consider refactor and rename infilQ to Infiltration_viaInfilEx_m to avoid maintaining two variables for same value
	//Infiltration_viaInfilEx_m (m) is renamed infilQ (m), which is fraction of Flux_to_Infiltration_viaInfilEx_m (m) that infiltrates, determined above by calling Green-Ampt routines below
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "infilQ");

	//Flux_InfilEx_to_Ponding_m (m) is Flux_to_Infiltration_viaInfilEx_m (m) - infilQ (m) sent to PerviousDepressionStorageCalc.cpp for ponding or runoff
	beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_InfilEx_to_Ponding_m") = Flux_to_Infiltration_viaInfilEx_m - beC->by_key(MapPixel_ID, DataFolder_ID, "infilQ");

	//Note: calcInfilEXQ function turned off to allow for ponding
	//Note: Consider refactor to place flag in HydroPlusConfig.xml to allow or disallow for infiltration of ponding
	//calcInfilEXQ(Flux_to_Infiltration_viaInfilEx_m, folder, input, beC, DataDrawer_ID, DataFolder_ID, timeStep);

	//Infiltrated water will be reduced by the depth of saturation excess runoff in DrainageToSaturationZone.cpp
	//satExcessInflow_m (m) is Flux_to_Infiltration_m reduced by macroPore_frac and 1-InfiltExcessGovernedArea_frac, regulated only later by saturation excess and not by Green Ampt infiltration
	double satExcessInflow_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_m");
	//infilExcessInflow_m (m) is Flux_to_Infiltration_m reduced by macroPore_frac and InfiltExcessGovernedArea_frac, regulated by Green Ampt infiltration and later by saturation excess
	double infilExcessInflow_m = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_m");
	//Infiltration_m (m) is sum of satExcessInflow_m and infilExcessInflow_m, weighted by InfiltExcessGovernedArea_frac in ::calculate function
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_m") = satExcessInflow_m + infilExcessInflow_m;
}

//Routine from USGS Topmodel Fortran Code from Wolock (USGS)/Bevin, subroutine titled "THERE IS NO PONDING IN THIS TIME STEP"
//Option: all Flux_to_Infiltration_m infiltrates, no ponding during infiltration; npR is Flux_to_Infiltration_m for no ponding situation
void InfiltrationExponentialDecay::noPonding(double npR, Inputs *input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//Set flag ponding = 0, which is false and indicates no ponding
	//beC->by_key(DataDrawer_ID, DataFolder_ID, "ponding") = 0;
	//infilQ (m) is depth of water infiltrated; if no ponding, this equals all Flux_to_Infiltration_m
	beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") = npR;
	//      TOPMODEL, DF = RINT*DT

	//Prohibit negative values for infilQ (m)
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") < 0.0) {
		beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") = 0.0;
	}
	if (!beC->has_var("Infiltration_viaInfilEx_cumulative_m")) { 
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") = 0.0;
	}
	//Update cumulative infiltration (m) depth
	double temp = beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m");
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") = temp;
}

//Note: calcInfilEXQ function turned off to allow for ponding
//Note: Consider refactor to place flag in HydroPlusConfig.xml to allow or disallow for infiltration of ponding
//Note: Routine from USGS Topmodel Fortran Code from Wolock (USGS)/Bevin, subroutine titled "C  INFILTRATION EXCESS CALCULATIONS"
//Note: Option: infiltration excess runoff; ifR is Flux_to_Infiltration_m for infiltration excess runoff situation
void InfiltrationExponentialDecay::calcInfilEXQ(double ifR, Inputs *input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//Runoff_InfilExcess_m (m) initiated to = 0 for pervious areas 
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_InfilExcess_m") = 0.0;
	//Runoff_InfilExcess_m (m) is residual of Flux_to_Infiltration_m (m) - infilQ (m), the water that infiltrated to SoilEvapZone; 
	//Note: No ponded water is retained on site, all is routed to overland runoff
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_InfilExcess_m") = ifR - beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ");
	//Note: TOPMODEL, REX = P - DF 

	//Prohibit negative values for Runoff_InfilExcess_m (m)
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_InfilExcess_m") < 0.0) {
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_InfilExcess_m") = 0.0;
	}
}

//Routine from USGS Topmodel Fortran Code from Wolock (USGS)/Bevin, subroutine titled "INITIAL ESTIMATE OF TIME TO PONDING & PONDING STARTS AT BEGINNING OF TIME STEP"
//Option: before ponding during infiltration; biT is timestep number, bR is Flux_to_Infiltration_m for no ponding situation, bDT is timestep duration
void InfiltrationExponentialDecay::beforePonding(int bIT,  double bR, double bDT, Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//K0 hydraulic conductivity (m / timestep) was adjusted to per time step with input->SimulationTimeStep_Duration_sec[timeStep] 
	double K0 = beC->by_key(DataDrawer_ID, DataFolder_ID, "HydraulicConductivity_Infiltration_mpdt");
	//      TOPMODEL, XKF = XK0; define hydraulic conductivity

	//mM is scaling parameter used to adjust infiltration rates from Config inputs; TOPMODEL theory; Beven et al. (1995a, 1995b)
	//m (m) is scaling parameter describing exponential decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
	//Beven, Lamb, et al. (1995a) 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 * 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. (1995b). 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"];
	//fF = 1/mM, scaling parameter for infiltration rates, taken directly from TOPMODEL code SZF = 1./SZM
	//fF could be thetaSoilMoisture/mMin, given in text below Eq 10 of Wang et al. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
	double fF = 1 / mM;
	//      TOPMODEL, SZF = 1./SZM
	//CCSF (m) coefficient is product of WFS_m (m), wetting front suction, and Green Ampt term DWSMI (delta soil moisture), Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3
	double CSSF = beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_WettingFront_Suction_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "DWSMI");
	//      TOPMODEL, CD=HF*DTH
	//bPP is Flux_to_Infiltration_m Rate (m / timestep) 
	double bPP = bR / bDT;
	//Note:     TOPMODEL, RINT = P/DT

	//F1 (m) is cumulative infiltration estimate at start 
	double F1 = 0.0;
	F1 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m");
	//      TOPMODEL, F1=CUMF

	//RR (m / timestep) value of rainfall rate when it equals potential infiltration rate at start of ponding, based on Green Ampt theory 
	//Theory from Eq. 2 to 8 in Beven (1984). Infiltration into a class of vertically non-uniform soils. Hydrological Sciences Journal 29: 425–434.
	//Note: 	Before ponding, potential infiltration rate > rainfall rate
	//Note: 	Potential infiltration rate declines from maximum when cumulative infiltration = 0 toward minimum capacity, defined as K0
	//Note: 	Rainfall rate has to exceed potential infiltration rate for ponding to occur
	double RR = 0.0;
	//RR is modified version of Eq 6 and 7 from Beven (1984) Infiltration paper; lifted directly from Beven TOPMODEL code. 
	//RR is represented in Eq 10 of Wang et al. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
	RR = - K0 * fF * (CSSF + F1) / (1 - exp(fF * F1));  
	//      TOPMODEL, R2=-XKF*SZF*(CD+F1)/(1-EXP(SZF*F1))

	//Ponding begins due to (Flux_to_Infiltration_m rate, bPP) > (rainfall rate when ponding begins, RR), and then infiltration excess runoff generated
	if(RR < bPP) {
	//Note: 	TOPMODEL, IF(R2.LT.RINT)THEN
		//pondTime (time units) is time ponding occurs = timestep number * timestep duration
		beC->by_key(DataDrawer_ID, DataFolder_ID, "pondTime") = bIT*bDT;             
		//Note: 	TOPMODEL, TP=(IT-1.)*DT
		if (beC->by_key(DataDrawer_ID, DataFolder_ID, "pondTime") < 0.0) {
			beC->by_key(DataDrawer_ID, DataFolder_ID, "pondTime") = 0.0;
		}

		//Flag ponding = true (1)
		//beC->by_key(DataDrawer_ID, DataFolder_ID, "ponding") = 1;
		//Note: 	TOPMODEL, IROF=1
		//F (m) is cumulative infiltration
		beC->by_key(DataDrawer_ID, DataFolder_ID, "F") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m");
		//Note: 	TOPMODEL, F=CUMF

		//prepPonding prepare parameters for infiltration during ponding
		prepPonding(bIT,bPP,bDT, input, beC, DataDrawer_ID, DataFolder_ID);
		//      TOPMODEL, GO TO 8 ... SET UP DEFINITE INTEGRAL CONSTANT USING FP
	}
	//Ponding does not begin, (Flux_to_Infiltration_m rate, bPP) < (rainfall rate when ponding begins, RR)
	else {
		//estimateTP to get estimate of when ponding will occur
		estimateTP(bIT, bR, bDT, input, beC, DataDrawer_ID, DataFolder_ID, timeStep);
		//Note: 	TOPMODEL, INITIAL ESTIMATE OF TIME TO PONDING
	}
}

//Routine from USGS Topmodel Fortran Code from Wolock (USGS)/Bevin, subroutine titled "INITIAL ESTIMATE OF TIME TO PONDING &     7 F2=CUMF+DT*RINT"
//Option: infiltrate all water if ponding has not occured, estimate time when ponding will occur based on current Flux_to_Infiltration_m rate; tpIT is timestep number, tpR is Flux_to_Infiltration_m for no ponding situation, tpDT is timestep duration
void InfiltrationExponentialDecay::estimateTP(int tpIT,  double tpR, double tpDT, Inputs *input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//K0 hydraulic conductivity (m / timestep) was adjusted to per time step with input->SimulationTimeStep_Duration_sec[timeStep] 
	double K0 = beC->by_key(DataDrawer_ID, DataFolder_ID, "HydraulicConductivity_Infiltration_mpdt");
	//      TOPMODEL, XKF = XK0; define hydraulic conductivity

	//mM is scaling parameter used to adjust infiltration rates from Config inputs; TOPMODEL theory; Beven et al. (1995a, 1995b)
	//m (m) is scaling parameter describing exponential decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
	//Beven, Lamb, et al. (1995a) 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 * 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. (1995b). 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"];
	//fF = 1/mM, scaling parameter for infiltration rates, taken directly from TOPMODEL code SZF = 1./SZM
	//fF could be thetaSoilMoisture/mMin, given in text below Eq 10 of Wang et al. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
	double fF = 1 / mM;
	//CCSF (m) coefficient is product of WFS_m (m), wetting front suction, and Green Ampt term DWSMI (delta soil moisture), Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3
	double CSSF = beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_WettingFront_Suction_m")*beC->by_key(DataDrawer_ID, DataFolder_ID, "DWSMI");

	//tpPP is Flux_to_Infiltration_m Rate (m / timestep) 
	double tpPP = tpR / tpDT;
	//FI is cumulative infiltration at start
	double F1 = 0.0;

	//RR (m / timestep) value of rainfall rate when it equals potential infiltration rate at start of ponding, based on Green Ampt theory 
	//Theory from Eq. 2 to 8 in Beven (1984). Infiltration into a class of vertically non-uniform soils. Hydrological Sciences Journal 29: 425–434.
	double RR = 0.0;
	//E (m) is epsilon error threhold for closure of iteration to find cumulative infiltration
	double E = 0.00001;
	//FF (m) is estimate of cumulative infiltration
	double FF = 0.0;
	//F2 (m) is cumulative infiltration + Flux_to_Infiltration_m
	double F2 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") + tpR;
	//Note: 	TOPMODEL, F2=CUMF+DT*RINT
	
	//Ponding will not occur during timestep, cumulative infiltration + Flux_to_Infiltration_m = 0, from TOPMODEL code as IF(F2.EQ.0.)GO TO 20
	if (F2 == 0.0) {
		//noPonding routine called, all Flux_to_Infiltration_m (tpR) is infiltrated as infilQ
		noPonding(tpR, input, beC, DataDrawer_ID, DataFolder_ID, timeStep);
	}

	//Ponding may occur during timestep, cumulative infiltration + runoffToIfil > 0
	else 
	{
		//RR (m / timestep) value of rainfall rate when it equals potential infiltration rate at start of ponding, based on Green Ampt theory 
		//Theory from Eq. 2 to 8 in Beven (1984). Infiltration into a class of vertically non-uniform soils. Hydrological Sciences Journal 29: 425–434.
		//Note: 	Before ponding, potential infiltration rate > rainfall rate
		//Note: 	Potential infiltration rate declines from maximum when cumulative infiltration = 0 toward minimum capacity, defined as K0
		//Note: 	Rainfall rate has to exceed potential infiltration rate for ponding to occur
		//RR is modified version of Eq 6 and 7 from Beven (1984) Infiltration paper; lifted directly from Beven TOPMODEL code. 
		//RR is represented in Eq 10 of Wang et al. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
		RR = -K0 * fF * (CSSF + F2) / (1 - exp(fF * F2));
		//      TOPMODEL, R2=-XKF*SZF*(CD+F2)/(1-EXP(SZF*F2))

		//Ponding does not begin due to (Flux_to_Infiltration_m rate, tpPP) < (rainfall rate when ponding begins, RR), from TOPMODEL code IF(R2.GT.RINT)GO TO 20
		if (RR > tpPP) {
			//noPonding routine called, all Flux_to_Infiltration_m (tpR) is infiltrated as infilQ
			noPonding(tpR, input, beC, DataDrawer_ID, DataFolder_ID, timeStep);
		}

		//Ponding does begin due to (Flux_to_Infiltration_m rate, tpPP) > (rainfall rate when ponding begins, RR), from TOPMODEL code IF(R2.GT.RINT)GO TO 20
		else
		{
			//initial estimate of F (m) cumulative infiltration when ponding occurs = cumulative infiltration + (rainfall rate when ponding begins) * timestep duration, from TOPMODEL code F=CUMF+R2*DT
			beC->by_key(DataDrawer_ID, DataFolder_ID, "F") =  beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") + tpDT*RR;
			//      TOMODEL, F=CUMF+R2*DT
			
			//for loop to find F (m) cumulative infiltration when ponding occurs, from TOPMODEL code = DO 9 I=1,20
			for(int i=0; i<20; i++)               
			{
				//RR (m / timestep) value of rainfall rate when it equals potential infiltration rate at start of ponding, based on Green Ampt theory 
				//RR is modified version of Eq 6 and 7 from Beven (1984) Infiltration paper; lifted directly from Beven TOPMODEL code. 
				RR =-K0 * fF * (CSSF+ beC->by_key(DataDrawer_ID, DataFolder_ID, "F")) / (1-exp(fF* beC->by_key(DataDrawer_ID, DataFolder_ID, "F")));
				//      TOPMODEL, R2=-XKF*SZF*(CD+F)/(1-EXP(SZF*F))
		
				//Ponding does not begin due to (Flux_to_Infiltration_m rate, tpPP) < (rainfall rate when ponding begins, RR)
				if(RR > tpPP) {
				//      TOPMODEL, IF(R2.GT.RINT)THEN
					//FI is cumulative infiltration = beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") + tpDT*RR, from TOPMODEL F1=F
					F1 = beC->by_key(DataDrawer_ID, DataFolder_ID, "F");
					//Note:   TOPMODEL, F1=F
					//F is average of two estimates of cumulative infiltration, F2 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") + tpR and beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") + tpDT*RR;
					beC->by_key(DataDrawer_ID, DataFolder_ID, "F") = 0.5 * (F2 + beC->by_key(DataDrawer_ID, DataFolder_ID, "F"));
					//Note:   TOPMODEL, F=(F2+F)*0.5
					//FF (m) is estimate of cumulative infiltration at ponding from prior timestep
					FF = F1;
				}

				//Ponding does begin due to (Flux_to_Infiltration_m rate, tpPP) < (rainfall rate when ponding begins, RR)
				else
				{
					//FI is cumulative infiltration = beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") + tpDT*RR, from TOPMODEL F1=F, from TOPMODEL F2=F
					F2 = beC->by_key(DataDrawer_ID, DataFolder_ID, "F");
					//Note:   TOPMODEL, F2=F
					//F is average of two estimates of cumulative infiltration, F1 = 0 and beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") + tpDT*RR;
					beC->by_key(DataDrawer_ID, DataFolder_ID, "F") =  0.5 * (F1 + beC->by_key(DataDrawer_ID, DataFolder_ID, "F"));
					//Note:   TOPMODEL, F=(F1+F)*0.5
					//FF (m) is estimate of cumulative infiltration at ponding from prior timestep
					FF = F2;
				}

				//close for loop to estimate ponding if change in estimated cumulative infiltration between timesteps < error threshold
				if(fabs(beC->by_key(DataDrawer_ID, DataFolder_ID, "F") - FF) < E) {
				//      TOPMODEL, IF(ABS(F-F1).LT.E)GO TO 11

					//pondTime (time units) is time ponding occurs = (timestep number * timestep duration) - interval of time for additional infiltration at given runoffToInfilRate
					beC->by_key(DataDrawer_ID, DataFolder_ID, "pondTime") = tpIT * tpDT + (beC->by_key(DataDrawer_ID, DataFolder_ID, "F") - (beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m"))) / tpPP;
					//      TOPMODEL, TP=(IT-1)*DT+(F-CUMF)/RINT

					//no ponding occurs in this timestep if pondTime occurs after the start of the next time step
					if(beC->by_key(DataDrawer_ID, DataFolder_ID, "pondTime") > ((tpIT + 1) * tpDT) ) {
					//      TOPMODEL, IF(TP.GT.IT*DT)GO TO 20
						//noPonding routine called, all Flux_to_Infiltration_m (tpR) is infiltrated as infilQ
						noPonding(tpR, input, beC, DataDrawer_ID, DataFolder_ID, timeStep);
						//End iteration
						break;
					}
			
					//ponding occurs in this timestep 
					//prepPonding prepare parameters for infiltration during ponding
					prepPonding(tpIT,tpPP,tpDT, input, beC, DataDrawer_ID, DataFolder_ID);
					//End iteration
					break;
				}
			
			}
		
		}
	
	}
}

//Routine from USGS Topmodel Fortran Code from Wolock (USGS)/Bevin, subroutine titled "NEWTON-RAPHSON SOLUTION FOR F(T)"
//Option: find infiltration rate during ponding; ipIT is timestep number, ipPP is Flux_to_Infiltration_m for no ponding situation, ipDT is timestep duration
void InfiltrationExponentialDecay::inPonding(int ipIT, double ipPP, double ipDT, Inputs *input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID)
{

	//K0 hydraulic conductivity (m / timestep) was adjusted to per time step with input->SimulationTimeStep_Duration_sec[timeStep] 
	double K0 = beC->by_key(DataDrawer_ID, DataFolder_ID, "HydraulicConductivity_Infiltration_mpdt");
	//      TOPMODEL, XKF = XK0; define hydraulic conductivity

	//mM is scaling parameter used to adjust infiltration rates from Config inputs; TOPMODEL theory; Beven et al. (1995a, 1995b)
	//m (m) is scaling parameter describing exponential decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
	//Beven, Lamb, et al. (1995a) 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 * 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. (1995b). 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"];
	//fF = 1/mM, scaling parameter for infiltration rates, taken directly from TOPMODEL code SZF = 1./SZM
	//fF could be thetaSoilMoisture/mMin, given in text below Eq 10 of Wang et al. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
	double fF = 1 / mM;
	//CCSF (m) coefficient is product of WFS_m (m), wetting front suction, and Green Ampt term DWSMI (delta soil moisture), Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3
	double CSSF = beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_WettingFront_Suction_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "DWSMI");

	//variables for Newton-Raphson algorithm
	int FAC ;
	double ADD = 0.0;
	double SUM ;
	//root function
	double FUNC = 0.0;
	//derivative of root function
	double DFUNC = 0.0;
	double FC = 0.0;
	//E (m) is epsilon error threhold for closure of iteration to find cumulative infiltration
	double E = 0.00001;


	//infilQ (m) initialized to 0
	beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") = 0.0;

		//Newton-Raphson algorithm to infiltration rate during ponding
		//for loop to find infiltration, from TOPMODEL code = DO 14 I = 1, 20
		for(int i=0; i<20; i++)
		{
			//FC (m) is used in Lambda function, cumulative infiltration + Green Ampt coefficient (wetting front suction * delta theta)
			FC = beC->by_key(DataDrawer_ID, DataFolder_ID, "F") + CSSF;
			//      FC=(F+CD)
			//SUM initialized to 0 on each loop 
			SUM = 0.0;
			//      SUM=0.
			//FAC initialized to 1 on each loop 
			FAC = 1;
			//      FAC=1.

			//for loop, from TOPMODEL code = DO 13 J=1,10
			for(int j = 1; j<10; j++)
			{
				//FAC is incremented by j factor for each loop
				FAC = FAC * j ;
				//      TOPMODEL, FAC=FAC*J
				//ADD is updated with j and FAC for each loop
				ADD = pow(FC * fF , j) / (j * FAC);
				//      TOPMODEL, ADD=(FC*SZF)**J/(J*FAC)
				//SUM (or constant) is updated each loop
				SUM = SUM + ADD ;
				//      TOPMODEL, SUM=SUM+ADD
			}
		
			//FUNC is Lambda defined below Eq 8 from Beven (1984) Infiltration paper; lifted directly from Beven TOPMODEL code
			//FUNC is Lambda in Eq 13 of Wang et al. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
			FUNC = -(log(FC) - (log(FC) + SUM) / exp(fF * CSSF) - beC->by_key(DataDrawer_ID, DataFolder_ID, "constant") ) / (K0 * fF)
				   - ((ipIT + 1) * ipDT - beC->by_key(DataDrawer_ID, DataFolder_ID, "pondTime")) ;
			//      TOPMODEL, FUNC=-(DLOG(FC)-(DLOG(FC)+SUM)/DEXP(SZF*CD)-CONST)/(XKF*SZF) - (IT*DT - TP)

			//Derivative of Lambda function for computing infiltration rate, from TOPMODEL, DFUNC=(EXP(SZF*F)-1)/(XKF*SZF*FC)
			DFUNC = ( exp(fF * beC->by_key(DataDrawer_ID, DataFolder_ID, "F")) - 1) / (K0 * fF * FC) ;
			//      TOPMODEL, DFUNC=(EXP(SZF*F)-1)/(XKF*SZF*FC)

			//Newton-Raphson (N-R) correction factor = (-FUNC/DFUNC) used to get InfilQ (m), incrementally updates cumulative infiltration to satisfy root
			beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") = (-FUNC / DFUNC);
			//      TOPMODEL, DF=-FUNC/DFUNC

			//N-R estimated F (m), cumulative infiltration, based on finding root that creates no additional infiltration to soil evapotranspiration zone
			beC->by_key(DataDrawer_ID, DataFolder_ID, "F") = beC->by_key(DataDrawer_ID, DataFolder_ID, "F") + beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ");
			//Note:      TOPMODEL, F=F+DF

			//close for loop to estimate cumulative infiltration if incremental improvement of cumulative infiltration estimate is less than threshold
			if (fabs(beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ")) < E) {
				//Note:      TOPMODEL, IF(ABS(DF).LE.E)GO TO 15
					//End iteration
				break;
			}
		}

	//Infiltration excess constraint when estimated N-R estimate of cumulative infiltration (m) < (prior cumulative infiltration + Flux_to_Infiltration_m)
	if(beC->by_key(DataDrawer_ID, DataFolder_ID, "F") < (beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") + ipPP) )
	//Note:      TOPMODEL, IF(F.LT.CUMF+RINT)THEN
	{
		//infilQ (m) is infiltration for timestep = maximum possible cumulative infiltration during timestep - cumulative infiltration from prior timestep
		beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") = (beC->by_key(DataDrawer_ID, DataFolder_ID, "F") - beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m"));
		//Note: 	TOPMODEL, DF=F-CUMF

		//Prohibit negative values for infilQ (m)
		if (beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") < 0.0) {
			//infilQ (m) is reset to 0 if negative; negative value may emerge during Newton-Raphson root finding, but should not emerge or cary forward from here; had been set to 0.00001
			beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") = 0.0;
		}

		//cumInfil (m) updated to new value at end of timestep
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_viaInfilEx_cumulative_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "F");
		//Note: 	TOPMODEL, CUMF=F
		//F (m) is updated to include the water infiltrated this timestep
		beC->by_key(DataDrawer_ID, DataFolder_ID, "F") = beC->by_key(DataDrawer_ID, DataFolder_ID, "F") + beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ");
		//Note: 	TOPMODEL, F=F+DF

	}
	//Prohibit negative values for infilQ (m)
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") < 0.0) {
		//infilQ (m) is reset to 0 if negative; negative value may emerge during Newton-Raphson root finding, but should not emerge or cary forward from here; had been set to 0.00001
		beC->by_key(DataDrawer_ID, DataFolder_ID, "infilQ") = 0.0;
	}
}

//Routine from USGS Topmodel Fortran Code from Wolock (USGS)/Bevin, subroutine titled "SET UP DEFINITE INTEGRAL CONSTANT USING FP"
//Option: prepare parameters for infiltration during ponding; prIT is timestep number, prPP is Flux_to_Infiltration_m for no ponding situation, prDT is timestep duration
void InfiltrationExponentialDecay::prepPonding(int prIT, double prPP, double prDT, Inputs *input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID)
{
	
	//mM is scaling parameter used to adjust infiltration rates from Config inputs; TOPMODEL theory; Beven et al. (1995a, 1995b)
	//m (m) is scaling parameter describing exponential decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
	//Beven, Lamb, et al. (1995a) 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 * 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. (1995b). 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"];

	//fF = 1/mM, scaling parameter for infiltration rates, taken directly from TOPMODEL code SZF = 1./SZM
	//fF could be thetaSoilMoisture/mMin, given in text below Eq 10 of Wang et al. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
	double fF = 1 / mM;
	//CCSF (m) coefficient is product of WFS_m (m), wetting front suction, and Green Ampt term DWSMI (delta soil moisture), Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3
	double CSSF = beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_WettingFront_Suction_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "DWSMI");
	//F (m) is cumulative infiltration 
	double F = beC->by_key(DataDrawer_ID, DataFolder_ID, "F");

	//variables for Newton-Raphson algorithm
	int FAC = 1;
	//Note: TOPMODEL, FAC=1
	double ADD = 0.0;
	beC->by_key(DataDrawer_ID, DataFolder_ID, "constant") = 0.0;
	//Note: TOPMODEL, CONST =0
	double FC = beC->by_key(DataDrawer_ID, DataFolder_ID, "F") + CSSF;
	//Note: TOPMODEL, FC=(F+CD)

	//for loop, from TOPMODEL code = DO 12 J=1,10
	for(int j = 1; j<=10; j++)
	{
		//FAC is incremented by j factor for each loop
		FAC = FAC * j;
		//Note: TOPMODEL, FAC=FAC*J

		//ADD is updated with j and FAC for each loop
		ADD = pow(FC * fF,j) / (j * FAC);
		//Note: TOPMODEL, ADD = (FC * SZF) * *J / (J * FAC)

		//constant (or SUM) is updated each loop
		beC->by_key(DataDrawer_ID, DataFolder_ID, "constant") = beC->by_key(DataDrawer_ID, DataFolder_ID, "constant") + ADD;
		//Note: TOPMODEL, CONST=CONST+ADD
	}
	
	//constant is Lambda defined below Eq 8 from Beven (1984) Infiltration paper; lifted directly from Beven TOPMODEL code
	//constant is Lambda in Eq 13 of Wang et al. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
	beC->by_key(DataDrawer_ID, DataFolder_ID, "constant") = log(FC) - (log(FC) + beC->by_key(DataDrawer_ID, DataFolder_ID, "constant")) / exp(fF * CSSF);
	//Note: TOPMODEL, CONST=DLOG(FC)-(DLOG(FC)+CONST)/DEXP(SZF*CD)

	//Flag ponding = true (1)
	//beC->by_key(DataDrawer_ID, DataFolder_ID, "ponding") = 1;
	//Note: TOPMODEL, IROF=1

	//F (m) cumulative infiltration is average of Flux_to_Infiltration_m at end of time step and at time of ponding
	beC->by_key(DataDrawer_ID, DataFolder_ID, "F") = beC->by_key(DataDrawer_ID, DataFolder_ID, "F") + 0.5 * prPP * ( (prIT + 1) * prDT -  beC->by_key(DataDrawer_ID, DataFolder_ID, "pondTime") );
	//TOPMODEL, F=F+0.5*RINT*(IT*DT-TP)

	//inPonding routine called
	inPonding(prIT, prPP, prDT, input, beC, DataDrawer_ID, DataFolder_ID);
	

	
}