#include "Inflow_StormwaterDevice.h"

#include "Runoff_StormwaterDevice.h"
#include "Groundwater_StormwaterDevice.h"
#include "Flux_Manager_StormwaterDevice.h"

//Green Infrastructure units may be comprised of 4 vertical storage layers 1) surface, 2) pavement, 3) soil and 4) vault.
//Note: The surface layer, number 1 above, is always at the top of the GI unit, and receives the precipitation and any runon
//Note: The sequence of layers below the surface is 2) pavement, 3) soil, and 4) vault, where the vault layer has porosity and a drain. 
//Note: The GI unit does not need to use all 3 layers below the surface layer, and the lowest layer exfiltrate to the catchment soils

//Transfers between GI and BulkArea and impacts on Catchment area (CA) water balance:
//Note: GI transfer of water to groundwater is through GI_Drainage_VadoseZone_m3, managed in StorageDeficit_VadoseZone::calculate function
//Note: GI transfer of water to runoff is through Runoff_Pervious_m3, Runoff_Impervious_m3, Runoff_Water_m3, and Runoff_Vault_m3, managed in RunoffSummation::calculate function 
//Note: CA groundwater presumed to be at equilibrium with GI device at start of simulation, and at end of every timeStep
//Note: CA groundwater potential releases via GI surface or vault drain not tallied in runoff, just as return flow is not tallied for BulkArea
//Note: CA groundwater level in GI with vault drain is unaffected by vault drain; groundwater would immediately replace any potentially drained
//Note: Variables are deliberately kept as beC (from HydroPlusConfig) or input->RepoDict, etc. 


//Similarities and differences between HydroPlus Green Infrastructure (GI) and EPA SWMM Low Impact Development (LID):
//Note: HydroPlus GI devices and functions were inspired by SWMM LID, and the initial HydroPlus GI code (e.g., Abdi PhD) took SWMM LID algorithms and connected them to TOPMODEL theory
//Note: HydroPlus GI algorithms are now different than the initial stage, having been rewritten (e.g., Endreny) to fit more cleanly into TOPMODEL theory and function within a water balance
//Note: Major differences include: a) use of one vertical flux sequence algorithm of all GI types; b) use of one vertical flux balance algorithm; c) use of porosity; d) initial saturation
//Note: GI variable for layer voids is porosity (phi=Vv/Vt, where Vv=volume voids, Vt=volume total). EPA SWMM used void ratios (e=Vv/Vs, where Vs=volume solids).
//Note: The use of porosity fits with computations. EPA SWMM had converted height of water, h, to layer depth, d, in SURF & STOR as if using porosity, d=h/e.
//Note: SWMM use of equation d=h/e is only correct if e=phi, which is mathematically incorrect when phi = e/(e+1) and e=phi/(1-phi).

//References:
//Chin, D. A. (2021). Water Resources Engineering, Fourth Edition. Hoboken, NJ: Pearson Education. ISBN-13: 9780135357750

