///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                    Module Library:                    //
//                      ihacres_cal                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//                    ihacres_cal.cpp                    //
//                                                       //
//                 Copyright (C) 2006 by                 //
//                     Stefan Liersch                    //
//-------------------------------------------------------//
//                                                       //
//    e-mail:     stefan.liersch@ufz.de                  //
//                stefan.liersch@gmail.com                   //
//                                                       //
//                     2006-07-27                        //
///////////////////////////////////////////////////////////
//                   
//-------------------------------------------------------//
// References:											 
// Jakeman, A.J. / Hornberger, G.M (1993).				 
//   How Much Complexity Is Warranted in a				 
//	 Rainfall-Runoff Model?								 
//	 Water Resources Research, (29), NO. 8 (2637-2649)   
// Kokkonen, T. S. et al. (2003).
//   Predicting daily flows in ungauged catchments:
//   model regionalization from catchment descriptors
//   at the Coweeta Hydrologic Laboratory, North Carolina
//   Hydrological Processes (17), 2219-2238
// Croke, B. F. W., W. S. Merritt, et al. (2004).
//   A dynamic model for predicting hydrologic response
//   to land cover changes in gauged and
//   ungauged catchments.
//   Journal Of Hydrology 291(1-2): 115-131.
//-------------------------------------------------------//

//---------------------------------------------------------
// ToDo
//---------------------------------------------------------
//	- what happens or what should happen if temperature data
//	  are not available?
//---------------------------------------------------------
#include <iostream> // used for textfile output (test only)
#include <fstream>  // used for textfile output (test only)

#include <stdlib.h> // random numbers
#include <time.h>	// random numbers

#include "ihacres_cal.h"
#include ".\..\model_tools.h"
#include "convert_sl.h"
//---------------------------------------------------------

// Constructor
Cihacres_cal::Cihacres_cal(void)
{
	//-----------------------------------------------------
	// 1. Info...
	Set_Name("IHACRES Calibration");

	Set_Author("copyrights (c) 2006 Stefan Liersch");

	Set_Description("Calibration Tool for the Model IHACRES \n \n \n"
		"Reference: \n \n"
		"Jakeman, A.J. / Hornberger, G.M. (1993). \n"
		"How Much Complexity Is Warranted in a Rainfall-Runoff Model? \n"
		"Water Resources Research, (29), NO. 8 (2637-2649) \n \n"   
		"Kokkonen, T. S. et al. (2003). \n"
		"Predicting daily flows in ungauged catchments:"
		"model regionalization from catchment descriptors"
		"at the Coweeta Hydrologic Laboratory, North Carolina \n "
		"Hydrological Processes (17), 2219-2238 \n \n"
		"Croke, B. F. W., W. S. Merritt, et al. (2004).\n"
		"A dynamic model for predicting hydrologic response"
		"to land cover changes in gauged and"
		"ungauged catchments. \n"
		"Journal Of Hydrology 291(1-2): 115-131."
		);

	//-----------------------------------------------------
	// Create First Module Dialog
	_CreateDialog1();
	//-----------------------------------------------------
}

//---------------------------------------------------------
// Destructor
Cihacres_cal::~Cihacres_cal(void)
{
	if (m_bSnowModule)
	{
		delete m_pSnowModule;
		delete[] m_pMeltRate;
	}
}

//---------------------------------------------------------
// This function is executed when the user is pressing the OK button
// in the first module dialog
bool Cihacres_cal::On_Execute(void)
{
	CSG_Table		*pTable, *pTable_out, *pTable_parms;
	int				dateField;	// table field numbers
	int				dischargeField;
	int				pcpField;
	int				tmpField;

	int				j, k;		//counter
	bool			bTMP;		// true, if temperature data are used to
								// to refine wetness index estimation.
	bool			bnl_only;	// if true, only non-linear module 
								// calibration will be performed

	int				IHAC_version; // Different versions of IHACRES exist, corresponding
								  // to the version...

	int				storconf;	// storage configuration
								// 0 = single,
								// 1 = two parallel,
								// 2 = two in series
	double			dev_eRainDis; // deviation of estimated effective rainfall
								  // to observed discharge in [%]
	double			dev_eRainPCP; // deviation of estimated effective rainfall
								  // to observed precipitation in [%]
	int				first;		// number of the record of first date
	int				last;		// number of the record of last date

	int				nTop = 0;
	//---------------------------------------------------------
	// Assign parameters from First Module Dialog
	pTable			= Parameters("TABLE")			->asTable();
	// Field numbers
	dateField		= Parameters("DATE_Field")		->asInt();
	dischargeField	= Parameters("DISCHARGE_Field")	->asInt();
	pcpField		= Parameters("PCP_Field")		->asInt();
	tmpField		= Parameters("TMP_Field")		->asInt();

	bTMP			= Parameters("USE_TMP")			->asBool();
	m_c_lb			= Parameters("CFAC_lb")			->asDouble();
	m_c_ub			= Parameters("CFAC_ub")			->asDouble();
	m_TwConst_lb	= Parameters("TwFAC_lb")		->asDouble();
	m_TwConst_ub	= Parameters("TwFAC_ub")		->asDouble();
	bnl_only		= Parameters("NL_ONLY")			->asBool();
	m_nsim			= Parameters("NSIM")			->asInt();
	m_area			= Parameters("AREA")			->asDouble();
	dev_eRainDis	= Parameters("DEV_ERAINDIS")	->asDouble();
	dev_eRainPCP	= Parameters("DEV_ERAINPCP")	->asDouble();
	storconf		= Parameters("STORAGE")			->asInt();
	IHAC_version	= Parameters("IHACVERS")		->asInt();
	m_bSnowModule		= Parameters("SNOW_MODULE")		->asBool();
	first = last	= 0;
	//---------------------------------------------------------
	m_counter = 0;
	//---------------------------------------------------------

	// After pressing OK in the first Dialog Dialog2 appears.
	// In the sec
	if (_CreateDialog2(bTMP,bnl_only,storconf, IHAC_version, m_bSnowModule ,pTable,dateField,dischargeField))
	//---------------------------------------------------------
	{
		//---------------------------------------------------------
		// Searching the first and the last record of the time range
		ihac.AssignFirstLastRec(*pTable, first, last, m_date1, m_date2, dateField);
		//---------------------------------------------------------

		//---------------------------------------------------------
		// Initialize arrays
		m_sizeAll = last - first + 1;
		_InitArrays_sizeAll(m_sizeAll,m_bSnowModule);
		_InitArrays_nsim(m_nsim, m_bSnowModule);
		if (m_bSnowModule)
		{
			m_pSnowModule = new CSnowModule(m_sizeAll);
		}
		//---------------------------------------------------------
	
		//---------------------------------------------------------
		// Initialize global variables / set variables to zero
		_ZeroVariables();
		//---------------------------------------------------------

		//---------------------------------------------------------
		// Assign selected time range to vector m_vec_date, 
		// discharge[], pcp[], and tmp[]
		for (j = 0, k = first; j < m_sizeAll, k < last + 1; j++, k++)
		{
			// put the date string into the vector
			m_vec_date[j].append(pTable->Get_Record(k)->asString(dateField));
			m_p_Q_obs_m3s[j]= pTable->Get_Record(k)->asDouble(dischargeField);
			m_pPCP[j]		= pTable->Get_Record(k)->asDouble(pcpField);
			m_pTMP[j]		= pTable->Get_Record(k)->asDouble(tmpField);
		}
		//---------------------------------------------------------

		//---------------------------------------------------------
		// Convert streamflow from m3/s to mm/day
		m_p_Q_obs_mmday = model_tools::m3s_to_mmday(m_p_Q_obs_m3s,m_p_Q_obs_mmday,m_sizeAll,m_area);
		// calculate sum of observed discharge in [mm]
		m_sum_obsDisMM = ihac.SumVector(m_p_Q_obs_mmday, m_sizeAll);

		//---------------------------------------------------------
		// Perform non-linear module calibration
		// NOTE: bool m_bSnowModule will always be set to 'false' here,
		// because the snow module will be calibrated later on!
		//---------------------------------------------------------
		_CalcNonLinearModule(bTMP, IHAC_version, false, m_nsim, dev_eRainDis, dev_eRainPCP);
		// store information about valid simulations in 'm_p_bValidSim'
		//---------------------------------------------------------

		if (!bnl_only)
		{
			//---------------------------------------------------------
			// Create Output Table for simulation parameters
			pTable_parms = SG_Create_Table();
			if (storconf == 0)
				_CreateTable1Storage(pTable_parms, m_bSnowModule);
			if (storconf == 1)
				_CreateTable2Parallel(pTable_parms, m_bSnowModule);
			//---------------------------------------------------------

			_Get_Best_NonLinear_Sims(5, nTop, IHAC_version);
			//---------------------------------------------------------
			// Perform linear module calculations
			_CalcLinearModule(bTMP, IHAC_version, storconf, m_bSnowModule, nTop, *pTable_parms);
			//---------------------------------------------------------

		} // end if (!bnl_only)

		//---------------------------------------------------------
		// Create output table and add it to SAGA workspace
		pTable_out = SG_Create_Table();

		_CreateTableCal(pTable_out, IHAC_version, m_bSnowModule, m_pNLP_random, m_nsim);

		// Add table to SAGA Workspace
		pTable_out->Set_Name("IHACRES_cal");
		Parameters("TABLEout")->Set_Value(pTable_out);

		//---------------------------------------------------------

		/////////////////////////////////////////
		// Delete Arrays and resize (zero)vectors
		//---------------------------------------
		if (nTop > 0)
		{
			delete m_pNLP_top;
			delete m_pNLP_random;
		}
		if (m_bSnowModule)
		{
			delete[] m_pT_Rain_rand;
			delete[] m_pT_Melt_rand;
			delete[] m_pDD_FAC_rand;
		}
		//delete[] m_p_Q_obs_m3s;
		//delete[] m_pPCP;
		//delete[] m_pTMP;
		//delete[] m_pTw;
		//delete[] m_pWI;
		//delete[] m_pExcessRain;
		//delete[] m_pQ_sim_m3s;
		//delete[] m_p_Q_obs_mmday;
		//delete[] m_pSum_eRainMM;
		//delete[] m_p_bValidSim;
		m_vec_date.resize(0);
		/////////////////////////////////////////

	} // end if _CreateDialog2
	//---------------------------------------------------------

	return( true );   // bool Cihacres_cal::On_Execute(void)
} // end function On_Execute()



