mo_eckhardt_filter.f90 6.59 KB
Newer Older
Sebastian Müller's avatar
Sebastian Müller committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
!> \file    mo_eckhardt_filter.f90
!> \copydoc mo_eckhardt_filter

!> \brief   Eckhardt filter for baseflow index calculation.
!> \details This module provides routines for the Eckardt filter to analyse discharge time series and extract the baseflow.
!!          The filter was proposed in Eckhardt (2008, doi: 10.1016/j.jhydrol.2008.01.005)
!> \version 0.1
!> \authors Sebastian Mueller
!> \authors Mariaines Di Dato
!> \date    Apr 2022
module mo_eckhardt_filter

  use mo_kind,       only: i4, dp
  use mo_moment,     only: mean
  use mo_percentile, only: percentile
  use mo_append,     only: append
  use mo_nelmin,     only: nelminrange
18
  use mo_message,    only: error_message
Sebastian Müller's avatar
Sebastian Müller committed
19
20
21
22
23
24
25
26
27
28

  implicit none
  private
  public :: fit_alpha
  public :: eckhardt_filter_fit
  public :: eckhardt_filter
  public :: weekly_average
  public :: BFI

  real(dp), allocatable :: temp_d7(:)
29
  logical, allocatable :: temp_qmin_mask(:), temp_mask(:)
Sebastian Müller's avatar
Sebastian Müller committed
30
31
32
33
34
35

contains

  !> \brief   Eckhardt filter for baseflow calculation from discharge time series with fitting.
  !> \details This filter was proposed in Eckhardt (2008, doi: 10.1016/j.jhydrol.2008.01.005)
  !> \return  baseflow
36
  function eckhardt_filter_fit(discharge, mask) result(baseflow)
Sebastian Müller's avatar
Sebastian Müller committed
37
38
39

    !> array with daily discharge
    real(dp), intent(in)  :: discharge(:)
40
41
    !> mask for daily discharge
    logical, intent(in), optional :: mask(:)
Sebastian Müller's avatar
Sebastian Müller committed
42
43
44
45
46

    real(dp), allocatable :: baseflow(:)

    real(dp) :: alpha

47
48
    alpha = fit_alpha(discharge, mask=mask)
    baseflow = eckhardt_filter(alpha, discharge, mask=mask)
Sebastian Müller's avatar
Sebastian Müller committed
49
50
51
52
53
54

  end function eckhardt_filter_fit

  !> \brief   Fitted alpha parameter for the Eckhardt filter.
  !> \details This filter was proposed in Eckhardt (2008, doi: 10.1016/j.jhydrol.2008.01.005)
  !> \return  alpha parameter for eckard filter
55
  real(dp) function fit_alpha(discharge, mask)
Sebastian Müller's avatar
Sebastian Müller committed
56
57
58

    !> array with daily discharge
    real(dp), intent(in)  :: discharge(:)
59
60
    !> mask for daily discharge
    logical, intent(in), optional :: mask(:)
Sebastian Müller's avatar
Sebastian Müller committed
61
62

    real(dp) :: alpha_out(1)
63
64
65
66
67
68
    real(dp), allocatable :: q_min(:), dummy(:)
    logical, dimension(size(discharge)) :: mask_
    integer(i4) :: i, j

    mask_(:) = .true.
    if ( present(mask) ) mask_ = mask
Sebastian Müller's avatar
Sebastian Müller committed
69

70
71
    temp_d7 = weekly_average(discharge, mask=mask)
    allocate(q_min(0))
Sebastian Müller's avatar
Sebastian Müller committed
72
    do i = 1, size(temp_d7), 365
73
74
75
76
      j = min(i+364, size(temp_d7))
      ! only use values in mask
      ! TODO: do we need a threshold for number in mask here?
      if ( any(mask_(i : j)) ) call append(q_min, minval(temp_d7(i : j), mask=mask_(i : j)))
Sebastian Müller's avatar
Sebastian Müller committed
77
    end do
78
    if ( size(q_min) < 2 ) call error_message("Eckhardt filter: Less than 2 years of discharge observations! (min. 10 recommended)")
Sebastian Müller's avatar
Sebastian Müller committed
79
80

    allocate(temp_qmin_mask(size(discharge)))
81
82
    allocate(temp_mask(size(discharge)))
    temp_mask = mask_
Sebastian Müller's avatar
Sebastian Müller committed
83
    temp_qmin_mask = (temp_d7 < percentile(q_min, 10.0_dp, mode_in=4))
84
85
86
87
88
89
90
91
92
93
94
    temp_qmin_mask = temp_qmin_mask .and. temp_mask

    ! set values outside of mask to 1 in d7
    allocate(dummy(count(.not.mask_)))
    dummy(:) = 1.0_dp  ! [1.0_dp, i=1,count(.not.mask_)]
    temp_d7 = unpack( &
      vector=dummy, &
      mask=.not.mask_, &
      field=temp_d7 &
    )
    deallocate(dummy)
Sebastian Müller's avatar
Sebastian Müller committed
95
96
97
98
99
100
101
102
103

    alpha_out = nelminrange( &
      func=func, &
      pstart=[0.9_dp], &
      prange=reshape([0._dp, 1._dp], [1, 2]) &
    )
    fit_alpha = alpha_out(1)

    deallocate(temp_qmin_mask)
104
    deallocate(temp_mask)
