Line data Source code
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
24 : module mo_mrm_net_startup
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
32 : PUBLIC :: L11_set_network_topology
33 : PUBLIC :: L11_routing_order
34 : PUBLIC :: L11_link_location
35 : PUBLIC :: L11_set_drain_outlet_gauges
36 : PUBLIC :: L11_stream_features
37 : PUBLIC :: L11_fraction_sealed_floodplain
38 : PUBLIC :: get_distance_two_lat_lon_points
39 : PUBLIC :: L11_flow_accumulation
40 : PUBLIC :: L11_calc_celerity
41 : contains
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 24 : subroutine L11_L1_mapping(iDomain)
62 :
63 : use mo_append, only : append
64 : use mo_common_constants, only : nodata_i4
65 : use mo_common_variables, only : level1
66 : use mo_mrm_global_variables, only : L11_L1_ID, L1_L11_ID, level11
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 24 : integer(i4), dimension(:, :), allocatable :: L11Id_on_L1
85 :
86 : ! mapping of L1 Id on L11
87 24 : integer(i4), dimension(:, :), allocatable :: L1Id_on_L11
88 :
89 : ! dummy ID
90 24 : integer(i4), dimension(:, :), allocatable :: dummy_2d_id
91 :
92 24 : real(dp) :: cellFactorRbyH
93 :
94 : integer(i4) :: cellFactorRbyH_inv
95 :
96 :
97 24 : nrows1 = level1(iDomain)%nrows
98 24 : nrows11 = level11(iDomain)%nrows
99 24 : ncols1 = level1(iDomain)%ncols
100 24 : ncols11 = level11(iDomain)%ncols
101 :
102 : ! allocate variables for mapping L11 Ids and L1 Ids
103 96 : allocate (L11Id_on_L1 (nrows1, ncols1))
104 96 : allocate (L1Id_on_L11 (nrows11, ncols11))
105 1983 : L11Id_on_L1(:, :) = nodata_i4
106 1812 : L1Id_on_L11(:, :) = nodata_i4
107 :
108 : ! set cell factor for routing
109 24 : 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 24 : if (cellFactorRbyH .lt. 1._dp) then
114 0 : allocate (dummy_2d_id (nrows1, ncols1))
115 0 : dummy_2d_id = unpack(level1(iDomain)%Id, level1(iDomain)%mask, nodata_i4)
116 0 : cellFactorRbyH_inv = int(1. / cellFactorRbyH, i4)
117 0 : kk = 0
118 0 : do jcc = 1, ncols1
119 0 : do icc = 1, nrows1
120 0 : if(.not. level1(iDomain)%mask(icc, jcc)) cycle
121 0 : kk = kk + 1
122 : !
123 0 : iu = (icc - 1) * cellFactorRbyH_inv + 1
124 0 : id = min(icc * cellFactorRbyH_inv, nrows11)
125 0 : jl = (jcc - 1) * cellFactorRbyH_inv + 1
126 0 : jr = min(jcc * cellFactorRbyH_inv, ncols11)
127 :
128 0 : 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 72 : allocate (dummy_2d_id (nrows11, ncols11))
133 24 : dummy_2d_id = unpack(level11(iDomain)%Id, level11(iDomain)%mask, nodata_i4)
134 :
135 24 : kk = 0
136 250 : do jcc = 1, ncols11
137 1812 : do icc = 1, nrows11
138 1562 : if(.not. level11(iDomain)%mask(icc, jcc)) cycle
139 936 : kk = kk + 1
140 :
141 : ! coord. of all corners L11 -> of finer scale level-1
142 936 : iu = (icc - 1) * nint(cellFactorRbyH, i4) + 1
143 936 : id = icc * nint(cellFactorRbyH, i4)
144 936 : jl = (jcc - 1) * nint(cellFactorRbyH, i4) + 1
145 936 : jr = jcc * nint(cellFactorRbyH, i4)
146 :
147 : ! constrain the range of up, down, left, and right boundaries
148 936 : if(iu < 1) iu = 1
149 936 : if(id > nrows1) id = nrows1
150 936 : if(jl < 1) jl = 1
151 936 : if(jr > ncols1) jr = ncols1
152 :
153 : ! Delimitation of level-11 cells on level-1 for L11 resolution lower than L1 resolution
154 3170 : 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 24 : call append(L1_L11_Id, pack (L11Id_on_L1(:, :), level1(iDomain)%mask))
161 : ! L11 data sets
162 24 : call append(L11_L1_Id, pack (L1Id_on_L11(:, :), level11(iDomain)%mask))
163 : ! free space
164 24 : deallocate(L11Id_on_L1, L1Id_on_L11, dummy_2d_id)
165 :
166 24 : 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 24 : subroutine L11_flow_direction(iDomain)
228 :
229 24 : use mo_append, only : append
230 : use mo_common_constants, only : nodata_i4
231 : use mo_common_types, only: Grid
232 : use mo_common_variables, only : domainMeta, level0
233 : use mo_mrm_global_variables, only : L0_draSC, L0_fAcc, L0_fDir, L0_l11_remap, L11_colOut, L11_fDir, &
234 : L11_nOutlets, L11_rowOut, domain_mrm, level11
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 24 : integer(i4), dimension(:, :), allocatable :: iD0
262 :
263 24 : integer(i4), dimension(:, :), allocatable :: fDir0
264 :
265 24 : integer(i4), dimension(:, :), allocatable :: fAcc0
266 :
267 24 : integer(i4), dimension(:, :), allocatable :: fDir11
268 :
269 : ! northing cell loc. of the Outlet
270 24 : integer(i4), dimension(:), allocatable :: rowOut
271 :
272 : ! easting cell loc. of the Outlet
273 24 : integer(i4), dimension(:), allocatable :: colOut
274 :
275 24 : integer(i4), dimension(:, :), allocatable :: draSC0
276 :
277 : ! output location in L0
278 24 : 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 24 : Noutlet = 0_i4
303 24 : 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 24 : nrows0 = level0_iDomain%nrows
317 24 : ncols0 = level0_iDomain%ncols
318 24 : nCells0 = level0_iDomain%ncells
319 24 : nrows11 = level11(iDomain)%nrows
320 24 : ncols11 = level11(iDomain)%ncols
321 24 : nNodes = level11(iDomain)%ncells
322 24 : s0 = level0_iDomain%iStart
323 24 : e0 = level0_iDomain%iEnd
324 :
325 96 : allocate (iD0 (nrows0, ncols0))
326 72 : allocate (fAcc0 (nrows0, ncols0))
327 72 : allocate (fDir0 (nrows0, ncols0))
328 72 : allocate (draSC0 (nrows0, ncols0))
329 96 : allocate (fDir11 (nrows11, ncols11))
330 72 : allocate (rowOut (nNodes))
331 48 : allocate (colOut (nNodes))
332 24 : allocate (oLoc (1, 2))
333 :
334 : ! initialize
335 2862360 : iD0(:, :) = nodata_i4
336 2862360 : fAcc0(:, :) = nodata_i4
337 2862360 : fDir0(:, :) = nodata_i4
338 2862360 : draSC0(:, :) = nodata_i4
339 1812 : fDir11(:, :) = nodata_i4
340 960 : rowOut(:) = nodata_i4
341 960 : colOut(:) = nodata_i4
342 120 : oLoc(:, :) = nodata_i4
343 :
344 :
345 : ! get iD, fAcc, fDir at L0
346 24 : iD0(:, :) = UNPACK(level0_iDomain%Id, level0_iDomain%mask, nodata_i4)
347 24 : fAcc0(:, :) = UNPACK(L0_fAcc (s0 : e0), level0_iDomain%mask, nodata_i4)
348 24 : fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), level0_iDomain%mask, nodata_i4)
349 :
350 : ! case where routing and input data scale is similar
351 24 : IF(nCells0 .EQ. nNodes) THEN
352 0 : oLoc(1, :) = maxloc(fAcc0, level0_iDomain%mask)
353 0 : kk = L0_L11_remap(iDomain)%lowres_id_on_highres(oLoc(1, 1), oLoc(1, 2))
354 : ! for a single node model run
355 0 : if(nCells0 .EQ. 1) then
356 0 : fDir11(1, 1) = fDir0(oLoc(1, 1), oLoc(1, 2))
357 : else
358 0 : fDir11(:, :) = fDir0(:, :)
359 : end if
360 0 : fDir11 (level11(iDomain)%CellCoor(kk, 1), level11(iDomain)%CellCoor(kk, 2)) = 0
361 : ! set location of main outlet in L11
362 0 : do kk = 1, nNodes
363 0 : ii = level11(iDomain)%CellCoor(kk, 1)
364 0 : jj = level11(iDomain)%CellCoor(kk, 2)
365 0 : rowOut(kk) = ii
366 0 : colOut(kk) = jj
367 : end do
368 0 : do kk = 1, ncells0
369 0 : ii = level0_iDomain%CellCoor(kk, 1)
370 0 : jj = level0_iDomain%CellCoor(kk, 2)
371 0 : 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 1066064 : do ii = 1, nCells0
380 1066040 : iRow = level0_iDomain%CellCoor(ii, 1)
381 1066040 : jCol = level0_iDomain%CellCoor(ii, 2)
382 1066040 : call moveDownOneCell(fDir0(iRow, jCol), iRow, jCol)
383 : ! check whether new location is inside bound
384 1066040 : is_outlet = .False.
385 : if ((iRow .le. 0_i4) .or. (iRow .gt. nrows0) .or. &
386 1066040 : (jCol .le. 0_i4) .or. (jCol .gt. ncols0)) then
387 : is_outlet = .True.
388 : else
389 1066040 : if (fdir0(iRow, jCol) .le. 0) is_outlet = .True.
390 : end if
391 : !
392 1066064 : if (is_outlet) then
393 24 : Noutlet = Noutlet + 1_i4
394 : ! cell is an outlet
395 24 : if (Noutlet .eq. 1) then
396 72 : oLoc(1, :) = level0_iDomain%CellCoor(ii, :)
397 : else
398 0 : call append(oLoc, level0_iDomain%CellCoor(ii : ii, :))
399 : end if
400 : ! drain this cell into corresponding L11 cell
401 24 : kk = L0_L11_remap(iDomain)%lowres_id_on_highres(oLoc(Noutlet, 1), oLoc(Noutlet, 2))
402 24 : draSC0(oLoc(Noutlet, 1), oLoc(Noutlet, 2)) = kk
403 : ! check whether cell has maximum flow accumulation
404 : ! coord. of all corners
405 24 : iu = l0_l11_remap(iDomain)%upper_bound (kk)
406 24 : id = l0_l11_remap(iDomain)%lower_bound (kk)
407 24 : jl = l0_l11_remap(iDomain)%left_bound (kk)
408 24 : jr = l0_l11_remap(iDomain)%right_bound(kk)
409 52992 : if (maxval(facc0(iu : id, jl : jr)) .eq. facc0(oLoc(Noutlet, 1), oLoc(Noutlet, 2))) then
410 : ! set location of outlet at L11
411 24 : rowOut(kk) = oLoc(Noutlet, 1)
412 24 : colOut(kk) = oLoc(Noutlet, 2)
413 24 : 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 960 : do kk = 1, nNodes
421 :
422 : ! exclude outlet L11
423 936 : if (rowOut(kk) > 0) cycle
424 :
425 912 : ic = level11(iDomain)%CellCoor(kk, 1)
426 912 : jc = level11(iDomain)%CellCoor(kk, 2)
427 :
428 : ! coord. of all corners
429 912 : iu = l0_l11_remap(iDomain)%upper_bound (kk)
430 912 : id = l0_l11_remap(iDomain)%lower_bound (kk)
431 912 : jl = l0_l11_remap(iDomain)%left_bound (kk)
432 912 : jr = l0_l11_remap(iDomain)%right_bound(kk)
433 :
434 912 : fAccMax = -9
435 912 : idMax = 0
436 912 : side = -1
437 : ! searching on side 4
438 39504 : do jj = jl, jr
439 38592 : if ((fAcc0(iu, jj) > fAccMax) .and. &
440 : (fDir0(iu, jj) == 32 .or. &
441 : fDir0(iu, jj) == 64 .or. &
442 912 : fDir0(iu, jj) == 128)) then
443 1858 : fAccMax = fAcc0(iu, jj)
444 1858 : idMax = id0(iu, jj)
445 1858 : side = 4
446 : end if
447 : end do
448 :
449 : ! searching on side 1
450 39504 : do ii = iu, id
451 38592 : if ((fAcc0(ii, jr) > fAccMax) .and. &
452 : (fDir0(ii, jr) == 1 .or. &
453 : fDir0(ii, jr) == 2 .or. &
454 912 : fDir0(ii, jr) == 128)) then
455 418 : fAccMax = fAcc0(ii, jr)
456 418 : idMax = id0(ii, jr)
457 418 : side = 1
458 : end if
459 : end do
460 :
461 : ! searching on side 2
462 39504 : do jj = jl, jr
463 38592 : if ((fAcc0(id, jj) > fAccMax) .and. &
464 : (fDir0(id, jj) == 2 .or. &
465 : fDir0(id, jj) == 4 .or. &
466 912 : fDir0(id, jj) == 8)) then
467 546 : fAccMax = fAcc0(id, jj)
468 546 : idMax = id0(id, jj)
469 546 : side = 2
470 : end if
471 : end do
472 :
473 : ! searching on side 3
474 39504 : do ii = iu, id
475 38592 : if ((fAcc0(ii, jl) > fAccMax) .and. &
476 : (fDir0(ii, jl) == 8 .or. &
477 : fDir0(ii, jl) == 16 .or. &
478 912 : fDir0(ii, jl) == 32)) then
479 628 : fAccMax = fAcc0(ii, jl)
480 628 : idMax = id0(ii, jl)
481 628 : side = 3
482 : end if
483 : end do
484 :
485 : ! set location of the cell-outlet (row, col) in L0
486 912 : ii = level0_iDomain%CellCoor(idMax, 1)
487 912 : jj = level0_iDomain%CellCoor(idMax, 2)
488 912 : rowOut(kk) = ii
489 912 : colOut(kk) = jj
490 912 : draSC0(ii, jj) = kk
491 :
492 : ! set fDir at L11
493 936 : if (ii == iu .and. jj == jl) then
494 4 : select case (fDir0(ii, jj))
495 : case (8, 16)
496 0 : fDir11(ic, jc) = 16
497 : case (32)
498 0 : fDir11(ic, jc) = 32
499 : case (64, 128)
500 4 : fDir11(ic, jc) = 64
501 : end select
502 908 : elseif (ii == iu .and. jj == jr) then
503 4 : select case (fDir0(ii, jj))
504 : case (32, 64)
505 2 : fDir11(ic, jc) = 64
506 : case (128)
507 0 : fDir11(ic, jc) = 128
508 : case (1, 2)
509 2 : fDir11(ic, jc) = 1
510 : end select
511 906 : elseif (ii == id .and. jj == jl) then
512 62 : select case (fDir0(ii, jj))
513 : case (2, 4)
514 6 : fDir11(ic, jc) = 4
515 : case (8)
516 2 : fDir11(ic, jc) = 8
517 : case (16, 32)
518 56 : fDir11(ic, jc) = 16
519 : end select
520 850 : elseif (ii == id .and. jj == jr) then
521 2 : select case (fDir0(ii, jj))
522 : case (128, 1)
523 0 : fDir11(ic, jc) = 1
524 : case (2)
525 0 : fDir11(ic, jc) = 2
526 : case (4, 8)
527 2 : fDir11(ic, jc) = 4
528 : end select
529 : else
530 : ! cell on one side
531 50 : select case (side)
532 : case (1)
533 50 : fDir11(ic, jc) = 1
534 : case (2)
535 130 : fDir11(ic, jc) = 4
536 : case (3)
537 354 : fDir11(ic, jc) = 16
538 : case (4)
539 314 : fDir11(ic, jc) = 64
540 : case default
541 872 : 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 24 : allocate(domain_mrm(iDomain)%L0_rowOutlet(1))
554 24 : allocate(domain_mrm(iDomain)%L0_colOutlet(1))
555 24 : domain_mrm(iDomain)%L0_Noutlet = nodata_i4
556 72 : domain_mrm(iDomain)%L0_rowOutlet = nodata_i4
557 72 : 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 24 : if (iDomain .eq. 1) then
563 12 : call append(L0_draSC, PACK (draSC0(:, :), level0_iDomain%mask))
564 12 : else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
565 9 : call append(L0_draSC, PACK (draSC0(:, :), level0_iDomain%mask))
566 : end if
567 :
568 24 : domain_mrm(iDomain)%L0_Noutlet = Noutlet
569 : ! set L0 outlet coordinates
570 24 : old_Noutlet = size(domain_mrm(iDomain)%L0_rowOutlet, dim = 1)
571 24 : if (Noutlet .le. old_Noutlet) then
572 48 : domain_mrm(iDomain)%L0_rowOutlet(: Noutlet) = oLoc(:, 1)
573 48 : domain_mrm(iDomain)%L0_colOutlet(: Noutlet) = oLoc(:, 2)
574 : else
575 : ! store up to size of old_Noutlet
576 0 : domain_mrm(iDomain)%L0_rowOutlet(: old_Noutlet) = oLoc(: old_Noutlet, 1)
577 0 : 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 0 : call append(domain_mrm(iDomain)%L0_rowOutlet, oLoc(old_Noutlet + 1 :, 1))
581 0 : call append(domain_mrm(iDomain)%L0_colOutlet, oLoc(old_Noutlet + 1 :, 2))
582 : end if
583 :
584 : ! L11 data sets
585 1812 : call append(L11_nOutlets, count(fdir11 .eq. 0_i4))
586 24 : call append(L11_fDir, PACK (fDir11(:, :), level11(iDomain)%mask))
587 24 : call append(L11_rowOut, rowOut(:))
588 24 : call append(L11_colOut, colOut(:))
589 :
590 : ! communicate
591 24 : call message('')
592 24 : call message(' Domain: ' // num2str(iDomain, '(i3)'))
593 24 : call message(' Number of outlets found at Level 0:.. ' // num2str(Noutlet, '(i7)'))
594 1812 : call message(' Number of outlets found at Level 11:. ' // num2str(count(fdir11 .eq. 0_i4), '(i7)'))
595 :
596 : ! free space
597 24 : deallocate(fDir0, fAcc0, fDir11, rowOut, colOut, draSC0)
598 :
599 24 : 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 24 : subroutine L11_set_network_topology(iDomain)
631 :
632 24 : use mo_append, only : append
633 : use mo_common_constants, only : nodata_i4
634 : use mo_mrm_global_variables, only : L11_fDir, L11_fromN, L11_toN, level11
635 :
636 : implicit none
637 :
638 : ! Domain Id
639 : integer(i4), intent(in) :: iDomain
640 :
641 24 : integer(i4), dimension(:, :), allocatable :: fDir11
642 :
643 24 : integer(i4), dimension(:, :), allocatable :: dummy_2d_id
644 :
645 : integer(i4) :: jj, kk, ic, jc
646 :
647 : integer(i4) :: fn, tn
648 :
649 24 : 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 72 : allocate (nLinkFromN (level11(iDomain)%nCells)) ! valid from (1 : nLinks)
656 48 : allocate (nLinkToN (level11(iDomain)%nCells)) ! "
657 96 : allocate (fDir11 (level11(iDomain)%nrows, level11(iDomain)%ncols))
658 72 : allocate (dummy_2d_id(level11(iDomain)%nrows, level11(iDomain)%ncols))
659 24 : dummy_2d_id = unpack(level11(iDomain)%Id, level11(iDomain)%mask, nodata_i4)
660 :
661 :
662 : ! initialize
663 960 : nLinkFromN(:) = nodata_i4
664 960 : nLinkToN(:) = nodata_i4
665 1812 : fDir11(:, :) = nodata_i4
666 :
667 : ! get grids of L11
668 24 : fDir11(:, :) = UNPACK(L11_fDir (level11(iDomain)%iStart : level11(iDomain)%iEnd), level11(iDomain)%mask, nodata_i4)
669 :
670 : ! ------------------------------------------------------------------
671 : ! network topology
672 : ! ------------------------------------------------------------------
673 :
674 24 : jj = 0
675 960 : do kk = 1, level11(iDomain)%nCells
676 936 : ic = level11(iDomain)%CellCoor(kk, 1)
677 936 : jc = level11(iDomain)%CellCoor(kk, 2)
678 936 : fn = kk
679 936 : call moveDownOneCell(fDir11(ic, jc), ic, jc)
680 936 : tn = dummy_2d_id(ic, jc)
681 936 : if (fn == tn) cycle
682 912 : jj = jj + 1
683 912 : nLinkFromN(jj) = fn
684 1872 : nLinkToN(jj) = tn
685 : end do
686 :
687 : !--------------------------------------------------------
688 : ! Start padding up local variables to global variables
689 : !--------------------------------------------------------
690 :
691 : ! L11 data sets
692 24 : call append(L11_fromN, nLinkFromN(:)) ! sinks are at the end
693 24 : call append(L11_toN, nLinkToN(:))
694 :
695 : ! free space
696 24 : deallocate (fDir11, nLinkFromN, nLinkToN)
697 :
698 24 : 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 24 : subroutine L11_routing_order(iDomain)
729 :
730 24 : use mo_append, only : append
731 : use mo_common_constants, only : nodata_i4
732 : use mo_mrm_global_variables, only : L11_fDir, L11_fromN, L11_label, L11_nOutlets, L11_netPerm, L11_rOrder, L11_sink, L11_toN, &
733 : level11, sink_cells
734 :
735 : implicit none
736 :
737 : ! Domain Id
738 : integer(i4), intent(in) :: iDomain
739 :
740 : integer(i4) :: nLinks
741 :
742 : ! from node
743 24 : integer(i4), dimension(:), allocatable :: nLinkFromN
744 :
745 : ! to node
746 24 : integer(i4), dimension(:), allocatable :: nLinkToN
747 :
748 : ! network routing order
749 24 : integer(i4), dimension(:), allocatable :: nLinkROrder
750 :
751 : ! label Id [0='', 1=HeadWater, 2=Sink]
752 24 : integer(i4), dimension(:), allocatable :: nLinkLabel
753 :
754 : ! == .true. if sink node reached
755 24 : logical, dimension(:), allocatable :: nLinkSink
756 :
757 : ! routing order (permutation)
758 24 : integer(i4), dimension(:), allocatable :: netPerm
759 :
760 : integer(i4) :: ii, jj, kk
761 :
762 : logical :: flag
763 :
764 :
765 24 : 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 72 : allocate (nLinkFromN (level11(iDomain)%nCells)) ! all vectors valid from (1 : nLinks)
771 48 : allocate (nLinkToN (level11(iDomain)%nCells))
772 48 : allocate (nLinkROrder (level11(iDomain)%nCells))
773 48 : allocate (nLinkLabel (level11(iDomain)%nCells))
774 48 : allocate (nLinkSink (level11(iDomain)%nCells))
775 48 : allocate (netPerm (level11(iDomain)%nCells))
776 : ! initialize
777 960 : nLinkFromN(:) = nodata_i4
778 960 : nLinkToN(:) = nodata_i4
779 960 : nLinkROrder(1 : nLinks) = 1
780 24 : nLinkROrder(level11(iDomain)%nCells) = nodata_i4
781 936 : nLinkLabel(1 : nLinks) = 0
782 24 : nLinkLabel(level11(iDomain)%nCells) = nodata_i4
783 960 : nLinkSink(:) = .FALSE.
784 960 : netPerm(:) = nodata_i4
785 :
786 : ! for a single node model run
787 24 : if(level11(iDomain)%nCells .GT. 1) then
788 : ! get network vectors of L11
789 960 : nLinkFromN(:) = L11_fromN (level11(iDomain)%iStart : level11(iDomain)%iEnd)
790 960 : nLinkToN(:) = L11_toN (level11(iDomain)%iStart : level11(iDomain)%iEnd)
791 :
792 936 : loop1 : do ii = 1, nLinks
793 31812 : loop2 : do jj = 1, nLinks
794 31474 : if (jj == ii) cycle loop2
795 30714 : if (nLinkFromN(ii) == nLinkToN(jj)) then
796 598 : nLinkROrder(ii) = -9
797 : end if
798 31028 : if (nLinkROrder(ii) == -9) cycle loop1
799 : end do loop2
800 : end do loop1
801 : ! counting headwaters
802 : kk = 0
803 936 : do ii = 1, nLinks
804 936 : if (nLinkROrder(ii) == 1) then
805 314 : kk = kk + 1
806 314 : nLinkROrder(ii) = kk
807 314 : nLinkLabel(ii) = 1 ! 'Head Water'
808 : end if
809 : end do
810 : ! counting downstream
811 11748 : do while (minval(nLinkROrder(1 : nLinks)) < 0)
812 : !! print *, count(nLinkROrder .lt. 0), minval(nLinkROrder)
813 10836 : loop3 : do ii = 1, nLinks
814 10584 : if (.NOT. nLinkROrder(ii) == -9) cycle loop3
815 76524 : flag = .TRUE.
816 76524 : loop4 : do jj = 1, nLinks
817 76524 : if (jj == ii .OR. nLinkFromN(ii) /= nLinkToN(jj)) then
818 : cycle loop4
819 3716 : 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 2508 : if (flag) then
827 598 : kk = kk + 1
828 598 : nLinkROrder(ii) = kk
829 : end if
830 : end do loop3
831 : end do
832 :
833 : ! identify sink cells
834 936 : do ii = 1, nLinks
835 936 : if (L11_fdir(level11(iDomain)%iStart + nLinkToN(ii) - 1_i4) .eq. 0_i4) then
836 26 : nlinksink(ii) = .True.
837 : ! add sink target cell to sink_cells if not already present (in case of multiple inflows (rare case))
838 98 : 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 960 : where(nlinksink) nLinkLabel = 2 ! 'Sink'
842 :
843 : ! keep routing order
844 936 : do ii = 1, nLinks
845 936 : 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 24 : call append(L11_rOrder, nLinkROrder(:))
856 24 : call append(L11_label, nLinkLabel(:))
857 24 : call append(L11_sink, nLinkSink(:))
858 24 : call append(L11_netPerm, netPerm(:))
859 :
860 : ! free space
861 24 : deallocate (nLinkFromN, nLinkToN, nLinkROrder, nLinkLabel, nLinkSink, netPerm)
862 :
863 24 : 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 24 : subroutine L11_link_location(iDomain)
892 :
893 24 : use mo_append, only : append
894 : use mo_common_constants, only : nodata_i4
895 : use mo_common_types, only: Grid
896 : use mo_common_variables, only : domainMeta, level0
897 : use mo_mrm_global_variables, only : L0_draSC, L0_fDir, L11_colOut, L11_fCol, L11_fRow, L11_fromN, &
898 : L11_nOutlets, L11_netPerm, L11_rowOut, L11_tCol, L11_tRow, domain_mrm, level11
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 24 : integer(i4), dimension(:), allocatable :: rowOut
910 :
911 : ! easting cell loc. of the Outlet
912 24 : integer(i4), dimension(:), allocatable :: colOut
913 :
914 24 : integer(i4), dimension(:), allocatable :: nLinkFromN
915 :
916 24 : integer(i4), dimension(:), allocatable :: netPerm
917 :
918 24 : integer(i4), dimension(:), allocatable :: nLinkFromRow
919 :
920 24 : integer(i4), dimension(:), allocatable :: nLinkFromCol
921 :
922 24 : integer(i4), dimension(:), allocatable :: nLinkToRow
923 :
924 24 : integer(i4), dimension(:), allocatable :: nLinkToCol
925 :
926 24 : integer(i4), dimension(:, :), allocatable :: fDir0
927 :
928 24 : 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 24 : 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 24 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
947 24 : s0 = level0_iDomain%iStart
948 24 : e0 = level0_iDomain%iEnd
949 :
950 24 : nOutlets = L11_nOutlets(iDomain)
951 :
952 24 : 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 72 : allocate (rowOut (level11(iDomain)%nCells))
958 48 : allocate (colOut (level11(iDomain)%nCells))
959 48 : allocate (nLinkFromN (level11(iDomain)%nCells)) ! all network vectors valid from (1 : nLinks)
960 48 : allocate (netPerm (level11(iDomain)%nCells))
961 48 : allocate (nLinkFromRow (level11(iDomain)%nCells))
962 48 : allocate (nLinkFromCol (level11(iDomain)%nCells))
963 48 : allocate (nLinkToRow (level11(iDomain)%nCells))
964 48 : allocate (nLinkToCol (level11(iDomain)%nCells))
965 96 : allocate (fDir0 (level0_iDomain%nrows, level0_iDomain%ncols))
966 72 : allocate (draSC0 (level0_iDomain%nrows, level0_iDomain%ncols))
967 :
968 : ! initialize
969 984 : rowOut = nodata_i4
970 984 : colOut = nodata_i4
971 984 : nLinkFromN = nodata_i4
972 984 : netPerm = nodata_i4
973 984 : nLinkFromRow = nodata_i4
974 984 : nLinkFromCol = nodata_i4
975 984 : nLinkToRow = nodata_i4
976 984 : nLinkToCol = nodata_i4
977 2862384 : fDir0 = nodata_i4
978 2862384 : draSC0 = nodata_i4
979 :
980 : ! for a single node model run
981 24 : if(level11(iDomain)%nCells .GT. 1) then
982 : ! get fDir at L0
983 24 : fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), level0_iDomain%mask, nodata_i4)
984 24 : draSC0(:, :) = UNPACK(L0_draSC (s0 : e0), level0_iDomain%mask, nodata_i4)
985 :
986 : ! get network vectors of L11
987 960 : nLinkFromN(:) = L11_fromN (level11(iDomain)%iStart : level11(iDomain)%iEnd)
988 960 : netPerm(:) = L11_netPerm (level11(iDomain)%iStart : level11(iDomain)%iEnd)
989 960 : rowOut(:) = L11_rowOut (level11(iDomain)%iStart : level11(iDomain)%iEnd)
990 960 : colOut(:) = L11_colOut (level11(iDomain)%iStart : level11(iDomain)%iEnd)
991 :
992 : ! finding main outlet (row, col) in L0
993 72 : allocate(oLoc(Noutlets, 2))
994 48 : oLoc(:, 1) = domain_mrm(iDomain)%L0_rowOutlet(: Noutlets)
995 48 : oLoc(:, 2) = domain_mrm(iDomain)%L0_colOutlet(: Noutlets)
996 :
997 : ! Location of the stream-joint cells (row, col)
998 936 : do rr = 1, nLinks
999 :
1000 912 : ii = netPerm(rr)
1001 912 : iNode = nLinkFromN(ii)
1002 912 : iRow = rowOut(iNode)
1003 912 : jCol = colOut(iNode)
1004 912 : call moveDownOneCell(fDir0(iRow, jcol), iRow, jcol)
1005 : ! set "from" cell
1006 912 : nLinkFromRow(ii) = iRow
1007 912 : nLinkFromCol(ii) = jCol
1008 :
1009 : ! check whether this location is an outlet
1010 912 : is_outlet = .False.
1011 1824 : do kk = 1, Noutlets
1012 1824 : if (iRow .eq. oLoc(kk, 1) .and. jCol .eq. oLoc(kk, 2)) is_outlet = .True.
1013 : end do
1014 :
1015 1848 : if (is_outlet) then
1016 :
1017 0 : nLinkToRow(ii) = iRow
1018 0 : nLinkToCol(ii) = jCol
1019 :
1020 : else
1021 :
1022 35952 : do while (.not. (draSC0(iRow, jCol) > 0))
1023 35040 : prevRow = iRow
1024 35040 : prevCol = jCol
1025 : call moveDownOneCell(fDir0(iRow, jcol), iRow, jCol)
1026 : ! check whether this location is an outlet and exit
1027 70056 : do kk = 1, Noutlets
1028 70056 : if (iRow .eq. oLoc(kk, 1) .and. jCol .eq. oLoc(kk, 2)) exit
1029 : end do
1030 35952 : 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 0 : num2str(jCol))
1034 : end if
1035 : end do
1036 : ! set "to" cell (when an outlet is reached)
1037 912 : nLinkToRow(ii) = iRow
1038 912 : 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 24 : call append(L11_fRow, nLinkFromRow(:))
1051 24 : call append(L11_fCol, nLinkFromCol(:))
1052 24 : call append(L11_tRow, nLinkToRow(:))
1053 24 : call append(L11_tCol, nLinkToCol(:))
1054 :
1055 : ! free space
1056 : deallocate (rowOut, colOut, nLinkFromN, netPerm, nLinkFromRow, &
1057 24 : nLinkFromCol, nLinkToRow, nLinkToCol, fDir0, draSC0)
1058 :
1059 24 : 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 :
1092 24 : subroutine L11_set_drain_outlet_gauges(iDomain)
1093 :
1094 24 : use mo_append, only : append
1095 : use mo_common_constants, only : nodata_i4
1096 : use mo_common_types, only: Grid
1097 : use mo_common_variables, only : domainMeta, level0
1098 : use mo_mrm_global_variables, only : L0_InflowgaugeLoc, L0_draCell, L0_draSC, L0_fDir, L0_gaugeLoc, domain_mrm, &
1099 : l0_l11_remap
1100 :
1101 : implicit none
1102 :
1103 : ! Domain Id
1104 : integer(i4), intent(in) :: iDomain
1105 :
1106 24 : integer(i4), dimension(:, :), allocatable :: draSC0
1107 :
1108 24 : integer(i4), dimension(:, :), allocatable :: fDir0
1109 :
1110 24 : integer(i4), dimension(:, :), allocatable :: gaugeLoc0
1111 :
1112 24 : integer(i4), dimension(:, :), allocatable :: InflowGaugeLoc0
1113 :
1114 24 : 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 24 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
1126 24 : s0 = level0_iDomain%iStart
1127 24 : e0 = level0_iDomain%iEnd
1128 :
1129 :
1130 : ! allocate
1131 96 : allocate (draSC0 (level0_iDomain%nrows, level0_iDomain%ncols))
1132 72 : allocate (fDir0 (level0_iDomain%nrows, level0_iDomain%ncols))
1133 72 : allocate (gaugeLoc0 (level0_iDomain%nrows, level0_iDomain%ncols))
1134 72 : allocate (InflowGaugeLoc0 (level0_iDomain%nrows, level0_iDomain%ncols))
1135 72 : allocate (draCell0 (level0_iDomain%nrows, level0_iDomain%ncols))
1136 :
1137 : ! initialize
1138 2862360 : draSC0(:, :) = nodata_i4
1139 2862360 : fDir0(:, :) = nodata_i4
1140 2862360 : gaugeLoc0(:, :) = nodata_i4
1141 2862360 : InflowGaugeLoc0(:, :) = nodata_i4
1142 2862360 : draCell0(:, :) = nodata_i4
1143 :
1144 0 : draSC0(:, :) = UNPACK(L0_draSC (s0 : e0), &
1145 24 : level0_iDomain%mask, nodata_i4)
1146 0 : fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), &
1147 24 : level0_iDomain%mask, nodata_i4)
1148 0 : gaugeLoc0(:, :) = UNPACK(L0_gaugeLoc (s0 : e0), &
1149 24 : level0_iDomain%mask, nodata_i4)
1150 0 : InflowGaugeLoc0(:, :) = UNPACK(L0_InflowgaugeLoc (s0 : e0), &
1151 24 : level0_iDomain%mask, nodata_i4)
1152 :
1153 1066064 : do kk = 1, level0_iDomain%nCells
1154 1066040 : ii = level0_iDomain%CellCoor(kk, 1)
1155 1066040 : jj = level0_iDomain%CellCoor(kk, 2)
1156 1066040 : iSc = draSC0(ii, jj)
1157 : ! find drainage path
1158 1066040 : iRow = ii
1159 1066040 : jCol = jj
1160 43644186 : do while (.NOT. iSC > 0)
1161 : ! move downstream
1162 42578146 : call moveDownOneCell(fDir0(iRow, jCol), iRow, jCol)
1163 42578146 : iSC = draSC0(iRow, jCol)
1164 : end do
1165 1066040 : 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 1066040 : if (gaugeLoc0(ii, jj) .NE. nodata_i4) then
1170 : ! evaluation gauges
1171 76 : do ll = 1, domain_mrm(iDomain)%nGauges
1172 : ! search for gaugeID in L0 grid and save ID on L11
1173 76 : if (domain_mrm(iDomain)%gaugeIdList(ll) .EQ. gaugeLoc0(ii, jj)) then
1174 21 : 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 1066064 : if (InflowGaugeLoc0(ii, jj) .NE. nodata_i4) then
1180 : ! inflow gauges
1181 39 : do ll = 1, domain_mrm(iDomain)%nInflowGauges
1182 : ! search for gaugeID in L0 grid and save ID on L11
1183 4 : if (domain_mrm(iDomain)%InflowGaugeIdList(ll) .EQ. InflowGaugeLoc0(ii, jj)) &
1184 37 : 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 24 : if (iDomain .eq. 1) then
1196 12 : call append(L0_draCell, PACK(draCell0(:, :), level0_iDomain%mask))
1197 12 : else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
1198 9 : call append(L0_draCell, PACK(draCell0(:, :), level0_iDomain%mask))
1199 : end if
1200 :
1201 : ! free space
1202 24 : deallocate (draSC0, fDir0, gaugeLoc0, draCell0)
1203 :
1204 24 : 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 24 : subroutine L11_stream_features(iDomain)
1238 :
1239 24 : use mo_append, only : append
1240 : use mo_common_constants, only : nodata_dp, nodata_i4
1241 : use mo_common_types, only: Grid
1242 : use mo_common_variables, only : domainMeta, L0_elev, iFlag_cordinate_sys, level0, processMatrix
1243 : use mo_mrm_global_variables, only : L0_fDir, &
1244 : L0_floodPlain, L0_streamNet, L11_aFloodPlain, L11_fCol, L11_fRow, L11_length, &
1245 : L11_nOutlets, L11_netPerm, L11_slope, L11_tCol, L11_tRow, level11
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 24 : integer(i4), dimension(:, :), allocatable :: iD0
1256 :
1257 24 : integer(i4), dimension(:, :), allocatable :: fDir0
1258 :
1259 24 : real(dp), dimension(:, :), allocatable :: elev0
1260 :
1261 24 : real(dp), dimension(:, :), allocatable :: cellarea0
1262 :
1263 24 : integer(i4), dimension(:, :), allocatable :: streamNet0
1264 :
1265 24 : integer(i4), dimension(:, :), allocatable :: floodPlain0
1266 :
1267 : ! routing order (permutation)
1268 24 : integer(i4), dimension(:), allocatable :: netPerm
1269 :
1270 24 : integer(i4), dimension(:), allocatable :: nLinkFromRow
1271 :
1272 24 : integer(i4), dimension(:), allocatable :: nLinkFromCol
1273 :
1274 24 : integer(i4), dimension(:), allocatable :: nLinkToRow
1275 :
1276 24 : integer(i4), dimension(:), allocatable :: nLinkToCol
1277 :
1278 24 : real(dp), dimension(:), allocatable :: nLinkLength
1279 :
1280 24 : real(dp), dimension(:), allocatable :: nLinkAFloodPlain
1281 :
1282 24 : 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 24 : integer(i4), dimension(:, :), allocatable :: stack, append_chunk
1291 :
1292 24 : integer(i4), dimension(:), allocatable :: dummy_1d
1293 :
1294 24 : real(dp) :: length
1295 :
1296 24 : integer(i4), dimension(:, :), allocatable :: nodata_i4_tmp
1297 :
1298 24 : real(dp), dimension(:, :), allocatable :: nodata_dp_tmp
1299 :
1300 : type(Grid), pointer :: level0_iDomain => null()
1301 :
1302 :
1303 24 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
1304 24 : s0 = level0_iDomain%iStart
1305 24 : e0 = level0_iDomain%iEnd
1306 24 : nLinks = level11(iDomain)%nCells - L11_nOutlets(iDomain)
1307 :
1308 :
1309 : ! allocate
1310 96 : allocate (iD0 (level0_iDomain%nrows, level0_iDomain%ncols))
1311 96 : allocate (elev0 (level0_iDomain%nrows, level0_iDomain%ncols))
1312 96 : allocate (fDir0 (level0_iDomain%nrows, level0_iDomain%ncols))
1313 96 : allocate (cellarea0 (level0_iDomain%nrows, level0_iDomain%ncols))
1314 96 : allocate (streamNet0 (level0_iDomain%nrows, level0_iDomain%ncols))
1315 72 : 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 96 : allocate (stack (level11(iDomain)%nCells, 2)) ! >> stack(nNodes, 2)
1320 24 : allocate (dummy_1d (2))
1321 24 : allocate (append_chunk (8, 2))
1322 72 : allocate (netPerm (level11(iDomain)%nCells))
1323 48 : allocate (nLinkFromRow (level11(iDomain)%nCells))
1324 48 : allocate (nLinkFromCol (level11(iDomain)%nCells))
1325 48 : allocate (nLinkToRow (level11(iDomain)%nCells))
1326 48 : allocate (nLinkToCol (level11(iDomain)%nCells))
1327 72 : allocate (nLinkLength (level11(iDomain)%nCells))
1328 48 : allocate (nLinkAFloodPlain (level11(iDomain)%nCells))
1329 48 : allocate (nLinkSlope (level11(iDomain)%nCells))
1330 :
1331 96 : allocate (nodata_i4_tmp (level0_iDomain%nrows, level0_iDomain%ncols))
1332 96 : allocate (nodata_dp_tmp (level0_iDomain%nrows, level0_iDomain%ncols))
1333 :
1334 : ! initialize
1335 2862360 : iD0(:, :) = nodata_i4
1336 2862360 : elev0(:, :) = nodata_dp
1337 2862360 : fDir0(:, :) = nodata_i4
1338 2862360 : cellarea0(:, :) = nodata_dp
1339 2862360 : streamNet0(:, :) = nodata_i4
1340 2862360 : floodPlain0(:, :) = nodata_i4
1341 :
1342 1944 : stack(:, :) = nodata_i4
1343 456 : append_chunk(:, :) = nodata_i4
1344 960 : netPerm(:) = nodata_i4
1345 960 : nLinkFromRow(:) = nodata_i4
1346 960 : nLinkFromCol(:) = nodata_i4
1347 960 : nLinkToRow(:) = nodata_i4
1348 960 : nLinkToCol(:) = nodata_i4
1349 960 : nLinkLength(:) = nodata_dp
1350 960 : nLinkAFloodPlain(:) = nodata_dp
1351 960 : nLinkSlope(:) = nodata_dp
1352 :
1353 2862360 : nodata_i4_tmp(:, :) = nodata_i4
1354 2862360 : nodata_dp_tmp(:, :) = nodata_dp
1355 :
1356 : ! for a single node model run
1357 24 : if(level11(iDomain)%nCells .GT. 1) then
1358 : ! get L0 fields
1359 24 : iD0(:, :) = UNPACK(level0_iDomain%Id, level0_iDomain%mask, nodata_i4_tmp)
1360 0 : elev0(:, :) = UNPACK(L0_elev (s0 : e0), &
1361 24 : level0_iDomain%mask, nodata_dp_tmp)
1362 0 : fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), &
1363 24 : level0_iDomain%mask, nodata_i4_tmp)
1364 24 : cellarea0(:, :) = UNPACK(level0_iDomain%CellArea, level0_iDomain%mask, nodata_dp_tmp)
1365 :
1366 : ! get network vectors of L11
1367 960 : netPerm(:) = L11_netPerm (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1368 960 : nLinkFromRow(:) = L11_fRow (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1369 960 : nLinkFromCol(:) = L11_fCol (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1370 960 : nLinkToRow(:) = L11_tRow (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1371 960 : nLinkToCol(:) = L11_tCol (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1372 :
1373 : ! Flood plains: stream network delineation
1374 2862360 : streamNet0(:, :) = nodata_i4
1375 2862360 : floodPlain0(:, :) = nodata_i4
1376 :
1377 936 : do rr = 1, nLinks
1378 :
1379 912 : ii = netPerm(rr)
1380 912 : frow = nLinkFromRow(ii)
1381 912 : fcol = nLinkFromCol(ii)
1382 :
1383 : ! Init
1384 912 : streamNet0(frow, fcol) = ii
1385 912 : floodPlain0(frow, fcol) = ii
1386 126872 : stack = 0
1387 18240 : append_chunk = 0
1388 912 : ns = 1
1389 912 : stack(ns, 1) = frow
1390 912 : stack(ns, 2) = fcol
1391 :
1392 912 : call cellLength(iDomain, fDir0(frow, fcol), fRow, fCol, iFlag_cordinate_sys, nLinkLength(ii))
1393 912 : nLinkSlope(ii) = elev0(frow, fcol)
1394 :
1395 912 : fId = iD0(frow, fcol)
1396 912 : tId = iD0(nLinkToRow(ii), nLinkToCol(ii))
1397 :
1398 35952 : 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 35040 : if (processMatrix(8,1) == 1) then
1404 4673474 : do while (ns > 0)
1405 4649970 : if (ns + 8 .gt. size(stack, 1)) then
1406 78 : call append(stack, append_chunk)
1407 : end if
1408 4649970 : call moveUp(elev0, fDir0, frow, fcol, stack, ns)
1409 4649970 : stack(1, 1) = 0
1410 4649970 : stack(1, 2) = 0
1411 : ! stack = cshift(stack, SHIFT = 1, DIM = 1)
1412 : ! substitute cshift <<<
1413 18599880 : dummy_1d = stack(1, :)
1414 1466468700 : stack(: size(stack, dim = 1) - 1, :) = stack(2 :, :)
1415 13949910 : stack(size(stack, dim = 1), :) = dummy_1d
1416 : ! substitute cshift >>>
1417 4649970 : if (stack(1, 1) > 0 .and. stack(1, 2) > 0) floodPlain0(stack(1, 1), stack(1, 2)) = ii
1418 742557794 : ns = count(stack > 0) / 2
1419 : end do
1420 : end if
1421 :
1422 : ! move downstream
1423 35040 : call moveDownOneCell(fDir0(frow, fcol), frow, fcol)
1424 35040 : streamNet0(frow, fcol) = ii
1425 35040 : floodPlain0(frow, fcol) = ii
1426 35040 : fId = iD0(frow, fcol)
1427 4519784 : stack = 0
1428 35040 : ns = 1
1429 35040 : stack(ns, 1) = frow
1430 35040 : stack(ns, 2) = fcol
1431 35040 : call cellLength(iDomain, fDir0(fRow, fCol), fRow, fCol, iFlag_cordinate_sys, length)
1432 35040 : nLinkLength(ii) = nLinkLength(ii) + length
1433 :
1434 : end do
1435 :
1436 : ! stream bed slope
1437 912 : nLinkSlope(ii) = (nLinkSlope(ii) - elev0(frow, fcol)) / nLinkLength(ii)
1438 :
1439 912 : if (nLinkSlope(ii) < 0.0001_dp) nLinkSlope(ii) = 0.0001_dp
1440 :
1441 : ! calculate area of floodplains (avoid overwriting)
1442 111450936 : 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 24 : if ((processMatrix(8, 1) .eq. 2) .or. (processMatrix(8, 1) .eq. 3)) then
1452 280 : if (count(nLinkLength(:) .ge. 0._dp) .gt. 2) then
1453 280 : length = percentile(pack(nLinkLength(:), nLinkLength(:) .ge. 0._dp), 40._dp)
1454 280 : 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 24 : if (iDomain .eq. 1) then
1466 12 : call append(L0_streamNet, PACK (streamNet0(:, :), level0_iDomain%mask))
1467 12 : call append(L0_floodPlain, PACK (floodPlain0(:, :), level0_iDomain%mask))
1468 12 : else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
1469 9 : call append(L0_streamNet, PACK (streamNet0(:, :), level0_iDomain%mask))
1470 9 : call append(L0_floodPlain, PACK (floodPlain0(:, :), level0_iDomain%mask))
1471 : end if
1472 :
1473 :
1474 : ! L11 network data sets
1475 24 : call append(L11_length, nLinkLength(:))
1476 24 : call append(L11_aFloodPlain, nLinkAFloodPlain(:))
1477 24 : call append(L11_slope, nLinkSlope(:))
1478 :
1479 : ! free space
1480 0 : deallocate (&
1481 0 : iD0, elev0, fDir0, streamNet0, floodPlain0, &
1482 0 : cellarea0, stack, netPerm, nLinkFromRow, nLinkFromCol, nLinkToRow, nLinkToCol, &
1483 24 : nLinkLength, nLinkAFloodPlain, nLinkSlope, dummy_1d)
1484 24 : deallocate(nodata_i4_tmp, nodata_dp_tmp)
1485 :
1486 24 : 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 13 : subroutine L11_fraction_sealed_floodplain(LCClassImp, do_init)
1520 :
1521 24 : use mo_append, only : append
1522 : use mo_common_constants, only : nodata_dp
1523 : use mo_common_types, only: Grid
1524 : use mo_common_variables, only : domainMeta, L0_LCover, level0, domainMeta, nLCoverScene
1525 : use mo_mrm_global_variables, only : L0_floodPlain, L11_aFloodPlain, &
1526 : L11_nLinkFracFPimp, L11_nOutlets, level11
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 13 : 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 38 : do iDomain = 1, domainMeta%nDomains
1548 100 : allocate(temp_array(level11(iDomain)%nCells, nLCoverScene))
1549 2040 : temp_array = nodata_dp
1550 25 : if (do_init) then
1551 16 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
1552 :
1553 16 : s0 = level0_iDomain%iStart
1554 16 : e0 = level0_iDomain%iEnd
1555 16 : nLinks = level11(iDomain)%nCells + 1 - L11_nOutlets(iDomain)
1556 16 : nLinkAFloodPlain => L11_aFloodPlain(level11(iDomain)%iStart : level11(iDomain)%iEnd)
1557 :
1558 48 : do iiLC = 1, nLCoverScene
1559 : ! for a single node model run
1560 48 : if(nLinks .GT. 0) then
1561 1360 : do ii = 1, nLinks
1562 2656 : temp_array(ii, iiLC) = sum(level0_iDomain%CellArea(:), &
1563 0 : mask = (L0_floodPlain(s0 : e0) == ii .and. L0_LCover(s0 : e0, iiLC) == LCClassImp)) &
1564 59876256 : / nLinkAFloodPlain(ii)
1565 : end do
1566 : end if
1567 : end do
1568 : end if
1569 25 : call append(L11_nLinkFracFPimp, temp_array(:,:))
1570 38 : deallocate(temp_array)
1571 : end do
1572 :
1573 13 : 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 4649970 : subroutine moveUp(elev0, fDir0, fi, fj, ss, nn)
1605 :
1606 13 : 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 4649970 : ii = ss(1, 1)
1628 4649970 : jj = ss(1, 2)
1629 4649970 : ip = ii + 1
1630 4649970 : im = ii - 1
1631 4649970 : jp = jj + 1
1632 4649970 : jm = jj - 1
1633 :
1634 4649970 : nrows = size(fDir0, 1)
1635 4649970 : ncols = size(fDir0, 2)
1636 :
1637 : !E
1638 4649970 : if (jp <= ncols) then
1639 9299940 : if ((fdir0(ii, jp) == 16) .and. &
1640 18599880 : (le((elev0(ii, jp) - elev0(fi, fj)), deltaH)) .and. &
1641 4649970 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1642 : ) then
1643 970916 : nn = nn + 1
1644 970916 : ss(nn, 1) = ii
1645 970916 : ss(nn, 2) = jp
1646 : !print *, i,jp
1647 : end if
1648 : end if
1649 :
1650 : !SE
1651 4649970 : if ((ip <= nrows) .and. &
1652 : (jp <= ncols)) then
1653 0 : if ((fdir0(ip, jp) == 32) .and. &
1654 9299940 : (le((elev0(ip, jp) - elev0(fi, fj)), deltaH)) .and. &
1655 4649970 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1656 : ) then
1657 321576 : nn = nn + 1
1658 321576 : ss(nn, 1) = ip
1659 321576 : ss(nn, 2) = jp
1660 : !print *, ip,jp
1661 : end if
1662 : end if
1663 :
1664 : !S
1665 4649970 : if ((ip <= nrows) .and. &
1666 : (jp <= ncols)) then
1667 0 : if ((fdir0(ip, jj) == 64) .and. &
1668 9299940 : (le((elev0(ip, jj) - elev0(fi, fj)), deltaH)) .and. &
1669 4649970 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1670 : ) then
1671 787164 : nn = nn + 1
1672 787164 : ss(nn, 1) = ip
1673 787164 : ss(nn, 2) = jj
1674 : !print *, ip,j
1675 : end if
1676 : end if
1677 :
1678 : !SW
1679 : if ((ip <= nrows) .and. &
1680 4649970 : (jp <= ncols) .and. &
1681 : (jm >= 1)) then
1682 0 : if ((fdir0(ip, jm) == 128) .and. &
1683 9299940 : (le((elev0(ip, jm) - elev0(fi, fj)), deltaH)) .and. &
1684 4649970 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1685 : ) then
1686 249878 : nn = nn + 1
1687 249878 : ss(nn, 1) = ip
1688 249878 : ss(nn, 2) = jm
1689 : !print *, ip,jm
1690 : end if
1691 : end if
1692 :
1693 : !W
1694 4649970 : if ((jm >= 1) .and. &
1695 : (jp <= ncols)) then
1696 0 : if ((fdir0(ii, jm) == 1) .and. &
1697 9299940 : (le((elev0(ii, jm) - elev0(fi, fj)), deltaH)) .and. &
1698 4649970 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1699 : ) then
1700 663334 : nn = nn + 1
1701 663334 : ss(nn, 1) = ii
1702 663334 : ss(nn, 2) = jm
1703 : !print *, i,jm
1704 : end if
1705 : end if
1706 :
1707 : !NW
1708 : if ((im >= 1) .and. &
1709 4649970 : (jp <= ncols) .and. &
1710 : (jm >= 1)) then
1711 0 : if ((fdir0(im, jm) == 2) .and. &
1712 9299940 : (le((elev0(im, jm) - elev0(fi, fj)), deltaH)) .and. &
1713 4649970 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1714 : ) then
1715 354012 : nn = nn + 1
1716 354012 : ss(nn, 1) = im
1717 354012 : ss(nn, 2) = jm
1718 : !print *, im,jm
1719 : end if
1720 : end if
1721 :
1722 : !N
1723 4649970 : if ((im >= 1) .and. &
1724 : (jp <= ncols)) then
1725 0 : if ((fdir0(im, jj) == 4) .and. &
1726 9299940 : (le((elev0(im, jj) - elev0(fi, fj)), deltaH)) .and. &
1727 4649970 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1728 : ) then
1729 822956 : nn = nn + 1
1730 822956 : ss(nn, 1) = im
1731 822956 : ss(nn, 2) = jj
1732 : !print *, im,j
1733 : end if
1734 : end if
1735 :
1736 : !NE
1737 4649970 : if ((im >= 1) .and. &
1738 : (jp <= ncols)) then
1739 0 : if ((fdir0(im, jp) == 8) .and. &
1740 9299940 : (le((elev0(im, jp) - elev0(fi, fj)), deltaH)) .and. &
1741 4649970 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1742 : ) then
1743 456630 : nn = nn + 1
1744 456630 : ss(nn, 1) = im
1745 456630 : ss(nn, 2) = jp
1746 : !print *, im,jp
1747 : end if
1748 : end if
1749 :
1750 4649970 : 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 43740628 : 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 45179838 : select case (fDir)
1786 : case(1) !E
1787 1439210 : jCol = jCol + 1
1788 : case(2) !SE
1789 1682387 : iRow = iRow + 1
1790 1682387 : jCol = jCol + 1
1791 : case(4) !S
1792 3804397 : iRow = iRow + 1
1793 : case(8) !SW
1794 5753487 : iRow = iRow + 1
1795 5753487 : jCol = jCol - 1
1796 : case(16) !W
1797 11095532 : jCol = jCol - 1
1798 : case(32) !NW
1799 10002909 : iRow = iRow - 1
1800 10002909 : jCol = jCol - 1
1801 : case(64) !N
1802 7261091 : iRow = iRow - 1
1803 : case(128) !NE
1804 2701591 : iRow = iRow - 1
1805 2701591 : jCol = jCol + 1
1806 : case default !sink
1807 : ! do nothing
1808 : end select
1809 :
1810 4649970 : 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 35952 : subroutine cellLength(iDomain, fDir, iRow, jCol, iCoorSystem, length)
1841 :
1842 43740628 : use mo_common_types, only: Grid
1843 : use mo_common_variables, only : domainMeta, level0
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 35952 : real(dp) :: lat_1, long_1, lat_2, long_2
1863 :
1864 : type(Grid), pointer :: level0_iDomain => null()
1865 :
1866 :
1867 35952 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
1868 :
1869 : ! regular X-Y cordinate system
1870 35952 : IF(iCoorSystem .EQ. 0) THEN
1871 :
1872 54786 : select case (fDir)
1873 : case(1, 4, 16, 64) ! E, S, W, N
1874 18834 : length = 1.0_dp
1875 : case(2, 8, 32, 128) ! SE, SW, NW, NE
1876 35952 : length = SQRT2_dp
1877 : end select
1878 35952 : length = length * level0_iDomain%cellsize
1879 :
1880 : ! regular lat-lon cordinate system
1881 0 : ELSE IF(iCoorSystem .EQ. 1) THEN
1882 0 : iRow_to = iRow
1883 0 : jCol_to = jCol
1884 :
1885 : ! move in the direction of flow
1886 0 : 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 : 0.5_dp * level0_iDomain%cellsize
1891 : long_1 = level0_iDomain%xllcorner + real((iRow - 1), dp) * level0_iDomain%cellsize + &
1892 0 : 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 : 0.5_dp * level0_iDomain%cellsize
1896 : long_2 = level0_iDomain%xllcorner + real((iRow_to - 1), dp) * level0_iDomain%cellsize + &
1897 0 : 0.5_dp * level0_iDomain%cellsize
1898 : ! get distance between two points
1899 0 : call get_distance_two_lat_lon_points(lat_1, long_1, lat_2, long_2, length)
1900 :
1901 : END IF
1902 : !
1903 35952 : 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 0 : subroutine get_distance_two_lat_lon_points(lat1, long1, lat2, long2, distance_out)
1936 :
1937 35952 : 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 0 : real(dp) :: theta1
1948 :
1949 0 : real(dp) :: phi1
1950 :
1951 0 : real(dp) :: theta2
1952 :
1953 0 : real(dp) :: phi2
1954 :
1955 : real(dp) :: dtor
1956 :
1957 0 : real(dp) :: term1
1958 :
1959 0 : real(dp) :: term2
1960 :
1961 0 : real(dp) :: term3
1962 :
1963 0 : real(dp) :: temp
1964 :
1965 :
1966 0 : dtor = TWOPI_dp / 360.0_dp
1967 0 : theta1 = dtor * long1
1968 0 : phi1 = dtor * lat1
1969 0 : theta2 = dtor * long2
1970 0 : phi2 = dtor * lat2
1971 :
1972 0 : term1 = cos(phi1) * cos(theta1) * cos(phi2) * cos(theta2)
1973 0 : term2 = cos(phi1) * sin(theta1) * cos(phi2) * sin(theta2)
1974 0 : term3 = sin(phi1) * sin(phi2)
1975 0 : temp = term1 + term2 + term3
1976 0 : if(temp .GT. 1.0_dp) temp = 1.0_dp
1977 :
1978 0 : distance_out = RadiusEarth_dp * acos(temp);
1979 :
1980 0 : 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 24 : subroutine L11_flow_accumulation(iDomain)
2032 :
2033 : use mo_mrm_global_variables, only: &
2034 0 : L11_fDir, & ! IN: flow direction at L11 (standard notation)
2035 : L11_fAcc ! OUT: flow accumulation at L11 [km^2]
2036 : use mo_mrm_global_variables, only : level11
2037 : use mo_common_constants, only : nodata_i4, nodata_dp
2038 : use mo_append, only : append
2039 :
2040 : implicit none
2041 :
2042 : ! local variables
2043 : integer(i4), intent(in) :: iDomain
2044 24 : 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 24 : integer(i4), dimension(:,:), allocatable :: fDir11 ! L11_fDir array
2049 24 : logical, dimension(:,:), allocatable :: mask11 ! Domain mask
2050 :
2051 : ! initialize Domain info
2052 24 : nrows11 = level11(iDomain)%nrows
2053 24 : ncols11 = level11(iDomain)%ncols
2054 24 : s11 = level11(iDomain)%iStart
2055 24 : e11 = level11(iDomain)%iEnd
2056 1860 : mask11 = level11(iDomain)%mask
2057 :
2058 : ! allocate arrays
2059 96 : allocate(fDir11 (nrows11, ncols11))
2060 96 : allocate(fAcc11 (nrows11, ncols11))
2061 :
2062 : ! initialize
2063 1812 : fDir11(:,:) = nodata_i4
2064 :
2065 : ! get data
2066 24 : fDir11(:,:) = UNPACK( L11_fDir(s11:e11), mask11, nodata_i4 )
2067 :
2068 : ! initialize fAcc11 with cell area
2069 960 : fAcc11 = UNPACK( level11(iDomain)%cellarea * 1.e-6, mask11, nodata_dp )
2070 :
2071 : ! For each sink call "calculate_L11_flow_accumulation"
2072 250 : do jj=1, ncols11
2073 1812 : do ii=1, nrows11
2074 1788 : 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 24 : ncol = ncols11)
2081 : end if
2082 : end do
2083 : end do
2084 :
2085 : ! Append
2086 24 : call append( L11_fAcc, pack(fAcc11(:,:),mask11))
2087 :
2088 : ! free space
2089 24 : deallocate(fDir11, fAcc11, mask11)
2090 :
2091 : contains
2092 :
2093 936 : 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 936 : if (jj+1 .le. ncol) then
2113 912 : if (fDir(ii,jj+1) .eq. 16_i4) then
2114 402 : call calculate_L11_flow_accumulation(fDir, fAcc, ii, jj+1, nrow, ncol)
2115 402 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii,jj+1)
2116 : end if
2117 : end if
2118 :
2119 936 : if ((ii+1 .le. nrow) .and. (jj+1 .le. ncol)) then
2120 798 : if (fDir(ii+1,jj+1) .eq. 32_i4) then
2121 0 : call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj+1, nrow, ncol)
2122 0 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj+1)
2123 : end if
2124 : end if
2125 :
2126 936 : if (ii+1 .le. nrow) then
2127 822 : if (fDir(ii+1,jj) .eq. 64_i4) then
2128 320 : call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj, nrow, ncol)
2129 320 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj)
2130 : end if
2131 : end if
2132 :
2133 936 : if ((ii+1 .le. nrow) .and. (jj-1 .ge. 1)) then
2134 770 : if (fDir(ii+1,jj-1) .eq. 128_i4) then
2135 0 : call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj-1, nrow, ncol)
2136 0 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj-1)
2137 : end if
2138 : end if
2139 :
2140 936 : if (jj-1 .ge. 1) then
2141 884 : if (fDir(ii,jj-1) .eq. 1_i4) then
2142 50 : call calculate_L11_flow_accumulation(fDir, fAcc, ii, jj-1, nrow, ncol)
2143 50 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii,jj-1)
2144 : end if
2145 : end if
2146 :
2147 936 : if ((ii-1 .ge. 1) .and. (jj-1 .ge. 1)) then
2148 842 : if (fDir(ii-1,jj-1) .eq. 2_i4) then
2149 0 : call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj-1, nrow, ncol)
2150 0 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii-1,jj-1)
2151 : end if
2152 : end if
2153 :
2154 936 : if (ii-1 .ge. 1) then
2155 892 : if (fDir(ii-1,jj) .eq. 4_i4) then
2156 138 : call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj, nrow, ncol)
2157 138 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii-1,jj)
2158 : end if
2159 : end if
2160 :
2161 936 : if ((ii-1 .ge. 1) .and. (jj+1 .le. ncol)) then
2162 868 : if (fDir11(ii-1,jj+1) .eq. 8_i4) then
2163 2 : call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj+1, nrow, ncol)
2164 2 : 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 24 : 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 17 : subroutine L11_calc_celerity(iDomain, param)
2222 : use mo_common_constants, only: nodata_i4, nodata_dp
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 17 : logical, dimension(:,:), allocatable :: mask0
2259 17 : integer(i4), dimension(:,:), allocatable :: iD0
2260 17 : integer(i4), dimension(:,:), allocatable :: fDir0
2261 17 : integer(i4), dimension(:,:), allocatable :: fAcc0
2262 17 : real(dp), dimension(:,:), allocatable :: slope0
2263 17 : real(dp), dimension(:), allocatable :: slope_tmp
2264 17 : real(dp), dimension(:,:), allocatable :: cellarea0
2265 17 : integer(i4), dimension(:), allocatable :: netPerm ! routing order (permutation)
2266 17 : integer(i4), dimension(:), allocatable :: nLinkFromRow
2267 17 : integer(i4), dimension(:), allocatable :: nLinkFromCol
2268 17 : integer(i4), dimension(:), allocatable :: nLinkToRow
2269 17 : integer(i4), dimension(:), allocatable :: nLinkToCol
2270 : integer(i4) :: ii, rr, ns
2271 : integer(i4) :: frow, fcol
2272 : integer(i4) :: fId, tId
2273 17 : real(dp), dimension(:), allocatable :: stack, append_chunk ! Stacks celerity along the L0 river-path
2274 17 : integer(i4), dimension(:), allocatable :: dummy_1d
2275 :
2276 17 : real(dp) :: L0_link_slope
2277 17 : real(dp), dimension(:), allocatable :: celerity11
2278 17 : real(dp), dimension(:,:), allocatable :: celerity0
2279 :
2280 17 : integer(i4), dimension(:,:), allocatable :: nodata_i4_tmp
2281 17 : real(dp), dimension(:,:), allocatable :: nodata_dp_tmp
2282 17 : logical, dimension(:), allocatable :: slopemask0
2283 :
2284 : type(Grid), pointer :: level0_iDomain
2285 :
2286 : ! level-0 information
2287 17 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
2288 17 : nrows0 = level0_iDomain%nrows
2289 17 : ncols0 = level0_iDomain%ncols
2290 17 : nCells0 = level0_iDomain%ncells
2291 17 : iStart0 = level0_iDomain%iStart
2292 17 : iEnd0 = level0_iDomain%iEnd
2293 2122467 : mask0 = level0_iDomain%mask
2294 :
2295 : ! level-11 information
2296 17 : iStart11 = level11(iDomain)%iStart
2297 17 : iEnd11 = level11(iDomain)%iEnd
2298 17 : nrows11 = level11(iDomain)%nrows
2299 17 : ncols11 = level11(iDomain)%ncols
2300 17 : nNodes = level11(iDomain)%ncells
2301 :
2302 17 : nLinks = nNodes - L11_nOutlets(iDomain)
2303 :
2304 : ! allocate
2305 68 : allocate ( iD0 ( nrows0, ncols0 ) )
2306 68 : allocate ( slope0 ( nrows0, ncols0 ) )
2307 51 : allocate ( fDir0 ( nrows0, ncols0 ) )
2308 51 : allocate ( fAcc0 ( nrows0, ncols0 ) )
2309 51 : allocate ( cellarea0 ( nrows0, ncols0 ) )
2310 51 : allocate ( celerity0 ( nrows0, ncols0 ) )
2311 51 : 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 17 : allocate ( stack ( 1 ) )
2316 17 : allocate ( append_chunk( 1 ) )
2317 17 : allocate ( dummy_1d ( 2 ))
2318 51 : allocate ( netPerm ( nNodes ) )
2319 34 : allocate ( nLinkFromRow( nNodes ) )
2320 34 : allocate ( nLinkFromCol( nNodes ) )
2321 34 : allocate ( nLinkToRow ( nNodes ) )
2322 34 : allocate ( nLinkToCol ( nNodes ) )
2323 51 : allocate ( celerity11 ( nNodes ) )
2324 34 : allocate ( slope_tmp ( nNodes ) )
2325 :
2326 51 : allocate (nodata_i4_tmp( nrows0, ncols0 ) )
2327 51 : allocate (nodata_dp_tmp( nrows0, ncols0 ) )
2328 :
2329 : ! initialize
2330 2122433 : iD0(:,:) = nodata_i4
2331 2122433 : fDir0(:,:) = nodata_i4
2332 2122433 : fAcc0(:,:) = nodata_i4
2333 2122433 : cellarea0(:,:) = nodata_dp
2334 2122433 : slope0(:,:) = nodata_dp
2335 :
2336 34 : stack(:) = nodata_dp
2337 34 : append_chunk(:) = nodata_dp
2338 595 : netPerm(:) = nodata_i4
2339 595 : nLinkFromRow(:) = nodata_i4
2340 595 : nLinkFromCol(:) = nodata_i4
2341 595 : nLinkToRow(:) = nodata_i4
2342 595 : nLinkToCol(:) = nodata_i4
2343 595 : celerity11(:) = nodata_dp
2344 2122433 : celerity0(:,:) = nodata_dp
2345 791282 : slopemask0(:) = .False.
2346 :
2347 2122433 : nodata_i4_tmp(:,:) = nodata_i4
2348 2122433 : nodata_dp_tmp(:,:) = nodata_dp
2349 :
2350 : ! for a single node model run
2351 :
2352 17 : if(nNodes .GT. 1) then
2353 : ! get L0 fields
2354 17 : iD0(:,:) = UNPACK(level0_iDomain%Id(1:nCells0), mask0, nodata_i4_tmp)
2355 17 : fDir0(:,:) = UNPACK(L0_fDir(iStart0:iEnd0), mask0, nodata_i4_tmp)
2356 17 : fAcc0(:,:) = UNPACK(L0_fAcc(iStart0:iEnd0), mask0, nodata_i4_tmp)
2357 17 : cellarea0(:,:) = UNPACK(level0_iDomain%cellarea(1:nCells0), mask0, nodata_dp_tmp)
2358 :
2359 : ! smoothing river slope
2360 791299 : slope_tmp = L0_slope(iStart0:iEnd0)
2361 791282 : where ( slope_tmp .lt. 0.1_dp ) slope_tmp = 0.1_dp
2362 :
2363 791282 : slopemask0(:) = (L0_streamNet(iStart0:iEnd0) .ne. nodata_i4)
2364 :
2365 : ! smooth river cells if there is more than one cell
2366 791282 : if( count(slopemask0) .GT. 1) then
2367 1582564 : slope_tmp = mad(arr = slope_tmp, z = 2.25_dp, mask = slopemask0, tout='u', mval=0.1_dp)
2368 : end if
2369 17 : slope0(:,:) = UNPACK(slope_tmp, mask0, nodata_dp_tmp )
2370 :
2371 : ! get network vectors of L11
2372 595 : netPerm(:) = L11_netPerm ( iStart11 : iEnd11 )
2373 595 : nLinkFromRow(:) = L11_fRow ( iStart11 : iEnd11 )
2374 595 : nLinkFromCol(:) = L11_fCol ( iStart11 : iEnd11 )
2375 595 : nLinkToRow(:) = L11_tRow ( iStart11 : iEnd11 )
2376 595 : nLinkToCol(:) = L11_tCol ( iStart11 : iEnd11 )
2377 :
2378 578 : do rr = 1, nLinks
2379 :
2380 561 : ii = netPerm(rr)
2381 561 : frow = nLinkFromRow(ii)
2382 561 : fcol = nLinkFromCol(ii)
2383 :
2384 : ! Init
2385 1122 : stack(:) = 0_dp
2386 561 : ns = 1
2387 :
2388 561 : fId = iD0( frow, fcol )
2389 561 : tId = iD0( nLinkToRow(ii) , nLinkToCol(ii) )
2390 24514 : do
2391 25075 : L0_link_slope = slope0(frow, fcol) / 100._dp
2392 :
2393 : ! celerity parametrization
2394 25075 : stack(ns) = param(1) * sqrt(L0_link_slope)
2395 :
2396 25075 : celerity0(frow, fcol) = stack(ns)
2397 25075 : ns = ns + 1
2398 25075 : fId = iD0(frow, fcol)
2399 25075 : if( .NOT. (fID == tID)) then
2400 24514 : call append(stack, append_chunk)
2401 : else
2402 : exit
2403 : end if
2404 : ! move downstream
2405 24514 : call moveDownOneCell( fDir0(frow,fcol), frow, fcol )
2406 : end do
2407 :
2408 25636 : celerity11(ii) = size(stack) / sum(1/stack(:))
2409 561 : deallocate(stack)
2410 578 : 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 0 : celerity11(:) = 1._dp
2419 :
2420 : end if
2421 :
2422 : ! Write celerity
2423 595 : L11_celerity(iStart11:iEnd11) = celerity11(:)
2424 17 : L0_celerity(iStart0:iEnd0) = PACK(celerity0(:,:), mask0)
2425 :
2426 : ! free space
2427 0 : deallocate (&
2428 0 : mask0, iD0, slope_tmp, slopemask0, &
2429 0 : slope0, fDir0, cellarea0, &
2430 17 : stack, netPerm, nLinkFromRow, nLinkFromCol, nLinkToRow, nLinkToCol)
2431 :
2432 17 : end subroutine L11_calc_celerity
2433 :
2434 : end module mo_mrm_net_startup
|