/////////////////////////////////////////////////////////////////////////////
//
//                         PRIVATE MEMBER FUNCTIONS
//
/////////////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------
void Cihacres_cal::_ZeroVariables()
{
	m_sum_obsDisMM		= 0.0;
	m_sum_eRainGTpcp	= 0.0;
}
//---------------------------------------------------------------------
void Cihacres_cal::_InitArrays_sizeAll(int size, bool bSnowModule)
{
	m_vec_date.resize(size);
	m_p_Q_obs_m3s	= new double[size];
	m_p_Q_obs_mmday = new double[size];
	m_pPCP			= new double[size];
	m_pTMP			= new double[size];
	m_pTw			= new double[size];
	m_pWI			= new double[size];
	m_pExcessRain	= new double[size];
	m_pQ_sim_m3s	= new double[size];
	m_pQ_sim_mmday	= new double[size];
	if (bSnowModule) m_pMeltRate = new double[size];
}
//---------------------------------------------------------------------
void Cihacres_cal::_InitArrays_nsim(int nsim, bool bSnowModule)
{
	// random nsim arrays
	m_pNLP_random	= new C_IHAC_NonLinearParms(nsim);
	m_pSum_eRainMM	= new double[nsim];
	m_p_bValidSim	= new bool[nsim];

	if (bSnowModule) // if snow module is active
	{
		m_pT_Rain_rand = new double[nsim];
		m_pT_Melt_rand = new double[nsim];
		m_pDD_FAC_rand = new double[nsim];
	}
	for (int i = 0; i < nsim; i++)
		m_p_bValidSim[i] = false;
}
//---------------------------------------------------------------------


//---------------------------------------------------------------------
//
//						Calculations of the Non-linear Module
//						Calculate ExcessRain time series
//
//---------------------------------------------------------------------
void Cihacres_cal::_CalcNonLinearModule(bool bTMP, int IHAC_version, bool bSnowModule, int nsim, 
										double dev_eRainDis, double dev_eRainPCP)
{
	// initialize random function
	srand((unsigned) time(NULL)); // using time.h
	for (int j = 0; j < nsim && Set_Progress(j,nsim); j++)
	{
		m_pNLP_random->mp_tw[j] = _Random_double(m_TwConst_lb, m_TwConst_ub);
		m_pNLP_random->mp_f[j]  = _Random_double(m_f_lb, m_f_ub);
		m_pNLP_random->mp_c[j]  = _Random_double(m_c_lb, m_c_ub);

		_RunNonLinearModule(bTMP,IHAC_version,m_pTMP,m_pPCP,m_sizeAll,m_pNLP_random,j,m_pExcessRain, bSnowModule);

		// calculate difference between discharge and excess rainfall
		m_pNLP_random->mp_eR_flow_dif[j] = m_sum_obsDisMM - m_pSum_eRainMM[j];
		if (m_pNLP_random->mp_eR_flow_dif[j] < 0.0) m_pNLP_random->mp_eR_flow_dif[j] *= -1;

		//---------------------------------------------------------
		// decide whether a simulation is valid or not
		//---------------------------------------------------------
		if ( m_sum_eRainGTpcp  / m_pSum_eRainMM[j] * 100 <= dev_eRainPCP &&
			 m_pNLP_random->mp_eR_flow_dif[j] / m_sum_obsDisMM   * 100 <= dev_eRainDis   )
		{
			m_p_bValidSim[j] = true;
		}
	} // end nsim calibration (non-linear module)
}
//---------------------------------------------------------------------

void Cihacres_cal::_RunNonLinearModule(bool bTMP, int IHAC_version,
									   double *p_TMP, double *p_PCP, int size,
									   C_IHAC_NonLinearParms *nl_parms, int j,
									   double *p_ExcessRain, bool dummy_bSnowModule)
{
	double eR_init = 0.0;
	if (bTMP)
	{
		switch(IHAC_version)
		{
		case 0: // Jakeman & Hornberger (1993)
			ihac.CalcWetnessTimeConst(p_TMP, m_pTw, nl_parms->mp_tw[j], nl_parms->mp_f[j],m_sizeAll);
			break;
		case 1: // Croke et al. (2005) Redesign
			ihac.CalcWetnessTimeConst_Redesign(p_TMP, m_pTw, nl_parms->mp_tw[j], nl_parms->mp_f[j], m_sizeAll);
			break;
		}
	}

	// calculate the Catchment Wetness Index
	switch(IHAC_version)
	{
	case 0: ihac.CalcWetnessIndex(m_pTw,p_PCP,p_TMP,m_pWI,0.5, nl_parms->mp_c[j],dummy_bSnowModule,-100.0,m_sizeAll);
		break;
	case 1: ihac.CalcWetnessIndex_Redesign(m_pTw, p_PCP, m_pWI, 0.5,dummy_bSnowModule,-100.0, m_sizeAll);
		break;
	}

	// calculate Effective Rainfall
	switch(IHAC_version)
	{
	case 0: 
		m_pSum_eRainMM[j] = ihac.CalcExcessRain(p_PCP, p_TMP, m_pWI,p_ExcessRain, eR_init, m_sum_eRainGTpcp,m_sizeAll, false, 0,0,0);
		break;
	case 1: 
		m_pSum_eRainMM[j] = ihac.CalcExcessRain_Redesign(p_PCP, p_TMP, m_pWI, p_ExcessRain, eR_init, m_sum_eRainGTpcp,m_sizeAll,
							nl_parms->mp_c[j], nl_parms->mp_l[j], nl_parms->mp_p[j], false, 0,0,0);
		break;
	}
}

