Line data Source code
1 : !> \file mo_mrm_pre_routing.f90
2 : !> \brief \copybrief mo_mrm_pre_routing
3 : !> \details \copydetails mo_mrm_pre_routing
4 :
5 : !> \brief Performs pre-processing for routing for mHM at level L11.
6 : !> \details This module performs runoff accumulation from L1 to L11 and inflow summation.
7 : !> \changelog
8 : !! - Stephan Thober Aug 2015
9 : !! - adapted to mRM
10 : !! - Sebastian Mueller Jun 2020
11 : !! - separate module for pre-processing
12 : !> \authors Luis Samaniego
13 : !> \date Dec 2012
14 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
15 : !! mHM is released under the LGPLv3+ license \license_note
16 : !> \ingroup f_mrm
17 : MODULE mo_mrm_pre_routing
18 :
19 : ! This module performs pre-processing for routing for mHM at level L11.
20 :
21 : ! Written Sebastian Mueller Jun 2020
22 :
23 : USE mo_kind, ONLY : i4, dp
24 :
25 : IMPLICIT NONE
26 :
27 : PRIVATE
28 :
29 : PUBLIC :: L11_runoff_acc
30 : PUBLIC :: add_inflow
31 : PUBLIC :: L11_meteo_acc
32 : PUBLIC :: calc_L1_runoff_E
33 :
34 : ! ------------------------------------------------------------------
35 :
36 : CONTAINS
37 :
38 : ! ------------------------------------------------------------------
39 :
40 : ! NAME
41 : ! L11_runoff_acc
42 :
43 : ! PURPOSE
44 : !> \brief total runoff accumulation at L11.
45 :
46 : !> \details Upscales runoff in space from L1 to L11 if routing resolution
47 : !> is higher than hydrology resolution (map_flag equals .true.) or
48 : !> downscales runoff from L1 to L11 if routing resolution is lower
49 : !> than hydrology resolution.
50 :
51 : ! INTENT(IN)
52 : !> \param[in] "real(dp), dimension(:) :: qall" total runoff L1 [mm TS-1]
53 : !> \param[in] "real(dp), dimension(:) :: efecarea" effective area in [km2] at Level 1
54 : !> \param[in] "integer(i4), dimension(:) :: L1_L11_Id" L11 Ids mapped on L1
55 : !> \param[in] "real(dp), dimension(:) :: L11_areacell" effective area in [km2] at Level 11
56 : !> \param[in] "integer(i4), dimension(:) :: L11_L1_Id" L1 Ids mapped on L11
57 : !> \param[in] "integer(i4) :: TS" time step in [s]
58 : !> \param[in] "logical :: map_flag" Flag indicating whether routing resolution is higher than
59 : !> hydrologic one
60 :
61 : ! INTENT(OUT)
62 : !> \param[out] "real(dp), dimension(:) :: qAcc" aggregated runoff at L11 [m3 s-1]
63 :
64 : ! HISTORY
65 : !> \authors Luis Samaniego
66 :
67 : !> \date Jan 2013
68 :
69 : ! Modifications:
70 : ! Matthias Zink Mar 2014 - added inflow from upstream areas
71 : ! Matthias Zink Dec 2014 - adopted inflow gauges to ignore headwater cells
72 : ! Stephan Thober Sep 2015 - included downscaling of runoff
73 : ! Stephan Thober Feb 2016 - refactored upscaling of discharge from L1 to L11
74 : ! Stephan Thober Feb 2016 - refactored downscaling of discharge from L1 to L11
75 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
76 :
77 1281576 : SUBROUTINE L11_runoff_acc(qAll, efecArea, L1_L11_Id, L11_areaCell, L11_L1_Id, TS, map_flag, qAcc)
78 :
79 : use mo_constants, only : HourSecs
80 : use mo_common_constants, only : nodata_dp
81 :
82 : implicit none
83 :
84 : ! total runoff L1 [mm TS-1]
85 : real(dp), intent(in), dimension(:) :: qall
86 : ! effective area in [km2] at Level 1
87 : real(dp), intent(in), dimension(:) :: efecarea
88 : ! L11 Ids mapped on L1
89 : integer(i4), intent(in), dimension(:) :: L1_L11_Id
90 : ! effective area in [km2] at Level 11
91 : real(dp), intent(in), dimension(:) :: L11_areacell
92 : ! L1 Ids mapped on L11
93 : integer(i4), intent(in), dimension(:) :: L11_L1_Id
94 : ! time step in [h]
95 : integer(i4), intent(in) :: TS
96 : ! Flag indicating whether routing resolution is higher than hydrologic one
97 : logical, intent(in) :: map_flag
98 : ! aggregated runoff at L11 [m3 s-1]
99 : real(dp), intent(out), dimension(:) :: qAcc
100 :
101 : integer(i4) :: k
102 : ! [s] time step
103 640788 : real(dp) :: TST
104 :
105 :
106 : ! ------------------------------------------------------------------
107 : ! ACCUMULATION OF DISCHARGE TO A ROUTING CELL
108 : ! ------------------------------------------------------------------
109 : ! Hydrologic timestep in seconds
110 640788 : TST = HourSecs * TS
111 :
112 640788 : if (map_flag) then
113 : ! Estimate specific runoff at L11
114 : ! NOTE:
115 : ! 1) Total discharge depth aggregated at L11 level [mm/TST]
116 : ! 2) Transform depth [mm/TST] to discharge [m3/s]
117 : ! Total runoff should be divided by total_area to get
118 : ! specific discharge at L11. Then, to transform specific
119 : ! discharge from [mm/TST] to [m3/s], it should be multiplied by
120 : ! total_area [km2]*10^3 and divided by TST.
121 : ! Therefore, in this operation total_area cancels out.
122 28745580 : qAcc = 0._dp
123 : ! loop over high-resolution cells (L1) and add discharge to
124 : ! corresponding low-resolution cells (L11)
125 30385380 : do k = 1, size(qAll, 1)
126 30385380 : qAcc(L1_L11_Id(k)) = qAcc(L1_L11_Id(k)) + qAll(k) * efecArea(k)
127 : end do
128 : ! factor 1000 for conversion of mm*km2 in m3
129 28745580 : qAcc = qAcc * 1000.0_dp / TST
130 : !
131 : else
132 : ! initialize qout
133 0 : qAcc = nodata_dp
134 0 : do k = 1, size(qAcc, 1)
135 : ! map flux from coarse L1 resolution to fine L11 resolution
136 0 : qAcc(k) = qAll(L11_L1_Id(k))
137 : end do
138 : ! adjust flux by area cell
139 : ! factor 1000 for conversion of mm*km2 in m3
140 0 : qAcc(:) = qAcc(:) * L11_areaCell(:) * 1000.0_dp / TST
141 : end if
142 :
143 640788 : END SUBROUTINE L11_runoff_acc
144 :
145 : ! ------------------------------------------------------------------
146 :
147 : ! NAME
148 : ! add_inflow
149 :
150 : ! PURPOSE
151 : !> \brief
152 : !> Adds inflow discharge to the runoff produced at the
153 : !> cell where the inflow is occurring.
154 : !> \details
155 : !> If a inflow gauge is given, then this routine is adding the
156 : !> values to the runoff produced at the grid cell where the
157 : !> inflow is happening. The values are not directly added to the
158 : !> river network. If this cell is not a headwater then the streamflow
159 : !> produced upstream will be neglected.
160 :
161 : ! INTENT(IN)
162 : !> \param[in] "integer(i4) :: nInflowGauges" [-] number of inflow points
163 : !> \param[in] "integer(i4), dimension(:) :: InflowIndexList" [-] index of inflow points
164 : !> \param[in] "logical, dimension(:) :: InflowHeadwater" [-] if to consider headwater cells of inflow gauge
165 : !> \param[in] "integer(i4), dimension(:) :: InflowNodeList" [-] L11 ID of inflow points
166 : !> \param[in] "real(dp), dimension(:) :: QInflow" [m3 s-1] inflowing water
167 :
168 : ! INTENT(INOUT)
169 : !> \param[inout] "real(dp), dimension(:) :: qOut" [m3 s-1] Series of attenuated runoff
170 :
171 : ! HISTORY
172 : !> \authors Stephan Thober & Matthias Zink
173 :
174 : !> \date Jul 2016
175 :
176 : ! Modifications:
177 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
178 :
179 1268472 : subroutine add_inflow(nInflowGauges, InflowIndexList, InflowHeadwater, InflowNodeList, QInflow, qOut)
180 :
181 640788 : use mo_kind, only : dp, i4
182 :
183 : implicit none
184 :
185 : ! [-] number of inflow points
186 : integer(i4), intent(in) :: nInflowGauges
187 : ! [-] index of inflow points
188 : integer(i4), intent(in), dimension(:) :: InflowIndexList
189 : ! [-] if to consider headwater cells of inflow gauge
190 : logical, intent(in), dimension(:) :: InflowHeadwater
191 : ! [-] L11 ID of inflow points
192 : integer(i4), intent(in), dimension(:) :: InflowNodeList
193 : ! [m3 s-1] inflowing water
194 : real(dp), intent(in), dimension(:) :: QInflow
195 : ! [m3 s-1] Series of attenuated runoff
196 : real(dp), intent(inout), dimension(:) :: qOut
197 :
198 : integer(i4) :: ii
199 :
200 :
201 : ! discharge for inflow gauges (e.g. for missing upstream catchments) is added here
202 : ! should be put after UH attenuation because it is measured runoff at this cell
203 634236 : if (nInflowGauges .gt. 0) then
204 52416 : do ii = 1, nInflowGauges
205 52416 : if (InflowHeadwater(ii)) then
206 : ! add inflowing water to water produced by upstream/headwater cells
207 13104 : qOut(InflowNodeList(ii)) = qOut(InflowNodeList(ii)) + QInflow(InflowIndexList(ii))
208 : else
209 : ! put only timeseries and cut upstream/headwater cells produced water for routing
210 13104 : qOut(InflowNodeList(ii)) = QInflow(InflowIndexList(ii))
211 : end if
212 : end do
213 : end if
214 634236 : end subroutine add_inflow
215 :
216 : ! ------------------------------------------------------------------
217 :
218 : ! NAME
219 : ! L11_E_acc
220 :
221 : ! PURPOSE
222 : !> \brief temperature energy accumulation at L11.
223 :
224 : !> \details Upscales energy in space from L1 to L11 if routing resolution
225 : !> is higher than hydrology resolution (map_flag equals .true.) or
226 : !> downscales runoff from L1 to L11 if routing resolution is lower
227 : !> than hydrology resolution.
228 :
229 : ! INTENT(IN)
230 : !> \param[in] "real(dp), dimension(:) :: qall" total runoff L1 [mm K TS-1]
231 : !> \param[in] "real(dp), dimension(:) :: efecarea" effective area in [km2] at Level 1
232 : !> \param[in] "integer(i4), dimension(:) :: L1_L11_Id" L11 Ids mapped on L1
233 : !> \param[in] "real(dp), dimension(:) :: L11_areacell" effective area in [km2] at Level 11
234 : !> \param[in] "integer(i4), dimension(:) :: L11_L1_Id" L1 Ids mapped on L11
235 : !> \param[in] "integer(i4) :: TS" time step in [s]
236 : !> \param[in] "logical :: map_flag" Flag indicating whether routing resolution is higher than
237 : !> hydrologic one
238 :
239 : ! INTENT(OUT)
240 : !> \param[out] "real(dp), dimension(:) :: qAcc" aggregated runoff at L11 [m3 K s-1]
241 :
242 : ! HISTORY
243 : !> \authors Luis Samaniego
244 :
245 : !> \date Jan 2013
246 :
247 : ! Modifications:
248 : ! Matthias Zink Mar 2014 - added inflow from upstream areas
249 : ! Matthias Zink Dec 2014 - adopted inflow gauges to ignore headwater cells
250 : ! Stephan Thober Sep 2015 - included downscaling of runoff
251 : ! Stephan Thober Feb 2016 - refactored upscaling of discharge from L1 to L11
252 : ! Stephan Thober Feb 2016 - refactored downscaling of discharge from L1 to L11
253 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
254 :
255 : SUBROUTINE L11_E_acc(qAll, efecArea, L1_L11_Id, L11_areaCell, L11_L1_Id, TS, map_flag, qAcc)
256 :
257 634236 : use mo_constants, only : HourSecs
258 : use mo_common_constants, only : nodata_dp
259 :
260 : implicit none
261 :
262 : ! total runoff L1 [mm TS-1]
263 : real(dp), intent(in), dimension(:) :: qall
264 : ! effective area in [km2] at Level 1
265 : real(dp), intent(in), dimension(:) :: efecarea
266 : ! L11 Ids mapped on L1
267 : integer(i4), intent(in), dimension(:) :: L1_L11_Id
268 : ! effective area in [km2] at Level 11
269 : real(dp), intent(in), dimension(:) :: L11_areacell
270 : ! L1 Ids mapped on L11
271 : integer(i4), intent(in), dimension(:) :: L11_L1_Id
272 : ! time step in [s]
273 : integer(i4), intent(in) :: TS
274 : ! Flag indicating whether routing resolution is higher than hydrologic one
275 : logical, intent(in) :: map_flag
276 : ! aggregated runoff at L11 [m3 s-1]
277 : real(dp), intent(out), dimension(:) :: qAcc
278 :
279 : integer(i4) :: k
280 : ! [s] time step
281 : real(dp) :: TST
282 :
283 :
284 : ! ------------------------------------------------------------------
285 : ! ACCUMULATION OF DISCHARGE TO A ROUTING CELL
286 : ! ------------------------------------------------------------------
287 : ! Hydrologic timestep in seconds
288 : TST = HourSecs * TS
289 :
290 : if (map_flag) then
291 : ! Estimate specific runoff at L11
292 : ! NOTE:
293 : ! 1) Total discharge depth aggregated at L11 level [mm/TST]
294 : ! 2) Transform depth [mm/TST] to discharge [m3/s]
295 : ! Total runoff should be divided by total_area to get
296 : ! specific discharge at L11. Then, to transform specific
297 : ! discharge from [mm/TST] to [m3/s], it should be multiplied by
298 : ! total_area [km2]*10^3 and divided by TST.
299 : ! Therefore, in this operation total_area cancels out.
300 : qAcc = 0._dp
301 : ! loop over high-resolution cells (L1) and add discharge to
302 : ! corresponding low-resolution cells (L11)
303 : do k = 1, size(qAll, 1)
304 : qAcc(L1_L11_Id(k)) = qAcc(L1_L11_Id(k)) + qAll(k) * efecArea(k)
305 : end do
306 : qAcc = qAcc * 1000.0_dp / TST
307 : !
308 : else
309 : ! initialize qout
310 : qAcc = nodata_dp
311 : do k = 1, size(qAcc, 1)
312 : ! map flux from coarse L1 resolution to fine L11 resolution
313 : qAcc(k) = qAll(L11_L1_Id(k))
314 : end do
315 : ! adjust flux by area cell
316 : qAcc(:) = qAcc(:) * L11_areaCell(:) * 1000.0_dp / TST
317 : end if
318 :
319 : END SUBROUTINE L11_E_acc
320 :
321 : ! ------------------------------------------------------------------
322 :
323 : ! NAME
324 : ! calc_L1_runoff_E
325 :
326 : ! PURPOSE
327 : !> \brief calculate lateral temperature energy from runoff components.
328 :
329 : !> \details calculate lateral temperature energy from runoff components.
330 :
331 : ! INTENT(IN)
332 : !> \param[in] "REAL(dp) :: fSealed_area_fraction" sealed area fraction [1]
333 : !> \param[in] "REAL(dp) :: fast_interflow" \f$ q_0 \f$ Fast runoff component [mm TS-1]
334 : !> \param[in] "REAL(dp) :: slow_interflow" \f$ q_1 \f$ Slow runoff component [mm TS-1]
335 : !> \param[in] "REAL(dp) :: baseflow" \f$ q_2 \f$ Baseflow [mm TS-1]
336 : !> \param[in] "REAL(dp) :: direct_runoff" \f$ q_D \f$ Direct runoff from impervious areas [mm TS-1]
337 : !> \param[in] "REAL(dp) :: temp_air" air temperature [K]
338 : !> \param[in] "REAL(dp) :: mean_temp_air" annual mean air temperature [K]
339 :
340 : ! INTENT(OUT)
341 : !> \param[out] "REAL(dp) :: lateral_E" \f$ E_T \f$ Generated runoff [K mm TS-1]
342 :
343 : ! HISTORY
344 : !> \authors Sebastian Mueller
345 :
346 : !> \date Jun 2020
347 :
348 13104 : SUBROUTINE calc_L1_runoff_E( &
349 26208 : fSealed_area_fraction, &
350 26208 : fast_interflow, &
351 13104 : slow_interflow, &
352 26208 : baseflow, &
353 13104 : direct_runoff, &
354 13104 : temp_air, &
355 13104 : mean_temp_air, &
356 13104 : lateral_E &
357 : )
358 :
359 : use mo_constants, only : T0_dp ! 273.15 - Celcius <-> Kelvin [K]
360 :
361 : implicit none
362 :
363 : ! sealed area fraction [1]
364 : REAL(dp), dimension(:), INTENT(IN) :: fSealed_area_fraction
365 : ! \f$ q_0 \f$ Fast runoff component [mm TS-1]
366 : REAL(dp), dimension(:), INTENT(IN) :: fast_interflow
367 : ! \f$ q_1 \f$ Slow runoff component [mm TS-1]
368 : REAL(dp), dimension(:), INTENT(IN) :: slow_interflow
369 : ! \f$ q_2 \f$ Baseflow [mm TS-1]
370 : REAL(dp), dimension(:), INTENT(IN) :: baseflow
371 : ! \f$ q_D \f$ Direct runoff from impervious areas [mm TS-1]
372 : REAL(dp), dimension(:), INTENT(IN) :: direct_runoff
373 : ! air temperature [degC]
374 : real(dp), dimension(:), intent(in) :: temp_air
375 : ! annual mean air temperature [degC]
376 : real(dp), dimension(:), intent(in) :: mean_temp_air
377 : ! \f$ E_T \f$ Generated lateral Energy [K mm TS-1]
378 : REAL(dp), dimension(:), INTENT(inout) :: lateral_E
379 :
380 : ! convert temperatures form [deg C] to [K]
381 : ! accumulate in [K mm TS-1] -> convert later to [K m3 s-1] on L11
382 : ! following Wanders et.al. 2019
383 : lateral_E = lateral_E + ( &
384 0 : (baseflow * max(T0_dp + 5.0_dp, mean_temp_air + T0_dp) &
385 0 : + (slow_interflow + fast_interflow) * max(T0_dp, temp_air + T0_dp)) &
386 0 : * (1.0_dp - fSealed_area_fraction) &
387 0 : + direct_runoff * max(T0_dp, temp_air + T0_dp - 1.5_dp) * fSealed_area_fraction &
388 471744 : )
389 :
390 13104 : END SUBROUTINE calc_L1_runoff_E
391 :
392 : ! ------------------------------------------------------------------
393 :
394 : ! NAME
395 : ! L11_meteo_acc
396 :
397 : ! PURPOSE
398 : !> \brief meteo forcing accumulation at L11 for temperature routing.
399 :
400 : !> \details Upscales meteo forcing in space from L1 to L11 if routing resolution
401 : !> is higher than hydrology resolution (map_flag equals .true.) or
402 : !> downscales meteo forcing from L1 to L11 if routing resolution is lower
403 : !> than hydrology resolution.
404 :
405 : ! INTENT(IN)
406 : !> \param[in] "real(dp), dimension(:) :: meteo_all" meteo forcing
407 : !> \param[in] "real(dp), dimension(:) :: efecarea" effective area in [km2] at Level 1
408 : !> \param[in] "integer(i4), dimension(:) :: L1_L11_Id" L11 Ids mapped on L1
409 : !> \param[in] "real(dp), dimension(:) :: L11_areacell" effective area in [km2] at Level 11
410 : !> \param[in] "integer(i4), dimension(:) :: L11_L1_Id" L1 Ids mapped on L11
411 : !> \param[in] "logical :: map_flag" Flag indicating whether routing resolution is higher than
412 : !> hydrologic one
413 :
414 : ! INTENT(OUT)
415 : !> \param[out] "real(dp), dimension(:) :: meteo_acc" aggregated meteo forcing
416 :
417 : ! HISTORY
418 : !> \authors Sebastian Mueller
419 :
420 : !> \date Jul 2020
421 :
422 19657 : SUBROUTINE L11_meteo_acc(meteo_all, efecArea, L1_L11_Id, L11_areaCell, L11_L1_Id, map_flag, meteo_acc)
423 :
424 13104 : use mo_common_constants, only : nodata_dp
425 :
426 : implicit none
427 :
428 : ! meteo forcing (state variable)
429 : real(dp), intent(in), dimension(:) :: meteo_all
430 : ! effective area in [km2] at Level 1
431 : real(dp), intent(in), dimension(:) :: efecarea
432 : ! L11 Ids mapped on L1
433 : integer(i4), intent(in), dimension(:) :: L1_L11_Id
434 : ! effective area in [km2] at Level 11
435 : real(dp), intent(in), dimension(:) :: L11_areacell
436 : ! L1 Ids mapped on L11
437 : integer(i4), intent(in), dimension(:) :: L11_L1_Id
438 : ! Flag indicating whether routing resolution is higher than hydrologic one
439 : logical, intent(in) :: map_flag
440 : ! aggregated meteo forcing
441 : real(dp), intent(out), dimension(:) :: meteo_acc
442 :
443 : integer(i4) :: k
444 :
445 19657 : if (map_flag) then
446 687995 : meteo_acc = 0._dp
447 : ! loop over high-resolution cells (L1) and add meteo forcing (weighted by area) to
448 : ! corresponding low-resolution cells (L11)
449 687995 : do k = 1, size(meteo_all, 1)
450 687995 : meteo_acc(L1_L11_Id(k)) = meteo_acc(L1_L11_Id(k)) + meteo_all(k) * efecArea(k)
451 : end do
452 : ! divide by L11 cell area to get weighted mean
453 707652 : meteo_acc = meteo_acc / L11_areaCell(:)
454 : else
455 : ! initialize
456 0 : meteo_acc = nodata_dp
457 0 : do k = 1, size(meteo_acc, 1)
458 : ! map from coarse L1 resolution to fine L11 resolution
459 0 : meteo_acc(k) = meteo_all(L11_L1_Id(k))
460 : end do
461 : end if
462 :
463 19657 : END SUBROUTINE L11_meteo_acc
464 :
465 : END MODULE mo_mrm_pre_routing
|
