Line data Source code
1 : !> \file mo_mrm_signatures.f90
2 : !> \brief \copybrief mo_mrm_signatures
3 : !> \details \copydetails mo_mrm_signatures
4 :
5 : !> \brief Module with calculations for several hydrological signatures.
6 : !> \details This module contains calculations for hydrological signatures.
7 : !!
8 : !! It contains:
9 : !! - Autocorrelation
10 : !! - Rising and declining limb densities
11 : !! - Flow duration curves
12 : !! - Peak distribution
13 : !> \authors Remko Nijzink,
14 : !> \date March 2014
15 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
16 : !! mHM is released under the LGPLv3+ license \license_note
17 : !> \ingroup f_mrm
18 : MODULE mo_mrm_signatures
19 :
20 : USE mo_kind, ONLY : i4, sp, dp
21 : use mo_message, only : message, error_message
22 :
23 : IMPLICIT NONE
24 :
25 : PUBLIC :: Autocorrelation ! Autocorrelation function
26 : PUBLIC :: FlowDurationCurve ! Flow duration curve (i.e. CDF of runoff)
27 : PUBLIC :: Limb_densities ! Rising and declining limb densities
28 : PUBLIC :: Moments ! Moments of data and log-transformed data, e.g. mean and standard deviation.
29 : PUBLIC :: PeakDistribution ! Peak distribution parameter
30 : PUBLIC :: RunoffRatio ! Runoff ratio (accumulated daily discharge [mm/d] / accumulated daily precipitation [mm/d])
31 : PUBLIC :: ZeroFlowRatio ! Ratio of zero flow days to total observation days
32 :
33 : ! ------------------------------------------------------------------
34 :
35 : CONTAINS
36 :
37 : !-------------------------------------------------------------------------------
38 : ! NAME
39 : ! Autocorrelation
40 :
41 : ! PURPOSE
42 : !> \brief Autocorrelation of a given data series.
43 :
44 : !> \details Calculates the autocorrelation of a data series at given lags.
45 : !> An optional argument for masking data points can be given.
46 : !> The function is basically a wrapper of the function autocorr
47 : !> from the module mo_corr.
48 : !> An optional mask of data points can be specified.
49 :
50 : !> ADDITIONAL INFORMATION
51 : !> Used as hydrologic signature with lag 1 in
52 : !> Euser, T., Winsemius, H. C., Hrachowitz, M., Fenicia, F., Uhlenbrook, S., & Savenije, H. H. G. (2013).
53 : !> A framework to assess the realism of model structures using hydrological signatures.
54 : !> Hydrology and Earth System Sciences, 17(5), 1893-1912. doi:10.5194/hess-17-1893-2013
55 :
56 : ! INTENT(IN)
57 : !> \param[in] "real(dp), dimension(:) :: data" Array of data
58 : !> \param[in] "integer(i4), dimension(:) :: lags" Array of lags where autocorrelation is requested
59 :
60 : ! INTENT(IN), OPTIONAL
61 : !> \param[in] "logical, dimension(size(data, 1)), optional :: mask" Mask for data points givenWorks only with 1d
62 : !> double precision input data.
63 :
64 : ! HISTORY
65 : !> \authors Juliane Mai
66 :
67 : !> \date Jun 2015
68 :
69 : ! Modifications:
70 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
71 :
72 0 : FUNCTION Autocorrelation(data, lags, mask)
73 :
74 : use mo_corr, only : autocorr
75 :
76 : implicit none
77 :
78 : ! Array of data
79 : real(dp), dimension(:), intent(in) :: data
80 :
81 : ! Array of lags where autocorrelation is requested
82 : integer(i4), dimension(:), intent(in) :: lags
83 :
84 : ! Mask for data points givenWorks only with 1d double precision input data.
85 : logical, dimension(size(data, 1)), optional, intent(in) :: mask
86 :
87 : ! Autocorrelation of data at given lags
88 : real(dp), dimension(size(lags, 1)) :: autocorrelation
89 :
90 :
91 0 : if (present(mask)) then
92 0 : autocorrelation = autocorr(data, lags, mask = mask)
93 : else
94 0 : autocorrelation = autocorr(data, lags)
95 : end if
96 :
97 0 : END FUNCTION Autocorrelation
98 :
99 : !-------------------------------------------------------------------------------
100 : ! NAME
101 : ! FlowDurationCurve
102 :
103 : ! PURPOSE
104 : !> \brief Flow duration curves.
105 :
106 : !> \details Calculates the flow duration curves for a given data vector. The Flow duration curve at a
107 : !> certain quantile x is the data point p where x% of the data points are above the value p.
108 : !> Hence the function percentile of the module mo_percentile is used. But percentile is
109 : !> determining the point p where x% of the data points are below that value. Therfore, the
110 : !> given quantiles are transformed by (1.0-quantile) to get the percentiles of exceedance probabilities.
111 :
112 : !> Optionally, the concavity index CI can be calculated [Zhang2014]. CI is defined by
113 : !> \f[ CI = \frac{q_{10\%}-q_{99\%}}{q_{1\%}-q_{99\%}} \f]
114 : !> where \f$ q_{x} \f$ is the data point where x% of the data points are above that value.
115 : !> Hence, exceedance probabilities are used.
116 :
117 : !> Optionally, the FDC mid-segment slope \f$FDC_{MSS}\f$ as used by Shafii et. al (2014) can be returned.
118 : !> The \f$FDC_{MSS}\f$ is defined as
119 : !> \f[ FDC_{MSS} = \log(q_{m_1})-\log(q_{m_2}) \f]
120 : !> where \f$ m_1 \f$ and \f$ m_2 \f$ are the lowest and highest flow exceedance probabilities within the
121 : !> midsegment of FDC. The settings \f$m_1=0.2\f$ and \f$0.7\f$ are used by Shafii et. al (2014) and are
122 : !> implemented like that.
123 :
124 : !> Optionally, the FDC medium high-segment volume \f$FDC_{MHSV}\f$ as used by Shafii et. al (2014) can be
125 : !> returned. The \f$FDC_{MHSV}\f$ is defined as
126 : !> \f[ FDC_{MHSV} = \sum_{h=1}^{H} q_h \f]
127 : !> where \f$ h=1,2,...,H \f$ are flow indeces located within the high-flow segment (exceedance probabilities
128 : !> lower than \f$m_1\f$). \f$H\f$ is the index of the maximum flow. The settings \f$m_1=0.2\f$ is used here
129 : !> to be consistent with the definitions of the low-segment (0.7-1.0) and the mid-segment (0.2-0.7).
130 :
131 : !> Optionally, the FDC high-segment volume \f$FDC_{HSV}\f$ as used by Shafii et. al (2014) can be returned.
132 : !> The \f$FDC_{HSV}\f$ is defined as
133 : !> \f[ FDC_{HSV} = \sum_{h=1}^{H} q_h \f]
134 : !> where \f$ h=1,2,...,H \f$ are flow indeces located within the high-flow segment (exceedance probabilities
135 : !> lower than \f$m_1\f$). \f$H\f$ is the index of the maximum flow. The settings \f$m_1=0.02\f$ is used by
136 : !> Shafii et. al (2014) and is implemented like that.
137 :
138 : !> Optionally, the FDC low-segment volume \f$FDC_{LSV}\f$ as used by Shafii et. al (2014) can be returned.
139 : !> The \f$FDC_{LSV}\f$ is defined as
140 : !> \f[ FDC_{LSV} = -\sum_{l=1}^{L} (\log(q_l) - \log(q_L)) \f]
141 : !> where \f$ l=1,2,...,L \f$ are flow indeces located within the low-flow segment (exceedance probabilities
142 : !> larger than \f$m_1\f$). \f$L\f$ is the index of the minimum flow. The settings \f$m_1=0.7\f$ is used by
143 : !> Shafii et. al (2014) and is implemented like that.
144 :
145 : !> An optional mask of data points can be specified.
146 : !> ADDITIONAL INFORMATION
147 :
148 : !> Thresholds in mid_segment_slope, mhigh_segment_volume, high_segment_volume, low_segment_volume are hard
149 : !> coded.
150 : !> FDC is used as hydrologic signature (quantiles not specified) in
151 : !> Euser, T., Winsemius, H. C., Hrachowitz, M., Fenicia, F., Uhlenbrook, S., & Savenije, H. H. G. (2013).
152 : !> A framework to assess the realism of model structures using hydrological signatures.
153 : !> Hydrology and Earth System Sciences, 17(5), 1893-1912. doi:10.5194/hess-17-1893-2013
154 : !> Concavity Index used as hydrologic signature in
155 : !> Zhang, Y., Vaze, J., Chiew, F. H. S., Teng, J., & Li, M. (2014).
156 : !> Predicting hydrological signatures in ungauged catchments using spatial interpolation, index model, and
157 : !> rainfall-runoff modelling.
158 : !> Journal of Hydrology, 517(C), 936-948. doi:10.1016/j.jhydrol.2014.06.032
159 : !> Concavity index is defined using exceedance probabilities by
160 : !> Sauquet, E., & Catalogne, C. (2011).
161 : !> Comparison of catchment grouping methods for flow duration curve estimation at ungauged sites in France.
162 : !> Hydrology and Earth System Sciences, 15(8), 2421-2435. doi:10.5194/hess-15-2421-2011
163 : !> mid_segment_slope, high_segment_volume, low_segment_volume used as hydrologic signature in
164 : !> Shafii, M., & Tolson, B. A. (2015).
165 : !> Optimizing hydrological consistency by incorporating hydrological signatures into model calibration
166 : !> objectives.
167 : !> Water Resources Research, 51(5), 3796-3814. doi:10.1002/2014WR016520
168 :
169 : ! INTENT(IN)
170 : !> \param[in] "real(dp), dimension(:) :: data" data series
171 : !> \param[in] "real(dp), dimension(:) :: quantiles" Percentages of exceedance (x-axis of FDC)
172 :
173 : ! INTENT(IN), OPTIONAL
174 : !> \param[in] "logical, dimension(:), optional :: mask" mask of data array
175 :
176 : ! INTENT(OUT), OPTIONAL
177 : !> \param[out] "real(dp), optional :: concavity_index" concavity index as defined by Sauquet et al. (2011)
178 : !> \param[out] "real(dp), optional :: mid_segment_slope" mid-segment slope as defined by Shafii et al. (2014)
179 : !> \param[out] "real(dp), optional :: mhigh_segment_volume" medium high-segment volume
180 : !> \param[out] "real(dp), optional :: high_segment_volume" high-segment volume as defined by Shafii et al.
181 : !> (2014)
182 : !> \param[out] "real(dp), optional :: low_segment_volume" low-segment volume as defined by Shafii et al.
183 : !> (2014)
184 :
185 : ! RETURN
186 : !> \return real(dp), dimension(size(quantiles,1)) :: FlowDurationCurve — Flow Duration Curve value at
187 : !> resp. quantile
188 :
189 : ! HISTORY
190 : !> \authors Remko Nijzink, Juliane Mai
191 :
192 : !> \date March 2014
193 :
194 : ! Modifications:
195 : ! Juliane Mai Jun 2015 - mask added
196 : ! - function instead of subroutine
197 : ! - use of percentile
198 : ! - add concavity_index
199 : ! Juliane Mai Jun 2015 - add mid_segment_slope, mhigh_segment_volume, high_segment_volume, low_segment_volume
200 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
201 :
202 0 : FUNCTION FlowDurationCurve(data, quantiles, mask, concavity_index, mid_segment_slope, mhigh_segment_volume, &
203 : high_segment_volume, low_segment_volume)
204 :
205 0 : use mo_percentile, only : percentile
206 : use mo_utils, only : ge, le
207 :
208 : implicit none
209 :
210 : ! data series
211 : real(dp), dimension(:), intent(in) :: data
212 :
213 : ! Percentages of exceedance (x-axis of FDC)
214 : real(dp), dimension(:), intent(in) :: quantiles
215 :
216 : ! mask of data array
217 : logical, dimension(:), optional, intent(in) :: mask
218 :
219 : ! concavity index as defined by Sauquet et al. (2011)
220 : real(dp), optional, intent(out) :: concavity_index
221 :
222 : ! mid-segment slope as defined by Shafii et al. (2014)
223 : real(dp), optional, intent(out) :: mid_segment_slope
224 :
225 : ! medium high-segment volume
226 : real(dp), optional, intent(out) :: mhigh_segment_volume
227 :
228 : ! high-segment volume as defined by Shafii et al. (2014)
229 : real(dp), optional, intent(out) :: high_segment_volume
230 :
231 : ! low-segment volume as defined by Shafii et al. (2014)
232 : real(dp), optional, intent(out) :: low_segment_volume
233 :
234 : ! data point where x% of the data points
235 : ! are above that value
236 : real(dp), dimension(size(quantiles, 1)) :: FlowDurationCurve
237 :
238 : ! mask for data points
239 0 : logical, dimension(size(data, 1)) :: maske
240 :
241 : ! minimal flow value
242 0 : real(dp) :: min_flow_value
243 :
244 : ! flow value at a threshold quantile
245 0 : real(dp) :: flow_value_thres
246 :
247 :
248 : ! checking optionals
249 0 : if (present(mask)) then
250 0 : maske = mask
251 : else
252 0 : maske = .true.
253 : end if
254 :
255 0 : FlowDurationCurve = percentile(data, (1._dp - quantiles) * 100._dp, mask = maske, mode_in = 5)
256 :
257 0 : if (present(concavity_index)) then
258 : concavity_index = &
259 : (percentile(data, (1._dp - 0.10_dp) * 100._dp, mask = maske, mode_in = 5) - &
260 : percentile(data, (1._dp - 0.99_dp) * 100._dp, mask = maske, mode_in = 5)) / &
261 : (percentile(data, (1._dp - 0.01_dp) * 100._dp, mask = maske, mode_in = 5) - &
262 0 : percentile(data, (1._dp - 0.99_dp) * 100._dp, mask = maske, mode_in = 5))
263 : end if
264 :
265 0 : if (present(mid_segment_slope)) then
266 : ! mid-flows are defined to be between 0.2 and 0.7 by Shafii et. al (2014)
267 : mid_segment_slope = &
268 : log(percentile(data, (1._dp - 0.2_dp) * 100._dp, mask = maske, mode_in = 5)) - &
269 0 : log(percentile(data, (1._dp - 0.7_dp) * 100._dp, mask = maske, mode_in = 5))
270 : end if
271 :
272 0 : if (present(mhigh_segment_volume)) then
273 : ! medium high-flows are defined to be between 0.0 and 0.2 as to be constistent
274 : ! with the mid-segment (0.2-0.7) and low-segment (0.7-1.0) definitions
275 0 : flow_value_thres = percentile(data, (1._dp - 0.2_dp) * 100._dp, mask = maske, mode_in = 5)
276 0 : mhigh_segment_volume = sum(data, mask = (maske .and. ge(data, flow_value_thres)))
277 : ! print*, 'flow_value_thres = ',flow_value_thres
278 : ! print*, 'mhigh_segment_volume = ',mhigh_segment_volume
279 : end if
280 :
281 0 : if (present(high_segment_volume)) then
282 : ! high-flows are defined to be between 0.0 and 0.02 by Shafii et. al (2014)
283 0 : flow_value_thres = percentile(data, (1._dp - 0.02_dp) * 100._dp, mask = maske, mode_in = 5)
284 0 : high_segment_volume = sum(data, mask = (maske .and. ge(data, flow_value_thres)))
285 : end if
286 :
287 0 : if (present(low_segment_volume)) then
288 : ! low-flows are defined to be between 0.7 and 1.0 by Shafii et. al (2014)
289 0 : min_flow_value = minval(data, mask = maske)
290 0 : flow_value_thres = percentile(data, (1._dp - 0.7) * 100._dp, mask = maske, mode_in = 5)
291 : low_segment_volume = -1.0_dp * &
292 0 : sum(log(data) - log(min_flow_value), mask = (maske .and. le(data, flow_value_thres)))
293 : end if
294 :
295 0 : END FUNCTION FlowDurationCurve
296 :
297 : !-------------------------------------------------------------------------------
298 : ! NAME
299 : ! Limb_densities
300 :
301 : ! PURPOSE
302 : !> \brief Calculates limb densities
303 :
304 : !> \details Calculates rising and declinging limb densities. The peaks of the given series are
305 : !> first determined by looking for points where preceding and subsequent datapoint are lower.
306 : !> Second, the number of datapoints with rising values (nrise) and declining values (ndecline)
307 : !> are counted basically by comparing neighbors.
308 : !> The duration the data increase (nrise) divided by the number of peaks (npeaks)
309 : !> gives the rising limb density RLD
310 : !> \f[ RLD=t_{rise}/n_{peak} \f]
311 : !> whereas the duration the data decrease (ndecline) divided by the number of peaks (npeaks)
312 : !> gives the declining limb density DLD
313 : !> \f[ DLD=t_{fall}/n_{peak}. \f]
314 : !> An optional mask of data points can be specified.
315 :
316 : !> ADDITIONAL INFORMATION
317 : !> Rising limb density used as hydrologic signature in
318 : !> Euser, T., Winsemius, H. C., Hrachowitz, M., Fenicia, F., Uhlenbrook, S., & Savenije, H. H. G. (2013).
319 : !> A framework to assess the realism of model structures using hydrological signatures.
320 : !> Hydrology and Earth System Sciences, 17(5), 1893-1912. doi:10.5194/hess-17-1893-2013
321 :
322 : ! INTENT(IN)
323 : !> \param[in] "real(dp), dimension(:) :: data" data series
324 :
325 : ! INTENT(IN), OPTIONAL
326 : !> \param[in] "logical, dimension(size(data, 1)), optional :: mask" mask for data series
327 :
328 : ! INTENT(OUT), OPTIONAL
329 : !> \param[out] "real(dp), optional :: RLD" rising limb density
330 : !> \param[out] "real(dp), optional :: DLD" declining limb density
331 :
332 : ! HISTORY
333 : !> \authors Remko Nijzink
334 :
335 : !> \date March 2014
336 :
337 : ! Modifications:
338 : ! Juliane Mai Jun 2015 - RLD and DLD as optional
339 : ! - optional mask for data can be given
340 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
341 :
342 0 : SUBROUTINE Limb_densities(data, mask, RLD, DLD)
343 :
344 : implicit none
345 :
346 : ! data series
347 : real(dp), dimension(:), intent(in) :: data
348 :
349 : ! mask for data series
350 : logical, dimension(size(data, 1)), optional, intent(in) :: mask
351 :
352 : ! rising limb density
353 : real(dp), optional, intent(out) :: RLD
354 :
355 : ! declining limb density
356 : real(dp), optional, intent(out) :: DLD
357 :
358 : ! mask for data points
359 0 : logical, dimension(size(data, 1)) :: maske
360 :
361 : ! Counter
362 : integer(i4) :: jj
363 :
364 : ! Number of peaks
365 : integer(i4) :: n_peak
366 :
367 : ! Number of declining data points
368 : integer(i4) :: n_decline
369 :
370 : ! Number of rising data points
371 : integer(i4) :: n_rise
372 :
373 : ! True if rising, False if declining or contains masked values
374 0 : logical, dimension(size(data, 1)) :: goes_up
375 :
376 : ! Threshold that is has to rise at least to be detected as rising value
377 0 : real(dp) :: thres_rise
378 :
379 :
380 : ! checking optionals
381 0 : if (present(mask)) then
382 0 : maske = mask
383 : else
384 0 : maske = .true.
385 : end if
386 :
387 0 : if ((.not. present(RLD)) .and. (.not. present(DLD))) then
388 0 : call error_message('mo_signatures: limb_densities: Neither RLD or DLD is specified in calling sequence.')
389 : end if
390 :
391 : ! initialize
392 0 : n_rise = 0_i4
393 0 : n_decline = 0_i4
394 0 : n_peak = 0_i4
395 :
396 0 : goes_up = .False.
397 0 : thres_rise = 1.0_dp
398 0 : do jj = 1, size(data, 1) - 1
399 0 : if (maske(jj) .and. maske(jj + 1)) then
400 0 : if (data(jj) < data(jj + 1) - thres_rise) then
401 0 : goes_up(jj) = .true.
402 : ! print*, jj, ' ', data(jj), ' ', data(jj+1)
403 : end if
404 : end if
405 : end do
406 0 : n_rise = count(goes_up)
407 0 : n_decline = count(maske) - count(goes_up)
408 :
409 : ! write(*,*) 'goes_up = ', goes_up(1:178)
410 :
411 : ! peak is where goes_up switches from true to false
412 0 : n_peak = 0_i4
413 0 : do jj = 1, size(data, 1) - 1
414 0 : if (maske(jj) .and. maske(jj + 1)) then
415 0 : if (goes_up(jj) .and. .not.(goes_up(jj + 1))) then
416 0 : n_peak = n_peak + 1_i4
417 : ! print*, jj
418 : end if
419 : end if
420 : end do
421 :
422 : ! do jj=2, size(data,1)-1
423 :
424 : ! ! check for peak
425 : ! if (maske(jj-1) .and. maske(jj) .and. maske(jj+1)) then
426 : ! if ( (data(jj-1) .lt. data(jj)-1.0_dp) .and. (data(jj+1) .lt. data(jj)-1.0_dp) ) then
427 : ! n_peak = n_peak+1_i4
428 : ! write(*,*) jj-1
429 : ! end if
430 : ! end if
431 :
432 : ! ! check if data has rised
433 : ! if (maske(jj-1) .and. maske(jj)) then
434 : ! if (data(jj-1) .lt. data(jj)-1.0_dp) then
435 : ! n_rise = n_rise+1_i4
436 : ! else
437 : ! n_decline = n_decline+1_i4
438 : ! end if
439 : ! end if
440 :
441 : ! ! ! check if data has declined
442 : ! ! if (maske(jj-1) .and. maske(jj)) then
443 : ! ! if (data(jj-1) .gt. data(jj)) then
444 : ! ! n_decline = n_decline+1_i4
445 : ! ! end if
446 : ! ! end if
447 :
448 : ! end do
449 :
450 : ! write(*,*) 'n_peak = ', n_peak
451 :
452 0 : if (present(RLD)) then
453 0 : if (n_peak>0_i4) then
454 0 : RLD = real(n_rise, dp) / real(n_peak, dp)
455 : else
456 0 : RLD = 0.0_dp
457 : end if
458 : end if
459 :
460 0 : if (present(DLD)) then
461 0 : if (n_peak>0_i4) then
462 0 : DLD = real(n_decline, dp) / real(n_peak, dp)
463 : else
464 0 : DLD = 0.0_dp
465 : end if
466 : end if
467 :
468 0 : END SUBROUTINE Limb_densities
469 :
470 : !-------------------------------------------------------------------------------
471 : ! NAME
472 : ! MaximumMonthlyFlow
473 :
474 : ! PURPOSE
475 : !> \brief Maximum of average flows per months.
476 :
477 : !> \details Maximum of average flow per month is defined as
478 : !> \f[ max_{monthly flow} = Max( F(i), i=1,..12 ) \f]
479 : !> where \$f F(i) $\f is the average flow of month i.
480 : !> ADDITIONAL INFORMATION
481 : !> used as hydrologic signature in
482 : !> Shafii, M., & Tolson, B. A. (2015).
483 : !> Optimizing hydrological consistency by incorporating hydrological signatures into model calibration
484 : !> objectives.
485 : !> Water Resources Research, 51(5), 3796-3814. doi:10.1002/2014WR016520
486 :
487 : ! INTENT(IN)
488 : !> \param[in] "real(dp), dimension(:) :: data" array of data
489 :
490 : ! INTENT(IN), OPTIONAL
491 : !> \param[in] "logical, dimension(size(data, 1)), optional :: mask" mask for data points given
492 : !> \param[in] "integer(i4), optional :: yr_start" year of date of first data point given
493 : !> \param[in] "integer(i4), optional :: mo_start" month of date of first data point given
494 : !> (default: 1)
495 : !> \param[in] "integer(i4), optional :: dy_start" month of date of first data point given
496 : !> (default: 1)
497 :
498 : ! RETURN
499 : !> \return real(dp) :: MaximumMonthlyFlow — Maximum of average flow per month
500 : !> Works only with 1d double precision input data.
501 : !> Assumes data are daily values.
502 :
503 : ! HISTORY
504 : !> \authors Juliane Mai
505 :
506 : !> \date Jun 2015
507 :
508 : ! Modifications:
509 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
510 :
511 0 : FUNCTION MaximumMonthlyFlow(data, mask, yr_start, mo_start, dy_start)
512 :
513 0 : use mo_julian, only : date2dec, dec2date
514 :
515 : implicit none
516 :
517 : ! array of data
518 : real(dp), dimension(:), intent(in) :: data
519 :
520 : ! mask for data points given
521 : logical, dimension(size(data, 1)), optional, intent(in) :: mask
522 :
523 : ! year of date of first data point given
524 : integer(i4), optional, intent(in) :: yr_start
525 :
526 : ! month of date of first data point given (default: 1)
527 : integer(i4), optional, intent(in) :: mo_start
528 :
529 : ! month of date of first data point given (default: 1)
530 : integer(i4), optional, intent(in) :: dy_start
531 :
532 : ! return: maximum of average monthly flow
533 : real(dp) :: MaximumMonthlyFlow
534 :
535 0 : logical, dimension(size(data, 1)) :: maske
536 :
537 : ! counter
538 : integer(i4) :: ii
539 :
540 : ! date variables
541 : integer(i4) :: yr, mo, dy, imo
542 :
543 : ! number of data points per month
544 : integer(i4), dimension(12) :: counter
545 :
546 : ! summed data points per months
547 0 : real(dp), dimension(12) :: flow_month
548 :
549 : ! julian day of one day before start day
550 0 : real(dp) :: ref_jul_day
551 :
552 :
553 0 : if (present(mask)) then
554 0 : maske = mask
555 : else
556 0 : maske = .true.
557 : end if
558 :
559 0 : if (.not. present(yr_start)) then
560 0 : call error_message('mo_signatures: MaximumMonthlyFlow: Year of of data point has to be given!')
561 : else
562 0 : yr = yr_start
563 : end if
564 :
565 0 : if (present(mo_start)) then
566 0 : mo = mo_start
567 : else
568 0 : mo = 1
569 : end if
570 :
571 0 : if (present(dy_start)) then
572 0 : dy = dy_start
573 : else
574 0 : dy = 1
575 : end if
576 :
577 0 : flow_month = 0.0_dp
578 0 : counter = 0_i4
579 0 : ref_jul_day = date2dec(yy = yr, mm = mo, dd = dy) - 1.0_dp
580 :
581 0 : do ii = 1, size(data, 1)
582 0 : if (maske(ii)) then
583 : ! determine current month
584 0 : call dec2date(ref_jul_day + real(ii, dp), mm = imo)
585 : ! add value
586 0 : counter(imo) = counter(imo) + 1
587 0 : flow_month(imo) = flow_month(imo) + data(ii)
588 : end if
589 : end do
590 :
591 0 : if (any(counter == 0_i4)) then
592 0 : call error_message('mo_signatures: MaximumMonthlyFlow: There are months with no data points!')
593 : end if
594 :
595 : ! average
596 0 : MaximumMonthlyFlow = maxval(flow_month / real(counter, dp))
597 :
598 0 : END FUNCTION MaximumMonthlyFlow
599 :
600 : !-------------------------------------------------------------------------------
601 : ! NAME
602 : ! Moments
603 :
604 : ! PURPOSE
605 : !> \brief Moments of data and log-transformed data, e.g. mean and standard deviation.
606 :
607 : !> \details Returns several moments of data series given, i.e.
608 : !> * mean of data
609 : !> * standard deviation of data
610 : !> * median of data
611 : !> * maximum/ peak of data
612 : !> * mean of log-transformed data
613 : !> * standard deviation of log-transformed data
614 : !> * median of log-transformed data
615 : !> * maximum/ peak of log-transformed data
616 : !> An optional mask of data points can be specified.
617 : !> ADDITIONAL INFORMATION
618 : !> mean_log and stddev_log used as hydrologic signature in
619 : !> Zhang, Y., Vaze, J., Chiew, F. H. S., Teng, J., & Li, M. (2014).
620 : !> Predicting hydrological signatures in ungauged catchments using spatial interpolation, index model, and
621 : !> rainfall-runoff modelling.
622 : !> Journal of Hydrology, 517(C), 936-948. doi:10.1016/j.jhydrol.2014.06.032
623 : !> mean_data, stddev_data, median_data, max_data, mean_log, and stddev_log used as hydrologic signature in
624 : !> Shafii, M., & Tolson, B. A. (2015).
625 : !> Optimizing hydrological consistency by incorporating hydrological signatures into model calibration
626 : !> objectives.
627 : !> Water Resources Research, 51(5), 3796-3814. doi:10.1002/2014WR016520
628 :
629 : ! INTENT(IN)
630 : !> \param[in] "real(dp), dimension(:) :: data" array of data
631 :
632 : ! INTENT(IN), OPTIONAL
633 : !> \param[in] "logical, dimension(size(data, 1)), optional :: mask" mask for data points given
634 :
635 : ! INTENT(OUT), OPTIONAL
636 : !> \param[out] "real(dp), optional :: mean_data" mean of data
637 : !> \param[out] "real(dp), optional :: stddev_data" standard deviation of data
638 : !> \param[out] "real(dp), optional :: median_data" median of data
639 : !> \param[out] "real(dp), optional :: max_data" maximum/ peak of data
640 : !> \param[out] "real(dp), optional :: mean_log" mean of log-transformed data
641 : !> \param[out] "real(dp), optional :: stddev_log" standard deviation of log-transformed data
642 : !> \param[out] "real(dp), optional :: median_log" median of log-transformed data
643 : !> \param[out] "real(dp), optional :: max_log" maximum/ peak of log-transformed dataWorks only with 1d
644 : !> double precision input data.
645 :
646 : ! HISTORY
647 : !> \authors Juliane Mai
648 :
649 : !> \date Jun 2015
650 :
651 : ! Modifications:
652 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
653 :
654 0 : SUBROUTINE Moments(data, mask, mean_data, stddev_data, median_data, max_data, mean_log, stddev_log, median_log, &
655 : max_log)
656 :
657 0 : use mo_moment, only : mean, stddev
658 : use mo_percentile, only : median
659 :
660 : implicit none
661 :
662 : ! array of data
663 : real(dp), dimension(:), intent(in) :: data
664 :
665 : ! mask for data points given
666 : logical, dimension(size(data, 1)), optional, intent(in) :: mask
667 :
668 : ! mean of data
669 : real(dp), optional, intent(out) :: mean_data
670 :
671 : ! standard deviation of data
672 : real(dp), optional, intent(out) :: stddev_data
673 :
674 : ! median of data
675 : real(dp), optional, intent(out) :: median_data
676 :
677 : ! maximum/ peak of data
678 : real(dp), optional, intent(out) :: max_data
679 :
680 : ! mean of log-transformed data
681 : real(dp), optional, intent(out) :: mean_log
682 :
683 : ! standard deviation of log-transformed data
684 : real(dp), optional, intent(out) :: stddev_log
685 :
686 : ! median of log-transformed data
687 : real(dp), optional, intent(out) :: median_log
688 :
689 : ! maximum/ peak of log-transformed dataWorks only with 1d double precision input data.
690 : real(dp), optional, intent(out) :: max_log
691 :
692 0 : logical, dimension(size(data, 1)) :: maske
693 :
694 0 : real(dp), dimension(size(data, 1)) :: logdata
695 :
696 :
697 0 : if (present(mask)) then
698 0 : maske = mask
699 : else
700 0 : maske = .true.
701 : end if
702 :
703 : if (.not.(present(mean_data)) .and. .not.(present(stddev_data)) .and. &
704 : .not.(present(median_data)) .and. .not.(present(max_data)) .and. &
705 : .not.(present(mean_log)) .and. .not.(present(stddev_log)) .and. &
706 0 : .not.(present(median_log)) .and. .not.(present(max_log))) then
707 0 : call error_message('mo_signatures: Moments: None of the optional output arguments is specified')
708 : end if
709 :
710 0 : if (present(mean_data)) mean_data = mean(data, mask = maske)
711 0 : if (present(stddev_data)) stddev_data = stddev(data, mask = maske)
712 0 : if (present(median_data)) median_data = median(data, mask = maske)
713 0 : if (present(max_data)) max_data = maxval(data, mask = maske)
714 :
715 0 : if (present(mean_log) .or. present(stddev_log)) then
716 0 : where (data > 0.0_dp)
717 0 : logdata = log(data)
718 : elsewhere
719 : logdata = -9999.0_dp ! will not be used, since mask is set to .false.
720 : maske = .false.
721 : end where
722 :
723 0 : if (present(mean_log)) mean_log = mean(logdata, mask = maske)
724 0 : if (present(stddev_log)) stddev_log = stddev(logdata, mask = maske)
725 0 : if (present(median_log)) median_log = median(logdata, mask = maske)
726 0 : if (present(max_log)) max_log = maxval(logdata, mask = maske)
727 : end if
728 :
729 0 : END SUBROUTINE Moments
730 :
731 : !-------------------------------------------------------------------------------
732 : ! NAME
733 : ! PeakDistribution
734 :
735 : ! PURPOSE
736 : !> \brief Calculates the peak distribution.
737 :
738 : !> \details First, the peaks of the time series given are identified. For the peak distribution
739 : !> only this subset of data points are considered. Second, the peak distribution at the
740 : !> quantiles given is calculated. Calculates the peak distribution at the quantiles given
741 : !> using mo_percentile. Since the exceedance probabilities are usually used in
742 : !> hydrology the function percentile is used with (1.0-quantiles).
743 :
744 : !> Optionally, the slope of the peak distribution between 10th and 50th percentile, i.e.
745 : !> \f[ slope = \frac{\mathrm{peak\_{data}}_{0.1}-\mathrm{peak\_{data}}_{0.5}}{0.9-0.5} \f]
746 : !> can be returned.
747 : !> An optional mask for the data points can be given.
748 : !> ADDITIONAL INFORMATION
749 : !> slope_peak_distribution used as hydrologic signature in
750 : !> Euser, T., Winsemius, H. C., Hrachowitz, M., Fenicia, F., Uhlenbrook, S., & Savenije, H. H. G. (2013).
751 : !> A framework to assess the realism of model structures using hydrological signatures.
752 : !> Hydrology and Earth System Sciences, 17(5), 1893-1912. doi:10.5194/hess-17-1893-2013
753 :
754 : ! INTENT(IN)
755 : !> \param[in] "real(dp), dimension(:) :: data" data array
756 : !> \param[in] "real(dp), dimension(:) :: quantiles" requested quantiles for distribution
757 :
758 : ! INTENT(IN), OPTIONAL
759 : !> \param[in] "logical, dimension(size(data, 1)), optional :: mask" mask of data array
760 :
761 : ! INTENT(OUT), OPTIONAL
762 : !> \param[out] "real(dp), optional :: slope_peak_distribution" slope of the Peak distribution between10th and
763 : !> 50th percentile
764 :
765 : ! RETURN
766 : !> \return real(dp), dimension(size(quantiles,1)) :: PeakDistribution — Distribution of peak values
767 : !> at resp. quantiles
768 :
769 : ! HISTORY
770 : !> \authors Remko Nijzink
771 :
772 : !> \date March 2014
773 :
774 : ! Modifications:
775 : ! Juliane Mai Jun 2015 - mask added
776 : ! - function instead of subroutine
777 : ! - use of percentile
778 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
779 :
780 0 : FUNCTION PeakDistribution(data, quantiles, mask, slope_peak_distribution)
781 :
782 0 : use mo_percentile, only : percentile
783 :
784 : implicit none
785 :
786 : ! data array
787 : real(dp), dimension(:), intent(in) :: data
788 :
789 : ! requested quantiles for distribution
790 : real(dp), dimension(:), intent(in) :: quantiles
791 :
792 : ! mask of data array
793 : logical, dimension(size(data, 1)), optional, intent(in) :: mask
794 :
795 : ! slope of the Peak distribution between10th and 50th percentile
796 : real(dp), optional, intent(out) :: slope_peak_distribution
797 :
798 : ! distribution of peaks in data
799 : ! returns values of distribution at
800 : ! given quantiles
801 : real(dp), dimension(size(quantiles, 1)) :: PeakDistribution
802 :
803 : ! counters
804 : integer(i4) :: ii, jj
805 :
806 : ! mask of data
807 0 : logical, dimension(size(data, 1)) :: maske
808 :
809 : ! Number of peaks
810 : integer(i4) :: n_peak
811 :
812 : ! array containing some quantiles
813 0 : real(dp), dimension(2) :: pp
814 :
815 : ! data points of quantiles pp
816 0 : real(dp), dimension(2) :: data_pp
817 :
818 : ! series containing only peak values of original series data
819 0 : real(dp), allocatable, dimension(:) :: data_peak
820 :
821 :
822 : ! checking optionals
823 0 : if (present(mask)) then
824 0 : maske = mask
825 : else
826 0 : maske = .true.
827 : end if
828 :
829 : ! count peaks
830 0 : n_peak = 0_i4
831 0 : do jj = 2, size(data, 1) - 1
832 0 : if (maske(jj - 1) .and. maske(jj) .and. maske(jj + 1)) then
833 0 : if ((data(jj - 1) .le. data(jj)) .and. (data(jj + 1) .le. data(jj))) then
834 0 : n_peak = n_peak + 1_i4
835 : end if
836 : end if
837 : end do
838 :
839 0 : allocate(data_peak(n_peak))
840 :
841 : ! find peaks
842 0 : jj = 0
843 0 : do ii = 2, size(data, 1) - 1
844 0 : if (maske(ii - 1) .and. maske(ii) .and. maske(ii + 1)) then
845 0 : if((data(ii - 1) .le. data(ii)) .and. (data(ii + 1) .le. data(ii))) then
846 0 : jj = jj + 1_i4
847 0 : data_peak(jj) = data(ii)
848 : end if
849 : end if
850 : end do
851 :
852 0 : if (present(slope_peak_distribution)) then
853 : ! calculate slope between 10% and 50% quantiles, per definition
854 0 : pp = (/ 0.1_dp, 0.5_dp /)
855 0 : data_pp = percentile(data_peak, (1.0_dp - pp) * 100._dp, mode_in = 5) ! (1-p) because exceedence probability is required
856 0 : slope_peak_distribution = (data_pp(1) - data_pp(2)) / (0.9_dp - 0.5_dp)
857 : end if
858 :
859 0 : PeakDistribution = percentile(data_peak, (1.0_dp - Quantiles) * 100._dp, mode_in = 5)
860 0 : deallocate(data_peak)
861 :
862 0 : END FUNCTION PeakDistribution
863 :
864 : !-------------------------------------------------------------------------------
865 : ! NAME
866 : ! RunoffRatio
867 :
868 : ! PURPOSE
869 : !> \brief Runoff ratio (accumulated daily discharge [mm/d] / accumulated daily precipitation [mm/d]).
870 :
871 : !> \details The runoff ratio is defined as
872 : !> \f[ runoff_ratio = \frac{\sum_{t=1}^{N} q_t}{\sum_{t=1}^{N} p_t}\f]
873 : !> where \f$p_t\f$ and \f$q_t\f$ are precipitation and discharge, respectively.
874 : !> Therefore, precipitation over the entire domain is required and both discharge and precipitation
875 : !> have to be converted to the same units [mm/d].
876 :
877 : !> Input discharge is given in [m**3/s] as this is mHM default while precipitation has to be given
878 : !> in [mm/km**2 / day].
879 :
880 : !> Either "precip_sum" or "precip_series" has to be specified.
881 : !> If "precip_series" is used the optional mask is also applied to precipitation values.
882 : !> The "precip_sum" is the accumulated "precip_series".
883 :
884 : !> Optionally, a mask for the data (=discharge) can be given. If optional "log_data" is set to .true.
885 : !> the runoff ratio will be calculated as
886 : !> \f[ runoff\_ratio = \frac{\sum_{t=1}^{N} \log(q_t)}{\sum_{t=1}^{N} p_t}\f]
887 : !> where \f$p_t\f$ and \f$q_t\f$ are precipitation and discharge, respectively.
888 : !> ADDITIONAL INFORMATION
889 : !> \return real(dp), dimension(size(lags,1)) :: RunoffRation — Ratio of discharge and precipitation
890 : !> Used as hydrologic signature in
891 : !> Shafii, M., & Tolson, B. A. (2015).
892 : !> Optimizing hydrological consistency by incorporating hydrological signatures into model calibration
893 : !> objectives.
894 : !> Water Resources Research, 51(5), 3796-3814. doi:10.1002/2014WR016520
895 :
896 : ! INTENT(IN)
897 : !> \param[in] "real(dp), dimension(:) :: data" array of data [m**3/s]
898 : !> \param[in] "real(dp) :: domain_area" area of domain [km**2]
899 :
900 : ! INTENT(IN), OPTIONAL
901 : !> \param[in] "logical, dimension(size(data, 1)), optional :: mask" mask for data points given
902 : !> \param[in] "real(dp), dimension(size(data, 1)), optional :: precip_series" daily precipitation values
903 : !> [mm/km**2 / day]
904 : !> \param[in] "real(dp), optional :: precip_sum" sum of daily precip. values of
905 : !> whole period[mm/km**2 / day]
906 : !> \param[in] "logical, optional :: log_data" ratio using logarithmic dataWorks
907 : !> only with 1d double precision input data.
908 :
909 : ! HISTORY
910 : !> \authors Juliane Mai
911 :
912 : !> \date Jun 2015
913 :
914 : ! Modifications:
915 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
916 :
917 0 : FUNCTION RunoffRatio(data, domain_area, mask, precip_series, precip_sum, log_data)
918 :
919 : implicit none
920 :
921 : ! array of data [m**3/s]
922 : real(dp), dimension(:), intent(in) :: data
923 :
924 : ! area of domain [km**2]
925 : real(dp), intent(in) :: domain_area
926 :
927 : ! mask for data points given
928 : logical, dimension(size(data, 1)), optional, intent(in) :: mask
929 :
930 : ! daily precipitation values [mm/km**2 / day]
931 : real(dp), dimension(size(data, 1)), optional, intent(in) :: precip_series
932 :
933 : ! sum of daily precip. values of whole period[mm/km**2 / day]
934 : real(dp), optional, intent(in) :: precip_sum
935 :
936 : ! ratio using logarithmic dataWorks only with 1d double precision input data.
937 : logical, optional, intent(in) :: log_data
938 :
939 : ! sum(data) / sum(precip)
940 : real(dp) :: RunoffRatio
941 :
942 : ! mask of data
943 0 : logical, dimension(size(data, 1)) :: maske
944 :
945 : ! if logarithmic data are used --> sum(log(data)) / sum(precip)
946 : logical :: log_dat
947 :
948 0 : real(dp) :: sum_discharge
949 :
950 0 : real(dp) :: sum_precip
951 :
952 :
953 : ! checking optionals
954 0 : if (present(mask)) then
955 0 : maske = mask
956 : else
957 0 : maske = .true.
958 : end if
959 :
960 0 : if (present(log_data)) then
961 0 : log_dat = log_data
962 : else
963 : log_dat = .false.
964 : end if
965 :
966 0 : if ((present(precip_series) .and. present(precip_sum)) .or. &
967 : (.not. present(precip_series) .and. .not. present(precip_sum))) then
968 0 : call error_message('mo_signatures: RunoffRatio: Exactly one precipitation information', raise=.false.)
969 0 : call error_message(' (precipitation series or sum of precipitation) ', raise=.false.)
970 0 : call error_message(' has to be specified!')
971 : end if
972 :
973 0 : if (present(mask) .and. present(precip_sum)) then
974 0 : call error_message('mo_signatures: RunoffRatio: Already aggregated precipitation (precip_sum) and', raise=.false.)
975 0 : call error_message(' mask can not be used together.', raise=.false.)
976 0 : call error_message(' Precip_series should be used instead!')
977 : end if
978 :
979 : ! mhm output [m**3/s] --> required [mm/d]
980 : ! [m**3/s] / [km**2] = [m**3/(s km**2)]
981 : ! => [m**3/(s km**2) * 60*60*24/1000**2] = [m/d]
982 : ! => [m**3/(s km**2) * 60*60*24*1000/1000**2] = [mm/d]
983 : ! => [m**3/(s km**2) * 86.4 ] = [mm/d]
984 : ! => discharge value [m**3/s] / catchment area [km**2] * 86.4 [km**2 s/m**3 * mm/d]
985 0 : if (log_dat) then
986 0 : sum_discharge = sum(log(data) * 86.4_dp / domain_area, mask = maske)
987 : else
988 0 : sum_discharge = sum(data * 86.4_dp / domain_area, mask = maske)
989 : end if
990 :
991 0 : if (present(precip_sum)) then
992 0 : sum_precip = precip_sum
993 : else
994 0 : sum_precip = sum(precip_series, mask = maske)
995 : end if
996 :
997 0 : RunoffRatio = sum_discharge / sum_precip
998 :
999 0 : END FUNCTION RunoffRatio
1000 :
1001 : !-------------------------------------------------------------------------------
1002 : ! NAME
1003 : ! ZeroFlowRatio
1004 :
1005 : ! PURPOSE
1006 : !> \brief Ratio of zero values to total number of data points.
1007 :
1008 : !> \details An optional mask of data points can be specified.
1009 : !> ADDITIONAL INFORMATION
1010 : !> \return real(dp), dimension(size(lags,1)) :: ZeroFlowRatio — Ratio of zero values to total number
1011 : !> of data points
1012 : !> Used as hydrologic signature in
1013 : !> Zhang, Y., Vaze, J., Chiew, F. H. S., Teng, J., & Li, M. (2014).
1014 : !> Predicting hydrological signatures in ungauged catchments using spatial interpolation, index model, and
1015 : !> rainfall-runoff modelling.
1016 : !> Journal of Hydrology, 517(C), 936-948. doi:10.1016/j.jhydrol.2014.06.032
1017 :
1018 : ! INTENT(IN)
1019 : !> \param[in] "real(dp), dimension(:) :: data" array of data
1020 :
1021 : ! INTENT(IN), OPTIONAL
1022 : !> \param[in] "logical, dimension(size(data, 1)), optional :: mask" mask for data points givenWorks only with 1d
1023 : !> double precision input data.
1024 :
1025 : ! HISTORY
1026 : !> \authors Juliane Mai
1027 :
1028 : !> \date Jun 2015
1029 :
1030 : ! Modifications:
1031 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1032 :
1033 0 : FUNCTION ZeroFlowRatio(data, mask)
1034 :
1035 0 : use mo_utils, only : eq
1036 :
1037 : implicit none
1038 :
1039 : ! array of data
1040 : real(dp), dimension(:), intent(in) :: data
1041 :
1042 : ! mask for data points givenWorks only with 1d double precision input data.
1043 : logical, dimension(size(data, 1)), optional, intent(in) :: mask
1044 :
1045 : ! Autocorrelation of data at given lags
1046 : real(dp) :: ZeroFlowRatio
1047 :
1048 : ! mask of data
1049 0 : logical, dimension(size(data, 1)) :: maske
1050 :
1051 : ! total number of data points
1052 : integer(i4) :: nall
1053 :
1054 : ! number of zero data points
1055 : integer(i4) :: nzero
1056 :
1057 :
1058 : ! checking optionals
1059 0 : if (present(mask)) then
1060 0 : maske = mask
1061 : else
1062 0 : maske = .true.
1063 : end if
1064 :
1065 0 : nall = count(maske)
1066 0 : nzero = count(maske .and. (eq(data, 0.0_dp)))
1067 :
1068 0 : if (nall > 0) then
1069 0 : ZeroFlowRatio = real(nzero, dp) / real(nall, dp)
1070 : else
1071 0 : call error_message('mo_signatures: ZeroFlowRatio: all data points are masked')
1072 : end if
1073 :
1074 0 : END FUNCTION ZeroFlowRatio
1075 :
1076 : END MODULE mo_mrm_signatures
|
