#include "Evaporation_StormwaterDevice.h"
#include "../Evaporation/RootZoneEvapoTranspirationCalc.h"
#include "../Evaporation/OpenWaterEvaporationCalc.h"

//GI_RateOfEvaporation function references ET rates computed externally, then scales to proper area; now used only for SURF and SOIL layers
void Evaporation_StormwaterDevice::GI_RateOfEvaporation(Inputs* input, CompactRagged* beC,
				map<string, function<void(Inputs*, int)> >& functors, int DataDrawer_ID, int DataFolder_ID, int timeStep) {

	double Storage_SurfaceLayer_m = 0.0;
	double Storage_PavementLayer_m = 0.0;
	double Storage_SoilLayer_m = 0.0;
	double Storage_VaultLayer_m = 0.0;
	double Depth_To_PavementWater_m = 0.0;
	double Depth_To_SoilWater_m = 0.0;
	double Depth_To_StorageWater_m = 0.0;
	double Storage_PavementLayer_In_EvapoTransZone_m = 0.0;
	double Storage_SoilLayer_In_EvapoTransZone_m = 0.0;
	double Storage_VaultLayer_In_EvapoTransZone_m = 0.0;
	double Storage_EvapoTransZone_PreEvapoTrans_m = 0.0;

	//Note: Multiplication of variables by fractional area below is to adjust depths from one analysis area to entire folder area
	//Rationale: 1. variables arrive as depths for an analysis area (e.g., pervious area)
	//Rationale: 2. multiplying these depths by their analysis area gives a volume, which is conserved when transfered between different areas
	//Rationale: 3. multiplying these depths by their fractional area is equivalent to multiplying by their analysis area and dividing by the folder area, getting depth for folder area
	//Note: Multiplication is equivalent to multiplying by ratio of analysis area to folder area to convert from /analysis_area to /folder_area

	//Note: Division of variables by fractional area below is to adjust depths from entire folder area to one analysis area
	//Rationale: 1. variables arrive as depths for the entire folder area
	//Rationale: 2. multiplying these depths by the folder area gives a volume
	//Rationale: 3. dividing these depths by the analysis area is equivalent to multiplying by the folder area and dividing by the analysis area, getting depth for analysis area
	//Note: Division is equivalent to multiplying by ratio of folder area to analysis area to convert from /folder_area to /analysis_area

	//If Storage_GI_Surface_Max_m3 > 0 then there is a berm on the device, and 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) {
		//Storage_SurfaceLayer_m (m) is Surface_Berm_Height_m (m) * (Storage_GI_Surface_m3 / Storage_GI_Surface_Max_m3) 
		//Note: Surface_Berm_Height_m (m) represents maximum possible depth relative to ratio (Storage_GI_Surface_m3 / Storage_GI_Surface_Max_m3) 
		//Note: This is a conservative estimate of depth, to ensure it is not too high and misrepresent water available for surface evaporation
		Storage_SurfaceLayer_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Surface_Berm_Height_m") * input->SafeDivide(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3"), beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_Max_m3"));
	}
	//Else Storage_GI_Surface_Max_m3 <= 0 and there is no berm in the GI, so no scaling needed for incremental saturation
	else {
		//Storage_SurfaceLayer_m (m) is Storage_GI_Surface_m3 (m3) / Area_GI_Unit_m2 (m2)
		Storage_SurfaceLayer_m = input->SafeDivide(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Surface_m3"), beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2"));
	}

	//Note: Storage_PavementLayer_m (m) is quotient of (Storage_GI_Pavement_m3 and Area_per_GIUnit_m2), divided by PerviousCover_frac
	Storage_PavementLayer_m = input->SafeDivide(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Pavement_m3"), beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2"));
	//Note: Storage_SoilLayer_m (m) is quotient of (Storage_GI_Soil_m3 and Area_per_GIUnit_m2), divided by PerviousCover_frac
	Storage_SoilLayer_m = input->SafeDivide(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Soil_m3"), beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2"));
	//Note: Storage_VaultLayer_m (m) is quotient of (Storage_GI_Vault_m3 and Area_per_GI_Unit_m2), not limiting vault to below pervious 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
	Storage_VaultLayer_m = input->SafeDivide(beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_GI_Vault_m3"), beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_m2"));

	//calculatePerDepEvap functor computes Evaporation_PerviousPondedWater_m (m) for pervious cover fraction
	functors["calculatePerDepEvap"](input, timeStep);
	//calculateImpDepEvap functor computes Evaporation_ImperviousPondedWater_m (m) for pervious cover fraction
	functors["calculateImpDepEvap"](input, timeStep);
	//calculateWaterDepEvap functor computes Evaporation_WaterPondedWater_m (m) for water cover fraction
	functors["calculateWaterDepEvap"](input, timeStep);

	//Storage_PerviousPondedWater_m (m) is Storage_SurfaceLayer_m (m) divided by PerviousCover_frac to get depth over pervious area
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_PerviousPondedWater_m") = Storage_SurfaceLayer_m * beC->by_key(DataDrawer_ID, DataFolder_ID, "PerviousCover_frac");
	//Storage_ImperviousPondedWater_m (m) is Storage_SurfaceLayer_m (m) divided by ImperviousCover_frac to get depth over impervious area
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_ImperviousPondedWater_m") = Storage_SurfaceLayer_m * beC->by_key(DataDrawer_ID, DataFolder_ID, "ImperviousCover_frac");
	//Storage_WaterPondedWater_m (m) is Storage_SurfaceLayer_m (m) divided by WaterCover_noTreeCanopy_frac to get depth over water area
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_WaterPondedWater_m") = Storage_SurfaceLayer_m * beC->by_key(DataDrawer_ID, DataFolder_ID, "WaterCover_noTreeCanopy_frac");

	//If Evaporation_PerviousPondedWater_m (m) > Storage_PerviousPondedWater_m (m) then set to Storage_PerviousPondedWater_m (m)
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_PerviousPondedWater_m")) {
		//Evaporation_PerviousPondedWater_m (m) is set to Storage_PerviousPondedWater_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_PerviousPondedWater_m");
	}
	//If Storage_Snow_noCanopy_m (m) > 0, there is snow on ground and liquid evaporation from storage is zero for this timestep
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") > 0.0) {
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") = 0.0;
	}
	//If Evaporation_ImperviousPondedWater_m (m) > Storage_ImperviousPondedWater_m (m) then set to Storage_ImperviousPondedWater_m (m)
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_ImperviousPondedWater_m")) {
		//Evaporation_ImperviousPondedWater_m (m) is set to Storage_ImperviousPondedWater_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_ImperviousPondedWater_m");
	}
	//If Storage_Snow_noCanopy_m (m) > 0, there is snow on ground and liquid evaporation from storage is zero for this timestep
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") > 0.0) {
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") = 0.0;
	}
	//If Evaporation_WaterPondedWater_m (m) > Storage_WaterPondedWater_m (m) then set to Storage_ImperviousPondedWater_m (m)
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") > beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_WaterPondedWater_m")) {
		//Evaporation_WaterPondedWater_m (m) is set to Storage_WaterPondedWater_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_WaterPondedWater_m");
	}
	//If Storage_Snow_noCanopy_m (m) > 0, there is snow on ground and liquid evaporation from storage is zero for this timestep
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_Snow_noCanopy_m") > 0.0) {
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") = 0.0;
	}

	//Evaporation_Surface_m3 (m3) sums Evaporation_PerviousPondedWater_m and Evaporation_ImperviousPondedWater_m as volumes for folder area
	//Note: Evaporation_PerviousPondedWater_m (m) for pervious cover  converted to volume by with Area_GI_Unit_Pervious_m2
	//Note: Evaporation_ImperviousPondedWater_m (m) for impervious cover converted to volume with Area_GI_Unit_Impervious_m2
	//Note: Evaporation_WaterPondedWater_m (m) for water cover converted to volume with Area_GI_Unit_Water_m2
	beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Surface_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2") +
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Impervious_m2") +
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Water_m2");

	//Storage_EvapoTransZone_m (m) reset to zero before recalculating its value, using the water from distinct GI layers in this zone
	//Note: HydroPlusConfig.xml element Evapotranspiration_Depth_m for GI DataFolder is used in comparison with Storage_EvapoTransZone_m (m)
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = 0;
	//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) reset to StorageDeficit_FieldCapacity_SoilEvapZone_Max_m before subtracting GI layer water
	beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_Max_m");

	//Depth_To_PavementWater_m (m) is depth from ground to top of water surface in pavement layer
	//Note: Defined as (Pavement_Thickness_m - Storage_PavementLayer_m)
	Depth_To_PavementWater_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") - Storage_PavementLayer_m;
	//If Depth_To_PavementWater_m (m) and Pavement_Thickness_m (m) are less than Evapotranspiration_Depth_m (m) then
	if (Depth_To_PavementWater_m < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") &&
				beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m")) {
		//Storage_SoilLayer_In_EvapoTransZone_m  (m) is pavement water within soil evapotranspiration zone depth, (Evapotranspiration_Depth_m - Depth_To_PavementWater_m)
		Storage_PavementLayer_In_EvapoTransZone_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") - Depth_To_PavementWater_m;
		//Storage_EvapoTransZone_m (m) is increased to include Storage_PavementLayer_In_EvapoTransZone_m  
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") + Storage_PavementLayer_In_EvapoTransZone_m;
		//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) reduced by Storage_PavementLayer_In_EvapoTransZone_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") - Storage_PavementLayer_In_EvapoTransZone_m;
	}
	//Else If Depth_To_PavementWater_m (m) is less than Evapotranspiration_Depth_m (m) and Pavement_Thickness_m (m) is greater than Evapotranspiration_Depth_m (m) then
	else if (Depth_To_PavementWater_m < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") &&
				beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") > beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m")) {
		//Storage_SoilLayer_In_EvapoTransZone_m  (m) is pavement water within soil evapotranspiration zone depth, (Evapotranspiration_Depth_m - Depth_To_PavementWater_m)
		Storage_PavementLayer_In_EvapoTransZone_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") - Depth_To_PavementWater_m;
		//Storage_EvapoTransZone_m (m) is increased to include Storage_PavementLayer_In_EvapoTransZone_m  
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") + Storage_PavementLayer_In_EvapoTransZone_m;
		//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) reduced by Storage_PavementLayer_In_EvapoTransZone_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") - Storage_PavementLayer_In_EvapoTransZone_m;
	}

	//Depth_To_SoilWater_m (m) is depth from ground to top of water surface in soil layer;
	//Note: Defined as sum of potential layers Pavement_Thickness_m and (Soil_Thickness_m - Storage_GI_Soil_m3 * Porosity_ThetaSat_ThetaWp_m3pm3)
	Depth_To_SoilWater_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") + (beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") - Storage_SoilLayer_m);
	//If Depth_To_SoilWater_m (m) and (Pavement_Thickness_m + Soil_Thickness_m) (m) < Evapotranspiration_Depth_m (m) then
	if (Depth_To_SoilWater_m < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") &&
				(beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m")) < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m")) {
		//Storage_SoilLayer_In_EvapoTransZone_m  (m) is soil water within soil depth, (Pavement_Thickness_m + Soil_Thickness_m) - Depth_To_SoilWater_m
		Storage_SoilLayer_In_EvapoTransZone_m = (beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m")) - Depth_To_SoilWater_m;
		//Storage_EvapoTransZone_m (m) is increased to include Storage_SoilLayer_In_EvapoTransZone_m
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") + Storage_SoilLayer_In_EvapoTransZone_m;
		//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) reduced by Storage_PavementLayer_In_EvapoTransZone_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") - Storage_SoilLayer_In_EvapoTransZone_m;
	}
	//Else If Depth_To_SoilWater_m (m) is less than Evapotranspiration_Depth_m (m) and (Pavement_Thickness_m + Soil_Thickness_m) (m) > Evapotranspiration_Depth_m (m) then
	else if (Depth_To_SoilWater_m < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") &&
				(beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m")) > beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m")) {
		//Storage_SoilLayer_In_EvapoTransZone_m  (m) is soil water within soil evapotranspiration zone depth, (Evapotranspiration_Depth_m - Depth_To_SoilWater_m)
		Storage_SoilLayer_In_EvapoTransZone_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") - Depth_To_SoilWater_m;
		//Storage_EvapoTransZone_m (m) is increased to include Storage_SoilLayer_In_EvapoTransZone_m
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") + Storage_SoilLayer_In_EvapoTransZone_m;
		//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) reduced by Storage_PavementLayer_In_EvapoTransZone_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") - Storage_SoilLayer_In_EvapoTransZone_m;
	}

	//Depth_To_StorageWater_m (m) is depth from ground to top of water surface in storage layer;
	//Note: Defined as sum of potential layers Pavement_Thickness_m and Soil_Thickness_m, and (Vault_Thickness_m - Storage_GI_Vault_m3 * Vault_Porosity_m3pm3)
	Depth_To_StorageWater_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") + (beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Thickness_m") - Storage_VaultLayer_m);
	//If Depth_To_StorageWater_m (m) and (Pavement_Thickness_m + Soil_Thickness_m + Vault_Thickness_m) (m) < Evapotranspiration_Depth_m (m) then
	if (Depth_To_StorageWater_m < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") &&
				(beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Thickness_m")) < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m")) {
		//Storage_VaultLayer_In_EvapoTransZone_m (m) is storage water within storage depth, (Pavement_Thickness_m + Soil_Thickness_m + Vault_Thickness_m) - Depth_To_StorageWater_m
		Storage_VaultLayer_In_EvapoTransZone_m = (beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Thickness_m")) - Depth_To_StorageWater_m;
		//Storage_EvapoTransZone_m (m) is increased to include Storage_VaultLayer_In_EvapoTransZone_m  
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") + Storage_VaultLayer_In_EvapoTransZone_m;
		//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) reduced by Storage_PavementLayer_In_EvapoTransZone_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") - Storage_VaultLayer_In_EvapoTransZone_m;
	}
	//Else If Depth_To_StorageWater_m (m) is less than Evapotranspiration_Depth_m (m) and (Pavement_Thickness_m + Soil_Thickness_m + Vault_Thickness_m) (m) > Evapotranspiration_Depth_m (m) then
	else if (Depth_To_StorageWater_m < beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") &&
				(beC->by_key(DataDrawer_ID, DataFolder_ID, "Pavement_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Soil_Thickness_m") + beC->by_key(DataDrawer_ID, DataFolder_ID, "Vault_Thickness_m")) >
				beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m")) {
		//Storage_VaultLayer_In_EvapoTransZone_m (m) is storage water within soil evapotranspiration zone depth, (Evapotranspiration_Depth_m - Depth_To_StorageWater_m)
		Storage_VaultLayer_In_EvapoTransZone_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Evapotranspiration_Depth_m") - Depth_To_StorageWater_m;
		//Storage_EvapoTransZone_m (m) is increased to include Storage_VaultLayer_In_EvapoTransZone_m  
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") + Storage_VaultLayer_In_EvapoTransZone_m;
		//StorageDeficit_FieldCapacity_SoilEvapZone_m (m) reduced by Storage_PavementLayer_In_EvapoTransZone_m (m)
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") - Storage_VaultLayer_In_EvapoTransZone_m;
	}

	//Storage_EvapoTransZone_PreEvapoTrans_m (m) depth over folder area saved for later analysis
	Storage_EvapoTransZone_PreEvapoTrans_m = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m");

	//If PerviousCover_frac > 0 then divide by PerviousCover_frac
	//Note: Water balance above was performed relative to pervious area, starting with volume divided by pervious area, and now is rescaled to folder
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "PerviousCover_frac") > 0) {
		//Storage_EvapoTransZone_m adjusted from folder area to pervious area with division by PerviousCover_frac
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") / beC->by_key(DataDrawer_ID, DataFolder_ID, "PerviousCover_frac");
		//Storage_EvapoTransZone_m adjusted from folder area to pervious area with division by PerviousCover_frac
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") / beC->by_key(DataDrawer_ID, DataFolder_ID, "PerviousCover_frac");
	}
	//Else avoid division by zero and define as same value, which is likely zero
	else {
		beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m");
		beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m");
	}

	//Ensure Storage_EvapoTransZone_m (m) is >= zero
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") < 0) { beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_EvapoTransZone_m") = 0.0; }
	//Ensure StorageDeficit_FieldCapacity_SoilEvapZone_m (m) is >= zero; 
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") < 0) { beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m") = 0.0; }

	//Storage_SoilEvapZone_m (m) is computed as difference of maximum storage deficit StorageDeficit_FieldCapacity_SoilEvapZone_Max_m (m) and actual storage deficit StorageDeficit_FieldCapacity_SoilEvapZone_m (m)
	//Note: Storage_SoilEvapZone_m (m) updated each time step here for GI, and updated in Flux_LateralGroundwater.cpp for BulkArea
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_SoilEvapZone_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_Max_m") - beC->by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_m");

	//TI_bin as current TI bin to get appropriate topographic index properties
	int TI_bin = 0;
	//Calculate function computes EvapoTranspiration_SoilEvapZone_m (m) for pervious cover fraction
	//Note: Call to RootZoneEvapoTranspirationCalc::Calculate obtains EvapoTranspiration_SoilEvapZone_m, needed for soil water balance
	RootZoneEvapoTranspirationCalc::Calculate(input, beC, TI_bin, DataDrawer_ID,DataFolder_ID, timeStep);

	//If Storage_EvapoTransZone_PreEvapoTrans_m > 0 then divide by Storage_EvapoTransZone_PreEvapoTrans_m to scale ET into different GI layers
	if (Storage_EvapoTransZone_PreEvapoTrans_m > 0) {
		//EvapoTrans_Pavement_m3 (m) is EvapoTranspiration_SoilEvapZone_TI_m (m) over pervious cover scaled to pavement layer fraction by ..
		// product with Storage_PavementLayer_In_EvapoTransZone_m / Storage_EvapoTransZone_PreEvapoTrans_m and scaled to folder volume by ..
		// product with Area_per_GIUnit_m2 and PerviousCover_frac
		beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Pavement_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2") *
			(Storage_PavementLayer_In_EvapoTransZone_m / Storage_EvapoTransZone_PreEvapoTrans_m);

		//EvapoTrans_Soil_m3 (m) is EvapoTranspiration_SoilEvapZone_TI_m (m) over pervious cover scaled to soil layer fraction by ..
		// product with Storage_SoilLayer_In_EvapoTransZone_m / Storage_EvapoTransZone_PreEvapoTrans_m and scaled to folder volume by ..
		// product with Area_per_GIUnit_m2 and PerviousCover_frac
		beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Soil_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2") *
			(Storage_SoilLayer_In_EvapoTransZone_m / Storage_EvapoTransZone_PreEvapoTrans_m);

		//EvapoTrans_Vault_m3 (m) is EvapoTranspiration_SoilEvapZone_TI_m (m) over pervious cover scaled to storage layer fraction by ..
		// product with Storage_VaultLayer_In_EvapoTransZone_m / Storage_EvapoTransZone_PreEvapoTrans_m and scaled to folder volume by ..
		// product with Area_per_GIUnit_m2 and 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(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Vault_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") * beC->by_key(DataDrawer_ID, DataFolder_ID, "Area_GI_Unit_Pervious_m2") *
			(Storage_VaultLayer_In_EvapoTransZone_m / Storage_EvapoTransZone_PreEvapoTrans_m);
	}

	//Else Storage_EvapoTransZone_PreEvapoTrans_m = 0 and avoid division
	else {
		//EvapoTrans_Pavement_m3 (m) is water depth for folder area at value of zero due to Storage_EvapoTransZone_PreEvapoTrans_m = 0 
		beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Pavement_m3") = 0;
		//EvapoTrans_Soil_m3 (m) is water depth for folder area at zero due to Storage_EvapoTransZone_PreEvapoTrans_m = 0 
		beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Soil_m3") = 0;
		//EvapoTrans_Vault_m3 (m) is water depth for folder area at zero due to Storage_EvapoTransZone_PreEvapoTrans_m = 0 
		beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Vault_m3") = 0;
	}

	//Control for negative values
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Surface_m3") < 0) { beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Surface_m3") = 0.0; }
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Pavement_m3") < 0) { beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Pavement_m3") = 0.0; }
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Soil_m3") < 0) { beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Soil_m3") = 0.0; }
	if (beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Vault_m3") < 0) { beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Vault_m3") = 0.0; }

	//EvapoTranspiration_GI_m3 (m3) is sum of all evaporation and transpiration components
	beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_GI_m3") = beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Surface_m3") + beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Pavement_m3") + 
		beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Soil_m3") + beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTrans_Vault_m3");
}
