#include "SimulationCoordinator.h"

#include "Interception/TreeInterceptionCalc.h"
#include "Interception/ShortVegInterceptionCalc.h"
#include "Evaporation/TreeEvaporationCalc.h"
#include "Evaporation/ShortVegEvaporationCalc.h"
#include "Evaporation/OpenWaterEvaporationCalc.h"
#include "Evaporation/PotentialEvaporationCalc.h"
#include "SnowOnGround/SnowAblationOpenAreaCalc.h"
#include "SnowOnGround/SnowAblationUnderShortVegCalc.h"
#include "SnowOnGround/SnowAblationUnderTreeCalc.h"
#include "InflowsToSurface/InflowsToSurfaceCalc.h"
#include "StorageSurface/DepressionImperviousCalc.h"
#include "StorageSurface/DepressionPerviousCalc.h"
#include "StorageSurface/DepressionWaterCalc.h"
#include "Infiltration/InfiltrationExponentialDecay.h"
#include "Infiltration/InfiltrationPowerDecay.h"
#include "LateralDistribution/Flux_LateralGroundWater.h"
#include "StorageSubsurface/StorageDeficit_VadoseZone.h"
#include "Runoff/RunoffSubsurface.h"
#include "Runoff/RunoffSummation.h"
#include "Aggregation/AggregateOutput.h"
#include "Runoff/RoutingDiffusion.h"
#include "WaterQuality/Loading_WaterQuality_EMC.h"
#include "WaterQuality/Removal_StormwaterPollutants.h"
#include "CoolAir/HeatFluxCal.h"
#include "CoolAir/BlockGroupHandler.h"
#include "CoolAir/HeatMetrics_Calc.h"
#include "GreenInfrastructure/Inflow_StormwaterDevice.h"
#include "CoolRiver/SolarCalculation.h"
#include "CoolRiver/StreamTemperature.h"
#include "Buffer/BufferDynamicCalc.h"
#include "Buffer/BufferSpatialCalc.h"
#include "CoolBuilding/EnergyLoad_Calc.h"
#include "Outputs/RoutedRunoff_TimeSeries.h"
#include "Outputs/TemperatureOutputWriter.h"
#include "Outputs/CoolRiverOutputWriter.h"
#include "Outputs/WriteOutputs.h"
#include "Inputs/WeatherProcessor.h"
#include <algorithm>
#include <limits>
#include <cstddef>
#include <string>
#include <initializer_list>
#include <filesystem>

void SimulationCoordinator::SimulationCoordinator_Launch(DataOrganizer *organizer, Inputs *input, CompactRagged* beC, WeatherProcessor* WeatherPro, string inputDirectory)
{
	//Print to Command Prompt Xterm the HydroPlusConfig.xml Model type 
	cout << "Model: " << input->SimulationStringParams["Model_Selection"] << endl;
	//If Model is SpatialTemperatureHydro or StatisticalHydro then Print to Command Prompt Algorithm_SoilKsatDecay and Algorithm_FlowDirection
	if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro" || input->SimulationStringParams["Model_Selection"] == "StatisticalHydro"|| input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {
		//Print to Command Prompt Xterm the HydroPlusConfig.xml Infiltration type 
		cout << "Infiltration: " << input->SimulationStringParams["Algorithm_SoilKsatDecay"] << endl;
		//Print to Command Prompt Xterm the HydroPlusConfig.xml FlowPathAlgorithm type 
		cout << "Flow Path Algorithm: " << input->SimulationStringParams["Algorithm_FlowDirection"] << endl;
	}
	//If Model is SpatialTemperatureHydro then Print to Command Prompt Flag_MultipleStations and Flag_HottestDaySimulated status
	if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") {
		//If SimulationStringParams Flag_MultipleStations == 1 then 
		if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
			//Print to Command Prompt Xterm the HydroPlusConfig.xml Flag_MultipleStations setting = Multiple
			cout << "Weather Stations: Multiple" << endl;
		}
		//Else if single weather station is active
		else {
			//Print to Command Prompt Xterm the HydroPlusConfig.xml Flag_MultipleStations setting = Single
			cout << "Weather Stations: Single" << endl;
		}
		//If SimulationStringParams Flag_HottestDaySimulated == 1 then 
		if (input->SimulationNumericalParams["Flag_HottestDaySimulated"] == 1) {
			//Print to Command Prompt Xterm the HydroPlusConfig.xml Flag_HottestDaySimulated setting = True
			cout << "Hottest Day Simulated: True" << endl;
		}
	}
	//Print to Command Prompt Xterm the computer is running the simulation 
	cout << "Running Simulation ..." << endl;

	//If Model_Selection is StatisticalHydro then call RunStatisticalModel
	if (input->SimulationStringParams["Model_Selection"] == "StatisticalHydro") {
		RunStatisticalModel(organizer, input, beC, WeatherPro, inputDirectory);
	}
	//If Model_Selection is SpatialTemperatureHydro then call RunTemperatureSpatialModel
	else if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") {
		RunTemperatureSpatialModel(organizer, input, beC, WeatherPro, inputDirectory);
	}
	//If Model_Selection is CoolRiver then call RunCoolRiverModel
	else if (input->SimulationStringParams["Model_Selection"] == "CoolRiver")	{
		RunCoolRiverModel(organizer, input, WeatherPro, inputDirectory);
	}
	//If Model_Selection is ECDynamic and run_ECDynamicModel is 1 then call BufferDynamicCalc::YearAverage then simulate
	//Note: When NLCD or discharge data is not provided, run_ECDynamicModel is set to 0 and ...
	//Note: ... the model only extracts and calculates EC and EMC values from the database 
	else if (input->SimulationStringParams["Model_Selection"] == "ECDynamic" && input->run_ECDynamicModel == 1) {
		BufferDynamicCalc::YearAverage(input);
	}
	//If Model_Selection is SpatialBuffer then call BufferModel_HotSpotMaps_DiffusePollution_main
	else if (input->SimulationStringParams["Model_Selection"] == "SpatialBuffer") {
		//Buffer:BufferModel_HotSpotMaps_DiffusePollution_main model called when Model_Selection = SpatialBuffer
		BufferSpatialCalc::BufferModel_HotSpotMaps_DiffusePollution_main(input);
	}
	//If Model_Selection is SpatialBufferwGI then call multiple models, 
	else if (input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {
		//for spatial GI, StatisticalModel will be run twice and using always use append to open the file. 
		//TimeSeries outputs need to be cleaned up before the firdt run
		for (const auto& entry : filesystem::directory_iterator(Inputs::SimulationStringParams["OutputFolder_Path"])) {
			const string filename = entry.path().filename().string();
			if (filename.rfind("TimeSeries", 0) == 0) {  // rfind returns 0 if "TimeSeries" is a prefix
				filesystem::remove_all(entry.path());
			}
		}
		//RunStatisticalModel model called when Model_Selection = SpatialBufferwGI to determine runoff generation
		RunStatisticalModel(organizer, input, beC, WeatherPro, inputDirectory); //get GI inflow Q
		//Buffer:BufferModel_HotSpotMaps_DiffusePollution_main model called when Model_Selection = SpatialBufferwGI
		//Note: Get pollutant C (EC*area*time/Q)
		BufferSpatialCalc::BufferModel_HotSpotMaps_DiffusePollution_main(input);
		//RunStatisticalModel model called when Model_Selection = SpatialBufferwGI to determine runoff trapping in GI
		RunStatisticalModel(organizer, input, beC, WeatherPro, inputDirectory);//run GI
	}				
	//If Model_Selection is CoolBuilding then call EnergyLoad
	else if (input->SimulationStringParams["Model_Selection"] == "CoolBuilding") {
		EnergyLoad_Calc::EnergyLoad(input, WeatherPro);
	}
}

