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
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) nlinksink(ii) = .True.
836 : end do
837 960 : where(nlinksink) nLinkLabel = 2 ! 'Sink'
838 :
839 : ! keep routing order
840 936 : do ii = 1, nLinks
841 936 : 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 24 : call append(L11_rOrder, nLinkROrder(:))
852 24 : call append(L11_label, nLinkLabel(:))
853 24 : call append(L11_sink, nLinkSink(:))
854 24 : call append(L11_netPerm, netPerm(:))
855 :
856 : ! free space
857 24 : deallocate (nLinkFromN, nLinkToN, nLinkROrder, nLinkLabel, nLinkSink, netPerm)
858 :
859 24 : 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 24 : subroutine L11_link_location(iDomain)
888 :
889 24 : use mo_append, only : append
890 : use mo_common_constants, only : nodata_i4
891 : use mo_common_types, only: Grid
892 : use mo_common_variables, only : domainMeta, level0
893 : use mo_mrm_global_variables, only : L0_draSC, L0_fDir, L11_colOut, L11_fCol, L11_fRow, L11_fromN, &
894 : L11_nOutlets, L11_netPerm, L11_rowOut, L11_tCol, L11_tRow, domain_mrm, level11
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 24 : integer(i4), dimension(:), allocatable :: rowOut
906 :
907 : ! easting cell loc. of the Outlet
908 24 : integer(i4), dimension(:), allocatable :: colOut
909 :
910 24 : integer(i4), dimension(:), allocatable :: nLinkFromN
911 :
912 24 : integer(i4), dimension(:), allocatable :: netPerm
913 :
914 24 : integer(i4), dimension(:), allocatable :: nLinkFromRow
915 :
916 24 : integer(i4), dimension(:), allocatable :: nLinkFromCol
917 :
918 24 : integer(i4), dimension(:), allocatable :: nLinkToRow
919 :
920 24 : integer(i4), dimension(:), allocatable :: nLinkToCol
921 :
922 24 : integer(i4), dimension(:, :), allocatable :: fDir0
923 :
924 24 : 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 24 : 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 24 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
943 24 : s0 = level0_iDomain%iStart
944 24 : e0 = level0_iDomain%iEnd
945 :
946 24 : nOutlets = L11_nOutlets(iDomain)
947 :
948 24 : 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 72 : allocate (rowOut (level11(iDomain)%nCells))
954 48 : allocate (colOut (level11(iDomain)%nCells))
955 48 : allocate (nLinkFromN (level11(iDomain)%nCells)) ! all network vectors valid from (1 : nLinks)
956 48 : allocate (netPerm (level11(iDomain)%nCells))
957 48 : allocate (nLinkFromRow (level11(iDomain)%nCells))
958 48 : allocate (nLinkFromCol (level11(iDomain)%nCells))
959 48 : allocate (nLinkToRow (level11(iDomain)%nCells))
960 48 : allocate (nLinkToCol (level11(iDomain)%nCells))
961 96 : allocate (fDir0 (level0_iDomain%nrows, level0_iDomain%ncols))
962 72 : allocate (draSC0 (level0_iDomain%nrows, level0_iDomain%ncols))
963 :
964 : ! initialize
965 984 : rowOut = nodata_i4
966 984 : colOut = nodata_i4
967 984 : nLinkFromN = nodata_i4
968 984 : netPerm = nodata_i4
969 984 : nLinkFromRow = nodata_i4
970 984 : nLinkFromCol = nodata_i4
971 984 : nLinkToRow = nodata_i4
972 984 : nLinkToCol = nodata_i4
973 2862384 : fDir0 = nodata_i4
974 2862384 : draSC0 = nodata_i4
975 :
976 : ! for a single node model run
977 24 : if(level11(iDomain)%nCells .GT. 1) then
978 : ! get fDir at L0
979 24 : fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), level0_iDomain%mask, nodata_i4)
980 24 : draSC0(:, :) = UNPACK(L0_draSC (s0 : e0), level0_iDomain%mask, nodata_i4)
981 :
982 : ! get network vectors of L11
983 960 : nLinkFromN(:) = L11_fromN (level11(iDomain)%iStart : level11(iDomain)%iEnd)
984 960 : netPerm(:) = L11_netPerm (level11(iDomain)%iStart : level11(iDomain)%iEnd)
985 960 : rowOut(:) = L11_rowOut (level11(iDomain)%iStart : level11(iDomain)%iEnd)
986 960 : colOut(:) = L11_colOut (level11(iDomain)%iStart : level11(iDomain)%iEnd)
987 :
988 : ! finding main outlet (row, col) in L0
989 72 : allocate(oLoc(Noutlets, 2))
990 48 : oLoc(:, 1) = domain_mrm(iDomain)%L0_rowOutlet(: Noutlets)
991 48 : oLoc(:, 2) = domain_mrm(iDomain)%L0_colOutlet(: Noutlets)
992 :
993 : ! Location of the stream-joint cells (row, col)
994 936 : do rr = 1, nLinks
995 :
996 912 : ii = netPerm(rr)
997 912 : iNode = nLinkFromN(ii)
998 912 : iRow = rowOut(iNode)
999 912 : jCol = colOut(iNode)
1000 912 : call moveDownOneCell(fDir0(iRow, jcol), iRow, jcol)
1001 : ! set "from" cell
1002 912 : nLinkFromRow(ii) = iRow
1003 912 : nLinkFromCol(ii) = jCol
1004 :
1005 : ! check whether this location is an outlet
1006 912 : is_outlet = .False.
1007 1824 : do kk = 1, Noutlets
1008 1824 : if (iRow .eq. oLoc(kk, 1) .and. jCol .eq. oLoc(kk, 2)) is_outlet = .True.
1009 : end do
1010 :
1011 1848 : if (is_outlet) then
1012 :
1013 0 : nLinkToRow(ii) = iRow
1014 0 : nLinkToCol(ii) = jCol
1015 :
1016 : else
1017 :
1018 35952 : do while (.not. (draSC0(iRow, jCol) > 0))
1019 35040 : prevRow = iRow
1020 35040 : prevCol = jCol
1021 : call moveDownOneCell(fDir0(iRow, jcol), iRow, jCol)
1022 : ! check whether this location is an outlet and exit
1023 70056 : do kk = 1, Noutlets
1024 70056 : if (iRow .eq. oLoc(kk, 1) .and. jCol .eq. oLoc(kk, 2)) exit
1025 : end do
1026 35952 : 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 0 : num2str(jCol))
1030 : end if
1031 : end do
1032 : ! set "to" cell (when an outlet is reached)
1033 912 : nLinkToRow(ii) = iRow
1034 912 : 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 24 : call append(L11_fRow, nLinkFromRow(:))
1047 24 : call append(L11_fCol, nLinkFromCol(:))
1048 24 : call append(L11_tRow, nLinkToRow(:))
1049 24 : call append(L11_tCol, nLinkToCol(:))
1050 :
1051 : ! free space
1052 : deallocate (rowOut, colOut, nLinkFromN, netPerm, nLinkFromRow, &
1053 24 : nLinkFromCol, nLinkToRow, nLinkToCol, fDir0, draSC0)
1054 :
1055 24 : 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 :
1088 24 : subroutine L11_set_drain_outlet_gauges(iDomain)
1089 :
1090 24 : use mo_append, only : append
1091 : use mo_common_constants, only : nodata_i4
1092 : use mo_common_types, only: Grid
1093 : use mo_common_variables, only : domainMeta, level0
1094 : use mo_mrm_global_variables, only : L0_InflowgaugeLoc, L0_draCell, L0_draSC, L0_fDir, L0_gaugeLoc, domain_mrm, &
1095 : l0_l11_remap
1096 :
1097 : implicit none
1098 :
1099 : ! Domain Id
1100 : integer(i4), intent(in) :: iDomain
1101 :
1102 24 : integer(i4), dimension(:, :), allocatable :: draSC0
1103 :
1104 24 : integer(i4), dimension(:, :), allocatable :: fDir0
1105 :
1106 24 : integer(i4), dimension(:, :), allocatable :: gaugeLoc0
1107 :
1108 24 : integer(i4), dimension(:, :), allocatable :: InflowGaugeLoc0
1109 :
1110 24 : 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 24 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
1122 24 : s0 = level0_iDomain%iStart
1123 24 : e0 = level0_iDomain%iEnd
1124 :
1125 :
1126 : ! allocate
1127 96 : allocate (draSC0 (level0_iDomain%nrows, level0_iDomain%ncols))
1128 72 : allocate (fDir0 (level0_iDomain%nrows, level0_iDomain%ncols))
1129 72 : allocate (gaugeLoc0 (level0_iDomain%nrows, level0_iDomain%ncols))
1130 72 : allocate (InflowGaugeLoc0 (level0_iDomain%nrows, level0_iDomain%ncols))
1131 72 : allocate (draCell0 (level0_iDomain%nrows, level0_iDomain%ncols))
1132 :
1133 : ! initialize
1134 2862360 : draSC0(:, :) = nodata_i4
1135 2862360 : fDir0(:, :) = nodata_i4
1136 2862360 : gaugeLoc0(:, :) = nodata_i4
1137 2862360 : InflowGaugeLoc0(:, :) = nodata_i4
1138 2862360 : draCell0(:, :) = nodata_i4
1139 :
1140 0 : draSC0(:, :) = UNPACK(L0_draSC (s0 : e0), &
1141 24 : level0_iDomain%mask, nodata_i4)
1142 0 : fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), &
1143 24 : level0_iDomain%mask, nodata_i4)
1144 0 : gaugeLoc0(:, :) = UNPACK(L0_gaugeLoc (s0 : e0), &
1145 24 : level0_iDomain%mask, nodata_i4)
1146 0 : InflowGaugeLoc0(:, :) = UNPACK(L0_InflowgaugeLoc (s0 : e0), &
1147 24 : level0_iDomain%mask, nodata_i4)
1148 :
1149 1066064 : do kk = 1, level0_iDomain%nCells
1150 1066040 : ii = level0_iDomain%CellCoor(kk, 1)
1151 1066040 : jj = level0_iDomain%CellCoor(kk, 2)
1152 1066040 : iSc = draSC0(ii, jj)
1153 : ! find drainage path
1154 1066040 : iRow = ii
1155 1066040 : jCol = jj
1156 43644186 : do while (.NOT. iSC > 0)
1157 : ! move downstream
1158 42578146 : call moveDownOneCell(fDir0(iRow, jCol), iRow, jCol)
1159 42578146 : iSC = draSC0(iRow, jCol)
1160 : end do
1161 1066040 : 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 1066040 : if (gaugeLoc0(ii, jj) .NE. nodata_i4) then
1166 : ! evaluation gauges
1167 76 : do ll = 1, domain_mrm(iDomain)%nGauges
1168 : ! search for gaugeID in L0 grid and save ID on L11
1169 76 : if (domain_mrm(iDomain)%gaugeIdList(ll) .EQ. gaugeLoc0(ii, jj)) then
1170 21 : 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 1066064 : if (InflowGaugeLoc0(ii, jj) .NE. nodata_i4) then
1176 : ! inflow gauges
1177 39 : do ll = 1, domain_mrm(iDomain)%nInflowGauges
1178 : ! search for gaugeID in L0 grid and save ID on L11
1179 4 : if (domain_mrm(iDomain)%InflowGaugeIdList(ll) .EQ. InflowGaugeLoc0(ii, jj)) &
1180 37 : 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 24 : if (iDomain .eq. 1) then
1192 12 : call append(L0_draCell, PACK(draCell0(:, :), level0_iDomain%mask))
1193 12 : else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
1194 9 : call append(L0_draCell, PACK(draCell0(:, :), level0_iDomain%mask))
1195 : end if
1196 :
1197 : ! free space
1198 24 : deallocate (draSC0, fDir0, gaugeLoc0, draCell0)
1199 :
1200 24 : 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 24 : subroutine L11_stream_features(iDomain)
1234 :
1235 24 : use mo_append, only : append
1236 : use mo_common_constants, only : nodata_dp, nodata_i4
1237 : use mo_common_types, only: Grid
1238 : use mo_common_variables, only : domainMeta, L0_elev, iFlag_cordinate_sys, level0, processMatrix
1239 : use mo_mrm_global_variables, only : L0_fDir, &
1240 : L0_floodPlain, L0_streamNet, L11_aFloodPlain, L11_fCol, L11_fRow, L11_length, &
1241 : L11_nOutlets, L11_netPerm, L11_slope, L11_tCol, L11_tRow, level11
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 24 : integer(i4), dimension(:, :), allocatable :: iD0
1252 :
1253 24 : integer(i4), dimension(:, :), allocatable :: fDir0
1254 :
1255 24 : real(dp), dimension(:, :), allocatable :: elev0
1256 :
1257 24 : real(dp), dimension(:, :), allocatable :: cellarea0
1258 :
1259 24 : integer(i4), dimension(:, :), allocatable :: streamNet0
1260 :
1261 24 : integer(i4), dimension(:, :), allocatable :: floodPlain0
1262 :
1263 : ! routing order (permutation)
1264 24 : integer(i4), dimension(:), allocatable :: netPerm
1265 :
1266 24 : integer(i4), dimension(:), allocatable :: nLinkFromRow
1267 :
1268 24 : integer(i4), dimension(:), allocatable :: nLinkFromCol
1269 :
1270 24 : integer(i4), dimension(:), allocatable :: nLinkToRow
1271 :
1272 24 : integer(i4), dimension(:), allocatable :: nLinkToCol
1273 :
1274 24 : real(dp), dimension(:), allocatable :: nLinkLength
1275 :
1276 24 : real(dp), dimension(:), allocatable :: nLinkAFloodPlain
1277 :
1278 24 : 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 24 : integer(i4), dimension(:, :), allocatable :: stack, append_chunk
1287 :
1288 24 : integer(i4), dimension(:), allocatable :: dummy_1d
1289 :
1290 24 : real(dp) :: length
1291 :
1292 24 : integer(i4), dimension(:, :), allocatable :: nodata_i4_tmp
1293 :
1294 24 : real(dp), dimension(:, :), allocatable :: nodata_dp_tmp
1295 :
1296 : type(Grid), pointer :: level0_iDomain => null()
1297 :
1298 :
1299 24 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
1300 24 : s0 = level0_iDomain%iStart
1301 24 : e0 = level0_iDomain%iEnd
1302 24 : nLinks = level11(iDomain)%nCells - L11_nOutlets(iDomain)
1303 :
1304 :
1305 : ! allocate
1306 96 : allocate (iD0 (level0_iDomain%nrows, level0_iDomain%ncols))
1307 96 : allocate (elev0 (level0_iDomain%nrows, level0_iDomain%ncols))
1308 96 : allocate (fDir0 (level0_iDomain%nrows, level0_iDomain%ncols))
1309 96 : allocate (cellarea0 (level0_iDomain%nrows, level0_iDomain%ncols))
1310 96 : allocate (streamNet0 (level0_iDomain%nrows, level0_iDomain%ncols))
1311 72 : 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 96 : allocate (stack (level11(iDomain)%nCells, 2)) ! >> stack(nNodes, 2)
1316 24 : allocate (dummy_1d (2))
1317 24 : allocate (append_chunk (8, 2))
1318 72 : allocate (netPerm (level11(iDomain)%nCells))
1319 48 : allocate (nLinkFromRow (level11(iDomain)%nCells))
1320 48 : allocate (nLinkFromCol (level11(iDomain)%nCells))
1321 48 : allocate (nLinkToRow (level11(iDomain)%nCells))
1322 48 : allocate (nLinkToCol (level11(iDomain)%nCells))
1323 72 : allocate (nLinkLength (level11(iDomain)%nCells))
1324 48 : allocate (nLinkAFloodPlain (level11(iDomain)%nCells))
1325 48 : allocate (nLinkSlope (level11(iDomain)%nCells))
1326 :
1327 96 : allocate (nodata_i4_tmp (level0_iDomain%nrows, level0_iDomain%ncols))
1328 96 : allocate (nodata_dp_tmp (level0_iDomain%nrows, level0_iDomain%ncols))
1329 :
1330 : ! initialize
1331 2862360 : iD0(:, :) = nodata_i4
1332 2862360 : elev0(:, :) = nodata_dp
1333 2862360 : fDir0(:, :) = nodata_i4
1334 2862360 : cellarea0(:, :) = nodata_dp
1335 2862360 : streamNet0(:, :) = nodata_i4
1336 2862360 : floodPlain0(:, :) = nodata_i4
1337 :
1338 1944 : stack(:, :) = nodata_i4
1339 456 : append_chunk(:, :) = nodata_i4
1340 960 : netPerm(:) = nodata_i4
1341 960 : nLinkFromRow(:) = nodata_i4
1342 960 : nLinkFromCol(:) = nodata_i4
1343 960 : nLinkToRow(:) = nodata_i4
1344 960 : nLinkToCol(:) = nodata_i4
1345 960 : nLinkLength(:) = nodata_dp
1346 960 : nLinkAFloodPlain(:) = nodata_dp
1347 960 : nLinkSlope(:) = nodata_dp
1348 :
1349 2862360 : nodata_i4_tmp(:, :) = nodata_i4
1350 2862360 : nodata_dp_tmp(:, :) = nodata_dp
1351 :
1352 : ! for a single node model run
1353 24 : if(level11(iDomain)%nCells .GT. 1) then
1354 : ! get L0 fields
1355 24 : iD0(:, :) = UNPACK(level0_iDomain%Id, level0_iDomain%mask, nodata_i4_tmp)
1356 0 : elev0(:, :) = UNPACK(L0_elev (s0 : e0), &
1357 24 : level0_iDomain%mask, nodata_dp_tmp)
1358 0 : fDir0(:, :) = UNPACK(L0_fDir (s0 : e0), &
1359 24 : level0_iDomain%mask, nodata_i4_tmp)
1360 24 : cellarea0(:, :) = UNPACK(level0_iDomain%CellArea, level0_iDomain%mask, nodata_dp_tmp)
1361 :
1362 : ! get network vectors of L11
1363 960 : netPerm(:) = L11_netPerm (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1364 960 : nLinkFromRow(:) = L11_fRow (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1365 960 : nLinkFromCol(:) = L11_fCol (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1366 960 : nLinkToRow(:) = L11_tRow (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1367 960 : nLinkToCol(:) = L11_tCol (level11(iDomain)%iStart : level11(iDomain)%iEnd)
1368 :
1369 : ! Flood plains: stream network delineation
1370 2862360 : streamNet0(:, :) = nodata_i4
1371 2862360 : floodPlain0(:, :) = nodata_i4
1372 :
1373 936 : do rr = 1, nLinks
1374 :
1375 912 : ii = netPerm(rr)
1376 912 : frow = nLinkFromRow(ii)
1377 912 : fcol = nLinkFromCol(ii)
1378 :
1379 : ! Init
1380 912 : streamNet0(frow, fcol) = ii
1381 912 : floodPlain0(frow, fcol) = ii
1382 144024 : stack = 0
1383 18240 : append_chunk = 0
1384 912 : ns = 1
1385 912 : stack(ns, 1) = frow
1386 912 : stack(ns, 2) = fcol
1387 :
1388 912 : call cellLength(iDomain, fDir0(frow, fcol), fRow, fCol, iFlag_cordinate_sys, nLinkLength(ii))
1389 912 : nLinkSlope(ii) = elev0(frow, fcol)
1390 :
1391 912 : fId = iD0(frow, fcol)
1392 912 : tId = iD0(nLinkToRow(ii), nLinkToCol(ii))
1393 :
1394 35952 : do while (.NOT. (fId == tId))
1395 : ! Search flood plain from point(frow,fcol) upwards, keep co-ordinates in STACK
1396 6659002 : do while (ns > 0)
1397 6623962 : if (ns + 8 .gt. size(stack, 1)) then
1398 118 : call append(stack, append_chunk)
1399 : end if
1400 6623962 : call moveUp(elev0, fDir0, frow, fcol, stack, ns)
1401 6623962 : stack(1, 1) = 0
1402 6623962 : stack(1, 2) = 0
1403 : ! stack = cshift(stack, SHIFT = 1, DIM = 1)
1404 : ! substitute cshift <<<
1405 26495848 : dummy_1d = stack(1, :)
1406 2009564780 : stack(: size(stack, dim = 1) - 1, :) = stack(2 :, :)
1407 19871886 : stack(size(stack, dim = 1), :) = dummy_1d
1408 : ! substitute cshift >>>
1409 6623962 : if (stack(1, 1) > 0 .and. stack(1, 2) > 0) floodPlain0(stack(1, 1), stack(1, 2)) = ii
1410 1018065354 : ns = count(stack > 0) / 2
1411 : end do
1412 :
1413 : ! move downstream
1414 35040 : call moveDownOneCell(fDir0(frow, fcol), frow, fcol)
1415 35040 : streamNet0(frow, fcol) = ii
1416 35040 : floodPlain0(frow, fcol) = ii
1417 35040 : fId = iD0(frow, fcol)
1418 5303016 : stack = 0
1419 35040 : ns = 1
1420 35040 : stack(ns, 1) = frow
1421 35040 : stack(ns, 2) = fcol
1422 35040 : call cellLength(iDomain, fDir0(fRow, fCol), fRow, fCol, iFlag_cordinate_sys, length)
1423 35040 : nLinkLength(ii) = nLinkLength(ii) + length
1424 :
1425 : end do
1426 :
1427 : ! stream bed slope
1428 912 : nLinkSlope(ii) = (nLinkSlope(ii) - elev0(frow, fcol)) / nLinkLength(ii)
1429 :
1430 912 : if (nLinkSlope(ii) < 0.0001_dp) nLinkSlope(ii) = 0.0001_dp
1431 :
1432 : ! calculate area of floodplains (avoid overwriting)
1433 111450936 : 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 24 : if ((processMatrix(8, 1) .eq. 2) .or. (processMatrix(8, 1) .eq. 3)) then
1443 280 : if (count(nLinkLength(:) .ge. 0._dp) .gt. 2) then
1444 280 : length = percentile(pack(nLinkLength(:), nLinkLength(:) .ge. 0._dp), 40._dp)
1445 280 : 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 24 : if (iDomain .eq. 1) then
1457 12 : call append(L0_streamNet, PACK (streamNet0(:, :), level0_iDomain%mask))
1458 12 : call append(L0_floodPlain, PACK (floodPlain0(:, :), level0_iDomain%mask))
1459 12 : else if (domainMeta%L0DataFrom(iDomain) == iDomain) then
1460 9 : call append(L0_streamNet, PACK (streamNet0(:, :), level0_iDomain%mask))
1461 9 : call append(L0_floodPlain, PACK (floodPlain0(:, :), level0_iDomain%mask))
1462 : end if
1463 :
1464 :
1465 : ! L11 network data sets
1466 24 : call append(L11_length, nLinkLength(:))
1467 24 : call append(L11_aFloodPlain, nLinkAFloodPlain(:))
1468 24 : call append(L11_slope, nLinkSlope(:))
1469 :
1470 : ! free space
1471 0 : deallocate (&
1472 0 : iD0, elev0, fDir0, streamNet0, floodPlain0, &
1473 0 : cellarea0, stack, netPerm, nLinkFromRow, nLinkFromCol, nLinkToRow, nLinkToCol, &
1474 24 : nLinkLength, nLinkAFloodPlain, nLinkSlope, dummy_1d)
1475 24 : deallocate(nodata_i4_tmp, nodata_dp_tmp)
1476 :
1477 24 : 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 13 : subroutine L11_fraction_sealed_floodplain(LCClassImp, do_init)
1511 :
1512 24 : use mo_append, only : append
1513 : use mo_common_constants, only : nodata_dp
1514 : use mo_common_types, only: Grid
1515 : use mo_common_variables, only : domainMeta, L0_LCover, level0, domainMeta, nLCoverScene
1516 : use mo_mrm_global_variables, only : L0_floodPlain, L11_aFloodPlain, &
1517 : L11_nLinkFracFPimp, L11_nOutlets, level11
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 13 : 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 38 : do iDomain = 1, domainMeta%nDomains
1539 100 : allocate(temp_array(level11(iDomain)%nCells, nLCoverScene))
1540 2040 : temp_array = nodata_dp
1541 25 : if (do_init) then
1542 16 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
1543 :
1544 16 : s0 = level0_iDomain%iStart
1545 16 : e0 = level0_iDomain%iEnd
1546 16 : nLinks = level11(iDomain)%nCells + 1 - L11_nOutlets(iDomain)
1547 16 : nLinkAFloodPlain => L11_aFloodPlain(level11(iDomain)%iStart : level11(iDomain)%iEnd)
1548 :
1549 48 : do iiLC = 1, nLCoverScene
1550 : ! for a single node model run
1551 48 : if(nLinks .GT. 0) then
1552 1360 : do ii = 1, nLinks
1553 2656 : temp_array(ii, iiLC) = sum(level0_iDomain%CellArea(:), &
1554 0 : mask = (L0_floodPlain(s0 : e0) == ii .and. L0_LCover(s0 : e0, iiLC) == LCClassImp)) &
1555 59876256 : / nLinkAFloodPlain(ii)
1556 : end do
1557 : end if
1558 : end do
1559 : end if
1560 25 : call append(L11_nLinkFracFPimp, temp_array(:,:))
1561 38 : deallocate(temp_array)
1562 : end do
1563 :
1564 13 : 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 6623962 : subroutine moveUp(elev0, fDir0, fi, fj, ss, nn)
1596 :
1597 13 : 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 6623962 : ii = ss(1, 1)
1619 6623962 : jj = ss(1, 2)
1620 6623962 : ip = ii + 1
1621 6623962 : im = ii - 1
1622 6623962 : jp = jj + 1
1623 6623962 : jm = jj - 1
1624 :
1625 6623962 : nrows = size(fDir0, 1)
1626 6623962 : ncols = size(fDir0, 2)
1627 :
1628 : !E
1629 6623962 : if (jp <= ncols) then
1630 13247924 : if ((fdir0(ii, jp) == 16) .and. &
1631 26495848 : (le((elev0(ii, jp) - elev0(fi, fj)), deltaH)) .and. &
1632 6623962 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1633 : ) then
1634 1415524 : nn = nn + 1
1635 1415524 : ss(nn, 1) = ii
1636 1415524 : ss(nn, 2) = jp
1637 : !print *, i,jp
1638 : end if
1639 : end if
1640 :
1641 : !SE
1642 6623962 : if ((ip <= nrows) .and. &
1643 : (jp <= ncols)) then
1644 0 : if ((fdir0(ip, jp) == 32) .and. &
1645 13247924 : (le((elev0(ip, jp) - elev0(fi, fj)), deltaH)) .and. &
1646 6623962 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1647 : ) then
1648 470784 : nn = nn + 1
1649 470784 : ss(nn, 1) = ip
1650 470784 : ss(nn, 2) = jp
1651 : !print *, ip,jp
1652 : end if
1653 : end if
1654 :
1655 : !S
1656 6623962 : if ((ip <= nrows) .and. &
1657 : (jp <= ncols)) then
1658 0 : if ((fdir0(ip, jj) == 64) .and. &
1659 13247924 : (le((elev0(ip, jj) - elev0(fi, fj)), deltaH)) .and. &
1660 6623962 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1661 : ) then
1662 1125660 : nn = nn + 1
1663 1125660 : ss(nn, 1) = ip
1664 1125660 : ss(nn, 2) = jj
1665 : !print *, ip,j
1666 : end if
1667 : end if
1668 :
1669 : !SW
1670 : if ((ip <= nrows) .and. &
1671 6623962 : (jp <= ncols) .and. &
1672 : (jm >= 1)) then
1673 0 : if ((fdir0(ip, jm) == 128) .and. &
1674 13247924 : (le((elev0(ip, jm) - elev0(fi, fj)), deltaH)) .and. &
1675 6623962 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1676 : ) then
1677 335846 : nn = nn + 1
1678 335846 : ss(nn, 1) = ip
1679 335846 : ss(nn, 2) = jm
1680 : !print *, ip,jm
1681 : end if
1682 : end if
1683 :
1684 : !W
1685 6623962 : if ((jm >= 1) .and. &
1686 : (jp <= ncols)) then
1687 0 : if ((fdir0(ii, jm) == 1) .and. &
1688 13247924 : (le((elev0(ii, jm) - elev0(fi, fj)), deltaH)) .and. &
1689 6623962 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1690 : ) then
1691 912382 : nn = nn + 1
1692 912382 : ss(nn, 1) = ii
1693 912382 : ss(nn, 2) = jm
1694 : !print *, i,jm
1695 : end if
1696 : end if
1697 :
1698 : !NW
1699 : if ((im >= 1) .and. &
1700 6623962 : (jp <= ncols) .and. &
1701 : (jm >= 1)) then
1702 0 : if ((fdir0(im, jm) == 2) .and. &
1703 13247924 : (le((elev0(im, jm) - elev0(fi, fj)), deltaH)) .and. &
1704 6623962 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1705 : ) then
1706 492132 : nn = nn + 1
1707 492132 : ss(nn, 1) = im
1708 492132 : ss(nn, 2) = jm
1709 : !print *, im,jm
1710 : end if
1711 : end if
1712 :
1713 : !N
1714 6623962 : if ((im >= 1) .and. &
1715 : (jp <= ncols)) then
1716 0 : if ((fdir0(im, jj) == 4) .and. &
1717 13247924 : (le((elev0(im, jj) - elev0(fi, fj)), deltaH)) .and. &
1718 6623962 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1719 : ) then
1720 1179748 : nn = nn + 1
1721 1179748 : ss(nn, 1) = im
1722 1179748 : ss(nn, 2) = jj
1723 : !print *, im,j
1724 : end if
1725 : end if
1726 :
1727 : !NE
1728 6623962 : if ((im >= 1) .and. &
1729 : (jp <= ncols)) then
1730 0 : if ((fdir0(im, jp) == 8) .and. &
1731 13247924 : (le((elev0(im, jp) - elev0(fi, fj)), deltaH)) .and. &
1732 6623962 : (ge((elev0(ii, jp) - elev0(fi, fj)), 0.0_dp)) &
1733 : ) then
1734 656846 : nn = nn + 1
1735 656846 : ss(nn, 1) = im
1736 656846 : ss(nn, 2) = jp
1737 : !print *, im,jp
1738 : end if
1739 : end if
1740 :
1741 6623962 : 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 43740628 : 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 45179838 : select case (fDir)
1777 : case(1) !E
1778 1439210 : jCol = jCol + 1
1779 : case(2) !SE
1780 1682387 : iRow = iRow + 1
1781 1682387 : jCol = jCol + 1
1782 : case(4) !S
1783 3804397 : iRow = iRow + 1
1784 : case(8) !SW
1785 5753487 : iRow = iRow + 1
1786 5753487 : jCol = jCol - 1
1787 : case(16) !W
1788 11095532 : jCol = jCol - 1
1789 : case(32) !NW
1790 10002909 : iRow = iRow - 1
1791 10002909 : jCol = jCol - 1
1792 : case(64) !N
1793 7261091 : iRow = iRow - 1
1794 : case(128) !NE
1795 2701591 : iRow = iRow - 1
1796 2701591 : jCol = jCol + 1
1797 : case default !sink
1798 : ! do nothing
1799 : end select
1800 :
1801 6623962 : 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 35952 : subroutine cellLength(iDomain, fDir, iRow, jCol, iCoorSystem, length)
1832 :
1833 43740628 : use mo_common_types, only: Grid
1834 : use mo_common_variables, only : domainMeta, level0
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 35952 : real(dp) :: lat_1, long_1, lat_2, long_2
1854 :
1855 : type(Grid), pointer :: level0_iDomain => null()
1856 :
1857 :
1858 35952 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
1859 :
1860 : ! regular X-Y cordinate system
1861 35952 : IF(iCoorSystem .EQ. 0) THEN
1862 :
1863 54786 : select case (fDir)
1864 : case(1, 4, 16, 64) ! E, S, W, N
1865 18834 : length = 1.0_dp
1866 : case(2, 8, 32, 128) ! SE, SW, NW, NE
1867 35952 : length = SQRT2_dp
1868 : end select
1869 35952 : length = length * level0_iDomain%cellsize
1870 :
1871 : ! regular lat-lon cordinate system
1872 0 : ELSE IF(iCoorSystem .EQ. 1) THEN
1873 0 : iRow_to = iRow
1874 0 : jCol_to = jCol
1875 :
1876 : ! move in the direction of flow
1877 0 : 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 : 0.5_dp * level0_iDomain%cellsize
1882 : long_1 = level0_iDomain%xllcorner + real((iRow - 1), dp) * level0_iDomain%cellsize + &
1883 0 : 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 : 0.5_dp * level0_iDomain%cellsize
1887 : long_2 = level0_iDomain%xllcorner + real((iRow_to - 1), dp) * level0_iDomain%cellsize + &
1888 0 : 0.5_dp * level0_iDomain%cellsize
1889 : ! get distance between two points
1890 0 : call get_distance_two_lat_lon_points(lat_1, long_1, lat_2, long_2, length)
1891 :
1892 : END IF
1893 : !
1894 35952 : 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 0 : subroutine get_distance_two_lat_lon_points(lat1, long1, lat2, long2, distance_out)
1927 :
1928 35952 : 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 0 : real(dp) :: theta1
1939 :
1940 0 : real(dp) :: phi1
1941 :
1942 0 : real(dp) :: theta2
1943 :
1944 0 : real(dp) :: phi2
1945 :
1946 : real(dp) :: dtor
1947 :
1948 0 : real(dp) :: term1
1949 :
1950 0 : real(dp) :: term2
1951 :
1952 0 : real(dp) :: term3
1953 :
1954 0 : real(dp) :: temp
1955 :
1956 :
1957 0 : dtor = TWOPI_dp / 360.0_dp
1958 0 : theta1 = dtor * long1
1959 0 : phi1 = dtor * lat1
1960 0 : theta2 = dtor * long2
1961 0 : phi2 = dtor * lat2
1962 :
1963 0 : term1 = cos(phi1) * cos(theta1) * cos(phi2) * cos(theta2)
1964 0 : term2 = cos(phi1) * sin(theta1) * cos(phi2) * sin(theta2)
1965 0 : term3 = sin(phi1) * sin(phi2)
1966 0 : temp = term1 + term2 + term3
1967 0 : if(temp .GT. 1.0_dp) temp = 1.0_dp
1968 :
1969 0 : distance_out = RadiusEarth_dp * acos(temp);
1970 :
1971 0 : 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 24 : subroutine L11_flow_accumulation(iDomain)
2023 :
2024 : use mo_mrm_global_variables, only: &
2025 0 : L11_fDir, & ! IN: flow direction at L11 (standard notation)
2026 : L11_fAcc ! OUT: flow accumulation at L11 [km^2]
2027 : use mo_mrm_global_variables, only : level11
2028 : use mo_common_constants, only : nodata_i4, nodata_dp
2029 : use mo_append, only : append
2030 :
2031 : implicit none
2032 :
2033 : ! local variables
2034 : integer(i4), intent(in) :: iDomain
2035 24 : 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 24 : integer(i4), dimension(:,:), allocatable :: fDir11 ! L11_fDir array
2040 24 : logical, dimension(:,:), allocatable :: mask11 ! Domain mask
2041 :
2042 : ! initialize Domain info
2043 24 : nrows11 = level11(iDomain)%nrows
2044 24 : ncols11 = level11(iDomain)%ncols
2045 24 : s11 = level11(iDomain)%iStart
2046 24 : e11 = level11(iDomain)%iEnd
2047 1860 : mask11 = level11(iDomain)%mask
2048 :
2049 : ! allocate arrays
2050 96 : allocate(fDir11 (nrows11, ncols11))
2051 96 : allocate(fAcc11 (nrows11, ncols11))
2052 :
2053 : ! initialize
2054 1812 : fDir11(:,:) = nodata_i4
2055 :
2056 : ! get data
2057 24 : fDir11(:,:) = UNPACK( L11_fDir(s11:e11), mask11, nodata_i4 )
2058 :
2059 : ! initialize fAcc11 with cell area
2060 960 : fAcc11 = UNPACK( level11(iDomain)%cellarea * 1.e-6, mask11, nodata_dp )
2061 :
2062 : ! For each sink call "calculate_L11_flow_accumulation"
2063 250 : do jj=1, ncols11
2064 1812 : do ii=1, nrows11
2065 1788 : 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 24 : ncol = ncols11)
2072 : end if
2073 : end do
2074 : end do
2075 :
2076 : ! Append
2077 24 : call append( L11_fAcc, pack(fAcc11(:,:),mask11))
2078 :
2079 : ! free space
2080 24 : deallocate(fDir11, fAcc11, mask11)
2081 :
2082 : contains
2083 :
2084 936 : 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 936 : if (jj+1 .le. ncol) then
2104 912 : if (fDir(ii,jj+1) .eq. 16_i4) then
2105 402 : call calculate_L11_flow_accumulation(fDir, fAcc, ii, jj+1, nrow, ncol)
2106 402 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii,jj+1)
2107 : end if
2108 : end if
2109 :
2110 936 : if ((ii+1 .le. nrow) .and. (jj+1 .le. ncol)) then
2111 798 : if (fDir(ii+1,jj+1) .eq. 32_i4) then
2112 0 : call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj+1, nrow, ncol)
2113 0 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj+1)
2114 : end if
2115 : end if
2116 :
2117 936 : if (ii+1 .le. nrow) then
2118 822 : if (fDir(ii+1,jj) .eq. 64_i4) then
2119 320 : call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj, nrow, ncol)
2120 320 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj)
2121 : end if
2122 : end if
2123 :
2124 936 : if ((ii+1 .le. nrow) .and. (jj-1 .ge. 1)) then
2125 770 : if (fDir(ii+1,jj-1) .eq. 128_i4) then
2126 0 : call calculate_L11_flow_accumulation(fDir, fAcc, ii+1, jj-1, nrow, ncol)
2127 0 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii+1,jj-1)
2128 : end if
2129 : end if
2130 :
2131 936 : if (jj-1 .ge. 1) then
2132 884 : if (fDir(ii,jj-1) .eq. 1_i4) then
2133 50 : call calculate_L11_flow_accumulation(fDir, fAcc, ii, jj-1, nrow, ncol)
2134 50 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii,jj-1)
2135 : end if
2136 : end if
2137 :
2138 936 : if ((ii-1 .ge. 1) .and. (jj-1 .ge. 1)) then
2139 842 : if (fDir(ii-1,jj-1) .eq. 2_i4) then
2140 0 : call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj-1, nrow, ncol)
2141 0 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii-1,jj-1)
2142 : end if
2143 : end if
2144 :
2145 936 : if (ii-1 .ge. 1) then
2146 892 : if (fDir(ii-1,jj) .eq. 4_i4) then
2147 138 : call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj, nrow, ncol)
2148 138 : fAcc(ii,jj) = fAcc(ii,jj) + fAcc(ii-1,jj)
2149 : end if
2150 : end if
2151 :
2152 936 : if ((ii-1 .ge. 1) .and. (jj+1 .le. ncol)) then
2153 868 : if (fDir11(ii-1,jj+1) .eq. 8_i4) then
2154 2 : call calculate_L11_flow_accumulation(fDir, fAcc, ii-1, jj+1, nrow, ncol)
2155 2 : 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 24 : 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 17 : subroutine L11_calc_celerity(iDomain, param)
2213 : use mo_common_constants, only: nodata_i4, nodata_dp
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 17 : logical, dimension(:,:), allocatable :: mask0
2250 17 : integer(i4), dimension(:,:), allocatable :: iD0
2251 17 : integer(i4), dimension(:,:), allocatable :: fDir0
2252 17 : integer(i4), dimension(:,:), allocatable :: fAcc0
2253 17 : real(dp), dimension(:,:), allocatable :: slope0
2254 17 : real(dp), dimension(:), allocatable :: slope_tmp
2255 17 : real(dp), dimension(:,:), allocatable :: cellarea0
2256 17 : integer(i4), dimension(:), allocatable :: netPerm ! routing order (permutation)
2257 17 : integer(i4), dimension(:), allocatable :: nLinkFromRow
2258 17 : integer(i4), dimension(:), allocatable :: nLinkFromCol
2259 17 : integer(i4), dimension(:), allocatable :: nLinkToRow
2260 17 : integer(i4), dimension(:), allocatable :: nLinkToCol
2261 : integer(i4) :: ii, rr, ns
2262 : integer(i4) :: frow, fcol
2263 : integer(i4) :: fId, tId
2264 17 : real(dp), dimension(:), allocatable :: stack, append_chunk ! Stacks celerity along the L0 river-path
2265 17 : integer(i4), dimension(:), allocatable :: dummy_1d
2266 :
2267 17 : real(dp) :: L0_link_slope
2268 17 : real(dp), dimension(:), allocatable :: celerity11
2269 17 : real(dp), dimension(:,:), allocatable :: celerity0
2270 :
2271 17 : integer(i4), dimension(:,:), allocatable :: nodata_i4_tmp
2272 17 : real(dp), dimension(:,:), allocatable :: nodata_dp_tmp
2273 17 : logical, dimension(:), allocatable :: slopemask0
2274 :
2275 : type(Grid), pointer :: level0_iDomain
2276 :
2277 : ! level-0 information
2278 17 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
2279 17 : nrows0 = level0_iDomain%nrows
2280 17 : ncols0 = level0_iDomain%ncols
2281 17 : nCells0 = level0_iDomain%ncells
2282 17 : iStart0 = level0_iDomain%iStart
2283 17 : iEnd0 = level0_iDomain%iEnd
2284 2122467 : mask0 = level0_iDomain%mask
2285 :
2286 : ! level-11 information
2287 17 : iStart11 = level11(iDomain)%iStart
2288 17 : iEnd11 = level11(iDomain)%iEnd
2289 17 : nrows11 = level11(iDomain)%nrows
2290 17 : ncols11 = level11(iDomain)%ncols
2291 17 : nNodes = level11(iDomain)%ncells
2292 :
2293 17 : nLinks = nNodes - L11_nOutlets(iDomain)
2294 :
2295 : ! allocate
2296 68 : allocate ( iD0 ( nrows0, ncols0 ) )
2297 68 : allocate ( slope0 ( nrows0, ncols0 ) )
2298 51 : allocate ( fDir0 ( nrows0, ncols0 ) )
2299 51 : allocate ( fAcc0 ( nrows0, ncols0 ) )
2300 51 : allocate ( cellarea0 ( nrows0, ncols0 ) )
2301 51 : allocate ( celerity0 ( nrows0, ncols0 ) )
2302 51 : 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 17 : allocate ( stack ( 1 ) )
2307 17 : allocate ( append_chunk( 1 ) )
2308 17 : allocate ( dummy_1d ( 2 ))
2309 51 : allocate ( netPerm ( nNodes ) )
2310 34 : allocate ( nLinkFromRow( nNodes ) )
2311 34 : allocate ( nLinkFromCol( nNodes ) )
2312 34 : allocate ( nLinkToRow ( nNodes ) )
2313 34 : allocate ( nLinkToCol ( nNodes ) )
2314 51 : allocate ( celerity11 ( nNodes ) )
2315 34 : allocate ( slope_tmp ( nNodes ) )
2316 :
2317 51 : allocate (nodata_i4_tmp( nrows0, ncols0 ) )
2318 51 : allocate (nodata_dp_tmp( nrows0, ncols0 ) )
2319 :
2320 : ! initialize
2321 2122433 : iD0(:,:) = nodata_i4
2322 2122433 : fDir0(:,:) = nodata_i4
2323 2122433 : fAcc0(:,:) = nodata_i4
2324 2122433 : cellarea0(:,:) = nodata_dp
2325 2122433 : slope0(:,:) = nodata_dp
2326 :
2327 34 : stack(:) = nodata_dp
2328 34 : append_chunk(:) = nodata_dp
2329 595 : netPerm(:) = nodata_i4
2330 595 : nLinkFromRow(:) = nodata_i4
2331 595 : nLinkFromCol(:) = nodata_i4
2332 595 : nLinkToRow(:) = nodata_i4
2333 595 : nLinkToCol(:) = nodata_i4
2334 595 : celerity11(:) = nodata_dp
2335 2122433 : celerity0(:,:) = nodata_dp
2336 791282 : slopemask0(:) = .False.
2337 :
2338 2122433 : nodata_i4_tmp(:,:) = nodata_i4
2339 2122433 : nodata_dp_tmp(:,:) = nodata_dp
2340 :
2341 : ! for a single node model run
2342 :
2343 17 : if(nNodes .GT. 1) then
2344 : ! get L0 fields
2345 17 : iD0(:,:) = UNPACK(level0_iDomain%Id(1:nCells0), mask0, nodata_i4_tmp)
2346 17 : fDir0(:,:) = UNPACK(L0_fDir(iStart0:iEnd0), mask0, nodata_i4_tmp)
2347 17 : fAcc0(:,:) = UNPACK(L0_fAcc(iStart0:iEnd0), mask0, nodata_i4_tmp)
2348 17 : cellarea0(:,:) = UNPACK(level0_iDomain%cellarea(1:nCells0), mask0, nodata_dp_tmp)
2349 :
2350 : ! smoothing river slope
2351 791299 : slope_tmp = L0_slope(iStart0:iEnd0)
2352 791282 : where ( slope_tmp .lt. 0.1_dp ) slope_tmp = 0.1_dp
2353 :
2354 791282 : slopemask0(:) = (L0_streamNet(iStart0:iEnd0) .ne. nodata_i4)
2355 :
2356 : ! smooth river cells if there is more than one cell
2357 791282 : if( count(slopemask0) .GT. 1) then
2358 1582564 : slope_tmp = mad(arr = slope_tmp, z = 2.25_dp, mask = slopemask0, tout='u', mval=0.1_dp)
2359 : end if
2360 17 : slope0(:,:) = UNPACK(slope_tmp, mask0, nodata_dp_tmp )
2361 :
2362 : ! get network vectors of L11
2363 595 : netPerm(:) = L11_netPerm ( iStart11 : iEnd11 )
2364 595 : nLinkFromRow(:) = L11_fRow ( iStart11 : iEnd11 )
2365 595 : nLinkFromCol(:) = L11_fCol ( iStart11 : iEnd11 )
2366 595 : nLinkToRow(:) = L11_tRow ( iStart11 : iEnd11 )
2367 595 : nLinkToCol(:) = L11_tCol ( iStart11 : iEnd11 )
2368 :
2369 578 : do rr = 1, nLinks
2370 :
2371 561 : ii = netPerm(rr)
2372 561 : frow = nLinkFromRow(ii)
2373 561 : fcol = nLinkFromCol(ii)
2374 :
2375 : ! Init
2376 1122 : stack(:) = 0_dp
2377 561 : ns = 1
2378 :
2379 561 : fId = iD0( frow, fcol )
2380 561 : tId = iD0( nLinkToRow(ii) , nLinkToCol(ii) )
2381 24514 : do
2382 25075 : L0_link_slope = slope0(frow, fcol) / 100._dp
2383 :
2384 : ! celerity parametrization
2385 25075 : stack(ns) = param(1) * sqrt(L0_link_slope)
2386 :
2387 25075 : celerity0(frow, fcol) = stack(ns)
2388 25075 : ns = ns + 1
2389 25075 : fId = iD0(frow, fcol)
2390 25075 : if( .NOT. (fID == tID)) then
2391 24514 : call append(stack, append_chunk)
2392 : else
2393 : exit
2394 : end if
2395 : ! move downstream
2396 24514 : call moveDownOneCell( fDir0(frow,fcol), frow, fcol )
2397 : end do
2398 :
2399 25636 : celerity11(ii) = size(stack) / sum(1/stack(:))
2400 561 : deallocate(stack)
2401 578 : 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 0 : celerity11(:) = 1._dp
2410 :
2411 : end if
2412 :
2413 : ! Write celerity
2414 595 : L11_celerity(iStart11:iEnd11) = celerity11(:)
2415 17 : L0_celerity(iStart0:iEnd0) = PACK(celerity0(:,:), mask0)
2416 :
2417 : ! free space
2418 0 : deallocate (&
2419 0 : mask0, iD0, slope_tmp, slopemask0, &
2420 0 : slope0, fDir0, cellarea0, &
2421 17 : stack, netPerm, nLinkFromRow, nLinkFromCol, nLinkToRow, nLinkToCol)
2422 :
2423 17 : end subroutine L11_calc_celerity
2424 :
2425 : end module mo_mrm_net_startup
|
