mo_multi_param_reg.f90

Go to the documentation of this file.

!> \file mo_multi_param_reg.f90
!> \brief \copybrief mo_multi_param_reg
!> \details \copydetails mo_multi_param_reg
 
!> \brief Multiscale parameter regionalization (MPR).
!> \details This module provides the routines for multiscale parameter regionalization (MPR).
!> \authors Stephan Thober, Rohini Kumar
!> \date Dec 2012
!> \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_mpr
MODULE mo_multi_param_reg
 
  use mo_kind, only : i4, dp
  use mo_common_constants, only : nodata_dp, nodata_i4
  use mo_message, only : message, error_message
 
  implicit none
 
  private
 
  PUBLIC :: mpr                     ! calculates effective regionalised parameters
  PUBLIC :: canopy_intercept_param  ! estimate effective max. canopy interception
 
contains
 
  !> \brief Regionalizing and Upscaling process parameters
  !> \details calculating process parameters at L0 scale (Regionalization), like:
  !! - Baseflow recession parameter
  !! - Soil moisture parameters
  !! - PET correction for aspect
  !!
  !! and upscale these parameters to retrieve effective parameters at scale L1.
  !! Further parameter regionalizations are done for:
  !! - snow accumulation and melting parameters
  !! - threshold parameter for runoff generation on impervious layer
  !! - karstic percolation loss
  !! - setting up the Regionalized Routing Parameters
  !!
  !> \changelog
  !! - Stephan Thober           Jan 2013
  !!   - updated calling sequence for upscaling operators
  !! - Luis Samaniego           Feb 2013
  !!   - calling sequence, initial CHECK, call mpr_runoff
  !! - Stephan Thober           Feb 2013
  !!   - added subroutine for karstic percolation loss removed L1_, L0_ in variable names
  !! - Stephan Thober           Aug 2015
  !!   - moved regionalization of routing to mRM
  !! - Rohini Kumar             Mar 2016
  !!   - changes for handling multiple soil database options
  !! - Zink M. & Demirel M.C.   Mar 2017
  !!   - Added Jarvis soil water stress function at SM process(3)
  !! - Demirel M.C. & S. Stisen Apr 2017
  !!   - Added FC dependency on root fraction coefficient at SM process(3)
  !! - Robert Schweppe          Dec 2017
  !!   - added loop over LCscenes inside MPR, renamed variables rewrite
  !! - Robert Schweppe          Jun 2018
  !!   - refactoring and reformatting
  !! - Demirel M.C. & S. Stisen Jun 2020
  !!   - Added Feddes and global FC dependency on root fraction coefficient at SM process(3)=4
  !! - Rohini Kumar             Oct 2021
  !!   - Added Neutron count module to mHM integrate into develop branch (5.11.2)
  !! - Sebastian Müller         Mar 2023
  !!   - made L1_alpha, L1_kSlowFlow, L1_kBaseFlow and L1_kPerco land cover dependent
  !> \authors Stephan Thober, Rohini Kumar
  !> \date Dec 2012