//---------------------------------------------------------------------
//				Find best non-linear module simulations
//---------------------------------------------------------------------
void Cihacres_cal::_Get_Best_NonLinear_Sims(int nTopSims, int &nTop, int IHAC_version)
{
	//---------------------------------------------------------
	// local variables
	//---------------------------------------------------------
	//nTopSims: number of top simulations from non-linear module
	//			used to calibrate the linear module
	//int nTop: number of Top simulations if the number
	//			of valid simulations is less than nTopSims
	int j,k;
	int n_valid = _CountNumberTrue(m_p_bValidSim,m_nsim); //count number of valid simulation

	double*	Tw_temp;
	double*	c_temp;
	double*	f_temp;
	double*	e_temp; // temporary array dif_eRainPCP
	double*	l_temp;
	double* p_temp;
	//---------------------------------------------------------

	//---------------------------------------------------------
	// TOP 5 non-linear module simulations for use
	// in the linear module
	//---------------------------------------------------------
	if (n_valid > 0)
	{
		// create temp arrays (c, Tw, f and dif_eRainDis) size = nvalid
		c_temp	= new double[n_valid];
		Tw_temp	= new double[n_valid];
		f_temp	= new double[n_valid];
		e_temp	= new double[n_valid]; // temporary array dif_eRainPCP
		if (IHAC_version == 1) // Croke et al. (2005) Redesign
		{
			l_temp = new double[n_valid];
			p_temp = new double[n_valid];
		}

		k = 0;
		for (j = 0; j < m_nsim; j++)
		{
			if (m_p_bValidSim[j])
			{
				c_temp[k]  = m_pNLP_random->mp_c[j];
				Tw_temp[k] = m_pNLP_random->mp_tw[j];
				f_temp[k]  = m_pNLP_random->mp_f[j];
				e_temp[k]  = m_pNLP_random->mp_eR_flow_dif[j];
				if (IHAC_version == 1)
				{
					l_temp[k] = m_pNLP_random->mp_l[j];
					p_temp[k] = m_pNLP_random->mp_p[j];
				}
				k++;
			}
		}

		// In case the number of valid simulations is less then 5
		// the number of Top-Simulations must be reduced to
		// the number of valid simulations.
		if (n_valid >= nTopSims) nTop = nTopSims;
		else nTop = n_valid;

		m_pNLP_top = new C_IHAC_NonLinearParms(nTop);

		int		*Top5Index	= new int[nTop];

		// Identify TOP 5 simulations
		// Assign TOP 5 parameter settings to temp1 arrays
		_FindTopIndex(e_temp, n_valid, Top5Index, nTop);
				
		for (j = 0; j < nTop; j++)
		{
			m_pNLP_top->mp_tw[j]	= Tw_temp[Top5Index[j]];
			m_pNLP_top->mp_c[j]	= c_temp[Top5Index[j]];
			m_pNLP_top->mp_f[j]	= f_temp[Top5Index[j]];
			m_pNLP_top->mp_eR_flow_dif[j] = e_temp[Top5Index[j]];
			if (IHAC_version == 1) // Croke et al. (2005) Redesign
			{
				m_pNLP_top->mp_l[j] = l_temp[Top5Index[j]];
				m_pNLP_top->mp_p[j] = p_temp[Top5Index[j]];
			}
		}

		// delete pointers
		delete[] c_temp;
		delete[] Tw_temp;
		delete[] f_temp;
		delete[] e_temp;
		if (IHAC_version == 1) // Croke et al. (2005) Redesign
		{
			delete[] l_temp;
			delete[] p_temp;
		}
	}
}

