Line data    Source code

       1             : !> \file    mo_mhm_bfi.f90
       2             : !> \brief   \copybrief mo_mhm_bfi
       3             : !> \details \copydetails mo_mhm_bfi
       4             : 
       5             : !> \brief   Module to calculate BFI form gauging stations in mHM.
       6             : !> \version 0.1
       7             : !> \authors Sebastian Mueller
       8             : !> \date    Apr 2022
       9             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      10             : !! mHM is released under the LGPLv3+ license \license_note
      11             : !> \ingroup f_mhm
      12             : module mo_mhm_bfi
      13             : 
      14             :   use mo_kind,       only: i4, dp
      15             : 
      16             :   implicit none
      17             :   private
      18             :   public :: calculate_BFI
      19             : 
      20             : contains
      21             : 
      22             :   !> \brief Calculate BFI from given discharge observation.
      23           0 :   subroutine calculate_BFI()
      24             :     use mo_common_mhm_mrm_variables, only : evalPer, nTstepDay, warmingDays
      25             :     use mo_common_variables, only : domainMeta
      26             :     use mo_global_variables, only : BFI_obs
      27             :     use mo_message, only : error_message, message
      28             :     use mo_string_utils, only : num2str
      29             :     use mo_utils, only : ge
      30             :     use mo_mrm_global_variables, only : gauge, nMeasPerDay
      31             :     use mo_orderpack, only: sort_index
      32             :     use mo_eckhardt_filter, only: eckhardt_filter_fit, BFI
      33             : 
      34             :     implicit none
      35           0 :     real(dp), allocatable :: baseflow(:), runoff_obs(:)
      36           0 :     logical, allocatable :: runoff_obs_mask(:)
      37           0 :     integer(i4), allocatable :: id_sort(:)
      38             :     integer(i4) :: iDomain, id, length, mask_length, tt
      39             : 
      40             :     ! check for daily timesteps
      41           0 :     if ( nMeasPerDay > 1 ) call error_message("BFI: calculation only possible with daily discharge values!")
      42             : 
      43             :     ! extract domain Id from gauge Id
      44           0 :     if (size(gauge%domainId) /= domainMeta%nDomains) call error_message("BFI: number of gauges and domains need to be equal!")
      45           0 :     allocate(id_sort(size(gauge%domainId)))
      46           0 :     id_sort = sort_index(gauge%domainId)
      47             : 
      48           0 :     call message()
      49           0 :     do iDomain = 1, domainMeta%nDomains
      50             :       ! skip current calculation if BFI is given
      51           0 :       if ( .not. (BFI_obs(iDomain) < 0.0_dp) ) then
      52             :         call message( &
      53             :           "  BFI: using given BFI value ", &
      54             :           trim(adjustl(num2str(BFI_obs(iDomain)))), &
      55             :           " for domain ", &
      56             :           trim(adjustl(num2str(iDomain))) &
      57           0 :         )
      58           0 :         cycle
      59             :       end if
      60             : 
      61           0 :       id = id_sort(iDomain)
      62           0 :       if (gauge%domainId(id) /= iDomain) call error_message("BFI: exactly one gauge per domain required!")
      63             : 
      64             :       ! get length of evaluation period times
      65           0 :       length = (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1) * nMeasPerDay
      66             : 
      67             :       ! extract measurements
      68           0 :       if (allocated(runoff_obs_mask)) deallocate(runoff_obs_mask)
      69           0 :       if (allocated(runoff_obs)) deallocate(runoff_obs)
      70           0 :       if (allocated(baseflow)) deallocate(baseflow)
      71           0 :       allocate(runoff_obs_mask(length))
      72           0 :       allocate(runoff_obs(length))
      73           0 :       runoff_obs = gauge%Q(1 : length, id)
      74             : 
      75             :       ! create mask of observed runoff
      76           0 :       forall(tt = 1 : length) runoff_obs_mask(tt) = ge(runoff_obs(tt), 0.0_dp)
      77             : 
      78             :       ! calculate BFI
      79           0 :       baseflow = eckhardt_filter_fit(runoff_obs, mask=runoff_obs_mask)
      80           0 :       BFI_obs(iDomain) = BFI(baseflow, runoff_obs, mask=runoff_obs_mask)
      81             : 
      82             :       call message( &
      83             :         "  BFI:  calculated BFI value ", &
      84           0 :         trim(adjustl(num2str(BFI_obs(iDomain)))), &
      85             :         " for domain ", &
      86             :         trim(adjustl(num2str(iDomain))) &
      87           0 :       )
      88             :     end do
      89             : 
      90           0 :   end subroutine calculate_BFI
      91             : 
      92             : end module mo_mhm_bfi