mo_mrm_mpr.f90

Go to the documentation of this file.

!> \file mo_mrm_mpr.f90
!> \brief \copybrief mo_mrm_mpr
!> \details \copydetails mo_mrm_mpr
 
!> \brief Perform Multiscale Parameter Regionalization on Routing Parameters
!> \details This module contains the subroutine for calculating the regionalized
!! routing parameters (beta-parameters) given the five global routing parameters (gamma) at the level 0 scale.
!> \authors Luis Samaniego, Stephan Thober
!> \date Aug 2015
!> \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_mrm
module mo_mrm_mpr
  use mo_kind, only : dp
  use mo_message, only : message
 
  implicit none
 
  public :: reg_rout
  public :: mrm_init_param
  public :: mrm_update_param
  private
contains
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        reg_rout
 
  !    PURPOSE
  !>       \brief Regionalized routing
 
  !>       \details sets up the Regionalized Routing parameters
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param(1) = muskingumTravelTime_constant
  !>       - param(2) = muskingumTravelTime_riverLength
  !>       - param(3) = muskingumTravelTime_riverSlope
  !>       - param(4) = muskingumTravelTime_impervious
  !>       - param(5) = muskingumAttenuation_riverSlope
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(5) :: param"  input parameter
  !>       \param[in] "real(dp), dimension(:) :: length" [m] total length
  !>       \param[in] "real(dp), dimension(:) :: slope"  average slope
  !>       \param[in] "real(dp), dimension(:) :: fFPimp" fraction of the flood plain with
  !>       impervious layer
  !>       \param[in] "real(dp) :: TS"                   - [h] time step in
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp), dimension(:) :: C1" routing parameter C1 (Chow, 25-41)
  !>       \param[out] "real(dp), dimension(:) :: C2" routing parameter C2 (")
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine reg_rout(param, length, slope, fFPimp, TS, C1, C2)
    implicit none
 
    ! input parameter
    real(dp), dimension(5), intent(in) :: param
 
    ! [m] total length
    real(dp), dimension(:), intent(in) :: length
 
    ! average slope
    real(dp), dimension(:), intent(in) :: slope
 
    ! fraction of the flood plain with
    ! impervious layer
    real(dp), dimension(:), intent(in) :: ffpimp
 
    ! - [h] time step in
    real(dp), intent(in) :: ts
 
    ! routing parameter C1 (Chow, 25-41)
    real(dp), dimension(:), intent(out) :: c1
 
    ! routing parameter C2 (")
    real(dp), dimension(:), intent(out) :: c2
 
    ! stream slope max
    real(dp) :: ssmax
 
    ! [d] Muskingum travel time parameter
    real(dp), dimension(size(fFPimp, 1)) :: k
 
    ! [1] Muskingum diffusion parameter (attenuation)
    real(dp), dimension(size(fFPimp, 1)) :: xi
 
 
    ! normalize stream bed slope
    ssmax = maxval(slope(:))
 
    ! New regional relationship; K = f(length, slope, & fFPimp)
    k = param(1) + param(2) * (length * 0.001_dp) &
            + param(3) * slope &
            + param(4) * ffpimp
 
    ! Xi = f(slope)
    xi = param(5) * (1.0_dp + slope / ssmax)
 
    ! constraints on Xi
    xi = merge(0.5_dp, xi, xi > 0.5_dp)
    xi = merge(0.005_dp, xi, xi < 0.005_dp)
 
    ! constrains on Ki
    k = merge(0.5_dp * ts / xi, k, k > 0.5_dp * ts / xi)
    k = merge(0.5_dp * ts / (1.0_dp - xi), k, k < 0.5_dp * ts / (1.0_dp - xi))
 
    ! Muskingum parameters
    c1 = ts / (k * (1.0_dp - xi) + 0.5_dp * ts)
    c2 = 1.0_dp - c1 * k / ts
 
  end subroutine reg_rout
 
  ! --------------------------------------------------------------------------
  ! L11 PARAMETERS
  ! --------------------------------------------------------------------------
  ! The parameters are set following Thober et al. 2017
  ! Modified:
  !    NAME
  !        mrm_init_param
 
  !    PURPOSE
  !>       \brief TODO: add description
 
  !>       \details TODO: add description
 
  !    INTENT(IN)
  !>       \param[in] "integer(i4) :: iDomain"           domain number
  !>       \param[in] "real(dp), dimension(:) :: param" input parameter (param(1) is celerity in m/s)
 
  !    HISTORY
  !>       \authors Stephan Thober, Matthias Kelbling
 
  !>       \date Jun 2018
 
  ! Modifications:
 
  subroutine mrm_init_param(iDomain, param)
 
    use mo_constants, only : hoursecs
    use mo_common_mhm_mrm_variables, only : resolutionrouting, timestep, optimize
    use mo_common_variables, only : iflag_cordinate_sys, domainmeta, processmatrix
    use mo_kind, only : dp, i4
    use mo_mrm_constants, only : given_ts
    use mo_mrm_global_variables, only : level11, l11_tsrout, domain_mrm, l11_celerity
    use mo_string_utils, only : num2str
    use mo_utils, only : locate, notequal
    use mo_mrm_net_startup, only : l11_calc_celerity
 
    implicit none
 
    ! domain number
    integer(i4), intent(in) :: idomain
 
    ! input parameter (param(1) is celerity in m/s)
    real(dp), dimension(:), intent(in) :: param
 
    ! index selected from given_TS
    integer(i4) :: index
 
    ! spatial routing resolution
    real(dp) :: deltax
 
    ! [s] wave travel time parameter
    real(dp) :: k
 
    ! start and end index at level11
    integer(i4) :: s11, e11
 
    ! Number of cells within
 
    ! initialize indices
    s11 = level11(idomain)%iStart
    e11 = level11(idomain)%iEnd
 
    ! temporal resolution of routing
    if (idomain .eq. 1 .and. .not. allocated(l11_tsrout)) then
      allocate(l11_tsrout(domainmeta%nDomains))
      l11_tsrout = 0._dp
    end if
 
    if (processmatrix(8, 1) .eq. 1) then
       l11_tsrout = timestep * hoursecs
 
       if ( notequal(mod(hoursecs * 24.0_dp, l11_tsrout(idomain)), 0.0_dp) .and. &
            (domain_mrm(idomain)%nInflowGauges .gt. 0)) then
          call message('***WARNING: routing timestep is not a multiple of 24 h.')
          call message('            Inflowgauge timeseries is averaged over values')
          call message('            of different days, small mismatches at')
          call message('            inflowgauge location might occur.')
       end if
 
    else
 
      ! called for initialization
      call mrm_update_param(idomain, param)
 
    end if
 
    call message('')
    call message('    domain: '//num2str(idomain, '(i3)'))
    call message('      routing resolution [s]:. '//num2str(l11_tsrout(idomain), '(f7.0)'))
    call message('      routing factor:......... '//num2str(l11_tsrout(idomain) / (timestep * hoursecs), '(f5.2)'))
 
    if ( notequal(mod(hoursecs * 24.0_dp, l11_tsrout(idomain)), 0.0_dp) .and. &
        (domain_mrm(idomain)%nInflowGauges .gt. 0)) then
       call message('***WARNING: routing timestep is not a multiple of 24 h.')
       call message('            Inflowgauge timeseries is averaged over values')
       call message('            of different days, small mismatches at')
       call message('            inflowgauge location might occur.')
    end if
 
  end subroutine mrm_init_param
 
  !    NAME
  !        mrm_update_param
 
  !    PURPOSE
  !>       \brief TODO: add description
 
  !>       \details TODO: add description
 
  !    INTENT(IN)
  !>       \param[in] "integer(i4) :: iDomain"           domain number
  !>       \param[in] "real(dp), dimension(1) :: param" celerity parameter [m s-1]
 
  !    HISTORY
  !>       \authors Stehpan Thober, Matthias Kelbling
 
  !>       \date Jun 2018
 
  ! Modifications:
 
  subroutine mrm_update_param(iDomain, param)
 
    use mo_kind, only: i4, dp
    use mo_common_variables, only: processmatrix, iflag_cordinate_sys
    use mo_mrm_global_variables, only: &
         ! input variable
         level11, &
         l11_tsrout, &
         l11_celerity, l11_noutlets, l11_length, &
         ! output variables
         l11_c1, l11_c2
    use mo_common_mhm_mrm_variables, only: resolutionrouting, optimize, timestep, &
         optimize
    use mo_mrm_constants, only: rout_space_weight, given_ts
    use mo_utils, only: locate
    use mo_mrm_net_startup, only: l11_calc_celerity
    use mo_mrm_constants, only: given_ts
    use mo_constants, only: hoursecs
    use mo_string_utils, only: num2str
    use mo_utils, only: locate
 
    implicit none
 
    ! domain number
    integer(i4), intent(in) :: idomain
 
    ! celerity parameter [m s-1]
    real(dp), intent(in), dimension(1) :: param
 
    integer(i4) :: s11
    integer(i4) :: e11
    integer(i4) :: nnodes
 
    ! index selected from given_TS
    integer(i4) :: ind
 
    ! spatial routing resolution
    real(dp) :: deltax
    real(dp), allocatable :: length(:)
 
    ! [s] wave travel time parameter
    real(dp), allocatable :: k(:)
 
    ! [1] Muskingum diffusion parameter (attenuation)
    real(dp) :: xi
 
    ! get domain information
    s11 = level11(idomain)%iStart
    e11 = level11(idomain)%iEnd
    nnodes = level11(idomain)%nCells
 
    allocate(k(nnodes))
 
    if (processmatrix(8, 1) .eq. 2_i4) then
 
      ! [s] wave travel time parameter
      k(:) = l11_length(s11: e11) / param(1)
 
    else if (processmatrix(8, 1) .eq. 3_i4) then
 
      ! [s] wave travel time parameter
      call l11_calc_celerity( idomain, param)
 
      ! Allocate and calculate K
      k(:) = l11_length(s11: e11) / l11_celerity(s11:e11)
 
    end if
 
    ! set time-weighting scheme
    xi = abs(rout_space_weight) ! set weighting factor to 0._dp
 
    ! determine routing timestep
    ind = locate(given_ts, minval(k(1:(nnodes-l11_noutlets(idomain)))))
 
    ! set min-wave traveltime to min given_TS
    if (ind .lt. 1) ind = 1
    l11_tsrout(idomain) = given_ts(ind)
 
    ! Muskingum parameters
    l11_c1(s11:e11) = l11_tsrout(idomain) / ( k(:) * (1.0_dp - xi) + 0.5_dp * l11_tsrout(idomain) )
    l11_c2(s11:e11) = 1.0_dp - l11_c1(s11:e11) * k(:) / l11_tsrout(idomain)
 
    deallocate(k)
 
    ! optional print
    ! print *, 'C1 Muskingum routing parameter: ', L11_C1(s11)
    ! print *, 'C2 Muskingum routing parameter: ', L11_C2(s11)
 
  end subroutine mrm_update_param
 
end module mo_mrm_mpr