#include "Runoff_StormwaterDevice.h"

//Runoff_Surface_Outlet is providing controlled outflow from GI
//Note: Refactor to allow option of orifice controlled discharge vs Manning n controlled discharge
double Runoff_StormwaterDevice::Runoff_Surface_Outlet(Inputs* input, CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep) {

	//Note: EPA SWMM LID Devices using getSurfaceOutflowRate are: biocellFluxRates, greenRoofFluxRates, trenchFluxRates, pavementFluxRates, roofFluxRates
	//Note: Other SWMM LID functions are: barrelFluxRates, swaleFluxRates

	//Head_GI_Outlet_m (m) is water for controlled release through GI outlet, which may be set to height of 0 m
	double Head_GI_Outlet_m = 0.0;
	double Storage_SurfaceLayer_Available_m = 0.0;
	double Storage_GI_SurfaceDepression_Max_m3 = 0.0;
	double Runoff_GI_Outlet_Available_m3 = 0.0;
	double Runoff_GI_Outlet_m3ps = 0.0;
	double Runoff_GI_Outlet_mps = 0.0;
	double Runoff_GI_Outlet_m3= 0.0;
	double Runoff_GI_Outlet_m = 0.0;
	double Area_GI_CrossSection_m2 = 0.0;
	double Perimeter_GI_CrossSection_m = 0.0;
	double TopWidth_GI_CrossSection_m = 0.0;
	double Area_GI_Unit_Wetted_m2 = 0.0;
	double HydraulicRadius_GI_m = 0.0;

	//If Storage_GI_Surface_Max_m3 > 0 then enter for division 
	//Note: Depth may not be stored in a rectangular GI unit due to Surface_SideSlope_HtoV_mpm, and storage volume ratios handle this ambiguity
	//Note: If Storage_GI_Surface_Max_m3 = 0, then Surface_SideSlope_HtoV_mpm, Surface_Berm_Height_m, Surface_Berm_Height_m, or Surface_Porosity_m3pm3 = 0
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_Max_m3") > 0 && beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm") > 0 && beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Length_GI_Unit_m") > 0) {

		//If (Storage_GI_Surface_m3 - Storage_GI_SurfaceDepression_m3) >= 0.75 * Storage_GI_Surface_Max_m3 then nearing top of channel 
		if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3") - beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_SurfaceDepression_m3") >= 0.75 * beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_Max_m3")) {
			//Storage_SurfaceLayer_Available_m (m) equals Surface_Berm_Height_m (m) * (Storage_GI_Surface_m3 - Storage_GI_SurfaceDepression_m3) / Storage_GI_Surface_Max_m3 (m3)
			//Note: Right hand side term is a ratio of volumes, and provides fractional scaling of Surface_Berm_Height_m (m)
			//Note: Equation works well at 75% volume or greater, but tends to underestimate at lower volumes
			//Note: Surface_Porosity_m3pm3 (m3/m3) incorporated into Storage_GI_Surface_Max_m3 (m3), and hence into Storage_SurfaceLayer_Available_m (m)
			Storage_SurfaceLayer_Available_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Berm_Height_m") * (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3") - beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_SurfaceDepression_m3")) / beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_Max_m3");
		}
		//Else (Storage_GI_Surface_m3 - Storage_GI_SurfaceDepression_m3) < 0.75 * Storage_GI_Surface_Max_m3 then in lower part of channel
		else {
			//Storage_SurfaceLayer_Available_m (m) derived from relation between depth and volume for trapezoidal channel
			//Note: V = (B*D + D^2*m) * L, then positive quadratic equation gives D = (-B+(B^2+4*m*V/L)^0.5)/(2*m) where B is bottom width, ...
			//Note: ... D is depth, V is volume, m is side slope, and L is length
			//Note: Equation works well at less than 75% volume, but tends to underestimate at greater volumes
			//Note: Equation does not adjut for Surface_Porosity_m3pm3, to simplify assumptions in estimates of multi-threaded channels through vegetation
			Storage_SurfaceLayer_Available_m = (-beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_BottomWidth_m")+ pow((pow(beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_BottomWidth_m"), 2) + 4 * beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm") * (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3") - beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_SurfaceDepression_m3")) / beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Length_GI_Unit_m")), 0.5)) / (2 * beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm"));
			
			//Storage_SurfaceLayer_Available_m (m) equals Storage_SurfaceLayer_Available_m (m) / Surface_Porosity_m3pm3 (m3/m3)
			//Note: Surface_Porosity_m3pm3 (m3/m3) incorporated into Storage_GI_Surface_Max_m3 (m3), and hence into Storage_SurfaceLayer_Available_m (m)
			//Note: Surface_Porosity_m3pm3 (m3/m3) has minimum value of Constant_1E_negative10, so no chance for division by zero
			Storage_SurfaceLayer_Available_m = Storage_SurfaceLayer_Available_m / beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Porosity_m3pm3");
		}
	}
	//Else Storage_GI_Surface_Max_m3 <= 0 and determine Flux_GI_SurfaceDepression_m3 (m3) using Area_GI_Unit_m2 (m2)
	else {
		//Storage_SurfaceLayer_Available_m (m) equals (Storage_GI_Surface_m3 - Storage_GI_SurfaceDepression_m3) / Area_GI_Unit_m2 (m2)
		Storage_SurfaceLayer_Available_m = (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3") - beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_SurfaceDepression_m3")) / beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2");

		//Storage_SurfaceLayer_Available_m (m) equals Storage_SurfaceLayer_Available_m (m) / Surface_Porosity_m3pm3 (m3/m3)
		//Note: Surface_Porosity_m3pm3 (m3/m3) incorporated into Storage_GI_Surface_Max_m3 (m3), and hence into Storage_SurfaceLayer_Available_m (m)
		//Note: Surface_Porosity_m3pm3 (m3/m3) has minimum value of Constant_1E_negative10, so no chance for division by zero
		Storage_SurfaceLayer_Available_m = Storage_SurfaceLayer_Available_m / beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Porosity_m3pm3");
	}

	//Head_GI_Outlet_m (m) is head above Storage_GI_SurfaceDepression_Max_m3 and above Surface_Outlet_Offset_m (m)
	//Note: Head_GI_Outlet_m is not adjusted upwards for Surface_Porosity_m3pm3; doing so without elaborate algorithm generates instability
	Head_GI_Outlet_m = Storage_SurfaceLayer_Available_m - beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Outlet_Offset_m");
	
	//If Head_GI_Outlet_m (m) > Surface_Berm_Height_m (m) then set Head_GI_Outlet_m (m) = Surface_Berm_Height_m (m), which approximates physical limit
	//Note: Head_GI_Outlet_m (m) is adjusted for Surface_Porosity_m3pm3 (m3/m3)
	if (Head_GI_Outlet_m > beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Berm_Height_m")) { Head_GI_Outlet_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Berm_Height_m"); }

	//If Head_GI_Outlet_m (m) is less than or equal to 0, then Runoff_GI_Outlet_m3 equals zero
	if (Head_GI_Outlet_m <= 0.0) {
		Runoff_GI_Outlet_m3 = 0.0;
		return Runoff_GI_Outlet_m3;
	}
	
	//If Surface_BottomWidth_m (m) is > 0 then allow for controlled outflow using Manning function
	//Note: For GI type swale, Surface_BottomWidth_m (m) set to max of HydroPlusConfig.xml input or Surface_Outlet_Width_Min_m in BuildDataOrganizer
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_BottomWidth_m")> 0) {

		//Area_GI_CrossSection_m2 (m2) is cross-sectional flow area for trapezoidal channel, for depth above Surface_Outlet_Offset_m; 
		//Note: Eq 5.3 from Chin (2013) A = b*y + m*y^2, where b = bottom width, y = depth, m = side slope
		//Note: EPA SWMM for Swale computes these terms using the top width of depression storage, but HydroPlus takes a different approach
		Area_GI_CrossSection_m2 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_BottomWidth_m")* Head_GI_Outlet_m +
			beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm") * pow(Head_GI_Outlet_m, 2.0);
		//Perimeter_GI_CrossSection_m (m) is cross-sectional flow wetted perimeter for trapezoidal channel
		//Note: Eq 5.4 from Chin (2013) P = b + 2*y(1+m^2)^0.5, where b = bottom width, y = depth, m = side slope
		Perimeter_GI_CrossSection_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_BottomWidth_m")+ 2.0 * Head_GI_Outlet_m *
			sqrt(1.0 + beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm") * beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm"));
		//HydraulicRadius_GI_m (m) is ratio of Area_GI_CrossSection_m2 (m2) and Perimeter_GI_CrossSection_m (m)
		HydraulicRadius_GI_m = Area_GI_CrossSection_m2 / Perimeter_GI_CrossSection_m;

		//TopWidth_GI_CrossSection_m (m) is cross-sectional flow wetted top width for trapezoidal channel
		//Note: Top width equation for trapezoidal channel from Fig 5.2 in Chin (2021); Calculates wetted top width
		//Note: Equation will work for rectangular (Surface_SideSlope_HtoV_mpm=0) or triangular (Surface_BottomWidth_m=0) if variable = zero
		TopWidth_GI_CrossSection_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_BottomWidth_m")+ 2 * beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm") * Head_GI_Outlet_m;

		//Area_GI_Unit_Wetted_m2 (m2) equals Surface_Length_GI_Unit_m (m) * TopWidth_GI_CrossSection_m (m)
		//Note: Surface_Porosity_m3pm3 (m3/m3) incorporated into Storage_SurfaceLayer_Available_m (m), not Area_GI_Unit_Wetted_m2 (m2)
		Area_GI_Unit_Wetted_m2 = (beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Length_GI_Unit_m") * TopWidth_GI_CrossSection_m);

		//If Surface_ManningRoughness greater than 0 then enter for division, GI device can release water with Manning Eq
		if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_ManningRoughness") > 0) {

			//If GI type Swale then use one form of Manning Eq
			if (beC->by_key_str(DataDrawer_ID, DataFolder_ID, "Type") == "Swale") {
				//Runoff_GI_Outlet_m3ps (m3/s) is computed with Manning Eq to obtain volumetric flow (m3/s)
				Runoff_GI_Outlet_m3ps = pow(HydraulicRadius_GI_m, (2.0 / 3.0)) * Area_GI_CrossSection_m2 * sqrt(beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_LongitudinalSlope_mpm")) / beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_ManningRoughness");
				//If Area_GI_Unit_Wetted_m2 (m2) < 0 then set to Area_GI_Unit_m2
				if (Area_GI_Unit_Wetted_m2 < 0) {
					//Area_GI_Unit_Wetted_m2 (m2) equals Area_GI_Unit_m2 (m2)
					//Note: Surface_Porosity_m3pm3 (m3/m3) incorporated into Storage_SurfaceLayer_Available_m, not Area_GI_Unit_Wetted_m2
					Area_GI_Unit_Wetted_m2 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2");
				}
				//Runoff_GI_Outlet_m (m) equals (Runoff_GI_Outlet_m3ps (m3/s) * SimulationTimeStep_Duration_sec[timeStep]) / Area_per_GI_Unit_m2 (m2) 
				//Note: Division by Area_GI_Unit_Wetted_m2 scales flow to its actual surface area, not Area_GI_Unit_m2
				//Note: Area_GI_Unit_Wetted_m2 (m2) can change across timesteps without violating the GI water balance
				//Note: Runoff_GI_Outlet_m (m) * Area_GI_Unit_m2 (m2) becoming Runoff_GI_Outlet_m3 (m3) maintains water balance
				Runoff_GI_Outlet_m = input->SafeDivide(Runoff_GI_Outlet_m3ps * input->SimulationTimeStep_Duration_sec[timeStep], Area_GI_Unit_Wetted_m2);
			}
			//Else if all other GI devices but swale, then use other form of Manning Eq
			else {

				//Runoff_GI_Outlet_m3ps (m3/s) based on Manning Eq from SWMM code lidproc.c and Eq 3.3 of SWMM Ref Manual Vol I - Hydrology (2016) 
				//Derivation:	Q (m/s) = [GI_Head^(5/3) * GI_Width * Slope^(1/2) / n] / GI_Area_planform; GI_Area_planform = GI_Width * GI_Length
				//Note: Q (m/s) is runoff depth across GI, obtained by dividing Q (m^3/s) by GI_Area_planform
				//Note: GI_Depth^(5/3) = R^(2/3) * GI_Area_xs ^ (3/3) if R = GI_Head && GI_Area_xs = GI_Head x GI_Width
				//Note: Q (m^3/s) = Head_GI_Outlet_m ^ (5/3) * Surface_BottomWidth_m * Surface_LongitudinalSlope_mpm ^ 0.5 / Surface_ManningRoughness
				//Note: Head_GI_Outlet_m = GI_Head
				Runoff_GI_Outlet_m3ps = pow(Head_GI_Outlet_m, 5.0 / 3.0) * beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_BottomWidth_m")* input->SafeDivide(sqrt(beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_LongitudinalSlope_mpm")), beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_ManningRoughness"));
				//If Area_GI_Unit_Wetted_m2 (m2) < 0 then set to Area_GI_Unit_m2
				if (Area_GI_Unit_Wetted_m2 < 0) {
					//Area_GI_Unit_Wetted_m2 (m2) equals Area_GI_Unit_m2 (m2)
					//Note: Surface_Porosity_m3pm3 (m3/m3) incorporated into Storage_SurfaceLayer_Available_m, not Area_GI_Unit_Wetted_m2
					Area_GI_Unit_Wetted_m2 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2");
				}
				//Runoff_GI_Outlet_m (m) equals (Runoff_GI_Outlet_m3ps (m3/s) * SimulationTimeStep_Duration_sec[timeStep]) / Area_per_GI_Unit_m2 (m2) 
				//Note: Division by Area_GI_Unit_Wetted_m2 scales flow to its actual surface area, not Area_GI_Unit_m2
				//Note: Area_GI_Unit_Wetted_m2 (m2) can change across timesteps without violating the GI water balance
				//Note: Runoff_GI_Outlet_m (m) * Area_GI_Unit_m2 (m2) becoming Runoff_GI_Outlet_m3 (m3) maintains water balance
				Runoff_GI_Outlet_m = input->SafeDivide(Runoff_GI_Outlet_m3ps * input->SimulationTimeStep_Duration_sec[timeStep], Area_GI_Unit_Wetted_m2);
			}

			//Runoff_GI_Outlet_m (m) is potentially controlled by rate of Manning Eq, minimum of Runoff_GI_Outlet_m and Head_GI_Outlet_m 
			Runoff_GI_Outlet_m = MIN(Runoff_GI_Outlet_m, Head_GI_Outlet_m);
		}
		//Else Surface_ManningRoughness <= 0 then do not use Manning Eq for calculation of Runoff_GI_Outlet_m
		else {
			//Runoff_GI_Outlet_m (m) equal to Head_GI_Outlet_m and all water immediately released
			Runoff_GI_Outlet_m = Head_GI_Outlet_m;
		}
	}
	//Else Surface_BottomWidth_m is <= 0 then Runoff_GI_Outlet_m equals zero; no controlled outflow
	//Note: Storage_GI_Surface_m3 (m3) may still exit with GI_RateOfSurfaceEmergencyOverFlow function
	else { Runoff_GI_Outlet_m = 0.0; }

	//Runoff_GI_Outlet_m3 (m3) is Runoff_GI_Outlet_m (m) over folder area multiplied by Area_GI_Unit_Wetted_m2 (m2)
	//Note: Runoff_GI_Outlet_m (m) * Area_GI_Unit_Wetted_m2 (m2) becoming Runoff_GI_Outlet_m3 (m3) maintains water balance
	Runoff_GI_Outlet_m3 = Runoff_GI_Outlet_m * Area_GI_Unit_Wetted_m2;

	//Velocity_GI_Outlet_Base_mps (m/s) set to 0.1, and then adjsuted below to allow variation across events
	double Velocity_GI_Outlet_Base_mps = 0.10;
	//If Runoff_GI_Outlet_m3 (m3) > 0 then and Type = Swale then enter to delay intial runoff, allowing arrival at outlet
	//Note: Delay of runoff only occurs at start of storm runoff, while Time_RunoffStarted_hr < time for initial runoff to arrive at outlet
	//Note: Delay is based on observed data for Platz et al. (2020) for UMD Bioswale shown in Figure 9 and in other shared data modeled by HydroPlus
	//Reference: Platz, M., M. Simon and M. Tryby (2020). "Testing of the Storm Water Management Model Low Impact Development Modules." JAWRA 56(2): 283-296, https://doi-org.esf.idm.oclc.org/10.1111/1752-1688.12832

	if (Runoff_GI_Outlet_m3 > 0 && beC->by_key_str(DataDrawer_ID, DataFolder_ID, "Type") == "Swale") {
		//Runoff_GI_Outlet_mps (m/s) equals Runoff_GI_Outlet_m3ps (m3/s) / Area_GI_CrossSection_m2 (m2)
		Runoff_GI_Outlet_mps = Runoff_GI_Outlet_m3ps / Area_GI_CrossSection_m2;
		//Time_RunoffStarted_hr (hr) is increased by Ratio_Hour_to_Second (hr/s) * SimulationTimeStep_Duration_sec[timeStep]
		//Note: Time_RunoffStarted_hr (hr) is not reset to 0 when runof ends, as it records first runoff of event
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Time_RunoffStarted_hr") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Time_RunoffStarted_hr") + Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep];
		//Velocity_GI_Outlet_Base_mps (m/s) adjusted by Runoff_GI_Outlet_m3ps / Area_GI_CrossSection_m2 (m/s) to allow variation across events
		//Note: Runoff_GI_Outlet_m3ps / Area_GI_CrossSection_m2 provides velocity estimate, but based on US EPA data this is too low
		Velocity_GI_Outlet_Base_mps = Velocity_GI_Outlet_Base_mps + Runoff_GI_Outlet_m3ps / Area_GI_CrossSection_m2;

		//If (Time_RunoffStarted_hr * Ratio_s_to_hr) < (Surface_Length_GI_Unit_m * Velocity_GI_Outlet_Base_mps) 
		if (((beC->by_key(DataDrawer_ID, DataFolder_ID, "Time_RunoffStarted_hr") * Ratio_s_to_hr) < beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Length_GI_Unit_m") / Velocity_GI_Outlet_Base_mps)) {
			//Runoff_GI_Outlet_m3 (m3) set to 0 for Type = Swale, retaining runoff for infiltration given the Swale is long and wave takes time to arrive
			Runoff_GI_Outlet_m3 = 0;
		}
	}

	//return Runoff_GI_Outlet_m3 (m3) to become Runoff_Surface_m3 (m3)
	return Runoff_GI_Outlet_m3;
}