subroutine mpr(mask0, geoUnit0, soilId0, Asp0, gridded_LAI0, LCover0, slope_emp0, y0, Id0, upper_bound1, lower_bound1, &
       left_bound1, right_bound1, n_subcells1, fSealed1, alpha1, degDayInc1, degDayMax1, degDayNoPre1, fAsp1, &
       HarSamCoeff1, PrieTayAlpha1, aeroResist1, surfResist1, fRoots1, kFastFlow1, kSlowFlow1, kBaseFlow1, &
       kPerco1, karstLoss1, soilMoistFC1, soilMoistSat1, soilMoistExp1, jarvis_thresh_c1, tempThresh1, &
       unsatThresh1, sealedThresh1, wiltingPoint1, maxInter1, petLAIcorFactor, &
       No_Count1, bulkDens1, latticeWater1, COSMICL31, &
       parameterset )
 
    use mo_common_variables, only : global_parameters, processmatrix
    use mo_mpr_smhorizons, only : mpr_smhorizons
    use mo_mpr_global_variables, only : horizondepth_mhm, fracsealed_cityarea, iflag_soildb, nsoilhorizons_mhm, &
         soildb
    use mo_mpr_pet, only : bulksurface_resistance, pet_correctbyasp, pet_correctbylai, priestley_taylor_alpha
    use mo_mpr_runoff, only : mpr_runoff
    use mo_mpr_soilmoist, only : mpr_sm
    use mo_upscaling_operators, only : l0_fractionalcover_in_lx, &
         upscale_arithmetic_mean
    use mo_mpr_neutrons,        only: mpr_neutrons
    implicit none
 
    !> mask at level 0 field
    logical, dimension(:, :), intent(in) :: mask0
    !> L0 geological units
    integer(i4), dimension(:), intent(in) :: geounit0
    !> soil Ids at level 0
    integer(i4), dimension(:, :), intent(in) :: soilid0
    !> [degree] Aspect at Level 0
    real(dp), dimension(:), intent(in) :: asp0
    !> LAI grid at level 0, with dim2 = time
    real(dp), dimension(:, :), intent(in) :: gridded_lai0
    !> land cover at level 0
    integer(i4), dimension(:, :), intent(in) :: lcover0
    !> Empirical quantiles of slope
    real(dp), dimension(:), intent(in) :: slope_emp0
    !> Cell ids at level 0
    integer(i4), dimension(:), intent(in) :: id0
    !> Upper row of hi res block
    integer(i4), dimension(:), intent(in) :: upper_bound1
    !> Lower row of hi res block
    integer(i4), dimension(:), intent(in) :: lower_bound1
    !> Left column of hi res block
    integer(i4), dimension(:), intent(in) :: left_bound1
    !> Right column of hi res block
    integer(i4), dimension(:), intent(in) :: right_bound1
    !> Number of L0 cells within a L1 cell
    integer(i4), dimension(:), intent(in) :: n_subcells1
    !> y0 at level 0
    real(dp), dimension(:), intent(in) :: y0
    !> [1] fraction of sealed area
    real(dp), dimension(:, :, :), intent(inout) :: fsealed1
    !> Parameter that determines the rel. contribution to SM, upscal. Bulk den.
    real(dp), dimension(:, :, :), intent(inout) :: soilmoistexp1
    !> [1] jarvis critical value for norm SWC
    real(dp), dimension(:, :, :), intent(inout) :: jarvis_thresh_c1
    !> [10^-3 m] depth of saturated SM
    real(dp), dimension(:, :, :), intent(inout) :: soilmoistsat1
    !> [10^-3 m] field capacity
    real(dp), dimension(:, :, :), intent(inout) :: soilmoistfc1
    !> [10^-3 m] permanent wilting point
    real(dp), dimension(:, :, :), intent(inout) :: wiltingpoint1
    !> fraction of roots in soil horizon
    real(dp), dimension(:, :, :), intent(inout) :: froots1
    !> [degreeC] threshold temperature for snow rain
    real(dp), dimension(:, :, :), intent(inout) :: tempthresh1
    !> [mm-1 degreeC-1] Degree-day factor with no precipitation
    real(dp), dimension(:, :, :), intent(inout) :: degdaynopre1
    !> [mm-1 degreeC-1] Maximum Degree-day factor
    real(dp), dimension(:, :, :), intent(inout) :: degdaymax1
    !> [d-1 degreeC-1]  Increase of the Degree-day factor per mm of increase in precipitation
    real(dp), dimension(:, :, :), intent(inout) :: degdayinc1
    !> [1]     PET correction for Aspect at level 1
    real(dp), dimension(:, :, :), intent(inout) :: fasp1
    !> [1]     PET Hargreaves Samani coeff. at level 1
    real(dp), dimension(:, :, :), intent(inout) :: harsamcoeff1
    !> [1]     PET Priestley Taylor coeff. at level 1
    real(dp), dimension(:, :, :), intent(inout) :: prietayalpha1
    !> [s m-1] PET aerodynamical resitance at level 1
    real(dp), dimension(:, :, :), intent(inout) :: aeroresist1
    !> [s m-1] PET bulk surface resitance at level 1
    real(dp), dimension(:, :, :), intent(inout) :: surfresist1
    !> threshold parameter
    real(dp), dimension(:, :, :), intent(inout) :: sealedthresh1
    !> [10^-3 m] Threshhold water depth in upper reservoir (for Runoff contribution)
    real(dp), dimension(:, :, :), intent(inout) :: unsatthresh1
    !> [10^-3 m] Recession coefficient of the upper reservoir, upper outlet
    real(dp), dimension(:, :, :), intent(inout) :: kfastflow1
    !> [10^-3 m] Recession coefficient of the upper reservoir, lower outlet
    real(dp), dimension(:, :, :), intent(inout) :: kslowflow1
    !> Level 1 baseflow recession
    real(dp), dimension(:, :, :), intent(inout) :: kbaseflow1
    !> [1] Exponent for the upper reservoir
    real(dp), dimension(:, :, :), intent(inout) :: alpha1
    !> [d-1] percolation coefficient
    real(dp), dimension(:, :, :), intent(inout) :: kperco1
    !> karstic percolation loss
    real(dp), dimension(:, :, :), intent(inout) :: karstloss1
    !> max interception
    real(dp), dimension(:, :, :), intent(inout) :: maxinter1
    !> pet cor factor at level-1
    real(dp), dimension(:, :, :), intent(inout) :: petlaicorfactor
    !> [-] inital neutron count
    real(dp), dimension(:, :, :), intent(inout) :: no_count1
    !> [gcm-3] bulk density
    real(dp), dimension(:, :, :), intent(inout) :: bulkdens1
    !> [mm/mm] lattice water content
    real(dp), dimension(:, :, :), intent(inout) :: latticewater1
    !> [-] cosmic L3 parameter
    real(dp), dimension(:, :, :), intent(inout) :: cosmicl31
 
    !> array of global parameters
    real(dp), dimension(:), intent(in), optional, target :: parameterset
 
    ! array of global parameters
    real(dp), dimension(:), pointer :: param
 
    real(dp), dimension(:, :, :), allocatable :: thetas_till
 
    real(dp), dimension(:, :, :), allocatable :: thetafc_till
 
    real(dp), dimension(:, :, :), allocatable :: thetapw_till
 
    ! saturated hydraulic conductivity
    real(dp), dimension(:, :, :), allocatable :: ks
 
    ! Bulk density
    real(dp), dimension(:, :, :), allocatable :: db
    real(dp), dimension(:, :), allocatable :: thetas
    real(dp), dimension(:, :), allocatable :: thetafc
    real(dp), dimension(:, :), allocatable :: thetapw
 
    ! neutron count
    real(dp), dimension(:,:,:), allocatable :: latwat_till
    real(dp), dimension(:,:,:), allocatable :: cosmic_l3_till
    real(dp), dimension(:,:), allocatable   :: latwat         ! lattice water
    real(dp), dimension(:,:), allocatable   :: cosmic_l3      ! COSMIC parameter L3
 
 
    ! relative variability of saturated
    ! hydraulic cound. for Horizantal flow
    real(dp), dimension(:), allocatable :: ksvar_h0
 
    ! relative variability of saturated
    ! hydraulic cound. for vertical flow
    real(dp), dimension(:), allocatable :: ksvar_v0
 
    ! soil mositure deficit from
    ! field cap. w.r.t to saturation
    real(dp), dimension(:), allocatable :: sms_fc0
 
    ! L0 baseflow parameter
    real(dp), dimension(size(Id0, 1)) :: k2_0
 
    ! L1 baseflow parameter
    real(dp), dimension(:), allocatable :: k2_1
 
    ! L0 Aspect
    real(dp), dimension(size(Id0, 1)) :: fasp0
 
    ! number of soil classes
    integer(i4) :: msoil
 
    ! maximum of number of Tillage horizons
    integer(i4) :: mtill
 
    ! maximum number of horizons
    integer(i4) :: mhor
 
    ! number of Landcover classes
    integer(i4) :: mlc
 
    ! indexing of parameter vector - start
    integer(i4) :: istart
 
    ! indexing of parameter vector - end
    integer(i4) :: iend
 
    ! 2nd indexing of parameter vector - start
    integer(i4) :: istart2
 
    ! 2nd indexing of parameter vector - end
    integer(i4) :: iend2
 
    ! counter for looping over LCscenes
    integer(i4) :: iilc
 
    ! [1]  Fraction of forest cover
    real(dp), dimension(size(fSealed1, dim = 1)) :: fforest1
 
    ! [1]  Fraction of permeable cover
    real(dp), dimension(size(fSealed1, dim = 1)) :: fperm1
 
 
    if (present(parameterset)) then
       param => parameterset
    else
       param => global_parameters(:, 3)
    end if
 
 
    ! loop over all LCover scenes
    do iilc = 1, size(lcover0, 2)
 
       ! estimate land cover fractions for dominant landcover class
       ! fSealed is intent inout, the rest only intent in
       fforest1(:) = l0_fractionalcover_in_lx(lcover0(:, iilc), 1, mask0, &
            upper_bound1, &
            lower_bound1, &
            left_bound1, &
            right_bound1, &
            n_subcells1)
       fsealed1(:, 1, iilc) = l0_fractionalcover_in_lx(lcover0(:, iilc), 2, mask0, &
            upper_bound1, &
            lower_bound1, &
            left_bound1, &
            right_bound1, &
            n_subcells1)
       fperm1(:) = l0_fractionalcover_in_lx(lcover0(:, iilc), 3, mask0, &
            upper_bound1, &
            lower_bound1, &
            left_bound1, &
            right_bound1, &
            n_subcells1)
       !---------------------------------------------------------
       ! Update fractions of sealed area fractions
       ! based on the sealing fraction[0-1] in cities
       !---------------------------------------------------------
       ! a factor is applied to the sealed area, effectively reducing it
       fsealed1(:, 1, iilc) = fracsealed_cityarea * fsealed1(:, 1, iilc)
       ! the forest area is kept constant, but the permeable area is increased so that the
       ! sum off all fractions equals 1 again
       fperm1(:) = 1.0_dp - fsealed1(:, 1, iilc) - fforest1(:)
 
       ! ------------------------------------------------------------------
       ! snow parameters
       ! ------------------------------------------------------------------
       select case(processmatrix(2,1))
       case(1)
 
          istart = processmatrix(2, 3) - processmatrix(2, 2) + 1
          iend = processmatrix(2, 3)
 
          call snow_acc_melt_param(param(istart : iend), & ! intent(in)
               fforest1, fsealed1(:, 1, iilc), fperm1, & ! intent(in)
               tempthresh1(:, 1, iilc), degdaynopre1(:, 1, iilc), & ! intent(out)
               degdayinc1(:, 1, iilc), degdaymax1(:, 1, iilc) & ! intent(out)
               )
       case DEFAULT
          call error_message('***ERROR: Process description for process "snow pack" does not exist! mo_multi_param_reg')
       end select
 
       ! ------------------------------------------------------------------
       ! Soil moisture parametrization
       ! ------------------------------------------------------------------
       msoil = size(soildb%is_present, 1)
       mlc = maxval(lcover0(:, iilc), (lcover0(:, iilc) .ne. nodata_i4))
 
       ! depending on which kind of soil database processing is to be performed
       if(iflag_soildb .eq. 0)then
          mtill = maxval(soildb%nTillHorizons, (soildb%nTillHorizons .ne. nodata_i4))
          mhor  = maxval(soildb%nHorizons,     (soildb%nHorizons     .ne. nodata_i4))
       else if(iflag_soildb .eq. 1) then
          ! here for each soil type both till and non-till soil hydraulic properties are to be estimated
          ! since a given soil type can lie in any horizon (till or non-till ones)
          ! adopt it in a way that it do not break the consistency of iFlag_soilDB = 0
          ! ** NOTE: SDB_nTillHorizons and SDB_nHorizons are also assigned in
          !          this flag option (see mo_soildatabase.f90 file - read_soil_LUT).
          !          But we are not using those variables here since in this case we have not
          !          varying number of soil horizons or either tillage horizons.
          !          So assigning them with a value = 1 is more than enough.
          mtill = 1
          mhor  = 1
       end if
 
       allocate(thetas_till(msoil, mtill,  mlc))
       allocate(thetafc_till(msoil, mtill, mlc))
       allocate(thetapw_till(msoil, mtill, mlc))
       allocate(thetas(msoil,  mhor))
       allocate(thetafc(msoil, mhor))
       allocate(thetapw(msoil, mhor))
       allocate(ks(msoil, mhor, mlc))
       allocate(db(msoil, mhor, mlc))
 
       ! neutron count related ones
       ! allocate and initalize here
       allocate(   latwat_till(msoil, mtill, mlc))
       allocate(cosmic_l3_till(msoil, mtill, mlc))
       allocate(        latwat(msoil, mhor      ))
       allocate(     cosmic_l3(msoil, mhor      ))
       latwat_till    = 0.000001_dp
       cosmic_l3_till = 0.000001_dp
       cosmic_l3      = 0.000001_dp
       latwat         = 0.000001_dp
 
 
 
       ! earlier these variables were allocated with  size(soilId0,1)
       ! in which the variable "soilId0" changes according to the iFlag_soilDB
       ! so better to use other variable which is common to both soilDB (0 AND 1) flags
       allocate(ksvar_h0(size(id0, 1)))
       allocate(ksvar_v0(size(id0, 1)))
       allocate( sms_fc0(size(id0, 1)))
 
       select case(processmatrix(3,1))
       case(1)
          ! first thirteen parameters go to this routine
          istart = processmatrix(3, 3) - processmatrix(3, 2) + 1
          iend = processmatrix(3, 3) - 4
 
          ! next four parameters go here
          ! (the first three for the fRoots and the fourth one for the beta)
          istart2 = processmatrix(3, 3) - 4 + 1
          iend2 = processmatrix(3, 3)
 
       case(2)
          ! first thirteen parameters go to this routine
          istart = processmatrix(3, 3) - processmatrix(3, 2) + 1
          iend = processmatrix(3, 3) - 5
 
          ! next four parameters go here
          ! (the first three for the fRoots and the fourth one for the beta)
          istart2 = processmatrix(3, 3) - 5 + 1
          iend2 = processmatrix(3, 3) - 1
          ! last parameter is jarvis parameter - no need to be regionalized
          jarvis_thresh_c1 = param(processmatrix(3, 3))
          !write(*,*) 'jarvis_thresh_c1 = ', jarvis_thresh_c1
          !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
 
       case(3)
          ! first thirteen parameters go to this routine
          istart = processmatrix(3, 3) - processmatrix(3, 2) + 1
          iend = processmatrix(3, 3) - 9
 
          ! next four parameters go here
          ! (the first three for the fRoots and the fourth one for the beta)
          istart2 = processmatrix(3, 3) - 8
          iend2 = processmatrix(3, 3) - 1
 
          ! last parameter is jarvis parameter - no need to be regionalized
          jarvis_thresh_c1 = param(processmatrix(3, 3))
 
          !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
          !write(*,*) 'jarvis_thresh_c1 = ', jarvis_thresh_c1
 
       case(4)
          ! first thirteen parameters go to this routine
          istart = processmatrix(3, 3) - processmatrix(3, 2) + 1
          iend = processmatrix(3, 3) - 8
 
          ! next four parameters go here
          ! (the first three for the fRoots and the fourth one for the beta)
          istart2 = processmatrix(3, 3) - 7
          iend2 = processmatrix(3, 3)
          !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
 
       case DEFAULT
          call error_message('***ERROR: Process description for process "soil moisture parametrization"', &
               'does not exist! mo_multi_param_reg')
       end select
 
       call mpr_sm(param(istart : iend), processmatrix, &
            soildb%is_present, soildb%nHorizons, soildb%nTillHorizons, &
            soildb%sand, soildb%clay, soildb%DbM, &
            id0, soilid0, lcover0(:, iilc), &
            thetas_till, thetafc_till, thetapw_till, thetas, &
            thetafc, thetapw, ks, db, ksvar_h0, ksvar_v0, sms_fc0)
 
       ! >> neutron count related parameters
       if ( processmatrix(10,1) .GT. 0 ) &
            call mpr_neutrons( processmatrix(10,1), &  ! IN: processmatrix case
            param( processmatrix(10,3)-processmatrix(10,2)+1:processmatrix(10,3) ) , & ! IN:  global parameter set
            soildb%is_present       , & ! IN:  flag indicating presence of soil
            soildb%nHorizons        , & ! IN:  Number of Horizons of Soiltype
            soildb%nTillHorizons    , & ! IN:  Number of tillage Horizons
            lcover0(:, iilc)        , & ! IN:  land cover ids at level 0
            soildb%clay             , & ! IN:  clay content
            soildb%DbM              , & ! IN:  mineral Bulk density
            db                      , & ! IN: Bulk density
            cosmic_l3_till          , & ! OUT: COSMIC parameter L3 tillage layer
            latwat_till             , & ! OUT: COSMIC parameter Lattice Water tillage layer
            cosmic_l3               , & ! OUT: COSMIC parameter L3
            latwat                    & ! OUT: COSMIC parameter Lattice Water
            )
 
       call mpr_smhorizons(param(istart2:iend2), processmatrix, &
            iflag_soildb, nsoilhorizons_mhm, horizondepth_mhm, &
            lcover0(:, iilc), soilid0, &
            soildb%nHorizons, soildb%nTillHorizons, &
            thetas_till, thetafc_till, thetapw_till, &
            thetas, thetafc, thetapw, &
            soildb%Wd, db, soildb%DbM, soildb%RZdepth, &
            mask0, id0, &
            upper_bound1, lower_bound1, left_bound1, right_bound1, n_subcells1, &
            soilmoistexp1(:, :, iilc), soilmoistsat1(:, :, iilc), soilmoistfc1(:, :, iilc), &
            wiltingpoint1(:, :, iilc), froots1(:, :, iilc), &
            ! >>>>>> neutron count
            latwat_till, cosmic_l3_till, latwat, cosmic_l3, &
            bulkdens1(:,:,iilc), latticewater1(:,:,iilc), cosmicl31(:,:,iilc) &
            )
 
       deallocate(thetas_till)
       deallocate(thetafc_till)
       deallocate(thetapw_till)
       deallocate(thetas)
       deallocate(thetafc)
       deallocate(thetapw)
       deallocate(ks)
       deallocate(db)
 
       ! neutron count
       deallocate( latwat_till    )
       deallocate( cosmic_l3_till )
       deallocate( latwat     )
       deallocate( cosmic_l3  )
 
       ! ------------------------------------------------------------------
       ! potential evapotranspiration (PET)
       ! ------------------------------------------------------------------
       ! Penman-Monteith method is only method that is LCscene dependent
       if (processmatrix(5, 1) == 3) then
          istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
          iend = processmatrix(5, 3)
          call aerodynamical_resistance(gridded_lai0, lcover0(:, iilc), param(istart : iend - 1), mask0, &
               id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
               aeroresist1(:, :, iilc))
       else if (processmatrix(5, 1) == -1) then
          istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
          iend = processmatrix(5, 3)
 
          call pet_correctbylai(param(istart : iend), nodata_dp, &
               lcover0(:, iilc), gridded_lai0, mask0, id0, &
               upper_bound1, lower_bound1, left_bound1, &
               right_bound1, n_subcells1, petlaicorfactor(:, :, iilc))
       end if
 
       ! ------------------------------------------------------------------
       ! interflow
       ! ------------------------------------------------------------------
       select case(processmatrix(6, 1))
       case (1)
          !
          istart = processmatrix(6, 3) - processmatrix(6, 2) + 1
          iend = processmatrix(6, 3)
          ! TODO: this subroutine should be split into each param (or at least extract kFastFlow1)
          ! because it is in the loop unnecessarily
          call mpr_runoff(lcover0(:, iilc), mask0, sms_fc0, slope_emp0, &
               ksvar_h0, param(istart : iend), id0, upper_bound1, lower_bound1, &
               left_bound1, right_bound1, n_subcells1, unsatthresh1(:, 1, 1), kfastflow1(:, 1, iilc), &
               kslowflow1(:, 1, iilc), alpha1(:, 1, iilc))
       case DEFAULT
          call error_message('***ERROR: Process description for process "interflow" does not exist! mo_multi_param_reg')
       END select
 
       ! ------------------------------------------------------------------
       ! percolation cofficient, karstic percolation loss
       ! ------------------------------------------------------------------
       select case(processmatrix(7, 1))
       case(1)
 
          istart = processmatrix(7, 3) - processmatrix(7, 2) + 1
          iend = processmatrix(7, 3)
          call karstic_layer(& ! In
               param(istart : iend), & ! In
               geounit0, mask0, & ! In
               sms_fc0, ksvar_v0, id0, & ! In
               n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, & ! In
               karstloss1(:, 1, 1), kperco1(:, 1, iilc)                            & ! Out
               )
 
       case DEFAULT
          call error_message('***ERROR: Process description for process "percolation" does not exist! mo_multi_param_reg')
       end select
 
       deallocate(ksvar_h0)
       deallocate(ksvar_v0)
       deallocate(sms_fc0)
 
    end do !! >>>>>>> LAND COVER SCENE LOOP
 
 
    ! ------------------------------------------------------------------
    ! sealed area threshold for runoff generation
    ! ------------------------------------------------------------------
    select case(processmatrix(4, 1))
    case (1)
       istart = processmatrix(4, 3) - processmatrix(4, 2) + 1
       iend = processmatrix(4, 3)
       call iper_thres_runoff(param(istart : iend), sealedthresh1)
    case DEFAULT
       call error_message('***ERROR: Process description for process "runoff_generation" does not exist! mo_multi_param_reg')
    end select
 
    ! ------------------------------------------------------------------
    ! potential evapotranspiration (PET)
    ! ------------------------------------------------------------------
    select case(processmatrix(5, 1))
    case(-1) ! LAI correction of input PET
       iend = -9999 ! dummy statement
    case(0) ! aspect correction of input PET
       istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
       iend = processmatrix(5, 3)
       call pet_correctbyasp(id0, y0, asp0, param(istart : iend), nodata_dp, fasp0)
       fasp1(:, 1, 1) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
            left_bound1, right_bound1, id0, mask0, nodata_dp, fasp0)
    case(1) ! Hargreaves-Samani method
       istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
       iend = processmatrix(5, 3)
       call pet_correctbyasp(id0, y0, asp0, param(istart : iend - 1), nodata_dp, fasp0)
       fasp1(:, 1, 1) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
            left_bound1, right_bound1, id0, mask0, nodata_dp, fasp0)
       harsamcoeff1 = param(iend)
    case(2) ! Priestley-Taylor Method
       istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
       iend = processmatrix(5, 3)
       call priestley_taylor_alpha(gridded_lai0, param(istart : iend), &
            mask0, nodata_dp, id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
            prietayalpha1(:, :, 1))
    case(3) ! Penman-Monteith method
       ! aerodynamic resistance is calculated inside LCscene loop
       istart = processmatrix(5, 3) - processmatrix(5, 2) + 1
       iend = processmatrix(5, 3)
       call bulksurface_resistance(gridded_lai0, param(iend), mask0, &
            nodata_dp, id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
            surfresist1(:, :, 1))
    case default
       call error_message('***ERROR: Process description for process "pet correction" does not exist! mo_multi_param_reg')
    end select
    ! ------------------------------------------------------------------
    ! baseflow recession parameter
    ! ------------------------------------------------------------------
    select case(processmatrix(9, 1))
      case(1)
        ! the number of process parameters, so the number in processMatrix(9,2) has
        ! to be equal to the size of geo_unit_list
        istart = processmatrix(9, 3) - processmatrix(9, 2) + 1
        iend = processmatrix(9, 3)
 
        call baseflow_param(param(istart : iend), geounit0, k2_0)
 
        ! Upscale by arithmetic mean
        allocate(k2_1(size(kbaseflow1, 1)))
        k2_1 = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
            left_bound1, right_bound1, id0, mask0, nodata_dp, k2_0)
        ! loop over all LCover scenes
        do iilc = 1, size(lcover0, 2)
          kbaseflow1(:, 1, iilc) = k2_1
        end do
        deallocate(k2_1)
 
        ! correction and unit conversion
        ! if percolation is ON: correct K2 such that it is at least k1
        ! since kSlowFlow1 is LCover dependent, kBaseFlow1 is too
        if (processmatrix(7, 1) .gt. 0) kbaseflow1 = merge(kslowflow1, kbaseflow1, kbaseflow1 .lt. kslowflow1)
 
      case DEFAULT
        call error_message('***ERROR: Process description for process "baseflow Recession" does not exist! mo_multi_param_reg')
    end select
 
    ! ------------------------------------------------------------------
    ! Neutron count related parameters
    ! >> only N0 parameter - others are defined above in soil parameters
    ! ------------------------------------------------------------------
    select case(processmatrix(10, 1))
    case(0)
       ! do nothing
    case(1)
       ! the number of process parameters, so the number in processMatrix(9,2) has
       istart = processmatrix(10, 3) - processmatrix(10, 2) + 1
       iend = processmatrix(10, 3)
       no_count1 = param(istart)  ! >> 1st parameter --> N0 parameter
    case(2)
       ! the number of process parameters, so the number in processMatrix(9,2) has
       istart = processmatrix(10, 3) - processmatrix(10, 2) + 1
       iend = processmatrix(10, 3)
       no_count1 = param(istart)  ! >> 1st parameter --> N0 parameter
    case DEFAULT
       call error_message('***ERROR: Process description for process "Neutron count" does not exist! mo_multi_param_reg')
    end select
 
 
    !-------------------------------------------------------------------
    ! call regionalization of parameters related to LAI
    ! it is now outside of mHM since LAI is now dynamic variable
    !-------------------------------------------------------------------
    call canopy_intercept_param(processmatrix, param(:), &
         gridded_lai0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, id0, mask0, &
         nodata_dp, maxinter1(:, :, 1))
 
  end subroutine mpr
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        baseflow_param
 
  !    PURPOSE
  !>       \brief baseflow recession parameter
 
  !>       \details This subroutine calculates the baseflow recession parameter
  !>       based on the geological units at the Level 0 scale. For each level 0
  !>       cell, it assigns the value specified in the parameter array param for the
  !>       geological unit in this cell.
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param(1) = GeoParam(1,:)
  !>       - param(2) = GeoParam(2,:)
  !>       - ...
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(:) :: param"       list of required parameters
  !>       \param[in] "integer(i4), dimension(:) :: geoUnit0" ids of geological units at L0
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp), dimension(:) :: k2_0" - baseflow recession parameter at Level 0
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Stephan Thober  Dec 2013 - changed intent(inout) to intent(out)
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine baseflow_param(param, geoUnit0, k2_0)
 
    use mo_mpr_global_variables, only : geounitlist
    !$ use omp_lib
 
    implicit none
 
    ! list of required parameters
    real(dp), dimension(:), intent(in) :: param
 
    ! ids of geological units at L0
    integer(i4), dimension(:), intent(in) :: geoUnit0
 
    ! - baseflow recession parameter at Level 0
    real(dp), dimension(:), intent(out) :: k2_0
 
    ! loop variable
    integer(i4) :: ii
 
    ! geo unit
    integer(i4), dimension(1) :: gg
 
 
    if (size(param) .ne. size(geounitlist)) &
          call error_message(' mo_multi_param_reg: baseflow_param: size mismatch, subroutine baseflow parameters ')
 
    k2_0 = nodata_dp
 
    !$OMP PARALLEL
    !$OMP DO PRIVATE(gg) SCHEDULE(STATIC)
    do ii = 1, size(k2_0)
       ! get parameter index in geoUnitList
       gg = minloc(abs(geounitlist - geounit0(ii)))
       k2_0(ii) = param(gg(1))
    end do
    !$OMP END DO
    !$OMP END PARALLEL
 
  end subroutine baseflow_param
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        snow_acc_melt_param
 
  !    PURPOSE
  !>       \brief Calculates the snow parameters.
 
  !>       \details This subroutine calculates the snow parameters
  !>       threshold temperature (TT), degree-day factor without precipitation (DD)
  !>       and maximum degree-day factor (DDmax) as well as increase of degree-day
  !>       factor per mm of increase in precipitation (IDDP).
 
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param(1) = snowTreshholdTemperature
  !>       - param(2) = degreeDayFactor_forest
  !>       - param(3) = degreeDayFactor_impervious
  !>       - param(4) = degreeDayFactor_pervious
  !>       - param(5) = increaseDegreeDayFactorByPrecip
  !>       - param(6) = maxDegreeDayFactor_forest
  !>       - param(7) = maxDegreeDayFactor_impervious
  !>       - param(8) = maxDegreeDayFactor_pervious
  !>       INTENT(IN)
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(8) :: param"    eight global parameters
  !>       \param[in] "real(dp), dimension(:) :: fForest1" [1] fraction of forest cover
  !>       \param[in] "real(dp), dimension(:) :: fIperm1"  [1] fraction of sealed area
  !>       \param[in] "real(dp), dimension(:) :: fPerm1"   [1] fraction of permeable area
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp), dimension(:) :: tempThresh1"  [degreeC] threshold temperature for snow rain
  !>       \param[out] "real(dp), dimension(:) :: degDayNoPre1" [mm-1 degreeC-1] Degree-day factor with
  !>       \param[out] "real(dp), dimension(:) :: degDayInc1"   [d-1 degreeC-1]  Increase of the Degree-day
  !>       \param[out] "real(dp), dimension(:) :: degDayMax1"   [mm-1 degreeC-1] Maximum Degree-day factor
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Juliane Mai    Oct 2013 - OLD parametrization --> param(1) = snowTreshholdTemperature
  !                                                --> param(2) = degreeDayFactor_forest
  !                                                --> param(3) = degreeDayFactor_impervious
  !                                                --> param(4) = degreeDayFactor_pervious
  !                                                --> param(5) = increaseDegreeDayFactorByPrecip
  !                                                --> param(6) = maxDegreeDayFactor_forest
  !                                                --> param(7) = maxDegreeDayFactor_impervious
  !                                                --> param(8) = maxDegreeDayFactor_pervious
  !                                             -------------------------------
  !                                             degreeDayFactor_impervious    = degreeDayFactor_forest + delta_1 + delta_2
  !                                             degreeDayFactor_pervious      = degreeDayFactor_forest + delta_1
  !                                             maxDegreeDayFactor_forest     = degreeDayFactor_forest + delta_3
  !                                             maxDegreeDayFactor_impervious = degreeDayFactor_impervious + delta_5
  !                                                                           = degreeDayFactor_forest + delta_1 + delta_2 + delta_5
  !                                             maxDegreeDayFactor_pervious   = degreeDayFactor_pervious + delta_4
  !                                                                           = degreeDayFactor_forest + delta_1 + delta_4
  !                                             -------------------------------
  !                                             NEW parametrization
  !                                                --> param(1) = snowTreshholdTemperature
  !                                                --> param(2) = degreeDayFactor_forest
  !                                                --> param(3) = delta_2
  !                                                --> param(4) = delta_1
  !                                                --> param(5) = increaseDegreeDayFactorByPrecip
  !                                                --> param(6) = delta_3
  !                                                --> param(7) = delta_5
  !                                                --> param(8) = delta_4
  ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine snow_acc_melt_param(param, fForest1, fIperm1, fPerm1, tempThresh1, degDayNoPre1, degDayInc1, &
       degDayMax1)
    implicit none
 
    ! eight global parameters
    real(dp), dimension(8), intent(in) :: param
 
    ! [1] fraction of forest cover
    real(dp), dimension(:), intent(in) :: fForest1
 
    ! [1] fraction of sealed area
    real(dp), dimension(:), intent(in) :: fIperm1
 
    ! [1] fraction of permeable area
    real(dp), dimension(:), intent(in) :: fPerm1
 
    ! [degreeC] threshold temperature for snow rain
    real(dp), dimension(:), intent(out) :: tempThresh1
 
    ! [mm-1 degreeC-1] Degree-day factor with
    real(dp), dimension(:), intent(out) :: degDayNoPre1
 
    ! [mm-1 degreeC-1] Maximum Degree-day factor
    real(dp), dimension(:), intent(out) :: degDayMax1
 
    ! [d-1 degreeC-1]  Increase of the Degree-day
    real(dp), dimension(:), intent(out) :: degDayInc1
 
    real(dp) :: tmp_degreeDayFactor_forest, tmp_degreeDayFactor_impervious, tmp_degreeDayFactor_pervious
 
    real(dp) :: tmp_maxDegreeDayFactor_forest, tmp_maxDegreeDayFactor_impervious, tmp_maxDegreeDayFactor_pervious
 
 
    tmp_degreedayfactor_forest = param(2)                                    ! OLD: param(2)
    tmp_degreedayfactor_impervious = param(2) + param(4) + param(3)              ! OLD: param(3)
    tmp_degreedayfactor_pervious = param(2) + param(4)                         ! OLD: param(4)
    tmp_maxdegreedayfactor_forest = param(2) + param(6)   ! OLD: param(6)
    tmp_maxdegreedayfactor_impervious = param(2) + param(4) + param(3) + param(7)   ! OLD: param(7)
    tmp_maxdegreedayfactor_pervious = param(2) + param(4) + param(8)   ! OLD: param(8)
 
    tempthresh1 = param(1)
    degdayinc1 = param(5)
 
    degdaynopre1 = (&
         tmp_degreedayfactor_forest * fforest1 + &
         tmp_degreedayfactor_impervious * fiperm1 + &
         tmp_degreedayfactor_pervious * fperm1)
    degdaymax1 = (&
         tmp_maxdegreedayfactor_forest * fforest1 + &
         tmp_maxdegreedayfactor_impervious * fiperm1 + &
         tmp_maxdegreedayfactor_pervious * fperm1)
 
  end subroutine snow_acc_melt_param
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        iper_thres_runoff
 
  !    PURPOSE
  !>       \brief sets the impervious layer threshold parameter for runoff generation
 
  !>       \details to be done by Kumar
  !>       ....
 
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param(1) = imperviousStorageCapacity
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(1) :: param" - given threshold parameter
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp), dimension(:, :, :) :: sealedThresh1"
 
  !    HISTORY
  !>       \authors Stephan Thober, Rohini Kumar
 
  !>       \date Dec 2012
 
  ! Modifications:
  ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine iper_thres_runoff(param, sealedThresh1)
    implicit none
 
    ! - given threshold parameter
    real(dp), dimension(1), intent(in) :: param
 
    real(dp), dimension(:, :, :), intent(out) :: sealedThresh1
 
 
    sealedthresh1 = param(1)
 
  end subroutine iper_thres_runoff
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        karstic_layer
 
  !    PURPOSE
  !>       \brief calculates the Karstic percolation loss
 
  !>       \details This subroutine calls first the karstic_fraction upscaling
  !>       routine for determine the karstic fraction area for every Level 1
  !>       cell. Then, the karstic percolation loss is estimated given two
  !>       shape parameters by
  !>       \f[ karstLoss1 = 1 + ( fKarArea * param(1)) *( (-1)**INT(param(2),i4) ) \f]
  !>       where \f$ karstLoss1 \f$ is the karstic percolation loss and \f$ fKarArea \f$
  !>       is the fraction of karstic area at level 1
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param(1) = rechargeCoefficient
  !>       - param(2) = rechargeFactor_karstic
  !>       - param(3) = gain_loss_GWreservoir_karstic
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(3) :: param"           parameters
  !>       \param[in] "integer(i4), dimension(:) :: geoUnit0"     id of the Karstic formation
  !>       \param[in] "logical, dimension(:, :) :: mask0"         mask at level 0
  !>       \param[in] "real(dp), dimension(:) :: SMs_FC0"         [-] soil mositure deficit from field
  !>       capacity w.r.t to saturation
  !>       \param[in] "real(dp), dimension(:) :: KsVar_V0"        [-] relative variability of saturated
  !>       \param[in] "integer(i4), dimension(:) :: Id0"          Cell ids of hi res field
  !>       \param[in] "integer(i4), dimension(:) :: n_subcells1"  number of l0 cells within a l1 cell
  !>       \param[in] "integer(i4), dimension(:) :: upper_bound1" upper row of a l1 cell in l0 grid
  !>       \param[in] "integer(i4), dimension(:) :: lower_bound1" lower row of a l1 cell in l0 grid
  !>       \param[in] "integer(i4), dimension(:) :: left_bound1"  left col of a l1 cell in l0 grid
  !>       \param[in] "integer(i4), dimension(:) :: right_bound1" right col of a l1 cell in l0 grid
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp), dimension(:) :: karstLoss1" [-]    Karstic percolation loss
  !>       \param[out] "real(dp), dimension(:) :: L1_Kp"      [d-1] percolation coefficient
 
  !    HISTORY
  !>       \authors Rohini Kumar, Stephan Thober
 
  !>       \date Feb 2013
 
  ! Modifications:
  ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
  ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine karstic_layer(param, geoUnit0, mask0, SMs_FC0, KsVar_V0, Id0, n_subcells1, upper_bound1, lower_bound1, &
       left_bound1, right_bound1, karstLoss1, L1_Kp)
 
    use mo_mpr_global_variables, only : geounitlist, geounitkar
    use mo_upscaling_operators, only : l0_fractionalcover_in_lx, upscale_arithmetic_mean
    !$ use omp_lib
 
    implicit none
 
    ! parameters
    real(dp), dimension(3), intent(in) :: param
 
    ! id of the Karstic formation
    integer(i4), dimension(:), intent(in) :: geoUnit0
 
    ! mask at level 0
    logical, dimension(:, :), intent(in) :: mask0
 
    ! [-] soil mositure deficit from field
    ! capacity w.r.t to saturation
    real(dp), dimension(:), intent(in) :: SMs_FC0
 
    ! [-] relative variability of saturated
    real(dp), dimension(:), intent(in) :: KsVar_V0
 
    ! Cell ids of hi res field
    integer(i4), dimension(:), intent(in) :: Id0
 
    ! number of l0 cells within a l1 cell
    integer(i4), dimension(:), intent(in) :: n_subcells1
 
    ! upper row of a l1 cell in l0 grid
    integer(i4), dimension(:), intent(in) :: upper_bound1
 
    ! lower row of a l1 cell in l0 grid
    integer(i4), dimension(:), intent(in) :: lower_bound1
 
    ! left col of a l1 cell in l0 grid
    integer(i4), dimension(:), intent(in) :: left_bound1
 
    ! right col of a l1 cell in l0 grid
    integer(i4), dimension(:), intent(in) :: right_bound1
 
    ! [-]    Karstic percolation loss
    real(dp), dimension(:), intent(out) :: karstLoss1
 
    ! [d-1] percolation coefficient
    real(dp), dimension(:), intent(out) :: L1_Kp
 
    ! fraction of karstic area
    real(dp), dimension(:), allocatable :: fKarArea
 
    ! temporal variable
    real(dp), dimension(size(SMs_FC0, 1)) :: tmp
 
    integer(i4) :: nGeoUnits
 
    integer(i4) :: i
 
 
    ! ------------------------------------------------------------------
    ! PERCOLATION; 1/Kp = f(Ks)
    ! ------------------------------------------------------------------
    ! Regionalise Kp with variability of last soil layer property
    !$OMP PARALLEL
    tmp = merge(param(1) * (1.0_dp + sms_fc0) / (1.0_dp + ksvar_v0), &
         nodata_dp, id0 .ne. nodata_i4)
    !$OMP END PARALLEL
 
    l1_kp = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
         left_bound1, right_bound1, id0, mask0, nodata_dp, tmp)
 
    ! minimum constrains
    l1_kp = merge(2.0_dp, l1_kp, l1_kp .lt. 2.0_dp)
 
    ngeounits = size(geounitlist, 1)
 
    ! 1st calculate fraction of Karstic area
    allocate(fkararea(size(karstloss1, 1)))
    fkararea = 0.0_dp
 
    do i = 1, ngeounits
       if(geounitkar(i) .eq. 0) cycle
       fkararea(:) = l0_fractionalcover_in_lx(geounit0, geounitlist(i), mask0, &
            upper_bound1, lower_bound1, left_bound1, right_bound1, n_subcells1)
    end do
 
    ! 2nd calculate karstLoss1
    karstloss1 = 1.0_dp - (fkararea * param(2))
 
    deallocate(fkararea)
 
  end subroutine karstic_layer
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        canopy_intercept_param
 
  !    PURPOSE
  !>       \brief estimate effective maximum interception capacity at L1
 
  !>       \details estimate effective maximum interception capacity at L1 for a given
  !>       Leaf Area Index field.
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       Process Case 1:
  !>       - param(1) = canopyInterceptionFactor
 
  !    INTENT(IN)
  !>       \param[in] "integer(i4), dimension(:, :) :: processMatrix" indicate processes
  !>       \param[in] "real(dp), dimension(:) :: param"               array of global parameters
  !>       \param[in] "real(dp), dimension(:, :) :: LAI0"             LAI at level-0(nCells0, time)
  !>       \param[in] "integer(i4), dimension(:) :: n_subcells1"      Number of L0 cells within a L1 cell
  !>       \param[in] "integer(i4), dimension(:) :: upper_bound1"     Upper row of high resolution block
  !>       \param[in] "integer(i4), dimension(:) :: lower_bound1"     Lower row of high resolution block
  !>       \param[in] "integer(i4), dimension(:) :: left_bound1"      Left column of high resolution block
  !>       \param[in] "integer(i4), dimension(:) :: right_bound1"     Right column of high resolution block
  !>       \param[in] "integer(i4), dimension(:) :: Id0"              Cell ids at level 0
  !>       \param[in] "logical, dimension(:, :) :: mask0"             mask at level 0 field
  !>       \param[in] "real(dp) :: nodata"                            - nodata value
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp), dimension(:, :) :: max_intercept1" max interception at level-1(nCells1, time)
 
  !    HISTORY
  !>       \authors Rohini Kumar
 
  !>       \date Aug. 2013
 
  ! Modifications:
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine canopy_intercept_param(processMatrix, param, LAI0, n_subcells1, upper_bound1, lower_bound1, left_bound1, &
       right_bound1, Id0, mask0, nodata, max_intercept1)
 
    use mo_string_utils, only : num2str
    use mo_upscaling_operators, only : upscale_arithmetic_mean
 
    implicit none
 
    ! indicate processes
    integer(i4), dimension(:, :), intent(in) :: processmatrix
 
    ! array of global parameters
    real(dp), dimension(:), intent(in) :: param
 
    ! LAI at level-0(nCells0, time)
    real(dp), dimension(:, :), intent(in) :: lai0
 
    ! Number of L0 cells within a L1 cell
    integer(i4), dimension(:), intent(in) :: n_subcells1
 
    ! Upper row of high resolution block
    integer(i4), dimension(:), intent(in) :: upper_bound1
 
    ! Lower row of high resolution block
    integer(i4), dimension(:), intent(in) :: lower_bound1
 
    ! Left column of high resolution block
    integer(i4), dimension(:), intent(in) :: left_bound1
 
    ! Right column of high resolution block
    integer(i4), dimension(:), intent(in) :: right_bound1
 
    ! Cell ids at level 0
    integer(i4), dimension(:), intent(in) :: id0
 
    ! mask at level 0 field
    logical, dimension(:, :), intent(in) :: mask0
 
    ! - nodata value
    real(dp), intent(in) :: nodata
 
    ! max interception at level-1(nCells1, time)
    real(dp), dimension(:, :), intent(out) :: max_intercept1
 
    integer(i4) :: istart, iend, it
 
    real(dp), dimension(:), allocatable :: max_intercept0
 
    real(dp), dimension(:), allocatable :: gamma_intercept
 
 
    ! ------------------------------------------------------------------
    ! Maximum interception parameter
    ! ------------------------------------------------------------------
    select case(processmatrix(1, 1))
    case(1)
       istart = processmatrix(1, 3) - processmatrix(1, 2) + 1
       iend = processmatrix(1, 3)
 
       ! allocate space
       allocate(gamma_intercept(iend - istart + 1))
       allocate(max_intercept0(size(id0, 1)))
 
       ! estimate max. intercept at Level-0
       gamma_intercept(:) = param(istart : iend)
 
       do it = 1, size(lai0, 2)
          !$OMP PARALLEL
          max_intercept0(:) = lai0(:, it) * gamma_intercept(1)
          !$OMP END PARALLEL
 
          ! Upscale by arithmetic mean
          max_intercept1(:, it) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, left_bound1, &
               right_bound1, id0, mask0, nodata, max_intercept0(:))
 
       end do
 
       deallocate(gamma_intercept)
       deallocate(max_intercept0)
    CASE DEFAULT
       call error_message('mo_multi_param_reg: This processMatrix=', num2str(processmatrix(1, 1)), ' is not implemented!')
    end select
 
  end subroutine canopy_intercept_param
 
 
  ! ----------------------------------------------------------------------------
 
  !    NAME
  !        aerodynamical_resistance
 
  !    PURPOSE
  !>       \brief Regionalization of aerodynamic resistance
 
  !>       \details estimation of aerodynamical resistance
  !>       Global parameters needed (see mhm_parameter.nml):
  !>       - param(1) = canopyheigth_forest
  !>       - param(2) = canopyheigth_impervious
  !>       - param(3) = canopyheigth_pervious
  !>       - param(4) = displacementheight_coeff
  !>       - param(5) = roughnesslength_momentum_coeff
  !>       - param(6) = roughnesslength_heat_coeff
 
  !    INTENT(IN)
  !>       \param[in] "real(dp), dimension(:, :) :: LAI0"         LAI at level-0
  !>       \param[in] "integer(i4), dimension(:) :: LCover0"      land cover field
  !>       \param[in] "real(dp), dimension(6) :: param"           input parameter
  !>       \param[in] "logical, dimension(:, :) :: mask0"         mask at level 0
  !>       \param[in] "integer(i4), dimension(:) :: Id0"          Cell ids of hi res field
  !>       \param[in] "integer(i4), dimension(:) :: n_subcells1"  number of l0 cells within a l1 cell
  !>       \param[in] "integer(i4), dimension(:) :: upper_bound1" upper row of a l1 cell in l0 grid
  !>       \param[in] "integer(i4), dimension(:) :: lower_bound1" lower row of a l1 cell in l0 grid
  !>       \param[in] "integer(i4), dimension(:) :: left_bound1"  left col of a l1 cell in l0 grid
  !>       \param[in] "integer(i4), dimension(:) :: right_bound1" right col of a l1 cell in l0 grid
 
  !    INTENT(OUT)
  !>       \param[out] "real(dp), dimension(:, :) :: aerodyn_resistance1" aerodynmaical resistance
 
  !    HISTORY
  !>       \authors Matthias Zink
 
  !>       \date Apr 2013
 
  ! Modifications:
  ! Matthias Zink Jun 2017 - moved from mo_multi_scale_param_reg.f90 to mo_mpr_pet.f90
  ! Robert Schweppe Jun 2018 - refactoring and reformatting
 
  subroutine aerodynamical_resistance(LAI0, LCover0, param, mask0, Id0, n_subcells1, upper_bound1, lower_bound1, &
       left_bound1, right_bound1, aerodyn_resistance1)
 
    use mo_common_constants, only : eps_dp
    use mo_mpr_constants, only : windmeasheight, karman
    use mo_upscaling_operators, only : upscale_arithmetic_mean
 
    implicit none
 
    ! LAI at level-0
    real(dp), dimension(:, :), intent(in) :: LAI0
 
    ! land cover field
    integer(i4), dimension(:), intent(in) :: LCover0
 
    ! input parameter
    real(dp), dimension(6), intent(in) :: param
 
    ! mask at level 0
    logical, dimension(:, :), intent(in) :: mask0
 
    ! Cell ids of hi res field
    integer(i4), dimension(:), intent(in) :: Id0
 
    ! number of l0 cells within a l1 cell
    integer(i4), dimension(:), intent(in) :: n_subcells1
 
    ! upper row of a l1 cell in l0 grid
    integer(i4), dimension(:), intent(in) :: upper_bound1
 
    ! lower row of a l1 cell in l0 grid
    integer(i4), dimension(:), intent(in) :: lower_bound1
 
    ! left col of a l1 cell in l0 grid
    integer(i4), dimension(:), intent(in) :: left_bound1
 
    ! right col of a l1 cell in l0 grid
    integer(i4), dimension(:), intent(in) :: right_bound1
 
    ! aerodynmaical resistance
    real(dp), dimension(:, :), intent(out) :: aerodyn_resistance1
 
    integer(i4) :: tt
 
    real(dp), dimension(:), allocatable :: maxLAI
 
    real(dp), dimension(:), allocatable :: zm
 
    real(dp), dimension(:), allocatable :: canopy_height0
 
    real(dp), dimension(:), allocatable :: zm_zero, zh_zero, displace
 
    ! dim 1 = number of cells on level 0,
    ! dim2=month of year
    real(dp), dimension(:, :), allocatable :: aerodyn_resistance0
 
 
    ! initialize some things
    allocate(zm(size(lcover0, dim = 1))) ; zm = nodata_dp
    allocate(zm_zero(size(lcover0, dim = 1))) ; zm_zero = nodata_dp
    allocate(zh_zero(size(lcover0, dim = 1))) ; zh_zero = nodata_dp
    allocate(displace(size(lcover0, dim = 1))) ; displace = nodata_dp
    allocate(canopy_height0(size(lcover0, dim = 1))) ; canopy_height0 = nodata_dp
    allocate(aerodyn_resistance0(size(lcover0, dim = 1), size(lai0, 2))) ; aerodyn_resistance0 = nodata_dp
    allocate(maxlai(size(lcover0, dim = 1))) ; maxlai = nodata_dp
 
    ! regionalization of canopy height
    ! substitute with canopy height
    canopy_height0 = merge(param(1), canopy_height0, lcover0 == 1)  ! forest
    canopy_height0 = merge(param(2), canopy_height0, lcover0 == 2)  ! impervious
 
    ! old implementation used values from LUT statically for all cells (Jan-Dec, 12 values):
    ! 7 Intensive-orchards: 2.0, 2.0,  2.0,  2.0,  3.0, 3.5,  4.0,  4.0,  4.0,  2.5,  2.0,  2.0
    maxlai = maxval(lai0, dim=2)
 
    do tt = 1, size(lai0, 2)
 
       ! pervious canopy height is scaled with LAI
       canopy_height0 = merge((param(3) * lai0(:, tt) / maxlai), canopy_height0, lcover0 == 3)  ! pervious
 
       ! estimation of the aerodynamic resistance on the lower level
       ! see FAO Irrigation and Drainage Paper No. 56 (p. 19 ff) for more information
       zm = windmeasheight
       ! correction: if wind measurement height is below canopy height loagarithm becomes negative
       zm = merge(canopy_height0 + zm, zm, ((abs(zm - nodata_dp) .GT. eps_dp) .AND. (zm .LT. canopy_height0)))
 
       ! zh       = zm
       displace = param(4) * canopy_height0
       zm_zero = param(5) * canopy_height0
       zh_zero = param(6) * zm_zero
       !
       ! calculate aerodynamic resistance (changes monthly)
       aerodyn_resistance0(:, tt) = log((zm - displace) / zm_zero) * log((zm - displace) / zh_zero) / (karman**2.0_dp)
       aerodyn_resistance1(:, tt) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
            left_bound1, right_bound1, id0, mask0, nodata_dp, aerodyn_resistance0(:, tt))
 
    end do
 
  end subroutine aerodynamical_resistance
 
END MODULE mo_multi_param_reg