//GI_InitializeAndLaunchWaterBalance function receives inflows and initializes the GI water balance
void Inflow_StormwaterDevice::GI_InitializeAndLaunchWaterBalance(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep,
	map<string, function<void(Inputs*, int)> >& functors)
{
	//Inflow_StormwaterDevice::ResetVariablesToZero function resets or initializes variables depending on time step
	Inflow_StormwaterDevice::ResetVariablesToZero(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
	
	//Flux_to_GI_Precipitation_m3 (m3) is 3 types of precipitation on GI, with no removal of precipitation by DCIA etc. 
	//Note: Flux_to_PerviousArea_Rain_SnowMelt_Irrigation_m (m) is only precipitation to GI pervious cover, multiplied by Area_GI_Unit_Pervious_m2 to get volume on GI (m3)
	//Note: Flux_to_ImperviousArea_Rain_SnowMelt_Irrigation_m (m) is only precipitation to impervious cover, multiplied by Area_GI_Unit_Impervious_m2 to get volume on GI (m3)
	//Note: Flux_to_WaterArea_Rain_SnowMelt_m (m) is only precipitation to water cover, multiplied by Area_GI_Unit_Water_m2 to get volume on GI (m3)
	//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
	beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_GI_Precipitation_m3") = (beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_PerviousArea_Rain_SnowMelt_Irrigation_m") * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2")) +
		(beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_ImperviousArea_Rain_SnowMelt_Irrigation_m") * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Impervious_m2")) +
		(beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_WaterArea_Rain_SnowMelt_m") * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Water_m2"));

	//Flux_BulkArea_to_GI_Runon_m3 (m3) is 3 types of runoff from Bulk Area to GI, divided by Count_GI_Units to get volume per GI unit
	//Note: Transfer water into GI from BulkArea; at timeStep = 0 typically zero, arriving after Flux_to_GI_Precipitation_m3 due to BulkArea delay
	//Note: Runoff_Impervious_BAtoGI_m3 (m3) was determined in DepressionImpervious.cpp
	//Note: Runoff_Water_BAtoGI_m3 (m3) was determined in DepressionWater.cpp
	//Note: Runoff_SatExcess_BAtoGI_m3 and Runoff_InfilExcess_BAtoGI_m3 (m3) was determined in Flux_LateralGroundWater.cpp
	//Note: All BulkArea transfers to GI presumed to go to GI pervious area, regardless of GI ImperviousCover_frac
	beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_BulkArea_to_GI_Runon_m3") = (input->RepoDict["Runoff_Impervious_BAtoGI_m3"] + input->RepoDict["Runoff_Water_BAtoGI_m3"] +
		input->RepoDict["Runoff_SatExcess_BAtoGI_m3"] + input->RepoDict["Runoff_InfilExcess_BAtoGI_m3"]) / beC->by_key(MapPixel_ID, DataFolder_ID, "Count_GI_Units");

	//Flux_GI_to_GI_Runon_m3 (m3) is defined by runon from upslope folder of GI devices if its Flag_Outflow_to_NextFolder = 1
	//Note: As done above, runon volume to GI is divided by Count_GI_Units to get volume per GI unit
	//Note: Flag_Outflow_to_NextFolder defined in HydroPlusConfig.xml file for GI device that will drain to next GI folder; undefined = 0
	beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_GI_to_GI_Runon_m3") = (input->RepoDict["Runoff_Impervious_GItoGI_m3"] + input->RepoDict["Runoff_Pervious_GItoGI_m3"] +
		input->RepoDict["Runoff_Vault_GItoGI_m3"] + input->RepoDict["Runoff_Water_GItoGI_m3"]) / beC->by_key(MapPixel_ID, DataFolder_ID, "Count_GI_Units");

	//Runoff_Impervious_to_Outlet_m3 (m3) is product of Flux_to_ImperviousArea_Rain_SnowMelt_Irrigation_m (m), impervious cover depth, Area_GI_Unit_Impervious_m2 and DCIA_frac
	//Note: Runoff_Impervious_to_Outlet_m3 (m3) allows for GI to generate runoff heading directly to outlet
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Impervious_to_Outlet_m3") = (beC->by_key(MapPixel_ID, DataFolder_ID, "Flux_to_ImperviousArea_Rain_SnowMelt_Irrigation_m") * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Impervious_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "DirectlyConnectedImperviousArea_frac"));

	//Call Flux_Manager_StormwaterDevice::Calculate function to simulate water movement in GI device
	Flux_Manager_StormwaterDevice::Calculate(input, beC, functors, MapPixel_ID, DataFolder_ID, timeStep);
}

