Line data Source code
1 : !> \file mo_mrm_mpr.f90
2 : !> \brief \copybrief mo_mrm_mpr
3 : !> \details \copydetails mo_mrm_mpr
4 :
5 : !> \brief Perform Multiscale Parameter Regionalization on Routing Parameters
6 : !> \details This module contains the subroutine for calculating the regionalized
7 : !! routing parameters (beta-parameters) given the five global routing parameters (gamma) at the level 0 scale.
8 : !> \authors Luis Samaniego, Stephan Thober
9 : !> \date Aug 2015
10 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
11 : !! mHM is released under the LGPLv3+ license \license_note
12 : !> \ingroup f_mrm
13 : module mo_mrm_mpr
14 : use mo_kind, only : dp
15 : use mo_message, only : message
16 :
17 : implicit none
18 :
19 : public :: reg_rout
20 : public :: mrm_init_param
21 : public :: mrm_update_param
22 : private
23 : contains
24 :
25 : ! ----------------------------------------------------------------------------
26 :
27 : ! NAME
28 : ! reg_rout
29 :
30 : ! PURPOSE
31 : !> \brief Regionalized routing
32 :
33 : !> \details sets up the Regionalized Routing parameters
34 : !> Global parameters needed (see mhm_parameter.nml):
35 : !> - param(1) = muskingumTravelTime_constant
36 : !> - param(2) = muskingumTravelTime_riverLength
37 : !> - param(3) = muskingumTravelTime_riverSlope
38 : !> - param(4) = muskingumTravelTime_impervious
39 : !> - param(5) = muskingumAttenuation_riverSlope
40 :
41 : ! INTENT(IN)
42 : !> \param[in] "real(dp), dimension(5) :: param" input parameter
43 : !> \param[in] "real(dp), dimension(:) :: length" [m] total length
44 : !> \param[in] "real(dp), dimension(:) :: slope" average slope
45 : !> \param[in] "real(dp), dimension(:) :: fFPimp" fraction of the flood plain with
46 : !> impervious layer
47 : !> \param[in] "real(dp) :: TS" - [h] time step in
48 :
49 : ! INTENT(OUT)
50 : !> \param[out] "real(dp), dimension(:) :: C1" routing parameter C1 (Chow, 25-41)
51 : !> \param[out] "real(dp), dimension(:) :: C2" routing parameter C2 (")
52 :
53 : ! HISTORY
54 : !> \authors Stephan Thober, Rohini Kumar
55 :
56 : !> \date Dec 2012
57 :
58 : ! Modifications:
59 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
60 :
61 990192 : subroutine reg_rout(param, length, slope, fFPimp, TS, C1, C2)
62 : implicit none
63 :
64 : ! input parameter
65 : real(dp), dimension(5), intent(in) :: param
66 :
67 : ! [m] total length
68 : real(dp), dimension(:), intent(in) :: length
69 :
70 : ! average slope
71 : real(dp), dimension(:), intent(in) :: slope
72 :
73 : ! fraction of the flood plain with
74 : ! impervious layer
75 : real(dp), dimension(:), intent(in) :: fFPimp
76 :
77 : ! - [h] time step in
78 : real(dp), intent(in) :: TS
79 :
80 : ! routing parameter C1 (Chow, 25-41)
81 : real(dp), dimension(:), intent(out) :: C1
82 :
83 : ! routing parameter C2 (")
84 : real(dp), dimension(:), intent(out) :: C2
85 :
86 : ! stream slope max
87 23151264 : real(dp) :: ssMax
88 :
89 : ! [d] Muskingum travel time parameter
90 23646360 : real(dp), dimension(size(fFPimp, 1)) :: K
91 :
92 : ! [1] Muskingum diffusion parameter (attenuation)
93 23646360 : real(dp), dimension(size(fFPimp, 1)) :: xi
94 :
95 :
96 : ! normalize stream bed slope
97 23646360 : ssMax = maxval(slope(:))
98 :
99 : ! New regional relationship; K = f(length, slope, & fFPimp)
100 495096 : K = param(1) + param(2) * (length * 0.001_dp) &
101 0 : + param(3) * slope &
102 23646360 : + param(4) * fFPimp
103 :
104 : ! Xi = f(slope)
105 23151264 : xi = param(5) * (1.0_dp + slope / ssMax)
106 :
107 : ! constraints on Xi
108 23151264 : xi = merge(0.5_dp, xi, xi > 0.5_dp)
109 23151264 : xi = merge(0.005_dp, xi, xi < 0.005_dp)
110 :
111 : ! constrains on Ki
112 23151264 : K = merge(0.5_dp * TS / xi, K, K > 0.5_dp * TS / xi)
113 23151264 : K = merge(0.5_dp * TS / (1.0_dp - xi), K, K < 0.5_dp * TS / (1.0_dp - xi))
114 :
115 : ! Muskingum parameters
116 23646360 : C1 = TS / (K * (1.0_dp - xi) + 0.5_dp * TS)
117 23646360 : C2 = 1.0_dp - C1 * K / TS
118 :
119 495096 : end subroutine reg_rout
120 :
121 : ! --------------------------------------------------------------------------
122 : ! L11 PARAMETERS
123 : ! --------------------------------------------------------------------------
124 : ! The parameters are set following Thober et al. 2017
125 : ! Modified:
126 : ! NAME
127 : ! mrm_init_param
128 :
129 : ! PURPOSE
130 : !> \brief TODO: add description
131 :
132 : !> \details TODO: add description
133 :
134 : ! INTENT(IN)
135 : !> \param[in] "integer(i4) :: iDomain" domain number
136 : !> \param[in] "real(dp), dimension(:) :: param" input parameter (param(1) is celerity in m/s)
137 :
138 : ! HISTORY
139 : !> \authors Stephan Thober, Matthias Kelbling
140 :
141 : !> \date Jun 2018
142 :
143 : ! Modifications:
144 :
145 25 : subroutine mrm_init_param(iDomain, param)
146 :
147 495096 : use mo_constants, only : HourSecs
148 : use mo_common_mHM_mRM_variables, only : resolutionRouting, timeStep, optimize
149 : use mo_common_variables, only : iFlag_cordinate_sys, domainMeta, processMatrix
150 : use mo_kind, only : dp, i4
151 : use mo_mrm_constants, only : given_TS
152 : use mo_mrm_global_variables, only : level11, L11_tsRout, domain_mrm, L11_celerity
153 : use mo_string_utils, only : num2str
154 : use mo_utils, only : locate, notequal
155 : use mo_mrm_net_startup, only : L11_calc_celerity
156 :
157 : implicit none
158 :
159 : ! domain number
160 : integer(i4), intent(in) :: iDomain
161 :
162 : ! input parameter (param(1) is celerity in m/s)
163 : real(dp), dimension(:), intent(in) :: param
164 :
165 : ! index selected from given_TS
166 : integer(i4) :: index
167 :
168 : ! spatial routing resolution
169 : real(dp) :: deltaX
170 :
171 : ! [s] wave travel time parameter
172 : real(dp) :: K
173 :
174 : ! start and end index at level11
175 : integer(i4) :: s11, e11
176 :
177 : ! Number of cells within
178 :
179 : ! initialize indices
180 25 : s11 = level11(iDomain)%iStart
181 25 : e11 = level11(iDomain)%iEnd
182 :
183 : ! temporal resolution of routing
184 25 : if (iDomain .eq. 1 .and. .not. allocated(L11_tsRout)) then
185 36 : allocate(L11_tsRout(domainMeta%nDomains))
186 36 : L11_TSrout = 0._dp
187 : end if
188 :
189 25 : if (processMatrix(8, 1) .eq. 1) then
190 68 : L11_tsRout = timestep * HourSecs
191 :
192 32 : if ( NOTEQUAL(mod(HourSecs * 24.0_dp, L11_tsRout(iDomain)), 0.0_dp) .and. &
193 16 : (domain_mrm(iDomain)%nInflowGauges .gt. 0)) then
194 0 : call message('***WARNING: routing timestep is not a multiple of 24 h.')
195 0 : call message(' Inflowgauge timeseries is averaged over values')
196 0 : call message(' of different days, small mismatches at')
197 0 : call message(' inflowgauge location might occur.')
198 : end if
199 :
200 : else
201 :
202 : ! called for initialization
203 9 : call mrm_update_param(iDomain, param)
204 :
205 : end if
206 :
207 25 : call message('')
208 25 : call message(' domain: '//num2str(iDomain, '(i3)'))
209 25 : call message(' routing resolution [s]:. '//num2str(L11_tsRout(iDomain), '(f7.0)'))
210 25 : call message(' routing factor:......... '//num2str(L11_tsRout(iDomain) / (timestep * HourSecs), '(f5.2)'))
211 :
212 50 : if ( NOTEQUAL(mod(HourSecs * 24.0_dp, L11_tsRout(iDomain)), 0.0_dp) .and. &
213 25 : (domain_mrm(iDomain)%nInflowGauges .gt. 0)) then
214 0 : call message('***WARNING: routing timestep is not a multiple of 24 h.')
215 0 : call message(' Inflowgauge timeseries is averaged over values')
216 0 : call message(' of different days, small mismatches at')
217 0 : call message(' inflowgauge location might occur.')
218 : end if
219 :
220 25 : end subroutine mrm_init_param
221 :
222 : ! NAME
223 : ! mrm_update_param
224 :
225 : ! PURPOSE
226 : !> \brief TODO: add description
227 :
228 : !> \details TODO: add description
229 :
230 : ! INTENT(IN)
231 : !> \param[in] "integer(i4) :: iDomain" domain number
232 : !> \param[in] "real(dp), dimension(1) :: param" celerity parameter [m s-1]
233 :
234 : ! HISTORY
235 : !> \authors Stehpan Thober, Matthias Kelbling
236 :
237 : !> \date Jun 2018
238 :
239 : ! Modifications:
240 :
241 19 : subroutine mrm_update_param(iDomain, param)
242 :
243 25 : use mo_kind, only: i4, dp
244 : use mo_common_variables, only: processMatrix, iFlag_cordinate_sys
245 : use mo_mrm_global_variables, only: &
246 : ! input variable
247 : level11, &
248 : L11_TSrout, &
249 : L11_celerity, L11_nOutlets, L11_length, &
250 : ! output variables
251 : L11_C1, L11_C2
252 : use mo_common_mHM_mRM_variables, only: resolutionRouting, optimize, timeStep, &
253 : optimize
254 : use mo_mrm_constants, only: rout_space_weight, given_TS
255 : use mo_utils, only: locate
256 : use mo_mrm_net_startup, only: L11_calc_celerity
257 : use mo_mrm_constants, only: given_TS
258 : use mo_constants, only: HourSecs
259 : use mo_string_utils, only: num2str
260 : use mo_utils, only: locate
261 :
262 : implicit none
263 :
264 : ! domain number
265 : integer(i4), intent(in) :: iDomain
266 :
267 : ! celerity parameter [m s-1]
268 : real(dp), intent(in), dimension(1) :: param
269 :
270 : integer(i4) :: s11
271 : integer(i4) :: e11
272 : integer(i4) :: nNodes
273 :
274 : ! index selected from given_TS
275 : integer(i4) :: ind
276 :
277 : ! spatial routing resolution
278 : real(dp) :: deltaX
279 : real(dp), allocatable :: length(:)
280 :
281 : ! [s] wave travel time parameter
282 19 : real(dp), allocatable :: K(:)
283 :
284 : ! [1] Muskingum diffusion parameter (attenuation)
285 19 : real(dp) :: xi
286 :
287 : ! get domain information
288 19 : s11 = level11(iDomain)%iStart
289 19 : e11 = level11(iDomain)%iEnd
290 19 : Nnodes = level11(iDomain)%nCells
291 :
292 57 : allocate(K(nNodes))
293 :
294 19 : if (ProcessMatrix(8, 1) .eq. 2_i4) then
295 :
296 : ! [s] wave travel time parameter
297 70 : K(:) = L11_length(s11: e11) / param(1)
298 :
299 17 : else if (ProcessMatrix(8, 1) .eq. 3_i4) then
300 :
301 : ! [s] wave travel time parameter
302 17 : call L11_calc_celerity( iDomain, param)
303 :
304 : ! Allocate and calculate K
305 595 : K(:) = L11_length(s11: e11) / L11_celerity(s11:e11)
306 :
307 : end if
308 :
309 : ! set time-weighting scheme
310 19 : xi = abs(rout_space_weight) ! set weighting factor to 0._dp
311 :
312 : ! determine routing timestep
313 684 : ind = locate(given_TS, minval(K(1:(nNodes-L11_nOutlets(iDomain)))))
314 :
315 : ! set min-wave traveltime to min given_TS
316 19 : if (ind .lt. 1) ind = 1
317 19 : L11_TSrout(iDomain) = given_TS(ind)
318 :
319 : ! Muskingum parameters
320 665 : L11_C1(s11:e11) = L11_TSrout(iDomain) / ( K(:) * (1.0_dp - xi) + 0.5_dp * L11_TSrout(iDomain) )
321 665 : L11_C2(s11:e11) = 1.0_dp - L11_C1(s11:e11) * K(:) / L11_TSrout(iDomain)
322 :
323 19 : deallocate(K)
324 :
325 : ! optional print
326 : ! print *, 'C1 Muskingum routing parameter: ', L11_C1(s11)
327 : ! print *, 'C2 Muskingum routing parameter: ', L11_C2(s11)
328 :
329 19 : end subroutine mrm_update_param
330 :
331 : end module mo_mrm_mpr
|