//---------------------------------------------------------------------
//				Calculations of the Linear Module
//---------------------------------------------------------------------
void Cihacres_cal::_CalcLinearModule(bool bTMP, int IHAC_version, int storconf, bool bSnowModule,
									 int nTop, CSG_Table &pTable_parms)
{
	//---------------------------------------------------------
	// local variables
	//---------------------------------------------------------
	double	eR_init			= 0.0;
	int		j,k;		// counter
	double	NSE				= 0.0; // Nash-Sutcliffe Efficiency
	double	NSE_highflow	= 0.0; // Nash-Sutcliffe Efficiency adapted to high flow
	double	NSE_lowflow		= 0.0; // Nash-Sutcliffe Efficiency adapted to low flow
	double	NSE_temp		= 0.0;
	// efficiency output parameters to status bar
	double	NSE_max			= -100.0;
	string	nse;
	string	nse_text;
	//---------------------------------------------------------

	// simulations with lowest difference between excess rain and streamflow volumes
	// initialize random function
	srand((unsigned) time(NULL)); // using time.h

	switch (storconf)
	{
	case 0: // single storage
		m_pLP_random = new C_IHAC_LinearParms(m_nsim, 1);

		for (j = 0; j < nTop; j++) // Top5 Simulations
		{
			//---------------------------------------------------------
			// non-linear module calculations
			//---------------------------------------------------------
			_RunNonLinearModule(bTMP,IHAC_version,m_pTMP,m_pPCP,m_sizeAll,m_pNLP_top,j,m_pExcessRain, false);
			//---------------------------------------------------------

			//---------------------------------------------------------
			// calculate linear module
			for (k = 0; k < m_nsim && Set_Progress(k,m_nsim); k++)
			{
				m_pLP_random->a[k] = _Random_double(m_a_lb, m_a_ub);
				m_pLP_random->b[k] = 1 + m_pLP_random->a[k];

				// calculate simulated streamflow
				ihac.SimStreamflowSingle(m_pExcessRain, m_p_Q_obs_mmday[0],
					m_pQ_sim_mmday, m_delay, m_pLP_random->a[k], m_pLP_random->b[k],m_sizeAll);
				//---------------------------------------------------------

				//---------------------------------------------------------
				// calculate Nash-Sutcliffe efficiency
				NSE = model_tools::CalcEfficiency(m_p_Q_obs_mmday, m_pQ_sim_mmday,m_sizeAll);
				NSE_highflow = model_tools::Calc_NSE_HighFlow(m_p_Q_obs_mmday, m_pQ_sim_mmday, m_sizeAll);
				NSE_lowflow  = model_tools::Calc_NSE_LowFlow(m_p_Q_obs_mmday, m_pQ_sim_mmday, m_sizeAll);
				//---------------------------------------------------------

				//---------------------------------------------------------
				// Write all simulations with required goodness in a table
				NSE_temp = ihac._Assign_NSE_temp(m_obj_func, NSE, NSE_highflow, NSE_lowflow);

				if (NSE_temp >= m_NSEmin)
				{
					if (NSE_temp > NSE_max)
					{
						NSE_max = NSE_temp;
						nse = convert_sl::Double2String(NSE_max);
						nse_text = "max. Nash Sutcliffe Efficiency: " + nse;
						Process_Set_Text(nse_text.c_str());
					}

					_WriteTable1Storage(&pTable_parms, bSnowModule, NSE,m_pNLP_top->mp_eR_flow_dif[j],
						m_pNLP_top->mp_tw[j],m_pNLP_top->mp_f[j],m_pNLP_top->mp_c[j],
						0,0,0,
						m_pLP_random->a[k], m_pLP_random->b[k]);
				}
				//---------------------------------------------------------

			} // end calibration (linear module)
			//---------------------------------------------------------
		} // end TOP 5 Simulations

		delete m_pLP_random;

		break;

	case 1: // two parallel storages
		m_pLP_random = new C_IHAC_LinearParms(m_nsim, 2);
		m_pLP_random->aq = new double[m_nsim];
		m_pLP_random->as = new double[m_nsim];
		m_pLP_random->bq = new double[m_nsim];
		m_pLP_random->bs = new double[m_nsim];

		for (j = 0; j < nTop; j++)
		{
			//---------------------------------------------------------
			// non-linear module calculations
			//---------------------------------------------------------
			_RunNonLinearModule(bTMP,IHAC_version,m_pTMP,m_pPCP,m_sizeAll,m_pNLP_top,j,m_pExcessRain, false);
			//---------------------------------------------------------

			//---------------------------------------------------------
			// calculate linear module
			double check_vq; // check whether parameters a and b are appropriate
			for (k = 0; k < m_nsim && Set_Progress(k,m_nsim); k++)
			{
				// calculate random a and b parameters as long as m_vq is valid
				do
				{
					m_pLP_random->aq[k] = _Random_double(m_aq_lb, m_aq_ub);
					m_pLP_random->as[k] = _Random_double(m_as_lb, m_as_ub);
					m_pLP_random->bq[k] = _Random_double(m_bq_lb, m_bq_ub);
					// Calculate parameter m_vq to check parms a and b
					// Equation after Jakeman & Hornberger (1993)
					check_vq = m_pLP_random->bq[k] / (1 + m_pLP_random->aq[k]);
				}
				while (check_vq < 0.0 || check_vq > 1.0);
				m_pLP_random->bs[k] = ihac.Calc_Parm_BS(m_pLP_random->aq[k], m_pLP_random->as[k], m_pLP_random->bq[k]);

				//---------------------------------------------------------
				// Simulate streamflow with two parallel storages								
				ihac.SimStreamflow2Parallel(m_pExcessRain, m_pQ_sim_mmday,m_p_Q_obs_mmday[0],
					m_pLP_random->aq[k],m_pLP_random->as[k], m_pLP_random->bq[k], m_pLP_random->bs[k],
					m_vq, m_vs,
					IHAC_version, m_sizeAll, m_delay);	
				//---------------------------------------------------------

				//---------------------------------------------------------
				// Calculate Nash-Sutcliffe efficiency
				NSE = model_tools::CalcEfficiency(m_p_Q_obs_mmday, m_pQ_sim_mmday,m_sizeAll);
				NSE_highflow = model_tools::Calc_NSE_HighFlow(m_p_Q_obs_mmday, m_pQ_sim_mmday, m_sizeAll);
				NSE_lowflow  = model_tools::Calc_NSE_LowFlow(m_p_Q_obs_mmday, m_pQ_sim_mmday, m_sizeAll);
				//---------------------------------------------------------

				//---------------------------------------------------------
				// Write all simulations with required goodness of fit in a table
				NSE_temp = ihac._Assign_NSE_temp(m_obj_func, NSE, NSE_highflow, NSE_lowflow);

				if (NSE_temp >= m_NSEmin)
				{
					if (NSE_temp > NSE_max)
					{
						if (bSnowModule)
						{
							m_snowcal_parms.aq = m_pLP_random->aq[k];
							m_snowcal_parms.as	= m_pLP_random->as[k];
							m_snowcal_parms.bq	= m_pLP_random->bq[k];
							m_snowcal_parms.bs	= m_pLP_random->bs[k];
							m_snowcal_parms.nl_index = j;
							m_snowcal_parms.tw = m_pNLP_top->mp_tw[j];
							m_snowcal_parms.f	 = m_pNLP_top->mp_f[j];
							m_snowcal_parms.c  = m_pNLP_top->mp_c[j];
						}

						NSE_max = NSE_temp;
						nse = convert_sl::Double2String(NSE_max);
						nse_text = "max. Nash Sutcliffe Efficiency: " + nse;
						Process_Set_Text(nse_text.c_str());
					}
					if (!bSnowModule)
					{
						_WriteTable2Parallel(&pTable_parms, bSnowModule,
							NSE, NSE_highflow,NSE_lowflow,m_pNLP_top->mp_eR_flow_dif[j],
							m_vq, m_vs, m_pNLP_top->mp_tw[j],m_pNLP_top->mp_f[j],m_pNLP_top->mp_c[j],
							0,0,0, // snow module parameters
							m_pLP_random->aq[k],m_pLP_random->as[k],m_pLP_random->bq[k],m_pLP_random->bs[k]);
					}
				}
				//---------------------------------------------------------
			} // End calibration (linear module)
		} // End TOP 5

		//----------------------------------------------------------------
		// Snow Module CALIBRATION
		//----------------------------------------------------------------
		if (bSnowModule)
		{
			Process_Set_Text("Calibrating Snow Module...");
			m_pMeltRate = new double[m_sizeAll];
			double *pwi = new double[m_sizeAll];

			// calculate wetness time constant
			if (IHAC_version == 1) // Croke et al. (2005) Redesign
			{
				ihac.CalcWetnessTimeConst_Redesign(m_pTMP, m_pTw, m_snowcal_parms.tw, m_snowcal_parms.f, m_sizeAll);
			} else {
				ihac.CalcWetnessTimeConst(m_pTMP, m_pTw, m_snowcal_parms.tw, m_snowcal_parms.f, m_sizeAll);
			}

			//// calculate pwi
			//if (IHAC_version == 1) // Croke et al. (2005) Redesign
			//{
			//	ihac.CalcWetnessIndex_Redesign(m_pTw, m_pPCP, pwi, m_bSnowModule,m_sizeAll);
			//} else {
			//	ihac.CalcWetnessIndex(m_pTw,m_pPCP,m_pTMP,pwi,m_snowcal_parms.c,m_bSnowModule,m_sizeAll);
			//}

			// simulate
			// initialize random function
			srand((unsigned) time(NULL)); // using time.h

			for (int sim = 0; sim < m_nsim && Set_Progress(sim,m_nsim); sim++)
			{
				// calculate random snow parameters
				m_pSnowModule->Set_T_Rain(_Random_double(m_T_Rain_lb, m_T_Rain_ub));
				m_pSnowModule->Set_T_Melt(_Random_double(m_T_Melt_lb, m_T_Melt_ub));
				m_pSnowModule->Set_DD_FAC(_Random_double(m_DD_FAC_lb, m_DD_FAC_ub));
				m_pT_Rain_rand[sim] = m_pSnowModule->Get_T_Rain();
				m_pT_Melt_rand[sim] = m_pSnowModule->Get_T_Melt();
				m_pDD_FAC_rand[sim] = m_pSnowModule->Get_DD_FAC();

				// calculate pwi
				if (IHAC_version == 1) // Croke et al. (2005) Redesign
				{
					ihac.CalcWetnessIndex_Redesign(m_pTw, m_pPCP, pwi, 0.5,m_bSnowModule, m_pSnowModule->Get_T_Rain(),m_sizeAll);
				} else {
					ihac.CalcWetnessIndex(m_pTw,m_pPCP,m_pTMP,pwi, 0.5, m_snowcal_parms.c,m_bSnowModule, m_pSnowModule->Get_T_Rain(),m_sizeAll);
				}

				m_pSnowModule->Calc_SnowModule(m_pTMP,m_pPCP,m_sizeAll,
											   m_pSnowModule->Get_T_Rain(),
											   m_pSnowModule->Get_T_Melt(),
											   m_pSnowModule->Get_DD_FAC());

				m_pMeltRate = m_pSnowModule->Get_MeltRate(m_pMeltRate,m_sizeAll);

				// calculate excess rainfall
				if (IHAC_version == 1) // Croke et al. (2005) Redesign
				{
					ihac.CalcExcessRain_Redesign(m_pPCP, m_pTMP, pwi, m_pExcessRain, eR_init, m_sum_eRainGTpcp,m_sizeAll,
						m_snowcal_parms.c, m_pNLP_random->mp_l[sim],m_pNLP_random->mp_p[sim],
						//m_snowcal_parms.c, 0.0, 1.4,
						true, m_pSnowModule->Get_T_Rain(),m_pSnowModule->Get_T_Melt(),m_pMeltRate);
					//m_nonlinear.c[j], 0.0, 1.4, bSnowModule, 0,0,0);
				} else {
					ihac.CalcExcessRain(m_pPCP, m_pTMP, pwi,m_pExcessRain,eR_init, m_sum_eRainGTpcp,m_sizeAll,
						true, m_pSnowModule->Get_T_Rain(),m_pSnowModule->Get_T_Melt(),m_pMeltRate);
						//bSnowModule, 0,0,0);
				}

				//---------------------------------------------------------
				// Simulate streamflow with two parallel storages								
				ihac.SimStreamflow2Parallel(m_pExcessRain, m_pQ_sim_mmday,m_p_Q_obs_mmday[0],
					m_snowcal_parms.aq,m_snowcal_parms.as, m_snowcal_parms.bq, m_snowcal_parms.bs,
					m_vq, m_vs,
					IHAC_version, m_sizeAll, m_delay);	
				//---------------------------------------------------------

				//---------------------------------------------------------
				// Calculate Nash-Sutcliffe efficiency
				NSE = model_tools::CalcEfficiency(m_p_Q_obs_mmday, m_pQ_sim_mmday,m_sizeAll);
				NSE_highflow = model_tools::Calc_NSE_HighFlow(m_p_Q_obs_mmday, m_pQ_sim_mmday, m_sizeAll);
				NSE_lowflow  = model_tools::Calc_NSE_LowFlow(m_p_Q_obs_mmday, m_pQ_sim_mmday, m_sizeAll);
				//---------------------------------------------------------
				// Write all simulations with required goodness of fit in a table
				NSE_temp = ihac._Assign_NSE_temp(m_obj_func, NSE, NSE_highflow, NSE_lowflow);

				if (NSE_temp >= m_NSEmin)
				{
					if (NSE_temp > NSE_max)
					{								
						NSE_max = NSE_temp;
						nse = convert_sl::Double2String(NSE_max);
						nse_text = "Calibrating Snow Module... NSE: " + nse;
						Process_Set_Text(nse_text.c_str());
					}

					_WriteTable2Parallel(&pTable_parms, bSnowModule,
						NSE, NSE_highflow,NSE_lowflow, m_pNLP_top->mp_eR_flow_dif[m_snowcal_parms.nl_index],
						m_vq, m_vs, m_snowcal_parms.tw,m_snowcal_parms.f,m_snowcal_parms.c,
						m_pSnowModule->Get_T_Rain(),m_pSnowModule->Get_T_Melt(),m_pSnowModule->Get_DD_FAC(),
						m_snowcal_parms.aq,m_snowcal_parms.as,m_snowcal_parms.bq,m_snowcal_parms.bs);
				}
			}
			delete[] pwi;
			//----------------------------------------------------------------
		} // End snow module calibration

		delete m_pLP_random;

		break;
	case 2: // two storages in series
		break;
	} // end switch

	// Show table in SAGA Workspace
	if (storconf == 0)
		pTable_parms.Set_Name("SimParms1Stor");
	if (storconf == 1)
		pTable_parms.Set_Name("SimParms2Stors");
	Parameters("TABLEparms")->Set_Value(&pTable_parms);

}