//Inflow_StormwaterDevice::ResetVariablesToZero function resets or initializes variables depending on time step
void Inflow_StormwaterDevice::ResetVariablesToZero(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//Reset variables to zero at start of each time step
	beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_Surface_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Surface_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Evaporation_Surface_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "EvapoTrans_Pavement_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Percolation_Pavement_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "EvapoTrans_Soil_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Percolation_Soil_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Exfiltration_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "EvapoTrans_Vault_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Pervious_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Impervious_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Water_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Vault_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "GI_Drainage_VadoseZone_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Ratio_RunoffInfilExcess_to_RunoffPervious") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Runoff_Surface_BermCrest_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Groundwater_PavementLayer_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Groundwater_SoilLayer_m3") = 0.0;
	beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Groundwater_VaultLayer_m3") = 0.0;

	//Note: Consider refactoring so variables transferring water from BulkArea to GI could be reset to zero after transfer in Inflow_StormwaterDevice.cpp
	//Note: letting terms below keep values may be advantageous and still avoid summing them twice per timeStep in AggregateOutput 
	//input->RepoDict["Runoff_Impervious_BAtoGI_m3"] = 0.0;
	//input->RepoDict["Runoff_SatExcess_BAtoGI_m3"] = 0.0;
	//input->RepoDict["Runoff_InfilExcess_BAtoGI_m3"] = 0.0;

	//If timeStep equals first
	if (timeStep == 0.0) {
		//Time_PrecipitationEnded_hr, Time_RunoffStarted_hr, Time_PeriodSincePavementRegenerated_hr (hr) set to 0 day
		//Note: Could refactor HydroPlusConfig.xml to define these terms
		beC->by_key(MapPixel_ID, DataFolder_ID, "Time_PrecipitationEnded_hr") = 0 * Ratio_Hour_to_Day;
		beC->by_key(MapPixel_ID, DataFolder_ID, "Time_RunoffStarted_hr") = 0 * Ratio_Hour_to_Day;
		beC->by_key(MapPixel_ID, DataFolder_ID, "Time_PeriodSincePavementRegenerated_hr") = 0 * Ratio_Hour_to_Day;

		//Total_GI_Drainage_VadoseZone_m3 (m3) is volume passing from vadose zone to catchment groundwater
		beC->by_key(MapPixel_ID, DataFolder_ID, "Total_GI_Drainage_VadoseZone_m3") = 0.0;
		beC->by_key(MapPixel_ID, DataFolder_ID, "Total_Drainage_GI_Pavement_m3") = 0.0;

		//Initialize GI storage layers to zero before updating to values given in HydroPlusConfig.xml at start of first time step
		//Storage_GI_Surface_m3 (m3) is surface layer of GI unit; all GI units have this layer to receive runon and or precipitation
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Surface_m3") = 0.0;
		//Storage_GI_Pavement_m3 (m3) is pavement layer of GI unit; if it exits, then this layer is always directly below the surface layer
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Pavement_m3") = 0.0;
		//Storage_GI_Soil_m3 (m3) is soil layer of GI unit; if it exits, then this layer is always directly below the pavement or surface layer 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Soil_m3") = 0.0;
		//Storage_GI_Vault_m3 (m3) is vault layer of GI unit; if it exits, then this layer is always directly below the soil, pavement or surface layer 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Vault_m3") = 0.0;
		//Storage_GroundwaterMound_m3 is saturated mounding below GI device that drains to catchment groundwater table
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GroundwaterMound_m3") = 0;

		//Area_GI_Unit_m2 (m2) is Area_m2 divided by Count_GI_Units
		beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_m2") = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2") / beC->by_key(MapPixel_ID, DataFolder_ID, "Count_GI_Units");
		//Area_GI_Unit_Pervious_m2 (m2) is Area_per_GIUnit_m2 * 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
		beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2") = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "PerviousCover_frac");
		//Area_GI_Unit_Impervious_m2 (m2) is Area_per_GIUnit_m2 * PerviousCover_frac
		beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Impervious_m2") = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousCover_frac");
		//Area_GI_Unit_Water_m2 (m2) is Area_per_GIUnit_m2 * WaterCover_noTreeCanopy_frac
		beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Water_m2") = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "WaterCover_noTreeCanopy_frac");

		//Porosity_ThetaSat_ThetaWp_m3pm3 is the pore space between Soil_SaturationPoint_m3pm3 and Soil_WiltingPoint_m3pm3
		beC->by_key(MapPixel_ID, DataFolder_ID, "Porosity_ThetaSat_ThetaWp_m3pm3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_SaturationPoint_m3pm3") -
			beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_WiltingPoint_m3pm3");
		//Porosity_ThetaSat_ThetaFc_m3pm3 is the pore space between Soil_SaturationPoint_m3pm3 and Soil_FieldCapacity_m3pm3
		beC->by_key(MapPixel_ID, DataFolder_ID, "Porosity_ThetaSat_ThetaFc_m3pm3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_SaturationPoint_m3pm3") -
			beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_FieldCapacity_m3pm3");
		//Porosity_ThetaFc_ThetaWp_m3pm3 is the pore space between Soil_FieldCapacity_m3pm3 and Soil_WiltingPoint_m3pm3
		beC->by_key(MapPixel_ID, DataFolder_ID, "Porosity_ThetaFc_ThetaWp_m3pm3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_FieldCapacity_m3pm3") -
			beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_WiltingPoint_m3pm3");

		//Note: Some beC variables renamed beC->by_key variables to allow for variable redefinition when zero and avoid division by zero
		beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Porosity_m3pm3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Porosity_m3pm3");
		beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Porosity_m3pm3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Porosity_m3pm3");
		beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Thickness_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Thickness_m");
		beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Thickness_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Thickness_m");
		beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Thickness_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Thickness_m");
		beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_StorageInitial_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_StorageInitial_m");
		beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Porosity_m3pm3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Porosity_m3pm3");

		//If key variables = zero, then avoid subsequent NaN problems with quotients by redefining variables = Constant_1E_negative10
		//Note: Location of if conditionals must precede use of their terms, such as:
		// e.g., Volume_GI_SoilLayer_m3 = Soil_Thickness_m ... And Storage_GI_Soil_Max_m3 = Porosity_ThetaSat_ThetaWp_m3pm3 ...
		//Note: Code for Inflow_StormwaterDevice could be written with if conditionals to avoid division by zero, but this simpler approach should work
		//Note: Zero values for other key variables are handled in BuildDataOrganizer::Folder_CoverData_Statistical function ...
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2") <= 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2") = Constant_1E_negative10; }
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Porosity_m3pm3") <= 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Porosity_m3pm3") = Constant_1E_negative10; }
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Thickness_m") <= 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Thickness_m") = Constant_1E_negative10; }
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Porosity_m3pm3") <= 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Porosity_m3pm3") = Constant_1E_negative10; }
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Thickness_m") <= 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Thickness_m") = Constant_1E_negative10; }
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Porosity_ThetaSat_ThetaWp_m3pm3") <= 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Porosity_ThetaSat_ThetaWp_m3pm3") = Constant_1E_negative10; }
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Thickness_m") <= 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Thickness_m") = Constant_1E_negative10; }
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Porosity_m3pm3") <= 0.0) { beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Porosity_m3pm3") = Constant_1E_negative10; }


		//Surface_TopWidth_Max_m (m) = Surface_BottomWidth_m (m) + 2 * Surface_SideSlope_HtoV_mpm * Surface_Berm_Height_m (m)
		//Note: Top width equation for trapezoidal channel from Fig 5.2 in Chin (2021); Calculates maxium top width, which may be > active top width 
		//Note: Equation will work for rectangular (Surface_SideSlope_HtoV_mpm=0) or triangular (Surface_BottomWidth_m=0) if variable = zero
		beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_TopWidth_Max_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_BottomWidth_m")+ 2 * beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm") * beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Berm_Height_m");
		//Area_GI_CrossSection_Max_m2 (m2) = (Surface_BottomWidth_m (m) * Surface_Berm_Height_m (m) + Surface_SideSlope_HtoV_mpm * pow(Surface_Berm_Height_m,2)) * Surface_Porosity_m3pm3 (m3/m3)
		//Note: Cross section equation for trapezoidal channel from Eq 5.3 in Chin (2021); Calculates maxium cross section, which may be > active cross section ...
		//Note: ... Modifies the equation by removing area occupied by vegetation or other objects, using Surface_Porosity_m3pm3 which is fraction open to water
		//Note: Equation will work for rectangular (Surface_SideSlope_HtoV_mpm=0) or triangular (Surface_BottomWidth_m=0) if variable = zero
		beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_CrossSection_Max_m2") = (beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_BottomWidth_m")* beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Berm_Height_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_SideSlope_HtoV_mpm") * pow(beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Berm_Height_m"), 2)) * beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Porosity_m3pm3");
		//Surface_Length_GI_Unit_m (m) = Area_GI_Unit_m2 (m2) / Surface_TopWidth_Max_m (m) based on EPA SWMM and logic of Swale type GI
		beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Length_GI_Unit_m") = input->SafeDivide(beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_m2"), beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_TopWidth_Max_m"));

		//Storage_GI_Surface_Max_m3 (m3) = Area_GI_CrossSection_Max_m2 (m) * Area_GI_CrossSection_Max_m2
		//Note: Area_GI_CrossSection_Max_m2 (m2) represents a reduced area when Surface_Porosity_m3pm3 (m3/m3) < 1, 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Surface_Max_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_CrossSection_Max_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Length_GI_Unit_m");

		//Storage_GI_SurfaceDepression_Max_m3 (m3) is sum of ((ImperviousDepressionStorage_mm * Ratio_m_to_mm * Area_GI_Unit_Impervious_m2) + ... 
		//Note: ... (WaterDepressionStorage_mm * Ratio_m_to_mm * Area_GI_Unit_Water_m2) + ...
		//Note: (PerviousDepressionStorage_mm * Ratio_m_to_mm * Area_GI_Unit_Pervious_m2)) * Surface_Porosity_m3pm3 (m3/m3)
		//Note: Surface depression storage scaled by Surface_Porosity_m3pm3 (due to plants, hardscape) in same was as surface storage
		//Note: Storage_Depression_Max_m3 (m3) is volume in depressions below grade in the surface layer, not part of available runoff
		//Note: EPA SWMM lidproc.c swaleFluxRates (Build 5.2.4) uses dStore = 0, but treats it as within the active channel and ...
		//Note: ... resizes flowArea by removing dStore, and computes hydRadius (wetted perimeter) by using dStore topwidth as channel bottom width 
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_SurfaceDepression_Max_m3") = ((beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousDepressionStorage_mm") * Ratio_m_to_mm * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Impervious_m2")) + (beC->by_key(MapPixel_ID, DataFolder_ID, "WaterDepressionStorage_mm") * Ratio_m_to_mm * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Water_m2")) + (beC->by_key(MapPixel_ID, DataFolder_ID, "PerviousDepressionStorage_mm") * Ratio_m_to_mm * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2"))) * beC->by_key(MapPixel_ID, DataFolder_ID, "Surface_Porosity_m3pm3");

		//Storage_GI_Soil_m3 (m3) is product of (Soil_MoistureInitial_m3pm3 - Soil_WiltingPoint_m3pm3), Soil_Thickness_m, Area_per_GIUnit_m2 and PerviousCover_frac
		//Note: Initial soil moisture water scaled to wilting point to represent gravitational and plant available water
		//Note: Soil water is scaled to wilting point to represent gravitational and plant available water, with a minimum = 0.
		//Note:	EPA SWMM Version 5.1 scales similarly, e.g., lidUnit->initSat* (LidProcs[k].soil.porosity - LidProcs[k].soil.wiltPoint).
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Soil_m3") = (beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_MoistureInitial_m3pm3") - beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_WiltingPoint_m3pm3")) *
			beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Thickness_m") * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");
		//Storage_GI_Soil_TS_0th_m3 (m3) recorded at start of simulation, product of Storage_GI_Soil_m3 (m3) and Count_GI_Units to get GI folder Area volume
		//Note: Storage_GI_Soil_m3 (m3) at end of 1st time step may misrepresent the initial state of the GI system
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Soil_TS_0th_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Soil_m3") * beC->by_key(MapPixel_ID, DataFolder_ID, "Count_GI_Units");
		//Storage_Soil_FieldCapacity_m3 (m3) is product of Porosity_ThetaFc_ThetaWp_m3pm3, Soil_Thickness_m (m) and Area_GI_Unit_Pervious_m2 (m2)
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Soil_FieldCapacity_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Porosity_ThetaFc_ThetaWp_m3pm3") * beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Thickness_m") *
			beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");
		//Volume_GI_SoilLayer_m3 (m3) is product of Soil_Thickness_m (m) and Area_GI_Unit_Pervious_m2 (m2)
		//Note: Volume_GI_SoilLayer_m3 (m3) is volume of voids and solids, the entire layer, used for deriving volumetric moisture
		beC->by_key(MapPixel_ID, DataFolder_ID, "Volume_GI_SoilLayer_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Thickness_m") * beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");

		//Area_GI_Unit_Functional_m2 (m2) initiated to zero outside of if else statement
		double Area_GI_Unit_Functional_m2 = 0;
		//If Type is RainBarrel or RoofDisconnect set the functional area to GI footprint
		if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "RainBarrel" || beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "RoofDisconnect") {
			//Area_GI_Unit_Functional_m2 (m2) is redefined to Area_GI_Unit_m2 (m2)
			Area_GI_Unit_Functional_m2 = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_m2");
		}
		//Else Type is not RainBarrel or RoofDisconnect or vault has thickness <= Constant_1E_negative10 then
		else {
			//Area_GI_Unit_Functional_m2 (m2) is defaulted to Area_GI_Unit_Pervious_m2 (m2)
			Area_GI_Unit_Functional_m2 = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2");
		}
		//Area_GI_Unit_Vault_Functional_m2 (m2) presumed to have area of Area_GI_Unit_m2 for all GI units with such storage
		double Area_GI_Unit_Vault_Functional_m2 = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_GI_Unit_m2");

		//Storage_GI_Pavement_Max_Potential_m3 (m3) is product of Pavement_Porosity_m3pm3, Pavement_Thickness_m and Area_GI_Unit_Pervious_m2
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Pavement_Max_Potential_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Porosity_m3pm3") * beC->by_key(MapPixel_ID, DataFolder_ID, "Pavement_Thickness_m") *
			Area_GI_Unit_Functional_m2;
		//Storage_GI_Soil_Max_Potential_m3 (m3) is product of Porosity_ThetaSat_ThetaWp_m3pm3, Soil_Thickness_m and Area_GI_Unit_Pervious_m2
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Soil_Max_Potential_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Porosity_ThetaSat_ThetaWp_m3pm3") * beC->by_key(MapPixel_ID, DataFolder_ID, "Soil_Thickness_m") *
			Area_GI_Unit_Functional_m2;
		//Storage_GI_Vault_Max_Potential_m3 (m3) is product of Vault_Porosity_m3pm3, Vault_Thickness_m and Area_GI_Unit_Pervious_m2
		//Note: Vault is presumed to have area of Area_GI_Unit_m2, with water entering from (and evaporating) above through Area_GI_Unit_Pervious_m2
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Vault_Max_Potential_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Porosity_m3pm3") * beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Thickness_m") *
			Area_GI_Unit_Vault_Functional_m2;
		//Threshold_VaultDrain_m3 (m3) is the threshold above which vault drain will release water; pervious folder area is dedicated to storage
		//Note: Vault is presumed to have area of Area_GI_Unit_m2, with water entering from (and evaporating) above through Area_GI_Unit_Pervious_m2
		beC->by_key(MapPixel_ID, DataFolder_ID, "Threshold_VaultDrain_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Outlet_Offset_m") * beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Porosity_m3pm3") *
			Area_GI_Unit_Vault_Functional_m2;
		//Storage_GI_Vault_m3 (m3) equals product ot Storage_GI_Vault_Max_Potential_m3 (m3) and ratio of Vault_StorageInitial_m (m) to Vault_Thickness_m (m)
		//Note: Set at start of first time step, then responds to water balance dynamics during simulation
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Vault_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Vault_Max_Potential_m3") *
			beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_StorageInitial_m") / beC->by_key(MapPixel_ID, DataFolder_ID, "Vault_Thickness_m");

		//Initialize actual maximum storage to potential storage, noting for devices without exposure to GW this will not change 
		//Storage_GI_Pavement_Max_m3 (m3) initialized to Storage_GI_Pavement_Max_Potential_m3
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Pavement_Max_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Pavement_Max_Potential_m3");
		//Storage_GI_Soil_Max_m3 (m3) is initialized to Storage_GI_Soil_Max_Potential_m3
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Soil_Max_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Soil_Max_Potential_m3");
		//Storage_GI_Vault_Max_m3 (m3) is initialized to Storage_GI_Vault_Max_Potential_m3
		beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Vault_Max_m3") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_GI_Vault_Max_Potential_m3");

	}
}