#include "BlockGroupHandler.h"
#include <algorithm>

void BlockGroupHandler::AddToBlockGroup(Inputs *input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID)
{
	bool folderIsInABlockGroupRange = PrintBlockGroupCheck(MapPixel_ID, input);
	if (input->BlockGroup_ID[MapPixel_ID] == Inputs::NODATA_code) {
		return;
	}
	if (input->TemperatureCalculationParams["PrintAllBlockGroups"] == 1 || folderIsInABlockGroupRange) {
		AddDEMLocationValuesToBlockGroup(input, beC, input->BlockGroup_ID[MapPixel_ID], MapPixel_ID, DataFolder_ID);
	}
}

//Note: For model = SpatialTemperatureHydro, DataDrawer_ID has maximum = maxRows * maxCols, and DataFolder_ID typically has maximum = 1
bool BlockGroupHandler::PrintBlockGroupCheck(int DataDrawer_ID, Inputs* input)
{
	for (int iter2 = 0; iter2 < input->StartingBlockIdsToPrintTempOutput.size(); ++iter2) {
		if (input->BlockGroup_ID[DataDrawer_ID] >= input->StartingBlockIdsToPrintTempOutput[iter2]
			&& input->BlockGroup_ID[DataDrawer_ID] <= input->EndingBlockIdsToPrintTempOutput[iter2]) {

			return true;
		}
	}
	return false;
}