/////////////////////////////////////////////////////////////////////////////
//
//							CREATE TABLES
//
/////////////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------
//              Create Table for non linear module output
//---------------------------------------------------------------------
void Cihacres_cal::_CreateTableCal(CSG_Table *pTable_out, int IHAC_version, bool bSnowModule, C_IHAC_NonLinearParms *nl_parms, int size)
{
	int i = 0; // used in function Get_Record(i)
	int x = 0; // field position
	CSG_Table_Record	*pRecord;

	// creating the column titles
	pTable_out->Add_Field("Tw",				 TABLE_FIELDTYPE_Double);
	pTable_out->Add_Field("f",				 TABLE_FIELDTYPE_Double);
	pTable_out->Add_Field("c",				 TABLE_FIELDTYPE_Double);
	if (IHAC_version == 1)
	{
		pTable_out->Add_Field("l",			TABLE_FIELDTYPE_Double);
		pTable_out->Add_Field("p",			TABLE_FIELDTYPE_Double);
	}
	if (bSnowModule)
	{
		pTable_out->Add_Field("T_Rain",		TABLE_FIELDTYPE_Double);
		pTable_out->Add_Field("T_Melt",		TABLE_FIELDTYPE_Double);
		pTable_out->Add_Field("DD_FAC",		TABLE_FIELDTYPE_Double);
	}
	pTable_out->Add_Field("DifExcessFlow",	 TABLE_FIELDTYPE_Double);

	for (int j = 0; j < size; j++)
	{
		x = 0;
		if (m_p_bValidSim[j])
		{
			// add a new record to the table
			pTable_out->Add_Record();
			pRecord	= pTable_out->Get_Record(i);

			// writing the data to the current row
			pRecord->Set_Value(0,nl_parms->mp_tw[j]);
			pRecord->Set_Value(1,nl_parms->mp_f[j]);
			pRecord->Set_Value(2,nl_parms->mp_c[j]);
			if (IHAC_version == 1)
			{
				pRecord->Set_Value(3,nl_parms->mp_l[j]);
				pRecord->Set_Value(4,nl_parms->mp_p[j]);
				x += 2;
			}
			if (bSnowModule)
			{
				pRecord->Set_Value(x+3, m_pT_Rain_rand[j]);
				pRecord->Set_Value(x+4, m_pT_Melt_rand[j]);
				pRecord->Set_Value(x+5, m_pDD_FAC_rand[j]);
				x += 3;
			}

			pRecord->Set_Value(x+3,nl_parms->mp_eR_flow_dif[j]);
			i++;
		}
	}
}

