Line data    Source code

       1             : !> \file    mo_eckhardt_filter.f90
       2             : !> \brief \copybrief mo_eckhardt_filter
       3             : !> \details \copydetails mo_eckhardt_filter
       4             : 
       5             : !> \brief   Eckhardt filter for baseflow index calculation.
       6             : !> \details This module provides routines for the Eckardt filter to analyse discharge time series and extract the baseflow.
       7             : !!          The filter was proposed in Eckhardt (2008, doi: 10.1016/j.jhydrol.2008.01.005)
       8             : !> \version 0.1
       9             : !> \authors Sebastian Mueller
      10             : !> \authors Mariaines Di Dato
      11             : !> \date    Apr 2022
      12             : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
      13             : !! FORCES is released under the LGPLv3+ license \license_note
      14             : module mo_eckhardt_filter
      15             : 
      16             :   use mo_kind,       only: i4, dp
      17             :   use mo_moment,     only: mean
      18             :   use mo_percentile, only: percentile
      19             :   use mo_append,     only: append
      20             :   use mo_nelmin,     only: nelminrange
      21             :   use mo_message,    only: error_message
      22             : 
      23             :   implicit none
      24             :   private
      25             :   public :: fit_alpha
      26             :   public :: eckhardt_filter_fit
      27             :   public :: eckhardt_filter
      28             :   public :: weekly_average
      29             :   public :: BFI
      30             : 
      31             :   real(dp), allocatable :: temp_d7(:)
      32             :   logical, allocatable :: temp_qmin_mask(:), temp_mask(:)
      33             : 
      34             : contains
      35             : 
      36             :   !> \brief   Eckhardt filter for baseflow calculation from discharge time series with fitting.
      37             :   !> \details This filter was proposed in Eckhardt (2008, doi: 10.1016/j.jhydrol.2008.01.005)
      38             :   !> \return  baseflow
      39           0 :   function eckhardt_filter_fit(discharge, mask) result(baseflow)
      40             : 
      41             :     !> array with daily discharge
      42             :     real(dp), intent(in)  :: discharge(:)
      43             :     !> mask for daily discharge
      44             :     logical, intent(in), optional :: mask(:)
      45             : 
      46             :     real(dp), allocatable :: baseflow(:)
      47             : 
      48           0 :     real(dp) :: alpha
      49             : 
      50           0 :     alpha = fit_alpha(discharge, mask=mask)
      51           0 :     baseflow = eckhardt_filter(alpha, discharge, mask=mask)
      52             : 
      53           0 :   end function eckhardt_filter_fit
      54             : 
      55             :   !> \brief   Fitted alpha parameter for the Eckhardt filter.
      56             :   !> \details This filter was proposed in Eckhardt (2008, doi: 10.1016/j.jhydrol.2008.01.005)
      57             :   !> \return  alpha parameter for eckard filter
      58           4 :   real(dp) function fit_alpha(discharge, mask)
      59             : 
      60             :     !> array with daily discharge
      61             :     real(dp), intent(in)  :: discharge(:)
      62             :     !> mask for daily discharge
      63             :     logical, intent(in), optional :: mask(:)
      64             : 
      65           4 :     real(dp) :: alpha_out(1)
      66           4 :     real(dp), allocatable :: q_min(:), dummy(:)
      67           2 :     logical, dimension(size(discharge)) :: mask_
      68             :     integer(i4) :: i, j
      69             : 
      70       14602 :     mask_(:) = .true.
      71        7303 :     if ( present(mask) ) mask_ = mask
      72             : 
      73       14605 :     temp_d7 = weekly_average(discharge, mask=mask)
      74           2 :     allocate(q_min(0))
      75          42 :     do i = 1, size(temp_d7), 365
      76          40 :       j = min(i+364, size(temp_d7))
      77             :       ! only use values in mask
      78             :       ! TODO: do we need a threshold for number in mask here?
      79       14742 :       if ( any(mask_(i : j)) ) call append(q_min, minval(temp_d7(i : j), mask=mask_(i : j)))
      80             :     end do
      81           2 :     if ( size(q_min) < 2 ) call error_message("Eckhardt filter: Less than 2 years of discharge observations! (min. 10 recommended)")
      82             : 
      83           8 :     allocate(temp_qmin_mask(size(discharge)))
      84           4 :     allocate(temp_mask(size(discharge)))
      85       14604 :     temp_mask = mask_
      86       14604 :     temp_qmin_mask = (temp_d7 < percentile(q_min, 10.0_dp, mode_in=4))
      87       14604 :     temp_qmin_mask = temp_qmin_mask .and. temp_mask
      88             : 
      89             :     ! set values outside of mask to 1 in d7
      90       14605 :     allocate(dummy(count(.not.mask_)))
      91        3652 :     dummy(:) = 1.0_dp  ! [1.0_dp, i=1,count(.not.mask_)]
      92             :     temp_d7 = unpack( &
      93             :       vector=dummy, &
      94             :       mask=.not.mask_, &
      95             :       field=temp_d7 &
      96       43804 :     )
      97           2 :     deallocate(dummy)
      98             : 
      99             :     alpha_out = nelminrange( &
     100             :       func=func, &
     101             :       pstart=[0.9_dp], &
     102             :       prange=reshape([0._dp, 1._dp], [1, 2]) &
     103           2 :     )
     104           2 :     fit_alpha = alpha_out(1)
     105             : 
     106           2 :     deallocate(temp_qmin_mask)
     107           2 :     deallocate(temp_mask)
     108           2 :     deallocate(temp_d7)
     109             : 
     110           2 :   end function fit_alpha
     111             : 
     112             :   !> \brief   Eckhardt filter for baseflow calculation from discharge time series.
     113             :   !> \details This filter was proposed in Eckhardt (2008, doi: 10.1016/j.jhydrol.2008.01.005)
     114             :   !> \return  baseflow
     115        1508 :   function eckhardt_filter(alpha, discharge, mask) result(baseflow)
     116             : 
     117             :     !> filter parameter
     118             :     real(dp), intent(in)  :: alpha
     119             :     !> array with daily discharge
     120             :     real(dp), intent(in)  :: discharge(:)
     121             :     !> mask for daily discharge
     122             :     logical, intent(in), optional :: mask(:)
     123             : 
     124             :     real(dp), allocatable :: baseflow(:)
     125             : 
     126         754 :     real(dp), allocatable :: d7(:), d7_perc(:), d_temp(:), b_temp(:)
     127         754 :     logical, dimension(size(discharge)) :: mask_
     128         754 :     real(dp) :: BFI_max
     129             :     integer(i4) :: i
     130             : 
     131     5504954 :     mask_(:) = .true.
     132     5497654 :     if ( present(mask) ) mask_ = mask
     133             : 
     134        2262 :     allocate(baseflow(size(discharge)))
     135        1508 :     allocate(d7_perc(size(discharge)))
     136             : 
     137             :     ! 20 percent percentile with Linear interpolation
     138     5505709 :     d7 = weekly_average(discharge, mask=mask)
     139     5504955 :     d7_perc(:) = percentile(d7, 20.0_dp, mask=mask, mode_in=4)
     140         754 :     BFI_max = BFI(d7_perc, discharge, mask=mask)
     141             : 
     142     5506462 :     allocate(b_temp(count(mask_)))
     143     5506462 :     allocate(d_temp(count(mask_)))
     144         754 :     d_temp = pack(discharge, mask=mask_)
     145             : 
     146             :     ! Applying the equation Eq. (6) from Eckhardt (2008) (only at mask)
     147         754 :     b_temp(1) = ((1 - alpha)*BFI_max * d_temp(1)) / (1 - alpha*BFI_max)
     148     4529650 :     do i = 2, size(d_temp)
     149     4529650 :       b_temp(i) = ((1 - BFI_max)*alpha*b_temp(i-1) + (1 - alpha)*BFI_max*d_temp(i)) / (1 - alpha*BFI_max)
     150             :     end do
     151     5504954 :     baseflow(:) = 0.0_dp
     152    11009908 :     baseflow = unpack(vector=b_temp, mask=mask_, field=baseflow)
     153             : 
     154         756 :   end function eckhardt_filter
     155             : 
     156             :   !> \brief   This function returns the 7days-averaged discharge.
     157             :   !> \return  array with weekly moving average
     158        1512 :   function weekly_average(discharge, mask) result(d7)
     159             : 
     160             :     !> array with daily discharge
     161             :     real(dp), intent(in) :: discharge(:)
     162             :     !> mask for daily discharge
     163             :     logical, intent(in), optional :: mask(:)
     164             : 
     165             :     real(dp), allocatable :: d7(:)
     166             : 
     167         756 :     logical, dimension(size(discharge)) :: mask_
     168             :     integer(i4) :: i, n, m
     169             : 
     170     5519556 :     mask_(:) = .true.
     171     5505710 :     if ( present(mask) ) mask_ = mask
     172             : 
     173        2268 :     allocate(d7(size(discharge)))
     174     5519556 :     d7(:) = 0.0_dp
     175             : 
     176     5519556 :     do i = 1, size(discharge)
     177     5518800 :       n = max(1,i-3)
     178     5518800 :       m = min(size(discharge),i+3)
     179             :       ! TODO: do we need a threshold for number in mask here?
     180    17878520 :       if ( any(mask_(n : m)) ) d7(i) = mean(discharge(n : m), mask=mask_(n : m))
     181             :     end do
     182             : 
     183         754 :   end function weekly_average
     184             : 
     185             :   !> \brief   Calculate the baseflow index as ratio between baseflow and discharge.
     186             :   !> \return  baseflow index
     187        1512 :   real(dp) function BFI(baseflow, discharge, mask)
     188             : 
     189             :     !> array with daily baseflow values
     190             :     real(dp), intent(in) :: baseflow(:)
     191             :     !> array with daily discharge
     192             :     real(dp), intent(in) :: discharge(:)
     193             :     !> mask for daily discharge
     194             :     logical, intent(in), optional :: mask(:)
     195             : 
     196         756 :     logical, dimension(size(discharge)) :: mask_
     197             : 
     198     5519556 :     mask_(:) = .true.
     199     5505710 :     if ( present(mask) ) mask_ = mask
     200             : 
     201    11039868 :     BFI = sum(baseflow, mask=mask_) / sum(discharge, mask=mask_)
     202             : 
     203         756 :   end function BFI
     204             : 
     205             :   !> \brief   Target function for fitting Eckhardt filter.
     206             :   !> \return  Objective value.
     207         752 :   function func(pp)
     208             : 
     209             :     real(dp), dimension(:), intent(in) :: pp  !< alpha (single value)
     210             : 
     211             :     real(dp) :: func
     212             : 
     213         752 :     real(dp), allocatable :: baseflow(:)
     214             : 
     215     5490352 :     baseflow = eckhardt_filter(alpha=pp(1), discharge=temp_d7, mask=temp_mask)
     216     5490352 :     func = mean((baseflow/temp_d7 - 1)**2, mask=temp_qmin_mask)
     217             : 
     218        1508 :   end function func
     219             : 
     220             : end module mo_eckhardt_filter