void BlockGroupHandler::TimeStepProcessing(Inputs *input, int timeStep)
{
	int Time_HH_int;
	//Flag_TS_LastForDay initiated to false
	bool Flag_TS_LastForDay = false;
	//Year_YYYY = input->SimulationDate_Output_GD[timeStep] / 10000; determine year for Month_Year_to_Days_in_Month
	int Year_YYYY = input->SimulationDate_Output_GD[timeStep] / 10000;
	///Flag_LeapYear = ((Year_YYYY % 4 == 0 && Year_YYYY % 100 != 0) || (Year_YYYY % 400 == 0))
	//Note: Leap years are divisible by 4 and not divisible by 100, or divisible by 400
	bool Flag_LeapYear = ((Year_YYYY % 4 == 0 && Year_YYYY % 100 != 0) || (Year_YYYY % 400 == 0));
	//Days_in_Year = Flag_LeapYear ? 366 : 365; determines days based on ternary
	int Days_in_Year = Flag_LeapYear ? 366 : 365;

	//Date_YYYYMMDDHH_str is string created from SimulationDate_Output_GDH (YYYYMMDDHH) of timeStep, vector defined in readWeatherFile
	string Date_YYYYMMDDHH_str = to_string(input->SimulationDate_Output_GDH[timeStep]);
	//Date_YYYYMMDD_str is string created from SimulationDate_Output_GD (YYYYMMDD) of timeStep, vector defined in readWeatherFile
	string Date_YYYYMMDD_str = to_string(input->SimulationDate_Output_GD[timeStep]);
	//Time_HH_int (HH) is integer extracted from Date_YYYYMMDDHH_str string in 8th position, 2 digits
	Time_HH_int = atoi(Date_YYYYMMDDHH_str.substr(8, 2).c_str());

	//Day_JulianDay is SimulationDate_Output_GD[timeStep] Julian Day extracted from Inputs::Date_YYYYMMDD_to_JulianDay function
	int Day_JulianDay = input->Date_YYYYMMDD_to_JulianDay(input->SimulationDate_Output_GD[timeStep]);
	//Day_JulianDay_TS_Next initiated as SimulationDate_Output_GD[timeStep] Julian Day extracted from Inputs::Date_YYYYMMDD_to_JulianDay function
	int Day_JulianDay_TS_Next = input->Date_YYYYMMDD_to_JulianDay(input->SimulationDate_Output_GD[timeStep]);
	//Day_JulianDay_TS_Prior initiated as SimulationDate_Output_GD[timeStep] Julian Day extracted from Inputs::Date_YYYYMMDD_to_JulianDay function
	int Day_JulianDay_TS_Prior = input->Date_YYYYMMDD_to_JulianDay(input->SimulationDate_Output_GD[timeStep]);

	//If timeStep + 1 does not exceed SimulationDate_Output_GD.size then extract the next date in vector
	if (timeStep + 1 < input->SimulationDate_Output_GD.size()) {
		//Day_JulianDay_TS_Next is SimulationDate_Output_GD[timeStep + 1] Julian Day extracted from Inputs::Date_YYYYMMDD_to_JulianDay function
		Day_JulianDay_TS_Next = input->Date_YYYYMMDD_to_JulianDay(input->SimulationDate_Output_GD[timeStep + 1]);
	}
	//Else timeStep + 1 exceeds SimulationDate_Output_GD.size then at last time step for simulation and hence day
	else {
		//Flag_TS_LastForDay set to true
		Flag_TS_LastForDay = true;
	}

	//SimulationDay_Duration_sec incremented by input->SimulationTimeStep_Duration_sec[timeStep] to contain all seconds in simulated day
	//Note: Some simulations will only simulate partial days, some will use irregular time step durations
	input->SimulationDay_Duration_sec = input->SimulationDay_Duration_sec + input->SimulationTimeStep_Duration_sec[timeStep];

	//If (timeStep - 1 >= 0) then vector index exists
	if ((timeStep - 1) >= 0) {
		//Day_JulianDay_TS_Prior is SimulationDate_Output_GD[timeStep - 1] Julian Day extracted from Inputs::Date_YYYYMMDD_to_JulianDay function
		Day_JulianDay_TS_Prior = input->Date_YYYYMMDD_to_JulianDay(input->SimulationDate_Output_GD[timeStep - 1]);
	}

	//If (Day_JulianDay_TS_Next - Day_JulianDay == 1) then 1 day separates timeSteps and this is last time step of day, in same year
	//Note: Flag_TS_LastForDay initiated to false
	if (Day_JulianDay_TS_Next - Day_JulianDay == 1) {
		//Flag_TS_LastForDay set to true
		Flag_TS_LastForDay = true;
	}
	//Else If (Day_JulianDay - Day_JulianDay_TS_Next >= (Days_in_Year - 1)) then New Year's Eve separates timeSteps and this is last time step in year
	else if (Day_JulianDay - Day_JulianDay_TS_Next >= (Days_in_Year - 1)) {
		//Flag_TS_LastForDay set to true
		Flag_TS_LastForDay = true;
	}

	//If (Day_JulianDay - Day_JulianDay_TS_Prior == 1) then 1 day separates timeSteps and this is first time step of new day, in same year
	if (Day_JulianDay - Day_JulianDay_TS_Prior == 1) {
		//SimulationDay_Duration_sec is input->SimulationTimeStep_Duration_sec[timeStep]
		input->SimulationDay_Duration_sec = input->SimulationTimeStep_Duration_sec[timeStep];
	}
	//Else If (Day_JulianDay_TS_Prior - Day_JulianDay >= (Days_in_Year - 1)) then New Year's Eve separates timeSteps and this is first time step of a new year
	else if (Day_JulianDay_TS_Prior - Day_JulianDay >= (Days_in_Year - 1)) {
		//SimulationDay_Duration_sec is input->SimulationTimeStep_Duration_sec[timeStep]
		input->SimulationDay_Duration_sec = input->SimulationTimeStep_Duration_sec[timeStep];
	}


	//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 auto blockGroupIter : BlockGroupIDToItsFolder; loop through all Unique_BlockGroup_IDs
	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 timeStep is 0 then
		if (timeStep == 0) {
			//ResetBlockGroupVar_to_zeroDailyValues function called to clear daily vectors
			ResetBlockGroupVar_to_zeroDailyValues(input, blockGroupNum);
		}

		//if (PrintBlockGroupDailyAndHourly) true then call writeToBlockGroupHourlyFile
		if (Flag_needHourlyBlockGroupFiles) {
			//Call writeToBlockGroupHourlyFile(input->SimulationStringParams["OutputFolder_Path"], Date_YYYYMMDDHH_str, blockGroupNum)
			TemperatureOutputWriter::writeToBlockGroupHourlyFile(input, input->SimulationStringParams["OutputFolder_Path"], Date_YYYYMMDDHH_str, blockGroupNum);
		}

		//DailyCal function accumulates blockgroup values into daily statistics; perhaps send counter
		//Note: Each call updates the MAX value as the current time step or prior
		DailyCal(input, blockGroupNum, timeStep);

		//ResetBlockGroupVar_to_zero function called to clear hourly vectors
		ResetBlockGroupVar_to_zero(input, blockGroupNum);

		//If Flag_TS_LastForDay and Flag_needDailyBlockGroupFiles true, then call writeToBlockGroupDailyFile
		//Note: writeToBlockGroupDailyFile is called for any last hour of the simulation, even if a day has not fully completed at hour 23
		if (Flag_TS_LastForDay && Flag_needDailyBlockGroupFiles) {
			//TemperatureOutputWriter::writeToBlockGroupDailyFile function called
			TemperatureOutputWriter::writeToBlockGroupDailyFile(input, input->SimulationStringParams["OutputFolder_Path"],
				Date_YYYYMMDD_str, blockGroupNum, timeStep);

			//ResetBlockGroupVar_to_zeroDailyValues function called to clear daily vectors
			ResetBlockGroupVar_to_zeroDailyValues(input, blockGroupNum);
		}
	}
}