void SimulationCoordinator::RunTemperatureSpatialModel(DataOrganizer *organizer, Inputs *input, CompactRagged* beC, WeatherProcessor* WeatherPro, string inputDirectory)
{
	//SolarCalculator solarCalculator(input) creates single instance of SolarCalculator which is accessed by pointer elsewhere
	SolarCalculator solarCalculator(input);

	//HeatFluxCal input structure from HeatFluxCal.h is referenced as heatFluxCalc object
	//Note: HeatFluxCal object heatFluxCalc created using a constructor that initializes the solarCalculator pointer
	//Note: Optionally, create a pointer without &solarCalculator, supported by the similar constructor
	HeatFluxCal heatFluxCalc(input, &solarCalculator);
	//BlockGroupHandler structure from BlockGroupHandler.h is referenced as handleBlockGroups object
	BlockGroupHandler handleBlockGroups(input);
	
	//TemperatureOutputWriter functions called before entering timeStep loop to prepare files for output at each time step
	//If Flag_PrintBlockGroupDailyAndHourly_str != "0" then process BlockGroup data
	if (input->Flag_PrintBlockGroupDailyAndHourly_str != "0") {
		//Call TemperatureOutputWriter::generate_BlockGroup_Files function to prepare Blockgroup files
		TemperatureOutputWriter::generate_BlockGroup_Files(input);
	}

	//Call TemperatureOutputWriter::generateRowColFiles function to prepare RowColHeat and RowColHydro files
	TemperatureOutputWriter::generateRowColFiles(organizer, input, beC);

	//Iteration_Maximum initialized to IterationCount from HydroPlusConfig.xml file
	int Iteration_Maximum = input->TemperatureCalculationParams["IterationCount"];
	//Flag_AdiabaticLapseRate = 1 to simulate effect of elevation on temperature; = 0 to ignore effect of elevation
	bool Flag_AdiabaticLapseRate = input->TemperatureCalculationParams["Flag_AdiabaticLapseRate"];

	//If Flag_MultipleStations equals 1 then launch algorithm to compute weights for weather stations
	if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
		//Print to Command Prompt Xterm the computer is obtaining multiple station records
		cout << "Processing meteorological data for multiple stations. This may take some time ..." << endl;
		//WeatherProcessor::WeatherDataReadandProcessing function will read in map for all pixels, including NODATA_code
		WeatherProcessor::WeatherDataReadandProcessing(input, organizer);
	}

	//Date_prior_YYYYMMDD initialized to -1 to allow for updates
	int Date_prior_YYYYMMDD = -1;
	//For (int timeStep = 0; timeStep < input->SimulationTimePeriod_timeSteps; ++timeStep) 
	//Note: HydroPlus requires loop through time outside of loop through pixels to implement topographic index water distribution
	//Note: SimulationTimePeriod_timeSteps WeatherProcessor.cpp; dynamically defined to allow for variable time steps
	for (int timeStep = 0; timeStep < input->SimulationTimePeriod_timeSteps; ++timeStep) {

		//If (input->SimulationDate_YYYYMMDD[timeStep] != Date_prior_YYYYMMDD); then current date does not equal prior date
		if (input->SimulationDate_YYYYMMDD[timeStep] != Date_prior_YYYYMMDD) {
			//Date_prior_YYYYMMDD = input->SimulationDate_YYYYMMDD[timeStep]; set prior date to current date
			Date_prior_YYYYMMDD = input->SimulationDate_YYYYMMDD[timeStep];
			//input->CleanTimesPerDayMap(input->IrrigationCountByDay, input->SimulationDate_YYYYMMDD[timeStep]); will clean out array each day
			input->CleanTimesPerDayMap(input->IrrigationCountByDay, input->SimulationDate_YYYYMMDD[timeStep]);
		}

		//If (input->SimulationTimePeriod_timeSteps > 3 && timeStep % int(0.25 * input->SimulationTimePeriod_timeSteps) == 0)
		//Note: Modulus timeStep and 25% of SimulationTimePeriod_timeSteps equals zero provides quarterly updates
		if (input->SimulationTimePeriod_timeSteps > 3 && timeStep % int(0.25 * input->SimulationTimePeriod_timeSteps) == 0) {
			cout << "Simulating time step " << timeStep + 1 << " of " << input->SimulationTimePeriod_timeSteps << endl;
		}
		//Else if (timeStep + 1 == input->SimulationTimePeriod_timeSteps) {
		//Note: Write out the final time step
		else if (timeStep + 1 == input->SimulationTimePeriod_timeSteps) {
			cout << "Simulating time step " << timeStep + 1 << " of " << input->SimulationTimePeriod_timeSteps << endl;
		}

		//StorageDeficit_VadoseZone_TS_Prior_m is constant during time step, used by all folders while loop through DataDrawers iteratively updates AveSMD in StorageDeficit_VadoseZone
		//Note:	StorageDeficit_VadoseZone_TS_Prior_m (m) avoids use of iteratively changing input->RepoDict["StorageDeficit_VadoseZone_m"], updated in StorageDeficit_VadoseZone.cpp for each DataDrawer
		input->InputXml["StorageDeficit_VadoseZone_TS_Prior_m"] = input->RepoDict["StorageDeficit_VadoseZone_m"];

		//Flag_simulateReferenceStation initialized to one at start of each timeStep
		bool Flag_simulateReferenceStation = 1;
		//Flag_Call_Load_TemperatureCalculationParams initialized to one at start of each timeStep
		bool Flag_Call_Load_TemperatureCalculationParams = 1;
		//Flag_Update_VadoseZone_TS_Prior_StandardFolder set to true to limit updates of vadose zone prior timeStep
		bool Flag_Update_VadoseZone_TS_Prior_StandardFolder = true;
		//Flag_Update_VadoseZone_TS_Prior_ReferenceFolder set to true to limit updates of vadose zone prior timeStep
		bool Flag_Update_VadoseZone_TS_Prior_ReferenceFolder = true;
		//Counter_SimulatedUnit is used to identify start of simulated folders, initialized to zero
		input->RepoDict["Counter_SimulatedUnit"] = 0;
		
		//For (int MapPixel_ID = 0; MapPixel_ID < beC->n_drawers(beC); ++MapPixel_ID)
		//Note: beC->n_drawers(beC) is maxRows * maxCols, defined in BuildDataOrganizer; not same as input->DataDrawers.size()
		//Note: Conversion with MapPixel_ID starting at 0, row & col starting at 0 for internal processing
		//Note: Eqs: row=MapPixel_ID / Inputs::nCols; col=MapPixel_ID % Inputs::nCols; MapPixel_ID = (row * nCols) + col
		//Note: Conversion with MapPixel_ID starting at 0, row & col starting at 1 for map input & output
		//Note: Eqs: row=MapPixel_ID / Inputs::nCols+1; col=MapPixel_ID % Inputs::nCols+1; MapPixel_ID = (row * nCols) + col
		for (int MapPixel_ID = 0; MapPixel_ID < beC->n_drawers(beC); ++MapPixel_ID) {

			// Folder_count = beC->folders_in_drawer(beC, MapPixel_ID);
			int Folder_count = beC->folders_in_drawer(beC, MapPixel_ID);

			//If if (Folder_count == 0) then a NODATA_code pixel with empty DataDrawer; continue to next MapPixel_ID
			if (Folder_count == 0) { continue; }

			//For (int DataFolder_ID = 0; DataFolder_ID < Folder_count; ++DataFolder_ID)
			//Note: Folder_count is number of DataFolder in DataDrawer, typically 1
			//Note: Each DataDrawer contains only 1 BulkArea DataFolder in HydroPlusConfig.xml for Model = SpatialTemperatureHydro 
			for (int DataFolder_ID = 0; DataFolder_ID < Folder_count; ++DataFolder_ID) {

				//Flag_Use_ReferenceStationFolder is true if Flag_InputXml_ReferenceStationFolder && isReferenceStationPixel(MapPixel_ID)
				bool Flag_Use_ReferenceStationFolder = (input->Flag_InputXml_ReferenceStationFolder && input->isReferenceStationPixel(MapPixel_ID));

				//If (Flag_Use_ReferenceStationFolder && Flag_Update_VadoseZone_TS_Prior_ReferenceFolder) then update StorageDeficit_VadoseZone_TS_Prior_m
				if (Flag_Use_ReferenceStationFolder && Flag_Update_VadoseZone_TS_Prior_ReferenceFolder) {
					input->InputXml_ReferenceStationFolder["StorageDeficit_VadoseZone_TS_Prior_m"] = input->RepoDict["StorageDeficit_VadoseZone_m"];
					//Flag_Update_VadoseZone_TS_Prior_ReferenceFolder = false until next time step
					Flag_Update_VadoseZone_TS_Prior_ReferenceFolder = false;
				}
				//Else If (not Flag_Use_ReferenceStationFolder && Flag_Update_VadoseZone_TS_Prior_StandardFolder) then update StorageDeficit_VadoseZone_TS_Prior_m
				else if (!Flag_Use_ReferenceStationFolder && Flag_Update_VadoseZone_TS_Prior_StandardFolder) {
					input->InputXml_StandardFolder["StorageDeficit_VadoseZone_TS_Prior_m"] = input->RepoDict["StorageDeficit_VadoseZone_m"];
					//Flag_Update_VadoseZone_TS_Prior_StandardFolder = false until next time step
					Flag_Update_VadoseZone_TS_Prior_StandardFolder = false;
				}
				//InputXml defined with ternary using Flag_InputXml_ReferenceStationFolder, InputXml_ReferenceStationFolder if true, else InputXml_StandardFolder
				//Note: InputXml replaces InputXml_StandardFolder and InputXml_ReferenceStationFolder 
				input->InputXml = Flag_Use_ReferenceStationFolder
					? input->InputXml_ReferenceStationFolder
					: input->InputXml_StandardFolder;
				//InputXml_string defined with ternary, InputXmlString_ReferenceStationFolder if true, else InputXmlString_StandardFolder
				//Note: InputXml_string replaces InputXmlString_StandardFolder and InputXmlString_ReferenceStationFolder 
				input->InputXml_string = Flag_Use_ReferenceStationFolder
					? input->InputXmlString_ReferenceStationFolder
					: input->InputXmlString_StandardFolder;

				//If Flag_DegreeHour_Dependent_AH_Flux_Qcr and Flag_MultipleStations are 1, then enter
				if (input->SimulationNumericalParams["Flag_DegreeHour_Dependent_AH_Flux_Qcr"] == 1 && input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
					//WeatherProcessor::ComputeSpatiallyWeighted_DegreeHour(input, MapPixel_ID) called to compute DegreeHour values
					WeatherProcessor::ComputeSpatiallyWeighted_DegreeHour(input, MapPixel_ID);
				}

				//Counter_SimulatedUnit is used to identify start of simulated folders, increased by one
				input->RepoDict["Counter_SimulatedUnit"] = input->RepoDict["Counter_SimulatedUnit"] + 1;

				//Ref_MapPixel_ID = input->MapPixel_ReferenceStation_IDs[0]; Ref_MapPixel_ID set to default
				int Ref_MapPixel_ID = input->MapPixel_ReferenceStation_IDs[0];

				//If Flag_simulateReferenceStation equals 1, then 1st calculate the reference station, even before MapPixel_ID 0 
				//Note: This conditional will generate Tair_Mesoscale_K and AbsHumidity_Mesoscale_K values needed for local energy balance
				//Note: Flag_simulateReferenceStation is set to 1 each timeStep, then set to 0 at end of this If statement
				if (Flag_simulateReferenceStation == 1) {

					//If Flag_MultipleStations is true then process multiple reference stations
					//Note: This conditional generates Tair_Mesoscale_byStationID_K and AbsHumidity_Mesoscale_byStationID_kg_p_m3 vectors
					if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {

						//For loop with to get Station_Info through all weather stations in simulation
						//Note: StationInfoMap contains all weather stations for simulation; from HydroPlusConfig.xml element RefWeatherLocationAttributeFile 
						for (const auto& Station_Info : WeatherPro->StationInfoMap) {
							//StationID_string is string of Station_Info->first, taken from WeatherPro->StationMapID[], unique identifier of Weather Stations
							const string& StationID_string = Station_Info.first;
							//StationID_int is integer of WeatherPro->StationMapID[], unique identifier of Weather Stations, using C++ stoi command
							int StationID_int = stoi(StationID_string);
							//Station_Index = input->MultipleStationObjectID_to_OrganizerDataDrawersIndex_map[StationID_string]; using correspondence map
							//Note: Station_Index identifies where in organizer->DataDrawers vector the 
							int Station_Index = input->MultipleStationObjectID_to_OrganizerDataDrawersIndex_map[StationID_string];

							//Ref_MapPixel_ID = input->MapPixel_ReferenceStation_IDs[Station_Index]; Ref_MapPixel_ID dynamically updated
							Ref_MapPixel_ID = input->MapPixel_ReferenceStation_IDs[Station_Index];

							//ComputeSpatiallyWeighted_Weather_Inputs function gets weighted meteorological end elevation data for this pixel 
							WeatherProcessor::ComputeSpatiallyWeighted_Weather_Inputs(input, Ref_MapPixel_ID, timeStep);

							//runCalculations function called, sending input pointer, DataFolder ID, and timeStep
							runCalculations(input, beC, WeatherPro, Ref_MapPixel_ID, DataFolder_ID, timeStep);

							//Load_TemperatureCalculationParams function called for ReferenceFolder to initialize parameters
							//Note: This could be called only once per simulation if values were not changing between ref and non-ref folders
							heatFluxCalc.Load_TemperatureCalculationParams(input, Ref_MapPixel_ID, DataFolder_ID, timeStep);
							//EnergyFlux_WaterFlux_DefineTerms function calculates energy balance once for each folder, before water balance 
							heatFluxCalc.EnergyFlux_WaterFlux_DefineTerms(input,beC, Flag_simulateReferenceStation, timeStep, WeatherPro, StationID_string, Ref_MapPixel_ID, DataFolder_ID);
							//Tair_AbsHumidity_LayerBalance function calls AerodynamicResistance_CanopyToMesoLayer and Tair_AbsHumidity_LayerBalance 
							heatFluxCalc.Tair_AbsHumidity_LayerBalance(input, beC, organizer, WeatherPro, Ref_MapPixel_ID, DataFolder_ID, timeStep);
							//Tair_mesoScale_K (K) saved to WeatherPro Tair_Mesoscale_byStationID_K 2D vector with StationID_int
							WeatherPro->Tair_Mesoscale_byStationID_K[StationID_int].push_back(heatFluxCalc.Tair_mesoScale_K);
							//AbsHumidity_mesoScale_kg_p_m3 (K) saved to WeatherPro AbsHumidity_Mesoscale_byStationID_kg_p_m3 2D vector 
							WeatherPro->AbsHumidity_Mesoscale_byStationID_kg_p_m3[StationID_int].push_back(heatFluxCalc.AbsHumidity_mesoScale_kg_p_m3);

							//RescaleVariables_CoolAir function rescales computes corresponding mesoScale values for Tair and AbsoluteHumidity
							heatFluxCalc.RescaleVariables_CoolAir(input, beC, organizer, Ref_MapPixel_ID, DataFolder_ID, timeStep);
							//Tair_mesoScale_final_byStationID_K (K) saved to WeatherPro Tair_mesoScale_final_byStationID_K 2D vector with StationID_int
							WeatherPro->Tair_mesoScale_final_byStationID_K[StationID_int].push_back(heatFluxCalc.Tair_mesoScale_final_K);
							//Tair_mesoScale_K (K) saved to WeatherPro Tair_Mesoscale_byStationID_K 2D vector with StationID_int
							WeatherPro->H_total_byStationID_W_p_m2[StationID_int].push_back(heatFluxCalc.H_total_W_p_m2);

							//ComputeSpatiallyWeighted_Tair_DEM_Variables function will compute spatially weighted meteorological variables
							WeatherProcessor::ComputeSpatiallyWeighted_Tair_DEM_Variables(input, StationID_int, Ref_MapPixel_ID, timeStep);
							//CollectVariables_CoolAir function will store all values to memory
							heatFluxCalc.CollectVariables_CoolAir(input, beC, &solarCalculator, Ref_MapPixel_ID, DataFolder_ID, timeStep);
						}
					}
					//Else If Flag_MultipleStations is false then process a single reference station
					else {
						//runCalculations function called, sending input pointer, DataFolder ID, and timeStep
						runCalculations(input, beC, WeatherPro, Ref_MapPixel_ID, DataFolder_ID, timeStep);
						//Load_TemperatureCalculationParams function called for ReferenceFolder to initialize parameters
						//Note: This could be called only once per simulation if values were not changing between ref and non-ref folders
						heatFluxCalc.Load_TemperatureCalculationParams(input, Ref_MapPixel_ID, DataFolder_ID, timeStep);
						//EnergyFlux_WaterFlux_DefineTerms function calculates energy balance once for each folder, before water balance 
						heatFluxCalc.EnergyFlux_WaterFlux_DefineTerms(input, beC, Flag_simulateReferenceStation, timeStep, WeatherPro, "NA", Ref_MapPixel_ID, DataFolder_ID);
						//Tair_AbsHumidity_LayerBalance function calls AerodynamicResistance_CanopyToMesoLayer and Tair_AbsHumidity_LayerBalance 
						heatFluxCalc.Tair_AbsHumidity_LayerBalance(input, beC, organizer, WeatherPro, Ref_MapPixel_ID, DataFolder_ID, timeStep);

						//RescaleVariables_CoolAir function will rescale latent and sensible energy fluxes, compute corresponding mesoScale values for Tair and AbsoluteHumidity
						heatFluxCalc.RescaleVariables_CoolAir(input, beC, organizer, Ref_MapPixel_ID, DataFolder_ID, timeStep);
						//CollectVariables_CoolAir function will store all values to memory
						heatFluxCalc.CollectVariables_CoolAir(input, beC, &solarCalculator, Ref_MapPixel_ID, DataFolder_ID, timeStep);
					}

					//Flag_simulateReferenceStation set to zero until next timeStep
					Flag_simulateReferenceStation = 0;
					//Flag_Call_Load_TemperatureCalculationParams set to one after definition of reference station values
					Flag_Call_Load_TemperatureCalculationParams = 1;
				}

				//If not isReferenceStationPixel(MapPixel_ID) then enter for non-reference stations
				if (!input->isReferenceStationPixel(MapPixel_ID)) {

					//If Flag_MultipleStations equals one then compute spatially weighted meteorological variables
					if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
						//ComputeSpatiallyWeighted_Weather_Inputs function gets weighted meteorological end elevation data for this pixel 
						WeatherProcessor::ComputeSpatiallyWeighted_Weather_Inputs(input, MapPixel_ID, timeStep);
						//ComputeSpatiallyWeighted_Tair_DEM_Variables function will compute spatially weighted meteorological variables
						//Note: Send -1 in place of StationID_int, which enables function to loop through all stations
						WeatherProcessor::ComputeSpatiallyWeighted_Tair_DEM_Variables(input, -1, MapPixel_ID, timeStep);
					}
					//runCalculations function called, sending input pointer, DataFolder ID, and timeStep
					runCalculations(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);

					//If Flag_Call_Load_TemperatureCalculationParams is one then enter
					if (Flag_Call_Load_TemperatureCalculationParams == 1) {
						//Load_TemperatureCalculationParams function called for Folder to initialize parameters
						//Note: This is only needed once per time step for the non-reference station pixels
						heatFluxCalc.Load_TemperatureCalculationParams(input, MapPixel_ID, DataFolder_ID, timeStep);
						//Flag_Call_Load_TemperatureCalculationParams set to zero until next timeStep or reference station
						Flag_Call_Load_TemperatureCalculationParams = 0;
					}

					//HeatFluxCal.cpp energy balance functions called for all map pixels 
					//EnergyFlux_WaterFlux_DefineTerms function calculates energy balance once for each folder, before water balance 
					heatFluxCalc.EnergyFlux_WaterFlux_DefineTerms(input, beC, Flag_simulateReferenceStation, timeStep, WeatherPro, "NA", MapPixel_ID, DataFolder_ID);
					//Tair_AbsHumidity_LayerBalance_NumericalMethod function calculates temperature and humidity based on iterative numerical method
					heatFluxCalc.Tair_AbsHumidity_LayerBalance_NumericalMethod(input, beC, organizer, WeatherPro, MapPixel_ID,DataFolder_ID, timeStep, Iteration_Maximum);
					//If Flag_AdiabaticLapseRate equals 1 then call AdiabaticLapseRates to simulate lapse rate effect on air temperature
					//Note: Even if Flag_MultipleStations=1 and WeatherProcessor lapse rate adjusted Tair_weighted_K, this adjustment is needed 
					if (Flag_AdiabaticLapseRate == 1) {
						//AdiabaticLapseRates function adjusts air and dew point temperatures based on elevation at environmental lapse rate
						heatFluxCalc.AdiabaticLapseRates(input, MapPixel_ID, DataFolder_ID, timeStep, WeatherPro);
					}
					//RescaleVariables_CoolAir function rescales latent and sensible energy fluxes, computes related mesoScale varaiables
					heatFluxCalc.RescaleVariables_CoolAir(input, beC, organizer, MapPixel_ID,DataFolder_ID, timeStep);
					//CollectVariables_CoolAir function stores all values to memory, screens for an extreme value Tair_K error
					heatFluxCalc.CollectVariables_CoolAir(input, beC, &solarCalculator, MapPixel_ID, DataFolder_ID, timeStep);
				}

				//RunoffSummation::SumAcrossCoverTypes function called to sum folders of non-routed baseflow, pervious flow, impervious flow
				RunoffSummation::SumAcrossCoverTypes(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
				//AggregateOutput::SumCatchmentTotals function called to sum catchment variables into a total for each time step
				AggregateOutput::SumCatchmentTotals(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
				//TemperatureOutputWriter rowColFileBatchProcessor function called to process row and column output 
				TemperatureOutputWriter::rowColFileBatchProcessor(timeStep, MapPixel_ID, DataFolder_ID, input->SimulationStringParams["OutputFolder_Path"], input, beC);

				//If Flag_PrintBlockGroupDailyAndHourly_str != "0" then process BlockGroup data
				if (input->Flag_PrintBlockGroupDailyAndHourly_str != "0") {
					//handleBlockGroups AddToBlockGroup function called to create block group output
					handleBlockGroups.AddToBlockGroup(input, beC, MapPixel_ID, DataFolder_ID);
				}					
			}
		}
		//If Flag_PrintBlockGroupDailyAndHourly_str != "0" then process BlockGroup data
		if (input->Flag_PrintBlockGroupDailyAndHourly_str != "0") {
			//handleBlockGroups TimeStepProcessing function called to process block group output, sending input and timeStep
			handleBlockGroups.TimeStepProcessing(input, timeStep);
		}
		//TemperatureOutputWriter called to write output with MapBasedResultsFileProcessor function
		TemperatureOutputWriter::MapBasedResultsFileProcessor(organizer, input->SimulationStringParams["OutputFolder_Path"], input, beC, WeatherPro, timeStep);
		//AggregateOutput::SumTimeStepTotals function called to sum catchment variables into a total for each time step
		AggregateOutput::SumTimeStepTotals(input, timeStep);
		//RoutingDiffusion::RoutingRunoff will route summed flows from BulkArea for each time step
		RoutingDiffusion::RoutingRunoff(input, timeStep);
		//WriteOutputs::TimeSeriesOutput function writes time series output
		WriteOutputs::TimeSeriesOutput(input, organizer, timeStep);
	}
	//Write weighted inputs
	if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
		WeatherProcessor::MapPixelWeightedWeatherInputsWriter(input);
	}
	//RoutedRunoff_TimeSeries::ScalingRunoffDepth_WritingOutput function scales routed flows to ensure their volume balances
	RoutedRunoff_TimeSeries::ScalingRunoffDepth_WritingOutput(input);

	//WriteOutputs::SummaryOutput function writes summary files 
	WriteOutputs::SummaryOutput(organizer, input);

	cout << "Writing output to: " << input->SimulationStringParams["OutputFolder_Path"] << endl;
	cout << "===" << endl;
}


