Merge branch 'opti_add_BFI_output' into 'develop'

BFI: add eckhardt filter and optimization option

Closes #32

See merge request !47

src/mo_eckhardt_filter.f90

+!> \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
+  use mo_message,    only: error_message
+
+  implicit none
+  private
+  public :: fit_alpha
+  public :: eckhardt_filter_fit
+  public :: eckhardt_filter
+  public :: weekly_average
+  public :: BFI
+
+  real(dp), allocatable :: temp_d7(:)
+  logical, allocatable :: temp_qmin_mask(:), temp_mask(:)
+
+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
+  function eckhardt_filter_fit(discharge, mask) result(baseflow)
+
+    !> array with daily discharge
+    real(dp), intent(in)  :: discharge(:)
+    !> mask for daily discharge
+    logical, intent(in), optional :: mask(:)
+
+    real(dp), allocatable :: baseflow(:)
+
+    real(dp) :: alpha
+
+    alpha = fit_alpha(discharge, mask=mask)
+    baseflow = eckhardt_filter(alpha, discharge, mask=mask)
+
+  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
+  real(dp) function fit_alpha(discharge, mask)
+
+    !> array with daily discharge
+    real(dp), intent(in)  :: discharge(:)
+    !> mask for daily discharge
+    logical, intent(in), optional :: mask(:)
+
+    real(dp) :: alpha_out(1)
+    real(dp), allocatable :: q_min(:), dummy(:)
+    logical, dimension(size(discharge)) :: mask_
+    integer(i4) :: i, j
+
+    mask_(:) = .true.
+    if ( present(mask) ) mask_ = mask
+
+    temp_d7 = weekly_average(discharge, mask=mask)
+    allocate(q_min(0))
+    do i = 1, size(temp_d7), 365
+      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)))
+    end do
+    if ( size(q_min) < 2 ) call error_message("Eckhardt filter: Less than 2 years of discharge observations! (min. 10 recommended)")
+
+    allocate(temp_qmin_mask(size(discharge)))
+    allocate(temp_mask(size(discharge)))
+    temp_mask = mask_
+    temp_qmin_mask = (temp_d7 < percentile(q_min, 10.0_dp, mode_in=4))
+    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)
+
+    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)
+    deallocate(temp_mask)
+    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
+  function eckhardt_filter(alpha, discharge, mask) result(baseflow)
+
+    !> filter parameter
+    real(dp), intent(in)  :: alpha
+    !> array with daily discharge
+    real(dp), intent(in)  :: discharge(:)
+    !> mask for daily discharge
+    logical, intent(in), optional :: mask(:)
+
+    real(dp), allocatable :: baseflow(:)
+
+    real(dp), allocatable :: d7(:), d7_perc(:), d_temp(:), b_temp(:)
+    logical, dimension(size(discharge)) :: mask_
+    real(dp) :: BFI_max
+    integer(i4) :: i
+
+    mask_(:) = .true.
+    if ( present(mask) ) mask_ = mask
+
+    allocate(baseflow(size(discharge)))
+    allocate(d7_perc(size(discharge)))
+
+    ! 20 percent percentile with Linear interpolation
+    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)
+    end do
+    baseflow(:) = 0.0_dp
+    baseflow = unpack(vector=b_temp, mask=mask_, field=baseflow)
+
+  end function eckhardt_filter
+
+  !> \brief   This function returns the 7days-averaged discharge.
+  !> \return  array with weekly moving average
+  function weekly_average(discharge, mask) result(d7)
+
+    !> array with daily discharge
+    real(dp), intent(in) :: discharge(:)
+    !> mask for daily discharge
+    logical, intent(in), optional :: mask(:)
+
+    real(dp), allocatable :: d7(:)
+
+    logical, dimension(size(discharge)) :: mask_
+    integer(i4) :: i, n, m
+
+    mask_(:) = .true.
+    if ( present(mask) ) mask_ = mask
+
+    allocate(d7(size(discharge)))
+    d7(:) = 0.0_dp
+
+    do i = 1, size(discharge)
+      n = max(1,i-3)
+      m = min(size(discharge),i+3)
+      ! TODO: do we need a threshold for number in mask here?
+      if ( any(mask_(n : m)) ) d7(i) = mean(discharge(n : m), mask=mask_(n : m))
+    end do
+
+  end function weekly_average
+
+  !> \brief   Calculate the baseflow index as ratio between baseflow and discharge.
+  !> \return  baseflow index
+  real(dp) function BFI(baseflow, discharge, mask)
+
+    !> array with daily baseflow values
+    real(dp), intent(in) :: baseflow(:)
+    !> array with daily discharge
+    real(dp), intent(in) :: discharge(:)
+    !> mask for daily discharge
+    logical, intent(in), optional :: mask(:)
+
+    logical, dimension(size(discharge)) :: mask_
+
+    mask_(:) = .true.
+    if ( present(mask) ) mask_ = mask
+
+    BFI = sum(baseflow, mask=mask_) / sum(discharge, mask=mask_)
+
+  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(:)
+
+    baseflow = eckhardt_filter(alpha=pp(1), discharge=temp_d7, mask=temp_mask)
+    func = mean((baseflow/temp_d7 - 1)**2, mask=temp_qmin_mask)
+
+  end function func
+
+end module mo_eckhardt_filter