BlockGroupHandler::BlockGroupHandler(Inputs* input)
{
	// Initialize mins to a limit for all BG drawers
	double limit = 10000;
	for (size_t i = 0; i < input->Unique_BlockGroup_IDs.size(); ++i) {
		const int bg_ID = input->Unique_BlockGroup_IDs[i];
		input->by_BlockGroup(bg_ID, "WindChill_blockGroupMin_K") = limit;
		input->by_BlockGroup(bg_ID, "Tair_blockGroupMin_K") = limit;
		// seed identifiers
		input->by_BlockGroup(bg_ID, "blockGroupIndex") = bg_ID;
		// pixels Count_of this bgID
		input->by_BlockGroup(bg_ID, "pixelsCount_of_bgID") = 0.0;
	}
	
	/*unordered_set<int> bg_print_allow;
	// Case 1: print ALL block groups
	if (input->TemperatureCalculationParams["PrintAllBlockGroups"] != 0) {
		for (int bg_ID : input->Unique_BlockGroup_IDs)
			bg_print_allow.insert(bg_ID);
	}
	// Case 2: print only specified ranges
	else {
		for (size_t r = 0; r < input->StartingBlockIdsToPrintTempOutput.size(); ++r) {

			int start = input->StartingBlockIdsToPrintTempOutput[r];
			int end = input->EndingBlockIdsToPrintTempOutput[r];

			for (int bg_ID : input->Unique_BlockGroup_IDs) {
				if (bg_ID >= start && bg_ID <= end)
					bg_print_allow.insert(bg_ID);
			}
		}
	}*/


	// Walk pixels once to (a) count membership, (b) optionally skip nodata
	for (int pix = 0; pix < (int)input->BlockGroup_ID.size(); ++pix) {
		const int bgID = input->BlockGroup_ID[pix];

		bool CountBlockGroup = (input->TemperatureCalculationParams["PrintAllBlockGroups"] && (bgID != Inputs::NODATA_code));

		if (!CountBlockGroup) {
			CountBlockGroup = PrintBlockGroupCheck(pix, input);
		}
		if (!CountBlockGroup) {
			continue;
		}

		// Increment per-pixel membership count
		input->by_BlockGroup(bgID, "pixelsCount_of_bgID") += 1.0;
	}
}

