
#include "RunoffSurfaceSaturationExcess.h"

void RunoffSurfaceSaturationExcess::calculate(Inputs *input, CompactRagged* beC, int TI_bin, int DataDrawer_ID, int DataFolder_ID, int timeStep)
{
	//Saturation excess and variable source area computations based on algorithm taken from TOPMODEL code developed by Keith Beven and used by USGS D.M. Wolock https://pubs.er.usgs.gov/publication/wri934124
	//Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124).
	//Note: The algorithm requires the topographic index (TI) bins (ia) are sorted from highest (wettest) to lowest TI values
	//Note: The algorithm requires when ia = 0, the TI bin TI_Value[0] should be the upper limit and its TI_Area_frac[0] = 0
	//Note: The algorithm presumes for each ia, the area of catchment between discrete TI values TI_Value[ia] and TI_Value[ia+1] is TI_Area_frac[ia]
	//Note: The algorithm requires for each ia, the area of catchment bounding TI_Value[ia] is = to (TI_Area_frac[ia] + TI_Area_frac[ia+1])/2 

	//Original TOPMODEL weights discrete TI catchment area fractions of histogram with bounding TI parameter values
	//Note: TOPMODEL code used: OF=AC(IA)*(EX(IB)+EX(IA))/2, when EX[IA] values and EX[IB] had saturation excess, 
	//Note: where OF is saturation excess runoff, AC(IA) is area of catchment between bins IA and IA+1, 
	//Note: and EX(IB) and EX(IA) area saturation excess runoff in bins IA-1 and IA 
	//Note: TOPMODEL code used: OF = ACF * EX(IB) / 2, and ACF = ACF * EX(IB) / (EX(IB) - EX(IA)), when only EX[IB] had saturation excess runoff
	//Note: Hence, TOPMODEL used discrete TI area with bounding TI values, and apparently had logical mistakes, such as:
	//Note: 	1) presuming AC(IA) extends from EX(IB) to EX(IA), when actually AC(IA) extends from EX(IA) to EX(IA+1), where IB = IA-1, not IA+1;
	//Note: 	2) presuming EX(IA) is negative when only EX(IB) is positive, when actually EX(IA) equals 0, making the weight EX(IB) / (EX(IB) - EX(IA)) = 1
	//Note: HydroPlus weights discrete TI parameter values, i.e, SoilParameter[ia], with TI_binArea_frac, catchment area bounding discrete value

	//ia = TI_bin as current TI bin to get appropriate topographic index properties
	int ia = TI_bin;

	//Runoff_SatExcess_m (m) is the saturation excess runoff, reset to 0.0
	double Runoff_SatExcess_m = 0.0;
	//Groundwater_surficial_m (m) initialized for use in TI bin loop
	//Note: Groundwater_surficial_m (m) used to track surficial expression of variable source area
	double Groundwater_surficial_m = 0.0;
	//Groundwater_surficial_frac (frac) initialized for use in TI bin loop
	//Note: Groundwater_surficial_frac (frac) used to track surficial expression of variable source area
	double Groundwater_surficial_frac = 0.0;
	
	//if Runoff_SatExcess_TI_m[ia] > 0 then saturation excess runoff occurs in TI_binArea_frac, fraction of catchment bounding TI parameter value
	if (beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "Runoff_SatExcess_TI_m")[ia] > 0.0) {
		//Runoff_SatExcess_m = beC->vec_by_key(MapPixel_ID, DataFolder_ID, "Runoff_SatExcess_TI_m")[ia] * beC->by_key(MapPixel_ID, DataFolder_ID, "TI_binArea_frac")
		//Note: Runoff_SatExcess_TI_m (m) is determined in DrainageToSaturationZone::calculate function, after trying to infiltrate and percolate water
		Runoff_SatExcess_m = beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "Runoff_SatExcess_TI_m")[ia] * beC->by_key(DataDrawer_ID, DataFolder_ID, "TI_binArea_frac");
	}

	//if StorageDeficit_VadoseZone_TI_temp_m [ia] < 0 then groundwater breached surface occurs in TI_binArea_frac, fraction of catchment bounding TI parameter value
	if (beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_temp_m")[ia] < 0.0) {
		//Groundwater_surficial_m (m) is computed using Eq 21 of Wolock (1993) substituted into TOPMODEL code algorithm for saturation excess runoff
		//Note: HydroPlus weights discrete TI parameter values, i.e, SoilParameter[ia], with TI_binArea_frac, catchment area bounding discrete value
		Groundwater_surficial_m = abs(beC->vec_by_key(DataDrawer_ID, DataFolder_ID, "StorageDeficit_VadoseZone_TI_temp_m")[ia]) * beC->by_key(DataDrawer_ID, DataFolder_ID, "TI_binArea_frac");
		//Groundwater_surficial_frac (fraction) is fraction of watershed area with surficial groundwater
		Groundwater_surficial_frac = beC->by_key(DataDrawer_ID, DataFolder_ID, "TI_binArea_frac");
 	}

	//Runoff_SatExcess_m (m) increased by Runoff_SatExcess_m
	//Note: Accumulates saturation excess runoff due to effective precipitation on surficial groundwater table across all TI bins
	//Note: Consider Refactor to check if Runoff_SatExcess_TI_m * TI_binArea_frac fits within available pervious depression storage
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_SatExcess_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_SatExcess_m") + Runoff_SatExcess_m;

	//Groundwater_surficial_m (m) increased by Groundwater_surficial_m
	//Note: Accumulates upwelling groundwater depth across all TI bins
	//Note: Record this before resetting negative StorageDeficit_VadoseZone_TI_m to zero
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Groundwater_surficial_m") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Groundwater_surficial_m") + Groundwater_surficial_m;
	//Groundwater_surficial_frac (fraction) increased by Groundwater_surficial_frac
	//Note: Accumulates upwelling groundwater fraction of catchment across all TI bins
	//Note: Record this before resetting negative StorageDeficit_VadoseZone_TI_m to zero
	beC->by_key(DataDrawer_ID, DataFolder_ID, "Groundwater_surficial_frac") = beC->by_key(DataDrawer_ID, DataFolder_ID, "Groundwater_surficial_frac") + Groundwater_surficial_frac;
}