src/mo_io.f90

+!> \file    mo_io.f90
+!> \copydoc mo_io
+
+!> \brief   File reading routines.
+!> \details This module provides routines to load a file into an array.
+!!          This is mainly taken from the Fortran stdlib: https://github.com/fortran-lang/stdlib
+!> \version 0.1
+!> \authors Sebastian Mueller
+!> \date    Apr 2022
+module mo_io
+
+  use mo_kind, only: i4, dp, sp
+  use mo_message, only: error_message
+  use mo_string_utils, only: is_blank
+
+  implicit none
+  private
+  public :: loadtxt
+  public :: number_of_columns
+  public :: number_of_rows
+
+  !> \brief   Read a file into a 2D array containing reals.
+  interface loadtxt
+    module procedure loadtxt_dp, loadtxt_sp
+  end interface loadtxt
+
+contains
+
+  !> \brief   Read a file into a 2D array containing reals.
+  subroutine loadtxt_dp(filename, d, skiprows, max_rows)
+
+    !> Filename to load the array from
+    character(len=*), intent(in) :: filename
+    !> The array 'd' will be automatically allocated with the correct dimensions
+    real(dp), allocatable, intent(out) :: d(:,:)
+    !> lines to skip at the begining
+    integer(i4), intent(in), optional :: skiprows
+    !> Read max_rows lines of content after skiprows lines. The default is to read all the lines (negative values).
+    integer(i4), intent(in), optional :: max_rows
+
+    integer(i4) :: u,  nrow, ncol, i, skiprows_, max_rows_
+
+    skiprows_ = 0
+    if ( present(skiprows) ) skiprows_ = max(skiprows, 0)
+
+    max_rows_ = -1
+    if ( present(max_rows) ) max_rows_ = max_rows
+
+    open(newunit=u, file=filename, action='read', position='asis', status='old', access='stream', form='formatted')
+
+    ! determine size
+    ncol = number_of_columns(u)
+    nrow = number_of_rows(u)
+    skiprows_ = min(skiprows_, nrow)
+    if ( max_rows_ < 0 .or. max_rows_ > (nrow - skiprows_)) max_rows_ = nrow - skiprows_
+
+    allocate(d(max_rows_, ncol))
+
+    do i = 1, skiprows_
+      read(u, *)
+    end do
+
+    do i = 1, max_rows_
+      read(u, "(*(es24.16e3,1x))") d(i, :)
+    end do
+
+    close(u)
+
+  end subroutine loadtxt_dp
+
+  !> \brief   Read a file into a 2D array containing reals.
+  subroutine loadtxt_sp(filename, d, skiprows, max_rows)
+
+    !> Filename to load the array from
+    character(len=*), intent(in) :: filename
+    !> The array 'd' will be automatically allocated with the correct dimensions
+    real(sp), allocatable, intent(out) :: d(:,:)
+    !> lines to skip at the begining
+    integer(i4), intent(in), optional :: skiprows
+    !> Read max_rows lines of content after skiprows lines. The default is to read all the lines (negative values).
+    integer(i4), intent(in), optional :: max_rows
+
+    integer(i4) :: u, nrow, ncol, i, skiprows_, max_rows_
+
+    skiprows_ = 0
+    if ( present(skiprows) ) skiprows_ = max(skiprows, 0)
+
+    max_rows_ = -1
+    if ( present(max_rows) ) max_rows_ = max_rows
+
+    open(newunit=u, file=filename, action='read', position='asis', status='old', access='stream', form='formatted')
+
+    ! determine size
+    ncol = number_of_columns(u)
+    nrow = number_of_rows(u)
+    skiprows_ = min(skiprows_, nrow)
+    if ( max_rows_ < 0 .or. max_rows_ > (nrow - skiprows_)) max_rows_ = nrow - skiprows_
+
+    allocate(d(max_rows_, ncol))
+
+    do i = 1, skiprows_
+      read(u, *)
+    end do
+
+    do i = 1, max_rows_
+      read(u, "(*(es15.8e2,1x))") d(i, :)
+    end do
+
+    close(u)
+
+  end subroutine loadtxt_sp
+
+  !> \brief   Determine number of columns in a file. The columns are assumed to be separated by spaces or tabs.
+  !> \return  Number of columns.
+  integer function number_of_columns(u)
+
+    integer(i4),intent(in) :: u !< unit of the open file
+    integer(i4) :: ios
+    character :: c
+    logical :: lastblank
+
+    rewind(u)
+    number_of_columns = 0
+    lastblank = .true.
+    do
+      read(u, '(a)', advance='no', iostat=ios) c
+      if (ios /= 0) exit
+      if (lastblank .and. .not. is_blank(c)) number_of_columns = number_of_columns + 1
+      lastblank = is_blank(c)
+    end do
+    rewind(u)
+
+  end function number_of_columns
+
+  !> \brief   Determine number of rows in a file.
+  !> \return  Number of rows.
+  integer function number_of_rows(u)
+
+    integer(i4), intent(in) :: u !< unit of the open file
+    integer(i4) :: ios
+
+    rewind(u)
+    number_of_rows = 0
+    do
+      read(u, *, iostat=ios)
+      if (ios /= 0) exit
+      number_of_rows = number_of_rows + 1
+    end do
+    rewind(u)
+
+  end function number_of_rows
+
+end module mo_io