/*BlockGroupHandler::~BlockGroupHandler()
{
	for (auto b : BlockGroupIDToItsFolder) {
		delete b.second;
	}

}*/



void BlockGroupHandler::AddDEMLocationValuesToBlockGroup(Inputs* input, CompactRagged* beC, int blockGroupID, int MapPixel_ID, int DataFolder_ID)
{
	input->by_BlockGroup(blockGroupID, "ImpNR_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "ImpNR_W_p_m2");
	input->by_BlockGroup(blockGroupID, "ImpH_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "ImpH_W_p_m2");
	input->by_BlockGroup(blockGroupID, "ImpLE_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "ImpLE_W_p_m2");
	input->by_BlockGroup(blockGroupID, "TreeNR_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "TreeNR_W_p_m2");
	input->by_BlockGroup(blockGroupID, "TreeH_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "TreeH_W_p_m2");
	input->by_BlockGroup(blockGroupID, "TreeLE_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLE_W_p_m2");
	input->by_BlockGroup(blockGroupID, "TreeLEE_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLEE_W_p_m2");
	input->by_BlockGroup(blockGroupID, "TreeLET_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "TreeLET_W_p_m2");
	input->by_BlockGroup(blockGroupID, "SVegNR_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "SVegNR_W_p_m2");
	input->by_BlockGroup(blockGroupID, "SVegH_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "SVegH_W_p_m2");
	input->by_BlockGroup(blockGroupID, "SVegLE_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLE_W_p_m2");
	input->by_BlockGroup(blockGroupID, "SVegLEE_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLEE_W_p_m2");
	input->by_BlockGroup(blockGroupID, "SVegLET_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "SVegLET_W_p_m2");
	input->by_BlockGroup(blockGroupID, "SoilNR_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "SoilNR_W_p_m2");
	input->by_BlockGroup(blockGroupID, "SoilH_W_p_m2") +=beC->by_key(MapPixel_ID, DataFolder_ID, "SoilH_W_p_m2");
	input->by_BlockGroup(blockGroupID, "SoilLE_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "SoilLE_W_p_m2");
	input->by_BlockGroup(blockGroupID, "WaterNR_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "WaterNR_W_p_m2");
	input->by_BlockGroup(blockGroupID, "WaterH_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "WaterH_W_p_m2");
	input->by_BlockGroup(blockGroupID, "WaterLE_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "WaterLE_W_p_m2");
	input->by_BlockGroup(blockGroupID, "H_total_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "H_total_W_p_m2");
	input->by_BlockGroup(blockGroupID, "LE_total_W_p_m2") += beC->by_key(MapPixel_ID, DataFolder_ID, "LE_total_W_p_m2");
	input->by_BlockGroup(blockGroupID, "Tair_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K");
	input->by_BlockGroup(blockGroupID, "Tdew_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "Tdew_K");
	input->by_BlockGroup(blockGroupID, "HeatIndex_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "HeatIndex_K");
	input->by_BlockGroup(blockGroupID, "Humidex_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "Humidex_K");
	input->by_BlockGroup(blockGroupID, "Twetbulbglobe_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "Twetbulbglobe_K");
	input->by_BlockGroup(blockGroupID, "Twetbulbnatural_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "Twetbulbnatural_K");
	input->by_BlockGroup(blockGroupID, "Twetbulbpsychrometric_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "Twetbulbpsychrometric_K");
	input->by_BlockGroup(blockGroupID, "UTCI_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "UTCI_K");
	input->by_BlockGroup(blockGroupID, "WindChill_K") += beC->by_key(MapPixel_ID, DataFolder_ID, "WindChill_K");
	input->by_BlockGroup(blockGroupID, "AbsHumidity_kg_p_m3") += beC->by_key(MapPixel_ID, DataFolder_ID, "AbsHumidity_kg_p_m3");
	input->by_BlockGroup(blockGroupID, "Tair_mesoScale_K") +=  beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_mesoScale_final_K");
	input->by_BlockGroup(blockGroupID, "AbsHumidity_mesoScale_kg_p_m3") += beC->by_key(MapPixel_ID, DataFolder_ID, "AbsHumidity_mesoScale_final_kg_p_m3");
}

void BlockGroupHandler::ResetBlockGroupVar_to_zeroDailyValues(Inputs* input, int blockGroupID)
{
	const double limit = 10000.0;

	// List of all daily-reset variables and their initial values
	const vector<pair<string, double>> vars = {
		{"Tair_blockGroupMax_K", 0.0},
		{"Tair_blockGroupMin_K", limit},
		{"Tair_blockGroupAverage_K", 0.0},
		{"Tdew_blockGroupAverage_K", 0.0},
		{"HeatIndex_blockGroupMax_K", 0.0},
		{"HeatIndex_blockGroupAverage_K", 0.0},
		{"Humidex_blockGroupMax_K", 0.0},
		{"Humidex_blockGroupAverage_K", 0.0},
		{"Twetbulbglobe_blockGroupMax_K", 0.0},
		{"Twetbulbglobe_blockGroupAverage_K", 0.0},
		{"Twetbulbnatural_blockGroupMax_K", 0.0},
		{"Twetbulbnatural_blockGroupAverage_K", 0.0},
		{"Twetbulbpsychrometric_blockGroupMax_K", 0.0},
		{"Twetbulbpsychrometric_blockGroupAverage_K", 0.0},
		{"UTCI_blockGroupMax_K", 0.0},
		{"UTCI_blockGroupAverage_K", 0.0},
		{"HeatIndex_Lin_blockGroupMax_K", 0.0},
		{"HeatIndex_Lin_blockGroupAverage_K", 0.0},
		{"WindChill_blockGroupMin_K", limit},
		{"WindChill_blockGroupAverage_K", 0.0}
	};

	// Reset all listed variables to their default daily values
	for (const auto& [name, value] : vars) {
		input->by_BlockGroup(blockGroupID, name) = value;
	}
}

void BlockGroupHandler::ResetBlockGroupVar_to_zero(Inputs* input, int blockGroupID)
{
	// List of all BG variables to reset
	const vector<string> vars = {
		"ImpNR_W_p_m2", "ImpH_W_p_m2", "ImpLE_W_p_m2",
		"TreeNR_W_p_m2", "TreeH_W_p_m2", "TreeLE_W_p_m2",
		"TreeLEE_W_p_m2", "TreeLET_W_p_m2",
		"SVegNR_W_p_m2", "SVegH_W_p_m2", "SVegLE_W_p_m2",
		"SVegLEE_W_p_m2", "SVegLET_W_p_m2",
		"SoilNR_W_p_m2", "SoilH_W_p_m2", "SoilLE_W_p_m2",
		"WaterNR_W_p_m2", "WaterH_W_p_m2", "WaterLE_W_p_m2",
		"H_total_W_p_m2", "LE_total_W_p_m2",
		"Tdew_K", "Tair_K", "HeatIndex_K", "Humidex_K",
		"Twetbulbglobe_K", "Twetbulbnatural_K", "Twetbulbpsychrometric_K",
		"UTCI_K", "WindChill_K",
		"AbsHumidity_kg_p_m3", "Tair_mesoScale_final_K",
		"AbsHumidity_mesoScale_final_kg_p_m3"
	};

	// Reset each variable to 0 (auto-creates column if missing)
	for (const auto& name : vars) {
		input->by_BlockGroup(blockGroupID, name) = 0.0;
	}
}
//DailyCal function accumulates blockgroup values into daily statistics
//Note: Called each timestep to find average value for blockgroup and update maximum value for day
void BlockGroupHandler::DailyCal(Inputs* input, int blockGroupID, int timeStep)
{
	//Tair_K (C) is derived as quotient of blockGroupID summed Tair_K (K) and blockGroupID folderCount
	double Tair_K = input->by_BlockGroup(blockGroupID, "Tair_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID");
	//Tdew_K (C) is derived as quotient of blockGroupID summed Tdew_K (K) and blockGroupID folderCount
	double Tdew_K = input->by_BlockGroup(blockGroupID, "Tdew_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID");
	//Tair_C (C) is derived as quotient of blockGroupID summed Tair_K (K) and blockGroupID folderCount, converted to degrees C by subtracting 273.15
	double Tair_C = input->by_BlockGroup(blockGroupID, "Tair_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") - 273.15;
	//Tdew_C (C) is derived as quotient of blockGroupID summed Tdew_K (K) and blockGroupID folderCount, converted to degrees C by subtracting 273.15
	double Tdew_C = input->by_BlockGroup(blockGroupID, "Tdew_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") - 273.15;

	//HeatIndex_Lin_C defined using Eq text on bottom of pg 739 in Lin et al. (2009) 
	//References: 
	//Lin, S., M. Luo, R. J. Walker, X. Liu, S. A. Hwang and R. Chinery (2009). "Extreme high temperatures and hospital admissions for respiratory and cardiovascular diseases." Epidemiology 20(5): 738-746.
	//Basu, R., W.Y.Feng and B.D.Ostro(2008). "Characterizing temperature and mortality in nine California counties." Epidemiology 19(1) : 138 - 145.
	double HeatIndex_Lin_C = -2.653 + (0.994 * Tair_C) + 0.0153 * pow(Tdew_C,2);

	//Tair_blockGroupAverage_K is summed to contain all values from Blockgroup, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "Tair_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "Tair_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];
	
	//Tdew_blockGroupAverage_K is summed to contain all values from Blockgroup, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "Tdew_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "Tdew_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];
	
	//HeatIndex_blockGroupAverage_K is summed to contain all values from Blockgroup for day, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "HeatIndex_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "HeatIndex_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];

	//Humidex_blockGroupAverage_K is summed to contain all values from Blockgroup for day, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "Humidex_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "Humidex_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];

	//Twetbulbglobe_blockGroupAverage_K is summed to contain all values from Blockgroup, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "Twetbulbglobe_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "Twetbulbglobe_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];

	//Twetbulbnatural_blockGroupAverage_K is summed to contain all values from Blockgroup, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "Twetbulbnatural_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "Twetbulbnatural_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];

	//Twetbulbpsychrometric_blockGroupAverage_K is summed to contain all values from Blockgroup, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "Twetbulbpsychrometric_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "Twetbulbpsychrometric_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];

	//UTCI_blockGroupAverage_K is summed to contain all values from Blockgroup, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "UTCI_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "UTCI_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];

	//WindChill_blockGroupAverage_K is summed to contain all values from Blockgroup, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "WindChill_blockGroupAverage_K") += input->by_BlockGroup(blockGroupID, "WindChill_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID") * input->SimulationTimeStep_Duration_sec[timeStep];

	//HeatIndex_Lin_blockGroupAverage_K is summed to contain all values from Blockgroup, weighted by * SimulationTimeStep_Duration_sec in this routine
	//Note: Weighting of daily averages divides by SimulationDay_Duration_sec in TemperatureOutputWriter to handle daily variation of timeStep count or duration
	input->by_BlockGroup(blockGroupID, "HeatIndex_Lin_blockGroupAverage_K") += (HeatIndex_Lin_C + 273.15) * input->SimulationTimeStep_Duration_sec[timeStep];

	//Tair_blockGroupMax_K is MAX of Tair_blockGroupMax_K or Tair_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is larger than the prior maximum value 
	input->by_BlockGroup(blockGroupID, "Tair_blockGroupMax_K") = MAX(input->by_BlockGroup(blockGroupID, "Tair_blockGroupMax_K"), input->by_BlockGroup(blockGroupID, "Tair_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));

	//Tair_blockGroupMin_K is MIN of Tair_blockGroupMin_K or Tair_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is lesser than the prior minimum value 
	input->by_BlockGroup(blockGroupID, "Tair_blockGroupMin_K") = MIN(input->by_BlockGroup(blockGroupID, "Tair_blockGroupMin_K"), input->by_BlockGroup(blockGroupID, "Tair_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));

	//HeatIndex_blockGroupMax_K is MAX of HeatIndex_blockGroupMax_K or HeatIndex_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is larger than the prior maximum value 
	input->by_BlockGroup(blockGroupID, "HeatIndex_blockGroupMax_K") = MAX(input->by_BlockGroup(blockGroupID, "HeatIndex_blockGroupMax_K"), input->by_BlockGroup(blockGroupID, "HeatIndex_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));

	//Humidex_blockGroupMax_K is MAX of Humidex_blockGroupMax_K or Humidex_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is larger than the prior maximum value 
	input->by_BlockGroup(blockGroupID, "Humidex_blockGroupMax_K") = MAX(input->by_BlockGroup(blockGroupID, "Humidex_blockGroupMax_K"), input->by_BlockGroup(blockGroupID, "Humidex_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));

	//Twetbulbglobe_blockGroupMax_K is MAX of Twetbulbglobe_blockGroupMax_K or Twetbulbglobe_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is larger than the prior maximum value 
	input->by_BlockGroup(blockGroupID, "Twetbulbglobe_blockGroupMax_K") = MAX(input->by_BlockGroup(blockGroupID, "Twetbulbglobe_blockGroupMax_K"), input->by_BlockGroup(blockGroupID, "Twetbulbglobe_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));

	//Twetbulbnatural_blockGroupMax_K is MAX of Twetbulbnatural_blockGroupMax_K or Twetbulbnatural_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is larger than the prior maximum value 
	input->by_BlockGroup(blockGroupID, "Twetbulbnatural_blockGroupMax_K") = MAX(input->by_BlockGroup(blockGroupID, "Twetbulbnatural_blockGroupMax_K"), input->by_BlockGroup(blockGroupID, "Twetbulbnatural_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));

	//Twetbulbpsychrometric_blockGroupMax_K is MAX of Twetbulbpsychrometric_blockGroupMax_K or Twetbulbpsychrometric_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is larger than the prior maximum value 
	input->by_BlockGroup(blockGroupID, "Twetbulbpsychrometric_blockGroupMax_K") = MAX(input->by_BlockGroup(blockGroupID, "Twetbulbpsychrometric_blockGroupMax_K"), input->by_BlockGroup(blockGroupID, "Twetbulbpsychrometric_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));

	//UTCI_blockGroupMax_K is MAX of UTCI_blockGroupMax_K or UTCI_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is larger than the prior maximum value 
	input->by_BlockGroup(blockGroupID, "UTCI_blockGroupMax_K") = MAX(input->by_BlockGroup(blockGroupID, "UTCI_blockGroupMax_K"), input->by_BlockGroup(blockGroupID, "UTCI_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));

	//HeatIndex_Lin_blockGroupMax_K is MAX of HeatIndex_Lin_blockGroupMax_K or HeatIndex_Lin_C / folderCount
	//Note: Lin provides separate method for HeatIndex that was used during a USFS NUCFAC research project with Sinha et al (2022)
	input->by_BlockGroup(blockGroupID, "HeatIndex_Lin_blockGroupMax_K") = MAX(input->by_BlockGroup(blockGroupID, "HeatIndex_Lin_blockGroupMax_K"), (HeatIndex_Lin_C + 273.15));
	
	//WindChill_blockGroupMin_K is MIN of WindChill_blockGroupMin_K or WindChill_K / folderCount
	//Note: Updated each timestep to determine if the new blockbgroup average is smaller than the prior minimum value 
	input->by_BlockGroup(blockGroupID, "WindChill_blockGroupMin_K") = MIN(input->by_BlockGroup(blockGroupID, "WindChill_blockGroupMin_K"), input->by_BlockGroup(blockGroupID, "WindChill_K") / input->by_BlockGroup(blockGroupID, "pixelsCount_of_bgID"));
}