//Runoff_Surface_Emergency is flow over top of berm
double Runoff_StormwaterDevice::Runoff_Surface_Emergency(CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID) {

	double Storage_SurfaceLayer_Available_m3 = 0.0;
	double Storage_GI_SurfaceDepression_Max_m3 = 0.0;
	double Runoff_GI_Berm_m3 = 0.0;

	//Storage_SurfaceLayer_Available_m3 (m3) equals Storage_GI_Surface_m3 (m3) - Storage_GI_SurfaceDepression_m3 (m3)
	Storage_SurfaceLayer_Available_m3 = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3") - beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_SurfaceDepression_m3");

	//If Storage_SurfaceLayer_Available_m3 (m3) > Storage_GI_Surface_Max_m3 (m3) then immediately remove any excess water via overtopping berm
	if (Storage_SurfaceLayer_Available_m3 > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_Max_m3")) {

		//Runoff_GI_Berm_m3 (m3) is volume released by Storage_SurfaceLayer_Available_m3 (m3) - Storage_GI_Surface_Max_m3 (m3)
		Runoff_GI_Berm_m3 = Storage_SurfaceLayer_Available_m3 - beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_Max_m3");
		//Storage_GI_Surface_m3 (m3) is reduced by Runoff_GI_Berm_m3 (m3)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3") - Runoff_GI_Berm_m3;
	}
	//Else Storage_SurfaceLayer_Available_m3 (m3) <= Storage_GI_Surface_Max_m3 (m3) then no runoff is generated, Runoff_GI_Berm_m3 (m3) = 0
	else {
		Runoff_GI_Berm_m3 = 0.0;
	}

	//return Runoff_GI_Berm_m3 (m3) to become Runoff_Surface_m3 (m3)
	return Runoff_GI_Berm_m3;
}