void SimulationCoordinator::RunStatisticalModel(DataOrganizer *organizer, Inputs *input, CompactRagged* beC, WeatherProcessor* WeatherPro, string inputDirectory)
{
	//If Flag_MultipleStations equals 1 then launch algorithm to compute weights for weather stations
	if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
		//Print to Command Prompt Xterm the computer is obtaining multiple station records
		cout << "Processing meteorological data for multiple stations. This may take some time ..." << endl;
		//WeatherProcessor::WeatherDataReadandProcessing function will read in map for all pixels, including NODATA_code
		WeatherProcessor::WeatherDataReadandProcessing(input, organizer);
	}
	//Loop through all time steps in simulation, where SimulationTimePeriod_timeSteps based on HydroPlusConfig StartDate_YYYYMMDD and StopDate_YYYYMMDD
	for (int timeStep = 0; timeStep < input->SimulationTimePeriod_timeSteps; ++timeStep) {

		//If Flag_MultipleStations equals one then compute spatially weighted meteorological variables
		//Note: Only the center pixel is sent for RunStatisticalModel to obtain weighted station values relative to centroid
		if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
			//ComputeSpatiallyWeighted_Weather_Inputs function computes weighted meteorological end elevation data for this pixel of this time step
			WeatherProcessor::ComputeSpatiallyWeighted_Weather_Inputs(input, WeatherProcessor::centerMapPixel_ID, timeStep);
		}
		//Note: This Xterm output code for regular use, and should not be removed
		//If modulus timeStep and 20% of SimulationTimePeriod_timeSteps equals zero, then write output to Xterm window to update user
		if (timeStep % int(0.25 * input->SimulationTimePeriod_timeSteps) == 0) {
			cout << "Simulating time step " << timeStep + 1 << " of " << input->SimulationTimePeriod_timeSteps << endl;
		}
		//If timeStep equals SimulationTimePeriod_timeSteps, using modulus operator, then write output to Xterm window to update user
		else if (timeStep + 1 == input->SimulationTimePeriod_timeSteps) {
			cout << "Simulating time step " << timeStep + 1 << " of " << input->SimulationTimePeriod_timeSteps << endl;
		}

		//StorageDeficit_VadoseZone_TS_Prior_m is constant during time step, used by all folders while loop through DataDrawers iteratively updates AveSMD in StorageDeficit_VadoseZone
		//Note:	StorageDeficit_VadoseZone_TS_Prior_m (m) avoids use of iteratively changing input->RepoDict["StorageDeficit_VadoseZone_m"], updated in StorageDeficit_VadoseZone.cpp for each DataDrawer
		input->InputXml["StorageDeficit_VadoseZone_TS_Prior_m"] = input->RepoDict["StorageDeficit_VadoseZone_m"];

		//Counter_SimulatedUnit is used to identify start of simulated folders, initialized to zero
		input->RepoDict["Counter_SimulatedUnit"] = 0;

		//Loop through the DataDrawers, where the upper bound is based on the general DataDrawer attribute .size()
		for (int MapPixel_ID = 0; MapPixel_ID < beC->n_drawers(beC); ++MapPixel_ID) {

			// Folder_count = beC->folders_in_drawer(beC, MapPixel_ID);
			int Folder_count = beC->folders_in_drawer(beC, MapPixel_ID);

			//Loop through the DataFolders in each DataDrawer[MapPixel_ID], upper bound indicated by .size()
			for (int DataFolder_ID = 0; DataFolder_ID < Folder_count; ++DataFolder_ID) {

				//Counter_SimulatedUnit is used to identify start of simulated folders, increased by one
				input->RepoDict["Counter_SimulatedUnit"] = input->RepoDict["Counter_SimulatedUnit"] + 1;

				//runCalculations function called, sending input pointer, DataFolder ID, and timeStep
				runCalculations(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);

				//RunoffSummation::SumAcrossCoverTypes function called to sum folders of non-routed baseflow, pervious flow, impervious flow, water flow for BulkArea
				RunoffSummation::SumAcrossCoverTypes(input, beC, MapPixel_ID, DataFolder_ID, timeStep);

				//Loading_WaterQuality_EMC::Pollutant_Load_via_EMC function computes loading via event mean concentration theory
				Loading_WaterQuality_EMC::Pollutant_Load_via_EMC(input, beC, MapPixel_ID, DataFolder_ID, timeStep);

				//If it's a non-bulkarea folder (GI folder), call Pollutant_Removal_via_GI
				if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") != "BulkArea") {
					//Removal_StormwaterPollutants::Pollutant_Removal_via_GI function computes loading via event mean concentration theory
					Removal_StormwaterPollutants::Pollutant_Removal_via_GI(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
				}

				//AggregateOutput::SumCatchmentTotals function called to sum catchment variables into a total for each time step
				AggregateOutput::SumCatchmentTotals(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
				
				//only need to run for the last timestep
				if (timeStep == (input->SimulationTimePeriod_timeSteps - 1)) {
					//if it's running SpatialBuffer with GI model and the flag_reduce_EC_EMC_after_GI tuened on and if HydroPlusConfig.xml Type equals GI
					if (input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI" && beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") != "BulkArea" && BufferSpatialCalc::flag_reduce_EC_EMC_after_GI) {
						//start the process of reducing EC ans EMC value for contributing area
						Removal_StormwaterPollutants::reduce_EC_EMC_after_GI(input);
						//recreate the contributing area weighted EC and EMC map
						BufferSpatialCalc::BufferModel_IntermediateOutputs_writer("total_ec_tn_wtd_reduced");
						BufferSpatialCalc::BufferModel_IntermediateOutputs_writer("total_ec_tp_wtd_reduced");
						BufferSpatialCalc::BufferModel_IntermediateOutputs_writer("total_ec_tss_wtd_reduced");
						BufferSpatialCalc::BufferModel_IntermediateOutputs_writer("total_emc_tn_wtd_reduced");
						BufferSpatialCalc::BufferModel_IntermediateOutputs_writer("total_emc_tp_wtd_reduced");
						BufferSpatialCalc::BufferModel_IntermediateOutputs_writer("total_emc_tss_wtd_reduced");
					}
				}
			}
		}
		
		//AggregateOutput::SumTimeStepTotals function called to sum catchment variables into a total for each time step
		AggregateOutput::SumTimeStepTotals(input, timeStep);
		//WriteOutputs::TimeSeriesOutput function writes time series output
		WriteOutputs::TimeSeriesOutput(input, organizer, timeStep);
		//RoutingDiffusion::RoutingRunoff will route summed flows from BulkArea for each time step
		RoutingDiffusion::RoutingRunoff(input, timeStep);
	}
	//if Flag_MultipleStations == "1", write out the weighted input
	if (input->SimulationNumericalParams["Flag_MultipleStations"] == 1) {
		WeatherProcessor::MapPixelWeightedWeatherInputsWriter(input);
	}
	//RoutedRunoff_TimeSeries::ScalingRunoffDepth_WritingOutput function scales routed flows to ensure their volume balances
	RoutedRunoff_TimeSeries::ScalingRunoffDepth_WritingOutput(input);
	//WriteOutputs::SummaryOutput function writes summary files 
	WriteOutputs::SummaryOutput(organizer, input);

	cout << "Writing output to: " << input->SimulationStringParams["OutputFolder_Path"] << endl;
}

