#include "TemperatureOutputWriter.h"
#include "stdlib.h"
#include <fstream>
#include <iomanip>
#include <sstream>
#include <vector>
#include <cstring>
#include <string>
#include <algorithm>

//TemperatureOutputWriter::MapBasedResultsFileProcessor function coordinates time-based map output, calling createMapBasedResultFile function
void TemperatureOutputWriter::MapBasedResultsFileProcessor(DataOrganizer* organizer, string OutputDirectory_String , Inputs* input, CompactRagged* beC, WeatherProcessor* WeatherPro, int timeStep)
{
	int SimulationDate_Output_GDH;

	//Flag_MapBasedOutputAllTimeSteps is defined by value of MapBasedOutputAllTimeSteps, 1 or 0; set in Inputs::ProcessTemperatureExecutionParams
	//Note: Flag_MapBasedOutputAllTimeSteps will control whether printing time-based output maps for all time steps or specific time steps
	bool Flag_MapBasedOutputAllTimeSteps = input->TemperatureCalculationParams["MapBasedOutputAllTimeSteps"];

	//SimulationDate_Output_GDH is string created from SimulationDate_Output_GDH (YYYYMMDDHH) of timeStep, vector defined in readWeatherFile
	SimulationDate_Output_GDH = input->SimulationDate_Output_GDH[timeStep];

	//SimulationDate_GDH_string is created from SimulationDate_Output_GDH (YYYYMMDDHH)
	string SimulationDate_GDH_string = to_string(SimulationDate_Output_GDH);

	//If Flag_MapBasedOutputAllTimeSteps is not 1 or true, then enter
	if (!Flag_MapBasedOutputAllTimeSteps) {
		//For loop with iter counter through size of StartDatesForTimeBasedTempOutput vector, which contains start dates for time-based output
		for (int iter = 0; iter < input->StartDatesForTimeBasedTempOutput.size(); ++iter) {
			//If SimulationDate_Output_GDH witin StartDatesForTimeBasedTempOutput and EndDatesForTimeBasedTempOutput then enter
			if (SimulationDate_Output_GDH >= input->StartDatesForTimeBasedTempOutput[iter] &&
				SimulationDate_Output_GDH <= input->EndDatesForTimeBasedTempOutput[iter]) {

				//Flag_MapBasedOutputAllTimeSteps redefined as true for this time step so it will pass to next conditional
				Flag_MapBasedOutputAllTimeSteps = true;
				//break from for loop before passing to next specific date for time-based output
				break;
			}
		}
	}

	//Call createMapBasedResultFile function to create file and write six rows of header values for ArcGIS ASCII map type
	//If Flag_MapBasedOutputAllTimeSteps is true and PrintAllTimeBasedResults not true, then only print the map type defined in HydroPlusConfig.xml
	//Note: Any HydroPlus variable can be entered as the PrintMapBasedResults variable for printing
	//Note: TemperatureExecutionParams string variables converted to TemperatureCalculationParams numeric variables
	if (Flag_MapBasedOutputAllTimeSteps && input->TemperatureCalculationParams["PrintAllTimeBasedResults"] != 1) {

		//string printMapBasedResults = input->TemperatureExecutionParams["PrintMapBasedResults"] to read string value
		//Note: printMapBasedResults takes value of HydroPlusConfig.xml PrintMapBasedResults, single variable item or CSV variable list
		string printMapBasedResults = input->TemperatureExecutionParams["PrintMapBasedResults"];
		//stringstream stringstream_all(printMapBasedResults) to convert to stringstream value 
		//Note: stringstream_all will be searched in while loop below, finding each stringstream_item
		stringstream stringstream_all(printMapBasedResults);
		string stringstream_item;
		//vector<string> PrintMapBasedResults_variableNames_vec, which will hold all stringstream_item v
		vector<string> PrintMapBasedResults_variableNames_vec;

		//while (getline(stringstream_all, stringstream_item, ',') to split the string based on comma
		//Note: the stringstream_item value must be identical to one of the variable map types defined below in the MapBasedResultsFileProcessor function
		while (getline(stringstream_all, stringstream_item, ',')) {
			//stringstream_item erase to remove or trim any leading or trailing whitespace
			stringstream_item.erase(stringstream_item.begin(), find_if(stringstream_item.begin(), stringstream_item.end(), [](unsigned char ch) { return !isspace(ch); }));
			//stringstream_item erase to remove or trim any leading or trailing whitespace
			stringstream_item.erase(find_if(stringstream_item.rbegin(), stringstream_item.rend(), [](unsigned char ch) { return !isspace(ch); }).base(), stringstream_item.end());
			//PrintMapBasedResults_variableNames_vec increased using push_back with stringstream_item
			PrintMapBasedResults_variableNames_vec.push_back(stringstream_item);
		}

		//For loop to call createMapBasedResultFile for each PrintMapBasedResults_variableName_item within PrintMapBasedResults_variableNames_vec
		for (const string& PrintMapBasedResults_variableName_item : PrintMapBasedResults_variableNames_vec) {
			//createMapBasedResultFile function called with PrintMapBasedResults_variableName_item, e.g., Tair_K, Tair_weighted_K
			createMapBasedResultFile(PrintMapBasedResults_variableName_item, SimulationDate_Output_GDH, input);
		}
	}

	//Call createMapBasedResultFile function to create file and write six rows of header values for ArcGIS ASCII map type
	//If Flag_MapBasedOutputAllTimeSteps equals 1 and PrintAllTimeBasedResults equals 1, then call createMapBasedResultFile function 
	//Note: Any HydroPlus variable can be entered as the PrintMapBasedResults variable for printing
	if (Flag_MapBasedOutputAllTimeSteps && input->TemperatureCalculationParams["PrintAllTimeBasedResults"] == 1) {

		//call createMapBasedResultFile function to write header rows for 1st term in (), which is a variable in HydroPlus.sln
		//Note: Other variables are possible, including Tair_weighted_K for multiple station simulations
		createMapBasedResultFile("Tair_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("Tdew_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("HeatIndex_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("Humidex_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("Twetbulbglobe_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("Twetbulbnatural_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("Twetbulbpsychrometric_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("UTCI_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("WindChill_K", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("AbsHumidity_kg_p_m3", SimulationDate_Output_GDH, input);

		createMapBasedResultFile("ImpNR_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("ImpH_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("ImpLE_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("ImpDeltaQ_W_p_m2", SimulationDate_Output_GDH, input);

		createMapBasedResultFile("TreeNR_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("TreeH_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("TreeLE_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("TreeLEE_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("TreeLET_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("TreeDeltaQ_W_p_m2", SimulationDate_Output_GDH, input);

		createMapBasedResultFile("SVegNR_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("SVegH_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("SVegLE_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("SVegLEE_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("SVegLET_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("SVegDeltaQ_W_p_m2", SimulationDate_Output_GDH, input);

		createMapBasedResultFile("SoilNR_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("SoilH_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("SoilLE_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("SoilDeltaQ_W_p_m2", SimulationDate_Output_GDH, input);

		createMapBasedResultFile("WaterNR_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("WaterH_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("WaterLE_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("WaterDeltaQ_W_p_m2", SimulationDate_Output_GDH, input);

		createMapBasedResultFile("NR_total_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("H_total_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("LE_total_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("DeltaQ_W_p_m2", SimulationDate_Output_GDH, input);
		createMapBasedResultFile("AH_total_W_p_m2", SimulationDate_Output_GDH, input);
	}

	//Call writeToMapBasedResultFile function to write variable values under six rows of header values for ArcGIS ASCII map type
	//If Flag_MapBasedOutputAllTimeSteps is true and PrintAllTimeBasedResults not true, then only print the map type defined in HydroPlusConfig.xml
	//Note: Any HydroPlus variable can be entered as the PrintMapBasedResults variable for printing
	//Note: TemperatureExecutionParams string variables converted to TemperatureCalculationParams numeric variables
	if (Flag_MapBasedOutputAllTimeSteps && input->TemperatureCalculationParams["PrintAllTimeBasedResults"] != 1) {

		string printMapBasedResults = input->TemperatureExecutionParams["PrintMapBasedResults"];
		vector<string> PrintMapBasedResults_variableNames_vec;
		stringstream stringstream_all(printMapBasedResults);
		string stringstream_item;

		// Split the string based on comma
		while (getline(stringstream_all, stringstream_item, ',')) {
			//stringstream_item erase to remove or trim any leading or trailing whitespace
			stringstream_item.erase(stringstream_item.begin(), find_if(stringstream_item.begin(), stringstream_item.end(), [](unsigned char ch) { return !isspace(ch); }));
			//stringstream_item erase to remove or trim any leading or trailing whitespace
			stringstream_item.erase(find_if(stringstream_item.rbegin(), stringstream_item.rend(), [](unsigned char ch) { return !isspace(ch); }).base(), stringstream_item.end());
			//PrintMapBasedResults_variableNames_vec increased using push_back with stringstream_item
			PrintMapBasedResults_variableNames_vec.push_back(stringstream_item);
		}

		//For loop to call writeToMapBasedResultFile for each PrintMapBasedResults_variableName_item within PrintMapBasedResults_variableNames_vec
		for (const string& PrintMapBasedResults_variableName_item : PrintMapBasedResults_variableNames_vec) {
			//writeToMapBasedResultFile function called to write values to map for 1st term in (), which is a variable in HydroPlus.sln
			writeToMapBasedResultFile(PrintMapBasedResults_variableName_item, SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		}
	}
	
	//Call writeToMapBasedResultFile function to write variable values under six rows of header values for ArcGIS ASCII map type
	//If Flag_MapBasedOutputAllTimeSteps equals 1 and PrintAllTimeBasedResults equals 1, then call writeToMapBasedResultFile function for each time-based map type listed below
	//Note: TemperatureExecutionParams string variables converted to TemperatureCalculationParams numeric variables
	if (Flag_MapBasedOutputAllTimeSteps && input->TemperatureCalculationParams["PrintAllTimeBasedResults"]) {

		//call writeToMapBasedResultFile function to write values to map for 1st term in (), which is a variable in HydroPlus.sln
		writeToMapBasedResultFile("Tair_K", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("Tdew_K", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("AbsHumidity_kg_p_m3", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("NR_total_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("ImpNR_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("TreeNR_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SVegNR_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SoilNR_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("WaterNR_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("ImpH_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("ImpLE_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("TreeH_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("TreeLE_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("TreeLEE_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("TreeLET_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SVegH_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SVegLE_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SVegLEE_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SVegLET_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SoilH_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SoilLE_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("WaterH_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("WaterLE_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("H_total_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("LE_total_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("DeltaQ_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("ImpDeltaQ_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("TreeDeltaQ_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SVegDeltaQ_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("SoilDeltaQ_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("WaterDeltaQ_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
		writeToMapBasedResultFile("AH_total_W_p_m2", SimulationDate_Output_GDH, input->SimulationStringParams["OutputFolder_Path"], organizer, input, beC);
	}
}

//Variables defined in TemperatureOutputWriter.h need to be declared in TemperatureOutputWriter.cpp; relates to static vs non-static function
int TemperatureOutputWriter::Folder_Type_ID;

//TemperatureOutputWriter::generateRowColFiles function to create cell-based values for Hydro and Heat variables
void TemperatureOutputWriter::generateRowColFiles(DataOrganizer* organizer, Inputs* input, CompactRagged* beC)
{
	//define folder type for any variables that vary between BulkArea and GI
	int Folder_BulkArea_ID = 0;
	int Folder_GI_ID = 1;
	int row, col;
	
	//Loop through the DataDrawers, where the upper bound is based on the general DataDrawer attribute .size()
	//Note: Conversion between MapPixel_ID and row & column pair, with MapPixel_ID starting at 0, row & col starting at 1 w/ ESRI
	//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);

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

		//Loop through the DataFolders in each Folder_count, upper bound indicated by .size()
		//Note: Each DataDrawer contains only 1 BulkArea DataFolder for in HydroPlusConfig.xml for Model = SpatialTemperatureHydro 
		for (int DataFolder_ID = 0; DataFolder_ID < Folder_count; ++DataFolder_ID) {

			//Note: Conversion uses 0-based indexing: row = MapPixel_ID / nCols; col = MapPixel_ID % nCols;
			//MapPixel_ID = (row * nCols) + col; where MapPixel_ID ranges from 0 to (nRows * nCols - 1)
			row = MapPixel_ID / Inputs::nCols;
			col = MapPixel_ID % Inputs::nCols;

			//locationDataIsToBeOutputted boolean is defined as 1 for true if HydroPlusConfig.xml TemperatureCalculationParams = PrintAllRowColOutput
			bool locationDataIsToBeOutputted = input->TemperatureCalculationParams["PrintAllRowColOutput"];

			//if locationDataIsToBeOutputted is 0 or false, then do not print all row and column pairs
			if (!locationDataIsToBeOutputted) {
				locationDataIsToBeOutputted = input->LocationIdsToPrintTempOutput.find(MapPixel_ID) != input->LocationIdsToPrintTempOutput.end();
			}

			//if locationDataIsToBeOutputted is	1 or true, then print all row and column pairs
			if (locationDataIsToBeOutputted) {

				//If HydroPlusConfig.xml Type equals BulkArea then
				//Note: organizer->DataDrawers[MapPixel_ID][DataFolder_ID] is path to folder
				if (beC->by_key_str(MapPixel_ID, DataFolder_ID, "Type") == "BulkArea") {
					//Folder_Type_ID equals zero, Folder_BulkArea_ID = 0
					Folder_Type_ID = Folder_BulkArea_ID;
				}
				//Else Type is GreenInfrastructure
				else {
					//Folder_Type_ID equals 1, Folder_GI_ID = 1
					Folder_Type_ID = Folder_GI_ID;
				}

				//createHeatRowColFile function creates Row#Col#Heat.csv file for writing during time step loop, aka RowColHeat or Row#Col#Heat
				createHeatRowColFile(row, col, input->SimulationStringParams["OutputFolder_Path"], input);
				//createHeatRowColFile function creates Row#Col#Hydro.csv file for writing during time step loop, aka RowColHydro or Row#Col#Hydro
				createHydroRowColFile(row, col, input->SimulationStringParams["OutputFolder_Path"], input);
			}
		}
	}
}

//TemperatureOutputWriter::rowColFileBatchProcessor function to coordinate writing cell-based values for Hydro and Heat variables
void TemperatureOutputWriter::rowColFileBatchProcessor(int timeStep, int MapPixel_ID, int DataFolder_ID, string OutputDirectory_String , Inputs* input, CompactRagged* beC) {

	//locationDataIsToBeOutputted defined as input->TemperatureCalculationParams["PrintAllRowColOutput"]
	bool locationDataIsToBeOutputted = input->TemperatureCalculationParams["PrintAllRowColOutput"];

	//If not locationDataIsToBeOutputted then enter to set locationDataIsToBeOutputted to value
	if (!locationDataIsToBeOutputted) {
		//locationDataIsToBeOutputted defined as True of False based on comparison of LocationIdsToPrintTempOutput
		locationDataIsToBeOutputted = input->LocationIdsToPrintTempOutput.find(MapPixel_ID) != input->LocationIdsToPrintTempOutput.end();
	}
	//If not locationDataIsToBeOutputted, then the above assignment failed and exit function
	if (!locationDataIsToBeOutputted) {
		return;
	}

	//Note: Conversion uses 0-based indexing: row = MapPixel_ID / nCols; col = MapPixel_ID % nCols;
	//MapPixel_ID = (row * nCols) + col; where MapPixel_ID ranges from 0 to (nRows * nCols - 1)
	int row = MapPixel_ID / Inputs::nCols;
	int col = MapPixel_ID % Inputs::nCols;
	bool printLocation = input->TemperatureCalculationParams["RowColOutputAllTimeSteps"];

	//control to write out 1st time step, connects with SimulationCoordinator.cpp statement
	if (timeStep >= 0) {

		int SimulationDate_Output_GDH;

		//SimulationDate_Output_GDH is string created from SimulationDate_Output_GDH (YYYYMMDDHH) of timeStep, vector defined in readWeatherFile
		SimulationDate_Output_GDH = input->SimulationDate_Output_GDH[timeStep];

		//SimulationDate_GDH_string is created from SimulationDate_Output_GDH (YYYYMMDDHH)
		string SimulationDate_GDH_string = to_string(SimulationDate_Output_GDH);

		//If Not printLocation then enter and see if printLocation should be changed to True based on dates
		if (!printLocation) {
			for (int iter = 0; iter < input->StartDatesForRowColTempOutput.size(); ++iter) {
				if (SimulationDate_Output_GDH >= input->StartDatesForRowColTempOutput[iter] && 
					SimulationDate_Output_GDH <= input->EndDatesForRowColTempOutput[iter]) {
					printLocation = true;
					break;
				}
			}
		}
		//If printLocation then enter and call writeToHeatRowColFile and writeToHydroRowColFile
		if (printLocation) {
			//writeToHeatRowColFile function writes to Row#Col#Heat.csv file during time step loop, aka RowColHeat or Row#Col#Heat
			writeToHeatRowColFile(input, beC, input->SimulationStringParams["OutputFolder_Path"], MapPixel_ID, DataFolder_ID, timeStep, SimulationDate_Output_GDH, row, col);
			//writeToHeatRowColFile function writes to Row#Col#Heat.csv file during time step loop, aka RowColHydro or Row#Col#Hydro
			writeToHydroRowColFile(input, beC, input->SimulationStringParams["OutputFolder_Path"], MapPixel_ID, DataFolder_ID, timeStep, SimulationDate_Output_GDH, row, col);
		}
	}
}

//TemperatureOutputWriter::createHydroRowColFile function writes cell-based Hydro output, aka RowColHydro or Row#Col#Hydro
void TemperatureOutputWriter::createHydroRowColFile(int row, int col, string OutputDirectory_String , Inputs* input) {

	// Construct the filename using string and to_string
	string filename = "Row" + to_string(row) + "Col" + to_string(col) + "Hydro.csv";
	string outputPath = OutputDirectory_String + filename;

	ofstream outfile(outputPath);

	if (!outfile.good()) {
		cout << "Warning: Output folder does not exist or an output file in that folder is open and cannot be overwritten." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Correction: Confirm the HydroPlusConfig.xml parameter OutputFolder_Path and files in that folder are closed." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//Note: Conversion between MapPixel_ID and row & column pair, with MapPixel_ID starting at 0, row & col starting at 1 w/ ESRI
	//Note: Eqs: row=MapPixel_ID / Inputs::nCols+1; col=MapPixel_ID % Inputs::nCols+1; MapPixel_ID = (row * nCols) + col
	int MapPixel_ID = row * Inputs::nCols + col;
	//outfile writes pixel specific data variables
	outfile << "Elevation_m," << input->Elevation_DEM_m[row * Inputs::nCols + col] << endl;
	outfile << "SlopeGround_rad," << input->SlopeGround_rad[row * Inputs::nCols + col] << endl;
	outfile << "AspectGround_N_0_rad," << input->AspectGround_N_0_rad[row * Inputs::nCols + col] << endl;
	outfile << "Topographic_Index," << input->TI_Value[row * Inputs::nCols + col] << endl;
	outfile << "TI_Basin_Average," << input->TI_Value_Average << endl;
	outfile << "NLCD_Class," << input->LandCover_NLCD_Class[row * Inputs::nCols + col] << endl;
	outfile << "Tree_Cover_fraction," << (input->TreeCoverOnPerviousCover_frac[row * Inputs::nCols + col] + input->TreeCoverOnImperviousCover_frac[row * Inputs::nCols + col]) << endl;
	outfile << "Impervious_Cover_fraction," << (input->ImperviousCoverNoTreeCover_frac[row * Inputs::nCols + col] + input->TreeCoverOnImperviousCover_frac[row * Inputs::nCols + col]) << endl;
	outfile << "Impervious_NoTree_fraction," << input->ImperviousCoverNoTreeCover_frac[row * Inputs::nCols + col] << endl;
	

	//outfile receives header names for time series variables
	outfile << "YYYYMMDDHH,";
	outfile << "Interception_RainSnow_TreeCanopy_m,";
	outfile << "Storage_RainSnow_TreeCanopy_m,";
	outfile << "EvaporationSublimation_TreeCanopy_m,";
	outfile << "ThroughFall_RainSnow_TreeCanopy_m,";
	outfile << "Interception_RainSnow_SVegCanopy_m,";
	outfile << "Storage_RainSnow_SVegCanopy_m,";
	outfile << "EvaporationSublimation_SVegCanopy_m,";
	outfile << "ThroughFall_RainSnow_SVegCanopy_m,";
	outfile << "Flux_to_ImperviousArea_Rain_SnowMelt_Irrigation_m,";
	outfile << "Flux_Impervious_to_Pervious_m,";
	outfile << "Flux_to_PerviousArea_Rain_SnowMelt_Irrigation_ImpRunon_m,";
	outfile << "Flux_to_WaterArea_Rain_SnowMelt_m,";
	outfile << "Evaporation_ImperviousPondedWater_m,";
	outfile << "Evaporation_PerviousPondedWater_m,";
	outfile << "Evaporation_WaterPondedWater_m,";
	outfile << "Storage_ImperviousPondedWater_m,";
	outfile << "Storage_PerviousPondedWater_m,";
	outfile << "Storage_WaterPondedWater_m,";
	outfile << "Flux_to_Infiltration_m,";
	outfile << "Drainage_macroPore_m,";
	outfile << "Infiltration_m,";
	//EvapoTranspiration_TS_Sum_m (m) has header EvapoTranspiration_SoilEvapZone_m, only between Soil_FieldCapacity_m3pm3 and Soil_WiltingPoint_m3pm3
	outfile << "EvapoTranspiration_SoilEvapZone_m,";
	//Storage_SoilEvapZone_TS_Sum_m3_p_m3 (m3/m3) has header Storage_SoilEvapZone_m3pm3, only between Soil_FieldCapacity_m3pm3 and Soil_WiltingPoint_m3pm3
	outfile << "Storage_SoilEvapZone_m3pm3,";
	//Storage_SoilEvapZone_TS_Sum_m (m) has header Storage_SoilEvapZone_m, only between Soil_FieldCapacity_m3pm3 and Soil_WiltingPoint_m3pm3
	outfile << "Storage_SoilEvapZone_m,";
	//Drainage_VadoseZone_TS_Sum_m (m3) has header Drainage_VadoseZone_m, only between Soil_SaturationPoint_m3pm3 and Soil_FieldCapacity_m3pm3
	outfile << "Drainage_VadoseZone_m,";
	//StorageDeficit_VadoseZone_TS_Sum_m (m) has header StorageDeficit_VadoseZone_m, is storage capacity for gravitational water across entire catchment vadose zone
	outfile << "StorageDeficit_VadoseZone_m,";
	outfile << "Groundwater_surficial_frac,";
	outfile << "Runoff_Impervious_m,";
	outfile << "Runoff_Pervious_m,";
	outfile << "Runoff_Water_m,";
	outfile << "Runoff_Surface_m,";
	outfile << "Runoff_Subsurface_m,";
	outfile << "Runoff_m";
	outfile << endl;
	outfile.close();
}

//TemperatureOutputWriter::writeToHydroRowColFile function writes cell-based Hydro time series for folder row col pair, aka RowColHydro or Row#Col#Hydro
void TemperatureOutputWriter::writeToHydroRowColFile(Inputs* input, CompactRagged* beC, string OutputDirectory_String , int DataDrawer_ID, int DataFolder_ID, int timeStep, int SimulationDate_Output_GDH, int row, int col) {

	// Construct the filename using string and to_string
	string filename = "Row" + to_string(row) + "Col" + to_string(col) + "Hydro.csv";
	string outputPath = OutputDirectory_String + filename;

	// Open the output file in append mode
	ofstream outfile(outputPath, ios::app);

	if (!outfile.good()) {
		cout << "Cannot write Row#Col#Hydro.csv file" << endl;
		return;
	}
	
	//Folder_Type_ID equal to zero, which is BulkArea folder
	int Folder_Type_ID = 0;

	outfile << SimulationDate_Output_GDH << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Interception_RainSnow_TreeCanopy_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_RainSnow_TreeCanopy_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "EvaporationSublimation_TreeCanopy_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "ThroughFall_RainSnow_TreeCanopy_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Interception_RainSnow_SVegCanopy_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_RainSnow_SVegCanopy_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "EvaporationSublimation_SVegCanopy_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "ThroughFall_RainSnow_SVegCanopy_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Flux_to_ImperviousArea_Rain_SnowMelt_Irrigation_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Flux_Impervious_to_Pervious_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Flux_to_PerviousArea_Rain_SnowMelt_Irrigation_ImpRunon_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Flux_to_WaterArea_Rain_SnowMelt_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_ImperviousPondedWater_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_PerviousPondedWater_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Evaporation_WaterPondedWater_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_ImperviousPondedWater_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_PerviousPondedWater_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_WaterPondedWater_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Flux_to_Infiltration_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Drainage_macroPore_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Infiltration_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "EvapoTranspiration_SoilEvapZone_m") << ",";
	//Storage_SoilEvapZone_TS_Sum_m3_p_m3 (m3/m3) is only between Theta_FieldCapacity and Theta_WiltingPoint
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_SoilEvapZone_m3_p_m3") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Storage_SoilEvapZone_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Drainage_VadoseZone_m") << ",";
	outfile << input->RepoDict["StorageDeficit_VadoseZone_m"] << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Groundwater_surficial_frac") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_Impervious_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_Pervious_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_Water_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_Surface_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_Subsurface_m") << ",";
	outfile << beC->by_key(DataDrawer_ID, DataFolder_ID, "Runoff_m");
	outfile << endl;
	//outfile close function
	outfile.close();
}

//createHeatRowColFile function creates Row#Col#Heat.csv file for writing during time step loop, aka RowColHeat or Row#Col#Heat
void TemperatureOutputWriter::createHeatRowColFile(int row, int col, string OutputDirectory_String , Inputs* input) {
	
	string filename = "Row" + to_string(row) + "Col" + to_string(col) + "Heat.csv";
	string outputPath = OutputDirectory_String + filename;

	ofstream outfile(outputPath);

	if (!outfile.good()) {
		cout << "Warning: Output folder does not exist or an output file in that folder is open and cannot be overwritten." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Correction: Confirm the HydroPlusConfig.xml parameter OutputFolder_Path and files in that folder are closed." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//Note: Conversion between MapPixel_ID and row & column pair, with MapPixel_ID starting at 0, row & col starting at 1 w/ ESRI
	//Note: Eqs: row=MapPixel_ID / Inputs::nCols+1; col=MapPixel_ID % Inputs::nCols+1; MapPixel_ID = (row * nCols) + col
	int MapPixel_ID = row * Inputs::nCols + col;
	//outfile writes pixel specific data variables
	outfile << "Elevation_m," << input->Elevation_DEM_m[row * Inputs::nCols + col] << endl;
	outfile << "SlopeGround_rad," << input->SlopeGround_rad[row * Inputs::nCols + col] << endl;
	outfile << "AspectGround_N_0_rad," << input->AspectGround_N_0_rad[row * Inputs::nCols + col] << endl;
	outfile << "Topographic_Index," << input->TI_Value[row * Inputs::nCols + col] << endl;
	outfile << "TI_Basin_Average," << input->TI_Value_Average << endl;
	outfile << "NLCD_Class," << input->LandCover_NLCD_Class[row * Inputs::nCols + col] << endl;
	outfile << "Tree_Cover_fraction," << (input->TreeCoverOnPerviousCover_frac[row * Inputs::nCols + col] + input->TreeCoverOnImperviousCover_frac[row * Inputs::nCols + col]) << endl;
	outfile << "Impervious_Cover_fraction," << (input->ImperviousCoverNoTreeCover_frac[row * Inputs::nCols + col] + input->TreeCoverOnImperviousCover_frac[row * Inputs::nCols + col]) << endl;
	outfile << "Impervious_NoTree_fraction," << input->ImperviousCoverNoTreeCover_frac[row * Inputs::nCols + col] << endl;

	//outfile receives header names for time series variables
	outfile << "YYYYMMDDHH,";
	outfile << "Tair_K,";
	outfile << "Tdew_K,";
	outfile << "HeatIndex_K,";
	outfile << "Humidex_K,";
	outfile << "Twetbulbglobe_K,";
	outfile << "Twetbulbnatural_K,";
	outfile << "Twetbulbpsychrometric_K,";
	outfile << "UTCI_K,";
	outfile << "Tglobe_K,";
	outfile << "WindChill_K,";
	outfile << "AbsHumidity_kg_p_m3,";
	outfile << "ImpNR_W_p_m2,";
	outfile << "ImpH_W_p_m2,";
	outfile << "ImpLE_W_p_m2,";
	outfile << "ImpDeltaQ_W_p_m2,";
	outfile << "Imp_AH_W_p_m2,";
	outfile << "TreeNR_W_p_m2,";
	outfile << "TreeH_W_p_m2,";
	outfile << "TreeLE_W_p_m2,";
	outfile << "TreeLEE_W_p_m2,";
	outfile << "TreeLET_W_p_m2,";
	outfile << "TreeDeltaQ_W_p_m2,";
	outfile << "Tree_AH_W_p_m2,";
	outfile << "SVegNR_W_p_m2,";
	outfile << "SVegH_W_p_m2,";
	outfile << "SVegLE_W_p_m2,";
	outfile << "SVegLEE_W_p_m2,";
	outfile << "SVegLET_W_p_m2,";
	outfile << "SVegDeltaQ_W_p_m2,";
	outfile << "SVeg_AH_W_p_m2,";
	outfile << "SoilNR_W_p_m2,";
	outfile << "SoilH_W_p_m2,";
	outfile << "SoilLE_W_p_m2,";
	outfile << "SoilDeltaQ_W_p_m2,";
	outfile << "Soil_AH_W_p_m2,";
	outfile << "WaterNR_W_p_m2,";
	outfile << "WaterH_W_p_m2,";
	outfile << "WaterLE_W_p_m2,";
	outfile << "WaterDeltaQ_W_p_m2,";
	outfile << "Water_AH_W_p_m2,";
	outfile << "NR_W_p_m2,";
	outfile << "H_total_W_p_m2,";
	outfile << "LE_total_W_p_m2,";
	outfile << "DeltaQ_W_p_m2,";
	outfile << "AnthropogenicHeat_total_W_p_m2,";
	outfile << "WindSpeed_canopy_to_mesolayer_m_p_s,";
	outfile << "WindSpeed_folder_m_p_s,";
	outfile << "Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m,";
	outfile << "Resistance_Aerodynamic_Impervious_s_p_m,";
	outfile << "Resistance_Aerodynamic_Tree_s_p_m,";
	outfile << "Resistance_Aerodynamic_SVeg_s_p_m,";
	outfile << "Resistance_Aerodynamic_Soil_s_p_m,";
	outfile << "Resistance_Aerodynamic_Water_s_p_m,";
	outfile << "Resistance_Surface_Tree_s_p_m,";
	outfile << "Resistance_Surface_SVeg_s_p_m,";
	outfile << "Resistance_Surface_Soil_s_p_m,";
	outfile << "Tair_mesoScale_K,";
	outfile << "AbsHumidity_mesoScale_kg_p_m3,";

	outfile << "SurfaceSlopeAngle_rad,";
	outfile << "SurfaceAzimuthAngle_rad,";
	outfile << "LatitudeAngle_rad,";
	outfile << "DeclinationAngle_rad,";
	outfile << "HourAngle_Solar_rad,";
	outfile << "ZenithAngle_rad,";
	outfile << "Cos_SolarAzimuthAngle_Reda,";
	outfile << "Cos_IncidenceAngle_Reda";
	outfile << endl;
	outfile.close();
}

//TemperatureOutputWriter::writeToHeatRowColFile function writes cell-based Heat time series for folder row col pair, aka RowColHeat or Row#Col#Heat
void TemperatureOutputWriter::writeToHeatRowColFile(Inputs* input, CompactRagged* beC, string OutputDirectory_String , int MapPixel_ID, int DataFolder_ID, int timeStep, int SimulationDate_Output_GDH, int row, int col)
{
	
	string filename = "Row" + to_string(row) + "Col" + to_string(col) + "Heat.csv";
	ofstream outfile(OutputDirectory_String + filename, ios::app);

	if (!outfile.good()) {
		cout << "Cannot write Row#Col#Heat.csv file" << endl;
		return;
	}

	//Note: Conversion between MapPixel_ID and row & column pair, with MapPixel_ID starting at 0, row & col starting at 1 w/ ESRI
	//Note: Eqs: row=MapPixel_ID / Inputs::nCols+1; col=MapPixel_ID % Inputs::nCols+1; MapPixel_ID = (row * nCols) + col
	//int MapPixel_ID = row * Inputs::nCols + col;

	//outfile receives variables to cell-based output
	outfile << SimulationDate_Output_GDH << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Tdew_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "HeatIndex_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Humidex_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Twetbulbglobe_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Twetbulbnatural_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Twetbulbpsychrometric_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "UTCI_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Tglobe_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "WindChill_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "AbsHumidity_kg_p_m3") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "ImpNR_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "ImpH_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "ImpLE_W_p_m2") << ",";
	outfile << HeatFluxCal::ImpDeltaQ_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << HeatFluxCal::Imp_AH_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "TreeNR_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "TreeH_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLE_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLEE_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLET_W_p_m2") << ",";
	outfile << HeatFluxCal::TreeDeltaQ_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << HeatFluxCal::Tree_AH_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "SVegNR_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "SVegH_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLE_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLEE_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLET_W_p_m2") << ",";
	outfile << HeatFluxCal::SVegDeltaQ_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << HeatFluxCal::SVeg_AH_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "SoilNR_W_p_m2") << ",";
	outfile <<beC->by_key(MapPixel_ID, DataFolder_ID, "SoilH_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "SoilLE_W_p_m2") << ",";
	outfile << HeatFluxCal::SoilDeltaQ_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << HeatFluxCal::Soil_AH_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "WaterNR_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "WaterH_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "WaterLE_W_p_m2") << ",";
	outfile << HeatFluxCal::WaterDeltaQ_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << HeatFluxCal::Water_AH_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << HeatFluxCal::NR_total_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "H_total_W_p_m2") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "LE_total_W_p_m2") << ",";
	outfile << HeatFluxCal::DeltaQ_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << HeatFluxCal::AH_total_W_p_m2_vec[MapPixel_ID] << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "WindSpeed_canopy_to_mesolayer_m_p_s") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "WindSpeed_folder_m_p_s") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Impervious_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Tree_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_SVeg_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Soil_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Aerodynamic_Water_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Surface_Tree_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Surface_SVeg_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Resistance_Surface_Soil_s_p_m") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_mesoScale_final_K") << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "AbsHumidity_mesoScale_final_kg_p_m3") << ",";
	outfile << input->SlopeGround_rad[MapPixel_ID] << ",";
	outfile << input->AspectGround_N_0_rad[MapPixel_ID] << ",";
	outfile << input->Latitude_rad << ",";
	outfile << input->DeclinationAngle_Solar_rad[timeStep] << ",";
	outfile << input->HourAngle_Solar_rad[timeStep] << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "ZenithAngle_Solar_rad") << ",";
	outfile << cos(beC->by_key(MapPixel_ID, DataFolder_ID, "AzimuthAngle_Solar_N_0_rad")) << ",";
	outfile << beC->by_key(MapPixel_ID, DataFolder_ID, "Cos_IncidenceAngle_Solar");
	outfile << endl;
	outfile.close();
}

//TemperatureOutputWriter::createMapBasedResultFile function writes time-based map output header rows suitable for ArcGIS ASCII map .asc file
void TemperatureOutputWriter::createMapBasedResultFile(string outputVarName, int ts, Inputs* input)
{
	//string fName created by combining OutputDirectory_String , outputVarName, _, and to_string function of ts, ending with .asc
	string fName = input->SimulationStringParams["OutputFolder_Path"] + outputVarName + "_" + to_string(ts) + ".asc";
	//ofstream opens outfile fName to write header information
	ofstream outfile(fName);

	if (!outfile.good()) {
		cout << "Warning: Output folder does not exist or an output file in that folder is open and cannot be overwritten." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Correction: Confirm the HydroPlusConfig.xml parameter OutputFolder_Path and files in that folder are closed." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//nCols written to extend to length of nCols with width of 20, using width of 20 
	outfile << "ncols" << setw(20) << Inputs::nCols << endl;
	//nRows written to extend to length of nCols with width of 20, using width of 20 
	outfile << "nrows" << setw(20) << Inputs::nRows << endl;
	//xllcorner written to extend to length of nCols with width of 20, using width of 16
	outfile << "xllcorner" << setw(16) << Inputs::xllcorner_m << endl;
	//yllcorner written to extend to length of nCols with width of 20, using width of 16 
	outfile << "yllcorner" << setw(16) << Inputs::yllcorner_m << endl;
	//cellsize written to extend to length of nCols with width of 20, using width of 17 
	outfile << "cellsize" << setw(17) << Inputs::Length_Pixel_Side_m << endl;
	//NODATA_value written to extend to length of nCols with width of 20, using width of 13 
	outfile << "NODATA_value" << setw(13) << Inputs::NODATA_code << endl;
	outfile.close();
}

//writeToMapBasedResultFile function writes time-based map output values below header suitable for ArcGIS ASCII map .asc file
void TemperatureOutputWriter::writeToMapBasedResultFile(string outputVarName, int date, string& OutputDirectory_String , DataOrganizer* organizer, Inputs* input, CompactRagged* beC)
{
	//string fName created by combining OutputDirectory_String , outputVarName, _, and to_string function of ts, ending with .asc
	//Note: The ASCII .asc file created with a header in createMapBasedResultFile function is re-opened here for writing values
	string fName = OutputDirectory_String  + outputVarName + "_" + to_string(date) + ".asc";
	//ofstream opens outfile fName with ios::app function to append
	ofstream outfile(fName, ios::app);
	//set the DataFolder_ID = 0; old version: folder defined as organizer->DataDrawers[MapPixel_ID][0]
	int DataFolder_ID = 0;

	if (!outfile.good()) {
		cout << "Cannot write file: " + fName << endl;
		return;
	}

	//mapPixel_Count = Inputs::nRows * Inputs::nCols
	int mapPixel_Count = Inputs::nRows * Inputs::nCols;

	//for (int MapPixel_ID = 0; MapPixel_ID < mapPixel_Count; ++MapPixel_ID)
	//Note: mapPixel_Count ensures map rows and columns are printed, while beC->n_drawers(beC) can exceed map size
	for (int MapPixel_ID = 0; MapPixel_ID < mapPixel_Count; ++MapPixel_ID) {

		//Note: Conversion uses 0-based indexing: row = MapPixel_ID / nCols; col = MapPixel_ID % nCols;
		//MapPixel_ID = (row * nCols) + col; where MapPixel_ID ranges from 0 to (nRows * nCols - 1)
		int col = MapPixel_ID % Inputs::nCols;
		//newLine is boolean variable set to false or true depending on outcome of ternary operator, testing condition col not equal to nCols
		bool newLine = ((col + 1) == Inputs::nCols);

		//If LandCover_NLCD_Class[MapPixel_ID] equals NODATA_code dem.asc then
		if (int(input->LandCover_NLCD_Class[MapPixel_ID]) == Inputs::NODATA_code) {
			//write NODATA_code to map
			outfile << Inputs::NODATA_code;
		}
		//Else If data value, then write variable to the map output
		else {
			//HeatFluxCal variables can be written without storage to folder
			if (outputVarName == "NR_total_W_p_m2") {
				outfile << HeatFluxCal::NR_total_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "DeltaQ_W_p_m2") {
				outfile << HeatFluxCal::DeltaQ_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "WaterDeltaQ_W_p_m2") {
				outfile << HeatFluxCal::WaterDeltaQ_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "SoilDeltaQ_W_p_m2") {
				outfile << HeatFluxCal::SoilDeltaQ_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "AH_total_W_p_m2") {
				outfile << HeatFluxCal::AH_total_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "Imp_AH_W_p_m2") {
				outfile << HeatFluxCal::Imp_AH_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "Tree_AH_W_p_m2") {
				outfile << HeatFluxCal::Tree_AH_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "SVeg_AH_W_p_m2") {
				outfile << HeatFluxCal::SVeg_AH_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "Soil_AH_W_p_m2") {
				outfile << HeatFluxCal::Soil_AH_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "Water_AH_W_p_m2") {
				outfile << HeatFluxCal::Water_AH_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "SVegDeltaQ_W_p_m2") {
				outfile << HeatFluxCal::SVegDeltaQ_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "TreeDeltaQ_W_p_m2") {
				outfile << HeatFluxCal::TreeDeltaQ_W_p_m2_vec[MapPixel_ID];
			}
			else if (outputVarName == "ImpDeltaQ_W_p_m2") {
				outfile << HeatFluxCal::ImpDeltaQ_W_p_m2_vec[MapPixel_ID];
			}
			//Else If another data value, then write folder value for outputVarName
			else {
				//outfile has appended beC->by_key outputVarName variable value 
				//Note: Any HydroPlus variable can be written to map
				outfile << beC->by_key(MapPixel_ID, DataFolder_ID, outputVarName);
			}
		}
		//If newLine boolean is true then col == nCol and line return needed
		if (newLine) {
			//outfile receives endl function for end of line return
			outfile << endl;
		}
		//Else If then newLine boolen is false
		else {
			//outfile receives space between variables
			outfile << " ";
		}
	}
	//outfile uses close function
	outfile.close();
}

//TemperatureOutputWriter::generate_BlockGroup_Files function coordinates block group averages
void TemperatureOutputWriter::generate_BlockGroup_Files(Inputs* input)
{

	//Flag_needDailyBlockGroupFiles is set to true if Flag_PrintBlockGroupDailyAndHourly_str.find is 1 or 2, not npos
	//Note: PrintBlockGroupDailyAndHourly options: 0=no output; 1=day and hour ouput; 2=day only output; 3=hour only output
	bool Flag_needDailyBlockGroupFiles =
		(input->Flag_PrintBlockGroupDailyAndHourly_str.find('1') != string::npos ||
			input->Flag_PrintBlockGroupDailyAndHourly_str.find('2') != string::npos);
	//Flag_needHourlyBlockGroupFiles is set to true if Flag_PrintBlockGroupDailyAndHourly_str.find is 1 or 3, not npos
	//Note: PrintBlockGroupDailyAndHourly options: 0=no output; 1=day and hour ouput; 2=day only output; 3=hour only output
	bool Flag_needHourlyBlockGroupFiles =
		(input->Flag_PrintBlockGroupDailyAndHourly_str.find('1') != string::npos ||
			input->Flag_PrintBlockGroupDailyAndHourly_str.find('3') != string::npos);


	//For Loop from blockGroupIDToItsLocation.begin() through blockGroupIDToItsLocation.end() using iter as counter
	for (int blockGroupNum : input->Unique_BlockGroup_IDs) {
		// Skip BGs that are not in any requested range when PrintAllBlockGroups == 0
		if (input->TemperatureCalculationParams["PrintAllBlockGroups"] == 0) {
			bool inRange = false;
			for (size_t k = 0; k < input->StartingBlockIdsToPrintTempOutput.size(); ++k) {
				if (blockGroupNum >= input->StartingBlockIdsToPrintTempOutput[k] && blockGroupNum <= input->EndingBlockIdsToPrintTempOutput[k]) {
					inRange = true;
					break;
				}
			}
			if (!inRange) {
				continue; // skip this BG entirely
			}
		}

		//If Flag_needHourlyBlockGroupFiles is true then
		if (Flag_needHourlyBlockGroupFiles) {
			//call createBlockGroupHourlyFile function
			createBlockGroupHourlyFile(input->SimulationStringParams["OutputFolder_Path"], blockGroupNum);
		}
		//If Flag_needDailyBlockGroupFiles is true then
		if (Flag_needDailyBlockGroupFiles) {
			//call createBlockGroupDailyFile function
			createBlockGroupDailyFile(input->SimulationStringParams["OutputFolder_Path"], blockGroupNum);
		}
	}
}

//TemperatureOutputWriter::createBlockGroupHourlyFile function writes hourly block group header
void TemperatureOutputWriter::createBlockGroupHourlyFile(string OutputDirectory_String , int blockGroupIndex)
{
	//FileName_str contains BlockGroup defined as string, appended to blockGroupIndex (converted to string) and .csv to create file name
	string FileName_str = "BlockGroup" + to_string(blockGroupIndex) + ".csv";
	ofstream outfile(OutputDirectory_String  + FileName_str);

	if (!outfile.good()) {
		cout << "Warning: Output folder does not exist or an output file in that folder is open and cannot be overwritten." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Correction: Confirm the HydroPlusConfig.xml parameter OutputFolder_Path and files in that folder are closed." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}
	//Header written as YYYYMMDDHH,Tair_K,Tdew_K,HeatIndex_K,Humidex_K,Twetbulbglobe_K,Twetbulbnatural_K,Twetbulbpsychrometric_K,UTCI_K,WindChill_K
	outfile << "YYYYMMDDHH,";
	outfile << "Tair_K,";
	outfile << "Tdew_K,";
	outfile << "HeatIndex_K,";
	outfile << "Humidex_K,";
	outfile << "Twetbulbglobe_K,";
	outfile << "Twetbulbnatural_K,";
	outfile << "Twetbulbpsychrometric_K,";
	outfile << "UTCI_K,";
	outfile << "WindChill_K";
	outfile << endl;
	outfile.close();
}

//TemperatureOutputWriter::writeToBlockGroupHourlyFile function writes hourly block group variable values
void TemperatureOutputWriter::writeToBlockGroupHourlyFile(Inputs* input, string OutputDirectory_String , string date, int blockGroupIndex)
{
	int blockGroup_ID = blockGroupIndex;
	//FileName_str contains BlockGroup defined as string, appended to blockGroup_ID (converted to string) and .csv to create file name
	string FileName_str = "BlockGroup" + to_string(blockGroup_ID) + ".csv";

	ofstream outfile(OutputDirectory_String  + FileName_str, ios::app);

	if (!outfile.good()) {
		cout << "Warning: Output folder does not exist or an output file in that folder is open and cannot be overwritten." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Correction: Confirm the HydroPlusConfig.xml parameter OutputFolder_Path and files in that folder are closed." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	outfile << date << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "Tair_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID") << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "Tdew_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID") << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "HeatIndex_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID") << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "Humidex_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID") << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "Twetbulbglobe_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID") << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "Twetbulbnatural_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID") << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "Twetbulbpsychrometric_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID") << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "UTCI_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID") << ",";
	outfile << input->by_BlockGroup(blockGroup_ID, "WindChill_K") / input->by_BlockGroup(blockGroup_ID, "pixelsCount_of_bgID");
	outfile << endl;
}

//TemperatureOutputWriter::createBlockGroupDailyFile function writes daily block group header
void TemperatureOutputWriter::createBlockGroupDailyFile(string OutputDirectory_String , int blockGroupIndex)
{
	//FileName_str contains BlockGroup defined as string, appended to blockGroupIndex (converted to string) and .csv to create file name
	string FileName_str = "BlockGroupDaily" + to_string(blockGroupIndex) + ".csv";
	ofstream outfile(OutputDirectory_String  + FileName_str);

	if (!outfile.good()) {
		cout << "Warning: Output folder does not exist or an output file in that folder is open and cannot be overwritten." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Correction: Confirm the HydroPlusConfig.xml parameter OutputFolder_Path and files in that folder are closed." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//HeatIndex terms are from Rothfusz, Lans P. National Weather Service Technical Attachment SR 90-23 (1990) 
	//HeatIndex_Lin terms are from Lin, S., M. Luo, R. J. Walker, X. Liu, S. A. Hwang and R. Chinery (2009)

	outfile << "YYYYMMDD,"; 
	outfile << "Tair_Max_K,"; 
	outfile << "Tair_Min_K,"; 
	outfile << "Tair_Avg_K,";
	outfile << "Tdew_Avg_K,";
	outfile << "HeatIndex_Max_K,";
	outfile << "HeatIndex_Avg_K,"; 
	outfile << "Humidex_Max_K,";
	outfile << "Humidex_Avg_K,";
	outfile << "Twetbulbglobe_Max_K,";
	outfile << "Twetbulbglobe_Avg_K,";
	outfile << "Twetbulbnatural_Max_K,";
	outfile << "Twetbulbnatural_Avg_K,";
	outfile << "Twetbulbpsychrometric_Max_K,";
	outfile << "Twetbulbpsychrometric_Avg_K,";
	outfile << "UTCI_Max_K,";
	outfile << "UTCI_Avg_K,";
	//Commented out but retained in case needed by future study, 2025
	//outfile << "HeatIndex_Lin_Max_K,";
	//outfile << "HeatIndex_Lin_Avg_K,";
	//outfile << "HeatIndex_Lin_AvgDaily_K,"; 
	outfile << "WindChill_Min_K,"; 
	outfile << "WindChill_Avg_K";
	outfile << endl;
	outfile.close();
}

//TemperatureOutputWriter::writeToBlockGroupDailyFile function writes daily block group variable value
void TemperatureOutputWriter::writeToBlockGroupDailyFile(Inputs* input, string OutputDirectory_String , string date, int blockGroupIndex, int timeStep)
{
	int blockGroup_ID = blockGroupIndex;
	//FileName_str contains BlockGroup defined as string, appended to blockGroup_ID (converted to string) and .csv to create file name
	string FileName_str = "BlockGroupDaily" + to_string(blockGroup_ID) + ".csv";
	ofstream outfile(OutputDirectory_String + FileName_str, ios::app);

	if (!outfile.good()) {
		cout << "Warning: Output folder does not exist or an output file in that folder is open and cannot be overwritten." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Correction: Confirm the HydroPlusConfig.xml parameter OutputFolder_Path and files in that folder are closed." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//HeatIndex_Lin_blockGroupAverageDaily_K (K) terms are from Lin et al. (2009) calculated with daily average Tair_K and Tdew_K
	//Note: HeatIndex_Lin_blockGroupAverage_K (K) by contrast uses average of hourly heat index values, within 0.1% of HeatIndex_Lin_blockGroupAverageDaily_K 
	//Note: Tair_blockGroupAverage_K divided by SimulationDay_Duration_sec to complete timeStep weighting of variable, initiated in BlockGroupHandler.cpp
	double HeatIndex_Lin_blockGroupAverageDaily_K = -2.653 + 0.994 * ((input->by_BlockGroup(blockGroup_ID, "Tair_blockGroupAverage_K") / input->SimulationDay_Duration_sec) - 273.15) + 0.0153 * pow(((input->by_BlockGroup(blockGroup_ID, "Tdew_blockGroupAverage_K") / input->SimulationDay_Duration_sec) - 273.15), 2);

	//BlockGroupDaily written, with averages (e.g., Tair_blockGroupAverage_K) divided by SimulationDay_Duration_sec, updated in BLockGroupHandler
	//Note: Tair_blockGroupAverage_K divided by SimulationDay_Duration_sec to complete timeStep weighting of variable, initiated in BlockGroupHandler.cpp
	outfile << date << ","
		<< input->by_BlockGroup(blockGroup_ID, "Tair_blockGroupMax_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "Tair_blockGroupMin_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "Tair_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		<< input->by_BlockGroup(blockGroup_ID, "Tdew_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		<< input->by_BlockGroup(blockGroup_ID, "HeatIndex_blockGroupMax_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "HeatIndex_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		<< input->by_BlockGroup(blockGroup_ID, "Humidex_blockGroupMax_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "Humidex_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		<< input->by_BlockGroup(blockGroup_ID, "Twetbulbglobe_blockGroupMax_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "Twetbulbglobe_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		<< input->by_BlockGroup(blockGroup_ID, "Twetbulbnatural_blockGroupMax_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "Twetbulbnatural_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		<< input->by_BlockGroup(blockGroup_ID, "Twetbulbpsychrometric_blockGroupMax_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "Twetbulbpsychrometric_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		<< input->by_BlockGroup(blockGroup_ID, "UTCI_blockGroupMax_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "UTCI_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		//Commented out but retained in case needed by future study, 2025
		//<< input->by_BlockGroup(blockGroup_ID, "HeatIndex_Lin_blockGroupMax_K") << ","
		//<< input->by_BlockGroup(blockGroup_ID, "HeatIndex_Lin_blockGroupAverage_K") / input->SimulationDay_Duration_sec << ","
		//<< (HeatIndex_Lin_blockGroupAverageDaily_K + 273.15) << ","
		<< input->by_BlockGroup(blockGroup_ID, "WindChill_blockGroupMin_K") << ","
		<< input->by_BlockGroup(blockGroup_ID, "WindChill_blockGroupAverage_K") / input->SimulationDay_Duration_sec << endl;

	outfile.close();
}