1 /**********************************************************
2 * Version $Id: ihacres_eq.cpp 911 2011-02-14 16:38:15Z reklov_w $
3 *********************************************************/
4 ///////////////////////////////////////////////////////////
5 // //
6 // ihacres_eq.cpp //
7 // //
8 // Copyright (C) 2006 by //
9 // Stefan Liersch //
10 //-------------------------------------------------------//
11 // //
12 // e-mail: stefan.liersch@ufz.de //
13 // stefan.liersch@gmail.com //
14 // //
15 // 2006-08-27 //
16 // //
17 //-------------------------------------------------------//
18
19 #include <iostream> // used for textfile output (test only)
20 #include <fstream> // used for textfile output (test only)
21
22 #include "ihacres_eq.h"
23 #include <math.h> // exp()
24 //---------------------------------------------------------
25
26 ///////////////////////////////////////////////////////////
27 //
28 // COMMENTS
29 //
30 //---------------------------------------------------------
31 // 2007-11-08
32 // added to function CalcWetnessIndex()
33 // if (WetnessIndex[i] > 1.0) WetnessIndex[i] = 1.0;
34 //---------------------------------------------------------
35 ///////////////////////////////////////////////////////////
36
37
38 ///////////////////////////////////////////////////////////////////////
39 //
40 // CONSTRUCTORS
41 //
42 ///////////////////////////////////////////////////////////////////////
43
44 //---------------------------------------------------------------------
45 // DEFAULT CONSTRUCTOR
46 //---------------------------------------------------------------------
Cihacres_eq()47 Cihacres_eq::Cihacres_eq()
48 {};
49
50 //---------------------------------------------------------------------
51 // two storages
52 // using vector<double> as input
53 // if no temperature data are available
Cihacres_eq(date_array date_in,vector_d streamflow,vector_d pcp,double TwConst,double f,double c,double l,double p,double aq,double as,double bq,double bs)54 Cihacres_eq::Cihacres_eq(date_array date_in,
55 vector_d streamflow,
56 vector_d pcp,
57 double TwConst, double f, double c,
58 double l, double p,
59 double aq, double as, double bq, double bs)
60 {
61 sizeAll = (int)streamflow.size();
62 date = date_in;
63 streamflow_obs = streamflow;
64 precipitation = pcp;
65 this->TwConst = TwConst;
66 this->f = f;
67 this->c = c;
68 this->l = l;
69 this->p = p;
70 this->aq = aq;
71 this->as = as;
72 this->bq = bq;
73 this->bs = bs;
74 // Initialize Vectors
75 _InitVectorsStart((int)streamflow_obs.size());
76 }
77 //---------------------------------------------------------------------
78 // two storages
79 // using vector<double> as input
80 // if temperature data are available
Cihacres_eq(date_array date_in,vector_d streamflow,vector_d pcp,vector_d tmp,double TwConst,double f,double c,double l,double p,double aq,double as,double bq,double bs,double area,bool TMP_data_exist,int IHAC_vers,int storconf,bool bSnowModule,CSnowModule * SnowMod,int delay)81 Cihacres_eq::Cihacres_eq(date_array date_in,
82 vector_d streamflow,
83 vector_d pcp,
84 vector_d tmp,
85 double TwConst, double f, double c,
86 double l, double p,
87 double aq, double as, double bq, double bs,
88 double area, bool TMP_data_exist,
89 int IHAC_vers,
90 int storconf,
91 bool bSnowModule,
92 CSnowModule *SnowMod,
93 //double T_Rain, double T_Melt, double DD_FAC,
94 int delay)
95 {
96 // Initialize Parameters and Vectors
97 sizeAll = (int)streamflow.size();
98 date = date_in;
99 streamflow_obs = streamflow;
100 precipitation = pcp;
101 temperature = tmp;
102 this->TwConst = TwConst;
103 this->f = f;
104 this->c = c;
105 this->l = l;
106 this->p = p;
107 this->aq = aq;
108 this->as = as;
109 this->bq = bq;
110 this->bs = bs;
111 this->delay = delay;
112 this->area = area;
113 IHAC_version = IHAC_vers;
114 this->bSnowModule = bSnowModule;
115 m_pSnowMod = SnowMod;
116
117 // Initialize Vectors containing calculated values
118 _InitVectorsStart(sizeAll);
119
120 // Convert Streamflow vector from m3/s*day-1 to mm/day
121 streamflowMM_obs = model_tools::m3s_to_mmday(streamflow_obs, streamflowMM_obs, area);
122
123 // perform simulation
124 if (bSnowModule)
125 {
126 RunNonLinearModule(TMP_data_exist, bSnowModule, m_pSnowMod->Get_T_Rain());
127 } else {
128 RunNonLinearModule(TMP_data_exist, bSnowModule, 0.0);
129 }
130 //switch (IHAC_version)
131 //{
132 //case 0: // Jakeman & Hornberger (1993)
133 // if (bSnowModule)
134 // {
135 // RunNonLinearModule(TMP_data_exist, bSnowModule, m_pSnowMod->Get_T_Rain());
136 // } else {
137 // RunNonLinearModule(TMP_data_exist, bSnowModule, 0.0);
138 // }
139 // break;
140 //case 1: // Croke et al. (2005) Redesign
141 // //RunNonLinearModule5Parms();
142 // break;
143 //}
144
145 switch(storconf)
146 {
147 case 0: // single storage
148 this->a = aq;
149 this->b = bq;
150 SimStreamflowSingle(excessRain, streamflowMM_obs[0], streamflow_sim, delay, a, b);
151 break;
152 case 1: // two storages in parallel
153 SimStreamflow2Parallel(excessRain, streamflow_sim,
154 streamflowMM_obs[0],
155 aq, as, bq, bs,
156 vq, vs,
157 IHAC_vers,
158 delay);
159 break;
160 } // end switch(storconf)
161
162 NSE = model_tools::CalcEfficiency(streamflowMM_obs, streamflow_sim);
163 }
164 //---------------------------------------------------------------------
165 // two storages
166 // using double arrays as input
167 // if no temperature data are available
Cihacres_eq(int size,date_array date_in,double * streamflow,double * pcp,double TwConst,double f,double c,double aq,double as,double bq,double bs)168 Cihacres_eq::Cihacres_eq(int size, // array size
169 date_array date_in,
170 double *streamflow,
171 double *pcp,
172 double TwConst, double f, double c,
173 double aq, double as, double bq, double bs)
174 {
175 // assign values to global parameters
176 sizeAll = size;
177 date = date_in;
178 streamflow_obs.resize(sizeAll);
179 precipitation.resize(sizeAll);
180 for (int i = 0; i < sizeAll; i++)
181 {
182 streamflow_obs[i] = streamflow[i];
183 precipitation[i] = pcp[i];
184 }
185 this->TwConst = TwConst;
186 this->f = f;
187 this->c = c;
188 this->aq = aq;
189 this->as = as;
190 this->bq = bq;
191 this->bs = bs;
192 // Initialize Vectors
193 _InitVectorsStart(sizeAll);
194 }
195 //---------------------------------------------------------------------
196 // two storages
197 // using double arrays as input
198 // if temperature data are available
Cihacres_eq(int size,date_array date_in,double * streamflow,double * pcp,double * tmp,double TwConst,double f,double c,double aq,double as,double bq,double bs)199 Cihacres_eq::Cihacres_eq(int size, // array size
200 date_array date_in,
201 double *streamflow,
202 double *pcp,
203 double *tmp,
204 double TwConst, double f, double c,
205 double aq, double as, double bq, double bs)
206 {
207 sizeAll = size;
208 date = date_in;
209 streamflow_obs.resize(size);
210 precipitation.resize(size);
211 temperature.resize(size);
212 for (int i = 0; i < size; i++)
213 {
214 streamflow_obs[i] = streamflow[i];
215 precipitation[i] = pcp[i];
216 temperature[i] = tmp[i];
217 }
218 this->TwConst = TwConst;
219 this->f = f;
220 this->c = c;
221 this->aq = aq;
222 this->as = as;
223 this->bq = bq;
224 this->bs = bs;
225 // Initialize Vectors
226 _InitVectorsStart((int)streamflow_obs.size());
227 }
228 //---------------------------------------------------------------------
229 // end constructors ///////////////////////////////////////////////////
230
231 // destructor
~Cihacres_eq(void)232 Cihacres_eq::~Cihacres_eq(void)
233 {
234 _ZeroAllVectors();
235 }
236
237 ///////////////////////////////////////////////////////////////////////
238 //
239 // PUBLIC MEMBER FUNCTIONS
240 //
241 ///////////////////////////////////////////////////////////////////////
242
243 //---------------------------------------------------------------------
244 // Run Non-Linear Tool
245 //---------------------------------------------------------------------
RunNonLinearModule(bool TMP_data_exist,bool bSnowModule,double T_Rain)246 void Cihacres_eq::RunNonLinearModule(bool TMP_data_exist, bool bSnowModule, double T_Rain)
247 {
248 double WI_init = 0.5;
249 double eR_init = 0.0;
250
251 switch (IHAC_version)
252 {
253 case 0: // Jakeman & Hornberger (1993)
254 // if temperature data are available (TMP_data_exist = true), then adjust the rate
255 // at which the catchment wetness declines in the absence of rainfall
256 // to daily temperature
257 if (TMP_data_exist){
258 CalcWetnessTimeConst(temperature, Tw, TwConst, f);
259 }
260 if (bSnowModule)
261 {
262 // calculate the catchment wetness index
263 CalcWetnessIndex(Tw, precipitation, temperature,
264 WetnessIndex, WI_init, c,
265 bSnowModule, m_pSnowMod->Get_T_Rain());
266 // calculate effective rainfall
267 sum_eRainMM = CalcExcessRain(precipitation, temperature, WetnessIndex, excessRain, eR_init,
268 sum_eRainGTpcp, bSnowModule, m_pSnowMod);
269 } else {
270 CalcWetnessIndex(Tw, precipitation, temperature, WetnessIndex, WI_init, c,
271 bSnowModule, 0.0);
272 sum_eRainMM = CalcExcessRain(precipitation, temperature, WetnessIndex, excessRain, eR_init,
273 sum_eRainGTpcp, bSnowModule, m_pSnowMod);
274 }
275 break;
276
277 case 1: // Croke et al. (2005)
278 if (TMP_data_exist) {
279 CalcWetnessTimeConst_Redesign(temperature, Tw, TwConst, f);
280 }
281 if (bSnowModule)
282 {
283 // calculate the catchment wetness index
284 CalcWetnessIndex_Redesign(Tw, precipitation, WetnessIndex, bSnowModule, m_pSnowMod->Get_T_Rain());
285 // calculate effective rainfall
286 sum_eRainMM = CalcExcessRain_Redesign(precipitation, temperature,
287 WetnessIndex, excessRain, eR_init, sum_eRainGTpcp,
288 c, l, p, bSnowModule, m_pSnowMod);
289 } else {
290 // calculate the catchment wetness index
291 CalcWetnessIndex_Redesign(Tw, precipitation, WetnessIndex, bSnowModule, 0.0);
292 // calculate effective rainfall
293 sum_eRainMM = CalcExcessRain_Redesign(precipitation, temperature,
294 WetnessIndex, excessRain, eR_init, sum_eRainGTpcp,
295 c, l, p, bSnowModule, m_pSnowMod);
296 }
297 break;
298 }
299 }
300
301 //---------------------------------------------------------------------
302 // Simulate Streamflow (single storage)
303 //---------------------------------------------------------------------
SimStreamflowSingle(vector_d & excessRain,double initVal,vector_d & streamflow_sim,int delay,double a,double b)304 void Cihacres_eq::SimStreamflowSingle(vector_d &excessRain, double initVal,
305 vector_d &streamflow_sim, int delay,
306 double a, double b)
307 {
308 int i;
309 int size = (int)streamflow_sim.size();
310 // using the first observed streamflow value as initial simulation value
311 for (i = 0; i < delay; i++)
312 streamflow_sim[i] = initVal;
313 // start calculation with second value
314 for (i = delay; i < size; i++)
315 {
316 streamflow_sim[i] = -a * streamflow_sim[i-1] + b * excessRain[i-delay];
317 }
318 }
319 //---------------------------------------------------------------------
320 // Simulate Streamflow (single storage)
321 //---------------------------------------------------------------------
SimStreamflowSingle(double * excessRain,double initVal,double * streamflow_sim,int delay,double a,double b,int size)322 void Cihacres_eq::SimStreamflowSingle(double *excessRain, double initVal,
323 double *streamflow_sim, int delay,
324 double a, double b, int size)
325 {
326 int i;
327 // using the first observed streamflow value as initial simulation value
328 for (i = 0; i < delay; i++)
329 streamflow_sim[i] = initVal;
330 // start calculation with second value
331 for (i = delay; i < size; i++)
332 {
333 streamflow_sim[i] = -a * streamflow_sim[i-1] + b * excessRain[i-delay];
334 }
335 }
336
337 //---------------------------------------------------------------------
338 // Simulate Streamflow (2 parallel storages)
339 //---------------------------------------------------------------------
SimStreamflow2Parallel(vector_d & excessRain,vector_d & streamflow_sim,double initVal,double aq,double as,double bq,double bs,double & vq,double & vs,int IHAC_vers,int delay)340 void Cihacres_eq::SimStreamflow2Parallel(vector_d &excessRain, vector_d &streamflow_sim,
341 double initVal,
342 double aq, double as, double bq, double bs,
343 double &vq, double &vs,
344 int IHAC_vers, int delay)
345 {
346 int i;
347 int size = (int)streamflow_sim.size();
348 double *sf_q = new double[size]; // quick streamflow component
349 double *sf_s = new double[size]; // slow streamflow component
350
351 // calculate the dependent b-value
352 // after Jakeman etc.
353 vq = bq / (1 + aq);
354 vs = 1 - vq;
355
356 // using the first observed streamflow value as initial simulation value
357 for (i = 0; i < delay; i++)
358 {
359 streamflow_sim[i] = initVal;
360 sf_q[i] = initVal * vq;
361 sf_s[i] = initVal * vs;
362 }
363
364 // using the first observed streamflow value as initial simulation value
365 //for (i = 0; i < delay; i++)
366 //{
367 // streamflow_sim[i] = initVal;
368 // sf_q[i] = initVal / 2;
369 // sf_s[i] = initVal / 2;
370 //}
371
372 //// calculate the dependent b-value
373 //if (IHAC_vers == 1) // after Kokkonen
374 //{
375 // vq = bq / (bq + bs);
376 // vs = 1 - vq;
377 //} else { // after Jakeman etc.
378 // vq = bq / (1 + aq);
379 // vs = 1 - vq;
380 //}
381
382 // calculate streamflow
383 for (i = delay; i < size; i++)
384 {
385 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
386 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
387 streamflow_sim[i] = sf_q[i] + sf_s[i];
388 }
389 /*
390 switch(IHAC_vers)
391 {
392 case 0 : // after Jakeman & Hornberger (1993)
393 vq = bq / (1 + aq);
394 vs = 1 - vq;
395 if (!b_freebee)
396 bs = vs * (1 + as);
397
398 // calculate quick and slow components
399 for (i = delay; i < size; i++)
400 {
401 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
402 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
403 streamflow_sim[i] = sf_q[i] + sf_s[i];
404 }
405 break;
406 case 1 : // after Kokkonen et al. (2003)
407 vq = bq / (bq + bs);
408 vs = 1 - vq;
409 if (!b_freebee)
410 bs = vs * (1 + as);
411 // calculate quick and slow components
412 for (i = delay; i < size; i++)
413 {
414 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
415 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
416 streamflow_sim[i] = sf_q[i] + sf_s[i];
417 }
418 break;
419 case 2: // Ahandere Mahn
420 vq = bq / (1 + aq);
421 vs = 1 - vq;
422 if (!b_freebee)
423 bs = vs * (1 + as);
424
425 // calculate quick and slow components
426 for (i = delay; i < size; i++)
427 {
428 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
429 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
430 streamflow_sim[i] = sf_q[i] + sf_s[i];
431 }
432 break;
433 case 3 : // after Croke et al. (2005)
434 vq = bq / (1 + aq);
435 vs = 1 - vq;
436 if (!b_freebee)
437 bs = vs * (1 + as);
438
439 // calculate quick and slow components
440 for (i = delay; i < size; i++)
441 {
442 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
443 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
444 streamflow_sim[i] = sf_q[i] + sf_s[i];
445 }
446 break;
447 }
448 */
449 delete[] sf_q;
450 delete[] sf_s;
451 }
452 //---------------------------------------------------------------------
453 // Simulate Streamflow (2 parallel storages)
454 //---------------------------------------------------------------------
SimStreamflow2Parallel(double * excessRain,double * streamflow_sim,double initVal,double aq,double as,double bq,double bs,double & vq,double & vs,int IHAC_vers,int size,int delay)455 void Cihacres_eq::SimStreamflow2Parallel(double *excessRain, double *streamflow_sim,
456 double initVal,
457 double aq, double as, double bq, double bs,
458 double &vq, double &vs,
459 int IHAC_vers, int size, int delay)
460 {
461 int i;
462 double *sf_q = new double[size]; // quick streamflow component
463 double *sf_s = new double[size]; // slow streamflow component
464
465 // calculate the dependent b-value
466 // after Jakeman etc.
467 vq = bq / (1 + aq);
468 vs = 1 - vq;
469
470 // using the first observed streamflow value as initial simulation value
471 for (i = 0; i < delay; i++)
472 {
473 streamflow_sim[i] = initVal;
474 sf_q[i] = initVal * vq;
475 sf_s[i] = initVal * vs;
476 }
477
478 //// using the first observed streamflow value as initial simulation value
479 //for (i = 0; i < delay; i++)
480 //{
481 // streamflow_sim[i] = initVal;
482 // sf_q[i] = initVal / 2;
483 // sf_s[i] = initVal / 2;
484 //}
485
486 //// calculate the dependent b-value
487 //if (IHAC_vers == 1) // after Kokkonen
488 //{
489 // vq = bq / (bq + bs);
490 // vs = 1 - vq;
491 //} else { // after Jakeman etc.
492 // vq = bq / (1 + aq);
493 // vs = 1 - vq;
494 //}
495
496 // calculate streamflow
497 for (i = delay; i < size; i++)
498 {
499 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
500 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
501 streamflow_sim[i] = sf_q[i] + sf_s[i];
502 }
503 /*
504 switch(IHAC_vers)
505 {
506 case 0 : // after Jakeman & Hornberger (1993)
507 vq = bq / (1 + aq);
508 vs = 1 - vq;
509 if (!b_freebee)
510 bs = vs * (1 + as);
511
512 // calculate quick and slow components
513 for (i = delay; i < size; i++)
514 {
515 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
516 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
517 streamflow_sim[i] = sf_q[i] + sf_s[i];
518 }
519 break;
520 case 1 : // after Kokkonen et al. (2003)
521 vq = bq / (bq + bs);
522 vs = 1 - vq;
523 if (!b_freebee)
524 bs = vs * (1 + as);
525 // calculate quick and slow components
526 for (i = delay; i < size; i++)
527 {
528 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
529 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
530 streamflow_sim[i] = sf_q[i] + sf_s[i];
531 }
532 break;
533 case 2: // Jakeman & SnowModule
534 vq = bq / (1 + aq);
535 vs = 1 - vq;
536 if (!b_freebee)
537 bs = vs * (1 + as);
538
539 // calculate quick and slow components
540 for (i = delay; i < size; i++)
541 {
542 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
543 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
544 streamflow_sim[i] = sf_q[i] + sf_s[i];
545 }
546 break;
547 case 3 : // after Croke et al. (2005)
548 vq = bq / (1 + aq);
549 vs = 1 - vq;
550 if (!b_freebee)
551 bs = vs * (1 + as);
552
553 // calculate quick and slow components
554 for (i = delay; i < size; i++)
555 {
556 sf_q[i] = -aq * sf_q[i-1] + bq * excessRain[i-delay];
557 sf_s[i] = -as * sf_s[i-1] + bs * excessRain[i-delay];
558 streamflow_sim[i] = sf_q[i] + sf_s[i];
559 }
560 break;
561 }
562 */
563 delete[] sf_q;
564 delete[] sf_s;
565 }
566
567 //---------------------------------------------------------------------
568 // Simulate Streamflow (2 parallel storages)
569 //---------------------------------------------------------------------
SimStreamflow2Parallel(double * excessRain,double * streamflow_sim,double initVal,C_IHAC_LinearParms * linparms,int index,double & vq,double & vs,int size,int delay)570 void Cihacres_eq::SimStreamflow2Parallel(double *excessRain, double *streamflow_sim,
571 double initVal, // first observed streamflow value
572 C_IHAC_LinearParms* linparms, int index,
573 double &vq, double &vs, int size, int delay)
574 {
575 int i;
576 double *sf_q = new double[size]; // quick streamflow component
577 double *sf_s = new double[size]; // slow streamflow component
578
579 // calculate the dependent b-value
580 // after Jakeman etc.
581 vq = linparms->bq[index] / (1 + linparms->aq[index]);
582 vs = 1 - vq;
583
584 // using the first observed streamflow value as initial simulation value
585 for (i = 0; i < delay; i++)
586 {
587 streamflow_sim[i] = initVal;
588 sf_q[i] = initVal * vq;
589 sf_s[i] = initVal * vs;
590 }
591
592 // calculate streamflow
593 for (i = delay; i < size; i++)
594 {
595 sf_q[i] = -linparms->aq[index] * sf_q[i-1] + linparms->bq[index] * excessRain[i-delay];
596 sf_s[i] = -linparms->as[index] * sf_s[i-1] + linparms->bs[index] * excessRain[i-delay];
597 streamflow_sim[i] = sf_q[i] + sf_s[i];
598 }
599 delete[] sf_q;
600 delete[] sf_s;
601 }
602
603
604 //---------------------------------------------------------------------
605 // Calculate Parameter b(q)
606 //---------------------------------------------------------------------
Calc_Parm_BS(double aq,double as,double bq)607 double Cihacres_eq::Calc_Parm_BS(double aq, double as, double bq)
608 {
609 return( (1 - (bq / (1 + aq))) * (1 + as) );
610 }
611 //---------------------------------------------------------------------
612
613 //---------------------------------------------------------------------
614 // Calculate time of decay (quick or slow)
615 //---------------------------------------------------------------------
Calc_TimeOfDecay(double a)616 double Cihacres_eq::Calc_TimeOfDecay(double a)
617 {
618 return( -1 / log(-a) );
619 }
620 //---------------------------------------------------------------------
621
622 //---------------------------------------------------------------------
623 // Calculate Wetness Time Constant
624 //---------------------------------------------------------------------
CalcWetnessTimeConst(vector_d & temperature,vector_d & Tw,double TwConst,double f)625 void Cihacres_eq::CalcWetnessTimeConst(vector_d &temperature,
626 vector_d &Tw,
627 double TwConst,
628 double f)
629 {
630 for (unsigned int i = 0; i < Tw.size(); i++)
631 {
632 Tw[i] = TwConst * exp((20.0 - temperature[i]) * f);
633 }
634 }
635 //---------------------------------------------------------------------
636 // Calculate Wetness Time Constant
637 //---------------------------------------------------------------------
CalcWetnessTimeConst(double * temperature,double * Tw,double TwConst,double f,int size)638 void Cihacres_eq::CalcWetnessTimeConst(double *temperature,
639 double *Tw,
640 double TwConst,
641 double f,
642 int size)
643 {
644 for (int i = 0; i < size; i++)
645 {
646 Tw[i] = TwConst * exp((20.0 - temperature[i]) * f);
647 }
648 }
649
650 //---------------------------------------------------------------------
651 // Calculate Wetness Time Constant
652 //---------------------------------------------------------------------
CalcWetnessTimeConst(double * temperature,double * Tw,C_IHAC_NonLinearParms * nonlinparms,int index,int size)653 void Cihacres_eq::CalcWetnessTimeConst(double *temperature, double *Tw,
654 C_IHAC_NonLinearParms* nonlinparms, int index,
655 int size)
656 {
657 for (int i = 0; i < size; i++)
658 {
659 Tw[i] = nonlinparms->mp_tw[index] * exp((20.0 - temperature[i]) * nonlinparms->mp_f[index]);
660 }
661 }
662
663 //---------------------------------------------------------------------
664 // Calculate Wetness Time Constant
665 // For ihacres_climate_scen
666 //---------------------------------------------------------------------
CalcWetnessTimeConst_scen(double * temperature,double * Tw,C_IHAC_NonLinearParms * nonlinparms,int index,int size)667 void Cihacres_eq::CalcWetnessTimeConst_scen(double *temperature, double *Tw,
668 C_IHAC_NonLinearParms* nonlinparms, int index,
669 int size)
670 {
671 Tw[0] = 0.0;
672 for (int i = 1; i < size; i++)
673 {
674 Tw[i] = nonlinparms->mp_tw[index] * exp((20.0 - temperature[i]) * nonlinparms->mp_f[index]);
675 }
676 }
677
678 //---------------------------------------------------------------------
679 // Calculate Wetness Time Constant (Croke et al. 2005)
680 //---------------------------------------------------------------------
CalcWetnessTimeConst_Redesign(vector_d & temperature,vector_d & Tw,double TwConst,double f)681 void Cihacres_eq::CalcWetnessTimeConst_Redesign(vector_d &temperature,
682 vector_d &Tw,
683 double TwConst,
684 double f)
685 {
686 double Tr = 20.0; // reference temperature
687 for (unsigned int i = 0; i < Tw.size(); i++)
688 {
689 Tw[i] = TwConst * exp(0.062 * f * (Tr - temperature[i]));
690 }
691 }
692
693 //---------------------------------------------------------------------
694 // Calculate Wetness Time Constant (Croke et al. 2005)
695 //---------------------------------------------------------------------
CalcWetnessTimeConst_Redesign(double * temperature,double * Tw,double TwConst,double f,int size)696 void Cihacres_eq::CalcWetnessTimeConst_Redesign(double *temperature,
697 double *Tw,
698 double TwConst,
699 double f,
700 int size)
701 {
702 double Tr = 20.0; // reference temperature
703 for (int i = 0; i < size; i++)
704 {
705 Tw[i] = TwConst * exp(0.062 * f * (Tr - temperature[i]));
706 }
707 }
708 //---------------------------------------------------------------------
709 // Calculate Wetness Time Constant (Croke et al. 2005)
710 //---------------------------------------------------------------------
CalcWetnessTimeConst_Redesign(double * temperature,double * Tw,C_IHAC_NonLinearParms * nonlinparms,int index,int size)711 void Cihacres_eq::CalcWetnessTimeConst_Redesign(double *temperature, double *Tw,
712 C_IHAC_NonLinearParms* nonlinparms, int index,
713 int size)
714 {
715 double Tr = 20.0; // reference temperature
716 for (int i = 0; i < size; i++)
717 {
718 Tw[i] = nonlinparms->mp_tw[index] * exp(0.062 * nonlinparms->mp_f[index] * (Tr - temperature[i]));
719 }
720 }
721
722 //---------------------------------------------------------------------
723 // Calculate Wetness Index
724 //---------------------------------------------------------------------
CalcWetnessIndex(vector_d & Tw,vector_d & precipitation,vector_d & temperature,vector_d & WetnessIndex,double WI_init,double c,bool bSnowModule,double T_Rain)725 void Cihacres_eq::CalcWetnessIndex(vector_d &Tw,
726 vector_d &precipitation, vector_d &temperature,
727 vector_d &WetnessIndex, double WI_init,
728 double c, bool bSnowModule, double T_Rain)
729 {
730 WetnessIndex[0] = WI_init;
731 // starting at the second value (i=1)
732 for (unsigned int i = 1; i < WetnessIndex.size(); i++)
733 {
734 if (bSnowModule && temperature[i] < T_Rain)
735 {
736 WetnessIndex[i] = (1 - (1 / Tw[i])) * WetnessIndex[i-1];
737 } else {
738 WetnessIndex[i] = c * precipitation[i] + (1 - (1 / Tw[i])) * WetnessIndex[i-1];
739 }
740 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
741 // if (WetnessIndex[i] > 1.0) WetnessIndex[i] = 1.0;
742 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
743 }
744 }
745 //---------------------------------------------------------------------
746 // Calculate Wetness Index
747 //---------------------------------------------------------------------
CalcWetnessIndex(double * Tw,double * precipitation,double * temperature,double * WetnessIndex,double WI_init,double c,bool bSnowModule,double T_Rain,int size)748 void Cihacres_eq::CalcWetnessIndex(double *Tw,
749 double *precipitation, double *temperature,
750 double *WetnessIndex, double WI_init,
751 double c, bool bSnowModule, double T_Rain,
752 int size)
753 {
754 //WetnessIndex[0] = 0.5;
755 WetnessIndex[0] = WI_init;
756 // starting at the second value (i=1)
757 for (int i = 1; i < size; i++)
758 {
759 if (bSnowModule && temperature[i] < T_Rain)
760 {
761 WetnessIndex[i] = (1 - (1 / Tw[i])) * WetnessIndex[i-1];
762 } else {
763 WetnessIndex[i] = c * precipitation[i] + (1 - (1 / Tw[i])) * WetnessIndex[i-1];
764 }
765 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
766 // if (WetnessIndex[i] > 1.0) WetnessIndex[i] = 1.0;
767 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
768 }
769 }
770
771 //---------------------------------------------------------------------
772 // Calculate Wetness Index after Croke et al. (2005)
773 //---------------------------------------------------------------------
CalcWetnessIndex_Redesign(vector_d & Tw,vector_d & precipitation,vector_d & WetnessIndex,bool bSnowModule,double T_Rain)774 void Cihacres_eq::CalcWetnessIndex_Redesign(vector_d &Tw,
775 vector_d &precipitation,
776 vector_d &WetnessIndex, bool bSnowModule,double T_Rain)
777 {
778 WetnessIndex[0] = 0.5;
779 // starting at the second value (i=1)
780 for (unsigned int i = 1; i < WetnessIndex.size(); i++)
781 {
782 WetnessIndex[i] = precipitation[i] + (1 - (1 / Tw[i])) * WetnessIndex[i-1];
783 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
784 // if (WetnessIndex[i] > 1.0) WetnessIndex[i] = 1.0;
785 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
786 }
787 }
788 //---------------------------------------------------------------------
789 // Calculate Wetness Index after Croke et al. (2005)
790 //---------------------------------------------------------------------
CalcWetnessIndex_Redesign(double * Tw,double * precipitation,double * WetnessIndex,double WI_init,bool bSnowModule,double T_Rain,int size)791 void Cihacres_eq::CalcWetnessIndex_Redesign(double *Tw,
792 double *precipitation,double *WetnessIndex, double WI_init,
793 bool bSnowModule,double T_Rain,
794 int size)
795 {
796 //WetnessIndex[0] = 0.5;
797 WetnessIndex[0] = WI_init;
798 // starting at the second value (i=1)
799 for (int i = 1; i < size; i++)
800 {
801 WetnessIndex[i] = precipitation[i] + (1 - (1 / Tw[i])) * WetnessIndex[i-1];
802 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
803 // if (WetnessIndex[i] > 1.0) WetnessIndex[i] = 1.0;
804 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
805 }
806 }
807 //---------------------------------------------------------------------
808 // Calculate Effective or Excess Rainfall
809 //---------------------------------------------------------------------
CalcExcessRain(vector_d & precipitation,vector_d & temperature,vector_d & WetnessIndex,vector_d & excessRain,double eR_init,double & sum_eRainGTpcp,bool bSnowModule,CSnowModule * pSnowModule)810 double Cihacres_eq::CalcExcessRain(vector_d &precipitation, vector_d &temperature, vector_d &WetnessIndex,
811 vector_d &excessRain, double eR_init, double &sum_eRainGTpcp,
812 bool bSnowModule, CSnowModule* pSnowModule)
813 //double T_Rain, double T_Melt, double* MeltRate)
814 {
815 double sum = 0.0; // sum of total ExcessRain of the period
816 sum_eRainGTpcp = 0.0;
817
818 excessRain[0] = eR_init;
819 if (precipitation[0] > 0.0) excessRain[0] = precipitation[0] / 2;
820 // starting at the second value (i=1)
821 for (unsigned int i = 1; i < excessRain.size(); i++)
822 {
823 // "excess" rainfall after Jakeman & Hornberger (1993)
824 // ExcessRain[i] = pcp[i] * WetnessIndex[i];
825
826 // "excess" rainfall after Croke et al. (2004)
827 excessRain[i] = precipitation[i] * ((WetnessIndex[i] + WetnessIndex[i-1]) / 2);
828
829 if (excessRain[i] > precipitation[i])
830 {
831 // if calculated excess rain volume is greater than observed precipitation,
832 // then summarize volume differences and set to current pcp value
833 sum_eRainGTpcp += excessRain[i] - precipitation[i];
834 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
835 // if excess rainfall is greater than observed precip, then
836 // it is reduced to precip!
837 //excessRain[i] = precipitation[i];
838 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
839 }
840
841 if (excessRain[i] < 0.0) excessRain[i] = 0.0;
842
843 // snow module ********************************************
844 if (bSnowModule)
845 {
846 if (temperature[i] < pSnowModule->Get_T_Rain()) excessRain[i] = 0.0;
847 if (temperature[i] > pSnowModule->Get_T_Melt()) excessRain[i] += pSnowModule->Get_MeltRate(i);
848 if ((temperature[i] < pSnowModule->Get_T_Melt()) && (temperature[i] > pSnowModule->Get_T_Rain()))
849 excessRain[i] += pSnowModule->Get_MeltRate(i);
850 }
851 // end snow module ****************************************
852
853 sum += excessRain[i];
854 }
855 sum += excessRain[0]; // add the initial value
856 return sum;
857 }
858 //---------------------------------------------------------------------
859 // Calculate Effective or Excess Rainfall
860 //---------------------------------------------------------------------
CalcExcessRain(double * precipitation,double * temperature,double * WetnessIndex,double * excessRain,double eR_init,double & sum_eRainGTpcp,int size,bool bSnowModule,double T_Rain,double T_Melt,double * MeltRate)861 double Cihacres_eq::CalcExcessRain(double *precipitation, double* temperature, double *WetnessIndex,
862 double *excessRain, double eR_init,
863 double &sum_eRainGTpcp, int size,
864 bool bSnowModule, double T_Rain, double T_Melt, double* MeltRate)
865 {
866 double sum = 0.0; // sum of total ExcessRain of the period
867 sum_eRainGTpcp = 0.0;
868
869 //excessRain[0] = 0.0;
870 excessRain[0] = eR_init;
871 //if (precipitation[0] > 0.0) excessRain[0] = precipitation[0] / 2;
872 // starting at the second value (i=1)
873 for (int i = 1; i < size; i++)
874 {
875 // "excess" rainfall after Jakeman & Hornberger (1993)
876 // ExcessRain[i] = pcp[i] * WetnessIndex[i];
877
878 // "excess" rainfall after Croke et al. (2004)
879 excessRain[i] = precipitation[i] * ((WetnessIndex[i] + WetnessIndex[i-1]) / 2);
880
881 if (excessRain[i] > precipitation[i])
882 {
883 // if calculated excess rain volume is greater than observed precipitation,
884 // then summarize volume differences and set to current pcp value
885 sum_eRainGTpcp += excessRain[i] - precipitation[i];
886 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
887 // if excess rainfall is greater than observed precip, then
888 // it is reduced to precip!
889 //excessRain[i] = precipitation[i];
890 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
891 }
892 if (excessRain[i] < 0.0) excessRain[i] = 0.0;
893
894 // snow module ********************************************
895 if (bSnowModule)
896 {
897 if (temperature[i] < T_Rain) excessRain[i] = 0.0;
898 if (temperature[i] > T_Melt) excessRain[i] += MeltRate[i];
899 if ((temperature[i] < T_Melt) && (temperature[i] > T_Rain))
900 excessRain[i] += MeltRate[i];
901 }
902 // end snow module ****************************************
903
904 sum += excessRain[i];
905 }
906 sum += excessRain[0]; // add the initial value
907 return sum;
908 }
909
910 //---------------------------------------------------------------------
911 // Calculate Effective or Excess Rainfall after Croke et al. (2005)
912 //---------------------------------------------------------------------
CalcExcessRain_Redesign(vector_d & precipitation,vector_d & temperature,vector_d & WetnessIndex,vector_d & excessRain,double eR_init,double & sum_eRainGTpcp,double c,double l,double p,bool bSnowModule,CSnowModule * pSnowMod)913 double Cihacres_eq::CalcExcessRain_Redesign(vector_d &precipitation, vector_d& temperature, vector_d &WetnessIndex,
914 vector_d &excessRain, double eR_init, double &sum_eRainGTpcp,
915 double c, double l, double p, bool bSnowModule, CSnowModule* pSnowMod)
916 {
917 double sum = 0.0; // sum of total ExcessRain of the period
918 sum_eRainGTpcp = 0.0;
919
920 excessRain[0] = eR_init;
921 if (precipitation[0] > 0.0) excessRain[0] = precipitation[0] / 2;
922 // starting at the second value (i=1)
923 for (unsigned int i = 1; i < excessRain.size(); i++)
924 {
925 // "excess" rainfall after Jakeman & Hornberger (1993)
926 // ExcessRain[i] = pcp[i] * WetnessIndex[i];
927
928 // "excess" rainfall after Croke et al. (2004)
929 // excessRain[i] = precipitation[i] * ((WetnessIndex[i] + WetnessIndex[i-1]) / 2);
930
931 // excess rainfall after Croke et al. (2005)
932 // The first equation is as described in the paper,
933 // but Barry Croke told me that parameter "c" must be outside the brackets !!!
934 // excessRain[i] = pow((c *(WetnessIndex[i] - l)),p) * precipitation[i];
935 if ( WetnessIndex[i] - l < 0.0 ) {
936 excessRain[i] = 0.0;
937 } else {
938 excessRain[i] = c * pow((WetnessIndex[i] - l),p) * precipitation[i];
939 }
940
941 if (excessRain[i] > precipitation[i])
942 {
943 // if calculated excess rain volume is greater than observed precipitation,
944 // then summarize volume differences and set to current pcp value
945 sum_eRainGTpcp += excessRain[i] - precipitation[i];
946 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
947 // if excess rainfall is greater than observed precip, then
948 // it is reduced to precip!
949 //excessRain[i] = precipitation[i];
950 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
951 }
952 if (excessRain[i] < 0.0) excessRain[i] = 0.0;
953
954 // snow module ********************************************
955 if (bSnowModule)
956 {
957 if (temperature[i] < pSnowMod->Get_T_Rain()) excessRain[i] = 0.0;
958 if (temperature[i] > pSnowMod->Get_T_Melt()) excessRain[i] += pSnowMod->Get_MeltRate(i);
959 if ((temperature[i] < pSnowMod->Get_T_Melt()) && (temperature[i] > pSnowMod->Get_T_Rain()))
960 excessRain[i] += pSnowMod->Get_MeltRate(i);
961 }
962 // end snow module ****************************************
963
964 sum += excessRain[i];
965 }
966 sum += excessRain[0]; // add the initial value
967 return sum;
968 }
969 //---------------------------------------------------------------------
970 // Calculate Effective or Excess Rainfall after Croke et al. (2005)
971 //---------------------------------------------------------------------
CalcExcessRain_Redesign(double * precipitation,double * temperature,double * WetnessIndex,double * excessRain,double eR_init,double & sum_eRainGTpcp,int size,double c,double l,double p,bool bSnowModule,double T_Rain,double T_Melt,double * MeltRate)972 double Cihacres_eq::CalcExcessRain_Redesign(double *precipitation, double* temperature, double *WetnessIndex,
973 double *excessRain, double eR_init,
974 double &sum_eRainGTpcp,
975 int size, double c, double l, double p,
976 bool bSnowModule, double T_Rain, double T_Melt, double* MeltRate)
977 {
978 double sum = 0.0; // sum of total ExcessRain of the period
979 sum_eRainGTpcp = 0.0;
980
981 //excessRain[0] = 0.0;
982 excessRain[0] = eR_init;
983 //if (precipitation[0] > 0.0) excessRain[0] = precipitation[0] / 2;
984 // starting at the second value (i=1)
985 for (int i = 1; i < size; i++)
986 {
987 // excess rainfall after Croke et al. (2005)
988 // The first equation is as described in the paper,
989 // but Barry Croke told me that parameter "c" must be outside the brackets !!!
990 // excessRain[i] = pow((c *(WetnessIndex[i] - l)),p) * precipitation[i];
991 excessRain[i] = c * pow((WetnessIndex[i] - l),p) * precipitation[i];
992
993 if (excessRain[i] > precipitation[i])
994 {
995 // if calculated excess rain volume is greater than observed precipitation,
996 // then summarize volume differences and set to current pcp value
997 sum_eRainGTpcp += excessRain[i] - precipitation[i];
998 excessRain[i] = precipitation[i];
999 }
1000 if (excessRain[i] < 0.0) excessRain[i] = 0.0;
1001
1002 // snow module ********************************************
1003 if (bSnowModule)
1004 {
1005 if (temperature[i] < T_Rain) excessRain[i] = 0.0;
1006 if (temperature[i] > T_Melt) excessRain[i] += MeltRate[i];
1007 if ((temperature[i] < T_Melt) && (temperature[i] > T_Rain))
1008 excessRain[i] += MeltRate[i];
1009 }
1010 // end snow module ****************************************
1011
1012 sum += excessRain[i];
1013 }
1014 sum += excessRain[0]; // add the initial value
1015 return sum;
1016 }
1017
1018
1019
1020 //---------------------------------------------------------------------
1021 // Assign first and last record (user selection)
1022 //---------------------------------------------------------------------
AssignFirstLastRec(CSG_Table & pTable,int & first,int & last,CSG_String date1,CSG_String date2,int dateField)1023 void Cihacres_eq::AssignFirstLastRec(CSG_Table &pTable, int &first, int &last,
1024 CSG_String date1, CSG_String date2, int dateField)
1025 {
1026 int j;
1027 ///////////////////////////////////////////////////////////////
1028 // searching the first and the last record of the time range
1029 ///////////////////////////////////////////////////////////////
1030
1031 for (j = 0; j < pTable.Get_Record_Count(); j++)
1032 {
1033 if (!date1.Cmp(pTable.Get_Record(j)->asString(dateField)))
1034 {
1035 first = j;
1036 }
1037 else if (!date2.Cmp(pTable.Get_Record(j)->asString(dateField)))
1038 {
1039 last = j;
1040 }
1041 }
1042 }
1043
1044 //---------------------------------------------------------------------
1045
Assign_nStorages(int storconf)1046 int Cihacres_eq::Assign_nStorages(int storconf)
1047 {
1048 switch(storconf)
1049 {
1050 case 0: return(1); break;
1051 case 1: return(2); break;
1052 case 2: return(2); break;
1053 default : return(1);
1054 }
1055 }
1056 //---------------------------------------------------------------------
1057
1058
1059 //---------------------------------------------------------------------
1060 // Summarize all Values in a Vector
1061 //---------------------------------------------------------------------
SumVector(vector_d & input)1062 double Cihacres_eq::SumVector(vector_d &input)
1063 {
1064 double sum = 0;
1065 for (unsigned int j = 0; j < input.size(); j++)
1066 {
1067 sum += input[j];
1068 }
1069 return sum;
1070 }
1071
1072 //---------------------------------------------------------------------
1073 // Summarize all Values in an Array
1074 //---------------------------------------------------------------------
SumVector(double * input,int size)1075 double Cihacres_eq::SumVector(double *input, int size)
1076 {
1077 double sum = 0;
1078 for (int j = 0; j < size; j++)
1079 {
1080 sum += input[j];
1081 }
1082 return sum;
1083 }
1084
1085
1086 //---------------------------------------------------------------------
1087 // Assign temporary Nash-Sutcliffe efficiency according to
1088 // objective function (NSE, NSE_highflow, NSE_lowflow)
1089 //---------------------------------------------------------------------
_Assign_NSE_temp(int obj_func,double NSE,double NSE_highflow,double NSE_lowflow)1090 double Cihacres_eq::_Assign_NSE_temp(int obj_func, double NSE, double NSE_highflow, double NSE_lowflow)
1091 {
1092 switch(obj_func)
1093 {
1094 case 0: // NSE
1095 return(NSE);
1096 case 1: // NSE high flow
1097 return(NSE_highflow);
1098 case 2: // NSE low flow (not yet implemented!)
1099 return(NSE_lowflow);
1100 break;
1101 default: return(NSE);
1102 } // end switch(m_obj_func)
1103 }
1104 //---------------------------------------------------------------------
1105
1106
1107
1108 ///////////////////////////////////////////////////////////////////////
1109 //
1110 // GET FUNCTIONS
1111 //
1112 ///////////////////////////////////////////////////////////////////////
get_vq()1113 double Cihacres_eq::get_vq()
1114 {
1115 // return (bq / (bq + bs));
1116 return(vq);
1117 }
1118 //---------------------------------------------------------------------
get_vs()1119 double Cihacres_eq::get_vs()
1120 {
1121 // return (1 - get_vq());
1122 return(vs);
1123 }
1124 //---------------------------------------------------------------------
1125 //---------------------------------------------------------------------
get_sum_streamflowMM_obs(int size)1126 double Cihacres_eq::get_sum_streamflowMM_obs(int size)
1127 {
1128 double sum_streamflowMM_obs = 0.0;
1129 // calculate sum of observed streamflow in [mm]
1130 for (int j = 0; j < size; j++)
1131 {
1132 sum_streamflowMM_obs += streamflowMM_obs[j];
1133 }
1134 return(sum_streamflowMM_obs);
1135 }
1136 //---------------------------------------------------------------------
get_sum_precipitation(int size)1137 double Cihacres_eq::get_sum_precipitation(int size)
1138 {
1139 double sum_pcp = 0.0;
1140 for (int i = 0; i < size; i++)
1141 {
1142 sum_pcp += precipitation[i];
1143 }
1144 return(sum_pcp);
1145 }
1146 //---------------------------------------------------------------------
get_NSE()1147 double Cihacres_eq::get_NSE()
1148 {
1149 return(NSE);
1150 }
1151 //---------------------------------------------------------------------
get_streamflow_sim()1152 vector_d Cihacres_eq::get_streamflow_sim()
1153 {
1154 return(streamflow_sim);
1155 }
1156 //---------------------------------------------------------------------
get_excessRain()1157 vector_d Cihacres_eq::get_excessRain()
1158 {
1159 return(excessRain);
1160 }
1161 //---------------------------------------------------------------------
get_WetnessIndex()1162 vector_d Cihacres_eq::get_WetnessIndex()
1163 {
1164 return(WetnessIndex);
1165 }
1166 //---------------------------------------------------------------------
get_Tw()1167 vector_d Cihacres_eq::get_Tw()
1168 {
1169 return(Tw);
1170 }
1171 //---------------------------------------------------------------------
1172
1173 //---------------------------------------------------------------------
1174
1175 ///////////////////////////////////////////////////////////////////////
1176 //
1177 // PRIVATE MEMBER FUNCTIONS
1178 //
1179 ///////////////////////////////////////////////////////////////////////
1180
1181 //---------------------------------------------------------------------
1182 // Initialize (Resize) global Vectors
1183 //---------------------------------------------------------------------
_InitVectorsStart(int size)1184 void Cihacres_eq::_InitVectorsStart(int size)
1185 {
1186 streamflow_sim.resize(size);
1187 excessRain.resize(size);
1188 WetnessIndex.resize(size);
1189 Tw.resize(size);
1190 streamflowMM_obs.resize(size);
1191 //if (IHAC_version == 2)
1192 //{
1193 // m_pSnowMod->SnowStorage = new double[size];
1194 // m_pSnowMod->MeltRate = new double[size];
1195 // for (int i = 0; i < size; i++)
1196 // {
1197 // m_pSnowMod->SnowStorage[i] = 0.0;
1198 // m_pSnowMod->MeltRate[i] = 0.0;
1199 // }
1200 // //SnowModule(m_pSnowMod->SnowStorage,m_pSnowMod->MeltRate,temperature,precipitation, m_pSnowMod->T_Rain,m_pSnowMod->T_Melt,m_pSnowMod->DD_FAC, size);
1201 //}
1202 }
1203
1204 //---------------------------------------------------------------------
1205 // Zero all Vectors
1206 //---------------------------------------------------------------------
_ZeroAllVectors()1207 void Cihacres_eq::_ZeroAllVectors()
1208 {
1209 // Resize all Vectors to Zero
1210 streamflow_sim.resize(0);
1211 excessRain.resize(0);
1212 WetnessIndex.resize(0);
1213 Tw.resize(0);
1214 date.resize(0);
1215 streamflow_obs.resize(0);
1216 precipitation.resize(0);
1217 temperature.resize(0);
1218 }
1219