src/mo_likelihood.f90

@@ -147,7 +147,7 @@ CONTAINS
 
   ! -------------------------------
   !> \brief A Model: p1*x^2 + p2*x + p3
-  subroutine model_dp(parameterset, opti_domain_indices, runoff, smOptiSim, neutronsOptiSim, etOptiSim, twsOptiSim)
+  subroutine model_dp(parameterset, opti_domain_indices, runoff, smOptiSim, neutronsOptiSim, etOptiSim, twsOptiSim, BFI)
 
     use mo_kind, only: dp
     use mo_optimization_types, only: optidata_sim
@@ -160,6 +160,7 @@ CONTAINS
     type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim ! dim1=ncells, dim2=time
     type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim       ! dim1=ncells, dim2=time
     type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim      ! dim1=ncells, dim2=time
+    real(dp),    dimension(:), allocatable, optional, intent(out) :: BFI         !< baseflow index, dim1=domainID
     integer(i4) :: i, n
     ! for OMP
     !! !$  integer(i4)                           :: n_threads, is_thread

src/mo_opt_functions.f90

@@ -5708,7 +5708,7 @@ CONTAINS
   end function griewank_objective
 
 
- subroutine eval_dummy(parameterset, opti_domain_indices, runoff, smOptiSim, neutronsOptiSim, etOptiSim, twsOptiSim)
+ subroutine eval_dummy(parameterset, opti_domain_indices, runoff, smOptiSim, neutronsOptiSim, etOptiSim, twsOptiSim, BFI)
     use mo_kind, only : dp
     use mo_optimization_types, only : optidata_sim, optidata
 
@@ -5721,6 +5721,7 @@ CONTAINS
     type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim ! dim1=ncells, dim2=time
     type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim       ! dim1=ncells, dim2=time
     type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim      ! dim1=ncells, dim2=time
+    real(dp),    dimension(:), allocatable, optional, intent(out) :: BFI         !< baseflow index, dim1=domainID
 
     type(optidata) :: dummyData
     integer(i4) :: i
@@ -5757,6 +5758,11 @@ CONTAINS
       runoff(:, :) = 0.0_dp
     end if
 
+    if (present(BFI)) then
+      allocate(BFI(1))
+      BFI(:) = 0.0_dp
+    end if
+
   end subroutine eval_dummy
 
 END MODULE mo_opt_functions

src/mo_optimization_utils.f90

@@ -10,20 +10,21 @@ module mo_optimization_utils
 
   !> \brief Interface for evaluation routine.
   abstract interface
