#include "DrainageToEvapotranspirationZone.h"

//Soil water is represented as three components, with average soil moisture deficit representing total storage between soil surface and groundwater table: 
//Note: 1) evapotranspiration zone with capillarity held water for ET, releasing gravitational water into unsaturated zone; it can saturate from below
//Note: 2) vadose zone or AveSMD representing gravitational water, draining into saturated zone; it can saturate from below into evapotranspiration zone
//Note: 3) saturated zone for subsurface runoff; when it extends into evapotranspiration zone it fills the gravitational voids
//Note: The evapotranspiration and vadose zones are conceptually situated alongside each other, each filling a different volumetric component
//Note: The evapotranspiration zone storage depth extends to the limit of evapotranspiration accessing capillarity held water, a constant for all TI bins
//Note: The vadose zone extends from the groundwater table to the ground, and accounts for all gravitational held water, including evapotranspiration zone
//Note: ... StorageDeficit_FieldCapacity_SoilEvapZone_Max_m (m) = Evapotranspiration_Depth_m * (Soil_FieldCapacity_m3pm3 - Soil_WiltingPoint_m3pm3)
//Note: ... StorageDeficit_VadoseZone_m (m) or AveSMD (m) = (Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3) * Depth_SoilColumn_Ave_m

void  DrainageToEvapotranspirationZone::calculate(Inputs *input, CompactRagged* beC, int TI_bin, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//ia = TI_bin as current TI bin to get appropriate topographic index properties
	int ia = TI_bin;

	//Initialize vectors at start of simulation 
	if (timeStep == 0) {
		//StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) defined as StorageDeficit_FieldCapacity_SoilEvapZone_TS_first_m (m) for 1st time step
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TS_first_m is defined in BuildDataOrganizer as Maximum of ...
		//Note: ... zero and [Evapotranspiration_Depth_m * (Soil_FieldCapacity_m3pm3 - Soil_MoistureInitial_m3pm3)]
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m can extend beneath the groundwater table, as it represents capillary pore storage, ...
		//Note: ... while StorageDeficit_VadoseZone_TI_m represents gravitational pore storage
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TS_first_m is constant across all TI bins ...
		//Note: ... regardless of TI bin variation in depth to groundwater, e.g., StorageDeficit_VadoseZone_TI_m / Porosity_Drainable_m3pm3
		//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m should be initialized in Flux_LateralGroundWater::calculate with 0 first,
		// ...then filled with StorageDeficit_FieldCapacity_SoilEvapZone_TS_first_m here, at the beginning of the simulation of each bulkarea folder
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia]= beC->by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TS_first_m");
	}

	//Vectors of infiltration and drainage for each TI bin are set to their initial scalar value equivalent; this is repeated for each pass through the TI loop
	//Infiltration_TI_m (m) is defined as sum of Infiltration_viaInfilEx_PreRunoffSatExcess_m and Infiltration_viaSatEx_PreRunoffSatExcess_m for this TI bin
	//Note: Infiltration_viaInfilEx_PreRunoffSatExcess_m (m) and Infiltration_viaSatEx_PreRunoffSatExcess_m (m) hold scalar length from start of TI loop
	beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Infiltration_TI_m")[ia] = beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaInfilEx_PreRunoffSatExcess_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Infiltration_viaSatEx_PreRunoffSatExcess_m");
	//Drainage_macropore_TI_m (m) is defined to equal Drainage_macroPore_PreRunoffSatExcess_m for this TI bin, the latter holding the scalar value at start of TI loop
	beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_macropore_TI_m")[ia] = beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_macroPore_PreRunoffSatExcess_m");

	
	//StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) is reduced by Infiltration_TI_m (m) into evapotranspiration zone from infiltration
	//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m increases due to EvapoTranspiration_SoilEvapZone_TI_m in RootZoneEvapoTranspirationCalc.cpp
	//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m is a deficit, aka capacity; all TI bins have same thickness
	beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] - beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Infiltration_TI_m")[ia];

	//Drainage routine, combining macropore flow with gravitational water within soil evapotranspiration zone 
	//Drainage_SoilEvapZone_TI_m (m) is initialized to Drainage_macroPore_m (m), and is gravitational drainage from evapotranspiration zone to unsaturated zone
	//Drainage_macroPore_m (m) is the water from Flux_to_Infiltration_Drainage_m that bypasses infiltration and evapotranspiration zone, going directly to subsurface flow
	beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_macropore_TI_m")[ia];

	//If StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) < 0 there is gravitational water beyond field capacity, and this drains to the unsaturated zone
	//Note: deficits < 0 exceed maximum storage capacity; deficits = 0 are at maximum storage capacity which is capillarity held water for ET; 
	//Note: StorageDeficit_FieldCapacity_SoilEvapZone_TI_m (m) capacity is only capillarity held water, and Storage_VadoseZone_TI_m (m) capacity is only gravitational water
	if (Inputs::isLessThan(beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia], 0.0)) {
		//Drainage_SoilEvapZone_TI_m is increased by excess water in StorageDeficit_FieldCapacity_SoilEvapZone_TI_m, where subtracting a negative is equivalent to adding StorageDeficit_FieldCapacity_SoilEvapZone_TI_m to Drainage_SoilEvapZone_TI_m
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia] = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia] - beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia];

		//StorageDeficit_FieldCapacity_SoilEvapZone_TI_m is reset to 0
		beC->vec_by_key(MapPixel_ID, DataFolder_ID, "StorageDeficit_FieldCapacity_SoilEvapZone_TI_m")[ia] = 0.0;
	}

	//Drainage_SoilEvapZone_PreRunoffSatExcess_m (m) scalar is computed in TI loop as weighted sum of Drainage_SoilEvapZone_TI_m (m) and TI_binArea_frac (fraction), prior to any updates due to Runoff_SatExcess_TI_m
	beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_PreRunoffSatExcess_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_PreRunoffSatExcess_m") + beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Drainage_SoilEvapZone_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac");
	//Drainage_SoilEvapZone_m (m) scalar will be updated in VadoseZoneDrainageCalc.cpp after any needed updates to Drainage_SoilEvapZone_TI_m due to Runoff_SatExcess_TI_m
}