#include "RoutingDiffusion.h"
#include <math.h>
#include <fstream>
#include <vector>
#include <iostream>
#include <numeric>
#include <string>

//Theory: The 2-parameter diffusion model has 2 time-based parameters, alpha and beta, representing watershed scale flow diffusivity and flow celerity
//Note: The 2-parameter diffusion hydrograph runoff routing developed by Yang and Endreny (2013) and Yang et al. (2015), replacing the time-area hydrograph developed by Wang et al. (2008)
//Note: This 2-parameter diffusion model has excellent representation of natural hydrograph characteristics, including the inflection points in the rising and falling limbs and time to peak.
//Note: The 2 parameters alpha and beta represent watershed and storm influences on runoff and hydrograph features
//Note: If surface flow and subsurface flow have same time constants alpha and beta, the model can simulate discharge with mixed surface flow and subsurface flow.
//Note: If surface flow and subsurface flow have different time constants alpha and beta, a parallel model that consists of two-parameter models can be constructed to simulate the hydrograph
//Note: The 2-parameter surface flow diffusion hydrograph model is used for PerviousFlow (PAQ_m) and Impervious Flow (DCIAQ_m).
//Note; The 1-parameter subsurface flow diffusion hydrograph model is used for Baseflow (SubQ_m), based on Criss and Winston (2008a).
//Note: The time input, j, is used in calculating the discharge from the precip (pulse input) at time j in the future

//Yang, Yang, Theodore A. Endreny. (2013) "Watershed hydrograph model based on surface flow diffusion", Water Resources Research Vol 49:507-516.
//Yang, Y., Endreny, T.A., & Nowak, D.J. (2015).Simulating the effect of flow path roughness to examine how green infrastructure restores urban runoff timing and magnitude.Urban Forestry& Urban Greening, 14(2), 361 - 367. doi:https://doi.org/10.1016/j.ufug.2015.03.004
//Wang, J., Endreny, T.A, & Nowak, D.J. (2008) "Mechanistic Simulation of Tree Effects in an Impervious Water Balance Model1." JAWRA Journal of the American Water Resources Association Vol(44)1: 75-85 
//Criss, and Winston (2008a). Properties of a diffusive hydrograph and the interpretation of its single parameter. Mathematical Geosciences, 40(3), 313-325. doi:10.1007/s11004-008-9145-9
//Criss andWinston (2008b). Discharge predictions of a rainfall-driven theoretical hydrograph compared to common models and observed data. Water Resources Research, 44(10). doi:10.1029/2007wr006415
 
//Note: Determining alpha and beta, from Yang and Endreny (2013):
/*
	A visual demonstration of the influence of alpha and beta on the normalized hydrograph is shown in Figure 4 of Yang and Endreny (2013).
	Setting beta = 1, and varying alpha from 0.1 and 10, alpha and beta together determine the rising limb, t_max, and falling limb when (Figure 4a);
	The smaller alpha, the faster the rising and falling limb and the smaller t_max.
	Setting beta = 1 and alpha >10, beta determines the hydrograph shape and t_max = 2(beta/3).
	Setting alpha = 1 and varying beta from 0.1 to 100, the normalized hydrographs series indicate that beta always contributes to t_max (Figure 4b).
	The larger beta, the longer the hydrograph duration.
	If alpha goes to infinity, kinematic flow celerity = 0 and the 2-parameter model takes the same form as the 1-parameter model of Criss and Winston (2008b)
*/