-    subroutine eval_interface(parameterset, opti_domain_indices, runoff, smOptiSim, neutronsOptiSim, etOptiSim, twsOptiSim)
+    subroutine eval_interface(parameterset, opti_domain_indices, runoff, smOptiSim, neutronsOptiSim, etOptiSim, twsOptiSim, BFI)
       use mo_kind, only : dp, i4
       use mo_optimization_types, only : optidata_sim
       real(dp),    dimension(:), intent(in) :: parameterset
       integer(i4), dimension(:),                 optional, intent(in)  :: opti_domain_indices
-      real(dp),    dimension(:, :), allocatable, optional, intent(out) :: runoff        !< dim1=time dim2=gauge
+      real(dp),    dimension(:, :), allocatable, optional, intent(out) :: runoff   !< dim1=time dim2=gauge
       type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim       !< dim1=ncells, dim2=time
       type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim !< dim1=ncells, dim2=time
       type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim       !< dim1=ncells, dim2=time
       type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim      !< dim1=ncells, dim2=time
+      real(dp),    dimension(:), allocatable, optional, intent(out) :: BFI         !< baseflow index, dim1=domainID
     end subroutine
   end interface
 
-    !> \brief Interface for objective function.
+  !> \brief Interface for objective function.
   interface
     function objective_interface (parameterset, eval, arg1, arg2, arg3)
       use mo_kind, only : dp
@@ -39,5 +40,3 @@ module mo_optimization_utils
   end interface
 
 end module mo_optimization_utils
-
-

src/mo_string_utils.f90

@@ -48,6 +48,7 @@ MODULE mo_string_utils
   PUBLIC :: Replace_Text  ! replaces all text occurences in string
   PUBLIC :: replace_word  ! replaces all word occurences in string
   PUBLIC :: index_word    ! yields starting position of word in string or 0
+  PUBLIC :: is_blank
 
   ! public :: numarray2str
 
@@ -137,6 +138,19 @@ MODULE mo_string_utils
 
 CONTAINS
 
+  !> \brief   Check for blank characters.
+  !> \details Checks whether or not `c` is a blank character, namely a space and tab character.
+  !> \return  Truth value if `c` is a blank.
+  pure logical function is_blank(c)
+
+    character(len=1), intent(in) :: c !< The character to test.
+    integer :: ic
+
+    ic = iachar(c)             ! TAB
+    is_blank = (c == ' ') .or. (ic == int(z'09'));
+
+  end function is_blank
+
   ! ------------------------------------------------------------------
 
   !>    \brief Remove white spaces

src/pf_tests/test_mo_eckhardt_filter.pf

+module test_mo_eckhardt_filter
+  use funit
+  use mo_eckhardt_filter, only: eckhardt_filter, fit_alpha, BFI
+  use mo_io, only: loadtxt
+  use mo_kind, only : dp
+
+  implicit none
+
+  private
+
+  public :: test_fit
+  public :: test_fit_mask
+
+contains
+
+  @test
+  subroutine test_fit()
+    real(dp), allocatable :: dat(:,:), baseflow(:)
+    real(dp) :: alpha_out, BFI_out
+
+    call loadtxt("./files/discharge.txt", dat, skiprows=2)
+
+    alpha_out = fit_alpha(dat(:,1))
+    baseflow = eckhardt_filter(alpha_out, dat(:,1))
+    BFI_out = BFI(baseflow, dat(:,1))
+
+    @assertEqual(0.99718_dp, alpha_out, tolerance=0.0001, message='Check alpha')
+    @assertEqual(0.40635_dp, BFI_out, tolerance=0.0001, message='Check BFI')
+
+  end subroutine test_fit
+
+  @test
+  subroutine test_fit_mask()
+    real(dp), allocatable :: dat(:,:), baseflow(:)
+    logical, allocatable :: mask(:)
+    real(dp) :: alpha_out, BFI_out
+
+    call loadtxt("./files/discharge.txt", dat, skiprows=2)
+    allocate(mask(size(dat(:,1))))
+    mask(:) = .false.
+    ! only first 10 years
+    mask(1 : 365*10) = .true.
+    alpha_out = fit_alpha(dat(:,1), mask=mask)
+    baseflow = eckhardt_filter(alpha_out, dat(:,1), mask=mask)
+    BFI_out = BFI(baseflow, dat(:,1), mask=mask)
+
+    @assertEqual(0.99634_dp, alpha_out, tolerance=0.0001, message='Check alpha')
+    @assertEqual(0.43507_dp, BFI_out, tolerance=0.0001, message='Check BFI')
+
+  end subroutine test_fit_mask
+
+end module test_mo_eckhardt_filter