Line data Source code
1 : !> \file mo_mrm_restart.f90
2 : !> \brief \copybrief mo_mrm_restart
3 : !> \details \copydetails mo_mrm_restart
4 :
5 : !> \brief Restart routines
6 : !> \details This module contains the subroutines for reading and writing routing related variables to file.
7 : !> \authors Stephan Thober
8 : !> \date Aug 2015
9 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
10 : !! mHM is released under the LGPLv3+ license \license_note
11 : !> \ingroup f_mrm
12 : module mo_mrm_restart
13 : use mo_kind, only : i4, dp
14 : use mo_message, only : message, error_message
15 :
16 : implicit none
17 :
18 : public :: mrm_write_restart
19 : public :: mrm_read_restart_states
20 : public :: mrm_read_restart_config
21 : contains
22 :
23 : ! ------------------------------------------------------------------
24 :
25 : ! NAME
26 : ! mrm_write_restart
27 :
28 : ! PURPOSE
29 : !> \brief write routing states and configuration
30 :
31 : !> \details write configuration and state variables to a given restart
32 : !> directory.
33 :
34 : ! INTENT(IN)
35 : !> \param[in] "integer(i4) :: iDomain" number of domain
36 : !> \param[in] "character(256), dimension(:) :: OutFile" list of Output paths per Domain
37 :
38 : ! HISTORY
39 : !> \authors Stephan Thober
40 :
41 : !> \date Aug 2015
42 :
43 : ! Modifications:
44 : ! Stephan Thober Sep 2015 - added all write restart commands in this subroutine
45 : ! Stephan Thober Sep 2015 - added L11_areaCell L1_ID and L1_L11_Id for routing
46 : ! resolution higher than hydrology resolution
47 : ! David Schaefer Nov 2015 - mo_netcdf
48 : ! Stephan Thober May 2016 - split L0_OutletCoord into L0_rowOutlet & L0_colOutlet
49 : ! because multiple outlets could exist
50 : ! Stephan Thober Nov 2016 - added L11_TSrout, ProcessMatrix
51 : ! Matthias Kelbling Aug 2017 - added L11_fAcc, L0_slope, L0_celerity, L11_celerity, L11_meandering
52 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
53 : ! Stephan Thober Jun 2018 - including varying celerity functionality
54 : ! Stephan Thober May 2019 - added L0 info required for Process 3
55 :
56 41 : subroutine mrm_write_restart(iDomain, domainID, OutFile)
57 :
58 : use mo_common_constants, only : nodata_dp, nodata_i4
59 : use mo_common_restart, only : write_grid_info
60 : use mo_common_variables, only : level0, level1, nLCoverScene, processMatrix, domainMeta, &
61 : LC_year_start, LC_year_end
62 : use mo_common_constants, only : landCoverPeriodsVarName
63 : use mo_mrm_constants, only : nRoutingStates
64 : use mo_mpr_global_variables, only : L0_slope
65 : use mo_mrm_global_variables, only : L0_fdir, L0_fAcc, L0_streamnet, &
66 : L1_L11_Id, &
67 : L11_C1, L11_C2, L11_K, L11_L1_Id, L11_Qmod, &
68 : L11_TSrout, L11_aFloodPlain, L11_colOut, L11_colOut, L11_fCol, L11_fDir, &
69 : L11_fAcc, L11_fRow, L11_fromN, L11_label, L11_length, L11_nLinkFracFPimp, &
70 : L11_netPerm, L11_qOUT, L11_qTIN, L11_qTR, L11_rOrder, L11_rowOut, L11_rowOut, &
71 : L11_sink, L11_slope, L11_tCol, L11_tRow, L11_toN, L11_xi, L11_celerity, &
72 : level11, domain_mrm
73 : use mo_netcdf, only : NcDataset, NcDimension, NcVariable
74 : use mo_string_utils, only : num2str
75 :
76 : implicit none
77 :
78 : ! number of domain
79 : integer(i4), intent(in) :: iDomain
80 :
81 : ! domain
82 : integer(i4), intent(in) :: domainID
83 :
84 : ! list of Output paths per Domain
85 : character(256), dimension(:), intent(in) :: OutFile
86 :
87 : character(256) :: Fname
88 :
89 : integer(i4) :: ii
90 :
91 : ! number of outlets at Level 0
92 : integer(i4) :: Noutlet
93 :
94 : ! start index at level 0
95 : integer(i4) :: s0
96 :
97 : ! end index at level 0
98 : integer(i4) :: e0
99 :
100 : ! mask at level 0
101 20 : logical, dimension(:, :), allocatable :: mask0
102 :
103 : ! start index at level 1
104 : integer(i4) :: s1
105 :
106 : ! end index at level 1
107 : integer(i4) :: e1
108 :
109 : ! mask at level 1
110 20 : logical, dimension(:, :), allocatable :: mask1
111 :
112 : ! start index at level 11
113 : integer(i4) :: s11
114 :
115 : ! end index at level 11
116 : integer(i4) :: e11
117 :
118 : ! number of colums at level 11
119 : integer(i4) :: ncols11
120 :
121 : ! number of rows at level 11
122 : integer(i4) :: nrows11
123 :
124 : ! dummy variable for writing L11_sink
125 20 : integer(i4), allocatable :: dummy(:)
126 :
127 : ! mask at level 11
128 20 : logical, dimension(:, :), allocatable :: mask11
129 :
130 : ! dummy variable
131 20 : real(dp), dimension(:, :, :), allocatable :: dummy_d3
132 20 : real(dp), dimension(:), allocatable :: dummy_d1
133 :
134 : type(NcDataset) :: nc
135 :
136 : type(NcDimension) :: rows0, cols0, rows1, cols1, rows11, cols11, it11, lcscenes, nout
137 :
138 : type(NcDimension) :: links, nts, nproc
139 :
140 : type(NcVariable) :: var
141 :
142 :
143 : ! get Level1 and Level11 information about the Domain
144 20 : noutlet = domain_mrm(domainMeta%L0DataFrom(iDomain))%L0_noutlet
145 20 : s0 = level0(domainMeta%L0DataFrom(iDomain))%iStart
146 20 : e0 = level0(domainMeta%L0DataFrom(iDomain))%iEnd
147 2363004 : mask0 = level0(domainMeta%L0DataFrom(iDomain))%mask
148 20 : s1 = level1(iDomain)%iStart
149 20 : e1 = level1(iDomain)%iEnd
150 1787 : mask1 = level1(iDomain)%mask
151 20 : s11 = level11(iDomain)%iStart
152 20 : e11 = level11(iDomain)%iEnd
153 1616 : mask11 = level11(iDomain)%mask
154 20 : ncols11 = level11(iDomain)%ncols
155 20 : nrows11 = level11(iDomain)%nrows
156 100 : allocate(dummy_d3(nrows11, ncols11, nRoutingStates))
157 :
158 : ! set restart file name
159 20 : Fname = trim(OutFile(iDomain))
160 :
161 20 : call message(' Writing mRM restart file to ' // trim(Fname) // ' ...')
162 :
163 20 : nc = NcDataset(fname, "w")
164 :
165 20 : call write_grid_info(level0(domainMeta%L0DataFrom(iDomain)), "0", nc)
166 20 : call write_grid_info(level1(iDomain), "1", nc)
167 20 : call write_grid_info(level11(iDomain), "11", nc)
168 :
169 20 : nout = nc%setDimension("Noutlet", Noutlet)
170 20 : rows0 = nc%getDimension("nrows0")
171 20 : cols0 = nc%getDimension("ncols0")
172 20 : rows1 = nc%getDimension("nrows1")
173 20 : cols1 = nc%getDimension("ncols1")
174 20 : rows11 = nc%getDimension("nrows11")
175 20 : cols11 = nc%getDimension("ncols11")
176 :
177 20 : it11 = nc%setDimension("nIT", nRoutingStates)
178 20 : links = nc%setDimension("nLinks", size(L11_length(s11 : e11)))
179 20 : nts = nc%setDimension("TS", 1)
180 20 : nproc = nc%setDimension("Nprocesses", size(processMatrix, dim = 1))
181 60 : allocate(dummy_d1(nLCoverScene+1))
182 60 : dummy_d1(1:nLCoverScene) = LC_year_start(:)
183 : ! this is done because bounds are always stored as real so e.g.
184 : ! 1981-1990,1991-2000 is thus saved as 1981.0-1991.0,1991.0-2001.0
185 : ! it is translated back into ints correctly during reading
186 20 : dummy_d1(nLCoverScene+1) = LC_year_end(nLCoverScene) + 1
187 20 : lcscenes = nc%setCoordinate(trim(landCoverPeriodsVarName), nLCoverScene, dummy_d1, 0_i4)
188 20 : deallocate(dummy_d1)
189 :
190 :
191 : ! add processMatrix
192 40 : var = nc%setVariable("ProcessMatrix", "i32", (/nproc/))
193 20 : call var%setFillValue(nodata_i4)
194 20 : call var%setData(processMatrix(:, 1))
195 20 : call var%setAttribute("long_name", "Process Matrix")
196 :
197 : ! add L0 variables if processmatrix is equal to 3
198 20 : if (processMatrix(8, 1) .eq. 3_i4) then
199 : ! add L0_fdir, L0_fAcc, L0_slope, L0_streamnet
200 9 : var = nc%setVariable("L0_fDir", "i32", (/rows0, cols0/))
201 3 : call var%setFillValue(nodata_i4)
202 3 : call var%setData(unpack(L0_fdir(s0:e0), mask0, nodata_i4))
203 3 : call var%setAttribute("long_name", "flow direction at level 0")
204 :
205 9 : var = nc%setVariable("L0_fAcc", "i32", (/rows0, cols0/))
206 3 : call var%setFillValue(nodata_i4)
207 3 : call var%setData(unpack(L0_fAcc(s0:e0), mask0, nodata_i4))
208 3 : call var%setAttribute("long_name", "flow accumulation at level 0")
209 :
210 9 : var = nc%setVariable("L0_slope", "f64", (/rows0, cols0/))
211 3 : call var%setFillValue(nodata_dp)
212 3 : call var%setData(unpack(L0_slope(s0:e0), mask0, nodata_dp))
213 3 : call var%setAttribute("long_name", "slope at level 0")
214 :
215 9 : var = nc%setVariable("L0_streamnet", "i32", (/rows0, cols0/))
216 3 : call var%setFillValue(nodata_i4)
217 3 : call var%setData(unpack(L0_streamnet(s0:e0), mask0, nodata_i4))
218 3 : call var%setAttribute("long_name", "streamnet at level 0")
219 : end if
220 :
221 60 : var = nc%setVariable("L1_Id", "i32", (/rows1, cols1/))
222 20 : call var%setFillValue(nodata_i4)
223 20 : call var%setData(unpack(level1(iDomain)%Id(1:e1-s1+1), mask1, nodata_i4))
224 20 : call var%setAttribute("long_name", "cell IDs at level 1")
225 :
226 60 : var = nc%setVariable("L1_L11_Id", "i32", (/rows1, cols1/))
227 20 : call var%setFillValue(nodata_i4)
228 20 : call var%setData(unpack(L1_L11_Id(s1 : e1), mask1, nodata_i4))
229 20 : call var%setAttribute("long_name", "Mapping of L1 Id on L11")
230 :
231 60 : var = nc%setVariable("L11_Qmod", "f64", (/rows11, cols11/))
232 20 : call var%setFillValue(nodata_dp)
233 20 : call var%setData(unpack(L11_Qmod(s11 : e11), mask11, nodata_dp))
234 20 : call var%setAttribute("long_name", "simulated discharge at each node at level 11")
235 :
236 60 : var = nc%setVariable("L11_qOUT", "f64", (/rows11, cols11/))
237 20 : call var%setFillValue(nodata_dp)
238 20 : call var%setData(unpack(L11_qOUT(s11 : e11), mask11, nodata_dp))
239 20 : call var%setAttribute("long_name", "Total outflow from cells L11 at time tt at level 11")
240 :
241 60 : do ii = 1, size(dummy_d3, 3)
242 60 : dummy_d3(:, :, ii) = unpack(L11_qTIN(s11 : e11, ii), mask11, nodata_dp)
243 : end do
244 80 : var = nc%setVariable("L11_qTIN", "f64", (/rows11, cols11, it11/))
245 20 : call var%setFillValue(nodata_dp)
246 20 : call var%setData(dummy_d3)
247 20 : call var%setAttribute("long_name", "Total discharge inputs at t-1 and t at level 11")
248 :
249 60 : do ii = 1, size(dummy_d3, 3)
250 60 : dummy_d3(:, :, ii) = unpack(L11_qTR(s11 : e11, ii), mask11, nodata_dp)
251 : end do
252 80 : var = nc%setVariable("L11_qTR", "f64", (/rows11, cols11, it11/))
253 20 : call var%setFillValue(nodata_dp)
254 20 : call var%setData(dummy_d3)
255 20 : call var%setAttribute("long_name", "Routed outflow leaving a node at level 11")
256 :
257 60 : var = nc%setVariable("L11_K", "f64", (/rows11, cols11/))
258 20 : call var%setFillValue(nodata_dp)
259 20 : call var%setData(unpack(L11_K(s11 : e11), mask11, nodata_dp))
260 20 : call var%setAttribute("long_name", "kappa: Muskingum travel time parameter at level 11")
261 :
262 60 : var = nc%setVariable("L11_xi", "f64", (/rows11, cols11/))
263 20 : call var%setFillValue(nodata_dp)
264 20 : call var%setData(unpack(L11_xi(s11 : e11), mask11, nodata_dp))
265 20 : call var%setAttribute("long_name", "xi: Muskingum diffusion parameter at level 11")
266 :
267 60 : var = nc%setVariable("L11_C1", "f64", (/rows11, cols11/))
268 20 : call var%setFillValue(nodata_dp)
269 20 : call var%setData(unpack(L11_C1(s11 : e11), mask11, nodata_dp))
270 20 : call var%setAttribute("long_name", "Routing parameter C1=f(K,xi, DT) (Chow, 25-41) at level 11")
271 :
272 60 : var = nc%setVariable("L11_C2", "f64", (/rows11, cols11/))
273 20 : call var%setFillValue(nodata_dp)
274 20 : call var%setData(unpack(L11_C2(s11 : e11), mask11, nodata_dp))
275 20 : call var%setAttribute("long_name", "Routing parameter C2=f(K,xi, DT) (Chow, 25-41) at level 11")
276 :
277 20 : deallocate(dummy_d3)
278 100 : allocate(dummy_d3(nrows11, ncols11, nLCoverScene))
279 80 : do ii = 1, size(dummy_d3, 3)
280 60 : dummy_d3(:, :, ii) = unpack(L11_nLinkFracFPimp(s11 : e11, ii), mask11, nodata_dp)
281 : end do
282 80 : var = nc%setVariable("L11_nLinkFracFPimp", "f64", (/rows11, cols11, lcscenes/))
283 20 : call var%setFillValue(nodata_dp)
284 20 : call var%setData(dummy_d3)
285 20 : call var%setAttribute("long_name", "Fraction of the flood plain with impervious cover at level 11")
286 :
287 : ! ----------------------------------------------------------
288 : ! L11 config set
289 : ! ----------------------------------------------------------
290 60 : var = nc%setVariable("L11_domain_Mask", "i32", (/rows11, cols11/))
291 20 : call var%setFillValue(nodata_i4)
292 1556 : call var%setData(merge(1_i4, 0_i4, mask11))
293 20 : call var%setAttribute("long_name", "Mask at Level 11")
294 :
295 40 : var = nc%setVariable("L11_TSrout", "i32", (/nts/))
296 20 : call var%setFillValue(nodata_i4)
297 20 : call var%setData(L11_TSrout(iDomain))
298 20 : call var%setAttribute("long_name", "routing resolution at Level 11")
299 20 : call var%setAttribute("units", "s")
300 :
301 60 : var = nc%setVariable("L11_Id", "i32", (/rows11, cols11/))
302 20 : call var%setFillValue(nodata_i4)
303 20 : call var%setData(unpack(level11(iDomain)%Id(1:e11-s11+1), mask11, nodata_i4))
304 20 : call var%setAttribute("long_name", "cell Ids at Level 11")
305 :
306 60 : var = nc%setVariable("L11_fDir", "i32", (/rows11, cols11/))
307 20 : call var%setFillValue(nodata_i4)
308 20 : call var%setData(unpack(L11_fDir(s11 : e11), mask11, nodata_i4))
309 20 : call var%setAttribute("long_name", "flow Direction at Level 11")
310 :
311 60 : var = nc%setVariable("L11_fAcc", "f64", (/rows11, cols11/))
312 20 : call var%setFillValue(nodata_dp)
313 20 : call var%setData(unpack(L11_fAcc(s11:e11), mask11, nodata_dp))
314 20 : call var%setAttribute("long_name", "flow accumulation at Level 11")
315 :
316 60 : var = nc%setVariable("L11_rowOut", "i32", (/rows11, cols11/))
317 20 : call var%setFillValue(nodata_i4)
318 20 : call var%setData(unpack(L11_rowOut(s11 : e11), mask11, nodata_i4))
319 20 : call var%setAttribute("long_name", "Grid vertical location of the Outlet at Level 11")
320 :
321 60 : var = nc%setVariable("L11_colOut", "i32", (/rows11, cols11/))
322 20 : call var%setFillValue(nodata_i4)
323 20 : call var%setData(unpack(L11_colOut(s11 : e11), mask11, nodata_i4))
324 20 : call var%setAttribute("long_name", "Grid horizontal location of the Outlet at Level 11")
325 :
326 40 : var = nc%setVariable("L11_fromN", "i32", (/links/))
327 20 : call var%setFillValue(nodata_i4)
328 20 : call var%setData(L11_fromN(s11 : e11))
329 20 : call var%setAttribute("long_name", "From Node")
330 :
331 40 : var = nc%setVariable("L11_toN", "i32", (/links/))
332 20 : call var%setFillValue(nodata_i4)
333 20 : call var%setData(L11_toN(s11 : e11))
334 20 : call var%setAttribute("long_name", "To Node")
335 :
336 40 : var = nc%setVariable("L11_rOrder", "i32", (/links/))
337 20 : call var%setFillValue(nodata_i4)
338 20 : call var%setData(L11_rOrder(s11 : e11))
339 20 : call var%setAttribute("long_name", "Network routing order at Level 11")
340 :
341 40 : var = nc%setVariable("L11_label", "i32", (/links/))
342 20 : call var%setFillValue(nodata_i4)
343 20 : call var%setData(L11_label(s11 : e11))
344 20 : call var%setAttribute("long_name", "Label Id [0='', 1=HeadWater, 2=Sink] at Level 11")
345 :
346 40 : var = nc%setVariable("L11_sink", "i32", (/links/))
347 20 : call var%setFillValue(nodata_i4)
348 80 : allocate(dummy(e11 - s11 + 1))
349 840 : dummy = 0_i4
350 820 : where(L11_sink(s11 : e11))
351 : dummy = 1_i4
352 : end where
353 : ! call var%setData(merge(1_i4, 0_i4, L11_sink(s11 : e11)))
354 20 : call var%setData(dummy)
355 20 : deallocate(dummy)
356 20 : call var%setAttribute("long_name", ".true. if sink node reached at Level 11")
357 :
358 40 : var = nc%setVariable("L11_netPerm", "i32", (/links/))
359 20 : call var%setFillValue(nodata_i4)
360 20 : call var%setData(L11_netPerm(s11 : e11))
361 20 : call var%setAttribute("long_name", "Routing sequence (permutation of L11_rOrder) at Level 11")
362 :
363 40 : var = nc%setVariable("L11_fRow", "i32", (/links/))
364 20 : call var%setFillValue(nodata_i4)
365 20 : call var%setData(L11_fRow(s11 : e11))
366 20 : call var%setAttribute("long_name", "From row in L0 grid at Level 11")
367 :
368 40 : var = nc%setVariable("L11_fCol", "i32", (/links/))
369 20 : call var%setFillValue(nodata_i4)
370 20 : call var%setData(L11_fCol(s11 : e11))
371 20 : call var%setAttribute("long_name", "From col in L0 grid at Level 11")
372 :
373 40 : var = nc%setVariable("L11_tRow", "i32", (/links/))
374 20 : call var%setFillValue(nodata_i4)
375 20 : call var%setData(L11_tRow(s11 : e11))
376 20 : call var%setAttribute("long_name", "To row in L0 grid at Level 11")
377 :
378 40 : var = nc%setVariable("L11_tCol", "i32", (/links/))
379 20 : call var%setFillValue(nodata_i4)
380 20 : call var%setData(L11_tCol(s11 : e11))
381 20 : call var%setAttribute("long_name", "To Col in L0 grid at Level 11")
382 :
383 40 : var = nc%setVariable("L11_length", "f64", (/links/))
384 20 : call var%setFillValue(nodata_dp)
385 20 : call var%setData(L11_length(s11 : e11))
386 20 : call var%setAttribute("long_name", "Total length of river link [m]")
387 :
388 40 : var = nc%setVariable("L11_aFloodPlain", "f64", (/links/))
389 20 : call var%setFillValue(nodata_dp)
390 20 : call var%setData(L11_aFloodPlain(s11 : e11))
391 20 : call var%setAttribute("long_name", "Area of the flood plain [m2]")
392 :
393 40 : var = nc%setVariable("L11_slope", "f64", (/links/))
394 20 : call var%setFillValue(nodata_dp)
395 20 : call var%setData(L11_slope(s11 : e11))
396 20 : call var%setAttribute("long_name", "Average slope of river link")
397 :
398 60 : var = nc%setVariable("L11_L1_Id", "i32", (/rows11, cols11/))
399 20 : call var%setFillValue(nodata_i4)
400 20 : call var%setData(unpack(L11_L1_Id(s11 : e11), mask11, nodata_i4))
401 20 : call var%setAttribute("long_name", "Mapping of L1 Id on L11")
402 :
403 : var = nc%setVariable("gaugeNodeList", "i32", &
404 0 : (/nc%setDimension("Ngauges", size(domain_mrm(iDomain)%gaugeNodeList(:)))/) &
405 40 : )
406 20 : call var%setFillValue(nodata_i4)
407 20 : call var%setData(domain_mrm(iDomain)%gaugeNodeList(:))
408 20 : call var%setAttribute("long_name", "cell ID of gauges")
409 :
410 : var = nc%setVariable("InflowGaugeNodeList", "i32", &
411 0 : (/nc%setDimension("nInflowGauges", size(domain_mrm(iDomain)%InflowGaugeNodeList(:)))/) &
412 40 : )
413 20 : call var%setFillValue(nodata_i4)
414 20 : call var%setData(domain_mrm(iDomain)%InflowGaugeNodeList(:))
415 20 : call var%setAttribute("long_name", "cell ID of gauges")
416 :
417 20 : if (processMatrix(8, 1) .eq. 3) then
418 6 : var = nc%setVariable("L11_celerity", "f64", (/links/)) ! celerity
419 3 : call var%setFillValue(nodata_dp)
420 3 : call var%setData(L11_celerity(s11:e11))
421 3 : call var%setAttribute("long_name", "celerity at Level 11")
422 : end if
423 :
424 : ! free dummy variables
425 20 : deallocate(dummy_d3)
426 20 : call nc%close()
427 :
428 80 : end subroutine mrm_write_restart
429 :
430 : ! ------------------------------------------------------------------
431 :
432 : ! NAME
433 : ! mrm_read_restart_states
434 :
435 : ! PURPOSE
436 : !> \brief read routing states
437 :
438 : !> \details This subroutine reads the routing states from
439 : !> mRM_states_<domain_id>.nc that has to be in the given
440 : !> path directory. This subroutine has to be called directly
441 : !> each time the mHM_eval or mRM_eval is called such that the
442 : !> the states are always the same at the first simulation time
443 : !> step, crucial for optimization.
444 :
445 : ! INTENT(IN)
446 : !> \param[in] "integer(i4) :: iDomain" number of domains
447 : !> \param[in] "character(256) :: InFile" Input Path including trailing slash
448 :
449 : ! HISTORY
450 : !> \authors Stephan Thober
451 :
452 : !> \date Sep 2015
453 :
454 : ! Modifications:
455 : ! David Schaefer Mar 2016 - mo_netcdf
456 : ! Stephan Thober May 2016 - split L0_OutletCoord into L0_rowOutlet & L0_colOutlet because multiple outlets could exist
457 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
458 :
459 1 : subroutine mrm_read_restart_states(iDomain, domainID, InFile)
460 :
461 20 : use mo_common_variables, only : nLCoverScene
462 : use mo_mrm_constants, only : nRoutingStates
463 : use mo_mrm_global_variables, only : L11_C1, L11_C2, L11_K, L11_Qmod, L11_Qout, L11_nLinkFracFPimp, L11_qTIN, L11_qTR, &
464 : L11_xi, level11
465 : use mo_netcdf, only : NcDataset, NcVariable
466 : use mo_string_utils, only : num2str
467 :
468 : implicit none
469 :
470 : ! number of Domain
471 : integer(i4), intent(in) :: iDomain
472 :
473 : integer(i4), intent(in) :: domainID
474 :
475 : ! Input Path including trailing slash
476 : character(256), intent(in) :: InFile
477 :
478 : integer(i4) :: ii
479 :
480 : ! start index at level 11
481 : integer(i4) :: s11
482 :
483 : ! end index at level 11
484 : integer(i4) :: e11
485 :
486 : ! mask at level 11
487 1 : logical, dimension(:, :), allocatable :: mask11
488 :
489 : ! dummy, 2 dimension
490 1 : real(dp), dimension(:, :), allocatable :: dummyD2
491 :
492 : ! dummy, 3 dimension
493 1 : real(dp), dimension(:, :, :), allocatable :: dummyD3
494 :
495 : character(256) :: fname
496 :
497 : type(NcDataset) :: nc
498 :
499 : type(NcVariable) :: var
500 :
501 : !TODO-RIV-TEMP: read/write restart for riv-temp process
502 :
503 : ! set file name
504 1 : fname = trim(InFile)
505 :
506 : ! get Domain info at L11 including ncells, start, end, and mask
507 1 : s11 = level11(iDomain)%iStart
508 1 : e11 = level11(iDomain)%iEnd
509 67 : mask11 = level11(iDomain)%mask
510 :
511 1 : nc = NcDataset(fname, "r")
512 :
513 : ! simulated discharge at each node
514 1 : var = nc%getVariable("L11_Qmod")
515 1 : call var%getData(dummyD2)
516 1 : L11_Qmod(s11 : e11) = pack(dummyD2, mask11)
517 :
518 : ! Total outflow from cells L11 at time tt
519 1 : var = nc%getVariable("L11_qOUT")
520 1 : call var%getData(dummyD2)
521 1 : L11_qOUT(s11 : e11) = pack(dummyD2, mask11)
522 :
523 : ! Total discharge inputs at t-1 and t
524 1 : var = nc%getVariable("L11_qTIN")
525 1 : call var%getData(dummyD3)
526 3 : do ii = 1, nRoutingStates
527 3 : L11_qTIN(s11 : e11, ii) = pack(dummyD3(:, :, ii), mask11)
528 : end do
529 :
530 : ! Routed outflow leaving a node
531 1 : var = nc%getVariable("L11_qTR")
532 1 : call var%getData(dummyD3)
533 3 : do ii = 1, nRoutingStates
534 3 : L11_qTR(s11 : e11, ii) = pack(dummyD3(:, :, ii), mask11)
535 : end do
536 :
537 : ! kappa: Muskingum travel time parameter.
538 1 : var = nc%getVariable("L11_K")
539 1 : call var%getData(dummyD2)
540 1 : L11_K(s11 : e11) = pack(dummyD2, mask11)
541 :
542 : ! xi: Muskingum diffusion parameter
543 1 : var = nc%getVariable("L11_xi")
544 1 : call var%getData(dummyD2)
545 1 : L11_xi(s11 : e11) = pack(dummyD2, mask11)
546 :
547 : ! Routing parameter C1=f(K,xi, DT) (Chow, 25-41)
548 1 : var = nc%getVariable("L11_C1")
549 1 : call var%getData(dummyD2)
550 1 : L11_C1(s11 : e11) = pack(dummyD2, mask11)
551 :
552 : ! Routing parameter C2 =f(K,xi, DT) (Chow, 25-41)
553 1 : var = nc%getVariable("L11_C2")
554 1 : call var%getData(dummyD2)
555 1 : L11_C2(s11 : e11) = pack(dummyD2, mask11)
556 :
557 : ! Fraction of the flood plain with impervious cover
558 1 : var = nc%getVariable("L11_nLinkFracFPimp")
559 1 : deallocate(dummyD3)
560 1 : call var%getData(dummyD3)
561 3 : do ii = 1, nLCoverScene
562 3 : L11_nLinkFracFPimp(s11 : e11, ii) = pack(dummyD3(:, :, ii), mask11)
563 : end do
564 :
565 : ! free memory
566 1 : deallocate(dummyD2, dummyD3)
567 :
568 1 : end subroutine mrm_read_restart_states
569 :
570 : ! ------------------------------------------------------------------
571 :
572 : ! NAME
573 : ! mrm_read_restart_config
574 :
575 : ! PURPOSE
576 : !> \brief reads Level 11 configuration from a restart directory
577 :
578 : !> \details read Level 11 configuration variables from a given restart
579 : !> directory and initializes all Level 11 configuration variables,
580 : !> that are initialized in L11_variable_init,
581 : !> contained in module mo_startup.
582 :
583 : ! INTENT(IN)
584 : !> \param[in] "integer(i4) :: iDomain" number of Domain
585 : !> \param[in] "character(256) :: InFile" Input Path including trailing slash
586 :
587 : ! HISTORY
588 : !> \authors Stephan Thober
589 :
590 : !> \date Apr 2013
591 :
592 : ! Modifications:
593 : ! Matthias Zink Apr 2014 - added inflow gauge
594 : ! Stephan Thober Aug 2015 - adapted for mRM
595 : ! Stephan Thober Sep 2015 - added all read restart commands in this subroutine
596 : ! Stephan Thober Sep 2015 - added L11_areaCell, L1_ID and L1_L11_Id for routing resolution higher than hydrology resolution
597 : ! David Schaefer Mar 2016 - mo_netcdf
598 : ! Stephan Thober Nov 2016 - added L11_TSrout, ProcessMatrix
599 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
600 : ! Stephan Thober May 2019 - added L0 info required for Process 3
601 :
602 1 : subroutine mrm_read_restart_config(iDomain, domainID, InFile)
603 :
604 1 : use mo_append, only : append
605 : use mo_common_constants, only : nodata_dp
606 : use mo_common_variables, only : level0, level1, domainMeta, processMatrix, domainMeta
607 : use mo_kind, only : dp, i4
608 : use mo_mpr_global_variables, only : L0_slope
609 : use mo_mrm_global_variables, only : L0_fdir, L0_fAcc, L0_streamnet, &
610 : L11_L1_Id, L11_TSrout, L11_aFloodPlain, L11_colOut, L11_fCol, &
611 : L11_fDir, L11_fAcc, L11_fRow, L11_fromN, L11_label, L11_length, L11_nOutlets, L11_netPerm, &
612 : L11_rOrder, L11_rowOut, L11_sink, L11_slope, L11_tCol, L11_tRow, L11_toN, &
613 : L1_L11_Id, domain_mrm, level11
614 : use mo_netcdf, only : NcDataset, NcVariable
615 : use mo_string_utils, only : num2str
616 :
617 : implicit none
618 :
619 : ! number of Domain
620 : integer(i4), intent(in) :: iDomain
621 :
622 : ! domain
623 : integer(i4), intent(in) :: domainID
624 :
625 : ! Input Path including trailing slash
626 : character(256), intent(in) :: InFile
627 :
628 : character(256) :: fname
629 :
630 : ! Mask at Level 0
631 1 : logical, allocatable, dimension(:, :) :: mask0
632 :
633 : ! Mask at Level 1
634 1 : logical, allocatable, dimension(:, :) :: mask1
635 :
636 : ! Mask at Level 11
637 1 : logical, allocatable, dimension(:, :) :: mask11
638 :
639 : ! dummy, 1 dimension I4
640 1 : integer(i4), allocatable, dimension(:) :: dummyI1
641 :
642 : ! dummy, 2 dimension I4
643 1 : integer(i4), allocatable, dimension(:, :) :: dummyI2
644 :
645 : ! dummy, DP
646 1 : real(dp), allocatable, dimension(:) :: dummyD1
647 1 : real(dp), allocatable, dimension(:, :) :: dummyD2
648 1 : real(dp), allocatable :: dummyD0
649 :
650 : type(NcDataset) :: nc
651 :
652 : type(NcVariable) :: var
653 :
654 :
655 : ! set file name
656 1 : fname = trim(InFile)
657 1 : call message(' Reading mRM restart file: ', trim(adjustl(Fname)), ' ...')
658 :
659 : ! get Domain info at L11 mask
660 124852 : mask0 = level0(domainMeta%L0DataFrom(iDomain))%mask
661 67 : mask1 = level1(iDomain)%mask
662 67 : mask11 = level11(iDomain)%mask
663 :
664 1 : nc = NcDataset(fname, "r")
665 :
666 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
667 : ! Read Process Matrix for check <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
668 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
669 1 : var = nc%getVariable("ProcessMatrix")
670 1 : allocate(dummyI1(size(processMatrix, dim = 1)))
671 1 : call var%getData(dummyI1)
672 1 : if (dummyI1(8) .ne. processMatrix(8, 1)) then
673 0 : call error_message('***ERROR: process description for routing', raise=.false.)
674 0 : call error_message('***ERROR: given in restart file does not match', raise=.false.)
675 0 : call error_message('***ERROR: that in namelist', raise=.false.)
676 0 : call error_message('***ERROR: restart file value:. ' // num2str(dummyI1(8), '(i2)'), raise=.false.)
677 0 : call error_message('***ERROR: namelist value:..... ' // num2str(processMatrix(8, 1), '(i2)'), raise=.false.)
678 0 : call error_message('ERROR: mrm_read_restart_config')
679 : end if
680 1 : deallocate(dummyI1)
681 :
682 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
683 : ! Read L0 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
684 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
685 1 : if (processMatrix(8, 1) .eq. 3_i4) then
686 : ! add L0_fdir, L0_fAcc, L0_slope, L0_streamnet
687 1 : var = nc%getVariable("L0_fDir")
688 1 : call var%getData(dummyI2)
689 1 : call append(L0_fdir, pack(dummyI2, mask0))
690 :
691 1 : var = nc%getVariable("L0_fAcc")
692 1 : call var%getData(dummyI2)
693 1 : call append(L0_fAcc, pack(dummyI2, mask0))
694 :
695 1 : var = nc%getVariable("L0_slope")
696 1 : call var%getData(dummyD2)
697 1 : call append(L0_Slope, pack(dummyD2, mask0))
698 :
699 1 : var = nc%getVariable("L0_streamnet")
700 1 : call var%getData(dummyI2)
701 1 : call append(L0_streamnet, pack(dummyI2, mask0))
702 : end if
703 :
704 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
705 : ! Read L1 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
706 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
707 1 : var = nc%getVariable("L1_L11_Id")
708 1 : call var%getData(dummyI2)
709 1 : call append(L1_L11_Id, pack(dummyI2, mask1))
710 :
711 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
712 : ! Read L11 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
713 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
714 : ! read L11_TSrout
715 1 : if (iDomain .eq. 1) then
716 3 : allocate(L11_TSrout(domainMeta%nDomains))
717 3 : L11_tsRout = nodata_dp
718 : end if
719 1 : var = nc%getVariable("L11_TSrout")
720 1 : call var%getData(dummyD0)
721 1 : L11_TSrout(iDomain) = dummyD0
722 :
723 : ! L11 data sets
724 : ! Mapping of L1 Ids on L1
725 1 : var = nc%getVariable("L11_L1_Id")
726 1 : call var%getData(dummyI2)
727 1 : call append(L11_L1_Id, pack(dummyI2, mask11))
728 :
729 : ! Flow direction (standard notation)
730 1 : var = nc%getVariable("L11_fDir")
731 1 : call var%getData(dummyI2)
732 1 : call append(L11_fDir, pack(dummyI2, mask11))
733 : ! append Number of Outlets at Level 11 (where facc == 0 )
734 64 : call append(L11_nOutlets, count((dummyI2 .eq. 0_i4)))
735 :
736 : ! Flow accumulation
737 1 : var = nc%getVariable("L11_fAcc")
738 1 : call var%getData(dummyD2)
739 1 : call append(L11_fAcc, pack(dummyD2, mask11))
740 :
741 : ! Grid vertical location of the Outlet
742 1 : var = nc%getVariable("L11_rowOut")
743 1 : call var%getData(dummyI2)
744 1 : call append(L11_rowOut, pack(dummyI2, mask11))
745 :
746 : ! Grid horizontal location of the Outlet
747 1 : var = nc%getVariable("L11_colOut")
748 1 : call var%getData(dummyI2)
749 1 : call append(L11_colOut, pack(dummyI2, mask11))
750 :
751 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
752 : ! From Node
753 1 : var = nc%getVariable("L11_fromN")
754 1 : call var%getData(dummyI1)
755 1 : call append(L11_fromN, dummyI1)
756 :
757 : ! To Node
758 1 : var = nc%getVariable("L11_toN")
759 1 : call var%getData(dummyI1)
760 1 : call append(L11_toN, dummyI1)
761 :
762 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
763 : ! Network routing order
764 1 : var = nc%getVariable("L11_rOrder")
765 1 : call var%getData(dummyI1)
766 1 : call append(L11_rOrder, dummyI1)
767 :
768 : ! Label Id [0='', 1=HeadWater, 2=Sink]
769 1 : var = nc%getVariable("L11_label")
770 1 : call var%getData(dummyI1)
771 1 : call append(L11_label, dummyI1)
772 :
773 : ! .true. if sink node reached
774 1 : var = nc%getVariable("L11_sink")
775 1 : call var%getData(dummyI1)
776 35 : call append(L11_sink, (dummyI1 .eq. 1_i4))
777 :
778 : ! Routing sequence (permutation of L11_rOrder)
779 1 : var = nc%getVariable("L11_netPerm")
780 1 : call var%getData(dummyI1)
781 1 : call append(L11_netPerm, dummyI1)
782 :
783 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
784 : ! From row in L0 grid
785 1 : var = nc%getVariable("L11_fRow")
786 1 : call var%getData(dummyI1)
787 1 : call append(L11_fRow, dummyI1)
788 :
789 : ! From col in L0 grid
790 1 : var = nc%getVariable("L11_fCol")
791 1 : call var%getData(dummyI1)
792 1 : call append(L11_fCol, dummyI1)
793 :
794 : ! To row in L0 grid
795 1 : var = nc%getVariable("L11_tRow")
796 1 : call var%getData(dummyI1)
797 1 : call append(L11_tRow, dummyI1)
798 :
799 : ! To col in L0 grid
800 1 : var = nc%getVariable("L11_tCol")
801 1 : call var%getData(dummyI1)
802 1 : call append(L11_tCol, dummyI1)
803 :
804 :
805 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
806 : ! read gaugenodelist
807 1 : var = nc%getVariable("gaugeNodeList")
808 1 : call var%getData(dummyI1)
809 2 : domain_mrm(iDomain)%gaugeNodeList(:) = dummyI1
810 :
811 : ! read InflowGaugeNodelist
812 1 : if (domain_mrm(iDomain)%nInflowGauges > 0) then
813 0 : var = nc%getVariable("InflowGaugeNodeList")
814 0 : call var%getData(dummyI1)
815 0 : domain_mrm(iDomain)%InflowgaugeNodeList(:) = dummyI1
816 : end if
817 :
818 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
819 : ! L11 network data sets
820 : ! [m] Total length of river link
821 1 : var = nc%getVariable("L11_length")
822 1 : call var%getData(dummyD1)
823 1 : call append(L11_length, dummyD1)
824 :
825 : ! [m2] Area of the flood plain
826 1 : var = nc%getVariable("L11_aFloodPlain")
827 1 : call var%getData(dummyD1)
828 1 : call append(L11_aFloodPlain, dummyD1)
829 :
830 : ! Average slope of river link
831 1 : var = nc%getVariable("L11_slope")
832 1 : call var%getData(dummyD1)
833 1 : call append(L11_slope, dummyD1)
834 :
835 1 : end subroutine mrm_read_restart_config
836 :
837 : end module mo_mrm_restart
|
