5.13.3-dev0
mHM
The mesoscale Hydrological Model
Loading...
Searching...
No Matches
mo_mrm_net_startup.f90
Go to the documentation of this file.
1!> \file mo_mrm_net_startup.f90
2!> \brief \copybrief mo_mrm_net_startup
3!> \details \copydetails mo_mrm_net_startup
4
5!> \brief Startup drainage network for mHM.
6!> \details This module initializes the drainage network at L11 in mHM.
7!! - Delineation of drainage network at level 11.
8!! - Setting network topology (i.e. nodes and link).
9!! - Determining routing order.
10!! - Determining cell locations for network links.
11!! - Find drainage outlet.
12!! - Determine stream (links) features.
13!!
14!> \changelog
15!! - Rohini Kumar May 2014
16!! - cell area calulation based on a regular lat-lon grid or on a regular X-Y coordinate system
17!! - Robert Schweppe Jun 2018
18!! - refactoring and reformatting
19!> \authors Luis Samaniego
20!> \date Dec 2012
21!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
22!! mHM is released under the LGPLv3+ license \license_note
23!> \ingroup f_mrm
25 use mo_kind, only : i4, dp
26 use mo_message, only : message, error_message
27
28 implicit none
29
30 PUBLIC :: l11_l1_mapping
31 PUBLIC :: l11_flow_direction
33 PUBLIC :: l11_routing_order
34 PUBLIC :: l11_link_location
36 PUBLIC :: l11_stream_features
39 PUBLIC :: l11_flow_accumulation
40 PUBLIC :: l11_calc_celerity
41contains
42
43 ! NAME
44 ! L11_L1_mapping
45
46 ! PURPOSE
47 !> \brief TODO: add description
48
49 !> \details TODO: add description
50
51 ! INTENT(IN)
52 !> \param[in] "integer(i4) :: iDomain" domain
53
54 ! HISTORY
55 !> \authors Robert Schweppe
56
57 !> \date Jun 2018
58
59 ! Modifications:
60
61 subroutine l11_l1_mapping(iDomain)
62
63 use mo_append, only : append
65 use mo_common_variables, only : level1
67
68 implicit none
69
70 ! domain
71 integer(i4), intent(in) :: idomain
72
73 integer(i4) :: nrows1, ncols1
74
75 integer(i4) :: nrows11, ncols11
76
77 integer(i4) :: kk
78
79 integer(i4) :: icc, jcc
80
81 integer(i4) :: iu, id, jl, jr
82
83 ! mapping of L11 Id on L1
84 integer(i4), dimension(:, :), allocatable :: l11id_on_l1
85
86 ! mapping of L1 Id on L11
87 integer(i4), dimension(:, :), allocatable :: l1id_on_l11
88
89 ! dummy ID
90 integer(i4), dimension(:, :), allocatable :: dummy_2d_id
91
92 real(dp) :: cellfactorrbyh
93
94 integer(i4) :: cellfactorrbyh_inv
95
96
97 nrows1 = level1(idomain)%nrows
98 nrows11 = level11(idomain)%nrows
99 ncols1 = level1(idomain)%ncols
100 ncols11 = level11(idomain)%ncols
101
102 ! allocate variables for mapping L11 Ids and L1 Ids
103 allocate (l11id_on_l1(nrows1, ncols1))
104 allocate (l1id_on_l11(nrows11, ncols11))
105 l11id_on_l1(:, :) = nodata_i4
106 l1id_on_l11(:, :) = nodata_i4
107
108 ! set cell factor for routing
109 cellfactorrbyh = level11(idomain)%cellsize / level1(idomain)%cellsize
110
111 ! set mapping
112 ! create mapping between L11 and L1 for L11 resolution higher than L1 resolution
113 if (cellfactorrbyh .lt. 1._dp) then
114 allocate (dummy_2d_id(nrows1, ncols1))
115 dummy_2d_id = unpack(level1(idomain)%Id, level1(idomain)%mask, nodata_i4)
116 cellfactorrbyh_inv = int(1. / cellfactorrbyh, i4)
117 kk = 0
118 do jcc = 1, ncols1
119 do icc = 1, nrows1
120 if(.not. level1(idomain)%mask(icc, jcc)) cycle
121 kk = kk + 1
122 !
123 iu = (icc - 1) * cellfactorrbyh_inv + 1
124 id = min(icc * cellfactorrbyh_inv, nrows11)
125 jl = (jcc - 1) * cellfactorrbyh_inv + 1
126 jr = min(jcc * cellfactorrbyh_inv, ncols11)
127
128 l1id_on_l11(iu : id, jl : jr) = merge(dummy_2d_id(icc, jcc), nodata_i4, level11(idomain)%mask(iu : id, jl : jr))
129 end do
130 end do
131 else
132 allocate (dummy_2d_id(nrows11, ncols11))
133 dummy_2d_id = unpack(level11(idomain)%Id, level11(idomain)%mask, nodata_i4)
134
135 kk = 0
136 do jcc = 1, ncols11
137 do icc = 1, nrows11
138 if(.not. level11(idomain)%mask(icc, jcc)) cycle
139 kk = kk + 1
140
141 ! coord. of all corners L11 -> of finer scale level-1
142 iu = (icc - 1) * nint(cellfactorrbyh, i4) + 1
143 id = icc * nint(cellfactorrbyh, i4)
144 jl = (jcc - 1) * nint(cellfactorrbyh, i4) + 1
145 jr = jcc * nint(cellfactorrbyh, i4)
146
147 ! constrain the range of up, down, left, and right boundaries
148 if(iu < 1) iu = 1
149 if(id > nrows1) id = nrows1
150 if(jl < 1) jl = 1
151 if(jr > ncols1) jr = ncols1
152
153 ! Delimitation of level-11 cells on level-1 for L11 resolution lower than L1 resolution
154 l11id_on_l1(iu : id, jl : jr) = merge(dummy_2d_id(icc, jcc), nodata_i4, level1(idomain)%mask(iu : id, jl : jr))
155 end do
156 end do
157 end if
158
159 ! L1 data sets
160 call append(l1_l11_id, pack(l11id_on_l1(:, :), level1(idomain)%mask))
161 ! L11 data sets
162 call append(l11_l1_id, pack(l1id_on_l11(:, :), level11(idomain)%mask))
163 ! free space
164 deallocate(l11id_on_l1, l1id_on_l11, dummy_2d_id)
165
166 end subroutine l11_l1_mapping
167 ! --------------------------------------------------------------------------
168
169 ! NAME
170 ! L11_flow_direction
171
172 ! PURPOSE
173 !> \brief Determine the flow direction of the upscaled river
174 !> network at level L11.
175
176 !> \details The hydrographs generated at each cell are routed
177 !> through the drainage network at level-11 towards their
178 !> outlets. The drainage network at level-11 is conceptualized as a
179 !> graph whose nodes are hypothetically located at the center of
180 !> each grid cell connected by links that represent the river
181 !> reaches. The flow direction of a link correspond to the
182 !> direction towards a neighboring cell in which the net flow
183 !> accumulation (outflows minus inflows) attains its maximum
184 !> value. The net flow accumulation across a cell's boundary at
185 !> level-11 is estimated based on flow direction and flow
186 !> accumulation obtained at level-0 (\ref fig_routing "Routing
187 !> Network"). Note: level-1 denotes the modeling level, whereas
188 !> level-L11 is at least as coarse as level-1. Experience has
189 !> shown that routing can be done at a coarser resolution as
190 !> level-1, hence the level-11 was introduced.
191 !> \image html routing.png "Upscaling routing network from L0 to L1 (or L11)"
192 !> \anchor fig_routing \image latex routing.pdf "Upscaling routing network from L0 to L1 (or L11)" width=14cm
193 !> The left panel depicts a schematic derivation of a drainage
194 !> network at the level-11 based on level-0 flow direction and
195 !> flow accumulation. The dotted line circle denotes the point
196 !> with the highest flow accumulation within a grid cell. The
197 !> topology of a tipical drainage routing network at level-11 is
198 !> shown in the right panel. Gray color areas denote the flood
199 !> plains estimated in mo_init_mrm, where the network
200 !> upscaling is also carried out.
201 !> For the sake of simplicity, it is assumed that all runoff leaving
202 !> a given cell would exit through a major direction.
203 !> Note that multiple outlets can exist within the modelling domain.
204 !> If a variable is added or removed here, then it also has to
205 !> be added or removed in the subroutine L11_config_set in
206 !> module mo_restart and in the subroutine set_L11_config in module
207 !> mo_set_netcdf_restart
208 !> ADDITIONAL INFORMATION
209 !> L11_flow_direction
210
211 ! INTENT(IN)
212 !> \param[in] "integer(i4) :: iDomain" Domain Id
213
214 ! HISTORY
215 !> \authors Luis Samaniego
216
217 !> \date Dec 2005
218
219 ! Modifications:
220 ! Luis Samaniego Jan 2013 - modular version
221 ! Rohini Kumar Apr 2014 - Case of L0 is same as L11 implemented
222 ! Stephan Thober Aug 2015 - ported to mRM
223 ! Stephan Thober Sep 2015 - create mapping between L11 and L1 if L11 resolution is higher than L1 resolution
224 ! Stephan Thober May 2016 - introducing multiple outlets
225 ! Robert Schweppe Jun 2018 - refactoring and reformatting
226
227 subroutine l11_flow_direction(iDomain)
228
229 use mo_append, only : append
231 use mo_common_types, only: grid
235 use mo_string_utils, only : num2str
236
237 implicit none
238
239 ! Domain Id
240 integer(i4), intent(in) :: idomain
241
242 integer(i4) :: ncells0
243
244 integer(i4) :: nrows0, ncols0
245
246 integer(i4) :: s0, e0
247
248 integer(i4) :: nrows11, ncols11
249
250 ! = ncells11
251 integer(i4) :: nnodes
252
253 integer(i4) :: ii, jj, kk, ic, jc
254
255 integer(i4) :: iu, id
256
257 integer(i4) :: jl, jr
258
259 integer(i4) :: irow, jcol
260
261 integer(i4), dimension(:, :), allocatable :: id0
262
263 integer(i4), dimension(:, :), allocatable :: fdir0
264
265 integer(i4), dimension(:, :), allocatable :: facc0
266
267 integer(i4), dimension(:, :), allocatable :: fdir11
268
269 ! northing cell loc. of the Outlet
270 integer(i4), dimension(:), allocatable :: rowout
271
272 ! easting cell loc. of the Outlet
273 integer(i4), dimension(:), allocatable :: colout
274
275 integer(i4), dimension(:, :), allocatable :: drasc0
276
277 ! output location in L0
278 integer(i4), dimension(:, :), allocatable :: oloc
279
280 integer(i4) :: side
281
282 integer(i4) :: faccmax, idmax
283
284 ! Number of outlet found
285 integer(i4) :: noutlet
286
287 ! Number of outlets before this Domain
288 integer(i4) :: old_noutlet
289
290 ! flag whether outlet is found
291 logical :: is_outlet
292
293 type(grid), pointer :: level0_idomain => null()
294
295
296 !--------------------------------------------------------
297 ! STEPS:
298 ! 1) Estimate each variable locally for a given Domain
299 ! 2) Pad each variable to its corresponding global one
300 !--------------------------------------------------------
301 ! initialize
302 noutlet = 0_i4
303 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
304 !------------------------------------------------------------------
305 ! Set Flow Direction at Level 11
306 ! Searching order
307 ! jl jr
308 ! iu + 4 +
309 ! 3 1
310 ! id + 2 +
311 !------------------------------------------------------------------
312
313 ! flow direction at level-11
314
315 ! allocate
316 nrows0 = level0_idomain%nrows
317 ncols0 = level0_idomain%ncols
318 ncells0 = level0_idomain%ncells
319 nrows11 = level11(idomain)%nrows
320 ncols11 = level11(idomain)%ncols
321 nnodes = level11(idomain)%ncells
322 s0 = level0_idomain%iStart
323 e0 = level0_idomain%iEnd
324
325 allocate (id0(nrows0, ncols0))
326 allocate (facc0(nrows0, ncols0))
327 allocate (fdir0(nrows0, ncols0))
328 allocate (drasc0(nrows0, ncols0))
329 allocate (fdir11(nrows11, ncols11))
330 allocate (rowout(nnodes))
331 allocate (colout(nnodes))
332 allocate (oloc(1, 2))
333
334 ! initialize
335 id0(:, :) = nodata_i4
336 facc0(:, :) = nodata_i4
337 fdir0(:, :) = nodata_i4
338 drasc0(:, :) = nodata_i4
339 fdir11(:, :) = nodata_i4
340 rowout(:) = nodata_i4
341 colout(:) = nodata_i4
342 oloc(:, :) = nodata_i4
343
344
345 ! get iD, fAcc, fDir at L0
346 id0(:, :) = unpack(level0_idomain%Id, level0_idomain%mask, nodata_i4)
347 facc0(:, :) = unpack(l0_facc(s0 : e0), level0_idomain%mask, nodata_i4)
348 fdir0(:, :) = unpack(l0_fdir(s0 : e0), level0_idomain%mask, nodata_i4)
349
350 ! case where routing and input data scale is similar
351 IF(ncells0 .EQ. nnodes) THEN
352 oloc(1, :) = maxloc(facc0, level0_idomain%mask)
353 kk = l0_l11_remap(idomain)%lowres_id_on_highres(oloc(1, 1), oloc(1, 2))
354 ! for a single node model run
355 if(ncells0 .EQ. 1) then
356 fdir11(1, 1) = fdir0(oloc(1, 1), oloc(1, 2))
357 else
358 fdir11(:, :) = fdir0(:, :)
359 end if
360 fdir11(level11(idomain)%CellCoor(kk, 1), level11(idomain)%CellCoor(kk, 2)) = 0
361 ! set location of main outlet in L11
362 do kk = 1, nnodes
363 ii = level11(idomain)%CellCoor(kk, 1)
364 jj = level11(idomain)%CellCoor(kk, 2)
365 rowout(kk) = ii
366 colout(kk) = jj
367 end do
368 do kk = 1, ncells0
369 ii = level0_idomain%CellCoor(kk, 1)
370 jj = level0_idomain%CellCoor(kk, 2)
371 drasc0(ii, jj) = kk
372 end do
373
374 ! case where routing and input data scale differs
375 ELSE
376 ! =======================================================================
377 ! ST: find all cells whose downstream cells are outside the domain
378 ! =======================================================================
379 do ii = 1, ncells0
380 irow = level0_idomain%CellCoor(ii, 1)
381 jcol = level0_idomain%CellCoor(ii, 2)
382 call movedownonecell(fdir0(irow, jcol), irow, jcol)
383 ! check whether new location is inside bound
384 is_outlet = .false.
385 if ((irow .le. 0_i4) .or. (irow .gt. nrows0) .or. &
386 (jcol .le. 0_i4) .or. (jcol .gt. ncols0)) then
387 is_outlet = .true.
388 else
389 if (fdir0(irow, jcol) .le. 0) is_outlet = .true.
390 end if
391 !
392 if (is_outlet) then
393 noutlet = noutlet + 1_i4
394 ! cell is an outlet
395 if (noutlet .eq. 1) then
396 oloc(1, :) = level0_idomain%CellCoor(ii, :)
397 else
398 call append(oloc, level0_idomain%CellCoor(ii : ii, :))
399 end if
400 ! drain this cell into corresponding L11 cell
401 kk = l0_l11_remap(idomain)%lowres_id_on_highres(oloc(noutlet, 1), oloc(noutlet, 2))
402 drasc0(oloc(noutlet, 1), oloc(noutlet, 2)) = kk
403 ! check whether cell has maximum flow accumulation
404 ! coord. of all corners
405 iu = l0_l11_remap(idomain)%upper_bound (kk)
406 id = l0_l11_remap(idomain)%lower_bound (kk)
407 jl = l0_l11_remap(idomain)%left_bound (kk)
408 jr = l0_l11_remap(idomain)%right_bound(kk)
409 if (maxval(facc0(iu : id, jl : jr)) .eq. facc0(oloc(noutlet, 1), oloc(noutlet, 2))) then
410 ! set location of outlet at L11
411 rowout(kk) = oloc(noutlet, 1)
412 colout(kk) = oloc(noutlet, 2)
413 fdir11(level11(idomain)%CellCoor(kk, 1), level11(idomain)%CellCoor(kk, 2)) = 0
414 end if
415 end if
416 end do
417
418 ! finding cell L11 outlets - using L0_fAcc
419
420 do kk = 1, nnodes
421
422 ! exclude outlet L11
423 if (rowout(kk) > 0) cycle
424
425 ic = level11(idomain)%CellCoor(kk, 1)
426 jc = level11(idomain)%CellCoor(kk, 2)
427
428 ! coord. of all corners
429 iu = l0_l11_remap(idomain)%upper_bound (kk)
430 id = l0_l11_remap(idomain)%lower_bound (kk)
431 jl = l0_l11_remap(idomain)%left_bound (kk)
432 jr = l0_l11_remap(idomain)%right_bound(kk)
433
434 faccmax = -9
435 idmax = 0
436 side = -1
437 ! searching on side 4
438 do jj = jl, jr
439 if ((facc0(iu, jj) > faccmax) .and. &
440 (fdir0(iu, jj) == 32 .or. &
441 fdir0(iu, jj) == 64 .or. &
442 fdir0(iu, jj) == 128)) then
443 faccmax = facc0(iu, jj)
444 idmax = id0(iu, jj)
445 side = 4
446 end if
447 end do
448
449 ! searching on side 1
450 do ii = iu, id
451 if ((facc0(ii, jr) > faccmax) .and. &
452 (fdir0(ii, jr) == 1 .or. &
453 fdir0(ii, jr) == 2 .or. &
454 fdir0(ii, jr) == 128)) then
455 faccmax = facc0(ii, jr)
456 idmax = id0(ii, jr)
457 side = 1
458 end if
459 end do
460
461 ! searching on side 2
462 do jj = jl, jr
463 if ((facc0(id, jj) > faccmax) .and. &
464 (fdir0(id, jj) == 2 .or. &
465 fdir0(id, jj) == 4 .or. &
466 fdir0(id, jj) == 8)) then
467 faccmax = facc0(id, jj)
468 idmax = id0(id, jj)
469 side = 2
470 end if
471 end do
472
473 ! searching on side 3
474 do ii = iu, id
475 if ((facc0(ii, jl) > faccmax) .and. &
476 (fdir0(ii, jl) == 8 .or. &
477 fdir0(ii, jl) == 16 .or. &
478 fdir0(ii, jl) == 32)) then
479 faccmax = facc0(ii, jl)
480 idmax = id0(ii, jl)
481 side = 3
482 end if
483 end do
484
485 ! set location of the cell-outlet (row, col) in L0
486 ii = level0_idomain%CellCoor(idmax, 1)
487 jj = level0_idomain%CellCoor(idmax, 2)
488 rowout(kk) = ii
489 colout(kk) = jj
490 drasc0(ii, jj) = kk
491
492 ! set fDir at L11
493 if (ii == iu .and. jj == jl) then
494 select case (fdir0(ii, jj))
495 case (8, 16)
496 fdir11(ic, jc) = 16
497 case (32)
498 fdir11(ic, jc) = 32
499 case (64, 128)
500 fdir11(ic, jc) = 64
501 end select
502 elseif (ii == iu .and. jj == jr) then
503 select case (fdir0(ii, jj))
504 case (32, 64)
505 fdir11(ic, jc) = 64
506 case (128)
507 fdir11(ic, jc) = 128
508 case (1, 2)
509 fdir11(ic, jc) = 1
510 end select
511 elseif (ii == id .and. jj == jl) then
512 select case (fdir0(ii, jj))
513 case (2, 4)
514 fdir11(ic, jc) = 4
515 case (8)
516 fdir11(ic, jc) = 8
517 case (16, 32)
518 fdir11(ic, jc) = 16
519 end select
520 elseif (ii == id .and. jj == jr) then
521 select case (fdir0(ii, jj))
522 case (128, 1)
523 fdir11(ic, jc) = 1
524 case (2)
525 fdir11(ic, jc) = 2
526 case (4, 8)
527 fdir11(ic, jc) = 4
528 end select
529 else
530 ! cell on one side
531 select case (side)
532 case (1)
533 fdir11(ic, jc) = 1
534 case (2)
535 fdir11(ic, jc) = 4
536 case (3)
537 fdir11(ic, jc) = 16
538 case (4)
539 fdir11(ic, jc) = 64
540 case default
541 call error_message('Error L11_flow_direction: side = -1')
542 end select
543 end if
544
545 end do
546
547 END IF
548 !--------------------------------------------------------
549 ! Start padding up local variables to global variables
550 !--------------------------------------------------------
551
552 ! allocate space for row and col Outlet
553 allocate(domain_mrm(idomain)%L0_rowOutlet(1))
554 allocate(domain_mrm(idomain)%L0_colOutlet(1))
555 domain_mrm(idomain)%L0_Noutlet = nodata_i4
556 domain_mrm(idomain)%L0_rowOutlet = nodata_i4
557 domain_mrm(idomain)%L0_colOutlet = nodata_i4
558
559 ! L0 data sets
560 ! check whether L0 data is shared
561 !ToDo: check if change is correct
562 if (idomain .eq. 1) then
563 call append(l0_drasc, pack(drasc0(:, :), level0_idomain%mask))
564 else if (domainmeta%L0DataFrom(idomain) == idomain) then
565 call append(l0_drasc, pack(drasc0(:, :), level0_idomain%mask))
566 end if
567
568 domain_mrm(idomain)%L0_Noutlet = noutlet
569 ! set L0 outlet coordinates
570 old_noutlet = size(domain_mrm(idomain)%L0_rowOutlet, dim = 1)
571 if (noutlet .le. old_noutlet) then
572 domain_mrm(idomain)%L0_rowOutlet(: noutlet) = oloc(:, 1)
573 domain_mrm(idomain)%L0_colOutlet(: noutlet) = oloc(:, 2)
574 else
575 ! store up to size of old_Noutlet
576 domain_mrm(idomain)%L0_rowOutlet(: old_noutlet) = oloc(: old_noutlet, 1)
577 domain_mrm(idomain)%L0_colOutlet(: old_noutlet) = oloc(: old_noutlet, 2)
578 ! enlarge rowOutlet and colOutlet in domain_mrm structure
579 !TODO: do other domains also need to be enlarged accordingly???
580 call append(domain_mrm(idomain)%L0_rowOutlet, oloc(old_noutlet + 1 :, 1))
581 call append(domain_mrm(idomain)%L0_colOutlet, oloc(old_noutlet + 1 :, 2))
582 end if
583
584 ! L11 data sets
585 call append(l11_noutlets, count(fdir11 .eq. 0_i4))
586 call append(l11_fdir, pack(fdir11(:, :), level11(idomain)%mask))
587 call append(l11_rowout, rowout(:))
588 call append(l11_colout, colout(:))
589
590 ! communicate
591 call message('')
592 call message(' Domain: ' // num2str(idomain, '(i3)'))
593 call message(' Number of outlets found at Level 0:.. ' // num2str(noutlet, '(i7)'))
594 call message(' Number of outlets found at Level 11:. ' // num2str(count(fdir11 .eq. 0_i4), '(i7)'))
595
596 ! free space
597 deallocate(fdir0, facc0, fdir11, rowout, colout, drasc0)
598
599 end subroutine l11_flow_direction
600
601 ! ------------------------------------------------------------------
602
603 ! NAME
604 ! L11_set_network_topology
605
606 ! PURPOSE
607 !> \brief Set network topology
608
609 !> \details Set network topology from and to node for all links
610 !> at level-11 (\ref fig_routing "Routing Network").
611 !> If a variable is added or removed here, then it also has to
612 !> be added or removed in the subroutine L11_config_set in
613 !> module mo_restart and in the subroutine set_L11_config in module
614 !> mo_set_netcdf_restart.
615
616 ! INTENT(IN)
617 !> \param[in] "integer(i4) :: iDomain" Domain Id
618
619 ! HISTORY
620 !> \authors Luis Samaniego
621
622 !> \date Dec 2005
623
624 ! Modifications:
625 ! Luis Samaniego Jan 2013 - modular version
626 ! Stephan Thober Aug 2015 - ported to mRM
627 ! Stephan Thober May 2016 - moved calculation of sink here
628 ! Robert Schweppe Jun 2018 - refactoring and reformatting
629
630 subroutine l11_set_network_topology(iDomain)
631
632 use mo_append, only : append
635
636 implicit none
637
638 ! Domain Id
639 integer(i4), intent(in) :: idomain
640
641 integer(i4), dimension(:, :), allocatable :: fdir11
642
643 integer(i4), dimension(:, :), allocatable :: dummy_2d_id
644
645 integer(i4) :: jj, kk, ic, jc
646
647 integer(i4) :: fn, tn
648
649 integer(i4), dimension(:), allocatable :: nlinkfromn, nlinkton
650
651
652 ! Routing network vectors have nNodes size instead of nLinks to
653 ! avoid the need of having two extra indices to identify a Domain.
654 ! allocate
655 allocate (nlinkfromn(level11(idomain)%nCells)) ! valid from (1 : nLinks)
656 allocate (nlinkton(level11(idomain)%nCells)) ! "
657 allocate (fdir11(level11(idomain)%nrows, level11(idomain)%ncols))
658 allocate (dummy_2d_id(level11(idomain)%nrows, level11(idomain)%ncols))
659 dummy_2d_id = unpack(level11(idomain)%Id, level11(idomain)%mask, nodata_i4)
660
661
662 ! initialize
663 nlinkfromn(:) = nodata_i4
664 nlinkton(:) = nodata_i4
665 fdir11(:, :) = nodata_i4
666
667 ! get grids of L11
668 fdir11(:, :) = unpack(l11_fdir(level11(idomain)%iStart : level11(idomain)%iEnd), level11(idomain)%mask, nodata_i4)
669
670 ! ------------------------------------------------------------------
671 ! network topology
672 ! ------------------------------------------------------------------
673
674 jj = 0
675 do kk = 1, level11(idomain)%nCells
676 ic = level11(idomain)%CellCoor(kk, 1)
677 jc = level11(idomain)%CellCoor(kk, 2)
678 fn = kk
679 call movedownonecell(fdir11(ic, jc), ic, jc)
680 tn = dummy_2d_id(ic, jc)
681 if (fn == tn) cycle
682 jj = jj + 1
683 nlinkfromn(jj) = fn
684 nlinkton(jj) = tn
685 end do
686
687 !--------------------------------------------------------
688 ! Start padding up local variables to global variables
689 !--------------------------------------------------------
690
691 ! L11 data sets
692 call append(l11_fromn, nlinkfromn(:)) ! sinks are at the end
693 call append(l11_ton, nlinkton(:))
694
695 ! free space
696 deallocate (fdir11, nlinkfromn, nlinkton)
697
698 end subroutine l11_set_network_topology
699
700 ! ------------------------------------------------------------------
701
702 ! NAME
703 ! L11_routing_order
704
705 ! PURPOSE
706 !> \brief Find routing order, headwater cells and sink
707
708 !> \details Find routing order, headwater cells and sink.
709 !> If a variable is added or removed here, then it also has to
710 !> be added or removed in the subroutine L11_config_set in
711 !> module mo_restart and in the subroutine set_L11_config in module
712 !> mo_set_netcdf_restart
713
714 ! INTENT(IN)
715 !> \param[in] "integer(i4) :: iDomain" Domain Id
716
717 ! HISTORY
718 !> \authors Luis Samaniego
719
720 !> \date Dec 2005
721
722 ! Modifications:
723 ! Luis Samaniego Jan 2013 - modular version
724 ! Sa. Ku. Jan 2015 - corrected initialization of nLinkSink
725 ! Stephan Thober Aug 2015 - ported to mRM
726 ! Robert Schweppe Jun 2018 - refactoring and reformatting
727
728 subroutine l11_routing_order(iDomain)
729
730 use mo_append, only : append
734
735 implicit none
736
737 ! Domain Id
738 integer(i4), intent(in) :: idomain
739
740 integer(i4) :: nlinks
741
742 ! from node
743 integer(i4), dimension(:), allocatable :: nlinkfromn
744
745 ! to node
746 integer(i4), dimension(:), allocatable :: nlinkton
747
748 ! network routing order
749 integer(i4), dimension(:), allocatable :: nlinkrorder
750
751 ! label Id [0='', 1=HeadWater, 2=Sink]
752 integer(i4), dimension(:), allocatable :: nlinklabel
753
754 ! == .true. if sink node reached
755 logical, dimension(:), allocatable :: nlinksink
756
757 ! routing order (permutation)
758 integer(i4), dimension(:), allocatable :: netperm
759
760 integer(i4) :: ii, jj, kk
761
762 logical :: flag
763
764
765 nlinks = level11(idomain)%nCells - l11_noutlets(idomain)
766 ! Routing network vectors have nNodes size instead of nLinks to
767 ! avoid the need of having two extra indices to identify a Domain.
768
769 ! allocate
770 allocate (nlinkfromn(level11(idomain)%nCells)) ! all vectors valid from (1 : nLinks)
771 allocate (nlinkton(level11(idomain)%nCells))
772 allocate (nlinkrorder(level11(idomain)%nCells))
773 allocate (nlinklabel(level11(idomain)%nCells))
774 allocate (nlinksink(level11(idomain)%nCells))
775 allocate (netperm(level11(idomain)%nCells))
776 ! initialize
777 nlinkfromn(:) = nodata_i4
778 nlinkton(:) = nodata_i4
779 nlinkrorder(1 : nlinks) = 1
780 nlinkrorder(level11(idomain)%nCells) = nodata_i4
781 nlinklabel(1 : nlinks) = 0
782 nlinklabel(level11(idomain)%nCells) = nodata_i4
783 nlinksink(:) = .false.
784 netperm(:) = nodata_i4
785
786 ! for a single node model run
787 if(level11(idomain)%nCells .GT. 1) then
788 ! get network vectors of L11
789 nlinkfromn(:) = l11_fromn(level11(idomain)%iStart : level11(idomain)%iEnd)
790 nlinkton(:) = l11_ton(level11(idomain)%iStart : level11(idomain)%iEnd)
791
792 loop1 : do ii = 1, nlinks
793 loop2 : do jj = 1, nlinks
794 if (jj == ii) cycle loop2
795 if (nlinkfromn(ii) == nlinkton(jj)) then
796 nlinkrorder(ii) = -9
797 end if
798 if (nlinkrorder(ii) == -9) cycle loop1
799 end do loop2
800 end do loop1
801 ! counting headwaters
802 kk = 0
803 do ii = 1, nlinks
804 if (nlinkrorder(ii) == 1) then
805 kk = kk + 1
806 nlinkrorder(ii) = kk
807 nlinklabel(ii) = 1 ! 'Head Water'
808 end if
809 end do
810 ! counting downstream
811 do while (minval(nlinkrorder(1 : nlinks)) < 0)
812 !! print *, count(nLinkROrder .lt. 0), minval(nLinkROrder)
813 loop3 : do ii = 1, nlinks
814 if (.NOT. nlinkrorder(ii) == -9) cycle loop3
815 flag = .true.
816 loop4 : do jj = 1, nlinks
817 if (jj == ii .OR. nlinkfromn(ii) /= nlinkton(jj)) then
818 cycle loop4
819 else if (.NOT. (nlinkfromn(ii) == nlinkton(jj) .AND. nlinkrorder(jj) > 0)) then
820 flag = .false.
821 exit loop4
822 else
823 end if
824 end do loop4
825
826 if (flag) then
827 kk = kk + 1
828 nlinkrorder(ii) = kk
829 end if
830 end do loop3
831 end do
832
833 ! identify sink cells
834 do ii = 1, nlinks
835 if (l11_fdir(level11(idomain)%iStart + nlinkton(ii) - 1_i4) .eq. 0_i4) then
836 nlinksink(ii) = .true.
837 ! add sink target cell to sink_cells if not already present (in case of multiple inflows (rare case))
838 if (.not.any(sink_cells(idomain)%ids == nlinkton(ii))) sink_cells(idomain)%ids = [sink_cells(idomain)%ids, nlinkton(ii)]
839 end if
840 end do
841 where(nlinksink) nlinklabel = 2 ! 'Sink'
842
843 ! keep routing order
844 do ii = 1, nlinks
845 netperm(nlinkrorder(ii)) = ii
846 end do
847
848 ! end of multi-node network design loop
849 end if
850
851 !--------------------------------------------------------
852 ! Start padding up local variables to global variables
853 !--------------------------------------------------------
854 ! L11 network data sets
855 call append(l11_rorder, nlinkrorder(:))
856 call append(l11_label, nlinklabel(:))
857 call append(l11_sink, nlinksink(:))
858 call append(l11_netperm, netperm(:))
859
860 ! free space
861 deallocate (nlinkfromn, nlinkton, nlinkrorder, nlinklabel, nlinksink, netperm)
862
863 end subroutine l11_routing_order
864
865 ! ------------------------------------------------------------------
866
867 ! NAME
868 ! L11_link_location
869
870 ! PURPOSE
871 !> \brief Estimate the LO (row,col) location for each routing link at level L11
872
873 !> \details If a variable is added or removed here, then it also has to
874 !> be added or removed in the subroutine L11_config_set in
875 !> module mo_restart and in the subroutine set_L11_config in module
876 !> mo_set_netcdf_restart
877
878 ! INTENT(IN)
879 !> \param[in] "integer(i4) :: iDomain" Domain Id
880
881 ! HISTORY
882 !> \authors Luis Samaniego
883
884 !> \date Dec 2005
885
886 ! Modifications:
887 ! Luis Samaniego Jan 2013 - modular version
888 ! Stephan Thober Aug 2015 - ported to mRM
889 ! Robert Schweppe Jun 2018 - refactoring and reformatting
890
891 subroutine l11_link_location(iDomain)
892
893 use mo_append, only : append
895 use mo_common_types, only: grid
899 use mo_string_utils, only : num2str
900
901 implicit none
902
903 ! Domain Id
904 integer(i4), intent(in) :: idomain
905
906 integer(i4) :: nlinks
907
908 ! northing cell loc. of the Outlet
909 integer(i4), dimension(:), allocatable :: rowout
910
911 ! easting cell loc. of the Outlet
912 integer(i4), dimension(:), allocatable :: colout
913
914 integer(i4), dimension(:), allocatable :: nlinkfromn
915
916 integer(i4), dimension(:), allocatable :: netperm
917
918 integer(i4), dimension(:), allocatable :: nlinkfromrow
919
920 integer(i4), dimension(:), allocatable :: nlinkfromcol
921
922 integer(i4), dimension(:), allocatable :: nlinktorow
923
924 integer(i4), dimension(:), allocatable :: nlinktocol
925
926 integer(i4), dimension(:, :), allocatable :: fdir0
927
928 integer(i4), dimension(:, :), allocatable :: drasc0
929
930 integer(i4) :: ii, rr, kk, s0, e0
931
932 integer(i4) :: inode, irow, jcol, prevrow, prevcol
933
934 ! output location in L0
935 integer(i4), dimension(:, :), allocatable :: oloc
936
937 ! number of outlets in Domain
938 integer(i4) :: noutlets
939
940 ! flag for finding outlet
941 logical :: is_outlet
942
943 type(grid), pointer :: level0_idomain => null()
944
945
946 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
947 s0 = level0_idomain%iStart
948 e0 = level0_idomain%iEnd
949
950 noutlets = l11_noutlets(idomain)
951
952 nlinks = level11(idomain)%nCells - noutlets
953
954 ! Routing network vectors have level11(iDomain)%nCells size instead of nLinks to
955 ! avoid the need of having two extra indices to identify a Domain.
956 ! allocate
957 allocate (rowout(level11(idomain)%nCells))
958 allocate (colout(level11(idomain)%nCells))
959 allocate (nlinkfromn(level11(idomain)%nCells)) ! all network vectors valid from (1 : nLinks)
960 allocate (netperm(level11(idomain)%nCells))
961 allocate (nlinkfromrow(level11(idomain)%nCells))
962 allocate (nlinkfromcol(level11(idomain)%nCells))
963 allocate (nlinktorow(level11(idomain)%nCells))
964 allocate (nlinktocol(level11(idomain)%nCells))
965 allocate (fdir0(level0_idomain%nrows, level0_idomain%ncols))
966 allocate (drasc0(level0_idomain%nrows, level0_idomain%ncols))
967
968 ! initialize
969 rowout = nodata_i4
970 colout = nodata_i4
971 nlinkfromn = nodata_i4
972 netperm = nodata_i4
973 nlinkfromrow = nodata_i4
974 nlinkfromcol = nodata_i4
975 nlinktorow = nodata_i4
976 nlinktocol = nodata_i4
977 fdir0 = nodata_i4
978 drasc0 = nodata_i4
979
980 ! for a single node model run
981 if(level11(idomain)%nCells .GT. 1) then
982 ! get fDir at L0
983 fdir0(:, :) = unpack(l0_fdir(s0 : e0), level0_idomain%mask, nodata_i4)
984 drasc0(:, :) = unpack(l0_drasc(s0 : e0), level0_idomain%mask, nodata_i4)
985
986 ! get network vectors of L11
987 nlinkfromn(:) = l11_fromn(level11(idomain)%iStart : level11(idomain)%iEnd)
988 netperm(:) = l11_netperm(level11(idomain)%iStart : level11(idomain)%iEnd)
989 rowout(:) = l11_rowout(level11(idomain)%iStart : level11(idomain)%iEnd)
990 colout(:) = l11_colout(level11(idomain)%iStart : level11(idomain)%iEnd)
991
992 ! finding main outlet (row, col) in L0
993 allocate(oloc(noutlets, 2))
994 oloc(:, 1) = domain_mrm(idomain)%L0_rowOutlet(: noutlets)
995 oloc(:, 2) = domain_mrm(idomain)%L0_colOutlet(: noutlets)
996
997 ! Location of the stream-joint cells (row, col)
998 do rr = 1, nlinks
999
1000 ii = netperm(rr)
1001 inode = nlinkfromn(ii)
1002 irow = rowout(inode)
1003 jcol = colout(inode)
1004 call movedownonecell(fdir0(irow, jcol), irow, jcol)
1005 ! set "from" cell
1006 nlinkfromrow(ii) = irow
1007 nlinkfromcol(ii) = jcol
1008
1009 ! check whether this location is an outlet
1010 is_outlet = .false.
1011 do kk = 1, noutlets
1012 if (irow .eq. oloc(kk, 1) .and. jcol .eq. oloc(kk, 2)) is_outlet = .true.
1013 end do
1014
1015 if (is_outlet) then
1016
1017 nlinktorow(ii) = irow
1018 nlinktocol(ii) = jcol
1019
1020 else
1021
1022 do while (.not. (drasc0(irow, jcol) > 0))
1023 prevrow = irow
1024 prevcol = jcol
1025 call movedownonecell(fdir0(irow, jcol), irow, jcol)
1026 ! check whether this location is an outlet and exit
1027 do kk = 1, noutlets
1028 if (irow .eq. oloc(kk, 1) .and. jcol .eq. oloc(kk, 2)) exit
1029 end do
1030 if (prevrow .eq. irow .and. prevcol .eq. jcol) then
1031 call error_message('Something went wrong during L11_link_location, ', &
1032 'movedownonecell got stuck in infinite loop at cell (', num2str(irow), ' ', &
1033 num2str(jcol))
1034 end if
1035 end do
1036 ! set "to" cell (when an outlet is reached)
1037 nlinktorow(ii) = irow
1038 nlinktocol(ii) = jcol
1039
1040 end if
1041 end do
1042
1043 ! end of multi-node network design loop
1044 end if
1045
1046 !--------------------------------------------------------
1047 ! Start padding up local variables to global variables
1048 !--------------------------------------------------------
1049 ! L11 network data sets
1050 call append(l11_frow, nlinkfromrow(:))
1051 call append(l11_fcol, nlinkfromcol(:))
1052 call append(l11_trow, nlinktorow(:))
1053 call append(l11_tcol, nlinktocol(:))
1054
1055 ! free space
1056 deallocate (rowout, colout, nlinkfromn, netperm, nlinkfromrow, &
1057 nlinkfromcol, nlinktorow, nlinktocol, fdir0, drasc0)
1058
1059 end subroutine l11_link_location
1060
1061 ! ------------------------------------------------------------------
1062
1063 ! NAME
1064 ! L11_set_drain_outlet_gauges
1065
1066 ! PURPOSE
1067 !> \brief Draining cell identification and Set gauging node
1068
1069 !> \details Perform the following tasks:
1070 !> - Draining cell identification (cell at L0 to draining cell outlet at L11).
1071 !> - Set gauging nodes
1072 !> If a variable is added or removed here, then it also has to
1073 !> be added or removed in the subroutine L11_config_set in
1074 !> module mo_restart and in the subroutine set_L11_config in module
1075 !> mo_set_netcdf_restart
1076
1077 ! INTENT(IN)
1078 !> \param[in] "integer(i4) :: iDomain" Domain Id
1079
1080 ! HISTORY
1081 !> \authors Luis Samaniego
1082
1083 !> \date Dec 2005
1084
1085 ! Modifications:
1086 ! Luis Samaniego Jan 2013 - modular version
1087 ! Matthias Zink Mar 2014 - bugfix, added inflow gauge
1088 ! Rohini Kumar Apr 2014 - variable index is changed to index_gauge
1089 ! Stephan Thober Aug 2015 - ported to mRM
1090 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1091
1093
1094 use mo_append, only : append
1095 use mo_common_constants, only : nodata_i4
1096 use mo_common_types, only: grid
1100
1101 implicit none
1102
1103 ! Domain Id
1104 integer(i4), intent(in) :: idomain
1105
1106 integer(i4), dimension(:, :), allocatable :: drasc0
1107
1108 integer(i4), dimension(:, :), allocatable :: fdir0
1109
1110 integer(i4), dimension(:, :), allocatable :: gaugeloc0
1111
1112 integer(i4), dimension(:, :), allocatable :: inflowgaugeloc0
1113
1114 integer(i4), dimension(:, :), allocatable :: dracell0
1115
1116 integer(i4) :: ii, jj, kk, ll, s0, e0
1117
1118 integer(i4) :: isc
1119
1120 integer(i4) :: irow, jcol
1121
1122 type(grid), pointer :: level0_idomain => null()
1123
1124
1125 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
1126 s0 = level0_idomain%iStart
1127 e0 = level0_idomain%iEnd
1128
1129
1130 ! allocate
1131 allocate (drasc0(level0_idomain%nrows, level0_idomain%ncols))
1132 allocate (fdir0(level0_idomain%nrows, level0_idomain%ncols))
1133 allocate (gaugeloc0(level0_idomain%nrows, level0_idomain%ncols))
1134 allocate (inflowgaugeloc0(level0_idomain%nrows, level0_idomain%ncols))
1135 allocate (dracell0(level0_idomain%nrows, level0_idomain%ncols))
1136
1137 ! initialize
1138 drasc0(:, :) = nodata_i4
1139 fdir0(:, :) = nodata_i4
1140 gaugeloc0(:, :) = nodata_i4
1141 inflowgaugeloc0(:, :) = nodata_i4
1142 dracell0(:, :) = nodata_i4
1143
1144 drasc0(:, :) = unpack(l0_drasc(s0 : e0), &
1145 level0_idomain%mask, nodata_i4)
1146 fdir0(:, :) = unpack(l0_fdir(s0 : e0), &
1147 level0_idomain%mask, nodata_i4)
1148 gaugeloc0(:, :) = unpack(l0_gaugeloc(s0 : e0), &
1149 level0_idomain%mask, nodata_i4)
1150 inflowgaugeloc0(:, :) = unpack(l0_inflowgaugeloc(s0 : e0), &
1151 level0_idomain%mask, nodata_i4)
1152
1153 do kk = 1, level0_idomain%nCells
1154 ii = level0_idomain%CellCoor(kk, 1)
1155 jj = level0_idomain%CellCoor(kk, 2)
1156 isc = drasc0(ii, jj)
1157 ! find drainage path
1158 irow = ii
1159 jcol = jj
1160 do while (.NOT. isc > 0)
1161 ! move downstream
1162 call movedownonecell(fdir0(irow, jcol), irow, jcol)
1163 isc = drasc0(irow, jcol)
1164 end do
1165 dracell0(ii, jj) = isc
1166
1167 ! find cell at L11 corresponding to gauges in Domain at L0 ! >> L11_on_L0 is Id of
1168 ! the routing cell at level-11
1169 if (gaugeloc0(ii, jj) .NE. nodata_i4) then
1170 ! evaluation gauges
1171 do ll = 1, domain_mrm(idomain)%nGauges
1172 ! search for gaugeID in L0 grid and save ID on L11
1173 if (domain_mrm(idomain)%gaugeIdList(ll) .EQ. gaugeloc0(ii, jj)) then
1174 domain_mrm(idomain)%gaugeNodeList(ll) = l0_l11_remap(idomain)%lowres_id_on_highres(ii, jj)
1175 end if
1176 end do
1177 end if
1178
1179 if (inflowgaugeloc0(ii, jj) .NE. nodata_i4) then
1180 ! inflow gauges
1181 do ll = 1, domain_mrm(idomain)%nInflowGauges
1182 ! search for gaugeID in L0 grid and save ID on L11
1183 if (domain_mrm(idomain)%InflowGaugeIdList(ll) .EQ. inflowgaugeloc0(ii, jj)) &
1184 domain_mrm(idomain)%InflowGaugeNodeList(ll) = l0_l11_remap(idomain)%lowres_id_on_highres(ii, jj)
1185 end do
1186 end if
1187 end do
1188
1189 !--------------------------------------------------------
1190 ! Start padding up local variables to global variables
1191 !--------------------------------------------------------
1192 ! L0 data sets
1193 ! check whether L0 data is shared
1194 ! ToDo: check if change is correct
1195 if (idomain .eq. 1) then
1196 call append(l0_dracell, pack(dracell0(:, :), level0_idomain%mask))
1197 else if (domainmeta%L0DataFrom(idomain) == idomain) then
1198 call append(l0_dracell, pack(dracell0(:, :), level0_idomain%mask))
1199 end if
1200
1201 ! free space
1202 deallocate (drasc0, fdir0, gaugeloc0, dracell0)
1203
1204 end subroutine l11_set_drain_outlet_gauges
1205
1206 ! ------------------------------------------------------------------
1207
1208 ! NAME
1209 ! L11_stream_features
1210
1211 ! PURPOSE
1212 !> \brief Stream features (stream network and floodplain)
1213
1214 !> \details Stream features (stream network and floodplain)
1215 !> If a variable is added or removed here, then it also has to
1216 !> be added or removed in the subroutine L11_config_set in
1217 !> module mo_restart and in the subroutine set_L11_config in module
1218 !> mo_set_netcdf_restart
1219
1220 ! INTENT(IN)
1221 !> \param[in] "integer(i4) :: iDomain" Domain Id
1222
1223 ! HISTORY
1224 !> \authors Luis Samaniego
1225
1226 !> \date Dec 2005
1227
1228 ! Modifications:
1229 ! Luis Samaniego Jan 2013 - modular version
1230 ! R. Kumar Oct 2013 - stack size increased from nNodes to 100
1231 ! Stephan Thober Aug 2015 - ported to mRM
1232 ! Stephan Thober Nov 2016 - only read flood plain area if processMatrix for routing equals 1
1233 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1234 ! Stephan Thober Jul 2018 - introduced cut off Length at 40 percentile to neglect short paths in headwaters for adaptive timesteps
1235 ! Stephan Thober, Pallav Kumar Shrestha, Sep 2020 - bug fix in cut off Length at 40 percentile, neglecting links with -9999. that occur if multiple outlets are present
1236
1237 subroutine l11_stream_features(iDomain)
1238
1239 use mo_append, only : append
1241 use mo_common_types, only: grid
1243 use mo_mrm_global_variables, only : l0_fdir, &
1246 use mo_percentile, only: percentile
1247
1248 implicit none
1249
1250 ! Domain Id
1251 integer(i4), intent(in) :: idomain
1252
1253 integer(i4) :: nlinks
1254
1255 integer(i4), dimension(:, :), allocatable :: id0
1256
1257 integer(i4), dimension(:, :), allocatable :: fdir0
1258
1259 real(dp), dimension(:, :), allocatable :: elev0
1260
1261 real(dp), dimension(:, :), allocatable :: cellarea0
1262
1263 integer(i4), dimension(:, :), allocatable :: streamnet0
1264
1265 integer(i4), dimension(:, :), allocatable :: floodplain0
1266
1267 ! routing order (permutation)
1268 integer(i4), dimension(:), allocatable :: netperm
1269
1270 integer(i4), dimension(:), allocatable :: nlinkfromrow
1271
1272 integer(i4), dimension(:), allocatable :: nlinkfromcol
1273
1274 integer(i4), dimension(:), allocatable :: nlinktorow
1275
1276 integer(i4), dimension(:), allocatable :: nlinktocol
1277
1278 real(dp), dimension(:), allocatable :: nlinklength
1279
1280 real(dp), dimension(:), allocatable :: nlinkafloodplain
1281
1282 real(dp), dimension(:), allocatable :: nlinkslope
1283
1284 integer(i4) :: ii, rr, ns, s0, e0
1285
1286 integer(i4) :: frow, fcol
1287
1288 integer(i4) :: fid, tid
1289
1290 integer(i4), dimension(:, :), allocatable :: stack, append_chunk
1291
1292 integer(i4), dimension(:), allocatable :: dummy_1d
1293
1294 real(dp) :: length
1295
1296 integer(i4), dimension(:, :), allocatable :: nodata_i4_tmp
1297
1298 real(dp), dimension(:, :), allocatable :: nodata_dp_tmp
1299
1300 type(grid), pointer :: level0_idomain => null()
1301
1302
1303 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
1304 s0 = level0_idomain%iStart
1305 e0 = level0_idomain%iEnd
1306 nlinks = level11(idomain)%nCells - l11_noutlets(idomain)
1307
1308
1309 ! allocate
1310 allocate (id0(level0_idomain%nrows, level0_idomain%ncols))
1311 allocate (elev0(level0_idomain%nrows, level0_idomain%ncols))
1312 allocate (fdir0(level0_idomain%nrows, level0_idomain%ncols))
1313 allocate (cellarea0(level0_idomain%nrows, level0_idomain%ncols))
1314 allocate (streamnet0(level0_idomain%nrows, level0_idomain%ncols))
1315 allocate (floodplain0(level0_idomain%nrows, level0_idomain%ncols))
1316
1317 ! Routing network vectors have nNodes size instead of nLinks to
1318 ! avoid the need of having two extra indices to identify a Domain.
1319 allocate (stack(level11(idomain)%nCells, 2)) ! >> stack(nNodes, 2)
1320 allocate (dummy_1d(2))
1321 allocate (append_chunk(8, 2))
1322 allocate (netperm(level11(idomain)%nCells))
1323 allocate (nlinkfromrow(level11(idomain)%nCells))
1324 allocate (nlinkfromcol(level11(idomain)%nCells))
1325 allocate (nlinktorow(level11(idomain)%nCells))
1326 allocate (nlinktocol(level11(idomain)%nCells))
1327 allocate (nlinklength(level11(idomain)%nCells))
1328 allocate (nlinkafloodplain(level11(idomain)%nCells))
1329 allocate (nlinkslope(level11(idomain)%nCells))
1330
1331 allocate (nodata_i4_tmp(level0_idomain%nrows, level0_idomain%ncols))
1332 allocate (nodata_dp_tmp(level0_idomain%nrows, level0_idomain%ncols))
1333
1334 ! initialize
1335 id0(:, :) = nodata_i4
1336 elev0(:, :) = nodata_dp
1337 fdir0(:, :) = nodata_i4
1338 cellarea0(:, :) = nodata_dp
1339 streamnet0(:, :) = nodata_i4
1340 floodplain0(:, :) = nodata_i4
1341
1342 stack(:, :) = nodata_i4
1343 append_chunk(:, :) = nodata_i4
1344 netperm(:) = nodata_i4
1345 nlinkfromrow(:) = nodata_i4
1346 nlinkfromcol(:) = nodata_i4
1347 nlinktorow(:) = nodata_i4
1348 nlinktocol(:) = nodata_i4
1349 nlinklength(:) = nodata_dp
1350 nlinkafloodplain(:) = nodata_dp
1351 nlinkslope(:) = nodata_dp
1352
1353 nodata_i4_tmp(:, :) = nodata_i4
1354 nodata_dp_tmp(:, :) = nodata_dp
1355
1356 ! for a single node model run
1357 if(level11(idomain)%nCells .GT. 1) then
1358 ! get L0 fields
1359 id0(:, :) = unpack(level0_idomain%Id, level0_idomain%mask, nodata_i4_tmp)
1360 elev0(:, :) = unpack(l0_elev(s0 : e0), &
1361 level0_idomain%mask, nodata_dp_tmp)
1362 fdir0(:, :) = unpack(l0_fdir(s0 : e0), &
1363 level0_idomain%mask, nodata_i4_tmp)
1364 cellarea0(:, :) = unpack(level0_idomain%CellArea, level0_idomain%mask, nodata_dp_tmp)
1365
1366 ! get network vectors of L11
1367 netperm(:) = l11_netperm(level11(idomain)%iStart : level11(idomain)%iEnd)
1368 nlinkfromrow(:) = l11_frow(level11(idomain)%iStart : level11(idomain)%iEnd)
1369 nlinkfromcol(:) = l11_fcol(level11(idomain)%iStart : level11(idomain)%iEnd)
1370 nlinktorow(:) = l11_trow(level11(idomain)%iStart : level11(idomain)%iEnd)
1371 nlinktocol(:) = l11_tcol(level11(idomain)%iStart : level11(idomain)%iEnd)
1372
1373 ! Flood plains: stream network delineation
1374 streamnet0(:, :) = nodata_i4
1375 floodplain0(:, :) = nodata_i4
1376
1377 do rr = 1, nlinks
1378
1379 ii = netperm(rr)
1380 frow = nlinkfromrow(ii)
1381 fcol = nlinkfromcol(ii)
1382
1383 ! Init
1384 streamnet0(frow, fcol) = ii
1385 floodplain0(frow, fcol) = ii
1386 stack = 0
1387 append_chunk = 0
1388 ns = 1
1389 stack(ns, 1) = frow
1390 stack(ns, 2) = fcol
1391
1392 call celllength(idomain, fdir0(frow, fcol), frow, fcol, iflag_cordinate_sys, nlinklength(ii))
1393 nlinkslope(ii) = elev0(frow, fcol)
1394
1395 fid = id0(frow, fcol)
1396 tid = id0(nlinktorow(ii), nlinktocol(ii))
1397
1398 do while (.NOT. (fid == tid))
1399
1400 ! Search flood plain from point(frow,fcol) upwards, keep co-ordinates in STACK
1401
1402 ! flood plain calculation is expensive and only required for processCase = 1
1403 if (processmatrix(8,1) == 1) then
1404 do while (ns > 0)
1405 if (ns + 8 .gt. size(stack, 1)) then
1406 call append(stack, append_chunk)
1407 end if
1408 call moveup(elev0, fdir0, frow, fcol, stack, ns)
1409 stack(1, 1) = 0
1410 stack(1, 2) = 0
1411 ! stack = cshift(stack, SHIFT = 1, DIM = 1)
1412 ! substitute cshift <<<
1413 dummy_1d = stack(1, :)
1414 stack(: size(stack, dim = 1) - 1, :) = stack(2 :, :)
1415 stack(size(stack, dim = 1), :) = dummy_1d
1416 ! substitute cshift >>>
1417 if (stack(1, 1) > 0 .and. stack(1, 2) > 0) floodplain0(stack(1, 1), stack(1, 2)) = ii
1418 ns = count(stack > 0) / 2
1419 end do
1420 end if
1421
1422 ! move downstream
1423 call movedownonecell(fdir0(frow, fcol), frow, fcol)
1424 streamnet0(frow, fcol) = ii
1425 floodplain0(frow, fcol) = ii
1426 fid = id0(frow, fcol)
1427 stack = 0
1428 ns = 1
1429 stack(ns, 1) = frow
1430 stack(ns, 2) = fcol
1431 call celllength(idomain, fdir0(frow, fcol), frow, fcol, iflag_cordinate_sys, length)
1432 nlinklength(ii) = nlinklength(ii) + length
1433
1434 end do
1435
1436 ! stream bed slope
1437 nlinkslope(ii) = (nlinkslope(ii) - elev0(frow, fcol)) / nlinklength(ii)
1438
1439 if (nlinkslope(ii) < 0.0001_dp) nlinkslope(ii) = 0.0001_dp
1440
1441 ! calculate area of floodplains (avoid overwriting)
1442 nlinkafloodplain(ii) = sum(cellarea0(:, :), mask = (floodplain0(:, :) == ii))
1443 ! old > real( count( floodPlain0(:,:,) == i), dp ) * cellarea0
1444
1445 end do
1446
1447 ! end of multi-node network design loop
1448 end if
1449
1450 ! cut off Length at 40 percentile to neglect short paths in headwaters
1451 if ((processmatrix(8, 1) .eq. 2) .or. (processmatrix(8, 1) .eq. 3)) then
1452 if (count(nlinklength(:) .ge. 0._dp) .gt. 2) then
1453 length = percentile(pack(nlinklength(:), nlinklength(:) .ge. 0._dp), 40._dp)
1454 nlinklength(:) = merge(nlinklength(:), length, (nlinklength(:) .gt. length))
1455 end if
1456 end if
1457
1458 !--------------------------------------------------------
1459 ! Start padding up local variables to global variables
1460 !--------------------------------------------------------
1461
1462 ! L0 data sets
1463 ! check whether L0 data is shared
1464 ! ToDo: check if change is correct
1465 if (idomain .eq. 1) then
1466 call append(l0_streamnet, pack(streamnet0(:, :), level0_idomain%mask))
1467 call append(l0_floodplain, pack(floodplain0(:, :), level0_idomain%mask))
1468 else if (domainmeta%L0DataFrom(idomain) == idomain) then
1469 call append(l0_streamnet, pack(streamnet0(:, :), level0_idomain%mask))
1470 call append(l0_floodplain, pack(floodplain0(:, :), level0_idomain%mask))
1471 end if
1472
1473
1474 ! L11 network data sets
1475 call append(l11_length, nlinklength(:))
1476 call append(l11_afloodplain, nlinkafloodplain(:))
1477 call append(l11_slope, nlinkslope(:))
1478
1479 ! free space
1480 deallocate (&
1481 id0, elev0, fdir0, streamnet0, floodplain0, &
1482 cellarea0, stack, netperm, nlinkfromrow, nlinkfromcol, nlinktorow, nlinktocol, &
1483 nlinklength, nlinkafloodplain, nlinkslope, dummy_1d)
1484 deallocate(nodata_i4_tmp, nodata_dp_tmp)
1485
1486 end subroutine l11_stream_features
1487
1488 ! ------------------------------------------------------------------
1489
1490 ! NAME
1491 ! L11_fraction_sealed_floodplain
1492
1493 ! PURPOSE
1494 !> \brief Fraction of the flood plain with impervious cover
1495
1496 !> \details Fraction of the flood plain with impervious cover (\ref fig_routing "Routing
1497 !> Network"). This proportion is used to regionalize the Muskingum parameters.
1498 !> Samaniego et al. \cite SB05 found out that this fraction is one of the statistically
1499 !> significant predictor variables of peak discharge in mesoscale Domains.
1500 !> If a variable is added or removed here, then it also has to
1501 !> be added or removed in the subroutine L11_config_set in
1502 !> module mo_restart and in the subroutine set_L11_config in module
1503 !> mo_set_netcdf_restart
1504
1505 ! INTENT(IN)
1506 !> \param[in] "integer(i4) :: LCClassImp" Impervious land cover class Id, e.g. = 2 (old code)
1507 !> \param[in] "logical :: do_init"
1508
1509 ! HISTORY
1510 !> \authors Luis Samaniego
1511
1512 !> \date Dec 2005
1513
1514 ! Modifications:
1515 ! Luis Samaniego Jan 2013 - modular version
1516 ! Stephan Thober Aug 2015 - ported to mRM
1517 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1518
1519 subroutine l11_fraction_sealed_floodplain(LCClassImp, do_init)
1520
1521 use mo_append, only : append
1522 use mo_common_constants, only : nodata_dp
1523 use mo_common_types, only: grid
1527
1528 implicit none
1529
1530 ! Impervious land cover class Id, e.g. = 2 (old code)
1531 integer(i4), intent(in) :: lcclassimp
1532
1533 logical, intent(in) :: do_init
1534
1535 integer(i4) :: nlinks
1536
1537 real(dp), dimension(:), pointer :: nlinkafloodplain => null()
1538
1539 real(dp), dimension(:,:), allocatable :: temp_array
1540
1541 integer(i4) :: ii, idomain, iilc, s0, e0
1542
1543 type(grid), pointer :: level0_idomain => null()
1544
1545
1546 ! initialization
1547 do idomain = 1, domainmeta%nDomains
1548 allocate(temp_array(level11(idomain)%nCells, nlcoverscene))
1549 temp_array = nodata_dp
1550 if (do_init) then
1551 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
1552
1553 s0 = level0_idomain%iStart
1554 e0 = level0_idomain%iEnd
1555 nlinks = level11(idomain)%nCells + 1 - l11_noutlets(idomain)
1556 nlinkafloodplain => l11_afloodplain(level11(idomain)%iStart : level11(idomain)%iEnd)
1557
1558 do iilc = 1, nlcoverscene
1559 ! for a single node model run
1560 if(nlinks .GT. 0) then
1561 do ii = 1, nlinks
1562 temp_array(ii, iilc) = sum(level0_idomain%CellArea(:), &
1563 mask = (l0_floodplain(s0 : e0) == ii .and. l0_lcover(s0 : e0, iilc) == lcclassimp)) &
1564 / nlinkafloodplain(ii)
1565 end do
1566 end if
1567 end do
1568 end if
1569 call append(l11_nlinkfracfpimp, temp_array(:,:))
1570 deallocate(temp_array)
1571 end do
1572
1573 end subroutine l11_fraction_sealed_floodplain
1574
1575 ! ------------------------------------------------------------------
1576 ! MOVE UPSTREAM FROM-TO
1577 ! ------------------------------------------------------------------
1578 ! NAME
1579 ! moveUp
1580
1581 ! PURPOSE
1582 !> \brief TODO: add description
1583
1584 !> \details TODO: add description
1585
1586 ! INTENT(IN)
1587 !> \param[in] "real(dp), dimension(:, :) :: elev0"
1588 !> \param[in] "integer(i4), dimension(:, :) :: fDir0"
1589 !> \param[in] "integer(i4) :: fi, fj" co-ordinate of the stream bed
1590 !> \param[in] "integer(i4) :: fi, fj" co-ordinate of the stream bed
1591
1592 ! INTENT(INOUT)
1593 !> \param[inout] "integer(i4), dimension(:, :) :: ss"
1594 !> \param[inout] "integer(i4) :: nn"
1595
1596 ! HISTORY
1597 !> \authors Robert Schweppe
1598
1599 !> \date Jun 2018
1600
1601 ! Modifications:
1602 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1603
1604 subroutine moveup(elev0, fDir0, fi, fj, ss, nn)
1605
1606 use mo_mrm_constants, only : deltah
1607 use mo_utils, only : ge, le
1608
1609 implicit none
1610
1611 real(dp), dimension(:, :), allocatable, intent(IN) :: elev0
1612
1613 integer(i4), dimension(:, :), allocatable, intent(IN) :: fDir0
1614
1615 ! co-ordinate of the stream bed
1616 integer(i4), intent(IN) :: fi, fj
1617
1618 integer(i4), dimension(:, :), intent(INOUT) :: ss
1619
1620 integer(i4), intent(INOUT) :: nn
1621
1622 integer(i4) :: ii, jj, ip, im, jp, jm
1623
1624 integer(i4) :: nrows, ncols
1625
1626
1627 ii = ss(1, 1)
1628 jj = ss(1, 2)
1629 ip = ii + 1
1630 im = ii - 1
1631 jp = jj + 1
1632 jm = jj - 1
1633
1634 nrows = size(fdir0, 1)
1635 ncols = size(fdir0, 2)
1636
1637 !E
1638 if (jp <= ncols) then
1639 if ((fdir0(ii, jp) == 16) .and. &
1640 (le((elev0(ii, jp) - elev0(fi, fj)), deltah)) .and. &
1641 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1642 ) then
1643 nn = nn + 1
1644 ss(nn, 1) = ii
1645 ss(nn, 2) = jp
1646 !print *, i,jp
1647 end if
1648 end if
1649
1650 !SE
1651 if ((ip <= nrows) .and. &
1652 (jp <= ncols)) then
1653 if ((fdir0(ip, jp) == 32) .and. &
1654 (le((elev0(ip, jp) - elev0(fi, fj)), deltah)) .and. &
1655 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1656 ) then
1657 nn = nn + 1
1658 ss(nn, 1) = ip
1659 ss(nn, 2) = jp
1660 !print *, ip,jp
1661 end if
1662 end if
1663
1664 !S
1665 if ((ip <= nrows) .and. &
1666 (jp <= ncols)) then
1667 if ((fdir0(ip, jj) == 64) .and. &
1668 (le((elev0(ip, jj) - elev0(fi, fj)), deltah)) .and. &
1669 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1670 ) then
1671 nn = nn + 1
1672 ss(nn, 1) = ip
1673 ss(nn, 2) = jj
1674 !print *, ip,j
1675 end if
1676 end if
1677
1678 !SW
1679 if ((ip <= nrows) .and. &
1680 (jp <= ncols) .and. &
1681 (jm >= 1)) then
1682 if ((fdir0(ip, jm) == 128) .and. &
1683 (le((elev0(ip, jm) - elev0(fi, fj)), deltah)) .and. &
1684 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1685 ) then
1686 nn = nn + 1
1687 ss(nn, 1) = ip
1688 ss(nn, 2) = jm
1689 !print *, ip,jm
1690 end if
1691 end if
1692
1693 !W
1694 if ((jm >= 1) .and. &
1695 (jp <= ncols)) then
1696 if ((fdir0(ii, jm) == 1) .and. &
1697 (le((elev0(ii, jm) - elev0(fi, fj)), deltah)) .and. &
1698 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1699 ) then
1700 nn = nn + 1
1701 ss(nn, 1) = ii
1702 ss(nn, 2) = jm
1703 !print *, i,jm
1704 end if
1705 end if
1706
1707 !NW
1708 if ((im >= 1) .and. &
1709 (jp <= ncols) .and. &
1710 (jm >= 1)) then
1711 if ((fdir0(im, jm) == 2) .and. &
1712 (le((elev0(im, jm) - elev0(fi, fj)), deltah)) .and. &
1713 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1714 ) then
1715 nn = nn + 1
1716 ss(nn, 1) = im
1717 ss(nn, 2) = jm
1718 !print *, im,jm
1719 end if
1720 end if
1721
1722 !N
1723 if ((im >= 1) .and. &
1724 (jp <= ncols)) then
1725 if ((fdir0(im, jj) == 4) .and. &
1726 (le((elev0(im, jj) - elev0(fi, fj)), deltah)) .and. &
1727 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1728 ) then
1729 nn = nn + 1
1730 ss(nn, 1) = im
1731 ss(nn, 2) = jj
1732 !print *, im,j
1733 end if
1734 end if
1735
1736 !NE
1737 if ((im >= 1) .and. &
1738 (jp <= ncols)) then
1739 if ((fdir0(im, jp) == 8) .and. &
1740 (le((elev0(im, jp) - elev0(fi, fj)), deltah)) .and. &
1741 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1742 ) then
1743 nn = nn + 1
1744 ss(nn, 1) = im
1745 ss(nn, 2) = jp
1746 !print *, im,jp
1747 end if
1748 end if
1749
1750 end subroutine moveup
1751
1752 ! ------------------------------------------------------------------
1753 ! MOVE DOWNSTREAM
1754 ! ------------------------------------------------------------------
1755 ! NAME
1756 ! moveDownOneCell
1757
1758 ! PURPOSE
1759 !> \brief TODO: add description
1760
1761 !> \details TODO: add description
1762
1763 ! INTENT(IN)
1764 !> \param[in] "integer(i4) :: fDir"
1765
1766 ! INTENT(INOUT)
1767 !> \param[inout] "integer(i4) :: iRow, jCol"
1768 !> \param[inout] "integer(i4) :: iRow, jCol"
1769
1770 ! HISTORY
1771 !> \authors Robert Schweppe
1772
1773 !> \date Jun 2018
1774
1775 ! Modifications:
1776
1777 subroutine movedownonecell(fDir, iRow, jCol)
1778 implicit none
1779
1780 integer(i4), intent(IN) :: fDir
1781
1782 integer(i4), intent(INOUT) :: iRow, jCol
1783
1784
1785 select case (fdir)
1786 case(1) !E
1787 jcol = jcol + 1
1788 case(2) !SE
1789 irow = irow + 1
1790 jcol = jcol + 1
1791 case(4) !S
1792 irow = irow + 1
1793 case(8) !SW
1794 irow = irow + 1
1795 jcol = jcol - 1
1796 case(16) !W
1797 jcol = jcol - 1
1798 case(32) !NW
1799 irow = irow - 1
1800 jcol = jcol - 1
1801 case(64) !N
1802 irow = irow - 1
1803 case(128) !NE
1804 irow = irow - 1
1805 jcol = jcol + 1
1806 case default !sink
1807 ! do nothing
1808 end select
1809
1810 end subroutine movedownonecell
1811
1812 ! ------------------------------------------------------------------
1813 ! CELL LENGTH
1814 ! ------------------------------------------------------------------
1815 ! NAME
1816 ! cellLength
1817
1818 ! PURPOSE
1819 !> \brief TODO: add description
1820
1821 !> \details TODO: add description
1822
1823 ! INTENT(IN)
1824 !> \param[in] "integer(i4) :: iDomain"
1825 !> \param[in] "integer(i4) :: fDir"
1826 !> \param[in] "integer(i4) :: iRow"
1827 !> \param[in] "integer(i4) :: jCol"
1828 !> \param[in] "integer(i4) :: iCoorSystem"
1829
1830 ! INTENT(OUT)
1831 !> \param[out] "real(dp) :: length"
1832
1833 ! HISTORY
1834 !> \authors Robert Schweppe
1835
1836 !> \date Jun 2018
1837
1838 ! Modifications:
1839
1840 subroutine celllength(iDomain, fDir, iRow, jCol, iCoorSystem, length)
1841
1842 use mo_common_types, only: grid
1844 use mo_constants, only : sqrt2_dp
1845
1846 implicit none
1847
1848 integer(i4), intent(IN) :: iDomain
1849
1850 integer(i4), intent(IN) :: fDir
1851
1852 integer(i4), intent(IN) :: iRow
1853
1854 integer(i4), intent(IN) :: jCol
1855
1856 integer(i4), intent(IN) :: iCoorSystem
1857
1858 real(dp), intent(OUT) :: length
1859
1860 integer(i4) :: iRow_to, jCol_to
1861
1862 real(dp) :: lat_1, long_1, lat_2, long_2
1863
1864 type(grid), pointer :: level0_iDomain => null()
1865
1866
1867 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
1868
1869 ! regular X-Y cordinate system
1870 IF(icoorsystem .EQ. 0) THEN
1871
1872 select case (fdir)
1873 case(1, 4, 16, 64) ! E, S, W, N
1874 length = 1.0_dp
1875 case(2, 8, 32, 128) ! SE, SW, NW, NE
1876 length = sqrt2_dp
1877 end select
1878 length = length * level0_idomain%cellsize
1879
1880 ! regular lat-lon cordinate system
1881 ELSE IF(icoorsystem .EQ. 1) THEN
1882 irow_to = irow
1883 jcol_to = jcol
1884
1885 ! move in the direction of flow
1886 call movedownonecell(fdir, irow_to, jcol_to)
1887
1888 ! estimate lat-lon points
1889 lat_1 = level0_idomain%yllcorner + real((level0_idomain%ncols - jcol), dp) * level0_idomain%cellsize + &
1890 0.5_dp * level0_idomain%cellsize
1891 long_1 = level0_idomain%xllcorner + real((irow - 1), dp) * level0_idomain%cellsize + &
1892 0.5_dp * level0_idomain%cellsize
1893
1894 lat_2 = level0_idomain%yllcorner + real((level0_idomain%ncols - jcol_to), dp) * level0_idomain%cellsize + &
1895 0.5_dp * level0_idomain%cellsize
1896 long_2 = level0_idomain%xllcorner + real((irow_to - 1), dp) * level0_idomain%cellsize + &
1897 0.5_dp * level0_idomain%cellsize
1898 ! get distance between two points
1899 call get_distance_two_lat_lon_points(lat_1, long_1, lat_2, long_2, length)
1900
1901 END IF
1902 !
1903 end subroutine celllength
1904
1905
1906 ! --------------------------------------------------------------------------
1907
1908 ! NAME
1909 ! get_distance_two_lat_lon_points
1910
1911 ! PURPOSE
1912 !> \brief estimate distance in [m] between two points in a lat-lon
1913
1914 !> \details estimate distance in [m] between two points in a lat-lon
1915 !> Code is based on one that is implemented in the VIC-3L model
1916
1917 ! INTENT(IN)
1918 !> \param[in] "real(dp) :: lat1, long1, lat2, long2" latitude of point-1
1919 !> \param[in] "real(dp) :: lat1, long1, lat2, long2" longitude of point-1
1920 !> \param[in] "real(dp) :: lat1, long1, lat2, long2" latitude of point-2
1921 !> \param[in] "real(dp) :: lat1, long1, lat2, long2" longitude of point-2
1922
1923 ! INTENT(OUT)
1924 !> \param[out] "real(dp) :: distance_out" distance between two points [m]
1925
1926 ! HISTORY
1927 !> \authors Rohini Kumar
1928
1929 !> \date May 2014
1930
1931 ! Modifications:
1932 ! Stephan Thober Aug 2015 - ported to mRM
1933 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1934
1935 subroutine get_distance_two_lat_lon_points(lat1, long1, lat2, long2, distance_out)
1936
1937 use mo_constants, only : radiusearth_dp, twopi_dp
1938
1939 implicit none
1940
1941 ! longitude of point-2
1942 real(dp), intent(in) :: lat1, long1, lat2, long2
1943
1944 ! distance between two points [m]
1945 real(dp), intent(out) :: distance_out
1946
1947 real(dp) :: theta1
1948
1949 real(dp) :: phi1
1950
1951 real(dp) :: theta2
1952
1953 real(dp) :: phi2
1954
1955 real(dp) :: dtor
1956
1957 real(dp) :: term1
1958
1959 real(dp) :: term2
1960
1961 real(dp) :: term3
1962
1963 real(dp) :: temp
1964
1965
1966 dtor = twopi_dp / 360.0_dp
1967 theta1 = dtor * long1
1968 phi1 = dtor * lat1
1969 theta2 = dtor * long2
1970 phi2 = dtor * lat2
1971
1972 term1 = cos(phi1) * cos(theta1) * cos(phi2) * cos(theta2)
1973 term2 = cos(phi1) * sin(theta1) * cos(phi2) * sin(theta2)
1974 term3 = sin(phi1) * sin(phi2)
1975 temp = term1 + term2 + term3
1976 if(temp .GT. 1.0_dp) temp = 1.0_dp
1977
1978 distance_out = radiusearth_dp * acos(temp);
1979
1980 end subroutine get_distance_two_lat_lon_points
1981
1982 ! --------------------------------------------------------------------------
1983
1984 ! NAME
1985 ! L11_flow_accumulation
1986
1987 ! PURPOSE
1988
1989 !> \brief Calculates L11 flow accumulation per grid cell
1990 !> \details Calculates L11 flow accumulation per grid cell using L11_fDir
1991 !> and L11_cellarea. L11_flow_accumulation contains the recursiv subroutine
1992 !> calculate_L11_flow_accumulation
1993
1994 ! INTENT(IN)
1995 !> iDomain
1996
1997 ! INTENT(INOUT)
1998 ! None
1999
2000 ! INTENT(OUT)
2001 ! L11_fAcc, L11_LinkIn_fAcc
2002
2003 ! INTENT(IN), OPTIONAL
2004 ! None
2005
2006 ! INTENT(INOUT), OPTIONAL
2007 ! None
2008
2009 ! INTENT(OUT), OPTIONAL
2010 ! None
2011
2012 ! RETURN
2013 ! None
2014
2015 ! RESTRICTIONS
2016 ! None
2017
2018 ! EXAMPLE
2019 ! None
2020
2021 ! LITERATURE
2022 ! None
2023
2024 ! HISTORY
2025 !> \author Matthias Kelbling
2026 !> \date Aug 2017
2027 ! Modified
2028 ! Stephan Thober, Jun 2018 - refactored to fit MPR_extract
2029
2030 ! --------------------------------------------------------------------------
2031 subroutine l11_flow_accumulation(iDomain)
2032
2033 use mo_mrm_global_variables, only: &
2034 l11_fdir, & ! IN: flow direction at L11 (standard notation)
2035 l11_facc ! OUT: flow accumulation at L11 [km^2]
2038 use mo_append, only : append
2039
2040 implicit none
2041
2042 ! local variables
2043 integer(i4), intent(in) :: idomain
2044 real(dp), dimension(:,:), allocatable :: facc11 ! L11_fAcc array
2045 integer(i4) :: ii, jj ! row and col index
2046 integer(i4) :: s11, e11 ! Vec. position iDomain - fAcc
2047 integer(i4) :: nrows11, ncols11 ! array size Domain
2048 integer(i4), dimension(:,:), allocatable :: fdir11 ! L11_fDir array
2049 logical, dimension(:,:), allocatable :: mask11 ! Domain mask
2050
2051 ! initialize Domain info
2052 nrows11 = level11(idomain)%nrows
2053 ncols11 = level11(idomain)%ncols
2054 s11 = level11(idomain)%iStart
2055 e11 = level11(idomain)%iEnd
2056 mask11 = level11(idomain)%mask
2057
2058 ! allocate arrays
2059 allocate(fdir11(nrows11, ncols11))
2060 allocate(facc11(nrows11, ncols11))
2061
2062 ! initialize
2063 fdir11(:,:) = nodata_i4
2064
2065 ! get data
2066 fdir11(:,:) = unpack( l11_fdir(s11:e11), mask11, nodata_i4 )
2067
2068 ! initialize fAcc11 with cell area
2069 facc11 = unpack( level11(idomain)%cellarea * 1.e-6, mask11, nodata_dp )
2070
2071 ! For each sink call "calculate_L11_flow_accumulation"
2072 do jj=1, ncols11
2073 do ii=1, nrows11
2074 if (fdir11(ii,jj) .eq. 0) then
2075 call calculate_l11_flow_accumulation(fdir = fdir11, &
2076 facc = facc11, &
2077 ii = ii, &
2078 jj = jj, &
2079 nrow = nrows11, &
2080 ncol = ncols11)
2081 end if
2082 end do
2083 end do
2084
2085 ! Append
2086 call append( l11_facc, pack(facc11(:,:),mask11))
2087
2088 ! free space
2089 deallocate(fdir11, facc11, mask11)
2090
2091 contains
2092
2093 recursive subroutine calculate_l11_flow_accumulation(fDir, fAcc, ii, jj, nrow, ncol)
2094
2095 implicit none
2096
2097 integer(i4), intent(in) :: fdir(:,:) ! flow Direction
2098 real(dp), intent (inout) :: facc(:,:) ! flow accumulation
2099 integer(i4), intent(in) :: ii, jj ! row and col index
2100 integer(i4), intent(in) :: nrow, ncol ! number of rows,cols in array
2101
2102 ! Scan order:
2103 !
2104 ! 6 7 8
2105 ! 5 x 1
2106 ! 4 3 2
2107 !
2108 ! Each:
2109 ! 1. Check if cell is inflow cell to current grid
2110 ! 2. If yes: Call calculate_subroutine and add result
2111
2112 if (jj+1 .le. ncol) then
2113 if (fdir(ii,jj+1) .eq. 16_i4) then
2114 call calculate_l11_flow_accumulation(fdir, facc, ii, jj+1, nrow, ncol)
2115 facc(ii,jj) = facc(ii,jj) + facc(ii,jj+1)
2116 end if
2117 end if
2118
2119 if ((ii+1 .le. nrow) .and. (jj+1 .le. ncol)) then
2120 if (fdir(ii+1,jj+1) .eq. 32_i4) then
2121 call calculate_l11_flow_accumulation(fdir, facc, ii+1, jj+1, nrow, ncol)
2122 facc(ii,jj) = facc(ii,jj) + facc(ii+1,jj+1)
2123 end if
2124 end if
2125
2126 if (ii+1 .le. nrow) then
2127 if (fdir(ii+1,jj) .eq. 64_i4) then
2128 call calculate_l11_flow_accumulation(fdir, facc, ii+1, jj, nrow, ncol)
2129 facc(ii,jj) = facc(ii,jj) + facc(ii+1,jj)
2130 end if
2131 end if
2132
2133 if ((ii+1 .le. nrow) .and. (jj-1 .ge. 1)) then
2134 if (fdir(ii+1,jj-1) .eq. 128_i4) then
2135 call calculate_l11_flow_accumulation(fdir, facc, ii+1, jj-1, nrow, ncol)
2136 facc(ii,jj) = facc(ii,jj) + facc(ii+1,jj-1)
2137 end if
2138 end if
2139
2140 if (jj-1 .ge. 1) then
2141 if (fdir(ii,jj-1) .eq. 1_i4) then
2142 call calculate_l11_flow_accumulation(fdir, facc, ii, jj-1, nrow, ncol)
2143 facc(ii,jj) = facc(ii,jj) + facc(ii,jj-1)
2144 end if
2145 end if
2146
2147 if ((ii-1 .ge. 1) .and. (jj-1 .ge. 1)) then
2148 if (fdir(ii-1,jj-1) .eq. 2_i4) then
2149 call calculate_l11_flow_accumulation(fdir, facc, ii-1, jj-1, nrow, ncol)
2150 facc(ii,jj) = facc(ii,jj) + facc(ii-1,jj-1)
2151 end if
2152 end if
2153
2154 if (ii-1 .ge. 1) then
2155 if (fdir(ii-1,jj) .eq. 4_i4) then
2156 call calculate_l11_flow_accumulation(fdir, facc, ii-1, jj, nrow, ncol)
2157 facc(ii,jj) = facc(ii,jj) + facc(ii-1,jj)
2158 end if
2159 end if
2160
2161 if ((ii-1 .ge. 1) .and. (jj+1 .le. ncol)) then
2162 if (fdir11(ii-1,jj+1) .eq. 8_i4) then
2163 call calculate_l11_flow_accumulation(fdir, facc, ii-1, jj+1, nrow, ncol)
2164 facc(ii,jj) = facc(ii,jj) + facc(ii-1,jj+1)
2165 end if
2166 end if
2167
2168 ! print *, ii, jj, nrow, ncol, fAcc(ii, jj)
2169
2170 end subroutine calculate_l11_flow_accumulation
2171
2172 end subroutine l11_flow_accumulation
2173
2174
2175 ! ------------------------------------------------------------------
2176
2177 ! NAME
2178 ! L11_calc_celerity
2179
2180 ! PURPOSE
2181 !> \brief L11 celerity based on L0 elevation and L0 fAcc
2182
2183 !> \details L11 celerity based on L0 elevation and L0 fAcc
2184
2185 ! INTENT(IN)
2186 !> \param[in]
2187
2188 ! INTENT(INOUT)
2189 ! None
2190
2191 ! INTENT(OUT)
2192 ! None
2193
2194 ! INTENT(IN), OPTIONAL
2195 ! None
2196
2197 ! INTENT(INOUT), OPTIONAL
2198 ! None
2199
2200 ! INTENT(OUT), OPTIONAL
2201 ! None
2202
2203 ! RETURN
2204 ! None
2205
2206 ! RESTRICTIONS
2207 ! None
2208
2209 ! EXAMPLE
2210 ! None
2211
2212 ! LITERATURE
2213 ! None
2214
2215 ! HISTORY
2216 !> \author Matthias Kelbling
2217 !> \date Oct 2017
2218
2219 ! ------------------------------------------------------------------
2220
2221 subroutine l11_calc_celerity(iDomain, param)
2223 use mo_mad, only: mad
2224 use mo_append, only: append
2225 use mo_mpr_global_variables, only: &
2226 l0_slope ! IN: slope [%]
2227 use mo_common_types, only: grid
2228 use mo_common_variables, only: &
2229 domainmeta, & ! IN: for L0 Domain indexer
2230 level0 ! IN: level 0 grid
2231 use mo_mrm_global_variables, only: &
2232 l0_fdir, & ! IN: flow direction (standard notation) L0
2233 l0_facc, & ! IN: flow accumulation (number of cells)?
2234 l0_streamnet, & ! IN: stream Network at Level 0
2235 level11, & ! IN: level 11 grid
2236 l11_frow, & ! IN: from row in L0 grid
2237 l11_fcol, & ! IN: from col in L0 grid
2238 l11_trow, & ! IN: to row in L0 grid
2239 l11_tcol, & ! IN: to col in L0 grid
2240 l11_netperm, & ! IN: routing order (permutation)
2241 l11_noutlets, & ! IN: Number of Outlets/Sinks
2242 l11_celerity, & ! INOUT: averaged celerity
2243 l0_celerity ! INOUT
2244
2245 implicit none
2246
2247 integer(i4), intent(in) :: idomain ! Domain
2248 real(dp), dimension(:), intent(in) :: param
2249
2250 ! local
2251 integer(i4) :: ncells0
2252 integer(i4) :: nnodes
2253 integer(i4) :: nlinks
2254 integer(i4) :: nrows0, ncols0
2255 integer(i4) :: istart0, iend0
2256 integer(i4) :: nrows11, ncols11
2257 integer(i4) :: istart11, iend11
2258 logical, dimension(:,:), allocatable :: mask0
2259 integer(i4), dimension(:,:), allocatable :: id0
2260 integer(i4), dimension(:,:), allocatable :: fdir0
2261 integer(i4), dimension(:,:), allocatable :: facc0
2262 real(dp), dimension(:,:), allocatable :: slope0
2263 real(dp), dimension(:), allocatable :: slope_tmp
2264 real(dp), dimension(:,:), allocatable :: cellarea0
2265 integer(i4), dimension(:), allocatable :: netperm ! routing order (permutation)
2266 integer(i4), dimension(:), allocatable :: nlinkfromrow
2267 integer(i4), dimension(:), allocatable :: nlinkfromcol
2268 integer(i4), dimension(:), allocatable :: nlinktorow
2269 integer(i4), dimension(:), allocatable :: nlinktocol
2270 integer(i4) :: ii, rr, ns
2271 integer(i4) :: frow, fcol
2272 integer(i4) :: fid, tid
2273 real(dp), dimension(:), allocatable :: stack, append_chunk ! Stacks celerity along the L0 river-path
2274 integer(i4), dimension(:), allocatable :: dummy_1d
2275
2276 real(dp) :: l0_link_slope
2277 real(dp), dimension(:), allocatable :: celerity11
2278 real(dp), dimension(:,:), allocatable :: celerity0
2279
2280 integer(i4), dimension(:,:), allocatable :: nodata_i4_tmp
2281 real(dp), dimension(:,:), allocatable :: nodata_dp_tmp
2282 logical, dimension(:), allocatable :: slopemask0
2283
2284 type(grid), pointer :: level0_idomain
2285
2286 ! level-0 information
2287 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
2288 nrows0 = level0_idomain%nrows
2289 ncols0 = level0_idomain%ncols
2290 ncells0 = level0_idomain%ncells
2291 istart0 = level0_idomain%iStart
2292 iend0 = level0_idomain%iEnd
2293 mask0 = level0_idomain%mask
2294
2295 ! level-11 information
2296 istart11 = level11(idomain)%iStart
2297 iend11 = level11(idomain)%iEnd
2298 nrows11 = level11(idomain)%nrows
2299 ncols11 = level11(idomain)%ncols
2300 nnodes = level11(idomain)%ncells
2301
2302 nlinks = nnodes - l11_noutlets(idomain)
2303
2304 ! allocate
2305 allocate ( id0( nrows0, ncols0 ) )
2306 allocate ( slope0( nrows0, ncols0 ) )
2307 allocate ( fdir0( nrows0, ncols0 ) )
2308 allocate ( facc0( nrows0, ncols0 ) )
2309 allocate ( cellarea0( nrows0, ncols0 ) )
2310 allocate ( celerity0( nrows0, ncols0 ) )
2311 allocate ( slopemask0( ncells0 ) )
2312
2313 ! Routing network vectors have nNodes size instead of nLinks to
2314 ! avoid the need of having two extra indices to identify a Domain.
2315 allocate ( stack( 1 ) )
2316 allocate ( append_chunk( 1 ) )
2317 allocate ( dummy_1d( 2 ))
2318 allocate ( netperm( nnodes ) )
2319 allocate ( nlinkfromrow( nnodes ) )
2320 allocate ( nlinkfromcol( nnodes ) )
2321 allocate ( nlinktorow( nnodes ) )
2322 allocate ( nlinktocol( nnodes ) )
2323 allocate ( celerity11( nnodes ) )
2324 allocate ( slope_tmp( nnodes ) )
2325
2326 allocate (nodata_i4_tmp( nrows0, ncols0 ) )
2327 allocate (nodata_dp_tmp( nrows0, ncols0 ) )
2328
2329 ! initialize
2330 id0(:,:) = nodata_i4
2331 fdir0(:,:) = nodata_i4
2332 facc0(:,:) = nodata_i4
2333 cellarea0(:,:) = nodata_dp
2334 slope0(:,:) = nodata_dp
2335
2336 stack(:) = nodata_dp
2337 append_chunk(:) = nodata_dp
2338 netperm(:) = nodata_i4
2339 nlinkfromrow(:) = nodata_i4
2340 nlinkfromcol(:) = nodata_i4
2341 nlinktorow(:) = nodata_i4
2342 nlinktocol(:) = nodata_i4
2343 celerity11(:) = nodata_dp
2344 celerity0(:,:) = nodata_dp
2345 slopemask0(:) = .false.
2346
2347 nodata_i4_tmp(:,:) = nodata_i4
2348 nodata_dp_tmp(:,:) = nodata_dp
2349
2350 ! for a single node model run
2351
2352 if(nnodes .GT. 1) then
2353 ! get L0 fields
2354 id0(:,:) = unpack(level0_idomain%Id(1:ncells0), mask0, nodata_i4_tmp)
2355 fdir0(:,:) = unpack(l0_fdir(istart0:iend0), mask0, nodata_i4_tmp)
2356 facc0(:,:) = unpack(l0_facc(istart0:iend0), mask0, nodata_i4_tmp)
2357 cellarea0(:,:) = unpack(level0_idomain%cellarea(1:ncells0), mask0, nodata_dp_tmp)
2358
2359 ! smoothing river slope
2360 slope_tmp = l0_slope(istart0:iend0)
2361 where ( slope_tmp .lt. 0.1_dp ) slope_tmp = 0.1_dp
2362
2363 slopemask0(:) = (l0_streamnet(istart0:iend0) .ne. nodata_i4)
2364
2365 ! smooth river cells if there is more than one cell
2366 if( count(slopemask0) .GT. 1) then
2367 slope_tmp = mad(arr = slope_tmp, z = 2.25_dp, mask = slopemask0, tout='u', mval=0.1_dp)
2368 end if
2369 slope0(:,:) = unpack(slope_tmp, mask0, nodata_dp_tmp )
2370
2371 ! get network vectors of L11
2372 netperm(:) = l11_netperm( istart11 : iend11 )
2373 nlinkfromrow(:) = l11_frow( istart11 : iend11 )
2374 nlinkfromcol(:) = l11_fcol( istart11 : iend11 )
2375 nlinktorow(:) = l11_trow( istart11 : iend11 )
2376 nlinktocol(:) = l11_tcol( istart11 : iend11 )
2377
2378 do rr = 1, nlinks
2379
2380 ii = netperm(rr)
2381 frow = nlinkfromrow(ii)
2382 fcol = nlinkfromcol(ii)
2383
2384 ! Init
2385 stack(:) = 0_dp
2386 ns = 1
2387
2388 fid = id0( frow, fcol )
2389 tid = id0( nlinktorow(ii) , nlinktocol(ii) )
2390 do
2391 l0_link_slope = slope0(frow, fcol) / 100._dp
2392
2393 ! celerity parametrization
2394 stack(ns) = param(1) * sqrt(l0_link_slope)
2395
2396 celerity0(frow, fcol) = stack(ns)
2397 ns = ns + 1
2398 fid = id0(frow, fcol)
2399 if( .NOT. (fid == tid)) then
2400 call append(stack, append_chunk)
2401 else
2402 exit
2403 end if
2404 ! move downstream
2405 call movedownonecell( fdir0(frow,fcol), frow, fcol )
2406 end do
2407
2408 celerity11(ii) = size(stack) / sum(1/stack(:))
2409 deallocate(stack)
2410 allocate(stack(1))
2411
2412 end do
2413
2414 else
2415
2416 ! There is only one cell, so no routing is taking place
2417 ! set dummy value of 1 m / s
2418 celerity11(:) = 1._dp
2419
2420 end if
2421
2422 ! Write celerity
2423 l11_celerity(istart11:iend11) = celerity11(:)
2424 l0_celerity(istart0:iend0) = pack(celerity0(:,:), mask0)
2425
2426 ! free space
2427 deallocate (&
2428 mask0, id0, slope_tmp, slopemask0, &
2429 slope0, fdir0, cellarea0, &
2430 stack, netperm, nlinkfromrow, nlinkfromcol, nlinktorow, nlinktocol)
2431
2432 end subroutine l11_calc_celerity
2433
2434end module mo_mrm_net_startup
recursive subroutine calculate_l11_flow_accumulation(fdir, facc, ii, jj, nrow, ncol)
Provides constants commonly used by mHM, mRM and MPR.
real(dp), parameter, public nodata_dp
integer(i4), parameter, public nodata_i4
Provides common types needed by mHM, mRM and/or mpr.
Provides structures needed by mHM, mRM and/or mpr.
type(domain_meta), public domainmeta
integer(i4), public nlcoverscene
real(dp), dimension(:), allocatable, public l0_elev
integer(i4), public iflag_cordinate_sys
integer(i4), dimension(:, :), allocatable, public l0_lcover
integer(i4), dimension(nprocesses, 3), public processmatrix
type(grid), dimension(:), allocatable, target, public level1
type(grid), dimension(:), allocatable, target, public level0
Global variables for mpr only.
real(dp), dimension(:), allocatable, public l0_slope
Provides mRM specific constants.
real(dp), parameter, public deltah
Global variables for mRM only.
type(gridremapper), dimension(:), allocatable, public l0_l11_remap
integer(i4), dimension(:), allocatable, public l0_floodplain
real(dp), dimension(:), allocatable, public l11_facc
integer(i4), dimension(:), allocatable, public l11_netperm
real(dp), dimension(:, :), allocatable, public l11_nlinkfracfpimp
integer(i4), dimension(:), allocatable, public l11_l1_id
integer(i4), dimension(:), allocatable, public l1_l11_id
real(dp), dimension(:), allocatable, public l0_celerity
real(dp), dimension(:), allocatable, target, public l11_afloodplain
integer(i4), dimension(:), allocatable, public l0_fdir
integer(i4), dimension(:), allocatable, public l11_fcol
type(sink_cells_t), dimension(:), allocatable sink_cells
sink cell ids for each domain
integer(i4), dimension(:), allocatable, public l0_dracell
integer(i4), dimension(:), allocatable, public l11_label
integer(i4), dimension(:), allocatable, public l0_streamnet
integer(i4), dimension(:), allocatable, public l11_fromn
real(dp), dimension(:), allocatable, public l11_length
integer(i4), dimension(:), allocatable, public l11_rowout
integer(i4), dimension(:), allocatable, public l11_ton
type(domaininfo_mrm), dimension(:), allocatable, target, public domain_mrm
integer(i4), dimension(:), allocatable, public l11_trow
type(grid), dimension(:), allocatable, target, public level11
integer(i4), dimension(:), allocatable, public l11_fdir
integer(i4), dimension(:), allocatable, public l0_facc
integer(i4), dimension(:), allocatable, public l11_tcol
integer(i4), dimension(:), allocatable, public l11_frow
real(dp), dimension(:), allocatable, public l11_slope
logical, dimension(:), allocatable, public l11_sink
integer(i4), dimension(:), allocatable, public l0_inflowgaugeloc
real(dp), dimension(:), allocatable, public l11_celerity
integer(i4), dimension(:), allocatable, public l11_noutlets
integer(i4), dimension(:), allocatable, public l0_drasc
integer(i4), dimension(:), allocatable, public l11_rorder
integer(i4), dimension(:), allocatable, public l0_gaugeloc
integer(i4), dimension(:), allocatable, public l11_colout
Startup drainage network for mHM.
subroutine, public l11_stream_features(idomain)
Stream features (stream network and floodplain)
subroutine, public l11_flow_direction(idomain)
Determine the flow direction of the upscaled river network at level L11.
subroutine celllength(idomain, fdir, irow, jcol, icoorsystem, length)
TODO: add description.
subroutine moveup(elev0, fdir0, fi, fj, ss, nn)
TODO: add description.
subroutine, public get_distance_two_lat_lon_points(lat1, long1, lat2, long2, distance_out)
estimate distance in [m] between two points in a lat-lon
subroutine, public l11_fraction_sealed_floodplain(lcclassimp, do_init)
Fraction of the flood plain with impervious cover.
subroutine, public l11_set_drain_outlet_gauges(idomain)
Draining cell identification and Set gauging node.
subroutine, public l11_set_network_topology(idomain)
Set network topology.
subroutine, public l11_routing_order(idomain)
Find routing order, headwater cells and sink.
subroutine movedownonecell(fdir, irow, jcol)
TODO: add description.
subroutine, public l11_link_location(idomain)
Estimate the LO (row,col) location for each routing link at level L11.
subroutine, public l11_calc_celerity(idomain, param)
L11 celerity based on L0 elevation and L0 fAcc.
subroutine, public l11_l1_mapping(idomain)
TODO: add description.
subroutine, public l11_flow_accumulation(idomain)
Calculates L11 flow accumulation per grid cell.