1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // sim_hydrology //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // overland_flow.cpp //
14 // //
15 // Copyright (C) 2020 by //
16 // Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "overland_flow.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
COverland_Flow(void)59 COverland_Flow::COverland_Flow(void)
60 {
61 Set_Name (_TL("Overland Flow"));
62
63 Set_Author ("O.Conrad (c) 2020");
64
65 Set_Description (_TW(
66 "A simple overland flow simulation."
67 ));
68
69 //-----------------------------------------------------
70 Parameters.Add_Grid("",
71 "DEM" , _TL("Elevation"),
72 _TL(""),
73 PARAMETER_INPUT
74 );
75
76 Parameters.Add_Grid_or_Const("",
77 "ROUGHNESS" , _TL("Roughness"),
78 _TL("Ks Strickler = 1/n Manning"),
79 20., 0., true
80 );
81
82 Parameters.Add_Choice("ROUGHNESS",
83 "STRICKLER" , _TL("Type"),
84 _TL("Ks Strickler = 1/n Gauckler-Manning"),
85 CSG_String::Format("%s|%s",
86 SG_T("Strickler Ks, [m^1/3 / s]"),
87 SG_T("Gauckler-Manning n, [s / m^1/3]")
88 ), 0
89 );
90
91 //-----------------------------------------------------
92 Parameters.Add_Grid_or_Const("",
93 "INTER_MAX" , _TL("Interception Capacity [mm]"),
94 _TL(""),
95 0., 0., true
96 );
97
98 Parameters.Add_Grid_or_Const("",
99 "POND_MAX" , _TL("Ponding Capacity [mm]"),
100 _TL(""),
101 0., 0., true
102 );
103
104 Parameters.Add_Grid_or_Const("",
105 "INFIL_MAX" , _TL("Infiltration Capacity [mm/h]"),
106 _TL(""),
107 0., 0., true
108 );
109
110 Parameters.Add_Grid("",
111 "INTERCEPT" , _TL("Interception [mm]"),
112 _TL(""),
113 PARAMETER_OUTPUT
114 );
115
116 Parameters.Add_Grid("",
117 "PONDING" , _TL("Ponding [mm]"),
118 _TL(""),
119 PARAMETER_OUTPUT
120 );
121
122 Parameters.Add_Grid("",
123 "INFILTRAT" , _TL("Infiltration [mm]"),
124 _TL(""),
125 PARAMETER_OUTPUT
126 );
127
128 Parameters.Add_Grid("",
129 "FLOW" , _TL("Flow [mm]"),
130 _TL(""),
131 PARAMETER_OUTPUT
132 );
133
134 Parameters.Add_Grid("",
135 "VELOCITY" , _TL("Velocity [m/h]"),
136 _TL(""),
137 PARAMETER_OUTPUT_OPTIONAL
138 );
139
140 //-----------------------------------------------------
141 Parameters.Add_Grid_System("",
142 "WEATHER" , _TL("Weather Grid System"),
143 _TL("")
144 );
145
146 Parameters.Add_Grid_or_Const("WEATHER",
147 "PRECIP" , _TL("Precipitation [mm/h]"),
148 _TL(""),
149 0., 0., true
150 );
151
152 Parameters.Add_Grid_or_Const("WEATHER",
153 "ET_POT" , _TL("Potential Evapotranspiration [mm/h]"),
154 _TL(""),
155 0., 0., true
156 );
157
158 Parameters.Add_Bool("",
159 "RESET" , _TL("Reset"),
160 _TL("If checked storages (flow, ponding, interception) and sinks (infiltration) will be set to zero."),
161 true
162 );
163
164 Parameters.Add_Double("",
165 "TIME_STOP" , _TL("Simulation Time [h]"),
166 _TL("Simulation time in hours."),
167 6., 0., true
168 );
169
170 Parameters.Add_Double("",
171 "TIME_STEP" , _TL("Time Step Adjustment"),
172 _TL("Choosing a lower value will result in a better numerical precision but also in a longer calculation time."),
173 0.5, 0.01, true, 1., true
174 );
175
176 // Parameters.Add_Double("",
177 // "V_MIN" , _TL("Minimum Velocity [m/h]"),
178 // _TL(""),
179 // 0., 0., true
180 // );
181
182 Parameters.Add_Bool("",
183 "FLOW_OUT" , _TL("Overland Flow Summary"),
184 _TL("Report the amount of overland flow that left the covered area."),
185 false
186 );
187
188 Parameters.Add_Table("",
189 "SUMMARY" , _TL("Overland Flow Summary"),
190 _TL(""),
191 PARAMETER_OUTPUT_OPTIONAL
192 );
193
194 Parameters.Add_Double("",
195 "TIME_UPDATE", _TL("Map Update Frequency [Minutes]"),
196 _TL("Map update frequency in minutes. Set to zero to update each simulation time step."),
197 1., 0., true
198 );
199
200 Parameters.Add_Bool("TIME_UPDATE", "UPDATE_FLOW_FIXED" , _TL("Fixed Color Stretch for Flow" ), _TL(""), false);
201 Parameters.Add_Range( "UPDATE_FLOW_FIXED", "UPDATE_FLOW_RANGE", _TL("Fixed Color Stretch" ), _TL(""), 0., 1., 0., true);
202
203 Parameters.Add_Bool("TIME_UPDATE", "UPDATE_VELO_FIXED" , _TL("Fixed Color Stretch for Velocity"), _TL(""), true);
204 Parameters.Add_Range( "UPDATE_VELO_FIXED", "UPDATE_VELO_RANGE", _TL("Fixed Color Stretch" ), _TL(""), 750., 1., 0., true);
205 }
206
207
208 ///////////////////////////////////////////////////////////
209 // //
210 ///////////////////////////////////////////////////////////
211
212 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)213 int COverland_Flow::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
214 {
215 if( pParameter->Cmp_Identifier("STRICKLER") && (*pParameters)("ROUGHNESS")->asDouble() > 0. )
216 {
217 pParameters->Set_Parameter("ROUGHNESS_DEFAULT", 1. / (*pParameters)("ROUGHNESS")->asDouble());
218 }
219
220 return( CSG_Tool_Grid::On_Parameter_Changed(pParameters, pParameter) );
221 }
222
223 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)224 int COverland_Flow::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
225 {
226 if( pParameter->Cmp_Identifier("INTER_MAX") )
227 {
228 pParameters->Set_Enabled("INTERCEPT", pParameter->asGrid() || pParameter->asDouble() > 0.);
229 }
230
231 if( pParameter->Cmp_Identifier("POND_MAX" ) )
232 {
233 pParameters->Set_Enabled("PONDING" , pParameter->asGrid() || pParameter->asDouble() > 0.);
234 }
235
236 if( pParameter->Cmp_Identifier("INFIL_MAX") )
237 {
238 pParameters->Set_Enabled("INFILTRAT", pParameter->asGrid() || pParameter->asDouble() > 0.);
239 }
240
241 if( pParameter->Cmp_Identifier("UPDATE_FLOW_FIXED") )
242 {
243 pParameters->Set_Enabled("UPDATE_FLOW_RANGE", pParameter->asBool());
244 }
245
246 if( pParameter->Cmp_Identifier("VELOCITY") )
247 {
248 pParameters->Set_Enabled("UPDATE_VELO_FIXED", pParameter->asDataObject() != NULL);
249 }
250
251 if( pParameter->Cmp_Identifier("UPDATE_VELO_FIXED") )
252 {
253 pParameters->Set_Enabled("UPDATE_VELO_RANGE", pParameter->asBool());
254 }
255
256 if( pParameter->Cmp_Identifier("FLOW_OUT") )
257 {
258 pParameters->Set_Enabled("SUMMARY", pParameter->asBool());
259 }
260
261 return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
262 }
263
264
265 ///////////////////////////////////////////////////////////
266 // //
267 ///////////////////////////////////////////////////////////
268
269 //---------------------------------------------------------
On_Execute(void)270 bool COverland_Flow::On_Execute(void)
271 {
272 if( !Initialize() )
273 {
274 Finalize();
275
276 return( false );
277 }
278
279 //-----------------------------------------------------
280 double Update_Last = 0., Update = Parameters("TIME_UPDATE")->asDouble() / 60.; // from minutes to hours
281
282 double Time, Time_Stop = Parameters("TIME_STOP")->asDouble();
283
284 for(Time=0.; Time<=Time_Stop && Set_Time_Stamp(Time); Time+=m_dTime)
285 {
286 SG_UI_ProgressAndMsg_Lock(true);
287
288 Do_Time_Step();
289
290 if( Time >= Update_Last )
291 {
292 if( Update > 0. )
293 {
294 Update_Last = Update * (1. + floor(Time / Update));
295 }
296
297 Do_Updates();
298 }
299
300 SG_UI_ProgressAndMsg_Lock(false);
301 }
302
303 //-----------------------------------------------------
304 double s = Time;
305 int h = (int)s; s = 60. * (s - h);
306 int m = (int)s; s = 60. * (s - m);
307
308 Message_Fmt("\n____\n%s: %02dh %02dm %02fs (= %g %s)\n", _TL("Simulation Time"), h, m, s, Time, _TL("hours"));
309
310 //-----------------------------------------------------
311 Finalize();
312
313 return( true );
314 }
315
316
317 ///////////////////////////////////////////////////////////
318 // //
319 ///////////////////////////////////////////////////////////
320
321 //---------------------------------------------------------
Initialize(void)322 bool COverland_Flow::Initialize(void)
323 {
324 m_pDEM = Parameters("DEM" )->asGrid ();
325
326 m_pRoughness = Parameters("ROUGHNESS")->asGrid ();
327 m_Roughness = Parameters("ROUGHNESS")->asDouble();
328
329 m_pPrecipitation = Parameters("PRECIP" )->asGrid ();
330 m_Precipitation = Parameters("PRECIP" )->asDouble();
331
332 m_pETpot = Parameters("ET_POT" )->asGrid ();
333 m_ETpot = Parameters("ET_POT" )->asDouble();
334
335 m_pIntercept_max = Parameters("INTER_MAX")->asGrid ();
336 m_Intercept_max = Parameters("INTER_MAX")->asDouble();
337 m_pIntercept = m_pIntercept_max || m_Intercept_max > 0.
338 ? Parameters("INTERCEPT")->asGrid () : NULL;
339
340 m_pPonding_max = Parameters("POND_MAX" )->asGrid ();
341 m_Ponding_max = Parameters("POND_MAX" )->asDouble();
342 m_pPonding = m_pPonding_max || m_Ponding_max > 0.
343 ? Parameters("PONDING" )->asGrid () : NULL;
344
345 m_pInfiltrat_max = Parameters("INFIL_MAX")->asGrid ();
346 m_Infiltrat_max = Parameters("INFIL_MAX")->asDouble();
347 m_pInfiltrat = m_pInfiltrat_max || m_Infiltrat_max > 0.
348 ? Parameters("INFILTRAT")->asGrid () : NULL;
349
350 m_pFlow = Parameters("FLOW" )->asGrid ();
351 m_pVelocity = Parameters("VELOCITY" )->asGrid ();
352
353 m_bStrickler = Parameters("STRICKLER")->asInt() == 0;
354
355 m_vMin = 0.; // = Parameters("V_MIN" )->asDouble();
356
357 m_bFlow_Out = Parameters("FLOW_OUT" )->asBool ();
358 m_Flow_Out = 0.;
359
360 //-----------------------------------------------------
361 if( Parameters("RESET")->asBool() )
362 {
363 #pragma omp parallel
364 for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
365 {
366 if( m_pDEM->is_NoData(x, y) )
367 {
368 if( m_pIntercept ) m_pIntercept->Set_NoData(x, y);
369 if( m_pPonding ) m_pPonding ->Set_NoData(x, y);
370 if( m_pInfiltrat ) m_pInfiltrat->Set_NoData(x, y);
371 if( m_pFlow ) m_pFlow ->Set_NoData(x, y);
372 }
373 else
374 {
375 if( m_pIntercept ) m_pIntercept->Set_Value(x, y, 0.);
376 if( m_pPonding ) m_pPonding ->Set_Value(x, y, 0.);
377 if( m_pInfiltrat ) m_pInfiltrat->Set_Value(x, y, 0.);
378 if( m_pFlow ) m_pFlow ->Set_Value(x, y, 0.);
379 }
380 }
381
382 CSG_Colors Colors(5, SG_COLORS_WHITE_BLUE); Colors.Set_Color(0, 240, 240, 240);
383
384 DataObject_Set_Colors(m_pIntercept, Colors);
385 DataObject_Set_Colors(m_pPonding , Colors);
386 DataObject_Set_Colors(m_pInfiltrat, Colors);
387 DataObject_Set_Colors(m_pFlow , Colors);
388 }
389
390 if( m_pVelocity )
391 {
392 m_pVelocity->Set_NoData_Value(0.);
393
394 CSG_Colors Colors(11, SG_COLORS_RAINBOW); Colors.Set_Color(0, 255, 255, 255);
395
396 DataObject_Set_Colors(m_pVelocity , Colors);
397 }
398
399 DataObject_Update(m_pFlow, SG_UI_DATAOBJECT_SHOW); // show in new map
400
401 //-----------------------------------------------------
402 m_Flow.Create(Get_System() , SG_DATATYPE_Float);
403 m_v .Create(Get_System(), 9, 0., SG_DATATYPE_Float);
404
405 return( true );
406 }
407
408 //---------------------------------------------------------
Finalize(void)409 bool COverland_Flow::Finalize(void)
410 {
411 m_Flow.Destroy();
412 m_v .Destroy();
413
414 if( !Process_Get_Okay() )
415 {
416 SG_UI_Process_Set_Okay();
417
418 if( m_pIntercept ) m_pIntercept->Update(true);
419 if( m_pPonding ) m_pPonding ->Update(true);
420 if( m_pInfiltrat ) m_pInfiltrat->Update(true);
421 if( m_pFlow ) m_pFlow ->Update(true);
422 }
423
424 if( m_bFlow_Out )
425 {
426 double Flow_Sum = 0., Inf_Sum = 0.;
427
428 for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
429 {
430 if( !m_pDEM->is_NoData(x, y) )
431 {
432 Flow_Sum += m_pFlow->asDouble(x, y);
433
434 if( m_pInfiltrat )
435 {
436 Inf_Sum += m_pInfiltrat->asDouble(x, y);
437 }
438 }
439 }
440
441 Message_Fmt("\n____\n%s\n%s:\t%g\n%s:\t%g\n%s:\t%g\n%s:\t%g\n", _TL("Overland Flow Summary"),
442 _TL("flow in area"), Flow_Sum,
443 _TL("flow out" ), m_Flow_Out,
444 _TL("infiltration"), Inf_Sum,
445 _TL("total" ), m_Flow_Out + Flow_Sum + Inf_Sum
446 );
447
448 CSG_Table *pSummary = Parameters("SUMMARY")->asTable();
449
450 if( pSummary )
451 {
452 pSummary->Destroy();
453 pSummary->Set_Name(_TL("Overland Flow Summary"));
454 pSummary->Add_Field(_TL("Parameter"), SG_DATATYPE_String);
455 pSummary->Add_Field(_TL("Value" ), SG_DATATYPE_Double);
456
457 #define Add_Entry(name, value) { CSG_Table_Record &r = *pSummary->Add_Record(); r.Set_Value(0, name); r.Set_Value(1, value); }
458
459 Add_Entry(_TL("flow in area"), Flow_Sum);
460 Add_Entry(_TL("flow out" ), m_Flow_Out);
461 Add_Entry(_TL("infiltration"), Inf_Sum);
462 Add_Entry(_TL("total" ), m_Flow_Out + Flow_Sum + Inf_Sum);
463 }
464 }
465
466 return( true );
467 }
468
469
470 ///////////////////////////////////////////////////////////
471 // //
472 ///////////////////////////////////////////////////////////
473
474 //---------------------------------------------------------
Do_Updates(void)475 bool COverland_Flow::Do_Updates(void)
476 {
477 DataObject_Update(m_pIntercept);
478 DataObject_Update(m_pInfiltrat);
479
480 if( Parameters("UPDATE_FLOW_FIXED")->asBool() == false )
481 {
482 DataObject_Update(m_pFlow);
483 }
484 else
485 {
486 DataObject_Update(m_pFlow,
487 Parameters("UPDATE_FLOW_RANGE.MIN")->asDouble(),
488 Parameters("UPDATE_FLOW_RANGE.MAX")->asDouble()
489 );
490 }
491
492 if( Parameters("UPDATE_VELO_FIXED")->asBool() == false )
493 {
494 DataObject_Update(m_pVelocity);
495 }
496 else
497 {
498 DataObject_Update(m_pVelocity,
499 Parameters("UPDATE_VELO_RANGE.MIN")->asDouble(),
500 Parameters("UPDATE_VELO_RANGE.MAX")->asDouble()
501 );
502 }
503
504 return( true );
505 }
506
507 //---------------------------------------------------------
Set_Time_Stamp(double Time)508 bool COverland_Flow::Set_Time_Stamp(double Time)
509 {
510 double s = Time;
511 int h = (int)s; s = 60. * (s - h);
512 int m = (int)s; s = 60. * (s - m);
513
514 Process_Set_Text(CSG_String::Format("%s: %02d:%02d:%02d %s: %.2f [min]",
515 _TL("Time"), h, m, (int)s,
516 _TL("Step"), 60. * m_dTime
517 ));
518
519 return( Process_Get_Okay() );
520 }
521
522
523 ///////////////////////////////////////////////////////////
524 // //
525 ///////////////////////////////////////////////////////////
526
527 //---------------------------------------------------------
Do_Time_Step(void)528 bool COverland_Flow::Do_Time_Step(void)
529 {
530 m_vMax = 0.;
531
532 #pragma omp parallel for
533 for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
534 {
535 Get_Velocity(x, y);
536 }
537
538 //-----------------------------------------------------
539 if( m_vMax > 0. )
540 {
541 m_dTime = Parameters("TIME_STEP")->asDouble() * Get_Cellsize() / m_vMax; // Courant�Friedrichs�Lewy (CFL) condition
542
543 #pragma omp parallel for
544 for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
545 {
546 Set_Flow_Lateral(x, y);
547 }
548 }
549 else
550 {
551 m_dTime = 1. / 60.; // 1 min
552 }
553
554 //-----------------------------------------------------
555 #pragma omp parallel for
556 for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
557 {
558 Set_Flow_Vertical(x, y);
559 }
560
561 //-----------------------------------------------------
562 return( true );
563 }
564
565
566 ///////////////////////////////////////////////////////////
567 // //
568 ///////////////////////////////////////////////////////////
569
570 //---------------------------------------------------------
571 #define GET_GRID_OR_CONST(g, c) { double v; if( !g || !g->Get_Value(Get_System().Get_Grid_to_World(x, y), v) ) { v = c; } return( v > 0. ? m_dTime * v : 0. ); }
572
573 //---------------------------------------------------------
Get_Precipitation(int x,int y)574 inline double COverland_Flow::Get_Precipitation(int x, int y)
575 {
576 GET_GRID_OR_CONST(m_pPrecipitation, m_Precipitation);
577 }
578
579 //---------------------------------------------------------
Get_ETpot(int x,int y)580 inline double COverland_Flow::Get_ETpot(int x, int y)
581 {
582 GET_GRID_OR_CONST(m_pETpot, m_ETpot);
583 }
584
585 //---------------------------------------------------------
586 #undef GET_GRID_OR_CONST
587
588
589 ///////////////////////////////////////////////////////////
590 // //
591 ///////////////////////////////////////////////////////////
592
593 //---------------------------------------------------------
594 #define GET_GRID_OR_CONST(g, c) (g && !g->is_NoData(x, y) ? g->asDouble(x, y) : c)
595
596 //---------------------------------------------------------
Get_Roughness(int x,int y)597 inline double COverland_Flow::Get_Roughness(int x, int y)
598 {
599 double Value = GET_GRID_OR_CONST(m_pRoughness, m_Roughness);
600
601 return( Value > 0. ? (m_bStrickler ? Value : 1. / Value) : 0. );
602 }
603
604 //---------------------------------------------------------
Get_Intercept_max(int x,int y)605 inline double COverland_Flow::Get_Intercept_max(int x, int y)
606 {
607 double Value = GET_GRID_OR_CONST(m_pIntercept_max, m_Intercept_max);
608
609 return( Value > 0. ? Value : 0. );
610 }
611
612 //---------------------------------------------------------
Get_Ponding(int x,int y)613 inline double COverland_Flow::Get_Ponding(int x, int y)
614 {
615 double Value = GET_GRID_OR_CONST(m_pPonding_max, m_Ponding_max);
616
617 return( Value > 0. ? Value : 0. );
618 }
619
620 //---------------------------------------------------------
Get_Infiltration(int x,int y)621 inline double COverland_Flow::Get_Infiltration(int x, int y)
622 {
623 double Value = GET_GRID_OR_CONST(m_pInfiltrat_max, m_Infiltrat_max);
624
625 return( Value > 0. ? m_dTime * Value : 0. );
626 }
627
628 //---------------------------------------------------------
629 #undef GET_GRID_OR_CONST
630
631
632 ///////////////////////////////////////////////////////////
633 // //
634 ///////////////////////////////////////////////////////////
635
636 //---------------------------------------------------------
Get_Surface(int x,int y)637 inline double COverland_Flow::Get_Surface(int x, int y)
638 {
639 return( m_pDEM->asDouble(x, y) + m_pFlow->asDouble(x, y) / 1000. );
640 }
641
642
643 ///////////////////////////////////////////////////////////
644 // //
645 ///////////////////////////////////////////////////////////
646
647 //---------------------------------------------------------
Get_Neighbour(int x,int y,int i,int & ix,int & iy)648 inline bool COverland_Flow::Get_Neighbour(int x, int y, int i, int &ix, int &iy)
649 {
650 return( m_pDEM->is_InGrid(ix = Get_xTo(i, x), iy = Get_yTo(i, y)) );
651 }
652
653 //---------------------------------------------------------
Get_Slope(int x,int y,int i)654 inline double COverland_Flow::Get_Slope(int x, int y, int i)
655 {
656 int ix, iy; double dz;
657
658 if ( Get_Neighbour(x, y, i , ix, iy) )
659 {
660 dz = Get_Surface( x, y) - Get_Surface(ix, iy);
661 }
662 else if( Get_Neighbour(x, y, i + 4, ix, iy) )
663 {
664 dz = Get_Surface(ix, iy) - Get_Surface( x, y);
665 }
666 else
667 {
668 return( 0. );
669 }
670
671 return( dz > 0. ? dz / Get_Length(i) : 0. ); // the tangens of slope
672 }
673
674
675 ///////////////////////////////////////////////////////////
676 // //
677 ///////////////////////////////////////////////////////////
678
679 //---------------------------------------------------------
Get_Velocity(double Flow,double Slope,double Roughness)680 inline double COverland_Flow::Get_Velocity(double Flow, double Slope, double Roughness)
681 {
682 // V = Ks * T^2/3 * S^1/2
683 // V : velocity [m/s]
684 // T : water depth [m]
685 // S : tangens of slope
686 // Ks: Strickler roughness factor (= 1/n, n: Manning roughness)
687
688 return( 3600 * Roughness * pow(Flow / 1000., 2. / 3.) * sqrt(Slope) ); // 3600 [m/s] => 1 [m/h]
689 }
690
691 //---------------------------------------------------------
Get_Velocity(int x,int y)692 bool COverland_Flow::Get_Velocity(int x, int y)
693 {
694 if( m_pDEM->is_NoData(x, y) )
695 {
696 return( false );
697 }
698
699 double Flow = m_pFlow->asDouble(x, y);
700
701 if( Flow > 0. )
702 {
703 double vMax = 0., vSum = 0.;
704
705 for(int i=0; i<8; i++)
706 {
707 double Slope = Get_Slope(x, y, i);
708
709 if( Slope > 0. )
710 {
711 double v = Get_Velocity(Flow, Slope, Get_Roughness(x, y)); if( v < m_vMin ) { v = m_vMin; }
712
713 if( vMax < v )
714 {
715 vMax = v;
716 }
717
718 vSum += v;
719
720 m_v[i].Set_Value(x, y, v);
721 }
722 else
723 {
724 m_v[i].Set_Value(x, y, 0.);
725 }
726 }
727
728 if( m_vMax < vMax )
729 {
730 #pragma omp critical
731 {
732 if( m_vMax < vMax ) // could have been changed by another thread after the comparison outside the critical section
733 {
734 m_vMax = vMax;
735 }
736 }
737 }
738
739 //-------------------------------------------------
740 m_v[8].Set_Value(x, y, vSum);
741
742 if( m_pVelocity )
743 {
744 m_pVelocity->Set_Value(x, y, vMax); // 1/3600 * [m/h] => [m/s]
745 }
746 }
747 else
748 {
749 if( m_pVelocity )
750 {
751 m_pVelocity->Set_Value(x, y, 0.);
752 }
753 }
754
755 return( true );
756 }
757
758
759 ///////////////////////////////////////////////////////////
760 // //
761 ///////////////////////////////////////////////////////////
762
763 //---------------------------------------------------------
Get_Flow_Lateral(int x,int y,int i,bool bInverse)764 inline double COverland_Flow::Get_Flow_Lateral(int x, int y, int i, bool bInverse)
765 {
766 if( bInverse )
767 {
768 if( !m_pDEM->is_InGrid(x = Get_xTo(i, x), y = Get_yTo(i, y)) )
769 {
770 return( 0. );
771 }
772
773 i = (i + 4) % 8;
774 }
775
776 double Flow, v;
777
778 if( (Flow = m_pFlow->asDouble(x, y)) > 0. && (v = m_v[i].asDouble(x, y)) > 0. )
779 {
780 Flow = Flow * v / m_v[8].asDouble(x, y) * m_dTime * v / Get_Length(i);
781
782 if( m_bFlow_Out && !bInverse && !is_InGrid(Get_xTo(i, x), Get_yTo(i, y)) )
783 {
784 #pragma omp atomic
785 m_Flow_Out += Flow;
786 }
787
788 return( Flow );
789 }
790
791 return( 0. );
792 }
793
794 //---------------------------------------------------------
Set_Flow_Lateral(int x,int y)795 bool COverland_Flow::Set_Flow_Lateral(int x, int y)
796 {
797 double iFlow, Flow = m_pFlow->asDouble(x, y);
798
799 for(int i=0; i<8; i++)
800 {
801 if ( (iFlow = Get_Flow_Lateral(x, y, i, false)) > 0. ) // downslope flow leaving cell
802 {
803 Flow -= iFlow;
804 }
805 else if( (iFlow = Get_Flow_Lateral(x, y, i, true)) > 0. ) // upslope flow entering cell
806 {
807 Flow += iFlow;
808 }
809 }
810
811 m_Flow.Set_Value(x, y, Flow > 0. ? Flow : 0.);
812
813 return( true );
814 }
815
816
817 ///////////////////////////////////////////////////////////
818 // //
819 ///////////////////////////////////////////////////////////
820
821 //---------------------------------------------------------
Set_Flow_Vertical(int x,int y)822 bool COverland_Flow::Set_Flow_Vertical(int x, int y)
823 {
824 if( m_pDEM->is_NoData(x, y) )
825 {
826 return( false );
827 }
828
829 //-----------------------------------------------------
830 double P = Get_Precipitation(x, y);
831
832 double I = m_pIntercept ? m_pIntercept->asDouble(x, y) : 0.;
833 double Imax = Get_Intercept_max(x, y);
834
835 if( Imax > 0. )
836 {
837 double Deficit = Imax - I;
838
839 if( Deficit > P )
840 {
841 I += P;
842 P = 0.;
843 }
844 else
845 {
846 I += Deficit;
847 P -= Deficit;
848 }
849 }
850
851 double Q = P + m_Flow.asDouble(x, y) + (m_pPonding ? m_pPonding->asDouble(x, y) : 0.);
852
853 //-----------------------------------------------------
854 if( Q > 0. )
855 {
856 double ETpot = Get_ETpot(x, y);
857
858 if( ETpot >= I + Q )
859 {
860 I = 0.;
861 Q = 0.;
862 }
863 else
864 {
865 double dIQ = I / Q;
866 double sIQ = (I + Q) - ETpot;
867
868 I = sIQ * ( dIQ);
869 Q = sIQ * (1. - dIQ);
870 }
871 }
872 else if( I > .0 )
873 {
874 double ETpot = Get_ETpot(x, y);
875
876 I = I > ETpot ? I - ETpot : 0.;
877 }
878
879 if( m_pIntercept )
880 {
881 m_pIntercept->Set_Value(x, y, I);
882 }
883
884 //-----------------------------------------------------
885 if( Q > 0. )
886 {
887 double Infiltration = Get_Infiltration(x, y);
888
889 if( Infiltration >= Q )
890 {
891 m_pInfiltrat->Add_Value(x, y, Q);
892
893 Q = 0.;
894 }
895 else if( Infiltration > 0. )
896 {
897 m_pInfiltrat->Add_Value(x, y, Infiltration);
898
899 Q -= Infiltration;
900 }
901 }
902
903 //-----------------------------------------------------
904 if( Q > 0. )
905 {
906 double Ponding_max = Get_Ponding(x, y);
907
908 if( Ponding_max >= Q )
909 {
910 m_pPonding->Set_Value(x, y, Q);
911
912 Q = 0.;
913 }
914 else if( Ponding_max > 0. )
915 {
916 m_pPonding->Set_Value(x, y, Ponding_max);
917
918 Q -= Ponding_max;
919 }
920 }
921
922 //-----------------------------------------------------
923 m_pFlow->Set_Value(x, y, Q);
924
925 return( true );
926 }
927
928
929 ///////////////////////////////////////////////////////////
930 // //
931 // //
932 // //
933 ///////////////////////////////////////////////////////////
934
935 //---------------------------------------------------------
936