//---------------------------------------------------------------------
//      Create Table (Simulation Parameters, 2 parallel storages)
//---------------------------------------------------------------------
void Cihacres_cal::_CreateTable2Parallel(CSG_Table *pTable_parms, bool bSnowModule)
{
	// creating the column titles
	pTable_parms->Add_Field("NSE",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("NSE_high",			 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("NSE_low",			 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("eR_Dif",			 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("v(q)",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("v(s)",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("Tw",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("f",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("c",				 TABLE_FIELDTYPE_Double);
	if (bSnowModule)
	{
		pTable_parms->Add_Field("T_Rain",		 TABLE_FIELDTYPE_Double);
		pTable_parms->Add_Field("T_Melt",		 TABLE_FIELDTYPE_Double);
		pTable_parms->Add_Field("DD_Fac",		 TABLE_FIELDTYPE_Double);
	}
	pTable_parms->Add_Field("a(q)",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("a(s)",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("b(q)",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("b(s)",				 TABLE_FIELDTYPE_Double);
}

//---------------------------------------------------------------------
//      Create Table (Simulation Parameters, 1 storage)
//---------------------------------------------------------------------
void Cihacres_cal::_CreateTable1Storage(CSG_Table *pTable_parms, bool bSnowModule)
{
	// creating the column titles
	pTable_parms->Add_Field("NSE",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("eR_Dif",			 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("Tw",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("f",				 TABLE_FIELDTYPE_Double);
	if (bSnowModule)
	{
		pTable_parms->Add_Field("T_Rain",			 TABLE_FIELDTYPE_Double);
		pTable_parms->Add_Field("T_Melt",			 TABLE_FIELDTYPE_Double);
		pTable_parms->Add_Field("DD_Fac",			 TABLE_FIELDTYPE_Double);
	}
	pTable_parms->Add_Field("c",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("a",				 TABLE_FIELDTYPE_Double);
	pTable_parms->Add_Field("b",				 TABLE_FIELDTYPE_Double);
}

//---------------------------------------------------------------------
//    Write data into table Parms2Parallel (Simulation Parameters)
//---------------------------------------------------------------------
void Cihacres_cal::_WriteTable2Parallel(CSG_Table *pTable_parms, bool bSnowModule,
										double NSE, double NSE_high, double NSE_low,
										double eR, double vq, double vs,
										double Tw, double f, double c,
									    double T_Rain, double T_Melt, double DD_Fac,
									    double aq, double as, double bq, double bs
										)
{
	//static int		j = 0;
	int count = 1;
	CSG_Table_Record	*pRecord;

	// add a new record to the table
	pTable_parms->Add_Record();
	pRecord	= pTable_parms->Get_Record(m_counter);

	// writing the data to the current row
	pRecord->Set_Value(0,NSE);
	pRecord->Set_Value(1,NSE_high);
	pRecord->Set_Value(2,NSE_low);
	pRecord->Set_Value(3,eR);
	pRecord->Set_Value(4,vq);
	pRecord->Set_Value(5,vs);
	pRecord->Set_Value(6,Tw);
	pRecord->Set_Value(7,f);
	pRecord->Set_Value(8,c);
	if (bSnowModule)
	{
		pRecord->Set_Value(8+count,T_Rain);
		count++;
		pRecord->Set_Value(8+count,T_Melt);
		count++;
		pRecord->Set_Value(8+count,DD_Fac);
		count++;
	}
	pRecord->Set_Value(8+count,aq);
	count++;
	pRecord->Set_Value(8+count,as);
	count++;
	pRecord->Set_Value(8+count,bq);
	count++;
	pRecord->Set_Value(8+count,bs);
	count++;
	m_counter++;
	//j += 1;
}

//---------------------------------------------------------------------
//    Write data into table Parms1Storage (Simulation Parameters)
//---------------------------------------------------------------------
void Cihacres_cal::_WriteTable1Storage(CSG_Table *pTable_parms, bool bSnowModule, double NSE, double eR,
									   double Tw, double f, double c,
									   double T_Rain, double T_Melt, double DD_Fac,
									   double a, double b)
{
	int count = 1;
	//static int		j = 0;
	CSG_Table_Record	*pRecord;

	// add a new record to the table
	pTable_parms->Add_Record();
	pRecord	= pTable_parms->Get_Record(m_counter);

	// writing the data to the current row
	pRecord->Set_Value(0,NSE);
	pRecord->Set_Value(1,eR);
	pRecord->Set_Value(2,Tw);
	pRecord->Set_Value(3,f);
	pRecord->Set_Value(4,c);
	if (bSnowModule)
	{
		pRecord->Set_Value(4+count,T_Rain);
		count++;
		pRecord->Set_Value(4+count,T_Melt);
		count++;
		pRecord->Set_Value(4+count,DD_Fac);
		count++;
	}
	pRecord->Set_Value(4+count,a);
	count++;
	pRecord->Set_Value(4+count,b);
	count++;
	m_counter++;
	//j += 1;
}

/////////////////////////////////////////////////////////////////////////////
//
//							HELPER
//
/////////////////////////////////////////////////////////////////////////////

//---------------------------------------------------------------------
//                       Calculate random number
//---------------------------------------------------------------------
double Cihacres_cal::_Random_double(double lb, double ub)
{
	// lb = lower bound, ub = upper bound
	double random;

	random = rand()/(1.+RAND_MAX);
	random = lb + (ub - lb) * random;

	return(random);
}

//---------------------------------------------------------------------
// Count number of true values in bool array
//---------------------------------------------------------------------
int Cihacres_cal::_CountNumberTrue(bool *array, int size)
{
	int counter = 0;
	// count number of true values
	for (int i = 0; i < size; i++)
	{
		if (array[i]) counter += 1;
	}
	return(counter);
}

//---------------------------------------------------------------------
// Find Top5 Simulation indices
//---------------------------------------------------------------------
void Cihacres_cal::_FindTopIndex(double *array, int size, int *TopIndex, int top)
{
	double min_temp = 99999999.0;
	double x		= -99999999.0;
	int index		= 0;

	for (int k = 0; k < top; k++)
	{
		for (int j = 0; j < size; j++)
		{
			if (array[j] < min_temp && array[j] > x)
			{
				min_temp = array[j];
				index = j;
			}
		}
		TopIndex[k] = index;
		x = min_temp;
		min_temp = 99999999.0;
	}
}

//---------------------------------------------------------------------
// Find Top x simulations depending on Nash-Sutcliffe Efficiency
//---------------------------------------------------------------------
void Cihacres_cal::_FindTopNSE(double *array, int size, int *TopNSE, int top)
{
	double max_temp = -10.0;
	double x		= 1.0;
	int index		= 0;

	for (int k = 0; k < top; k++)
	{
		for (int j = 0; j < size; j++)
		{
			if (array[j] > max_temp && array[j] < x)
			{
				max_temp = array[j];
				index = j;
			}
		}
		TopNSE[k] = index;
		x = max_temp;
		max_temp = -10.0;
	}
}




/////////////////////////////////////////////////////////////////////////////
//
//									DIALOGS
//
/////////////////////////////////////////////////////////////////////////////

void Cihacres_cal::_CreateDialog1()
{
	CSG_Parameter	*pNode;
	CSG_String s;

	pNode = Parameters.Add_Table(
		NULL	, "TABLE"	, _TL("Table"),
		"",
		PARAMETER_INPUT
	);

	Parameters.Add_Table_Field(
		pNode	, "DATE_Field"	, _TL("Date Column"),
		"Select the column containing the Date"
	);

	Parameters.Add_Table_Field(
		pNode	, "DISCHARGE_Field"	, _TL("Streamflow Column"),
		"Select the Column containing Discharge Values"
	);

	Parameters.Add_Table_Field(
		pNode	, "PCP_Field"	, _TL("Precipitation Column"),
		"Select the Column containing precipitation Values"
	);

	Parameters.Add_Table_Field(
		pNode	, "TMP_Field"	, _TL("Temperature Column"),
		"Select the Column containing Temperature Values"
	);

	Parameters.Add_Value(
		pNode,	"USE_TMP",	_TL("Using temperature data?"),
		_TL("If checked, then temperature data are used."),
		PARAMETER_TYPE_Bool, true
	);

	pNode = Parameters.Add_Value(
		NULL,	"NSIM",	_TL("Number of Simulations"),
		_TL("Number of Simulations for Calibration"),
		PARAMETER_TYPE_Int,
		1000, 1, true, 10000000, true
	);

	pNode = Parameters.Add_Value(
		NULL,	"AREA",	_TL("Area of the Watershed in [km2]"),
		_TL("Area of the Watershed in [km2] used for unit conversion"),
		PARAMETER_TYPE_Double,
		100.0, 0.00001, true, 0.0, false
	);

	s.Printf("Node2", 2);
		pNode = Parameters.Add_Node(NULL,s,"Non-Linear Module","");

	Parameters.Add_Value(
		pNode,	"NL_ONLY",	_TL("Calibrate non-linear module only"),
		_TL("If checked, then only the non-linear module will be calibrated"),
		PARAMETER_TYPE_Bool, true
	);

	Parameters.Add_Value(
		pNode,	"CFAC_lb",	_TL("Parameter (c) [lower bound]"),
		_TL("Parameter (c) to fit streamflow volume"),
		PARAMETER_TYPE_Double,
		0.001, 0.0, true, 1.0, true
	);
	Parameters.Add_Value(
		pNode,	"CFAC_ub",	_TL("Parameter (c) [upper bound]"),
		_TL("Parameter (c) to fit streamflow volume"),
		PARAMETER_TYPE_Double,
		0.01, 0.0, true, 1.0, true
	);

	Parameters.Add_Value(
		pNode,	"TwFAC_lb",	_TL("(Tw) wetness decline time constant [lower bound]"),
		_TL("Tw is approximately the time constant, or inversely,"
		"the rate at which the catchment wetness declines in the absence of rainfall"),
		PARAMETER_TYPE_Double,
		1.0, 0.01, true, 150.0, true
	);
	Parameters.Add_Value(
		pNode,	"TwFAC_ub",	_TL("(Tw) wetness decline time constant [upper bound]"),
		_TL("Tw is approximately the time constant, or inversely,"
		"the rate at which the catchment wetness declines in the absence of rainfall"),
		PARAMETER_TYPE_Double,
		1.0, 0.01, true, 150.0, true
		);

	Parameters.Add_Value(
		pNode,	"DEV_ERAINDIS",	_TL("Maximum Deviation: ExcessRain-Streamflow [%]"),
		_TL("Maximum Deviation between estimated excess rainfall and observed streamflow volumes"),
		PARAMETER_TYPE_Double,
		5.0, 0.0, true, 50.0, true
	);

	Parameters.Add_Value(
		pNode,	"DEV_ERAINPCP",	_TL("Maximum Overestimation of ExcessRain [%]"),
		_TL("Maximum Sum of overestimated ExcessRain to total ExcessRain in [%]"),
		PARAMETER_TYPE_Double,
		5.0, 0.0, true, 50.0, true
	);

	s.Printf("Node3", 3);
	pNode = Parameters.Add_Node(NULL,s,"Storage Configuration","");

	Parameters.Add_Choice(
		pNode	, "STORAGE"		, "Storage",
		"",
		CSG_String::Format("%s|%s|%s|",
			_TL("Single Storage"),			//  0
			_TL("Two Parallel Storages"),	//  1 
			_TL("Two Storages in Series")	//  2 
		)
	);

	s.Printf("Node4", 4);
	pNode = Parameters.Add_Node(NULL,s,"IHACRES Version","");

	Parameters.Add_Choice(
		pNode	, "IHACVERS"		, "IHACRES Version",
		"",
		CSG_String::Format("%s|%s|",
			_TL("Jakeman & Hornberger (1993)"),	//  0
			_TL("Croke et al. (2005)")			//	1
		)
	);

	Parameters.Add_Value(
		pNode,	"SNOW_MODULE",	_TL("Snow Module on/off"),
		_TL("If checked the snow module is active"),
		PARAMETER_TYPE_Bool, true
	);

	Parameters.Add_Table_Output(
		NULL	, "TABLEout"	, _TL("Table"),
		""
	);
	Parameters.Add_Table_Output(
		NULL	, "TABLEparms"	, _TL("Table"),
		""
	);
}
//-------------------------------------------------------------------

bool Cihacres_cal::_CreateDialog2(bool bTMP, bool bnl_only, int storconf, int IHAC_version, bool bSnowModule,
								 CSG_Table *pTable, int dateField, int dischargeField)
{
	CSG_String		s;
	CSG_Parameters	P;  // used to add Parameters in the second dialog
	CSG_Parameter	*pNode, *pNode1;


	//	Dialog design
	P.Set_Name(_TL("Choose Time Range"));

	s.Printf("Node1", 1);
	pNode = P.Add_Node(NULL,s,"Time Range","");

	s.Printf("FDAY" , 1-1);
	P.Add_String(pNode,s,_TL("First Day"),_TL(""),
				 pTable->Get_Record(0)->asString(dateField));

	s.Printf("LDAY" , 1-2);
	P.Add_String(pNode,s,_TL("Last Day"),_TL(""),
				 pTable->Get_Record(pTable->Get_Record_Count()-1)->asString(dateField));

	if (bTMP) // if observed temperature data are used
	{
		s.Printf("Node2", 2);
		pNode = P.Add_Node(NULL,s,"Non-Linear Module","");

		s.Printf("Node2-1", 2-1);
		pNode1 = P.Add_Node(pNode,s,"Temperature Modulation Factor","");

		P.Add_Value(
			pNode1	, "TFac_lb"		, "Temperature Modulation Factor (f) [lower bound]",
			_TL("Temperature Modulation Factor f"),
			PARAMETER_TYPE_Double,
			0.05, 0.0001, true, 5.0, true
		);
		P.Add_Value(
			pNode1	, "TFac_ub"		, "Temperature Modulation Factor (f) [upper bound]",
			_TL("Temperature Modulation Factor f"),
			PARAMETER_TYPE_Double,
			0.5, 0.0001, true, 5.0, true
		);
	} // End if (bTMP)

	if (IHAC_version == 1) // Croke et al. (2005) Redesign
	{
		s.Printf("Node2-2", 2-2);
		pNode1 = P.Add_Node(pNode,s,"Soil moisture index threshold","");

		P.Add_Value(
			pNode1	, "SMI_L_lb"		, "Soil moisture index (l) [lower bound]",
			_TL("Soil moisture index threshold (l), lower bound"),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 5.0, true
		);
		P.Add_Value(
			pNode1	, "SMI_L_ub"		, "Soil moisture index (l) [upper bound]",
			_TL("Soil moisture index threshold (l), upper bound"),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 5.0, true
		);
		s.Printf("Node2-3", 2-3);
		pNode1 = P.Add_Node(pNode,s,"Power on soil moisture","");

		P.Add_Value(
			pNode1	, "SM_P_lb"		, "Power on soil (p) [lower bound]",
			_TL("Power on soil moisture (p), lower bound"),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 5.0, true
		);
		P.Add_Value(
			pNode1	, "SM_P_ub"		, "Power on soil (p) [upper bound]",
			_TL("Power on soil moisture (p), upper bound"),
			PARAMETER_TYPE_Double,
			0.0, 0.0, true, 5.0, true
		);
	}

	if (bSnowModule) // if snow module is active
	{
		s.Printf("Node2-4", 2-4);
		pNode1 = P.Add_Node(pNode,s,"Snow Module","");

		P.Add_Value(
			pNode1,	"T_RAIN_lb",	"Temperature Threshold for Rainfall (lower bound)",
			"Below this threshold precipitation will fall as snow",
			PARAMETER_TYPE_Double,
			-1.0, -10.0, true, 10.0, true
		);
		P.Add_Value(
			pNode1,	"T_RAIN_ub",	"Temperature Threshold for Rainfall (upper bound)",
			"Below this threshold precipitation will fall as snow",
			PARAMETER_TYPE_Double,
			1.0, -10.0, true, 10.0, true
		);
		P.Add_Value(
			pNode1,	"T_MELT_lb",	"Temperature Threshold for Melting (lower bound)",
			"Above this threshold snow will start to melt",
			PARAMETER_TYPE_Double,
			-1.0, -5.0, true, 10.0, true
		);
		P.Add_Value(
			pNode1,	"T_MELT_ub",	"Temperature Threshold for Melting (upper bound)",
			"Above this threshold snow will start to melt",
			PARAMETER_TYPE_Double,
			1.0, -5.0, true, 10.0, true
		);
		P.Add_Value(
			pNode1,	"DD_FAC_lb",	"Day-Degree Factor (lower bound)",
			"Day-Degree Factor depends on catchment characteristics",
			PARAMETER_TYPE_Double,
			0.7, 0.7, true, 9.2, true
		);
		P.Add_Value(
			pNode1,	"DD_FAC_ub",	"Day-Degree Factor (upper bound)",
			"Day-Degree Factor depends on catchment characteristics",
			PARAMETER_TYPE_Double,
			9.2, 0.7, true, 9.2, true
		);
	}

	if (!bnl_only) // if not only non-linear module calibration
	{
		switch(storconf)
		{
		case 0: // single storage
			s.Printf("Node3", 3);
			pNode = P.Add_Node(NULL,s,"Linear Module","");

			s.Printf("Node3-1", 3-1);
			pNode1 = P.Add_Node(pNode,s,"Parameter (a)","");

			P.Add_Value(
				pNode1,	"AFAC_lb",	_TL("(a) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.5, -0.99, true, -0.01, true
			);
			P.Add_Value(
				pNode1,	"AFAC_ub",	_TL("(a) [upper bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.9, -0.99, true, -0.01, true
				);

			s.Printf("Node3-2", 3-2);
			pNode1 = P.Add_Node(pNode,s,"Parameter (b)","");

			P.Add_Value(
				pNode1,	"BFAC_lb",	_TL("(b) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				0.0, 0.001, true, 1.0, true
			);
			P.Add_Value(
				pNode1,	"BFAC_ub",	_TL("(b) [upper bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				1.0, 0.001, true, 1.0, true
			);
			break;

		case 1: // two parallel storages
			s.Printf("Node3", 3);
			pNode = P.Add_Node(NULL,s,"Linear Module","");

			// Parameter a
			s.Printf("Node3-1", 3-1);
			pNode1 = P.Add_Node(pNode,s,"Parameter (a)","");

			P.Add_Value(
				pNode1,	"AQ_lb",	_TL("a(q) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.5, -0.99, true, -0.01, true
			);
			P.Add_Value(
				pNode1,	"AQ_ub",	_TL("a(q) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.9, -0.99, true, -0.01, true
			);
			P.Add_Value(
				pNode1,	"AS_lb",	_TL("a(s) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.5, -0.99, true, -0.01, true
			);
			P.Add_Value(
				pNode1,	"AS_ub",	_TL("a(s) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				-0.9, -0.99, true, -0.01, true
			);

			// Parameter b
			s.Printf("Node3-2", 3-2);
			pNode1 = P.Add_Node(pNode,s,"Parameter (b)","");

			P.Add_Value(
				pNode1,	"BQ_lb",	_TL("b(q) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				0.0, 0.0, true, 1.0, true
			);
			P.Add_Value(
				pNode1,	"BQ_ub",	_TL("b(q) [lower bound]"),
				_TL(""),
				PARAMETER_TYPE_Double,
				1.0, 0.0, true, 1.0, true
			);
			//P.Add_Value(
			//	pNode1,	"BS_lb",	_TL("b(s) [lower bound]"),
			//	_TL(""),
			//	PARAMETER_TYPE_Double,
			//	0.0, 0.0, true, 1.0, true
			//);
			//P.Add_Value(
			//	pNode1,	"BS_ub",	_TL("b(s) [lower bound]"),
			//	_TL(""),
			//	PARAMETER_TYPE_Double,
			//	1.0, 0.0, true, 1.0, true
			//);

			break;
		case 2: // two storages in series
			break;
		} // end switch (storconf)

		s.Printf("Node4", 4);
		pNode1 = P.Add_Node(pNode,s,"Time Delay after Start of Rainfall (INTEGER)","");
		
		P.Add_Value(
			pNode1,	"DELAY",	"Time Delay (Rain-Runoff)",
			"The delay after the start of rainfall, before the discharge starts to rise.",
			PARAMETER_TYPE_Int,
			0, 1, true, 100, true
		);

		s.Printf("Node6", 6);
		pNode1 = P.Add_Node(pNode,s,"Nash-Sutcliffe Efficiency","");

		P.Add_Choice(
			pNode1, "OBJ_FUNC"		, "Objective Function",
			"",
			"NSE|NSE high flow|NSE low flow"
		);

		P.Add_Value(
			pNode1,	"NSEMIN",	"Minimum Nash-Sutcliffe Efficiency",
			"Minimum Nash-Sutcliffe Efficiency required to print simulation to calibration table",
			PARAMETER_TYPE_Double,
			0.7, 0.1, true, 1.0, true
		);

	} // End if (!bnl_only)
	// End of second dialog design


	if( SG_UI_Dlg_Parameters(&P, _TL("Choose Time Range")) )
	{
		// assign data from second dialog

		m_date1		= P(CSG_String::Format("FDAY",
										dateField).c_str())->asString();
		m_date2		= P(CSG_String::Format("LDAY",
										dischargeField).c_str())->asString();
		if (bTMP) // if observed temperature data are used
		{
			m_f_lb	= P("TFac_lb")			->asDouble();
			m_f_ub	= P("TFac_ub")			->asDouble();
		}
		if (IHAC_version == 1)
		{
			m_l_lb	= P("SMI_L_lb")			->asDouble();
			m_l_ub	= P("SMI_L_ub")			->asDouble();
			m_p_lb	= P("SM_P_lb")			->asDouble();
			m_p_ub	= P("SM_P_ub")			->asDouble();
		}

		if (!bnl_only) // if not only non-linear module calibration
					   // (if non-linear and linear modulde calibration)
		{
			switch (storconf)
			{
			case 0:
				m_a_lb	= P("AFAC_lb")			->asDouble();
				m_a_ub	= P("AFAC_ub")			->asDouble();
				m_b_lb	= P("BFAC_lb")			->asDouble();
				m_b_ub	= P("BFAC_ub")			->asDouble();
				m_delay	= P("DELAY")			->asInt();
				m_NSEmin	= P("NSEMIN")			->asDouble();
				break;
			case 1: // two parallel storages
				m_aq_lb	= P("AQ_lb")			->asDouble();
				m_aq_ub	= P("AQ_ub")			->asDouble();
				m_as_lb	= P("AS_lb")			->asDouble();
				m_as_ub	= P("AS_ub")			->asDouble();
				m_bq_lb	= P("BQ_lb")			->asDouble();
				m_bq_ub	= P("BQ_ub")			->asDouble();
				//m_bs_lb	= P("BS_lb")			->asDouble();
				//m_bs_ub	= P("BS_ub")			->asDouble();
				m_delay		= P("DELAY")		->asInt();
				m_obj_func	= P("OBJ_FUNC")		->asInt();
				m_NSEmin	= P("NSEMIN")		->asDouble();
				break;
			case 2: // two storages in series
				break;
			} // end switch

		} // end if (!bnl_only)
		if (bSnowModule)
		{
			m_T_Rain_lb	= P("T_RAIN_lb")		->asDouble();
			m_T_Rain_ub	= P("T_RAIN_ub")		->asDouble();
			m_T_Melt_lb	= P("T_MELT_lb")		->asDouble();
			m_T_Melt_ub	= P("T_MELT_ub")		->asDouble();
			m_DD_FAC_lb	= P("DD_FAC_lb")		->asDouble();
			m_DD_FAC_ub	= P("DD_FAC_ub")		->asDouble();
		}
		return(true);
	} else {
		return(false);
	}
}




