5.13.2-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
733 level11
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) nlinksink(ii) = .true.
836 end do
837 where(nlinksink) nlinklabel = 2 ! 'Sink'
838
839 ! keep routing order
840 do ii = 1, nlinks
841 netperm(nlinkrorder(ii)) = ii
842 end do
843
844 ! end of multi-node network design loop
845 end if
846
847 !--------------------------------------------------------
848 ! Start padding up local variables to global variables
849 !--------------------------------------------------------
850 ! L11 network data sets
851 call append(l11_rorder, nlinkrorder(:))
852 call append(l11_label, nlinklabel(:))
853 call append(l11_sink, nlinksink(:))
854 call append(l11_netperm, netperm(:))
855
856 ! free space
857 deallocate (nlinkfromn, nlinkton, nlinkrorder, nlinklabel, nlinksink, netperm)
858
859 end subroutine l11_routing_order
860
861 ! ------------------------------------------------------------------
862
863 ! NAME
864 ! L11_link_location
865
866 ! PURPOSE
867 !> \brief Estimate the LO (row,col) location for each routing link at level L11
868
869 !> \details If a variable is added or removed here, then it also has to
870 !> be added or removed in the subroutine L11_config_set in
871 !> module mo_restart and in the subroutine set_L11_config in module
872 !> mo_set_netcdf_restart
873
874 ! INTENT(IN)
875 !> \param[in] "integer(i4) :: iDomain" Domain Id
876
877 ! HISTORY
878 !> \authors Luis Samaniego
879
880 !> \date Dec 2005
881
882 ! Modifications:
883 ! Luis Samaniego Jan 2013 - modular version
884 ! Stephan Thober Aug 2015 - ported to mRM
885 ! Robert Schweppe Jun 2018 - refactoring and reformatting
886
887 subroutine l11_link_location(iDomain)
888
889 use mo_append, only : append
891 use mo_common_types, only: grid
895 use mo_string_utils, only : num2str
896
897 implicit none
898
899 ! Domain Id
900 integer(i4), intent(in) :: idomain
901
902 integer(i4) :: nlinks
903
904 ! northing cell loc. of the Outlet
905 integer(i4), dimension(:), allocatable :: rowout
906
907 ! easting cell loc. of the Outlet
908 integer(i4), dimension(:), allocatable :: colout
909
910 integer(i4), dimension(:), allocatable :: nlinkfromn
911
912 integer(i4), dimension(:), allocatable :: netperm
913
914 integer(i4), dimension(:), allocatable :: nlinkfromrow
915
916 integer(i4), dimension(:), allocatable :: nlinkfromcol
917
918 integer(i4), dimension(:), allocatable :: nlinktorow
919
920 integer(i4), dimension(:), allocatable :: nlinktocol
921
922 integer(i4), dimension(:, :), allocatable :: fdir0
923
924 integer(i4), dimension(:, :), allocatable :: drasc0
925
926 integer(i4) :: ii, rr, kk, s0, e0
927
928 integer(i4) :: inode, irow, jcol, prevrow, prevcol
929
930 ! output location in L0
931 integer(i4), dimension(:, :), allocatable :: oloc
932
933 ! number of outlets in Domain
934 integer(i4) :: noutlets
935
936 ! flag for finding outlet
937 logical :: is_outlet
938
939 type(grid), pointer :: level0_idomain => null()
940
941
942 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
943 s0 = level0_idomain%iStart
944 e0 = level0_idomain%iEnd
945
946 noutlets = l11_noutlets(idomain)
947
948 nlinks = level11(idomain)%nCells - noutlets
949
950 ! Routing network vectors have level11(iDomain)%nCells size instead of nLinks to
951 ! avoid the need of having two extra indices to identify a Domain.
952 ! allocate
953 allocate (rowout(level11(idomain)%nCells))
954 allocate (colout(level11(idomain)%nCells))
955 allocate (nlinkfromn(level11(idomain)%nCells)) ! all network vectors valid from (1 : nLinks)
956 allocate (netperm(level11(idomain)%nCells))
957 allocate (nlinkfromrow(level11(idomain)%nCells))
958 allocate (nlinkfromcol(level11(idomain)%nCells))
959 allocate (nlinktorow(level11(idomain)%nCells))
960 allocate (nlinktocol(level11(idomain)%nCells))
961 allocate (fdir0(level0_idomain%nrows, level0_idomain%ncols))
962 allocate (drasc0(level0_idomain%nrows, level0_idomain%ncols))
963
964 ! initialize
965 rowout = nodata_i4
966 colout = nodata_i4
967 nlinkfromn = nodata_i4
968 netperm = nodata_i4
969 nlinkfromrow = nodata_i4
970 nlinkfromcol = nodata_i4
971 nlinktorow = nodata_i4
972 nlinktocol = nodata_i4
973 fdir0 = nodata_i4
974 drasc0 = nodata_i4
975
976 ! for a single node model run
977 if(level11(idomain)%nCells .GT. 1) then
978 ! get fDir at L0
979 fdir0(:, :) = unpack(l0_fdir(s0 : e0), level0_idomain%mask, nodata_i4)
980 drasc0(:, :) = unpack(l0_drasc(s0 : e0), level0_idomain%mask, nodata_i4)
981
982 ! get network vectors of L11
983 nlinkfromn(:) = l11_fromn(level11(idomain)%iStart : level11(idomain)%iEnd)
984 netperm(:) = l11_netperm(level11(idomain)%iStart : level11(idomain)%iEnd)
985 rowout(:) = l11_rowout(level11(idomain)%iStart : level11(idomain)%iEnd)
986 colout(:) = l11_colout(level11(idomain)%iStart : level11(idomain)%iEnd)
987
988 ! finding main outlet (row, col) in L0
989 allocate(oloc(noutlets, 2))
990 oloc(:, 1) = domain_mrm(idomain)%L0_rowOutlet(: noutlets)
991 oloc(:, 2) = domain_mrm(idomain)%L0_colOutlet(: noutlets)
992
993 ! Location of the stream-joint cells (row, col)
994 do rr = 1, nlinks
995
996 ii = netperm(rr)
997 inode = nlinkfromn(ii)
998 irow = rowout(inode)
999 jcol = colout(inode)
1000 call movedownonecell(fdir0(irow, jcol), irow, jcol)
1001 ! set "from" cell
1002 nlinkfromrow(ii) = irow
1003 nlinkfromcol(ii) = jcol
1004
1005 ! check whether this location is an outlet
1006 is_outlet = .false.
1007 do kk = 1, noutlets
1008 if (irow .eq. oloc(kk, 1) .and. jcol .eq. oloc(kk, 2)) is_outlet = .true.
1009 end do
1010
1011 if (is_outlet) then
1012
1013 nlinktorow(ii) = irow
1014 nlinktocol(ii) = jcol
1015
1016 else
1017
1018 do while (.not. (drasc0(irow, jcol) > 0))
1019 prevrow = irow
1020 prevcol = jcol
1021 call movedownonecell(fdir0(irow, jcol), irow, jcol)
1022 ! check whether this location is an outlet and exit
1023 do kk = 1, noutlets
1024 if (irow .eq. oloc(kk, 1) .and. jcol .eq. oloc(kk, 2)) exit
1025 end do
1026 if (prevrow .eq. irow .and. prevcol .eq. jcol) then
1027 call error_message('Something went wrong during L11_link_location, ', &
1028 'movedownonecell got stuck in infinite loop at cell (', num2str(irow), ' ', &
1029 num2str(jcol))
1030 end if
1031 end do
1032 ! set "to" cell (when an outlet is reached)
1033 nlinktorow(ii) = irow
1034 nlinktocol(ii) = jcol
1035
1036 end if
1037 end do
1038
1039 ! end of multi-node network design loop
1040 end if
1041
1042 !--------------------------------------------------------
1043 ! Start padding up local variables to global variables
1044 !--------------------------------------------------------
1045 ! L11 network data sets
1046 call append(l11_frow, nlinkfromrow(:))
1047 call append(l11_fcol, nlinkfromcol(:))
1048 call append(l11_trow, nlinktorow(:))
1049 call append(l11_tcol, nlinktocol(:))
1050
1051 ! free space
1052 deallocate (rowout, colout, nlinkfromn, netperm, nlinkfromrow, &
1053 nlinkfromcol, nlinktorow, nlinktocol, fdir0, drasc0)
1054
1055 end subroutine l11_link_location
1056
1057 ! ------------------------------------------------------------------
1058
1059 ! NAME
1060 ! L11_set_drain_outlet_gauges
1061
1062 ! PURPOSE
1063 !> \brief Draining cell identification and Set gauging node
1064
1065 !> \details Perform the following tasks:
1066 !> - Draining cell identification (cell at L0 to draining cell outlet at L11).
1067 !> - Set gauging nodes
1068 !> If a variable is added or removed here, then it also has to
1069 !> be added or removed in the subroutine L11_config_set in
1070 !> module mo_restart and in the subroutine set_L11_config in module
1071 !> mo_set_netcdf_restart
1072
1073 ! INTENT(IN)
1074 !> \param[in] "integer(i4) :: iDomain" Domain Id
1075
1076 ! HISTORY
1077 !> \authors Luis Samaniego
1078
1079 !> \date Dec 2005
1080
1081 ! Modifications:
1082 ! Luis Samaniego Jan 2013 - modular version
1083 ! Matthias Zink Mar 2014 - bugfix, added inflow gauge
1084 ! Rohini Kumar Apr 2014 - variable index is changed to index_gauge
1085 ! Stephan Thober Aug 2015 - ported to mRM
1086 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1087
1089
1090 use mo_append, only : append
1091 use mo_common_constants, only : nodata_i4
1092 use mo_common_types, only: grid
1096
1097 implicit none
1098
1099 ! Domain Id
1100 integer(i4), intent(in) :: idomain
1101
1102 integer(i4), dimension(:, :), allocatable :: drasc0
1103
1104 integer(i4), dimension(:, :), allocatable :: fdir0
1105
1106 integer(i4), dimension(:, :), allocatable :: gaugeloc0
1107
1108 integer(i4), dimension(:, :), allocatable :: inflowgaugeloc0
1109
1110 integer(i4), dimension(:, :), allocatable :: dracell0
1111
1112 integer(i4) :: ii, jj, kk, ll, s0, e0
1113
1114 integer(i4) :: isc
1115
1116 integer(i4) :: irow, jcol
1117
1118 type(grid), pointer :: level0_idomain => null()
1119
1120
1121 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
1122 s0 = level0_idomain%iStart
1123 e0 = level0_idomain%iEnd
1124
1125
1126 ! allocate
1127 allocate (drasc0(level0_idomain%nrows, level0_idomain%ncols))
1128 allocate (fdir0(level0_idomain%nrows, level0_idomain%ncols))
1129 allocate (gaugeloc0(level0_idomain%nrows, level0_idomain%ncols))
1130 allocate (inflowgaugeloc0(level0_idomain%nrows, level0_idomain%ncols))
1131 allocate (dracell0(level0_idomain%nrows, level0_idomain%ncols))
1132
1133 ! initialize
1134 drasc0(:, :) = nodata_i4
1135 fdir0(:, :) = nodata_i4
1136 gaugeloc0(:, :) = nodata_i4
1137 inflowgaugeloc0(:, :) = nodata_i4
1138 dracell0(:, :) = nodata_i4
1139
1140 drasc0(:, :) = unpack(l0_drasc(s0 : e0), &
1141 level0_idomain%mask, nodata_i4)
1142 fdir0(:, :) = unpack(l0_fdir(s0 : e0), &
1143 level0_idomain%mask, nodata_i4)
1144 gaugeloc0(:, :) = unpack(l0_gaugeloc(s0 : e0), &
1145 level0_idomain%mask, nodata_i4)
1146 inflowgaugeloc0(:, :) = unpack(l0_inflowgaugeloc(s0 : e0), &
1147 level0_idomain%mask, nodata_i4)
1148
1149 do kk = 1, level0_idomain%nCells
1150 ii = level0_idomain%CellCoor(kk, 1)
1151 jj = level0_idomain%CellCoor(kk, 2)
1152 isc = drasc0(ii, jj)
1153 ! find drainage path
1154 irow = ii
1155 jcol = jj
1156 do while (.NOT. isc > 0)
1157 ! move downstream
1158 call movedownonecell(fdir0(irow, jcol), irow, jcol)
1159 isc = drasc0(irow, jcol)
1160 end do
1161 dracell0(ii, jj) = isc
1162
1163 ! find cell at L11 corresponding to gauges in Domain at L0 ! >> L11_on_L0 is Id of
1164 ! the routing cell at level-11
1165 if (gaugeloc0(ii, jj) .NE. nodata_i4) then
1166 ! evaluation gauges
1167 do ll = 1, domain_mrm(idomain)%nGauges
1168 ! search for gaugeID in L0 grid and save ID on L11
1169 if (domain_mrm(idomain)%gaugeIdList(ll) .EQ. gaugeloc0(ii, jj)) then
1170 domain_mrm(idomain)%gaugeNodeList(ll) = l0_l11_remap(idomain)%lowres_id_on_highres(ii, jj)
1171 end if
1172 end do
1173 end if
1174
1175 if (inflowgaugeloc0(ii, jj) .NE. nodata_i4) then
1176 ! inflow gauges
1177 do ll = 1, domain_mrm(idomain)%nInflowGauges
1178 ! search for gaugeID in L0 grid and save ID on L11
1179 if (domain_mrm(idomain)%InflowGaugeIdList(ll) .EQ. inflowgaugeloc0(ii, jj)) &
1180 domain_mrm(idomain)%InflowGaugeNodeList(ll) = l0_l11_remap(idomain)%lowres_id_on_highres(ii, jj)
1181 end do
1182 end if
1183 end do
1184
1185 !--------------------------------------------------------
1186 ! Start padding up local variables to global variables
1187 !--------------------------------------------------------
1188 ! L0 data sets
1189 ! check whether L0 data is shared
1190 ! ToDo: check if change is correct
1191 if (idomain .eq. 1) then
1192 call append(l0_dracell, pack(dracell0(:, :), level0_idomain%mask))
1193 else if (domainmeta%L0DataFrom(idomain) == idomain) then
1194 call append(l0_dracell, pack(dracell0(:, :), level0_idomain%mask))
1195 end if
1196
1197 ! free space
1198 deallocate (drasc0, fdir0, gaugeloc0, dracell0)
1199
1200 end subroutine l11_set_drain_outlet_gauges
1201
1202 ! ------------------------------------------------------------------
1203
1204 ! NAME
1205 ! L11_stream_features
1206
1207 ! PURPOSE
1208 !> \brief Stream features (stream network and floodplain)
1209
1210 !> \details Stream features (stream network and floodplain)
1211 !> If a variable is added or removed here, then it also has to
1212 !> be added or removed in the subroutine L11_config_set in
1213 !> module mo_restart and in the subroutine set_L11_config in module
1214 !> mo_set_netcdf_restart
1215
1216 ! INTENT(IN)
1217 !> \param[in] "integer(i4) :: iDomain" Domain Id
1218
1219 ! HISTORY
1220 !> \authors Luis Samaniego
1221
1222 !> \date Dec 2005
1223
1224 ! Modifications:
1225 ! Luis Samaniego Jan 2013 - modular version
1226 ! R. Kumar Oct 2013 - stack size increased from nNodes to 100
1227 ! Stephan Thober Aug 2015 - ported to mRM
1228 ! Stephan Thober Nov 2016 - only read flood plain area if processMatrix for routing equals 1
1229 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1230 ! Stephan Thober Jul 2018 - introduced cut off Length at 40 percentile to neglect short paths in headwaters for adaptive timesteps
1231 ! 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
1232
1233 subroutine l11_stream_features(iDomain)
1234
1235 use mo_append, only : append
1237 use mo_common_types, only: grid
1239 use mo_mrm_global_variables, only : l0_fdir, &
1242 use mo_percentile, only: percentile
1243
1244 implicit none
1245
1246 ! Domain Id
1247 integer(i4), intent(in) :: idomain
1248
1249 integer(i4) :: nlinks
1250
1251 integer(i4), dimension(:, :), allocatable :: id0
1252
1253 integer(i4), dimension(:, :), allocatable :: fdir0
1254
1255 real(dp), dimension(:, :), allocatable :: elev0
1256
1257 real(dp), dimension(:, :), allocatable :: cellarea0
1258
1259 integer(i4), dimension(:, :), allocatable :: streamnet0
1260
1261 integer(i4), dimension(:, :), allocatable :: floodplain0
1262
1263 ! routing order (permutation)
1264 integer(i4), dimension(:), allocatable :: netperm
1265
1266 integer(i4), dimension(:), allocatable :: nlinkfromrow
1267
1268 integer(i4), dimension(:), allocatable :: nlinkfromcol
1269
1270 integer(i4), dimension(:), allocatable :: nlinktorow
1271
1272 integer(i4), dimension(:), allocatable :: nlinktocol
1273
1274 real(dp), dimension(:), allocatable :: nlinklength
1275
1276 real(dp), dimension(:), allocatable :: nlinkafloodplain
1277
1278 real(dp), dimension(:), allocatable :: nlinkslope
1279
1280 integer(i4) :: ii, rr, ns, s0, e0
1281
1282 integer(i4) :: frow, fcol
1283
1284 integer(i4) :: fid, tid
1285
1286 integer(i4), dimension(:, :), allocatable :: stack, append_chunk
1287
1288 integer(i4), dimension(:), allocatable :: dummy_1d
1289
1290 real(dp) :: length
1291
1292 integer(i4), dimension(:, :), allocatable :: nodata_i4_tmp
1293
1294 real(dp), dimension(:, :), allocatable :: nodata_dp_tmp
1295
1296 type(grid), pointer :: level0_idomain => null()
1297
1298
1299 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
1300 s0 = level0_idomain%iStart
1301 e0 = level0_idomain%iEnd
1302 nlinks = level11(idomain)%nCells - l11_noutlets(idomain)
1303
1304
1305 ! allocate
1306 allocate (id0(level0_idomain%nrows, level0_idomain%ncols))
1307 allocate (elev0(level0_idomain%nrows, level0_idomain%ncols))
1308 allocate (fdir0(level0_idomain%nrows, level0_idomain%ncols))
1309 allocate (cellarea0(level0_idomain%nrows, level0_idomain%ncols))
1310 allocate (streamnet0(level0_idomain%nrows, level0_idomain%ncols))
1311 allocate (floodplain0(level0_idomain%nrows, level0_idomain%ncols))
1312
1313 ! Routing network vectors have nNodes size instead of nLinks to
1314 ! avoid the need of having two extra indices to identify a Domain.
1315 allocate (stack(level11(idomain)%nCells, 2)) ! >> stack(nNodes, 2)
1316 allocate (dummy_1d(2))
1317 allocate (append_chunk(8, 2))
1318 allocate (netperm(level11(idomain)%nCells))
1319 allocate (nlinkfromrow(level11(idomain)%nCells))
1320 allocate (nlinkfromcol(level11(idomain)%nCells))
1321 allocate (nlinktorow(level11(idomain)%nCells))
1322 allocate (nlinktocol(level11(idomain)%nCells))
1323 allocate (nlinklength(level11(idomain)%nCells))
1324 allocate (nlinkafloodplain(level11(idomain)%nCells))
1325 allocate (nlinkslope(level11(idomain)%nCells))
1326
1327 allocate (nodata_i4_tmp(level0_idomain%nrows, level0_idomain%ncols))
1328 allocate (nodata_dp_tmp(level0_idomain%nrows, level0_idomain%ncols))
1329
1330 ! initialize
1331 id0(:, :) = nodata_i4
1332 elev0(:, :) = nodata_dp
1333 fdir0(:, :) = nodata_i4
1334 cellarea0(:, :) = nodata_dp
1335 streamnet0(:, :) = nodata_i4
1336 floodplain0(:, :) = nodata_i4
1337
1338 stack(:, :) = nodata_i4
1339 append_chunk(:, :) = nodata_i4
1340 netperm(:) = nodata_i4
1341 nlinkfromrow(:) = nodata_i4
1342 nlinkfromcol(:) = nodata_i4
1343 nlinktorow(:) = nodata_i4
1344 nlinktocol(:) = nodata_i4
1345 nlinklength(:) = nodata_dp
1346 nlinkafloodplain(:) = nodata_dp
1347 nlinkslope(:) = nodata_dp
1348
1349 nodata_i4_tmp(:, :) = nodata_i4
1350 nodata_dp_tmp(:, :) = nodata_dp
1351
1352 ! for a single node model run
1353 if(level11(idomain)%nCells .GT. 1) then
1354 ! get L0 fields
1355 id0(:, :) = unpack(level0_idomain%Id, level0_idomain%mask, nodata_i4_tmp)
1356 elev0(:, :) = unpack(l0_elev(s0 : e0), &
1357 level0_idomain%mask, nodata_dp_tmp)
1358 fdir0(:, :) = unpack(l0_fdir(s0 : e0), &
1359 level0_idomain%mask, nodata_i4_tmp)
1360 cellarea0(:, :) = unpack(level0_idomain%CellArea, level0_idomain%mask, nodata_dp_tmp)
1361
1362 ! get network vectors of L11
1363 netperm(:) = l11_netperm(level11(idomain)%iStart : level11(idomain)%iEnd)
1364 nlinkfromrow(:) = l11_frow(level11(idomain)%iStart : level11(idomain)%iEnd)
1365 nlinkfromcol(:) = l11_fcol(level11(idomain)%iStart : level11(idomain)%iEnd)
1366 nlinktorow(:) = l11_trow(level11(idomain)%iStart : level11(idomain)%iEnd)
1367 nlinktocol(:) = l11_tcol(level11(idomain)%iStart : level11(idomain)%iEnd)
1368
1369 ! Flood plains: stream network delineation
1370 streamnet0(:, :) = nodata_i4
1371 floodplain0(:, :) = nodata_i4
1372
1373 do rr = 1, nlinks
1374
1375 ii = netperm(rr)
1376 frow = nlinkfromrow(ii)
1377 fcol = nlinkfromcol(ii)
1378
1379 ! Init
1380 streamnet0(frow, fcol) = ii
1381 floodplain0(frow, fcol) = ii
1382 stack = 0
1383 append_chunk = 0
1384 ns = 1
1385 stack(ns, 1) = frow
1386 stack(ns, 2) = fcol
1387
1388 call celllength(idomain, fdir0(frow, fcol), frow, fcol, iflag_cordinate_sys, nlinklength(ii))
1389 nlinkslope(ii) = elev0(frow, fcol)
1390
1391 fid = id0(frow, fcol)
1392 tid = id0(nlinktorow(ii), nlinktocol(ii))
1393
1394 do while (.NOT. (fid == tid))
1395 ! Search flood plain from point(frow,fcol) upwards, keep co-ordinates in STACK
1396 do while (ns > 0)
1397 if (ns + 8 .gt. size(stack, 1)) then
1398 call append(stack, append_chunk)
1399 end if
1400 call moveup(elev0, fdir0, frow, fcol, stack, ns)
1401 stack(1, 1) = 0
1402 stack(1, 2) = 0
1403 ! stack = cshift(stack, SHIFT = 1, DIM = 1)
1404 ! substitute cshift <<<
1405 dummy_1d = stack(1, :)
1406 stack(: size(stack, dim = 1) - 1, :) = stack(2 :, :)
1407 stack(size(stack, dim = 1), :) = dummy_1d
1408 ! substitute cshift >>>
1409 if (stack(1, 1) > 0 .and. stack(1, 2) > 0) floodplain0(stack(1, 1), stack(1, 2)) = ii
1410 ns = count(stack > 0) / 2
1411 end do
1412
1413 ! move downstream
1414 call movedownonecell(fdir0(frow, fcol), frow, fcol)
1415 streamnet0(frow, fcol) = ii
1416 floodplain0(frow, fcol) = ii
1417 fid = id0(frow, fcol)
1418 stack = 0
1419 ns = 1
1420 stack(ns, 1) = frow
1421 stack(ns, 2) = fcol
1422 call celllength(idomain, fdir0(frow, fcol), frow, fcol, iflag_cordinate_sys, length)
1423 nlinklength(ii) = nlinklength(ii) + length
1424
1425 end do
1426
1427 ! stream bed slope
1428 nlinkslope(ii) = (nlinkslope(ii) - elev0(frow, fcol)) / nlinklength(ii)
1429
1430 if (nlinkslope(ii) < 0.0001_dp) nlinkslope(ii) = 0.0001_dp
1431
1432 ! calculate area of floodplains (avoid overwriting)
1433 nlinkafloodplain(ii) = sum(cellarea0(:, :), mask = (floodplain0(:, :) == ii))
1434 ! old > real( count( floodPlain0(:,:,) == i), dp ) * cellarea0
1435
1436 end do
1437
1438 ! end of multi-node network design loop
1439 end if
1440
1441 ! cut off Length at 40 percentile to neglect short paths in headwaters
1442 if ((processmatrix(8, 1) .eq. 2) .or. (processmatrix(8, 1) .eq. 3)) then
1443 if (count(nlinklength(:) .ge. 0._dp) .gt. 2) then
1444 length = percentile(pack(nlinklength(:), nlinklength(:) .ge. 0._dp), 40._dp)
1445 nlinklength(:) = merge(nlinklength(:), length, (nlinklength(:) .gt. length))
1446 end if
1447 end if
1448
1449 !--------------------------------------------------------
1450 ! Start padding up local variables to global variables
1451 !--------------------------------------------------------
1452
1453 ! L0 data sets
1454 ! check whether L0 data is shared
1455 ! ToDo: check if change is correct
1456 if (idomain .eq. 1) then
1457 call append(l0_streamnet, pack(streamnet0(:, :), level0_idomain%mask))
1458 call append(l0_floodplain, pack(floodplain0(:, :), level0_idomain%mask))
1459 else if (domainmeta%L0DataFrom(idomain) == idomain) then
1460 call append(l0_streamnet, pack(streamnet0(:, :), level0_idomain%mask))
1461 call append(l0_floodplain, pack(floodplain0(:, :), level0_idomain%mask))
1462 end if
1463
1464
1465 ! L11 network data sets
1466 call append(l11_length, nlinklength(:))
1467 call append(l11_afloodplain, nlinkafloodplain(:))
1468 call append(l11_slope, nlinkslope(:))
1469
1470 ! free space
1471 deallocate (&
1472 id0, elev0, fdir0, streamnet0, floodplain0, &
1473 cellarea0, stack, netperm, nlinkfromrow, nlinkfromcol, nlinktorow, nlinktocol, &
1474 nlinklength, nlinkafloodplain, nlinkslope, dummy_1d)
1475 deallocate(nodata_i4_tmp, nodata_dp_tmp)
1476
1477 end subroutine l11_stream_features
1478
1479 ! ------------------------------------------------------------------
1480
1481 ! NAME
1482 ! L11_fraction_sealed_floodplain
1483
1484 ! PURPOSE
1485 !> \brief Fraction of the flood plain with impervious cover
1486
1487 !> \details Fraction of the flood plain with impervious cover (\ref fig_routing "Routing
1488 !> Network"). This proportion is used to regionalize the Muskingum parameters.
1489 !> Samaniego et al. \cite SB05 found out that this fraction is one of the statistically
1490 !> significant predictor variables of peak discharge in mesoscale Domains.
1491 !> If a variable is added or removed here, then it also has to
1492 !> be added or removed in the subroutine L11_config_set in
1493 !> module mo_restart and in the subroutine set_L11_config in module
1494 !> mo_set_netcdf_restart
1495
1496 ! INTENT(IN)
1497 !> \param[in] "integer(i4) :: LCClassImp" Impervious land cover class Id, e.g. = 2 (old code)
1498 !> \param[in] "logical :: do_init"
1499
1500 ! HISTORY
1501 !> \authors Luis Samaniego
1502
1503 !> \date Dec 2005
1504
1505 ! Modifications:
1506 ! Luis Samaniego Jan 2013 - modular version
1507 ! Stephan Thober Aug 2015 - ported to mRM
1508 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1509
1510 subroutine l11_fraction_sealed_floodplain(LCClassImp, do_init)
1511
1512 use mo_append, only : append
1513 use mo_common_constants, only : nodata_dp
1514 use mo_common_types, only: grid
1518
1519 implicit none
1520
1521 ! Impervious land cover class Id, e.g. = 2 (old code)
1522 integer(i4), intent(in) :: lcclassimp
1523
1524 logical, intent(in) :: do_init
1525
1526 integer(i4) :: nlinks
1527
1528 real(dp), dimension(:), pointer :: nlinkafloodplain => null()
1529
1530 real(dp), dimension(:,:), allocatable :: temp_array
1531
1532 integer(i4) :: ii, idomain, iilc, s0, e0
1533
1534 type(grid), pointer :: level0_idomain => null()
1535
1536
1537 ! initialization
1538 do idomain = 1, domainmeta%nDomains
1539 allocate(temp_array(level11(idomain)%nCells, nlcoverscene))
1540 temp_array = nodata_dp
1541 if (do_init) then
1542 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
1543
1544 s0 = level0_idomain%iStart
1545 e0 = level0_idomain%iEnd
1546 nlinks = level11(idomain)%nCells + 1 - l11_noutlets(idomain)
1547 nlinkafloodplain => l11_afloodplain(level11(idomain)%iStart : level11(idomain)%iEnd)
1548
1549 do iilc = 1, nlcoverscene
1550 ! for a single node model run
1551 if(nlinks .GT. 0) then
1552 do ii = 1, nlinks
1553 temp_array(ii, iilc) = sum(level0_idomain%CellArea(:), &
1554 mask = (l0_floodplain(s0 : e0) == ii .and. l0_lcover(s0 : e0, iilc) == lcclassimp)) &
1555 / nlinkafloodplain(ii)
1556 end do
1557 end if
1558 end do
1559 end if
1560 call append(l11_nlinkfracfpimp, temp_array(:,:))
1561 deallocate(temp_array)
1562 end do
1563
1564 end subroutine l11_fraction_sealed_floodplain
1565
1566 ! ------------------------------------------------------------------
1567 ! MOVE UPSTREAM FROM-TO
1568 ! ------------------------------------------------------------------
1569 ! NAME
1570 ! moveUp
1571
1572 ! PURPOSE
1573 !> \brief TODO: add description
1574
1575 !> \details TODO: add description
1576
1577 ! INTENT(IN)
1578 !> \param[in] "real(dp), dimension(:, :) :: elev0"
1579 !> \param[in] "integer(i4), dimension(:, :) :: fDir0"
1580 !> \param[in] "integer(i4) :: fi, fj" co-ordinate of the stream bed
1581 !> \param[in] "integer(i4) :: fi, fj" co-ordinate of the stream bed
1582
1583 ! INTENT(INOUT)
1584 !> \param[inout] "integer(i4), dimension(:, :) :: ss"
1585 !> \param[inout] "integer(i4) :: nn"
1586
1587 ! HISTORY
1588 !> \authors Robert Schweppe
1589
1590 !> \date Jun 2018
1591
1592 ! Modifications:
1593 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1594
1595 subroutine moveup(elev0, fDir0, fi, fj, ss, nn)
1596
1597 use mo_mrm_constants, only : deltah
1598 use mo_utils, only : ge, le
1599
1600 implicit none
1601
1602 real(dp), dimension(:, :), allocatable, intent(IN) :: elev0
1603
1604 integer(i4), dimension(:, :), allocatable, intent(IN) :: fDir0
1605
1606 ! co-ordinate of the stream bed
1607 integer(i4), intent(IN) :: fi, fj
1608
1609 integer(i4), dimension(:, :), intent(INOUT) :: ss
1610
1611 integer(i4), intent(INOUT) :: nn
1612
1613 integer(i4) :: ii, jj, ip, im, jp, jm
1614
1615 integer(i4) :: nrows, ncols
1616
1617
1618 ii = ss(1, 1)
1619 jj = ss(1, 2)
1620 ip = ii + 1
1621 im = ii - 1
1622 jp = jj + 1
1623 jm = jj - 1
1624
1625 nrows = size(fdir0, 1)
1626 ncols = size(fdir0, 2)
1627
1628 !E
1629 if (jp <= ncols) then
1630 if ((fdir0(ii, jp) == 16) .and. &
1631 (le((elev0(ii, jp) - elev0(fi, fj)), deltah)) .and. &
1632 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1633 ) then
1634 nn = nn + 1
1635 ss(nn, 1) = ii
1636 ss(nn, 2) = jp
1637 !print *, i,jp
1638 end if
1639 end if
1640
1641 !SE
1642 if ((ip <= nrows) .and. &
1643 (jp <= ncols)) then
1644 if ((fdir0(ip, jp) == 32) .and. &
1645 (le((elev0(ip, jp) - elev0(fi, fj)), deltah)) .and. &
1646 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1647 ) then
1648 nn = nn + 1
1649 ss(nn, 1) = ip
1650 ss(nn, 2) = jp
1651 !print *, ip,jp
1652 end if
1653 end if
1654
1655 !S
1656 if ((ip <= nrows) .and. &
1657 (jp <= ncols)) then
1658 if ((fdir0(ip, jj) == 64) .and. &
1659 (le((elev0(ip, jj) - elev0(fi, fj)), deltah)) .and. &
1660 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1661 ) then
1662 nn = nn + 1
1663 ss(nn, 1) = ip
1664 ss(nn, 2) = jj
1665 !print *, ip,j
1666 end if
1667 end if
1668
1669 !SW
1670 if ((ip <= nrows) .and. &
1671 (jp <= ncols) .and. &
1672 (jm >= 1)) then
1673 if ((fdir0(ip, jm) == 128) .and. &
1674 (le((elev0(ip, jm) - elev0(fi, fj)), deltah)) .and. &
1675 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1676 ) then
1677 nn = nn + 1
1678 ss(nn, 1) = ip
1679 ss(nn, 2) = jm
1680 !print *, ip,jm
1681 end if
1682 end if
1683
1684 !W
1685 if ((jm >= 1) .and. &
1686 (jp <= ncols)) then
1687 if ((fdir0(ii, jm) == 1) .and. &
1688 (le((elev0(ii, jm) - elev0(fi, fj)), deltah)) .and. &
1689 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1690 ) then
1691 nn = nn + 1
1692 ss(nn, 1) = ii
1693 ss(nn, 2) = jm
1694 !print *, i,jm
1695 end if
1696 end if
1697
1698 !NW
1699 if ((im >= 1) .and. &
1700 (jp <= ncols) .and. &
1701 (jm >= 1)) then
1702 if ((fdir0(im, jm) == 2) .and. &
1703 (le((elev0(im, jm) - elev0(fi, fj)), deltah)) .and. &
1704 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1705 ) then
1706 nn = nn + 1
1707 ss(nn, 1) = im
1708 ss(nn, 2) = jm
1709 !print *, im,jm
1710 end if
1711 end if
1712
1713 !N
1714 if ((im >= 1) .and. &
1715 (jp <= ncols)) then
1716 if ((fdir0(im, jj) == 4) .and. &
1717 (le((elev0(im, jj) - elev0(fi, fj)), deltah)) .and. &
1718 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1719 ) then
1720 nn = nn + 1
1721 ss(nn, 1) = im
1722 ss(nn, 2) = jj
1723 !print *, im,j
1724 end if
1725 end if
1726
1727 !NE
1728 if ((im >= 1) .and. &
1729 (jp <= ncols)) then
1730 if ((fdir0(im, jp) == 8) .and. &
1731 (le((elev0(im, jp) - elev0(fi, fj)), deltah)) .and. &
1732 (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1733 ) then
1734 nn = nn + 1
1735 ss(nn, 1) = im
1736 ss(nn, 2) = jp
1737 !print *, im,jp
1738 end if
1739 end if
1740
1741 end subroutine moveup
1742
1743 ! ------------------------------------------------------------------
1744 ! MOVE DOWNSTREAM
1745 ! ------------------------------------------------------------------
1746 ! NAME
1747 ! moveDownOneCell
1748
1749 ! PURPOSE
1750 !> \brief TODO: add description
1751
1752 !> \details TODO: add description
1753
1754 ! INTENT(IN)
1755 !> \param[in] "integer(i4) :: fDir"
1756
1757 ! INTENT(INOUT)
1758 !> \param[inout] "integer(i4) :: iRow, jCol"
1759 !> \param[inout] "integer(i4) :: iRow, jCol"
1760
1761 ! HISTORY
1762 !> \authors Robert Schweppe
1763
1764 !> \date Jun 2018
1765
1766 ! Modifications:
1767
1768 subroutine movedownonecell(fDir, iRow, jCol)
1769 implicit none
1770
1771 integer(i4), intent(IN) :: fDir
1772
1773 integer(i4), intent(INOUT) :: iRow, jCol
1774
1775
1776 select case (fdir)
1777 case(1) !E
1778 jcol = jcol + 1
1779 case(2) !SE
1780 irow = irow + 1
1781 jcol = jcol + 1
1782 case(4) !S
1783 irow = irow + 1
1784 case(8) !SW
1785 irow = irow + 1
1786 jcol = jcol - 1
1787 case(16) !W
1788 jcol = jcol - 1
1789 case(32) !NW
1790 irow = irow - 1
1791 jcol = jcol - 1
1792 case(64) !N
1793 irow = irow - 1
1794 case(128) !NE
1795 irow = irow - 1
1796 jcol = jcol + 1
1797 case default !sink
1798 ! do nothing
1799 end select
1800
1801 end subroutine movedownonecell
1802
1803 ! ------------------------------------------------------------------
1804 ! CELL LENGTH
1805 ! ------------------------------------------------------------------
1806 ! NAME
1807 ! cellLength
1808
1809 ! PURPOSE
1810 !> \brief TODO: add description
1811
1812 !> \details TODO: add description
1813
1814 ! INTENT(IN)
1815 !> \param[in] "integer(i4) :: iDomain"
1816 !> \param[in] "integer(i4) :: fDir"
1817 !> \param[in] "integer(i4) :: iRow"
1818 !> \param[in] "integer(i4) :: jCol"
1819 !> \param[in] "integer(i4) :: iCoorSystem"
1820
1821 ! INTENT(OUT)
1822 !> \param[out] "real(dp) :: length"
1823
1824 ! HISTORY
1825 !> \authors Robert Schweppe
1826
1827 !> \date Jun 2018
1828
1829 ! Modifications:
1830
1831 subroutine celllength(iDomain, fDir, iRow, jCol, iCoorSystem, length)
1832
1833 use mo_common_types, only: grid
1835 use mo_constants, only : sqrt2_dp
1836
1837 implicit none
1838
1839 integer(i4), intent(IN) :: iDomain
1840
1841 integer(i4), intent(IN) :: fDir
1842
1843 integer(i4), intent(IN) :: iRow
1844
1845 integer(i4), intent(IN) :: jCol
1846
1847 integer(i4), intent(IN) :: iCoorSystem
1848
1849 real(dp), intent(OUT) :: length
1850
1851 integer(i4) :: iRow_to, jCol_to
1852
1853 real(dp) :: lat_1, long_1, lat_2, long_2
1854
1855 type(grid), pointer :: level0_iDomain => null()
1856
1857
1858 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
1859
1860 ! regular X-Y cordinate system
1861 IF(icoorsystem .EQ. 0) THEN
1862
1863 select case (fdir)
1864 case(1, 4, 16, 64) ! E, S, W, N
1865 length = 1.0_dp
1866 case(2, 8, 32, 128) ! SE, SW, NW, NE
1867 length = sqrt2_dp
1868 end select
1869 length = length * level0_idomain%cellsize
1870
1871 ! regular lat-lon cordinate system
1872 ELSE IF(icoorsystem .EQ. 1) THEN
1873 irow_to = irow
1874 jcol_to = jcol
1875
1876 ! move in the direction of flow
1877 call movedownonecell(fdir, irow_to, jcol_to)
1878
1879 ! estimate lat-lon points
1880 lat_1 = level0_idomain%yllcorner + real((level0_idomain%ncols - jcol), dp) * level0_idomain%cellsize + &
1881 0.5_dp * level0_idomain%cellsize
1882 long_1 = level0_idomain%xllcorner + real((irow - 1), dp) * level0_idomain%cellsize + &
1883 0.5_dp * level0_idomain%cellsize
1884
1885 lat_2 = level0_idomain%yllcorner + real((level0_idomain%ncols - jcol_to), dp) * level0_idomain%cellsize + &
1886 0.5_dp * level0_idomain%cellsize
1887 long_2 = level0_idomain%xllcorner + real((irow_to - 1), dp) * level0_idomain%cellsize + &
1888 0.5_dp * level0_idomain%cellsize
1889 ! get distance between two points
1890 call get_distance_two_lat_lon_points(lat_1, long_1, lat_2, long_2, length)
1891
1892 END IF
1893 !
1894 end subroutine celllength
1895
1896
1897 ! --------------------------------------------------------------------------
1898
1899 ! NAME
1900 ! get_distance_two_lat_lon_points
1901
1902 ! PURPOSE
1903 !> \brief estimate distance in [m] between two points in a lat-lon
1904
1905 !> \details estimate distance in [m] between two points in a lat-lon
1906 !> Code is based on one that is implemented in the VIC-3L model
1907
1908 ! INTENT(IN)
1909 !> \param[in] "real(dp) :: lat1, long1, lat2, long2" latitude of point-1
1910 !> \param[in] "real(dp) :: lat1, long1, lat2, long2" longitude of point-1
1911 !> \param[in] "real(dp) :: lat1, long1, lat2, long2" latitude of point-2
1912 !> \param[in] "real(dp) :: lat1, long1, lat2, long2" longitude of point-2
1913
1914 ! INTENT(OUT)
1915 !> \param[out] "real(dp) :: distance_out" distance between two points [m]
1916
1917 ! HISTORY
1918 !> \authors Rohini Kumar
1919
1920 !> \date May 2014
1921
1922 ! Modifications:
1923 ! Stephan Thober Aug 2015 - ported to mRM
1924 ! Robert Schweppe Jun 2018 - refactoring and reformatting
1925
1926 subroutine get_distance_two_lat_lon_points(lat1, long1, lat2, long2, distance_out)
1927
1928 use mo_constants, only : radiusearth_dp, twopi_dp
1929
1930 implicit none
1931
1932 ! longitude of point-2
1933 real(dp), intent(in) :: lat1, long1, lat2, long2
1934
1935 ! distance between two points [m]
1936 real(dp), intent(out) :: distance_out
1937
1938 real(dp) :: theta1
1939
1940 real(dp) :: phi1
1941
1942 real(dp) :: theta2
1943
1944 real(dp) :: phi2
1945
1946 real(dp) :: dtor
1947
1948 real(dp) :: term1
1949
1950 real(dp) :: term2
1951
1952 real(dp) :: term3
1953
1954 real(dp) :: temp
1955
1956
1957 dtor = twopi_dp / 360.0_dp
1958 theta1 = dtor * long1
1959 phi1 = dtor * lat1
1960 theta2 = dtor * long2
1961 phi2 = dtor * lat2
1962
1963 term1 = cos(phi1) * cos(theta1) * cos(phi2) * cos(theta2)
1964 term2 = cos(phi1) * sin(theta1) * cos(phi2) * sin(theta2)
1965 term3 = sin(phi1) * sin(phi2)
1966 temp = term1 + term2 + term3
1967 if(temp .GT. 1.0_dp) temp = 1.0_dp
1968
1969 distance_out = radiusearth_dp * acos(temp);
1970
1971 end subroutine get_distance_two_lat_lon_points
1972
1973 ! --------------------------------------------------------------------------
1974
1975 ! NAME
1976 ! L11_flow_accumulation
1977
1978 ! PURPOSE
1979
1980 !> \brief Calculates L11 flow accumulation per grid cell
1981 !> \details Calculates L11 flow accumulation per grid cell using L11_fDir
1982 !> and L11_cellarea. L11_flow_accumulation contains the recursiv subroutine
1983 !> calculate_L11_flow_accumulation
1984
1985 ! INTENT(IN)
1986 !> iDomain
1987
1988 ! INTENT(INOUT)
1989 ! None
1990
1991 ! INTENT(OUT)
1992 ! L11_fAcc, L11_LinkIn_fAcc
1993
1994 ! INTENT(IN), OPTIONAL
1995 ! None
1996
1997 ! INTENT(INOUT), OPTIONAL
1998 ! None
1999
2000 ! INTENT(OUT), OPTIONAL
2001 ! None
2002
2003 ! RETURN
2004 ! None
2005
2006 ! RESTRICTIONS
2007 ! None
2008
2009 ! EXAMPLE
2010 ! None
2011
2012 ! LITERATURE
2013 ! None
2014
2015 ! HISTORY
2016 !> \author Matthias Kelbling
2017 !> \date Aug 2017
2018 ! Modified
2019 ! Stephan Thober, Jun 2018 - refactored to fit MPR_extract
2020
2021 ! --------------------------------------------------------------------------
2022 subroutine l11_flow_accumulation(iDomain)
2023
2024 use mo_mrm_global_variables, only: &
2025 l11_fdir, & ! IN: flow direction at L11 (standard notation)
2026 l11_facc ! OUT: flow accumulation at L11 [km^2]
2029 use mo_append, only : append
2030
2031 implicit none
2032
2033 ! local variables
2034 integer(i4), intent(in) :: idomain
2035 real(dp), dimension(:,:), allocatable :: facc11 ! L11_fAcc array
2036 integer(i4) :: ii, jj ! row and col index
2037 integer(i4) :: s11, e11 ! Vec. position iDomain - fAcc
2038 integer(i4) :: nrows11, ncols11 ! array size Domain
2039 integer(i4), dimension(:,:), allocatable :: fdir11 ! L11_fDir array
2040 logical, dimension(:,:), allocatable :: mask11 ! Domain mask
2041
2042 ! initialize Domain info
2043 nrows11 = level11(idomain)%nrows
2044 ncols11 = level11(idomain)%ncols
2045 s11 = level11(idomain)%iStart
2046 e11 = level11(idomain)%iEnd
2047 mask11 = level11(idomain)%mask
2048
2049 ! allocate arrays
2050 allocate(fdir11(nrows11, ncols11))
2051 allocate(facc11(nrows11, ncols11))
2052
2053 ! initialize
2054 fdir11(:,:) = nodata_i4
2055
2056 ! get data
2057 fdir11(:,:) = unpack( l11_fdir(s11:e11), mask11, nodata_i4 )
2058
2059 ! initialize fAcc11 with cell area
2060 facc11 = unpack( level11(idomain)%cellarea * 1.e-6, mask11, nodata_dp )
2061
2062 ! For each sink call "calculate_L11_flow_accumulation"
2063 do jj=1, ncols11
2064 do ii=1, nrows11
2065 if (fdir11(ii,jj) .eq. 0) then
2066 call calculate_l11_flow_accumulation(fdir = fdir11, &
2067 facc = facc11, &
2068 ii = ii, &
2069 jj = jj, &
2070 nrow = nrows11, &
2071 ncol = ncols11)
2072 end if
2073 end do
2074 end do
2075
2076 ! Append
2077 call append( l11_facc, pack(facc11(:,:),mask11))
2078
2079 ! free space
2080 deallocate(fdir11, facc11, mask11)
2081
2082 contains
2083
2084 recursive subroutine calculate_l11_flow_accumulation(fDir, fAcc, ii, jj, nrow, ncol)
2085
2086 implicit none
2087
2088 integer(i4), intent(in) :: fdir(:,:) ! flow Direction
2089 real(dp), intent (inout) :: facc(:,:) ! flow accumulation
2090 integer(i4), intent(in) :: ii, jj ! row and col index
2091 integer(i4), intent(in) :: nrow, ncol ! number of rows,cols in array
2092
2093 ! Scan order:
2094 !
2095 ! 6 7 8
2096 ! 5 x 1
2097 ! 4 3 2
2098 !
2099 ! Each:
2100 ! 1. Check if cell is inflow cell to current grid
2101 ! 2. If yes: Call calculate_subroutine and add result
2102
2103 if (jj+1 .le. ncol) then
2104 if (fdir(ii,jj+1) .eq. 16_i4) then
2105 call calculate_l11_flow_accumulation(fdir, facc, ii, jj+1, nrow, ncol)
2106 facc(ii,jj) = facc(ii,jj) + facc(ii,jj+1)
2107 end if
2108 end if
2109
2110 if ((ii+1 .le. nrow) .and. (jj+1 .le. ncol)) then
2111 if (fdir(ii+1,jj+1) .eq. 32_i4) then
2112 call calculate_l11_flow_accumulation(fdir, facc, ii+1, jj+1, nrow, ncol)
2113 facc(ii,jj) = facc(ii,jj) + facc(ii+1,jj+1)
2114 end if
2115 end if
2116
2117 if (ii+1 .le. nrow) then
2118 if (fdir(ii+1,jj) .eq. 64_i4) then
2119 call calculate_l11_flow_accumulation(fdir, facc, ii+1, jj, nrow, ncol)
2120 facc(ii,jj) = facc(ii,jj) + facc(ii+1,jj)
2121 end if
2122 end if
2123
2124 if ((ii+1 .le. nrow) .and. (jj-1 .ge. 1)) then
2125 if (fdir(ii+1,jj-1) .eq. 128_i4) then
2126 call calculate_l11_flow_accumulation(fdir, facc, ii+1, jj-1, nrow, ncol)
2127 facc(ii,jj) = facc(ii,jj) + facc(ii+1,jj-1)
2128 end if
2129 end if
2130
2131 if (jj-1 .ge. 1) then
2132 if (fdir(ii,jj-1) .eq. 1_i4) then
2133 call calculate_l11_flow_accumulation(fdir, facc, ii, jj-1, nrow, ncol)
2134 facc(ii,jj) = facc(ii,jj) + facc(ii,jj-1)
2135 end if
2136 end if
2137
2138 if ((ii-1 .ge. 1) .and. (jj-1 .ge. 1)) then
2139 if (fdir(ii-1,jj-1) .eq. 2_i4) then
2140 call calculate_l11_flow_accumulation(fdir, facc, ii-1, jj-1, nrow, ncol)
2141 facc(ii,jj) = facc(ii,jj) + facc(ii-1,jj-1)
2142 end if
2143 end if
2144
2145 if (ii-1 .ge. 1) then
2146 if (fdir(ii-1,jj) .eq. 4_i4) then
2147 call calculate_l11_flow_accumulation(fdir, facc, ii-1, jj, nrow, ncol)
2148 facc(ii,jj) = facc(ii,jj) + facc(ii-1,jj)
2149 end if
2150 end if
2151
2152 if ((ii-1 .ge. 1) .and. (jj+1 .le. ncol)) then
2153 if (fdir11(ii-1,jj+1) .eq. 8_i4) then
2154 call calculate_l11_flow_accumulation(fdir, facc, ii-1, jj+1, nrow, ncol)
2155 facc(ii,jj) = facc(ii,jj) + facc(ii-1,jj+1)
2156 end if
2157 end if
2158
2159 ! print *, ii, jj, nrow, ncol, fAcc(ii, jj)
2160
2161 end subroutine calculate_l11_flow_accumulation
2162
2163 end subroutine l11_flow_accumulation
2164
2165
2166 ! ------------------------------------------------------------------
2167
2168 ! NAME
2169 ! L11_calc_celerity
2170
2171 ! PURPOSE
2172 !> \brief L11 celerity based on L0 elevation and L0 fAcc
2173
2174 !> \details L11 celerity based on L0 elevation and L0 fAcc
2175
2176 ! INTENT(IN)
2177 !> \param[in]
2178
2179 ! INTENT(INOUT)
2180 ! None
2181
2182 ! INTENT(OUT)
2183 ! None
2184
2185 ! INTENT(IN), OPTIONAL
2186 ! None
2187
2188 ! INTENT(INOUT), OPTIONAL
2189 ! None
2190
2191 ! INTENT(OUT), OPTIONAL
2192 ! None
2193
2194 ! RETURN
2195 ! None
2196
2197 ! RESTRICTIONS
2198 ! None
2199
2200 ! EXAMPLE
2201 ! None
2202
2203 ! LITERATURE
2204 ! None
2205
2206 ! HISTORY
2207 !> \author Matthias Kelbling
2208 !> \date Oct 2017
2209
2210 ! ------------------------------------------------------------------
2211
2212 subroutine l11_calc_celerity(iDomain, param)
2214 use mo_mad, only: mad
2215 use mo_append, only: append
2216 use mo_mpr_global_variables, only: &
2217 l0_slope ! IN: slope [%]
2218 use mo_common_types, only: grid
2219 use mo_common_variables, only: &
2220 domainmeta, & ! IN: for L0 Domain indexer
2221 level0 ! IN: level 0 grid
2222 use mo_mrm_global_variables, only: &
2223 l0_fdir, & ! IN: flow direction (standard notation) L0
2224 l0_facc, & ! IN: flow accumulation (number of cells)?
2225 l0_streamnet, & ! IN: stream Network at Level 0
2226 level11, & ! IN: level 11 grid
2227 l11_frow, & ! IN: from row in L0 grid
2228 l11_fcol, & ! IN: from col in L0 grid
2229 l11_trow, & ! IN: to row in L0 grid
2230 l11_tcol, & ! IN: to col in L0 grid
2231 l11_netperm, & ! IN: routing order (permutation)
2232 l11_noutlets, & ! IN: Number of Outlets/Sinks
2233 l11_celerity, & ! INOUT: averaged celerity
2234 l0_celerity ! INOUT
2235
2236 implicit none
2237
2238 integer(i4), intent(in) :: idomain ! Domain
2239 real(dp), dimension(:), intent(in) :: param
2240
2241 ! local
2242 integer(i4) :: ncells0
2243 integer(i4) :: nnodes
2244 integer(i4) :: nlinks
2245 integer(i4) :: nrows0, ncols0
2246 integer(i4) :: istart0, iend0
2247 integer(i4) :: nrows11, ncols11
2248 integer(i4) :: istart11, iend11
2249 logical, dimension(:,:), allocatable :: mask0
2250 integer(i4), dimension(:,:), allocatable :: id0
2251 integer(i4), dimension(:,:), allocatable :: fdir0
2252 integer(i4), dimension(:,:), allocatable :: facc0
2253 real(dp), dimension(:,:), allocatable :: slope0
2254 real(dp), dimension(:), allocatable :: slope_tmp
2255 real(dp), dimension(:,:), allocatable :: cellarea0
2256 integer(i4), dimension(:), allocatable :: netperm ! routing order (permutation)
2257 integer(i4), dimension(:), allocatable :: nlinkfromrow
2258 integer(i4), dimension(:), allocatable :: nlinkfromcol
2259 integer(i4), dimension(:), allocatable :: nlinktorow
2260 integer(i4), dimension(:), allocatable :: nlinktocol
2261 integer(i4) :: ii, rr, ns
2262 integer(i4) :: frow, fcol
2263 integer(i4) :: fid, tid
2264 real(dp), dimension(:), allocatable :: stack, append_chunk ! Stacks celerity along the L0 river-path
2265 integer(i4), dimension(:), allocatable :: dummy_1d
2266
2267 real(dp) :: l0_link_slope
2268 real(dp), dimension(:), allocatable :: celerity11
2269 real(dp), dimension(:,:), allocatable :: celerity0
2270
2271 integer(i4), dimension(:,:), allocatable :: nodata_i4_tmp
2272 real(dp), dimension(:,:), allocatable :: nodata_dp_tmp
2273 logical, dimension(:), allocatable :: slopemask0
2274
2275 type(grid), pointer :: level0_idomain
2276
2277 ! level-0 information
2278 level0_idomain => level0(domainmeta%L0DataFrom(idomain))
2279 nrows0 = level0_idomain%nrows
2280 ncols0 = level0_idomain%ncols
2281 ncells0 = level0_idomain%ncells
2282 istart0 = level0_idomain%iStart
2283 iend0 = level0_idomain%iEnd
2284 mask0 = level0_idomain%mask
2285
2286 ! level-11 information
2287 istart11 = level11(idomain)%iStart
2288 iend11 = level11(idomain)%iEnd
2289 nrows11 = level11(idomain)%nrows
2290 ncols11 = level11(idomain)%ncols
2291 nnodes = level11(idomain)%ncells
2292
2293 nlinks = nnodes - l11_noutlets(idomain)
2294
2295 ! allocate
2296 allocate ( id0( nrows0, ncols0 ) )
2297 allocate ( slope0( nrows0, ncols0 ) )
2298 allocate ( fdir0( nrows0, ncols0 ) )
2299 allocate ( facc0( nrows0, ncols0 ) )
2300 allocate ( cellarea0( nrows0, ncols0 ) )
2301 allocate ( celerity0( nrows0, ncols0 ) )
2302 allocate ( slopemask0( ncells0 ) )
2303
2304 ! Routing network vectors have nNodes size instead of nLinks to
2305 ! avoid the need of having two extra indices to identify a Domain.
2306 allocate ( stack( 1 ) )
2307 allocate ( append_chunk( 1 ) )
2308 allocate ( dummy_1d( 2 ))
2309 allocate ( netperm( nnodes ) )
2310 allocate ( nlinkfromrow( nnodes ) )
2311 allocate ( nlinkfromcol( nnodes ) )
2312 allocate ( nlinktorow( nnodes ) )
2313 allocate ( nlinktocol( nnodes ) )
2314 allocate ( celerity11( nnodes ) )
2315 allocate ( slope_tmp( nnodes ) )
2316
2317 allocate (nodata_i4_tmp( nrows0, ncols0 ) )
2318 allocate (nodata_dp_tmp( nrows0, ncols0 ) )
2319
2320 ! initialize
2321 id0(:,:) = nodata_i4
2322 fdir0(:,:) = nodata_i4
2323 facc0(:,:) = nodata_i4
2324 cellarea0(:,:) = nodata_dp
2325 slope0(:,:) = nodata_dp
2326
2327 stack(:) = nodata_dp
2328 append_chunk(:) = nodata_dp
2329 netperm(:) = nodata_i4
2330 nlinkfromrow(:) = nodata_i4
2331 nlinkfromcol(:) = nodata_i4
2332 nlinktorow(:) = nodata_i4
2333 nlinktocol(:) = nodata_i4
2334 celerity11(:) = nodata_dp
2335 celerity0(:,:) = nodata_dp
2336 slopemask0(:) = .false.
2337
2338 nodata_i4_tmp(:,:) = nodata_i4
2339 nodata_dp_tmp(:,:) = nodata_dp
2340
2341 ! for a single node model run
2342
2343 if(nnodes .GT. 1) then
2344 ! get L0 fields
2345 id0(:,:) = unpack(level0_idomain%Id(1:ncells0), mask0, nodata_i4_tmp)
2346 fdir0(:,:) = unpack(l0_fdir(istart0:iend0), mask0, nodata_i4_tmp)
2347 facc0(:,:) = unpack(l0_facc(istart0:iend0), mask0, nodata_i4_tmp)
2348 cellarea0(:,:) = unpack(level0_idomain%cellarea(1:ncells0), mask0, nodata_dp_tmp)
2349
2350 ! smoothing river slope
2351 slope_tmp = l0_slope(istart0:iend0)
2352 where ( slope_tmp .lt. 0.1_dp ) slope_tmp = 0.1_dp
2353
2354 slopemask0(:) = (l0_streamnet(istart0:iend0) .ne. nodata_i4)
2355
2356 ! smooth river cells if there is more than one cell
2357 if( count(slopemask0) .GT. 1) then
2358 slope_tmp = mad(arr = slope_tmp, z = 2.25_dp, mask = slopemask0, tout='u', mval=0.1_dp)
2359 end if
2360 slope0(:,:) = unpack(slope_tmp, mask0, nodata_dp_tmp )
2361
2362 ! get network vectors of L11
2363 netperm(:) = l11_netperm( istart11 : iend11 )
2364 nlinkfromrow(:) = l11_frow( istart11 : iend11 )
2365 nlinkfromcol(:) = l11_fcol( istart11 : iend11 )
2366 nlinktorow(:) = l11_trow( istart11 : iend11 )
2367 nlinktocol(:) = l11_tcol( istart11 : iend11 )
2368
2369 do rr = 1, nlinks
2370
2371 ii = netperm(rr)
2372 frow = nlinkfromrow(ii)
2373 fcol = nlinkfromcol(ii)
2374
2375 ! Init
2376 stack(:) = 0_dp
2377 ns = 1
2378
2379 fid = id0( frow, fcol )
2380 tid = id0( nlinktorow(ii) , nlinktocol(ii) )
2381 do
2382 l0_link_slope = slope0(frow, fcol) / 100._dp
2383
2384 ! celerity parametrization
2385 stack(ns) = param(1) * sqrt(l0_link_slope)
2386
2387 celerity0(frow, fcol) = stack(ns)
2388 ns = ns + 1
2389 fid = id0(frow, fcol)
2390 if( .NOT. (fid == tid)) then
2391 call append(stack, append_chunk)
2392 else
2393 exit
2394 end if
2395 ! move downstream
2396 call movedownonecell( fdir0(frow,fcol), frow, fcol )
2397 end do
2398
2399 celerity11(ii) = size(stack) / sum(1/stack(:))
2400 deallocate(stack)
2401 allocate(stack(1))
2402
2403 end do
2404
2405 else
2406
2407 ! There is only one cell, so no routing is taking place
2408 ! set dummy value of 1 m / s
2409 celerity11(:) = 1._dp
2410
2411 end if
2412
2413 ! Write celerity
2414 l11_celerity(istart11:iend11) = celerity11(:)
2415 l0_celerity(istart0:iend0) = pack(celerity0(:,:), mask0)
2416
2417 ! free space
2418 deallocate (&
2419 mask0, id0, slope_tmp, slopemask0, &
2420 slope0, fdir0, cellarea0, &
2421 stack, netperm, nlinkfromrow, nlinkfromcol, nlinktorow, nlinktocol)
2422
2423 end subroutine l11_calc_celerity
2424
2425end 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
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.