//RoutingDiffusion::RoutingRunoff function hydrologically routes runoff depths from overland impervious, pervious, and water areas, and subsurface pervious area. 
//Note: Incoming runoff depths, PAQi_m, DCIAQi_m, WAQi_m, SubQi_m were transformed to folder area (e.g. bulk area) in FlowSumCalc.cpp or earlier. 
void RoutingDiffusion::RoutingRunoff(Inputs* input, int timeStep)
{
	//routing_flag is option to create and use alternative routing or no routing; flag = 1 use 2-parameter diffusion hydrograph, flag = 0 skip routing
	bool routing_flag = 1;	
	//timeStep_routed is timestep + time steps advanced by convolution for loop
	int timeStep_routed = 0;	
	auto const stride = 3;

	//TotalTimeSteps defined from SimulationTimePeriod_timeSteps, based on HydroPlusConfig StartDate_YYYYMMDD and StopDate_YYYYMMDD
	int TotalTimeSteps = input->SimulationTimePeriod_timeSteps;

	//If timeStep equals 0, or first time step, then dimensionalize vectors
	if (timeStep == 0) {
		//Dimensionalize size for catchment runoff depth (m) vectors as time distributed vector values
		//Note: These input-> vector double variables were created in Inputs.h
		//Pervious area runoff, including runon from impervious area (m)
		input->PAQi_Routed_m.insert(input->PAQi_Routed_m.end(), TotalTimeSteps, 0.0);
		//Directly connected impervious area runoff (m)
		input->DCIAQi_Routed_m.insert(input->DCIAQi_Routed_m.end(), TotalTimeSteps, 0.0);
		//Water area runoff, which can represent hydrologically connected waters (m)
		input->WAQi_Routed_m.insert(input->WAQi_Routed_m.end(), TotalTimeSteps, 0.0);
		//Surface runoff (m)
		input->SurfQi_Routed_m.insert(input->SurfQi_Routed_m.end(), TotalTimeSteps, 0.0);
		//Subsurface runoff (m)
		input->SubQi_Routed_m.insert(input->SubQi_Routed_m.end(), TotalTimeSteps, 0.0);
		//Total runoff (m)
		input->Qi_Routed_m.insert(input->Qi_Routed_m.end(), TotalTimeSteps, 0.0);
	}

	//Define pervious, impervious, and baseflow runoff terms using unrouted runoff depths from FLowSumCalc.cpp
	double PAQ_m = input->PAQi_m[timeStep];
	double DCIAQ_m = input->DCIAQi_m[timeStep];
	double WAQ_m = input->WAQi_m[timeStep];
	double SubQ_m = input->SubQi_m[timeStep];

	if (routing_flag == 1) {

		//PAQ_RT_A_h (h) is alpha time parameter of the two parameter surface flow diffusion hydrograph model
		double PAQ_RT_A_h = input->InputXml["Routing_PerviousArea_Alpha_h"];		
		//PAQ_RT_B_h (h) is beta time parameter of the two parameter surface flow diffusion hydrograph model
		double PAQ_RT_B_h = input->InputXml["Routing_PerviousArea_Beta_h"];
		//DCIAQ_RT_A_h (h) is alpha time parameter of the two parameter surface flow diffusion hydrograph model
		double DCIAQ_RT_A_h = input->InputXml["Routing_ImperviousArea_Alpha_h"];
		//DCIAQ_RT_B_h  (h) is beta time parameter of the two parameter surface flow diffusion hydrograph model
		double DCIAQ_RT_B_h = input->InputXml["Routing_ImperviousArea_Beta_h"];
		//SSQ_RT_h (h) is beta time parameter of the one parameter subsurface flow diffusion hydrograph model
		double SSQ_RT_h = input->InputXml["Routing_Subsurface_Beta_h"];

		//Note: t_max is computed in order to estimate the upper duration of the hydrograph for use in a convolution for loop
		//Time_PA_Max_h from Equation 6 of Yang et al. (2015) and Eq 3 of Yang et al. (2013), generating t_max time 
		double Time_PA_Max_h = (-0.75 * PAQ_RT_A_h) + (0.75 * PAQ_RT_A_h) * sqrt(1 + (16.0 / 9.0) * (PAQ_RT_B_h / PAQ_RT_A_h));
		//Time_DCIA_Max_h from Equation 6 of Yang et al. (2015) and Eq 3 of Yang et al. (2013), generating t_max time 
		double Time_DCIA_Max_h = (-0.75 * DCIAQ_RT_A_h) + (0.75 * DCIAQ_RT_A_h) * sqrt(1 + (16.0 / 9.0) * (DCIAQ_RT_B_h / DCIAQ_RT_A_h));
		//Time_Sub_Max_h from text below Equation 10 of Yang et al. (2013), where SSQ_RT_h is beta coefficient, generating t_max time 
		double Time_Sub_Max_h = (2.0 / 3.0) * SSQ_RT_h;

		//Time_PA_upper_h (hour, not time steps) is estimator of how far into future to route hydrograph and complete the falling or receding limb
		//Note: Yang derived upper time limits as 10*(b * 7.435 + 4.73), where any decimals are dropped when converted to integer variable
		int Time_PA_upper_h = static_cast<int>(10 * ((Time_PA_Max_h * 7.435) + 4.73));
		int Time_DCIA_upper_h = static_cast<int>(10 * ((Time_DCIA_Max_h * 7.435) + 4.73));
		int Time_Sub_upper_h = static_cast<int>(10 * ((Time_Sub_Max_h * 7.435) + 4.73));

		//Variables advance time through routed hydrograph
		double PAQ_it = 0.1;
		double DCIAQ_it = 0.1;
		double WAQ_it = 0.1;
		double SubQ_it = 0.1;
		//Variables sum the three temporal components of routed hydrograph
		double Delay_PAQ_TS_sum = 0.0;
		double Delay_DCIAQ_TS_sum = 0.0;
		double Delay_WAQ_TS_sum = 0.0;
		double Delay_SubQ_TS_sum = 0.0;

		//PAQ_m (m) is pervious area non-routed runoff, considered effective precipitation
		if (PAQ_m > 0.0) {
			//For loop serves as a hydrograph convolution integral, to distribute effective precipitation from time t_now as runoff at times t_future
			for (int j = 0; j < Time_PA_upper_h && timeStep + j < TotalTimeSteps; ++j) {
				//timeStep_routed is time for hydrograph, 
				timeStep_routed = timeStep + j;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, PAQ_m as Precipitation effective
				double first = PAQ_m / sqrt(3.14) * sqrt(PAQ_RT_B_h) * pow(PAQ_it, -1.5) * exp(-PAQ_RT_B_h / PAQ_it) * exp(-PAQ_it / PAQ_RT_A_h) * exp(2 * sqrt(PAQ_RT_B_h / PAQ_RT_A_h));
				//If j equals 0 time is set to zero
				if (j == 0) { PAQ_it = 0.0; }
				//time increased as part of convolution integral
				PAQ_it += 1;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, PAQ_m as Precipitation effective
				double second = PAQ_m / sqrt(3.14) * sqrt(PAQ_RT_B_h) * pow(PAQ_it, -1.5) * exp(-PAQ_RT_B_h / PAQ_it) * exp(-PAQ_it / PAQ_RT_A_h) * exp(2 * sqrt(PAQ_RT_B_h / PAQ_RT_A_h));
				//time increased as part of convolution integral
				PAQ_it += 1;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, PAQ_m as Precipitation effective
				double third = PAQ_m / sqrt(3.14) * sqrt(PAQ_RT_B_h) * pow(PAQ_it, -1.5) * exp(-PAQ_RT_B_h / PAQ_it) * exp(-PAQ_it / PAQ_RT_A_h) * exp(2 * sqrt(PAQ_RT_B_h / PAQ_RT_A_h));
				//time increased as part of convolution integral
				PAQ_it += 1;
				//Delay_PAQ_TS_sum sums three above terms to construct convolution integral
				Delay_PAQ_TS_sum = first + second + third;
				//PAQi_Routed_m vector at time timeStep_routed is stored
				input->PAQi_Routed_m[timeStep_routed] += Delay_PAQ_TS_sum;
			}
		}

		//DCIAQ_m (m) is impervious area non-routed runoff, considered effective precipitation
		if (DCIAQ_m > 0.0) {
			//For loop serves as a hydrograph convolution integral, to distribute effective precipitation from time t_now as runoff at times t_future
			for (int j = 0; j < Time_DCIA_upper_h && timeStep + j < TotalTimeSteps; ++j)	{
				//timeStep_routed is time for hydrograph, 
				timeStep_routed = timeStep + j;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, DICAQ_m as Precipitation effective
				double first = DCIAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(DCIAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / DCIAQ_it) * exp(-DCIAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				//If j equals 0 time is set to zero
				if (j == 0) { DCIAQ_it = 0; }
				//time increased as part of convolution integral
				DCIAQ_it += 1;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, DICAQ_m as Precipitation effective
				double second = DCIAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(DCIAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / DCIAQ_it) * exp(-DCIAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				//time increased as part of convolution integral
				DCIAQ_it += 1;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, DICAQ_m as Precipitation effective
				double third = DCIAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(DCIAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / DCIAQ_it) * exp(-DCIAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				//time increased as part of convolution integral
				DCIAQ_it += 1;
				//Delay_DCIAQ_TS_sum sums three above terms to construct convolution integral
				Delay_DCIAQ_TS_sum = first + second + third;
				//DCIAQi_Routed_m vector at time timeStep_routed is stored
				input->DCIAQi_Routed_m[timeStep_routed] += Delay_DCIAQ_TS_sum;
			}
		}

		//WAQ_m (m) is water area non-routed runoff, considered effective precipitation, using DCIAQ_RT_A_h and DCIAQ_RT_B_h 
		if (WAQ_m > 0.0) {
			//For loop serves as a hydrograph convolution integral, to distribute effective precipitation from time t_now as runoff at times t_future
			for (int j = 0; j < Time_PA_upper_h && timeStep + j < TotalTimeSteps; ++j) {
				//timeStep_routed is time for hydrograph, 
				timeStep_routed = timeStep + j;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, WAQ_m as Precipitation effective
				double first = WAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(WAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / WAQ_it) * exp(-WAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				//If j equals 0 time is set to zero
				if (j == 0) { WAQ_it = 0.0; }
				//time increased as part of convolution integral
				WAQ_it += 1;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, WAQ_m as Precipitation effective
				double second = WAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(WAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / WAQ_it) * exp(-WAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				//time increased as part of convolution integral
				WAQ_it += 1;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, WAQ_m as Precipitation effective
				double third = WAQ_m / sqrt(3.14) * sqrt(DCIAQ_RT_B_h) * pow(WAQ_it, -1.5) * exp(-DCIAQ_RT_B_h / WAQ_it) * exp(-WAQ_it / DCIAQ_RT_A_h) * exp(2 * sqrt(DCIAQ_RT_B_h / DCIAQ_RT_A_h));
				//time increased as part of convolution integral
				WAQ_it += 1;
				//Delay_DCIAQ_TS_sum sums three above terms to construct convolution integral
				Delay_WAQ_TS_sum = first + second + third;
				//DCIAQi_Routed_m vector at time timeStep_routed is stored
				input->WAQi_Routed_m[timeStep_routed] += Delay_WAQ_TS_sum;
			}
		}

		//SubQ_m (m) is pervious subsurface area non-routed runoff, using version of Equation 10 in Yang and Endreny (2013) based on Criss and Winston (2003)
		//Note: Consider refactor if it is important that routed runoff subsurface has greater similarity to unrouted
		if (SubQ_m > 0.0) {
			//For loop serves as a hydrograph convolution integral, to distribute effective precipitation from time t_now as runoff at times t_future
			for (int j = 0; j < Time_Sub_upper_h && timeStep + j < TotalTimeSteps; ++j) {
				//timeStep_routed is time for hydrograph, 
				timeStep_routed = timeStep + j;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, SubQ_m as Precipitation effective, dropping last 2 terms
				double first = SubQ_m / sqrt(3.14) * sqrt(SSQ_RT_h) * pow(SubQ_it, -1.5) * exp(-SSQ_RT_h / SubQ_it);
				if (j == 0) { SubQ_it = 0; }
				//time increased as part of convolution integral
				SubQ_it += 1;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, SubQ_m as Precipitation effective, dropping last 2 terms
				double second = SubQ_m / sqrt(3.14) * sqrt(SSQ_RT_h) * pow(SubQ_it, -1.5) * exp(-SSQ_RT_h / SubQ_it);
				//time increased as part of convolution integral
				SubQ_it += 1;
				//Equation 5 in Yang et al. (2015) Flow Path Roughness, with A removed in order to give runoff as depth, SubQ_m as Precipitation effective, dropping last 2 terms
				double third = SubQ_m / sqrt(3.14) * sqrt(SSQ_RT_h) * pow(SubQ_it, -1.5) * exp(-SSQ_RT_h / SubQ_it);
				//time increased as part of convolution integral
				SubQ_it += 1;
				//Delay_SubQ_TS_sum sums three above terms to construct convolution integral
				Delay_SubQ_TS_sum = first + second + third;
				//SubQi_Routed_m vector at time timeStep_routed is stored
				input->SubQi_Routed_m[timeStep_routed] = input->SubQi_Routed_m[timeStep_routed] + Delay_SubQ_TS_sum;
			}
		}
	}

	//flag = 0, unrouted flow
	else {
		timeStep_routed = timeStep;
		//PAQi_Routed_m increased by PAQ_m; Pervious area runoff, including runon from impervious area (m)
		input->PAQi_Routed_m[timeStep_routed] = input->PAQi_Routed_m[timeStep_routed] + PAQ_m;
		//DCIAQi_Routed_m increased by DCIAQ_m; Directly connected impervious area runoff (m)
		input->DCIAQi_Routed_m[timeStep_routed] = input->DCIAQi_Routed_m[timeStep_routed] + DCIAQ_m;
		//WAQi_Routed_m increased by WAQ_m; Water surface area runoff (m)
		input->WAQi_Routed_m[timeStep_routed] = input->WAQi_Routed_m[timeStep_routed] + WAQ_m;
		//SubQi_Routed_m increased by SubQ_m; Subsurface runoff (m)
		input->SubQi_Routed_m[timeStep_routed] = input->SubQi_Routed_m[timeStep_routed] + SubQ_m;
	}

	double Qcheck_m = 0;
	if (timeStep > 0) {
		//Remove any anomalously low pervious surface runoff generated by Yang routing algorithm 
		if (input->PAQi_Routed_m[timeStep] > 0 && input->PAQi_Routed_m[timeStep] < Constant_1E_negative7) {
			input->PAQi_Routed_m[timeStep] = 0;
		}
		//Set a taper of runoff by dividing it in half
		if (input->PAQi_Routed_m[timeStep] >= Constant_1E_negative7 && input->PAQi_Routed_m[timeStep] < Constant_1E_negative5) {
			Qcheck_m = input->PAQi_Routed_m[timeStep - 1] / 2;
			if (Qcheck_m < Constant_1E_negative7) {
				input->PAQi_Routed_m[timeStep] = 0;
			}
			else {
				input->PAQi_Routed_m[timeStep] = Qcheck_m;
			}
		}
		//Remove any anomalously low impervious surface runoff generated by Yang routing algorithm 
		if (input->DCIAQi_Routed_m[timeStep] > 0 && input->DCIAQi_Routed_m[timeStep] < Constant_1E_negative7) {
			input->DCIAQi_Routed_m[timeStep] = 0;
		}
		//Set a taper of runoff by dividing it in half
		if (input->DCIAQi_Routed_m[timeStep] > Constant_1E_negative7 && input->DCIAQi_Routed_m[timeStep] < Constant_1E_negative5) {
			Qcheck_m = input->DCIAQi_Routed_m[timeStep - 1] / 2;
			if (Qcheck_m < Constant_1E_negative7) {
				input->DCIAQi_Routed_m[timeStep] = 0;
			}
			else {
				input->DCIAQi_Routed_m[timeStep] = Qcheck_m;
			}
		}
		//Remove any anomalously low water surface runoff generated by Yang routing algorithm 
		if (input->WAQi_Routed_m[timeStep] > 0 && input->WAQi_Routed_m[timeStep] < Constant_1E_negative7) {
			input->WAQi_Routed_m[timeStep] = 0;
		}
		//Set a taper of runoff by dividing it in half
		if (input->WAQi_Routed_m[timeStep] >= Constant_1E_negative7 && input->WAQi_Routed_m[timeStep] < Constant_1E_negative5) {
			Qcheck_m = input->WAQi_Routed_m[timeStep - 1] / 2;
			if (Qcheck_m < Constant_1E_negative7) {
				input->WAQi_Routed_m[timeStep] = 0;
			}
			else {
				input->WAQi_Routed_m[timeStep] = Qcheck_m;
			}
		}
	}
	//If timeStep greater than first, then total runoff is sum of all runoff components
	if (timeStep > 0) {
		//Total runoff is sum of subsurface, impervious surface, pervious surface runoff
		input->Qi_Routed_m[timeStep] = input->SubQi_Routed_m[timeStep] + input->DCIAQi_Routed_m[timeStep] + input->PAQi_Routed_m[timeStep] + input->WAQi_Routed_m[timeStep];
	}
	//Else timeStep equals 0 then first timeStep of runoff includes the initial starting baseflow value
	else {
		//Total runoff is sum of initial subsurface runoff, impervious surface, pervious surface runoff
		//Note: Q0_mph variable is from HydroPlusConfig.xml, and with To_m2ph is used to derive AveSMD
		input->Qi_Routed_m[timeStep] = input->InputXml["Discharge_Subsurface_Initial_mph"] * Ratio_Hour_to_Second * input->SimulationTimeStep_Duration_sec[timeStep] + input->DCIAQi_Routed_m[timeStep] + input->PAQi_Routed_m[timeStep];
	}

	//Store catchment runoff depth (m) as routed (e.g.,time distributed) vector values
	//Runoff_Impervious_R_m (m) increased by DCIAQi_Routed_m
	input->RepoDict["Runoff_Impervious_R_m"] = input->RepoDict["Runoff_Impervious_R_m"] + input->DCIAQi_Routed_m[timeStep];
	//Runoff_Pervious_R_m (m) increased by PAQi_Routed_m
	input->RepoDict["Runoff_Pervious_R_m"] = input->RepoDict["Runoff_Pervious_R_m"] + input->PAQi_Routed_m[timeStep];
	//Runoff_Water_R_m (m) increased by PAQi_Routed_m
	input->RepoDict["Runoff_Water_R_m"] = input->RepoDict["Runoff_Water_R_m"] + input->WAQi_Routed_m[timeStep];
	//Runoff_Subsurface_R_m (m) increased by SubQi_Routed_m
	input->RepoDict["Runoff_Subsurface_R_m"] = input->RepoDict["Runoff_Subsurface_R_m"] + input->SubQi_Routed_m[timeStep];
	//Runoff_R_m (m) increased by Qi_Routed_m
	input->RepoDict["Runoff_R_m"] = input->RepoDict["Runoff_R_m"] + input->Qi_Routed_m[timeStep];
}