LCOV - code coverage report
Current view: top level - mRM - mo_mrm_mpr.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 55 64 85.9 %
Date: 2024-04-30 08:53:32 Functions: 3 3 100.0 %

          Line data    Source code
       1             : !> \file mo_mrm_mpr.f90
       2             : !> \brief \copybrief mo_mrm_mpr
       3             : !> \details \copydetails mo_mrm_mpr
       4             : 
       5             : !> \brief Perform Multiscale Parameter Regionalization on Routing Parameters
       6             : !> \details This module contains the subroutine for calculating the regionalized
       7             : !! routing parameters (beta-parameters) given the five global routing parameters (gamma) at the level 0 scale.
       8             : !> \authors Luis Samaniego, Stephan Thober
       9             : !> \date Aug 2015
      10             : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
      11             : !! mHM is released under the LGPLv3+ license \license_note
      12             : !> \ingroup f_mrm
      13             : module mo_mrm_mpr
      14             :   use mo_kind, only : dp
      15             :   use mo_message, only : message
      16             : 
      17             :   implicit none
      18             : 
      19             :   public :: reg_rout
      20             :   public :: mrm_init_param
      21             :   public :: mrm_update_param
      22             :   private
      23             : contains
      24             : 
      25             :   ! ----------------------------------------------------------------------------
      26             : 
      27             :   !    NAME
      28             :   !        reg_rout
      29             : 
      30             :   !    PURPOSE
      31             :   !>       \brief Regionalized routing
      32             : 
      33             :   !>       \details sets up the Regionalized Routing parameters
      34             :   !>       Global parameters needed (see mhm_parameter.nml):
      35             :   !>       - param(1) = muskingumTravelTime_constant
      36             :   !>       - param(2) = muskingumTravelTime_riverLength
      37             :   !>       - param(3) = muskingumTravelTime_riverSlope
      38             :   !>       - param(4) = muskingumTravelTime_impervious
      39             :   !>       - param(5) = muskingumAttenuation_riverSlope
      40             : 
      41             :   !    INTENT(IN)
      42             :   !>       \param[in] "real(dp), dimension(5) :: param"  input parameter
      43             :   !>       \param[in] "real(dp), dimension(:) :: length" [m] total length
      44             :   !>       \param[in] "real(dp), dimension(:) :: slope"  average slope
      45             :   !>       \param[in] "real(dp), dimension(:) :: fFPimp" fraction of the flood plain with
      46             :   !>       impervious layer
      47             :   !>       \param[in] "real(dp) :: TS"                   - [h] time step in
      48             : 
      49             :   !    INTENT(OUT)
      50             :   !>       \param[out] "real(dp), dimension(:) :: C1" routing parameter C1 (Chow, 25-41)
      51             :   !>       \param[out] "real(dp), dimension(:) :: C2" routing parameter C2 (")
      52             : 
      53             :   !    HISTORY
      54             :   !>       \authors Stephan Thober, Rohini Kumar
      55             : 
      56             :   !>       \date Dec 2012
      57             : 
      58             :   ! Modifications:
      59             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
      60             : 
      61      990192 :   subroutine reg_rout(param, length, slope, fFPimp, TS, C1, C2)
      62             :     implicit none
      63             : 
      64             :     ! input parameter
      65             :     real(dp), dimension(5), intent(in) :: param
      66             : 
      67             :     ! [m] total length
      68             :     real(dp), dimension(:), intent(in) :: length
      69             : 
      70             :     ! average slope
      71             :     real(dp), dimension(:), intent(in) :: slope
      72             : 
      73             :     ! fraction of the flood plain with
      74             :     ! impervious layer
      75             :     real(dp), dimension(:), intent(in) :: fFPimp
      76             : 
      77             :     ! - [h] time step in
      78             :     real(dp), intent(in) :: TS
      79             : 
      80             :     ! routing parameter C1 (Chow, 25-41)
      81             :     real(dp), dimension(:), intent(out) :: C1
      82             : 
      83             :     ! routing parameter C2 (")
      84             :     real(dp), dimension(:), intent(out) :: C2
      85             : 
      86             :     ! stream slope max
      87    23151264 :     real(dp) :: ssMax
      88             : 
      89             :     ! [d] Muskingum travel time parameter
      90    23646360 :     real(dp), dimension(size(fFPimp, 1)) :: K
      91             : 
      92             :     ! [1] Muskingum diffusion parameter (attenuation)
      93    23646360 :     real(dp), dimension(size(fFPimp, 1)) :: xi
      94             : 
      95             : 
      96             :     ! normalize stream bed slope
      97    23646360 :     ssMax = maxval(slope(:))
      98             : 
      99             :     ! New regional relationship; K = f(length, slope, & fFPimp)
     100      495096 :     K = param(1) + param(2) * (length * 0.001_dp) &
     101           0 :             + param(3) * slope &
     102    23646360 :             + param(4) * fFPimp
     103             : 
     104             :     ! Xi = f(slope)
     105    23151264 :     xi = param(5) * (1.0_dp + slope / ssMax)
     106             : 
     107             :     ! constraints on Xi
     108    23151264 :     xi = merge(0.5_dp, xi, xi > 0.5_dp)
     109    23151264 :     xi = merge(0.005_dp, xi, xi < 0.005_dp)
     110             : 
     111             :     ! constrains on Ki
     112    23151264 :     K = merge(0.5_dp * TS / xi, K, K > 0.5_dp * TS / xi)
     113    23151264 :     K = merge(0.5_dp * TS / (1.0_dp - xi), K, K < 0.5_dp * TS / (1.0_dp - xi))
     114             : 
     115             :     ! Muskingum parameters
     116    23646360 :     C1 = TS / (K * (1.0_dp - xi) + 0.5_dp * TS)
     117    23646360 :     C2 = 1.0_dp - C1 * K / TS
     118             : 
     119      495096 :   end subroutine reg_rout
     120             : 
     121             :   ! --------------------------------------------------------------------------
     122             :   ! L11 PARAMETERS
     123             :   ! --------------------------------------------------------------------------
     124             :   ! The parameters are set following Thober et al. 2017
     125             :   ! Modified:
     126             :   !    NAME
     127             :   !        mrm_init_param
     128             : 
     129             :   !    PURPOSE
     130             :   !>       \brief TODO: add description
     131             : 
     132             :   !>       \details TODO: add description
     133             : 
     134             :   !    INTENT(IN)
     135             :   !>       \param[in] "integer(i4) :: iDomain"           domain number
     136             :   !>       \param[in] "real(dp), dimension(:) :: param" input parameter (param(1) is celerity in m/s)
     137             : 
     138             :   !    HISTORY
     139             :   !>       \authors Stephan Thober, Matthias Kelbling
     140             : 
     141             :   !>       \date Jun 2018
     142             : 
     143             :   ! Modifications:
     144             : 
     145          25 :   subroutine mrm_init_param(iDomain, param)
     146             : 
     147      495096 :     use mo_constants, only : HourSecs
     148             :     use mo_common_mHM_mRM_variables, only : resolutionRouting, timeStep, optimize
     149             :     use mo_common_variables, only : iFlag_cordinate_sys, domainMeta, processMatrix
     150             :     use mo_kind, only : dp, i4
     151             :     use mo_mrm_constants, only : given_TS
     152             :     use mo_mrm_global_variables, only : level11, L11_tsRout, domain_mrm, L11_celerity
     153             :     use mo_string_utils, only : num2str
     154             :     use mo_utils, only : locate, notequal
     155             :     use mo_mrm_net_startup, only : L11_calc_celerity
     156             : 
     157             :     implicit none
     158             : 
     159             :     ! domain number
     160             :     integer(i4), intent(in) :: iDomain
     161             : 
     162             :     ! input parameter (param(1) is celerity in m/s)
     163             :     real(dp), dimension(:), intent(in) :: param
     164             : 
     165             :     ! index selected from given_TS
     166             :     integer(i4) :: index
     167             : 
     168             :     ! spatial routing resolution
     169             :     real(dp) :: deltaX
     170             : 
     171             :     ! [s] wave travel time parameter
     172             :     real(dp) :: K
     173             : 
     174             :     ! start and end index at level11
     175             :     integer(i4) :: s11, e11
     176             : 
     177             :     ! Number of cells within
     178             : 
     179             :     ! initialize indices
     180          25 :     s11 = level11(iDomain)%iStart
     181          25 :     e11 = level11(iDomain)%iEnd
     182             : 
     183             :     ! temporal resolution of routing
     184          25 :     if (iDomain .eq. 1 .and. .not. allocated(L11_tsRout)) then
     185          36 :       allocate(L11_tsRout(domainMeta%nDomains))
     186          36 :       L11_TSrout = 0._dp
     187             :     end if
     188             : 
     189          25 :     if (processMatrix(8, 1) .eq. 1) then
     190          68 :        L11_tsRout = timestep * HourSecs
     191             : 
     192          32 :        if ( NOTEQUAL(mod(HourSecs * 24.0_dp, L11_tsRout(iDomain)), 0.0_dp) .and. &
     193          16 :             (domain_mrm(iDomain)%nInflowGauges .gt. 0)) then
     194           0 :           call message('***WARNING: routing timestep is not a multiple of 24 h.')
     195           0 :           call message('            Inflowgauge timeseries is averaged over values')
     196           0 :           call message('            of different days, small mismatches at')
     197           0 :           call message('            inflowgauge location might occur.')
     198             :        end if
     199             : 
     200             :     else
     201             : 
     202             :       ! called for initialization
     203           9 :       call mrm_update_param(iDomain, param)
     204             : 
     205             :     end if
     206             : 
     207          25 :     call message('')
     208          25 :     call message('    domain: '//num2str(iDomain, '(i3)'))
     209          25 :     call message('      routing resolution [s]:. '//num2str(L11_tsRout(iDomain), '(f7.0)'))
     210          25 :     call message('      routing factor:......... '//num2str(L11_tsRout(iDomain) / (timestep * HourSecs), '(f5.2)'))
     211             : 
     212          50 :     if ( NOTEQUAL(mod(HourSecs * 24.0_dp, L11_tsRout(iDomain)), 0.0_dp) .and. &
     213          25 :         (domain_mrm(iDomain)%nInflowGauges .gt. 0)) then
     214           0 :        call message('***WARNING: routing timestep is not a multiple of 24 h.')
     215           0 :        call message('            Inflowgauge timeseries is averaged over values')
     216           0 :        call message('            of different days, small mismatches at')
     217           0 :        call message('            inflowgauge location might occur.')
     218             :     end if
     219             : 
     220          25 :   end subroutine mrm_init_param
     221             : 
     222             :   !    NAME
     223             :   !        mrm_update_param
     224             : 
     225             :   !    PURPOSE
     226             :   !>       \brief TODO: add description
     227             : 
     228             :   !>       \details TODO: add description
     229             : 
     230             :   !    INTENT(IN)
     231             :   !>       \param[in] "integer(i4) :: iDomain"           domain number
     232             :   !>       \param[in] "real(dp), dimension(1) :: param" celerity parameter [m s-1]
     233             : 
     234             :   !    HISTORY
     235             :   !>       \authors Stehpan Thober, Matthias Kelbling
     236             : 
     237             :   !>       \date Jun 2018
     238             : 
     239             :   ! Modifications:
     240             : 
     241          19 :   subroutine mrm_update_param(iDomain, param)
     242             : 
     243          25 :     use mo_kind, only: i4, dp
     244             :     use mo_common_variables, only: processMatrix, iFlag_cordinate_sys
     245             :     use mo_mrm_global_variables, only: &
     246             :          ! input variable
     247             :          level11, &
     248             :          L11_TSrout, &
     249             :          L11_celerity, L11_nOutlets, L11_length, &
     250             :          ! output variables
     251             :          L11_C1, L11_C2
     252             :     use mo_common_mHM_mRM_variables, only: resolutionRouting, optimize, timeStep, &
     253             :          optimize
     254             :     use mo_mrm_constants, only: rout_space_weight, given_TS
     255             :     use mo_utils, only: locate
     256             :     use mo_mrm_net_startup, only: L11_calc_celerity
     257             :     use mo_mrm_constants, only: given_TS
     258             :     use mo_constants, only: HourSecs
     259             :     use mo_string_utils, only: num2str
     260             :     use mo_utils, only: locate
     261             : 
     262             :     implicit none
     263             : 
     264             :     ! domain number
     265             :     integer(i4), intent(in) :: iDomain
     266             : 
     267             :     ! celerity parameter [m s-1]
     268             :     real(dp), intent(in), dimension(1) :: param
     269             : 
     270             :     integer(i4) :: s11
     271             :     integer(i4) :: e11
     272             :     integer(i4) :: nNodes
     273             : 
     274             :     ! index selected from given_TS
     275             :     integer(i4) :: ind
     276             : 
     277             :     ! spatial routing resolution
     278             :     real(dp) :: deltaX
     279             :     real(dp), allocatable :: length(:)
     280             : 
     281             :     ! [s] wave travel time parameter
     282          19 :     real(dp), allocatable :: K(:)
     283             : 
     284             :     ! [1] Muskingum diffusion parameter (attenuation)
     285          19 :     real(dp) :: xi
     286             : 
     287             :     ! get domain information
     288          19 :     s11 = level11(iDomain)%iStart
     289          19 :     e11 = level11(iDomain)%iEnd
     290          19 :     Nnodes = level11(iDomain)%nCells
     291             : 
     292          57 :     allocate(K(nNodes))
     293             : 
     294          19 :     if (ProcessMatrix(8, 1) .eq. 2_i4) then
     295             : 
     296             :       ! [s] wave travel time parameter
     297          70 :       K(:) = L11_length(s11: e11) / param(1)
     298             : 
     299          17 :     else if (ProcessMatrix(8, 1) .eq. 3_i4) then
     300             : 
     301             :       ! [s] wave travel time parameter
     302          17 :       call L11_calc_celerity( iDomain, param)
     303             : 
     304             :       ! Allocate and calculate K
     305         595 :       K(:) = L11_length(s11: e11) / L11_celerity(s11:e11)
     306             : 
     307             :     end if
     308             : 
     309             :     ! set time-weighting scheme
     310          19 :     xi = abs(rout_space_weight) ! set weighting factor to 0._dp
     311             : 
     312             :     ! determine routing timestep
     313         684 :     ind = locate(given_TS, minval(K(1:(nNodes-L11_nOutlets(iDomain)))))
     314             : 
     315             :     ! set min-wave traveltime to min given_TS
     316          19 :     if (ind .lt. 1) ind = 1
     317          19 :     L11_TSrout(iDomain) = given_TS(ind)
     318             : 
     319             :     ! Muskingum parameters
     320         665 :     L11_C1(s11:e11) = L11_TSrout(iDomain) / ( K(:) * (1.0_dp - xi) + 0.5_dp * L11_TSrout(iDomain) )
     321         665 :     L11_C2(s11:e11) = 1.0_dp - L11_C1(s11:e11) * K(:) / L11_TSrout(iDomain)
     322             : 
     323          19 :     deallocate(K)
     324             : 
     325             :     ! optional print
     326             :     ! print *, 'C1 Muskingum routing parameter: ', L11_C1(s11)
     327             :     ! print *, 'C2 Muskingum routing parameter: ', L11_C2(s11)
     328             : 
     329          19 :   end subroutine mrm_update_param
     330             : 
     331             : end module mo_mrm_mpr

Generated by: LCOV version 1.16