//Runoff_Vault_Drain function receives depths adjusted upwards by layer drainable porosity
//Note: Refactor to allow option of orifice controlled discharge vs Manning n controlled discharge
double Runoff_StormwaterDevice::Runoff_Vault_Drain(CompactRagged* beC, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	double Storage_SurfaceLayer_m = 0.0;
	double Head_SurfaceLayer_m = 0.0;
	double Storage_PavementLayer_m = 0.0;
	double Head_PavementLayer_m = 0.0;
	double Storage_SoilLayer_m = 0.0;
	double Head_SoilLayer_m = 0.0;
	double Storage_VaultLayer_m = 0.0;
	double Head_VaultLayer_m = 0.0;
	double Runoff_VaultDrain_m3 = 0.0;
	double Runoff_VaultDrain_m = 0.0;
	
	//If Vault_Outlet_Discharge_Coeff <= zero then return Runoff_VaultDrain_m3 as zero
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Outlet_Discharge_Coeff") <= 0) {
		Runoff_VaultDrain_m3 = 0.0;
		return Runoff_VaultDrain_m3;
	}

	//Storage volumes transformed to storage depths in pervious area
	//Note: Storage_SurfaceLayer_m (m) is quotient of (Storage_GI_Surface_m3 and Area_per_GIUnit_m2), divided by PerviousCover_frac
	Storage_SurfaceLayer_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3") / beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2");
	
	//Head_SurfaceLayer_m (m) initialized as Storage_SurfaceLayer_m (m) divided by Surface_Porosity_m3pm3 (m3/m3), increasing head due to solids
	//Note: Head_SurfaceLayer_m (m) includes water held in depression storage, as this head influences subsurface drainage
	//Note: Surface_Porosity_m3pm3 (m3/m3) has minimum value of Constant_1E_negative10, so no chance for division by zero
	Head_SurfaceLayer_m = Storage_SurfaceLayer_m / beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Porosity_m3pm3");

	//Note: Storage_PavementLayer_m (m) is quotient of (Storage_GI_Surface_m3 and Area_per_GIUnit_m2), divided by PerviousCover_frac
	Storage_PavementLayer_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Pavement_m3") / beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");
	//Head_PavementLayer_m (m) initialized as Storage_PavementLayer_m (m) divided by Pavement_Porosity_m3pm3 (m3/m3), increasing head due to solids
	//Note: Pavement_Porosity_m3pm3 (m3/m3) has minimum value of Constant_1E_negative10, so no chance for division by zero
	Head_PavementLayer_m = Storage_PavementLayer_m / beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Porosity_m3pm3");

	//Note: Storage_SoilLayer_m (m) is quotient of (Storage_GI_Surface_m3 and Area_per_GIUnit_m2), divided by PerviousCover_frac
	Storage_SoilLayer_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_m3") / beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");
	//Head_SoilLayer_m (m) initialized as Storage_SoilLayer_m (m) divided by Porosity_ThetaSat_ThetaFc_m3pm3 (m3/m3), increasing head due to solids
	Head_SoilLayer_m = Storage_SoilLayer_m / beC->by_key(DataDrawer_ID, DataFolder_ID, "Porosity_ThetaSat_ThetaFc_m3pm3");

	//Note: Storage_VaultLayer_m (m) is quotient of (Storage_GI_Surface_m3 and Area_per_GIUnit_m2), divided by PerviousCover_frac
	//Note: Vault is presumed to have area of Area_GI_Unit_m2 for exfiltration, with water otherwise infiltrating, percolating, and evaporating through Area_GI_Unit_Pervious_m2
	Storage_VaultLayer_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Vault_m3") / beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2");
	//Head_VaultLayer_m (m) initialized as Storage_VaultLayer_m (m) divided by Vault_Porosity_m3pm3 (m3/m3), increasing head due to solids
	//Note: Vault_Porosity_m3pm3 (m3/m3) has minimum value of Constant_1E_negative10, so no chance for division by zero
	Head_VaultLayer_m = Storage_VaultLayer_m / beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Porosity_m3pm3");


	//If Storage_GI_Vault_m3 (m3) is > Storage_GI_Vault_Max_m3 (m3) then layer is saturated and upper layers might also drain
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Vault_m3") >= beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Vault_Max_m3")) {
		//If Soil_Thickness_m (m) is > Constant_1E_negative10 then a soil layer exists above storage layer
		if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") > Constant_1E_negative10) {
			//If Storage_GI_Soil_m3 (m3) is > Storage_Soil_FieldCapacity_m3 (m3) then soil can drain
			if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_m3") > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Soil_FieldCapacity_m3")) {
				//Head_VaultLayer_m (m) is increased by Head_SoilLayer_m (m) drainage beyond Soil_FieldCapacity_m3pm3 (m3/m3)
				Head_VaultLayer_m = Head_VaultLayer_m + Head_SoilLayer_m;
				//If Storage_GI_Soil_m3 (m) is > Storage_GI_Soil_Max_m3 (m), then layer is saturated and upper layers might also drain
				if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_m3") >= beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_Max_m3")) {
					//If Pavement_Thickness_m (m) > Constant_1E_negative10 when Soil_Thickness_m (m) > Constant_1E_negative10 then get depth of pavement water above saturated soil
					if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") > Constant_1E_negative10) {
						//Head_VaultLayer_m (m) is increased by Head_PavementLayer_m (m)
						Head_VaultLayer_m = Head_VaultLayer_m + Head_PavementLayer_m;
						//If Storage_GI_Pavement_m3 (m3) > Storage_GI_Pavement_Max_m3 (m3), then layer is saturated and upper layers might also drain
						if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Pavement_m3") >= beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Pavement_Max_m3")) {
							//Head_VaultLayer_m (m) is increased by Head_SurfaceLayer_m (m)
							Head_VaultLayer_m = Head_VaultLayer_m + Head_SurfaceLayer_m;
						}
					}
					//Else Storage_GI_Soil_m3 (m) >= Storage_GI_Soil_Max_m3 (m) and Pavement_Thickness_m (m) equals 0, then get depth of surface water above saturated pavement
					else {
						//Head_VaultLayer_m (m) is increased by Head_SurfaceLayer_m (m)
						Head_VaultLayer_m = Head_VaultLayer_m + Head_SurfaceLayer_m;
					}
				}
			}
		}

		//If Pavement_Thickness_m (m) > Constant_1E_negative10, Soil_Thickness_m (m) <= 0 and Storage_GI_Vault_m3 >= Storage_GI_Vault_Max_m3, then add pavement head
		else if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") > Constant_1E_negative10 && beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") <= 0.0) {
			//Head_VaultLayer_m (m) is increased by Head_PavementLayer_m (m)
			Head_VaultLayer_m = Head_VaultLayer_m + Head_PavementLayer_m;
			//If Storage_GI_Pavement_m3 (m) > Pavement_Thickness_m (m), then get depth of surface water above saturated pavement
			if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Pavement_m3") >= beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Pavement_Max_m3")) {
				//Head_VaultLayer_m (m) is increased by Storage_GI_Surface_m3 (m)
				Head_VaultLayer_m = Head_VaultLayer_m + Head_SurfaceLayer_m;
			}
		}
		//If Storage_GI_Vault_m3 >= Storage_GI_Vault_Max_m3 AND
		//If Pavement_Thickness_m (m) and Soil_Thickness_m (m) <= 0, and Storage_GI_Vault_m3 >= Storage_GI_Vault_Max_m3
		//Else Storage Layer is directly below the Surface Layer, then add Head_SurfaceLayer_m
		else {
			//Head_VaultLayer_m (m) is increased by Head_PavementLayer_m (m)
			Head_VaultLayer_m = Head_VaultLayer_m + Head_SurfaceLayer_m;
		}
	}

	//Head_VaultLayer_m (m) is reduced by Vault_Outlet_Offset_m (m), to remove water below the pipe invert
	//Note: Vault_Outlet_Offset_m (m) is defined in HydroPlusConfig.xml as drain invert height above Storage Layer base
	Head_VaultLayer_m = Head_VaultLayer_m - beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Outlet_Offset_m");

	//If Head_VaultLayer_m (m) is > 0, then compute runoffStorageDrain_m3 (m3) from underdrain using EPA SWMM equation 
	if (Head_VaultLayer_m > 0) {
		//Runoff_VaultDrain_m (m) is product of Vault_Outlet_Discharge_Coeff and Head_VaultLayer_m to the power of Vault_Outlet_Discharge_Exponent
		//Note: Runoff equation is from EPA SWMM Version 5.1 lidproc.c, using SWMM GUI for Flow Coeff* variable (SI or BG units); Drain Advisor interface for guidance
		//Note: According to EPA SWMM, for SI uses -> Vault_Outlet_Discharge_Coeff (mm/hr) and Head_VaultLayer_m (m) with same Vault_Outlet_Discharge_Exponent
		//Note: According to EPA SWMM, for BG uses -> Vault_Outlet_Discharge_Coeff (in/hr) and Head_VaultLayer_ft (ft) with same Vault_Outlet_Discharge_Exponent
		//Note: EPA SWMM lidproc.c outflow = theLidProc->drain.coeff * pow(head, theLidProc->drain.expon), where head is units m based on function UCF(int u) and FlowUnits
		//Note: And there is no unit balance in the drain outflow equation used by EPA SWMM (and replicated here), making use of units such as mm/hr or in/hr meaningless
		//Note: And simulating Fig 14 of Platz et al. (2020), the Vault_Outlet_Discharge_Coeff was found to have same magnitude for SI and BG; if BG=0.35585, then SI=0.35585
		//Note: If Vault_Outlet_Discharge_Coeff = 0.35585 and Vault_Outlet_Discharge_Exponent = 1.06, then examine the Runoff_VaultDrain_m and Runoff_VaultDrain_ft ...
		//Note: ... when Head_VaultLayer_ft = [2, 1.5, 1.0], Runoff_VaultDrain_ft = [0.7436, 0.5476, 0.3558], then converting Head_VaultLayer_ft to Head_VaultLayer_m ...
		//Note: ... when Head_VaultLayer_m = [0.6096, 0.4572, 0.3048], Runoff_VaultDrain_m = [0.2105, 0.1552, 0.101], converted to Runoff_VaultDrain_m_to_ft = [0.6906, 0.5091, 0.3313]
		//Note: ... While there are notable differences in magnitude between Runoff_VaultDrain_m_to_ft and Runoff_VaultDrain_ft, they are close, linear, and highly correlated. ...
		//Note: ... Therefore, Vault_Outlet_Discharge_Coeff can be slightly adjusted to ensure Runoff_VaultDrain_m_to_ft equals Runoff_VaultDrain_ft.
		//Reference: Platz, M., M. Simon and M. Tryby (2020). "Testing of the Storm Water Management Model Low Impact Development Modules." JAWRA 56(2): 283-296, https://doi-org.esf.idm.oclc.org/10.1111/1752-1688.12832
		//Note: Consider refactor to optionally allow use of Vault_Outlet_ManningRoughness for vault outlet flow control, which will require slope and diameter of of pipe
		Runoff_VaultDrain_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Outlet_Discharge_Coeff") * pow(Head_VaultLayer_m, beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Outlet_Discharge_Exponent"));
		//Runoff_VaultDrain_m3 (m3) is product of Runoff_VaultDrain_ft (m) and Area_GI_Unit_m2 (m2) to get volume for folder area
		//Note: Vault is presumed to have area of Area_GI_Unit_m2 for exfiltration, with water otherwise infiltrating, percolating, and evaporating through Area_GI_Unit_Pervious_m2
		Runoff_VaultDrain_m3 = Runoff_VaultDrain_m * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2");
	}
	
	//Control for negative flows
	Runoff_VaultDrain_m3 = Inputs::forceZeroIfLessThanOrAlmostEqualZero(Runoff_VaultDrain_m3, Epsilon_Tolerance_1E_negative15);

	//Return Runoff_VaultDrain_m3 (m3) as folder volume
	return Runoff_VaultDrain_m3;
}
