mo_mhm_bfi.f90

Go to the documentation of this file.

!> \file    mo_mhm_bfi.f90
!> \brief   \copybrief mo_mhm_bfi
!> \details \copydetails mo_mhm_bfi
 
!> \brief   Module to calculate BFI form gauging stations in mHM.
!> \version 0.1
!> \authors Sebastian Mueller
!> \date    Apr 2022
!> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
!! mHM is released under the LGPLv3+ license \license_note
!> \ingroup f_mhm
module mo_mhm_bfi
 
  use mo_kind,       only: i4, dp
 
  implicit none
  private
  public :: calculate_bfi
 
contains
 
  !> \brief Calculate BFI from given discharge observation.
  subroutine calculate_bfi()
    use mo_common_mhm_mrm_variables, only : evalper, ntstepday, warmingdays
    use mo_common_variables, only : domainmeta
    use mo_global_variables, only : bfi_obs
    use mo_message, only : error_message, message
    use mo_string_utils, only : num2str
    use mo_utils, only : ge
    use mo_mrm_global_variables, only : gauge, nmeasperday
    use mo_orderpack, only: sort_index
    use mo_eckhardt_filter, only: eckhardt_filter_fit, bfi
 
    implicit none
    real(dp), allocatable :: baseflow(:), runoff_obs(:)
    logical, allocatable :: runoff_obs_mask(:)
    integer(i4), allocatable :: id_sort(:)
    integer(i4) :: idomain, id, length, mask_length, tt
 
    ! check for daily timesteps
    if ( nmeasperday > 1 ) call error_message("BFI: calculation only possible with daily discharge values!")
 
    ! extract domain Id from gauge Id
    if (size(gauge%domainId) /= domainmeta%nDomains) call error_message("BFI: number of gauges and domains need to be equal!")
    allocate(id_sort(size(gauge%domainId)))
    id_sort = sort_index(gauge%domainId)
 
    call message()
    do idomain = 1, domainmeta%nDomains
      ! skip current calculation if BFI is given
      if ( .not. (bfi_obs(idomain) < 0.0_dp) ) then
        call message( &
          "  BFI: using given BFI value ", &
          trim(adjustl(num2str(bfi_obs(idomain)))), &
          " for domain ", &
          trim(adjustl(num2str(idomain))) &
        )
        cycle
      end if
 
      id = id_sort(idomain)
      if (gauge%domainId(id) /= idomain) call error_message("BFI: exactly one gauge per domain required!")
 
      ! get length of evaluation period times
      length = (evalper(idomain)%julEnd - evalper(idomain)%julStart + 1) * nmeasperday
 
      ! extract measurements
      if (allocated(runoff_obs_mask)) deallocate(runoff_obs_mask)
      if (allocated(runoff_obs)) deallocate(runoff_obs)
      if (allocated(baseflow)) deallocate(baseflow)
      allocate(runoff_obs_mask(length))
      allocate(runoff_obs(length))
      runoff_obs = gauge%Q(1 : length, id)
 
      ! create mask of observed runoff
      forall(tt = 1 : length) runoff_obs_mask(tt) = ge(runoff_obs(tt), 0.0_dp)
 
      ! calculate BFI
      baseflow = eckhardt_filter_fit(runoff_obs, mask=runoff_obs_mask)
      bfi_obs(idomain) = bfi(baseflow, runoff_obs, mask=runoff_obs_mask)
 
      call message( &
        "  BFI:  calculated BFI value ", &
        trim(adjustl(num2str(bfi_obs(idomain)))), &
        " for domain ", &
        trim(adjustl(num2str(idomain))) &
      )
    end do
 
  end subroutine calculate_bfi
 
end module mo_mhm_bfi