Line data Source code
1 : !> \file mo_mrm_routing.f90
2 : !> \brief \copybrief mo_mrm_routing
3 : !> \details \copydetails mo_mrm_routing
4 :
5 : !> \brief Performs runoff routing for mHM at level L11.
6 : !> \details This module performs flood routing at a given time step through the stream network at level L11 to the sink cell.
7 : !! The Muskingum flood routing algorithm is used.
8 : !> \changelog
9 : !! - Stephan Thober Aug 2015
10 : !! - adapted to mRM
11 : !! - Sebastian Mueller Jun 2020
12 : !! - outsourcing helper functions
13 : !> \authors Luis Samaniego
14 : !> \date Dec 2012
15 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
16 : !! mHM is released under the LGPLv3+ license \license_note
17 : !> \ingroup f_mrm
18 : MODULE mo_mrm_routing
19 :
20 : ! This module performs runoff flood routing for mHM.
21 :
22 : ! Written Luis Samaniego, Dec 2012
23 :
24 : USE mo_kind, ONLY : i4, dp
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : PUBLIC :: mRM_routing
31 :
32 : ! ------------------------------------------------------------------
33 :
34 : CONTAINS
35 :
36 : ! ------------------------------------------------------------------
37 :
38 : ! NAME
39 : ! mRM_routing
40 :
41 : ! PURPOSE
42 : !> \brief route water given runoff
43 :
44 : !> \details This routine first performs mpr for the routing variables
45 : !> if required, then accumulates the runoff to the routing resolution
46 : !> and eventually routes the water in a third step. The last step is
47 : !> repeated multiple times if the routing timestep is smaller than
48 : !> the timestep of the hydrological timestep
49 :
50 : ! INTENT(IN)
51 : !> \param[in] "logical :: read_states" whether states are derived from restart file
52 : !> \param[in] "integer(i4) :: processCase" Process switch for routing
53 : !> \param[in] "real(dp), dimension(:) :: global_routing_param" routing parameters
54 : !> \param[in] "real(dp), dimension(:) :: L1_total_runoff" total runoff from L1 grid cells
55 : !> \param[in] "real(dp), dimension(:) :: L1_areaCell" L1 cell area
56 : !> \param[in] "integer(i4), dimension(:) :: L1_L11_Id" L1 cell ids on L11
57 : !> \param[in] "real(dp), dimension(:) :: L11_areaCell" L11 cell area
58 : !> \param[in] "integer(i4), dimension(:) :: L11_L1_Id" L11 cell ids on L1
59 : !> \param[in] "integer(i4), dimension(:) :: L11_netPerm" L11 routing order
60 : !> \param[in] "integer(i4), dimension(:) :: L11_fromN" L11 source grid cell order
61 : !> \param[in] "integer(i4), dimension(:) :: L11_toN" L11 target grid cell order
62 : !> \param[in] "integer(i4) :: L11_nOutlets" L11 number of outlets/sinks
63 : !> \param[in] "integer(i4) :: timestep" simulation timestep in [h]
64 : !> \param[in] "real(dp) :: tsRoutFactor" factor between routing timestep and
65 : !> hydrological timestep
66 : !> \param[in] "integer(i4) :: nNodes" number of nodes
67 : !> \param[in] "integer(i4) :: nInflowGauges" number of inflow gauges
68 : !> \param[in] "integer(i4), dimension(:) :: InflowGaugeIndexList" index list of inflow gauges
69 : !> \param[in] "logical, dimension(:) :: InflowGaugeHeadwater" flag for headwater cell of inflow gauge
70 : !> \param[in] "integer(i4), dimension(:) :: InflowGaugeNodeList" gauge node list at L11
71 : !> \param[in] "real(dp), dimension(:) :: InflowDischarge" inflowing discharge at discharge gauge at
72 : !> current day
73 : !> \param[in] "integer(i4) :: nGauges" number of recording gauges
74 : !> \param[in] "integer(i4), dimension(:) :: gaugeIndexList" index list for outflow gauges
75 : !> \param[in] "integer(i4), dimension(:) :: gaugeNodeList" gauge node list at L11
76 : !> \param[in] "logical :: map_flag" flag indicating whether routing resolution
77 : !> iscoarser than hydrologic resolution
78 : !> \param[in] "real(dp), dimension(:) :: L11_length" L11 link length
79 : !> \param[in] "real(dp), dimension(:) :: L11_slope" L11 slope
80 : !> \param[in] "real(dp), dimension(:) :: L11_FracFPimp" L11 fraction of flood plain with impervios
81 : !> cover
82 :
83 : ! INTENT(INOUT)
84 : !> \param[inout] "real(dp), dimension(:) :: L11_C1" L11 muskingum parameter 1
85 : !> \param[inout] "real(dp), dimension(:) :: L11_C2" L11 muskingum parameter 2
86 : !> \param[inout] "real(dp), dimension(:) :: L11_qOut" total runoff from L11 grid cells
87 : !> \param[inout] "real(dp), dimension(:, :) :: L11_qTIN" L11 inflow to the reach
88 : !> \param[inout] "real(dp), dimension(:, :) :: L11_qTR" L11 routed outflow
89 : !> \param[inout] "real(dp), dimension(:) :: L11_qMod" modelled discharge at each grid cell
90 : !> \param[inout] "real(dp), dimension(:) :: GaugeDischarge" modelled discharge at each gauge
91 :
92 : ! HISTORY
93 : !> \authors Stephan Thober
94 :
95 : !> \date Aug 2015
96 :
97 : ! Modifications:
98 : ! Stephan Thober Sep 2015 - using arguments instead of global variables
99 : ! Stephan Thober Sep 2015 - added variables for routing resolution higher than hydrologic resolution
100 : ! Stephan Thober May 2016 - added check whether gauge is actually inside modelling domain before copying simulated runoff
101 : ! Stephan Thober Nov 2016 - implemented second routing process i.e. adaptive timestep
102 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
103 :
104 634236 : subroutine mRM_routing( &
105 634236 : read_states, processCase, global_routing_param, L1_total_runoff, L1_areaCell, L1_L11_Id, &
106 1268472 : L11_areaCell, L11_L1_Id, L11_netPerm, L11_fromN, L11_toN, L11_nOutlets, timestep, tsRoutFactor, &
107 634236 : nNodes, nInflowGauges, InflowGaugeIndexList, InflowGaugeHeadwater, InflowGaugeNodeList, &
108 1902708 : InflowDischarge, nGauges, gaugeIndexList, gaugeNodeList, map_flag, L11_length, L11_slope, &
109 2536944 : L11_FracFPimp, sink_cells, L11_C1, L11_C2, L11_qOut, L11_qTIN, L11_qTR, L11_qMod, GaugeDischarge &
110 : )
111 :
112 : use mo_constants, only : T0_dp
113 : use mo_mrm_global_variables, only : riv_temp_pcs, is_start
114 : use mo_mrm_mpr, only : reg_rout
115 : use mo_mrm_pre_routing, only : L11_runoff_acc, add_inflow
116 : ! use mo_mrm_riv_temp_class, only : riv_temp_type
117 :
118 : implicit none
119 :
120 : ! whether states are derived from restart file
121 : logical, intent(in) :: read_states
122 : ! Process switch for routing
123 : integer(i4), intent(in) :: processCase
124 : ! routing parameters
125 : real(dp), dimension(:), intent(in) :: global_routing_param
126 : ! total runoff from L1 grid cells
127 : real(dp), dimension(:), intent(in) :: L1_total_runoff
128 : ! L1 cell area
129 : real(dp), dimension(:), intent(in) :: L1_areaCell
130 : ! L1 cell ids on L11
131 : integer(i4), dimension(:), intent(in) :: L1_L11_Id
132 : ! L11 cell area
133 : real(dp), dimension(:), intent(in) :: L11_areaCell
134 : ! L11 cell ids on L1
135 : integer(i4), dimension(:), intent(in) :: L11_L1_Id
136 : ! L11 routing order
137 : integer(i4), dimension(:), intent(in) :: L11_netPerm
138 : ! L11 source grid cell order
139 : integer(i4), dimension(:), intent(in) :: L11_fromN
140 : ! L11 target grid cell order
141 : integer(i4), dimension(:), intent(in) :: L11_toN
142 : ! L11 number of outlets/sinks
143 : integer(i4), intent(in) :: L11_nOutlets
144 : ! simulation timestep in [h]
145 : integer(i4), intent(in) :: timestep
146 : ! factor between routing timestep and hydrological timestep
147 : real(dp), intent(in) :: tsRoutFactor
148 : ! number of nodes
149 : integer(i4), intent(in) :: nNodes
150 : ! number of inflow gauges
151 : integer(i4), intent(in) :: nInflowGauges
152 : ! index list of inflow gauges
153 : integer(i4), dimension(:), intent(in) :: InflowGaugeIndexList
154 : ! flag for headwater cell of inflow gauge
155 : logical, dimension(:), intent(in) :: InflowGaugeHeadwater
156 : ! gauge node list at L11
157 : integer(i4), dimension(:), intent(in) :: InflowGaugeNodeList
158 : ! inflowing discharge at discharge gauge at current day
159 : real(dp), dimension(:), intent(in) :: InflowDischarge
160 : ! number of recording gauges
161 : integer(i4), intent(in) :: nGauges
162 : ! index list for outflow gauges
163 : integer(i4), dimension(:), intent(in) :: gaugeIndexList
164 : ! gauge node list at L11
165 : integer(i4), dimension(:), intent(in) :: gaugeNodeList
166 : ! flag indicating whether routing resolution iscoarser than hydrologic resolution
167 : logical, intent(in) :: map_flag
168 : ! L11 link length
169 : real(dp), dimension(:), intent(in) :: L11_length
170 : ! L11 slope
171 : real(dp), dimension(:), intent(in) :: L11_slope
172 : ! L11 fraction of flood plain with impervios cover
173 : real(dp), dimension(:), intent(in) :: L11_FracFPimp
174 : ! [-] sink nodes
175 : integer(i4), dimension(:), intent(in) :: sink_cells
176 : ! L11 muskingum parameter 1
177 : real(dp), dimension(:), intent(inout) :: L11_C1
178 : ! L11 muskingum parameter 2
179 : real(dp), dimension(:), intent(inout) :: L11_C2
180 : ! total runoff from L11 grid cells
181 : real(dp), dimension(:), intent(inout) :: L11_qOut
182 : ! L11 inflow to the reach
183 : real(dp), dimension(:, :), intent(inout) :: L11_qTIN
184 : ! L11 routed outflow
185 : real(dp), dimension(:, :), intent(inout) :: L11_qTR
186 : ! modelled discharge at each grid cell
187 : real(dp), dimension(:), intent(inout) :: L11_qMod
188 : ! modelled discharge at each gauge
189 : real(dp), dimension(:), intent(inout) :: GaugeDischarge
190 :
191 : integer(i4) :: s11, e11 ! only for riv temp routing
192 : integer(i4) :: gg
193 : integer(i4) :: tt
194 : ! number of routing loops
195 : integer(i4) :: rout_loop
196 : ! variable for accumulation
197 28516260 : real(dp), dimension(size(L11_qMod, dim = 1)) :: L11_qAcc
198 634236 : real(dp), dimension(:), allocatable :: L11_E_Acc
199 :
200 634236 : if ( riv_temp_pcs%active ) then
201 : ! allocate accumulated temperature energy
202 26208 : allocate(L11_E_Acc(size(L11_qMod, dim = 1)))
203 229320 : L11_E_Acc = 0._dp ! init to zero
204 : ! get shortcuts for start end ending of current L11-domain
205 6552 : s11 = riv_temp_pcs%s11
206 6552 : e11 = riv_temp_pcs%e11
207 : end if
208 :
209 634236 : if (is_start) is_start = .false.
210 :
211 : ! this is using the sealed fraction for determining the routing parameters
212 : ! MPR has already been done
213 634236 : if (processCase .eq. 1_i4 .AND. (.not. read_states)) then
214 : ! for a single node model run
215 495096 : if (nNodes .GT. 1) then
216 0 : call reg_rout(global_routing_param, &
217 495096 : L11_length, L11_slope, L11_FracFPimp(: nNodes - L11_nOutlets), &
218 990192 : real(timeStep, dp), L11_C1(: nNodes - L11_nOutlets), L11_C2(: nNodes - L11_nOutlets))
219 : end if
220 : end if
221 :
222 : ! =====================================================================
223 : ! NOW, EXECUTE ROUTING
224 : ! ====================================================================
225 : ! calculate number of routing loops
226 634236 : rout_loop = max(1_i4, nint(1._dp / tsRoutFactor))
227 :
228 : ! runoff accumulation from L1 to L11 level
229 : call L11_runoff_acc( &
230 : L1_total_runoff, L1_areaCell, L1_L11_Id, &
231 : L11_areaCell, L11_L1_Id, timeStep, & ! Intent IN
232 : map_flag, & ! Intent IN
233 : L11_qOut & ! Intent OUT
234 634236 : )
235 : ! add inflow
236 : call add_inflow( &
237 : nInflowGauges, &
238 : InflowGaugeIndexList, &
239 : InflowGaugeHeadwater, &
240 : InflowGaugeNodeList, &
241 : InflowDischarge, & ! Intent IN
242 : L11_qOUT & ! Intent INOUT
243 634236 : )
244 : ! for a single node model run
245 634236 : if(nNodes .GT. 1) then
246 28516260 : L11_qAcc = 0._dp
247 : ! routing multiple times if timestep is smaller than 1
248 1268472 : do tt = 1, rout_loop
249 : ! routing of water within river reaches
250 : call L11_routing( &
251 : nNodes, &
252 : nNodes - L11_nOutlets, &
253 : L11_netPerm, &
254 : L11_fromN, & ! Intent IN
255 : L11_toN, & ! Intent IN
256 : L11_C1, & ! Intent IN
257 : L11_C2, & ! Intent IN
258 : L11_qOut, & ! Intent IN
259 : nInflowGauges, & ! Intent IN
260 : InflowGaugeHeadwater, & ! Intent IN
261 : InflowGaugeNodeList, & ! Intent IN
262 : sink_cells, & ! Intent IN
263 : L11_qTIN, & ! Intent INOUT
264 : L11_qTR, & ! Intent INOUT
265 0 : L11_Qmod & ! Intent OUT
266 634236 : )
267 : ! accumulate values of individual subtimesteps
268 28516260 : L11_qAcc = L11_qAcc + L11_qMod
269 : ! do the temperature routing
270 1268472 : if ( riv_temp_pcs%active ) then
271 : call riv_temp_pcs%L11_routing_E( &
272 : nNodes - L11_nOutlets, &
273 : L11_netPerm, &
274 : L11_fromN, & ! Intent IN
275 : L11_toN, & ! Intent IN
276 : L11_C1, & ! Intent IN
277 : L11_C2, & ! Intent IN
278 : nInflowGauges, & ! Intent IN
279 : InflowGaugeHeadwater, & ! Intent IN
280 : InflowGaugeNodeList, & ! Intent IN
281 : sink_cells, & ! Intent IN
282 : L11_qTR(:, 1), & ! Intent IN
283 : L11_Qmod & ! Intent IN
284 6552 : )
285 : end if
286 : end do
287 : ! calculate mean over routing period (timestep)
288 28516260 : L11_qMod = L11_qAcc / real(rout_loop, dp)
289 : else
290 0 : L11_Qmod = L11_qOUT
291 0 : if ( riv_temp_pcs%active ) riv_temp_pcs%river_temp(s11 : e11) = &
292 0 : max(riv_temp_pcs%delta_T, riv_temp_pcs%netNode_E_out(s11 : e11) / L11_qOUT - T0_dp)
293 : end if
294 :
295 : !----------------------------------------------------------------------
296 : ! FOR STORING the optional arguments
297 : !
298 : ! FOR RUNOFF
299 : ! NOTE:: Node ID for a given gauging station is stored at gaugeindex's
300 : ! index in runoff. In consequence the gauges in runoff are
301 : ! ordered corresponing to gauge%Q(:,:)
302 : !----------------------------------------------------------------------
303 1294704 : do gg = 1, nGauges
304 1294704 : GaugeDischarge(gaugeIndexList(gg)) = L11_Qmod(gaugeNodeList(gg))
305 : end do
306 :
307 634236 : end subroutine mRM_routing
308 :
309 : ! ------------------------------------------------------------------
310 :
311 : ! NAME
312 : ! L11_routing
313 :
314 : ! PURPOSE
315 : !> \brief Performs runoff routing for mHM at L11 upscaled network
316 : !> (\ref fig_routing "Routing Network").
317 : !> \details
318 : !> Hydrograph routing is carried out with the Muskingum algorithm
319 : !> \cite CMM1988. This simplification of the St. Venant
320 : !> equations is justified in mHM because the potential areas of
321 : !> application of this model would hardly exhibit abruptly
322 : !> changing hydrographs with supercritical flows. The discharge
323 : !> leaving the river reach located on cell \f$ i \f$ \f$
324 : !> Q_{i}^{1}(t) \f$ at time step \f$ t \f$ can be determined by
325 : !> \f[ Q_{i}^{1}(t) = Q_{i}^{1}(t-1)
326 : !> + c_{1} \left( Q_{i}^{0}(t-1) - Q_{i}^{1}(t-1) \right)
327 : !> + c_{2} \left( Q_{i}^{0}(t) - Q_{i}^{0}(t-1) \right) \f]
328 : !> with
329 : !> \f[ Q_{i}^{0}(t) = Q_{i'}(t) + Q_{i'}^{1}(t) \f]
330 : !> \f[ c_{1}= \frac{\Delta t} { \kappa (1- \xi ) + \frac{\Delta t}{2} } \f]
331 : !> \f[ c_{2}= \frac{ \frac{\Delta t}{2} - \kappa \xi} { \kappa (1- \xi)
332 : !> + \frac{\Delta t}{2} } \f]
333 : !> where
334 : !> \f$ Q_{i}^{0} \f$ and \f$ Q_{i}^{1} \f$ denote the discharge
335 : !> entering and leaving the river reach located on cell \f$ i \f$
336 : !> respectively.
337 : !> \f$ Q_{i'} \f$ is the contribution from the upstream cell \f$
338 : !> i'\f$.
339 : !> \f$ \kappa \f$ Muskingum travel time parameter.
340 : !> \f$ \xi \f$ Muskingum attenuation parameter.
341 : !> \f$ \Delta t \f$ time interval in hours.
342 : !> \f$ t \f$ Time index for each \f$ \Delta t \f$ interval.
343 : !> To improve performance, a routing sequence "netPerm" is
344 : !> required. This permutation is determined in the mo_init_mrm
345 : !> routine.
346 :
347 : !> \details TODO: add description
348 :
349 : ! INTENT(IN)
350 : !> \param[in] "integer(i4) :: nNodes" number of network nodes = nCells1
351 : !> \param[in] "integer(i4) :: nLinks" number of stream segment (reaches)
352 : !> \param[in] "integer(i4), dimension(:) :: netPerm" routing order of a given domain (permutation)
353 : !> \param[in] "integer(i4), dimension(:) :: netLink_fromN" from node
354 : !> \param[in] "integer(i4), dimension(:) :: netLink_toN" to node
355 : !> \param[in] "real(dp), dimension(:) :: netLink_C1" routing parameter C1 (\cite CMM1988 p. 25-41)
356 : !> \param[in] "real(dp), dimension(:) :: netLink_C2" routing parameters C2 (id)
357 : !> \param[in] "real(dp), dimension(:) :: netNode_qOUT" Total outflow from cells (given domain) L11 at time
358 : !> tt in [m3 s-1]
359 : !> \param[in] "integer(i4) :: nInflowGauges" [-] number of inflow points
360 : !> \param[in] "logical, dimension(:) :: InflowHeadwater" [-] if to consider headwater cells of inflow
361 : !> gauge
362 : !> \param[in] "integer(i4), dimension(:) :: InflowNodeList" [-] L11 ID of inflow points
363 :
364 : ! INTENT(INOUT)
365 : !> \param[inout] "real(dp), dimension(:, :) :: netNode_qTIN" [m3 s-1] Total inputs at t-1 and t
366 : !> \param[inout] "real(dp), dimension(:, :) :: netNode_qTR" [m3 s-1] Transformed outflow leaving
367 : !> node I (Muskingum)
368 :
369 : ! INTENT(OUT)
370 : !> \param[out] "real(dp), dimension(nNodes) :: netNode_Qmod" [m3 s-1] Simulated routed discharge
371 :
372 : ! HISTORY
373 : !> \authors Luis Samaniego
374 :
375 : !> \date Dec 2005
376 :
377 : ! Modifications:
378 : ! Luis Samaniego Feb 2008 - routing module (cells)
379 : ! Rohini Kumar Aug 2011 - vector version of mHM-UFZ
380 : ! Nov 2011 - parallel version
381 : ! Luis Samaniego Jan 2013 - modularization, documentation
382 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
383 :
384 1902708 : subroutine L11_routing(nNodes, nLinks, netPerm, netLink_fromN, netLink_toN, netLink_C1, netLink_C2, netNode_qOUT, &
385 1902708 : nInflowGauges, InflowHeadwater, InflowNodeList, sink_cells, netNode_qTIN, netNode_qTR, netNode_Qmod)
386 : implicit none
387 :
388 : ! number of network nodes = nCells1
389 : integer(i4), intent(in) :: nNodes
390 : ! number of stream segment (reaches)
391 : integer(i4), intent(in) :: nLinks
392 : ! routing order of a given domain (permutation)
393 : integer(i4), dimension(:), intent(in) :: netPerm
394 : ! from node
395 : integer(i4), dimension(:), intent(in) :: netLink_fromN
396 : ! to node
397 : integer(i4), dimension(:), intent(in) :: netLink_toN
398 : ! routing parameter C1 (\cite CMM1988 p. 25-41)
399 : real(dp), dimension(:), intent(in) :: netLink_C1
400 : ! routing parameters C2 (id)
401 : real(dp), dimension(:), intent(in) :: netLink_C2
402 : ! Total outflow from cells (given domain) L11 at time tt in [m3 s-1]
403 : real(dp), dimension(:), intent(in) :: netNode_qOUT
404 : ! [-] number of inflow points
405 : integer(i4), intent(in) :: nInflowGauges
406 : ! [-] if to consider headwater cells of inflow gauge
407 : logical, dimension(:), intent(in) :: InflowHeadwater
408 : ! [-] L11 ID of inflow points
409 : integer(i4), dimension(:), intent(in) :: InflowNodeList
410 : ! [-] sink nodes
411 : integer(i4), dimension(:), intent(in) :: sink_cells
412 : ! [m3 s-1] Total inputs at t-1 and t
413 : real(dp), dimension(:, :), intent(inout) :: netNode_qTIN
414 : ! [m3 s-1] Transformed outflow leaving node I (Muskingum)
415 : real(dp), dimension(:, :), intent(inout) :: netNode_qTR
416 : ! [m3 s-1] Simulated routed discharge
417 : real(dp), dimension(nNodes), intent(out) :: netNode_Qmod
418 :
419 : integer(i4) :: g, i, k, iNode, tNode
420 : ! current routing state (2)
421 : integer(i4), parameter :: IT = 2
422 : ! past routing state (1)
423 : integer(i4), parameter :: IT1 = 1
424 :
425 : ! Entry value for the auxiliary vectors
426 : ! netNode_qTIN(iNode,:)
427 : ! netNode_qTR(iNode,:)
428 : ! which store current and past states of
429 : ! incoming and outgoing of discharge at iNode
430 : !--------------------------------------------------------------------------
431 : ! Muskingum Flood Routing
432 : !--------------------------------------------------------------------------
433 : ! initialize total input at point time IT in all nodes
434 28516260 : netNode_qTIN(:, IT) = 0.0_dp
435 : !--------------------------------------------------------------------------
436 : ! Links in sequential mode .... with single node
437 : !--------------------------------------------------------------------------
438 : ! ST - decent parallelization has to be done!!!
439 : !!$OMP parallel
440 : !!$OMP do private(g, i, inode, tnode)
441 27882024 : do k = 1, nLinks
442 : ! get LINK routing order -> i
443 27247788 : i = netPerm(k)
444 27247788 : iNode = netLink_fromN(i)
445 27247788 : tNode = netLink_toN(i)
446 :
447 : ! accumulate all inputs in iNode
448 27247788 : netNode_qTIN(iNode, IT) = netNode_qTIN(iNode, IT) + netNode_qOUT(iNode)
449 :
450 : ! routing iNode
451 27247788 : netNode_qTR(iNode, IT) = netNode_qTR(iNode, IT1) &
452 27247788 : + netLink_C1(i) * (netNode_qTIN(iNode, IT1) - netNode_qTR (iNode, IT1)) &
453 27247788 : + netLink_C2(i) * (netNode_qTIN(iNode, IT) - netNode_qTIN(iNode, IT1))
454 :
455 : ! check if the inflow from upstream cells should be deactivated
456 27247788 : if (nInflowGauges .GT. 0) then
457 1729728 : do g = 1, nInflowGauges
458 : ! check if downstream Node (tNode) is inflow gauge and headwaters should be ignored
459 1729728 : if ((tNode == InflowNodeList(g)) .AND. (.NOT. InflowHeadwater(g))) netNode_qTR(iNode, IT) = 0.0_dp
460 : end do
461 : end if
462 :
463 : ! add routed water to downstream node
464 27882024 : netNode_qTIN(tNode, IT) = netNode_qTIN(tNode, IT) + netNode_qTR(iNode, IT)
465 : end do
466 : !!$OMP end do
467 : !!$OMP end parallel
468 :
469 : !--------------------------------------------------------------------------
470 : ! add runoff to sinks
471 : !--------------------------------------------------------------------------
472 1268472 : do k = 1, size(sink_cells)
473 634236 : tNode = sink_cells(k)
474 1268472 : netNode_qTIN(tNode, IT) = netNode_qTIN(tNode, IT) + netNode_qOUT(tNode)
475 : end do
476 : !--------------------------------------------------------------------------
477 : ! save modeled discharge at time step tt then shift flow storages
478 : ! (NOTE aggregation to daily values to be done outside)
479 : !--------------------------------------------------------------------------
480 : ! !!$OMP parallel
481 : ! store generated discharge
482 29150496 : netNode_Qmod(1 : nNodes) = netNode_qTIN(1 : nNodes, IT)
483 : ! backflow t-> t-1
484 29150496 : netNode_qTR(1 : nNodes, IT1) = netNode_qTR(1 : nNodes, IT)
485 29150496 : netNode_qTIN(1 : nNodes, IT1) = netNode_qTIN(1 : nNodes, IT)
486 : ! !!$OMP end parallel
487 :
488 634236 : end subroutine L11_routing
489 :
490 : END MODULE mo_mrm_routing
|