Sebastian Müller's avatar
Sebastian Müller committed
105
106
107
108
109
110
111
    deallocate(temp_d7)

  end function fit_alpha

  !> \brief   Eckhardt filter for baseflow calculation from discharge time series.
  !> \details This filter was proposed in Eckhardt (2008, doi: 10.1016/j.jhydrol.2008.01.005)
  !> \return  baseflow
112
  function eckhardt_filter(alpha, discharge, mask) result(baseflow)
Sebastian Müller's avatar
Sebastian Müller committed
113
114
115
116
117

    !> filter parameter
    real(dp), intent(in)  :: alpha
    !> array with daily discharge
    real(dp), intent(in)  :: discharge(:)
118
119
    !> mask for daily discharge
    logical, intent(in), optional :: mask(:)
Sebastian Müller's avatar
Sebastian Müller committed
120
121
122

    real(dp), allocatable :: baseflow(:)

123
124
    real(dp), allocatable :: d7(:), d7_perc(:), d_temp(:), b_temp(:)
    logical, dimension(size(discharge)) :: mask_
Sebastian Müller's avatar
Sebastian Müller committed
125
126
127
    real(dp) :: BFI_max
    integer(i4) :: i

128
129
130
    mask_(:) = .true.
    if ( present(mask) ) mask_ = mask

Sebastian Müller's avatar
Sebastian Müller committed
131
132
133
134
    allocate(baseflow(size(discharge)))
    allocate(d7_perc(size(discharge)))

    ! 20 percent percentile with Linear interpolation
135
136
137
138
139
140
141
142
143
144
145
146
    d7 = weekly_average(discharge, mask=mask)
    d7_perc(:) = percentile(d7, 20.0_dp, mask=mask, mode_in=4)
    BFI_max = BFI(d7_perc, discharge, mask=mask)

    allocate(b_temp(count(mask_)))
    allocate(d_temp(count(mask_)))
    d_temp = pack(discharge, mask=mask_)

    ! Applying the equation Eq. (6) from Eckhardt (2008) (only at mask)
    b_temp(1) = ((1 - alpha)*BFI_max * d_temp(1)) / (1 - alpha*BFI_max)
    do i = 2, size(d_temp)
      b_temp(i) = ((1 - BFI_max)*alpha*b_temp(i-1) + (1 - alpha)*BFI_max*d_temp(i)) / (1 - alpha*BFI_max)
Sebastian Müller's avatar
Sebastian Müller committed
147
    end do
148
149
    baseflow(:) = 0.0_dp
    baseflow = unpack(vector=b_temp, mask=mask_, field=baseflow)
Sebastian Müller's avatar
Sebastian Müller committed
150
151
152
153
154

  end function eckhardt_filter

  !> \brief   This function returns the 7days-averaged discharge.
  !> \return  array with weekly moving average
155
  function weekly_average(discharge, mask) result(d7)
Sebastian Müller's avatar
Sebastian Müller committed
156
157
158

    !> array with daily discharge
    real(dp), intent(in) :: discharge(:)
159
160
    !> mask for daily discharge
    logical, intent(in), optional :: mask(:)
Sebastian Müller's avatar
Sebastian Müller committed
161
162
163

    real(dp), allocatable :: d7(:)

164
165
166
167
168
    logical, dimension(size(discharge)) :: mask_
    integer(i4) :: i, n, m

    mask_(:) = .true.
    if ( present(mask) ) mask_ = mask
Sebastian Müller's avatar
Sebastian Müller committed
169
170

    allocate(d7(size(discharge)))
171
    d7(:) = 0.0_dp
Sebastian Müller's avatar
Sebastian Müller committed
172
173

    do i = 1, size(discharge)
174
175
176
      n = max(1,i-3)
      m = min(size(discharge),i+3)
      ! TODO: do we need a threshold for number in mask here?
177
      if ( any(mask_(n : m)) ) d7(i) = mean(discharge(n : m), mask=mask_(n : m))
Sebastian Müller's avatar
Sebastian Müller committed
178
179
180
181
182
183
    end do

  end function weekly_average

  !> \brief   Calculate the baseflow index as ratio between baseflow and discharge.
  !> \return  baseflow index
184
  real(dp) function BFI(baseflow, discharge, mask)
Sebastian Müller's avatar
Sebastian Müller committed
185
186
187
188
189

    !> array with daily baseflow values
    real(dp), intent(in) :: baseflow(:)
    !> array with daily discharge
    real(dp), intent(in) :: discharge(:)
190
191
    !> mask for daily discharge
    logical, intent(in), optional :: mask(:)
Sebastian Müller's avatar
Sebastian Müller committed
192

193
    BFI = sum(baseflow, mask=mask) / sum(discharge, mask=mask)
Sebastian Müller's avatar
Sebastian Müller committed
194
195
196
197
198
199
200
201
202
203
204
205
206

  end function BFI

  !> \brief   Target function for fitting Eckhardt filter.
  !> \return  Objective value.
  function func(pp)

    real(dp), dimension(:), intent(in) :: pp  !< alpha (single value)

    real(dp) :: func

    real(dp), allocatable :: baseflow(:)

207
    baseflow = eckhardt_filter(alpha=pp(1), discharge=temp_d7, mask=temp_mask)
Sebastian Müller's avatar
Sebastian Müller committed
208
209
210
211
212
    func = mean((baseflow/temp_d7 - 1)**2, mask=temp_qmin_mask)

  end function func

end module mo_eckhardt_filter