#include "RunoffSubsurface.h"
#include <math.h>
#include <fstream>
#include <vector>
#include <iostream>
#include <numeric>
#include <string>

//Theory of exponential decay for computing Runof_Subsurface_m (m) emerged with Beven and Kirby (1979) and is documented by Wolock (1993)
//Theory of power function decay for computing Runof_Subsurface_m (m) emerged from exponential decay by Wolock (1993) and Beven and Kirby (1979)
//Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124). 
//Beven, K. J., & Kirkby, M. J. (1979). A physically based, variable contributing area model of basin hydrology. Hydrological Sciences Bulletin, 24(1), 43-69. doi:10.1080/02626667909491834

//Note: Theory of power function decay requires m <= S (or AveSMD) based on constraint governing decay with depth of transmissivity, Tz = T0*(1-S/m)^n  
//S (m) is soil moisture deficit between top of soil and groundwater table, defined as 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, Editor
//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(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

//RunoffSubsurface::calculate function computes Runof_Subsurface_m (m), depth of saturated zone runoff from pervious area, using exponential or power decay approach
//Note: StorageDeficit_VadoseZone_TS_Prior_m (m) used to avoid use of input->RepoDict["StorageDeficit_VadoseZone_m"], which changes throughout the time step
//Note: input->InputXml["StorageDeficit_VadoseZone_TS_Prior_m"] (m) is defined in SimulationCoordinator at start of each For Loop through DataDrawers, 
void RunoffSubsurface::calculate(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//T0_m2pdt transmissivity (m^2 / timestep) when AveSMD = 0 (watershed fully saturated) is adjusted to per time step 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];
	//m (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 * 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.
	//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(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.
	//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. 
	double mM = input->InputXml["Parameter_m_KsatExpDecay"];
	//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"];
	//Runoff_Subsurface_Max_m (m) initialized to 0.0
	double Runoff_Subsurface_Max_m = 0.0;

	//If SimulationStringParams["Algorithm_SoilKsatDecay" is ExponentialDecay then
	if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "ExponentialDecay") {
		//Runoff_Subsurface_Max_m (m) is maximum subsurface runoff from pervious watershed when soil fully saturated; multiply by watershed area to get volume (m3); runoff is per time step
		//Eq 32 from Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124). 
		Runoff_Subsurface_Max_m = T0_m2pdt * exp(-input->TI_Value_Average);
		//TOPMODEL Code: SZQ = EXP(T0DT - TL), where TL is average topographic index, and T0DT is natural log of transmissivity, ln(T0)

		//Runoff_Subsurface_m (m), pervious area depth, is saturation zone runoff from pervious area for actual level of saturation, using AveSMD and TOPMODEL relationship
		//Eq 28 from Wolock (1993), with 1st RHS term swapped out by Eq 31 (which replaces Eq 29 for spatial aggregation). Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124). 
		//Eq 3 from Beven and Kirby (1979) gives qb = q0 exp (S/m), where q0 is flow when S = 0 or storage deficit is 0
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_Subsurface_m") = Runoff_Subsurface_Max_m * exp(-input->InputXml["StorageDeficit_VadoseZone_TS_Prior_m"] / mM);
		//TOPMODEL Code: QB = SZQ * EXP(-SBAR / SZM), where SZQ is SatZoneQ (m), SBAR is AveSMD, SZM is m
	}
	//If SimulationStringParams["Algorithm_SoilKsatDecay" is PowerDecay then
	else if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") {
		//Runoff_Subsurface_Max_m (m) is maximum subsurface runoff from pervious watershed when soil fully saturated; multiply by watershed area to get volume (m3); runoff is per time step
		//Eq 5, 1st RHS term from 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. 
		Runoff_Subsurface_Max_m = T0_m2pdt * pow(input->TI_Value_Average, -nN);

		//Theory of power function decay requires m <= AveSMD or S based on constraint governing decay with depth of transmissivity, Tz = T0*(1-S/m)^n  
		//S (m) is soil moisture deficit, defined as Theta_sat-Theta_fc (based on Beven and Wood, 1983).
		//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(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
		//Note: StorageDeficit_VadoseZone_TS_Prior_m (m) is defined in SimulationCoordinator at start of each For Loop through DataDrawers, 
		//Note:	StorageDeficit_VadoseZone_TS_Prior_m (m) avoids use of iteratively changing input->RepoDict["StorageDeficit_VadoseZone_m"], updated in StorageDeficit_VadoseZone for each DataDrawer
		if (input->InputXml["StorageDeficit_VadoseZone_TS_Prior_m"] > mM) {input->InputXml["StorageDeficit_VadoseZone_TS_Prior_m"] = (mM);}

		//Runoff_Subsurface_m (m), pervious area depth, is saturation zone runoff from pervious area for actual level of saturation, using AveSMD and TOPMODEL relationship
		//Eq 5, 2nd RHS term from 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. 
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_Subsurface_m") = Runoff_Subsurface_Max_m * pow((1 - input->InputXml["StorageDeficit_VadoseZone_TS_Prior_m"] / mM), nN);
	}
}