//SimulationCoordinator::runCalculations function defines functors based on model and then sequentially goes through water balance
void SimulationCoordinator::runCalculations(Inputs* input, CompactRagged* beC, WeatherProcessor* WeatherPro, int MapPixel_ID, int DataFolder_ID, int timeStep)
{
	//Note: Functors are defined below with functors.emplace("NameOfFunctor", FunctionCpp::FunctionName)
	//Note: Funtors are called by passing path, e.g., FunctionName(Inputs* input, map<string, function<void(Inputs*, int) >& functors, int timeStep)
	// A uniform callable that call sites use:
	using StepFn = function<void(Inputs*, int)>;
	map<string, StepFn> functors;

	//bindPix wraps the “extended - args” function into a lambda that matches the StepFn signature(void(Inputs*, int))..
	//..Inside that lambda,bcapture beC, MapPixel_ID, and DataFolder_ID from the outer scope and pass them to the real function
	//Register “extended-arg” functions with lambdas that capture the extra params:
	auto bindPix = [beC, MapPixel_ID, DataFolder_ID](void(*fn)(Inputs*, CompactRagged*, int, int, int)) -> StepFn {
		return [=](Inputs* input, int timeStep) {
			fn(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
			};
		};


	if (input->SimulationStringParams["Model_Selection"] == "StatisticalHydro"|| input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {
		functors.emplace("calculateTreeEvaporation", bindPix(&TreeEvaporationCalc::calculateTreeEvaporationStatistical));
		functors.emplace("calculateShortVegEvaporation", bindPix(&ShortVegEvaporationCalc::calculateShortVegEvaporationStatistical));
		//Defining functor for StatisticalHydro for pervious ponded water called calculatePerDepEvap in Evaporation folder 
		functors.emplace("calculatePerDepEvap", bindPix(&OpenWaterEvaporationCalc::calculatePerDepEvapStatistical));
		//Defining functor for StatisticalHydro for impervious ponded water called calculateImpDepEvapTemperature in Evaporation folder 
		functors.emplace("calculateImpDepEvap", bindPix(&OpenWaterEvaporationCalc::calculateImpDepEvapStatistical));
		//Defining functor for StatisticalHydro for water cover ponded water called calculateWaterDepEvapTemperature in Evaporation folder 
		functors.emplace("calculateWaterDepEvap", bindPix(&OpenWaterEvaporationCalc::calculateWaterDepEvapStatistical));
		//Defining functor for StatisticalHydro for snow on ground sublimation called calculateGroundSublimationStatistical in Evaporation folder 
		functors.emplace("calculateGroundSublimation", bindPix(&OpenWaterEvaporationCalc::calculateGroundSublimationStatistical));
	}

	if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") {
		functors.emplace("calculateTreeEvaporation", bindPix(&TreeEvaporationCalc::calculateTreeEvaporationTemperature));
		functors.emplace("calculateShortVegEvaporation", bindPix(&ShortVegEvaporationCalc::calculateShortVegEvaporationTemperature));
		//Defining functor for SpatialTemperatureHydro for pervious ponded water called calculatePerDepEvap in Evaporation folder 
		functors.emplace("calculatePerDepEvap", bindPix(&OpenWaterEvaporationCalc::calculatePerDepEvapTemperature));
		//Defining functor for SpatialTemperatureHydro for impervious ponded water called calculateImpDepEvapTemperature in Evaporation folder 
		functors.emplace("calculateImpDepEvap", bindPix(&OpenWaterEvaporationCalc::calculateImpDepEvapTemperature));
		//Defining functor for SpatialTemperatureHydro for water cover ponded water called calculateWaterDepEvapTemperature in Evaporation folder 
		functors.emplace("calculateWaterDepEvap", bindPix(&OpenWaterEvaporationCalc::calculateWaterDepEvapTemperature));
		//Defining functor for SpatialTemperatureHydro for snow on ground sublimation called calculateGroundSublimationTemperature in Evaporation folder 
		functors.emplace("calculateGroundSublimation", bindPix(&OpenWaterEvaporationCalc::calculateGroundSublimationTemperature));
	}
	
	functors.emplace("calculateInfiltration", (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay")
		? bindPix(&InfiltrationPowerDecay::calculate)
		: bindPix(&InfiltrationExponentialDecay::calculate)
	);

	//BulkArea DataFolder processes
	//If "Type" not equal to BulkArea then allow BulkArea DataFolder processes
	if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "BulkArea") {
		//If SimulationStringParams["Flag_ReadEvaporationData"] not (!=) 1 then call PotentialEvaporationCalc::Method_PenmanMonteith
		if (input->SimulationNumericalParams["Flag_ReadEvaporationData"] != 1) {
			PotentialEvaporationCalc::Method_PenmanMonteith(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);
		}
		TreeInterceptionCalc::TreeInterceptionManager(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		ShortVegInterceptionCalc::ShortVegInterceptionManager(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		functors["calculateTreeEvaporation"](input, timeStep);
		functors["calculateShortVegEvaporation"](input, timeStep);
		functors["calculateGroundSublimation"](input, timeStep);
		SnowAblationOpenAreaCalc::SnowAblationOnGround(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);
		SnowAblationUnderTreeCalc::SnowAblationUnderTree(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);
		SnowAblationUnderShortVegCalc::SnowAblationUnderShortVeg(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);
		InflowsToSurfaceCalc::InflowsToSurfaceImpervious(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		DepressionImperviousCalc::DepressionImpervious(input, beC, MapPixel_ID, DataFolder_ID, timeStep, functors);
		InflowsToSurfaceCalc::InflowsToSurfacePervious(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		functors["calculateInfiltration"](input, timeStep);
		DepressionPerviousCalc::DepressionPervious(input, beC, MapPixel_ID, DataFolder_ID, timeStep, functors);
		InflowsToSurfaceCalc::InflowsToSurfaceWater(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		DepressionWaterCalc::DepressionWater(input, beC, MapPixel_ID, DataFolder_ID, timeStep, functors);
		//Flux_LateralGroundWater::calculate is only called for BulkArea, not for Green Infrastructure (GI)
		//Note: Bulk Area represents uses topographic index (TI) bins with unique soil moisture deficits to represent variable source area theory
		//Note: Green Infrastructure uses average TI and soil moisture deficit to determine average Depth_to_Groundwater_Table_m, and ...
		//Note: ... If Depth_to_Groundwater_Table_m is below the GI layer (e.g., soil, vault) then that layer can drain to vadose zone
		Flux_LateralGroundWater::calculate(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		RunoffSubsurface::calculate(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		StorageDeficit_VadoseZone::calculate(input, beC, MapPixel_ID, DataFolder_ID, timeStep);

		//If GI is operational in HydroPlusConfig.xml then create RepoDict variables to share in GI
		if (input->RepoDict["Flag_GI_Simulated"] == 1) {
			input->RepoDict["Area_BulkArea_m2"] = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2");
			input->RepoDict["Area_BulkArea_Pervious_m2"] = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "PerviousCover_frac");
			input->RepoDict["Area_BulkArea_Impervious_m2"] = beC->by_key(MapPixel_ID, DataFolder_ID, "Area_m2") * beC->by_key(MapPixel_ID, DataFolder_ID, "ImperviousCover_frac");
		}
	}

	//GreenInfrastructure DataFolder processes
	//If "Type" not equal to BulkArea then allow GreenInfrastructure DataFolder processes
	//Note: GI Types are Bioretention, InfilTrench, PermeablePavement, RainGarden, Swale, GreenRoof, RainBarrel, RoofDisconnect
	if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") != "BulkArea") {
		//If SimulationStringParams["Flag_ReadEvaporationData"] not (!=) 1 then call PotentialEvaporationCalc::Method_PenmanMonteith
		if (input->SimulationNumericalParams["Flag_ReadEvaporationData"] != 1) {
			PotentialEvaporationCalc::Method_PenmanMonteith(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);
		}
		TreeInterceptionCalc::TreeInterceptionManager(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		ShortVegInterceptionCalc::ShortVegInterceptionManager(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		functors["calculateTreeEvaporation"](input, timeStep);
		functors["calculateShortVegEvaporation"](input, timeStep);
		functors["calculateGroundSublimation"](input, timeStep);
		SnowAblationOpenAreaCalc::SnowAblationOnGround(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);
		SnowAblationUnderTreeCalc::SnowAblationUnderTree(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);
		SnowAblationUnderShortVegCalc::SnowAblationUnderShortVeg(input, beC, WeatherPro, MapPixel_ID, DataFolder_ID, timeStep);
		InflowsToSurfaceCalc::InflowsToSurfaceImpervious(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		InflowsToSurfaceCalc::InflowsToSurfacePervious(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		InflowsToSurfaceCalc::InflowsToSurfaceWater(input, beC, MapPixel_ID, DataFolder_ID, timeStep);

		//GI_InitializeAndLaunchWaterBalance function initializes variables and launches the GI water balance, calling Flux_Manager_StormwaterDevice
		//Note: BulkArea transfer of water to GI is through Flux_BulkArea_to_GI_Runon_m3 (m3), managed in Inflow_StormwaterDevice.cpp
		//Note: Flux_BulkArea_to_GI_Runon_m3 includes Runoff_Impervious_BAtoGI_m3, Runoff_SatExcess_BAtoGI_m3, Runoff_InfilExcess_BAtoGI_m3.
		//Note: GI transfer of water to groundwater is through GI_Drainage_VadoseZone_m3, managed in StorageDeficit_VadoseZone::calculate function ...
		//Note: ... GI uses average TI and soil moisture deficit to determine average Depth_to_Groundwater_Table_m, and ...
		//Note: ... If Depth_to_Groundwater_Table_m is below the GI layer (e.g., soil, vault) then that layer can drain to vadose zone
		//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 
		Inflow_StormwaterDevice::GI_InitializeAndLaunchWaterBalance(input, beC, MapPixel_ID, DataFolder_ID, timeStep, functors);

		//If Flag_Exfiltration_to_Catchment equals one, then allow for subsurface runoff from saturated zone, and vadose zone drainage to saturated zone
		//Note: When Flag_Exfiltration_to_Catchment != 1, mounding does not occur below GI
		if (beC->by_key(MapPixel_ID, DataFolder_ID, "Flag_Exfiltration_to_Catchment")== 1) {
			RunoffSubsurface::calculate(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
			StorageDeficit_VadoseZone::calculate(input, beC, MapPixel_ID, DataFolder_ID, timeStep);
		}
	}
}

void SimulationCoordinator::RunCoolRiverModel(DataOrganizer* organizer, Inputs* input, WeatherProcessor* WeatherPro, string inputDirectory)
{
	vector<double> solarAzimuth;
	vector<double> solarAltitude;

	//SolarCalculation pointer to instance sc constructed from SolarCalculation class
	SolarCalculation* sc = new SolarCalculation(input);

	//For Loop timeStep through all input->SimulationTimePeriod_timeSteps based on HydroPlusConfig StartDate_YYYYMMDD and StopDate_YYYYMMDD
	for (int timeStep = 0; timeStep < input->SimulationTimePeriod_timeSteps; ++timeStep) {
		//Call SolarCalculation::SolarHour function sending timeStep variable
		sc->SolarHour(input, timeStep);
		//Call SolarCalculation::getSolarAltitude function to populate solarAltitude vector
		solarAltitude.push_back(sc->getSolarAltitude());
		//Call SolarCalculation::getSolarAzimuth function to populate solarAltitude vector
		solarAzimuth.push_back(sc->getSolarAzimuth());
	}
	//StreamTemperature pointer to instance st constructed from StreamTemperature class
	StreamTemperature* st = new StreamTemperature(input, sc);

	//Call StreamTemperature::avgTemperature function
	//Note: avgTemperature function manages most of the CoolRiver functions, including interpolation of spatial and temporal variables
	//Note: Consider refactor avgTemperature to more cleanly divide total distance into river stations that include start (0m) and finish (1980m)
	//Note: Consider refactor avgTemperature to write river station 2D files with actual river station distance; check why it is ordered DS to US
	//Note: Consider refactor avgTemperature to only clarify difference between TotalTimeMod_CN, TotalTimeUns, TotalTime_HR 
	st->avgTemperature();

	//TotalTimeSteps_cnt (hr) is defined by SimulationTimePeriod_timeSteps
	int TotalTimeSteps_cnt = input->SimulationTimePeriod_timeSteps;
	int LengthSimulated_River_m = input->CoolRiverVariables["CR_totDist"];
	//dx_Interval_m (m) is computational interval in streamwise distance x
	int dx_Interval_m = input->CoolRiverVariables["CR_interval_dx"];
	int calcMethod = input->CoolRiverVariables["CR_calcMethod"];

	//TimeSimulated_hr (hr) presumes time steps are in hour; Consider Refactor to allow for other assumptions
	int TimeSimulated_hr = TotalTimeSteps_cnt;
	//CrossSectionsSimulated_cnt represents the LengthSimulated_River_m (m) / dx_Interval_m (m), as computational nodes
	int CrossSectionsSimulated_cnt = (LengthSimulated_River_m / dx_Interval_m);

	// Crank-Nicolson method
	int TotalTimeMod_CN = 1;
	int TotalTimeUns = 1;
	int TotalTime_HR = 1;

	if (calcMethod == 1) {
		TotalTimeMod_CN = TotalTimeSteps_cnt;
	}

	//Call Calculations for calculationMethod1 
	CoolRiverOutputWriter::calculationMethod1(input, inputDirectory, st, calcMethod, LengthSimulated_River_m, TotalTimeMod_CN);

	if (calcMethod == 2) {
		//TotalTimeUns (min) converts TotalTimeSteps_cnt from hours to minutes
		TotalTimeUns = TotalTimeSteps_cnt * 60;
	}

	//Call Calculations for calculationMethod2 
	CoolRiverOutputWriter::calculationMethod2(input, inputDirectory, st, calcMethod, LengthSimulated_River_m, TotalTimeUns, CrossSectionsSimulated_cnt, TotalTimeSteps_cnt);

	if (calcMethod == 3) {
		//TotalTimeUns (min) converts TimeSimulated_hr from hours to minutes
		TotalTime_HR = (TimeSimulated_hr - 1) * 60;
	}

	//Call Calculations for calculationMethod3 
	CoolRiverOutputWriter::calculationMethod3(input, inputDirectory, st, calcMethod, LengthSimulated_River_m, TotalTime_HR, CrossSectionsSimulated_cnt);

	//switch function checks calcMethod value
	switch (calcMethod) {
	//case 1 when HydroPlusConfig.xml calcMethod = 1, steady
	case 1:

		if (input->CoolRiverVariables["CR_Minute_output"] == 0) {
			//When CR_Minute_output is 0, then we are setting values in output for calcMethod1
			CoolRiverOutputWriter::swtichCalMethod1True(input, st, LengthSimulated_River_m, TotalTimeSteps_cnt);
		}
		else {
			//When CR_Minute_output is other than 0, then we are setting values in output for calcMethod1
			CoolRiverOutputWriter::swtichCalMethod1False(input, st, CrossSectionsSimulated_cnt, TotalTimeMod_CN);
		}
		break;

	//case 3 when HydroPlusConfig.xml calcMethod = 3, HEC-RAS steady
	case 3:
		if (input->CoolRiverVariables["CR_Minute_output"] == 0) {
			if (input->CoolRiverVariables["CR_Extended"] == 1) {
				//When CR_Minute_output is 0 and CR_Extended is 1, then we are setting values in output for calcMethod3
				CoolRiverOutputWriter::swtichCalMethod3ExtTrue(input, st, CrossSectionsSimulated_cnt, TimeSimulated_hr);

			}
			//When CR_Minute_output is 0, then we are setting values in output for calcMethod3
			CoolRiverOutputWriter::swtichCalMethod3True(input, st, CrossSectionsSimulated_cnt, TimeSimulated_hr);
		}
		else {
			if (input->CoolRiverVariables["CR_Extended"] == 1) {
				//When CR_Minute_output is other than 0 and CR_Extended is 1, then we are setting values in output for calcMethod3
				CoolRiverOutputWriter::swtichCalMethod3FalseExtTrue(input, st, CrossSectionsSimulated_cnt, TotalTime_HR);

			}
			//When CR_Minute_output is other than 0, then we are setting values in output for calcMethod3
			CoolRiverOutputWriter::swtichCalMethod3False(input, st, CrossSectionsSimulated_cnt, TotalTime_HR);
		}
		break;
	//case 2 when HydroPlusConfig.xml calcMethod = 2, unsteady
	case 2:
		if (input->CoolRiverVariables["CR_Minute_output"] == 0) {
			if (input->CoolRiverVariables["CR_Extended"] == 1) {
				//When CR_Minute_output is 0 and CR_Extended is 1, then we are setting values in output for calcMethod2
				CoolRiverOutputWriter::swtichCalMethod2ExtTrue(input, st, CrossSectionsSimulated_cnt, TotalTimeSteps_cnt);
			}
			//When CR_Minute_output is other than 0 , then we are setting values in output for calcMethod2
			CoolRiverOutputWriter::swtichCalMethod2True(input, st, CrossSectionsSimulated_cnt, TotalTimeSteps_cnt);
		}
		else {
			if (input->CoolRiverVariables["CR_Extended"] == 1) {
				//When CR_Minute_output other than 0 and CR_Extended is 1, then we are setting values in output for calcMethod2
				CoolRiverOutputWriter::swtichCalMethod2FalseExtTrue(input, st, CrossSectionsSimulated_cnt, TotalTimeSteps_cnt);

			}
			//When CR_Minute_output is other than 0, then we are setting values in output for calcMethod2
			CoolRiverOutputWriter::swtichCalMethod2False(input, st, CrossSectionsSimulated_cnt, TotalTimeSteps_cnt);
		}
		break;

	//default is when switch does not find case = 1, 2, or 3
	default:
		cout << "Warning: HydroPlusConfig.xml element CR_calcMethod not set to allowable options 1, 2, or 3." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model uses CR_calcMethod to select appropriate algorithms." << endl;
		cout << "Correction: Set CR_calcMethod to 1, 2, or 3 in the HydroPlusConfig.xml CoolRiver_Variables section." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	delete sc;

	cout << "Writing output to: " << input->SimulationStringParams["OutputFolder_Path"] << endl;
	cout << "===" << endl;
}
