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, 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 : ! L11 muskingum parameter 1
175 : real(dp), dimension(:), intent(inout) :: L11_C1
176 : ! L11 muskingum parameter 2
177 : real(dp), dimension(:), intent(inout) :: L11_C2
178 : ! total runoff from L11 grid cells
179 : real(dp), dimension(:), intent(inout) :: L11_qOut
180 : ! L11 inflow to the reach
181 : real(dp), dimension(:, :), intent(inout) :: L11_qTIN
182 : ! L11 routed outflow
183 : real(dp), dimension(:, :), intent(inout) :: L11_qTR
184 : ! modelled discharge at each grid cell
185 : real(dp), dimension(:), intent(inout) :: L11_qMod
186 : ! modelled discharge at each gauge
187 : real(dp), dimension(:), intent(inout) :: GaugeDischarge
188 :
189 : integer(i4) :: s11, e11 ! only for riv temp routing
190 : integer(i4) :: gg
191 : integer(i4) :: tt
192 : ! number of routing loops
193 : integer(i4) :: rout_loop
194 : ! variable for accumulation
195 28516260 : real(dp), dimension(size(L11_qMod, dim = 1)) :: L11_qAcc
196 634236 : real(dp), dimension(:), allocatable :: L11_E_Acc
197 :
198 634236 : if ( riv_temp_pcs%active ) then
199 : ! allocate accumulated temperature energy
200 26208 : allocate(L11_E_Acc(size(L11_qMod, dim = 1)))
201 229320 : L11_E_Acc = 0._dp ! init to zero
202 : ! get shortcuts for start end ending of current L11-domain
203 6552 : s11 = riv_temp_pcs%s11
204 6552 : e11 = riv_temp_pcs%e11
205 : end if
206 :
207 634236 : if (is_start) is_start = .false.
208 :
209 : ! this is using the sealed fraction for determining the routing parameters
210 : ! MPR has already been done
211 634236 : if (processCase .eq. 1_i4 .AND. (.not. read_states)) then
212 : ! for a single node model run
213 495096 : if (nNodes .GT. 1) then
214 0 : call reg_rout(global_routing_param, &
215 495096 : L11_length, L11_slope, L11_FracFPimp(: nNodes - L11_nOutlets), &
216 990192 : real(timeStep, dp), L11_C1(: nNodes - L11_nOutlets), L11_C2(: nNodes - L11_nOutlets))
217 : end if
218 : end if
219 :
220 : ! =====================================================================
221 : ! NOW, EXECUTE ROUTING
222 : ! ====================================================================
223 : ! calculate number of routing loops
224 634236 : rout_loop = max(1_i4, nint(1._dp / tsRoutFactor))
225 :
226 : ! runoff accumulation from L1 to L11 level
227 : call L11_runoff_acc( &
228 : L1_total_runoff, L1_areaCell, L1_L11_Id, &
229 : L11_areaCell, L11_L1_Id, timeStep, & ! Intent IN
230 : map_flag, & ! Intent IN
231 : L11_qOut & ! Intent OUT
232 634236 : )
233 : ! add inflow
234 : call add_inflow( &
235 : nInflowGauges, &
236 : InflowGaugeIndexList, &
237 : InflowGaugeHeadwater, &
238 : InflowGaugeNodeList, &
239 : InflowDischarge, & ! Intent IN
240 : L11_qOUT & ! Intent INOUT
241 634236 : )
242 : ! for a single node model run
243 634236 : if(nNodes .GT. 1) then
244 28516260 : L11_qAcc = 0._dp
245 : ! routing multiple times if timestep is smaller than 1
246 1268472 : do tt = 1, rout_loop
247 : ! routing of water within river reaches
248 : call L11_routing( &
249 : nNodes, &
250 : nNodes - L11_nOutlets, &
251 : L11_netPerm, &
252 : L11_fromN, & ! Intent IN
253 : L11_toN, & ! Intent IN
254 : L11_C1, & ! Intent IN
255 : L11_C2, & ! Intent IN
256 : L11_qOut, & ! Intent IN
257 : nInflowGauges, & ! Intent IN
258 : InflowGaugeHeadwater, & ! Intent IN
259 : InflowGaugeNodeList, & ! Intent IN
260 : L11_qTIN, & ! Intent INOUT
261 : L11_qTR, & ! Intent INOUT
262 0 : L11_Qmod & ! Intent OUT
263 634236 : )
264 : ! accumulate values of individual subtimesteps
265 28516260 : L11_qAcc = L11_qAcc + L11_qMod
266 : ! do the temperature routing
267 1268472 : if ( riv_temp_pcs%active ) then
268 : call riv_temp_pcs%L11_routing_E( &
269 : nNodes - L11_nOutlets, &
270 : L11_netPerm, &
271 : L11_fromN, & ! Intent IN
272 : L11_toN, & ! Intent IN
273 : L11_C1, & ! Intent IN
274 : L11_C2, & ! Intent IN
275 : nInflowGauges, & ! Intent IN
276 : InflowGaugeHeadwater, & ! Intent IN
277 : InflowGaugeNodeList, & ! Intent IN
278 : L11_qTR(:, 1), & ! Intent IN
279 : L11_Qmod & ! Intent IN
280 6552 : )
281 : end if
282 : end do
283 : ! calculate mean over routing period (timestep)
284 28516260 : L11_qMod = L11_qAcc / real(rout_loop, dp)
285 : else
286 0 : L11_Qmod = L11_qOUT
287 0 : if ( riv_temp_pcs%active ) riv_temp_pcs%river_temp(s11 : e11) = &
288 0 : max(riv_temp_pcs%delta_T, riv_temp_pcs%netNode_E_out(s11 : e11) / L11_qOUT - T0_dp)
289 : end if
290 :
291 : !----------------------------------------------------------------------
292 : ! FOR STORING the optional arguments
293 : !
294 : ! FOR RUNOFF
295 : ! NOTE:: Node ID for a given gauging station is stored at gaugeindex's
296 : ! index in runoff. In consequence the gauges in runoff are
297 : ! ordered corresponing to gauge%Q(:,:)
298 : !----------------------------------------------------------------------
299 1294704 : do gg = 1, nGauges
300 1294704 : GaugeDischarge(gaugeIndexList(gg)) = L11_Qmod(gaugeNodeList(gg))
301 : end do
302 :
303 634236 : end subroutine mRM_routing
304 :
305 : ! ------------------------------------------------------------------
306 :
307 : ! NAME
308 : ! L11_routing
309 :
310 : ! PURPOSE
311 : !> \brief Performs runoff routing for mHM at L11 upscaled network
312 : !> (\ref fig_routing "Routing Network").
313 : !> \details
314 : !> Hydrograph routing is carried out with the Muskingum algorithm
315 : !> \cite CMM1988. This simplification of the St. Venant
316 : !> equations is justified in mHM because the potential areas of
317 : !> application of this model would hardly exhibit abruptly
318 : !> changing hydrographs with supercritical flows. The discharge
319 : !> leaving the river reach located on cell \f$ i \f$ \f$
320 : !> Q_{i}^{1}(t) \f$ at time step \f$ t \f$ can be determined by
321 : !> \f[ Q_{i}^{1}(t) = Q_{i}^{1}(t-1)
322 : !> + c_{1} \left( Q_{i}^{0}(t-1) - Q_{i}^{1}(t-1) \right)
323 : !> + c_{2} \left( Q_{i}^{0}(t) - Q_{i}^{0}(t-1) \right) \f]
324 : !> with
325 : !> \f[ Q_{i}^{0}(t) = Q_{i'}(t) + Q_{i'}^{1}(t) \f]
326 : !> \f[ c_{1}= \frac{\Delta t} { \kappa (1- \xi ) + \frac{\Delta t}{2} } \f]
327 : !> \f[ c_{2}= \frac{ \frac{\Delta t}{2} - \kappa \xi} { \kappa (1- \xi)
328 : !> + \frac{\Delta t}{2} } \f]
329 : !> where
330 : !> \f$ Q_{i}^{0} \f$ and \f$ Q_{i}^{1} \f$ denote the discharge
331 : !> entering and leaving the river reach located on cell \f$ i \f$
332 : !> respectively.
333 : !> \f$ Q_{i'} \f$ is the contribution from the upstream cell \f$
334 : !> i'\f$.
335 : !> \f$ \kappa \f$ Muskingum travel time parameter.
336 : !> \f$ \xi \f$ Muskingum attenuation parameter.
337 : !> \f$ \Delta t \f$ time interval in hours.
338 : !> \f$ t \f$ Time index for each \f$ \Delta t \f$ interval.
339 : !> To improve performance, a routing sequence "netPerm" is
340 : !> required. This permutation is determined in the mo_init_mrm
341 : !> routine.
342 :
343 : !> \details TODO: add description
344 :
345 : ! INTENT(IN)
346 : !> \param[in] "integer(i4) :: nNodes" number of network nodes = nCells1
347 : !> \param[in] "integer(i4) :: nLinks" number of stream segment (reaches)
348 : !> \param[in] "integer(i4), dimension(:) :: netPerm" routing order of a given domain (permutation)
349 : !> \param[in] "integer(i4), dimension(:) :: netLink_fromN" from node
350 : !> \param[in] "integer(i4), dimension(:) :: netLink_toN" to node
351 : !> \param[in] "real(dp), dimension(:) :: netLink_C1" routing parameter C1 (\cite CMM1988 p. 25-41)
352 : !> \param[in] "real(dp), dimension(:) :: netLink_C2" routing parameters C2 (id)
353 : !> \param[in] "real(dp), dimension(:) :: netNode_qOUT" Total outflow from cells (given domain) L11 at time
354 : !> tt in [m3 s-1]
355 : !> \param[in] "integer(i4) :: nInflowGauges" [-] number of inflow points
356 : !> \param[in] "logical, dimension(:) :: InflowHeadwater" [-] if to consider headwater cells of inflow
357 : !> gauge
358 : !> \param[in] "integer(i4), dimension(:) :: InflowNodeList" [-] L11 ID of inflow points
359 :
360 : ! INTENT(INOUT)
361 : !> \param[inout] "real(dp), dimension(:, :) :: netNode_qTIN" [m3 s-1] Total inputs at t-1 and t
362 : !> \param[inout] "real(dp), dimension(:, :) :: netNode_qTR" [m3 s-1] Transformed outflow leaving
363 : !> node I (Muskingum)
364 :
365 : ! INTENT(OUT)
366 : !> \param[out] "real(dp), dimension(nNodes) :: netNode_Qmod" [m3 s-1] Simulated routed discharge
367 :
368 : ! HISTORY
369 : !> \authors Luis Samaniego
370 :
371 : !> \date Dec 2005
372 :
373 : ! Modifications:
374 : ! Luis Samaniego Feb 2008 - routing module (cells)
375 : ! Rohini Kumar Aug 2011 - vector version of mHM-UFZ
376 : ! Nov 2011 - parallel version
377 : ! Luis Samaniego Jan 2013 - modularization, documentation
378 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
379 :
380 1902708 : subroutine L11_routing(nNodes, nLinks, netPerm, netLink_fromN, netLink_toN, netLink_C1, netLink_C2, netNode_qOUT, &
381 1902708 : nInflowGauges, InflowHeadwater, InflowNodeList, netNode_qTIN, netNode_qTR, netNode_Qmod)
382 : implicit none
383 :
384 : ! number of network nodes = nCells1
385 : integer(i4), intent(in) :: nNodes
386 : ! number of stream segment (reaches)
387 : integer(i4), intent(in) :: nLinks
388 : ! routing order of a given domain (permutation)
389 : integer(i4), dimension(:), intent(in) :: netPerm
390 : ! from node
391 : integer(i4), dimension(:), intent(in) :: netLink_fromN
392 : ! to node
393 : integer(i4), dimension(:), intent(in) :: netLink_toN
394 : ! routing parameter C1 (\cite CMM1988 p. 25-41)
395 : real(dp), dimension(:), intent(in) :: netLink_C1
396 : ! routing parameters C2 (id)
397 : real(dp), dimension(:), intent(in) :: netLink_C2
398 : ! Total outflow from cells (given domain) L11 at time tt in [m3 s-1]
399 : real(dp), dimension(:), intent(in) :: netNode_qOUT
400 : ! [-] number of inflow points
401 : integer(i4), intent(in) :: nInflowGauges
402 : ! [-] if to consider headwater cells of inflow gauge
403 : logical, dimension(:), intent(in) :: InflowHeadwater
404 : ! [-] L11 ID of inflow points
405 : integer(i4), dimension(:), intent(in) :: InflowNodeList
406 : ! [m3 s-1] Total inputs at t-1 and t
407 : real(dp), dimension(:, :), intent(inout) :: netNode_qTIN
408 : ! [m3 s-1] Transformed outflow leaving node I (Muskingum)
409 : real(dp), dimension(:, :), intent(inout) :: netNode_qTR
410 : ! [m3 s-1] Simulated routed discharge
411 : real(dp), dimension(nNodes), intent(out) :: netNode_Qmod
412 :
413 : integer(i4) :: g, i, k, iNode, tNode
414 : ! current routing state (2)
415 : integer(i4), parameter :: IT = 2
416 : ! past routing state (1)
417 : integer(i4), parameter :: IT1 = 1
418 :
419 : ! Entry value for the auxiliary vectors
420 : ! netNode_qTIN(iNode,:)
421 : ! netNode_qTR(iNode,:)
422 : ! which store current and past states of
423 : ! incoming and outgoing of discharge at iNode
424 : !--------------------------------------------------------------------------
425 : ! Muskingum Flood Routing
426 : !--------------------------------------------------------------------------
427 : ! initialize total input at point time IT in all nodes
428 28516260 : netNode_qTIN(:, IT) = 0.0_dp
429 : !--------------------------------------------------------------------------
430 : ! Links in sequential mode .... with single node
431 : !--------------------------------------------------------------------------
432 : ! ST - decent parallelization has to be done!!!
433 : !!$OMP parallel
434 : !!$OMP do private(g, i, inode, tnode)
435 27882024 : do k = 1, nLinks
436 : ! get LINK routing order -> i
437 27247788 : i = netPerm(k)
438 27247788 : iNode = netLink_fromN(i)
439 27247788 : tNode = netLink_toN(i)
440 :
441 : ! accumulate all inputs in iNode
442 27247788 : netNode_qTIN(iNode, IT) = netNode_qTIN(iNode, IT) + netNode_qOUT(iNode)
443 :
444 : ! routing iNode
445 27247788 : netNode_qTR(iNode, IT) = netNode_qTR(iNode, IT1) &
446 27247788 : + netLink_C1(i) * (netNode_qTIN(iNode, IT1) - netNode_qTR (iNode, IT1)) &
447 27247788 : + netLink_C2(i) * (netNode_qTIN(iNode, IT) - netNode_qTIN(iNode, IT1))
448 :
449 : ! check if the inflow from upstream cells should be deactivated
450 27247788 : if (nInflowGauges .GT. 0) then
451 1729728 : do g = 1, nInflowGauges
452 : ! check if downstream Node (tNode) is inflow gauge and headwaters should be ignored
453 1729728 : if ((tNode == InflowNodeList(g)) .AND. (.NOT. InflowHeadwater(g))) netNode_qTR(iNode, IT) = 0.0_dp
454 : end do
455 : end if
456 :
457 : ! add routed water to downstream node
458 27882024 : netNode_qTIN(tNode, IT) = netNode_qTIN(tNode, IT) + netNode_qTR(iNode, IT)
459 : end do
460 : !!$OMP end do
461 : !!$OMP end parallel
462 :
463 : !--------------------------------------------------------------------------
464 : ! Accumulate all inputs in tNode (netNode_qOUT) ONLY for last link
465 : !--------------------------------------------------------------------------
466 634236 : tNode = netLink_toN(netPerm(nLinks))
467 634236 : netNode_qTIN(tNode, IT) = netNode_qTIN(tNode, IT) + netNode_qOUT(tNode)
468 :
469 : !--------------------------------------------------------------------------
470 : ! save modeled discharge at time step tt then shift flow storages
471 : ! (NOTE aggregation to daily values to be done outside)
472 : !--------------------------------------------------------------------------
473 : ! !!$OMP parallel
474 : ! store generated discharge
475 29150496 : netNode_Qmod(1 : nNodes) = netNode_qTIN(1 : nNodes, IT)
476 : ! backflow t-> t-1
477 29150496 : netNode_qTR(1 : nNodes, IT1) = netNode_qTR(1 : nNodes, IT)
478 29150496 : netNode_qTIN(1 : nNodes, IT1) = netNode_qTIN(1 : nNodes, IT)
479 : ! !!$OMP end parallel
480 :
481 634236 : end subroutine L11_routing
482 :
483 : END MODULE mo_mrm_routing
|