LCOV - code coverage report
Current view: top level - mRM - mo_mrm_signatures.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 0 178 0.0 %
Date: 2024-04-30 08:53:32 Functions: 0 8 0.0 %

          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 &mdash; 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 &mdash; 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 &mdash; 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 &mdash; 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

Generated by: LCOV version 1.16