LCOV - code coverage report
Current view: top level - MPR - mo_multi_param_reg.f90 (source / functions) Hit Total Coverage
Test: mHM coverage Lines: 256 278 92.1 %
Date: 2024-04-30 08:53:32 Functions: 7 7 100.0 %

          Line data    Source code
       1             : !> \file mo_multi_param_reg.f90
       2             : !> \brief \copybrief mo_multi_param_reg
       3             : !> \details \copydetails mo_multi_param_reg
       4             : 
       5             : !> \brief Multiscale parameter regionalization (MPR).
       6             : !> \details This module provides the routines for multiscale parameter regionalization (MPR).
       7             : !> \authors Stephan Thober, Rohini Kumar
       8             : !> \date Dec 2012
       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_mpr
      12             : MODULE mo_multi_param_reg
      13             : 
      14             :   use mo_kind, only : i4, dp
      15             :   use mo_common_constants, only : nodata_dp, nodata_i4
      16             :   use mo_message, only : message, error_message
      17             : 
      18             :   implicit none
      19             : 
      20             :   private
      21             : 
      22             :   PUBLIC :: mpr                     ! calculates effective regionalised parameters
      23             :   PUBLIC :: canopy_intercept_param  ! estimate effective max. canopy interception
      24             : 
      25             : contains
      26             : 
      27             :   !> \brief Regionalizing and Upscaling process parameters
      28             :   !> \details calculating process parameters at L0 scale (Regionalization), like:
      29             :   !! - Baseflow recession parameter
      30             :   !! - Soil moisture parameters
      31             :   !! - PET correction for aspect
      32             :   !!
      33             :   !! and upscale these parameters to retrieve effective parameters at scale L1.
      34             :   !! Further parameter regionalizations are done for:
      35             :   !! - snow accumulation and melting parameters
      36             :   !! - threshold parameter for runoff generation on impervious layer
      37             :   !! - karstic percolation loss
      38             :   !! - setting up the Regionalized Routing Parameters
      39             :   !!
      40             :   !> \changelog
      41             :   !! - Stephan Thober           Jan 2013
      42             :   !!   - updated calling sequence for upscaling operators
      43             :   !! - Luis Samaniego           Feb 2013
      44             :   !!   - calling sequence, initial CHECK, call mpr_runoff
      45             :   !! - Stephan Thober           Feb 2013
      46             :   !!   - added subroutine for karstic percolation loss removed L1_, L0_ in variable names
      47             :   !! - Stephan Thober           Aug 2015
      48             :   !!   - moved regionalization of routing to mRM
      49             :   !! - Rohini Kumar             Mar 2016
      50             :   !!   - changes for handling multiple soil database options
      51             :   !! - Zink M. & Demirel M.C.   Mar 2017
      52             :   !!   - Added Jarvis soil water stress function at SM process(3)
      53             :   !! - Demirel M.C. & S. Stisen Apr 2017
      54             :   !!   - Added FC dependency on root fraction coefficient at SM process(3)
      55             :   !! - Robert Schweppe          Dec 2017
      56             :   !!   - added loop over LCscenes inside MPR, renamed variables rewrite
      57             :   !! - Robert Schweppe          Jun 2018
      58             :   !!   - refactoring and reformatting
      59             :   !! - Demirel M.C. & S. Stisen Jun 2020
      60             :   !!   - Added Feddes and global FC dependency on root fraction coefficient at SM process(3)=4
      61             :   !! - Rohini Kumar             Oct 2021
      62             :   !!   - Added Neutron count module to mHM integrate into develop branch (5.11.2)
      63             :   !! - Sebastian Müller         Mar 2023
      64             :   !!   - made L1_alpha, L1_kSlowFlow, L1_kBaseFlow and L1_kPerco land cover dependent
      65             :   !> \authors Stephan Thober, Rohini Kumar
      66             :   !> \date Dec 2012
      67         621 : subroutine mpr(mask0, geoUnit0, soilId0, Asp0, gridded_LAI0, LCover0, slope_emp0, y0, Id0, upper_bound1, lower_bound1, &
      68        1035 :        left_bound1, right_bound1, n_subcells1, fSealed1, alpha1, degDayInc1, degDayMax1, degDayNoPre1, fAsp1, &
      69         621 :        HarSamCoeff1, PrieTayAlpha1, aeroResist1, surfResist1, fRoots1, kFastFlow1, kSlowFlow1, kBaseFlow1, &
      70         621 :        kPerco1, karstLoss1, soilMoistFC1, soilMoistSat1, soilMoistExp1, jarvis_thresh_c1, tempThresh1, &
      71         621 :        unsatThresh1, sealedThresh1, wiltingPoint1, maxInter1, petLAIcorFactor, &
      72         414 :        No_Count1, bulkDens1, latticeWater1, COSMICL31, &
      73         207 :        parameterset )
      74             : 
      75             :     use mo_common_variables, only : global_parameters, processMatrix
      76             :     use mo_mpr_SMhorizons, only : mpr_SMhorizons
      77             :     use mo_mpr_global_variables, only : HorizonDepth_mHM, fracSealed_CityArea, iFlag_soilDB, nSoilHorizons_mHM, &
      78             :          soilDB
      79             :     use mo_mpr_pet, only : bulksurface_resistance, pet_correctbyASP, pet_correctbyLAI, priestley_taylor_alpha
      80             :     use mo_mpr_runoff, only : mpr_runoff
      81             :     use mo_mpr_soilmoist, only : mpr_sm
      82             :     use mo_upscaling_operators, only : L0_fractionalCover_in_Lx, &
      83             :          upscale_arithmetic_mean
      84             :     use mo_mpr_neutrons,        only: mpr_neutrons
      85             :     implicit none
      86             : 
      87             :     !> mask at level 0 field
      88             :     logical, dimension(:, :), intent(in) :: mask0
      89             :     !> L0 geological units
      90             :     integer(i4), dimension(:), intent(in) :: geoUnit0
      91             :     !> soil Ids at level 0
      92             :     integer(i4), dimension(:, :), intent(in) :: soilId0
      93             :     !> [degree] Aspect at Level 0
      94             :     real(dp), dimension(:), intent(in) :: Asp0
      95             :     !> LAI grid at level 0, with dim2 = time
      96             :     real(dp), dimension(:, :), intent(in) :: gridded_LAI0
      97             :     !> land cover at level 0
      98             :     integer(i4), dimension(:, :), intent(in) :: LCOVER0
      99             :     !> Empirical quantiles of slope
     100             :     real(dp), dimension(:), intent(in) :: slope_emp0
     101             :     !> Cell ids at level 0
     102             :     integer(i4), dimension(:), intent(in) :: Id0
     103             :     !> Upper row of hi res block
     104             :     integer(i4), dimension(:), intent(in) :: upper_bound1
     105             :     !> Lower row of hi res block
     106             :     integer(i4), dimension(:), intent(in) :: lower_bound1
     107             :     !> Left column of hi res block
     108             :     integer(i4), dimension(:), intent(in) :: left_bound1
     109             :     !> Right column of hi res block
     110             :     integer(i4), dimension(:), intent(in) :: right_bound1
     111             :     !> Number of L0 cells within a L1 cell
     112             :     integer(i4), dimension(:), intent(in) :: n_subcells1
     113             :     !> y0 at level 0
     114             :     real(dp), dimension(:), intent(in) :: y0
     115             :     !> [1] fraction of sealed area
     116             :     real(dp), dimension(:, :, :), intent(inout) :: fSealed1
     117             :     !> Parameter that determines the rel. contribution to SM, upscal. Bulk den.
     118             :     real(dp), dimension(:, :, :), intent(inout) :: soilMoistExp1
     119             :     !> [1] jarvis critical value for norm SWC
     120             :     real(dp), dimension(:, :, :), intent(inout) :: jarvis_thresh_c1
     121             :     !> [10^-3 m] depth of saturated SM
     122             :     real(dp), dimension(:, :, :), intent(inout) :: soilMoistSat1
     123             :     !> [10^-3 m] field capacity
     124             :     real(dp), dimension(:, :, :), intent(inout) :: soilMoistFC1
     125             :     !> [10^-3 m] permanent wilting point
     126             :     real(dp), dimension(:, :, :), intent(inout) :: wiltingPoint1
     127             :     !> fraction of roots in soil horizon
     128             :     real(dp), dimension(:, :, :), intent(inout) :: fRoots1
     129             :     !> [degreeC] threshold temperature for snow rain
     130             :     real(dp), dimension(:, :, :), intent(inout) :: tempThresh1
     131             :     !> [mm-1 degreeC-1] Degree-day factor with no precipitation
     132             :     real(dp), dimension(:, :, :), intent(inout) :: degDayNoPre1
     133             :     !> [mm-1 degreeC-1] Maximum Degree-day factor
     134             :     real(dp), dimension(:, :, :), intent(inout) :: degDayMax1
     135             :     !> [d-1 degreeC-1]  Increase of the Degree-day factor per mm of increase in precipitation
     136             :     real(dp), dimension(:, :, :), intent(inout) :: degDayInc1
     137             :     !> [1]     PET correction for Aspect at level 1
     138             :     real(dp), dimension(:, :, :), intent(inout) :: fAsp1
     139             :     !> [1]     PET Hargreaves Samani coeff. at level 1
     140             :     real(dp), dimension(:, :, :), intent(inout) :: HarSamCoeff1
     141             :     !> [1]     PET Priestley Taylor coeff. at level 1
     142             :     real(dp), dimension(:, :, :), intent(inout) :: PrieTayAlpha1
     143             :     !> [s m-1] PET aerodynamical resitance at level 1
     144             :     real(dp), dimension(:, :, :), intent(inout) :: aeroResist1
     145             :     !> [s m-1] PET bulk surface resitance at level 1
     146             :     real(dp), dimension(:, :, :), intent(inout) :: surfResist1
     147             :     !> threshold parameter
     148             :     real(dp), dimension(:, :, :), intent(inout) :: sealedThresh1
     149             :     !> [10^-3 m] Threshhold water depth in upper reservoir (for Runoff contribution)
     150             :     real(dp), dimension(:, :, :), intent(inout) :: unsatThresh1
     151             :     !> [10^-3 m] Recession coefficient of the upper reservoir, upper outlet
     152             :     real(dp), dimension(:, :, :), intent(inout) :: kFastFlow1
     153             :     !> [10^-3 m] Recession coefficient of the upper reservoir, lower outlet
     154             :     real(dp), dimension(:, :, :), intent(inout) :: kSlowFlow1
     155             :     !> Level 1 baseflow recession
     156             :     real(dp), dimension(:, :, :), intent(inout) :: kBaseFlow1
     157             :     !> [1] Exponent for the upper reservoir
     158             :     real(dp), dimension(:, :, :), intent(inout) :: alpha1
     159             :     !> [d-1] percolation coefficient
     160             :     real(dp), dimension(:, :, :), intent(inout) :: kPerco1
     161             :     !> karstic percolation loss
     162             :     real(dp), dimension(:, :, :), intent(inout) :: karstLoss1
     163             :     !> max interception
     164             :     real(dp), dimension(:, :, :), intent(inout) :: maxInter1
     165             :     !> pet cor factor at level-1
     166             :     real(dp), dimension(:, :, :), intent(inout) :: petLAIcorFactor
     167             :     !> [-] inital neutron count
     168             :     real(dp), dimension(:, :, :), intent(inout) :: No_Count1
     169             :     !> [gcm-3] bulk density
     170             :     real(dp), dimension(:, :, :), intent(inout) :: bulkDens1
     171             :     !> [mm/mm] lattice water content
     172             :     real(dp), dimension(:, :, :), intent(inout) :: latticeWater1
     173             :     !> [-] cosmic L3 parameter
     174             :     real(dp), dimension(:, :, :), intent(inout) :: COSMICL31
     175             : 
     176             :     !> array of global parameters
     177             :     real(dp), dimension(:), intent(in), optional, target :: parameterset
     178             : 
     179             :     ! array of global parameters
     180         207 :     real(dp), dimension(:), pointer :: param
     181             : 
     182         207 :     real(dp), dimension(:, :, :), allocatable :: thetaS_till
     183             : 
     184         207 :     real(dp), dimension(:, :, :), allocatable :: thetaFC_till
     185             : 
     186         207 :     real(dp), dimension(:, :, :), allocatable :: thetaPW_till
     187             : 
     188             :     ! saturated hydraulic conductivity
     189         207 :     real(dp), dimension(:, :, :), allocatable :: Ks
     190             : 
     191             :     ! Bulk density
     192         207 :     real(dp), dimension(:, :, :), allocatable :: Db
     193         207 :     real(dp), dimension(:, :), allocatable :: thetaS
     194         207 :     real(dp), dimension(:, :), allocatable :: thetaFC
     195         207 :     real(dp), dimension(:, :), allocatable :: thetaPW
     196             : 
     197             :     ! neutron count
     198         207 :     real(dp), dimension(:,:,:), allocatable :: latWat_till
     199         207 :     real(dp), dimension(:,:,:), allocatable :: COSMIC_L3_till
     200         207 :     real(dp), dimension(:,:), allocatable   :: latWat         ! lattice water
     201         207 :     real(dp), dimension(:,:), allocatable   :: COSMIC_L3      ! COSMIC parameter L3
     202             : 
     203             : 
     204             :     ! relative variability of saturated
     205             :     ! hydraulic cound. for Horizantal flow
     206         207 :     real(dp), dimension(:), allocatable :: KsVar_H0
     207             : 
     208             :     ! relative variability of saturated
     209             :     ! hydraulic cound. for vertical flow
     210         207 :     real(dp), dimension(:), allocatable :: KsVar_V0
     211             : 
     212             :     ! soil mositure deficit from
     213             :     ! field cap. w.r.t to saturation
     214         207 :     real(dp), dimension(:), allocatable :: SMs_FC0
     215             : 
     216             :     ! L0 baseflow parameter
     217     9583982 :     real(dp), dimension(size(Id0, 1)) :: k2_0
     218             : 
     219             :     ! L1 baseflow parameter
     220         207 :     real(dp), dimension(:), allocatable :: k2_1
     221             : 
     222             :     ! L0 Aspect
     223     9584189 :     real(dp), dimension(size(Id0, 1)) :: fAsp0
     224             : 
     225             :     ! number of soil classes
     226             :     integer(i4) :: mSoil
     227             : 
     228             :     ! maximum of number of Tillage horizons
     229             :     integer(i4) :: mTill
     230             : 
     231             :     ! maximum number of horizons
     232             :     integer(i4) :: mHor
     233             : 
     234             :     ! number of Landcover classes
     235             :     integer(i4) :: mLC
     236             : 
     237             :     ! indexing of parameter vector - start
     238             :     integer(i4) :: iStart
     239             : 
     240             :     ! indexing of parameter vector - end
     241             :     integer(i4) :: iEnd
     242             : 
     243             :     ! 2nd indexing of parameter vector - start
     244             :     integer(i4) :: iStart2
     245             : 
     246             :     ! 2nd indexing of parameter vector - end
     247             :     integer(i4) :: iEnd2
     248             : 
     249             :     ! counter for looping over LCscenes
     250             :     integer(i4) :: iiLC
     251             : 
     252             :     ! [1]  Fraction of forest cover
     253        8022 :     real(dp), dimension(size(fSealed1, dim = 1)) :: fForest1
     254             : 
     255             :     ! [1]  Fraction of permeable cover
     256        8022 :     real(dp), dimension(size(fSealed1, dim = 1)) :: fPerm1
     257             : 
     258             : 
     259         207 :     if (present(parameterset)) then
     260         207 :        param => parameterset
     261             :     else
     262           0 :        param => global_parameters(:, 3)
     263             :     end if
     264             : 
     265             : 
     266             :     ! loop over all LCover scenes
     267         621 :     do iiLC = 1, size(LCover0, 2)
     268             : 
     269             :        ! estimate land cover fractions for dominant landcover class
     270             :        ! fSealed is intent inout, the rest only intent in
     271             :        fForest1(:) = L0_fractionalCover_in_Lx(LCover0(:, iiLC), 1, mask0, &
     272             :             upper_bound1, &
     273             :             lower_bound1, &
     274             :             left_bound1, &
     275             :             right_bound1, &
     276         414 :             n_subcells1)
     277         828 :        fSealed1(:, 1, iiLC) = L0_fractionalCover_in_Lx(LCover0(:, iiLC), 2, mask0, &
     278             :             upper_bound1, &
     279             :             lower_bound1, &
     280             :             left_bound1, &
     281             :             right_bound1, &
     282       16458 :             n_subcells1)
     283             :        fPerm1(:) = L0_fractionalCover_in_Lx(LCover0(:, iiLC), 3, mask0, &
     284             :             upper_bound1, &
     285             :             lower_bound1, &
     286             :             left_bound1, &
     287             :             right_bound1, &
     288         414 :             n_subcells1)
     289             :        !---------------------------------------------------------
     290             :        ! Update fractions of sealed area fractions
     291             :        ! based on the sealing fraction[0-1] in cities
     292             :        !---------------------------------------------------------
     293             :        ! a factor is applied to the sealed area, effectively reducing it
     294       15630 :        fSealed1(:, 1, iiLC) = fracSealed_CityArea * fSealed1(:, 1, iiLC)
     295             :        ! the forest area is kept constant, but the permeable area is increased so that the
     296             :        ! sum off all fractions equals 1 again
     297       16044 :        fPerm1(:) = 1.0_dp - fSealed1(:, 1, iiLC) - fForest1(:)
     298             : 
     299             :        ! ------------------------------------------------------------------
     300             :        ! snow parameters
     301             :        ! ------------------------------------------------------------------
     302         828 :        select case(processMatrix(2,1))
     303             :        case(1)
     304             : 
     305         414 :           iStart = processMatrix(2, 3) - processMatrix(2, 2) + 1
     306         414 :           iEnd = processMatrix(2, 3)
     307             : 
     308           0 :           call snow_acc_melt_param(param(iStart : iEnd), & ! intent(in)
     309             :                fForest1, fSealed1(:, 1, iiLC), fPerm1, & ! intent(in)
     310             :                tempThresh1(:, 1, iiLC), degDayNoPre1(:, 1, iiLC), & ! intent(out)
     311             :                degDayInc1(:, 1, iiLC), degDayMax1(:, 1, iiLC) & ! intent(out)
     312         414 :                )
     313             :        case DEFAULT
     314         414 :           call error_message('***ERROR: Process description for process "snow pack" does not exist! mo_multi_param_reg')
     315             :        end select
     316             : 
     317             :        ! ------------------------------------------------------------------
     318             :        ! Soil moisture parametrization
     319             :        ! ------------------------------------------------------------------
     320         414 :        msoil = size(soilDB%is_present, 1)
     321    19167964 :        mLC = maxval(LCover0(:, iiLC), (LCover0(:, iiLC) .ne. nodata_i4))
     322             : 
     323             :        ! depending on which kind of soil database processing is to be performed
     324         414 :        if(iFlag_soilDB .eq. 0)then
     325      611064 :           mtill = maxval(soilDB%nTillHorizons, (soilDB%nTillHorizons .ne. nodata_i4))
     326      611064 :           mHor  = maxval(soilDB%nHorizons,     (soilDB%nHorizons     .ne. nodata_i4))
     327           0 :        else if(iFlag_soilDB .eq. 1) then
     328             :           ! here for each soil type both till and non-till soil hydraulic properties are to be estimated
     329             :           ! since a given soil type can lie in any horizon (till or non-till ones)
     330             :           ! adopt it in a way that it do not break the consistency of iFlag_soilDB = 0
     331             :           ! ** NOTE: SDB_nTillHorizons and SDB_nHorizons are also assigned in
     332             :           !          this flag option (see mo_soildatabase.f90 file - read_soil_LUT).
     333             :           !          But we are not using those variables here since in this case we have not
     334             :           !          varying number of soil horizons or either tillage horizons.
     335             :           !          So assigning them with a value = 1 is more than enough.
     336           0 :           mtill = 1
     337           0 :           mHor  = 1
     338             :        end if
     339             : 
     340        2070 :        allocate(thetaS_till(msoil, mtill,  mLC))
     341        1656 :        allocate(thetaFC_till(msoil, mtill, mLC))
     342        1656 :        allocate(thetaPW_till(msoil, mtill, mLC))
     343        1656 :        allocate(thetaS(msoil,  mHor))
     344        1242 :        allocate(thetaFC(msoil, mHor))
     345        1242 :        allocate(thetaPW(msoil, mHor))
     346        2070 :        allocate(Ks(msoil, mHor, mLC))
     347        1656 :        allocate(Db(msoil, mHor, mLC))
     348             : 
     349             :        ! neutron count related ones
     350             :        ! allocate and initalize here
     351        1656 :        allocate(   latWat_till(msoil, mtill, mLC))
     352        1656 :        allocate(COSMIC_L3_till(msoil, mtill, mLC))
     353        1242 :        allocate(        latWat(msoil, mHor      ))
     354        1242 :        allocate(     COSMIC_L3(msoil, mHor      ))
     355     1835262 :        latWat_till    = 0.000001_dp
     356     1835262 :        COSMIC_L3_till = 0.000001_dp
     357     1834020 :        COSMIC_L3      = 0.000001_dp
     358     1834020 :        latWat         = 0.000001_dp
     359             : 
     360             : 
     361             : 
     362             :        ! earlier these variables were allocated with  size(soilId0,1)
     363             :        ! in which the variable "soilId0" changes according to the iFlag_soilDB
     364             :        ! so better to use other variable which is common to both soilDB (0 AND 1) flags
     365        1656 :        allocate(KsVar_H0(size(Id0, 1)))
     366         828 :        allocate(KsVar_V0(size(Id0, 1)))
     367         828 :        allocate( SMs_FC0(size(Id0, 1)))
     368             : 
     369         822 :        select case(processMatrix(3,1))
     370             :        case(1)
     371             :           ! first thirteen parameters go to this routine
     372         408 :           iStart = processMatrix(3, 3) - processMatrix(3, 2) + 1
     373         408 :           iEnd = processMatrix(3, 3) - 4
     374             : 
     375             :           ! next four parameters go here
     376             :           ! (the first three for the fRoots and the fourth one for the beta)
     377         408 :           iStart2 = processMatrix(3, 3) - 4 + 1
     378         408 :           iEnd2 = processMatrix(3, 3)
     379             : 
     380             :        case(2)
     381             :           ! first thirteen parameters go to this routine
     382           2 :           iStart = processMatrix(3, 3) - processMatrix(3, 2) + 1
     383           2 :           iEnd = processMatrix(3, 3) - 5
     384             : 
     385             :           ! next four parameters go here
     386             :           ! (the first three for the fRoots and the fourth one for the beta)
     387           2 :           iStart2 = processMatrix(3, 3) - 5 + 1
     388           2 :           iEnd2 = processMatrix(3, 3) - 1
     389             :           ! last parameter is jarvis parameter - no need to be regionalized
     390          74 :           jarvis_thresh_c1 = param(processMatrix(3, 3))
     391             :           !write(*,*) 'jarvis_thresh_c1 = ', jarvis_thresh_c1
     392             :           !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
     393             : 
     394             :        case(3)
     395             :           ! first thirteen parameters go to this routine
     396           2 :           iStart = processMatrix(3, 3) - processMatrix(3, 2) + 1
     397           2 :           iEnd = processMatrix(3, 3) - 9
     398             : 
     399             :           ! next four parameters go here
     400             :           ! (the first three for the fRoots and the fourth one for the beta)
     401           2 :           iStart2 = processMatrix(3, 3) - 8
     402           2 :           iEnd2 = processMatrix(3, 3) - 1
     403             : 
     404             :           ! last parameter is jarvis parameter - no need to be regionalized
     405          74 :           jarvis_thresh_c1 = param(processMatrix(3, 3))
     406             : 
     407             :           !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
     408             :           !write(*,*) 'jarvis_thresh_c1 = ', jarvis_thresh_c1
     409             : 
     410             :        case(4)
     411             :           ! first thirteen parameters go to this routine
     412           2 :           iStart = processMatrix(3, 3) - processMatrix(3, 2) + 1
     413           2 :           iEnd = processMatrix(3, 3) - 8
     414             : 
     415             :           ! next four parameters go here
     416             :           ! (the first three for the fRoots and the fourth one for the beta)
     417           2 :           iStart2 = processMatrix(3, 3) - 7
     418           2 :           iEnd2 = processMatrix(3, 3)
     419             :           !write(*,*) 'iStart, iEnd, iStart2, iEnd2 = ', iStart, iEnd, iStart2, iEnd2
     420             : 
     421             :        case DEFAULT
     422             :           call error_message('***ERROR: Process description for process "soil moisture parametrization"', &
     423         414 :                'does not exist! mo_multi_param_reg')
     424             :        end select
     425             : 
     426           0 :        call mpr_sm(param(iStart : iEnd), processMatrix, &
     427             :             soilDB%is_present, soilDB%nHorizons, soilDB%nTillHorizons, &
     428             :             soilDB%sand, soilDB%clay, soilDB%DbM, &
     429             :             Id0, soilId0, LCover0(:, iiLC), &
     430             :             thetaS_till, thetaFC_till, thetaPW_till, thetaS, &
     431         414 :             thetaFC, thetaPW, Ks, Db, KsVar_H0, KsVar_V0, SMs_FC0)
     432             : 
     433             :        ! >> neutron count related parameters
     434         414 :        if ( processMatrix(10,1) .GT. 0 ) &
     435             :             call mpr_neutrons( processMatrix(10,1), &  ! IN: processmatrix case
     436           0 :             param( processMatrix(10,3)-processMatrix(10,2)+1:processMatrix(10,3) ) , & ! IN:  global parameter set
     437             :             soilDB%is_present       , & ! IN:  flag indicating presence of soil
     438             :             soilDB%nHorizons        , & ! IN:  Number of Horizons of Soiltype
     439             :             soilDB%nTillHorizons    , & ! IN:  Number of tillage Horizons
     440             :             LCover0(:, iiLC)        , & ! IN:  land cover ids at level 0
     441             :             soilDB%clay             , & ! IN:  clay content
     442             :             soilDB%DbM              , & ! IN:  mineral Bulk density
     443             :             Db                      , & ! IN: Bulk density
     444             :             COSMIC_L3_till          , & ! OUT: COSMIC parameter L3 tillage layer
     445             :             latWat_till             , & ! OUT: COSMIC parameter Lattice Water tillage layer
     446             :             COSMIC_L3               , & ! OUT: COSMIC parameter L3
     447             :             latWat                    & ! OUT: COSMIC parameter Lattice Water
     448          12 :             )
     449             : 
     450           0 :        call mpr_SMhorizons(param(iStart2:iEnd2), processMatrix, &
     451             :             iFlag_soilDB, nSoilHorizons_mHM, HorizonDepth_mHM, &
     452             :             LCover0(:, iiLC), soilId0, &
     453             :             soilDB%nHorizons, soilDB%nTillHorizons, &
     454             :             thetaS_till, thetaFC_till, thetaPW_till, &
     455             :             thetaS, thetaFC, thetaPW, &
     456             :             soilDB%Wd, Db, soilDB%DbM, soilDB%RZdepth, &
     457             :             mask0, Id0, &
     458             :             upper_bound1, lower_bound1, left_bound1, right_bound1, n_subcells1, &
     459             :             soilMoistExp1(:, :, iiLC), soilMoistSat1(:, :, iiLC), soilMoistFC1(:, :, iiLC), &
     460             :             wiltingPoint1(:, :, iiLC), fRoots1(:, :, iiLC), &
     461             :             ! >>>>>> neutron count
     462             :             latWat_till, COSMIC_L3_till, latWat, COSMIC_L3, &
     463             :             bulkDens1(:,:,iiLC), latticeWater1(:,:,iiLC), COSMICL31(:,:,iiLC) &
     464         414 :             )
     465             : 
     466         414 :        deallocate(thetaS_till)
     467         414 :        deallocate(thetaFC_till)
     468         414 :        deallocate(thetaPW_till)
     469         414 :        deallocate(thetaS)
     470         414 :        deallocate(thetaFC)
     471         414 :        deallocate(thetaPW)
     472         414 :        deallocate(Ks)
     473         414 :        deallocate(Db)
     474             : 
     475             :        ! neutron count
     476         414 :        deallocate( latWat_till    )
     477         414 :        deallocate( COSMIC_L3_till )
     478         414 :        deallocate( latWat     )
     479         414 :        deallocate( COSMIC_L3  )
     480             : 
     481             :        ! ------------------------------------------------------------------
     482             :        ! potential evapotranspiration (PET)
     483             :        ! ------------------------------------------------------------------
     484             :        ! Penman-Monteith method is only method that is LCscene dependent
     485         414 :        if (processMatrix(5, 1) == 3) then
     486          12 :           iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
     487          12 :           iEnd = processMatrix(5, 3)
     488           0 :           call aerodynamical_resistance(gridded_LAI0, LCover0(:, iiLC), param(iStart : iEnd - 1), mask0, &
     489             :                Id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
     490          12 :                aeroResist1(:, :, iiLC))
     491         402 :        else if (processMatrix(5, 1) == -1) then
     492           4 :           iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
     493           4 :           iEnd = processMatrix(5, 3)
     494             : 
     495           0 :           call pet_correctbyLAI(param(iStart : iEnd), nodata_dp, &
     496             :                LCover0(:, iiLC), gridded_LAI0, mask0, Id0, &
     497             :                upper_bound1, lower_bound1, left_bound1, &
     498           4 :                right_bound1, n_subcells1, petLAIcorFactor(:, :, iiLC))
     499             :        end if
     500             : 
     501             :        ! ------------------------------------------------------------------
     502             :        ! interflow
     503             :        ! ------------------------------------------------------------------
     504         828 :        select case(processMatrix(6, 1))
     505             :        case (1)
     506             :           !
     507         414 :           iStart = processMatrix(6, 3) - processMatrix(6, 2) + 1
     508         414 :           iEnd = processMatrix(6, 3)
     509             :           ! TODO: this subroutine should be split into each param (or at least extract kFastFlow1)
     510             :           ! because it is in the loop unnecessarily
     511             :           call mpr_runoff(LCover0(:, iiLC), mask0, SMs_FC0, slope_emp0, &
     512           0 :                KsVar_H0, param(iStart : iEnd), Id0, upper_bound1, lower_bound1, &
     513             :                left_bound1, right_bound1, n_subcells1, unsatThresh1(:, 1, 1), kFastFlow1(:, 1, iiLC), &
     514         414 :                kSlowFlow1(:, 1, iiLC), alpha1(:, 1, iiLC))
     515             :        case DEFAULT
     516         414 :           call error_message('***ERROR: Process description for process "interflow" does not exist! mo_multi_param_reg')
     517             :        END select
     518             : 
     519             :        ! ------------------------------------------------------------------
     520             :        ! percolation cofficient, karstic percolation loss
     521             :        ! ------------------------------------------------------------------
     522         828 :        select case(processMatrix(7, 1))
     523             :        case(1)
     524             : 
     525         414 :           iStart = processMatrix(7, 3) - processMatrix(7, 2) + 1
     526         414 :           iEnd = processMatrix(7, 3)
     527             :           call karstic_layer(& ! In
     528           0 :                param(iStart : iEnd), & ! In
     529             :                geoUnit0, mask0, & ! In
     530             :                SMs_FC0, KsVar_V0, Id0, & ! In
     531             :                n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, & ! In
     532             :                karstLoss1(:, 1, 1), kPerco1(:, 1, iiLC)                            & ! Out
     533         414 :                )
     534             : 
     535             :        case DEFAULT
     536         414 :           call error_message('***ERROR: Process description for process "percolation" does not exist! mo_multi_param_reg')
     537             :        end select
     538             : 
     539         414 :        deallocate(KsVar_H0)
     540         414 :        deallocate(KsVar_V0)
     541         621 :        deallocate(SMs_FC0)
     542             : 
     543             :     end do !! >>>>>>> LAND COVER SCENE LOOP
     544             : 
     545             : 
     546             :     ! ------------------------------------------------------------------
     547             :     ! sealed area threshold for runoff generation
     548             :     ! ------------------------------------------------------------------
     549         414 :     select case(processMatrix(4, 1))
     550             :     case (1)
     551         207 :        iStart = processMatrix(4, 3) - processMatrix(4, 2) + 1
     552         207 :        iEnd = processMatrix(4, 3)
     553         207 :        call iper_thres_runoff(param(iStart : iEnd), sealedThresh1)
     554             :     case DEFAULT
     555         207 :        call error_message('***ERROR: Process description for process "runoff_generation" does not exist! mo_multi_param_reg')
     556             :     end select
     557             : 
     558             :     ! ------------------------------------------------------------------
     559             :     ! potential evapotranspiration (PET)
     560             :     ! ------------------------------------------------------------------
     561         207 :     select case(processMatrix(5, 1))
     562             :     case(-1) ! LAI correction of input PET
     563         404 :        iEnd = -9999 ! dummy statement
     564             :     case(0) ! aspect correction of input PET
     565         197 :        iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
     566         197 :        iEnd = processMatrix(5, 3)
     567         197 :        call pet_correctbyASP(Id0, y0, Asp0, param(iStart : iEnd), nodata_dp, fAsp0)
     568         197 :        fAsp1(:, 1, 1) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
     569        7662 :             left_bound1, right_bound1, Id0, mask0, nodata_dp, fAsp0)
     570             :     case(1) ! Hargreaves-Samani method
     571           1 :        iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
     572           1 :        iEnd = processMatrix(5, 3)
     573           1 :        call pet_correctbyASP(Id0, y0, Asp0, param(iStart : iEnd - 1), nodata_dp, fAsp0)
     574           1 :        fAsp1(:, 1, 1) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
     575          36 :             left_bound1, right_bound1, Id0, mask0, nodata_dp, fAsp0)
     576          37 :        HarSamCoeff1 = param(iEnd)
     577             :     case(2) ! Priestley-Taylor Method
     578           1 :        iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
     579           1 :        iEnd = processMatrix(5, 3)
     580           0 :        call priestley_taylor_alpha(gridded_LAI0, param(iStart : iEnd), &
     581             :             mask0, nodata_dp, Id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
     582           1 :             PrieTayAlpha1(:, :, 1))
     583             :     case(3) ! Penman-Monteith method
     584             :        ! aerodynamic resistance is calculated inside LCscene loop
     585           6 :        iStart = processMatrix(5, 3) - processMatrix(5, 2) + 1
     586           6 :        iEnd = processMatrix(5, 3)
     587           0 :        call bulksurface_resistance(gridded_LAI0, param(iEnd), mask0, &
     588             :             nodata_dp, Id0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, &
     589           6 :             surfResist1(:, :, 1))
     590             :     case default
     591         207 :        call error_message('***ERROR: Process description for process "pet correction" does not exist! mo_multi_param_reg')
     592             :     end select
     593             :     ! ------------------------------------------------------------------
     594             :     ! baseflow recession parameter
     595             :     ! ------------------------------------------------------------------
     596         414 :     select case(processMatrix(9, 1))
     597             :       case(1)
     598             :         ! the number of process parameters, so the number in processMatrix(9,2) has
     599             :         ! to be equal to the size of geo_unit_list
     600         207 :         iStart = processMatrix(9, 3) - processMatrix(9, 2) + 1
     601         207 :         iEnd = processMatrix(9, 3)
     602             : 
     603         207 :         call baseflow_param(param(iStart : iEnd), geoUnit0, k2_0)
     604             : 
     605             :         ! Upscale by arithmetic mean
     606         621 :         allocate(k2_1(size(kBaseFlow1, 1)))
     607             :         k2_1 = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
     608         207 :             left_bound1, right_bound1, Id0, mask0, nodata_dp, k2_0)
     609             :         ! loop over all LCover scenes
     610         621 :         do iiLC = 1, size(LCover0, 2)
     611       15837 :           kBaseFlow1(:, 1, iiLC) = k2_1
     612             :         end do
     613         207 :         deallocate(k2_1)
     614             : 
     615             :         ! correction and unit conversion
     616             :         ! if percolation is ON: correct K2 such that it is at least k1
     617             :         ! since kSlowFlow1 is LCover dependent, kBaseFlow1 is too
     618       16458 :         if (processMatrix(7, 1) .gt. 0) kBaseFlow1 = merge(kSlowFlow1, kBaseFlow1, kBaseFlow1 .lt. kSlowFlow1)
     619             : 
     620             :       case DEFAULT
     621         207 :         call error_message('***ERROR: Process description for process "baseflow Recession" does not exist! mo_multi_param_reg')
     622             :     end select
     623             : 
     624             :     ! ------------------------------------------------------------------
     625             :     ! Neutron count related parameters
     626             :     ! >> only N0 parameter - others are defined above in soil parameters
     627             :     ! ------------------------------------------------------------------
     628         213 :     select case(processMatrix(10, 1))
     629             :     case(0)
     630             :        ! do nothing
     631             :     case(1)
     632             :        ! the number of process parameters, so the number in processMatrix(9,2) has
     633           6 :        iStart = processMatrix(10, 3) - processMatrix(10, 2) + 1
     634           6 :        iEnd = processMatrix(10, 3)
     635         672 :        No_Count1 = param(iStart)  ! >> 1st parameter --> N0 parameter
     636             :     case(2)
     637             :        ! the number of process parameters, so the number in processMatrix(9,2) has
     638           0 :        iStart = processMatrix(10, 3) - processMatrix(10, 2) + 1
     639           0 :        iEnd = processMatrix(10, 3)
     640           0 :        No_Count1 = param(iStart)  ! >> 1st parameter --> N0 parameter
     641             :     case DEFAULT
     642         207 :        call error_message('***ERROR: Process description for process "Neutron count" does not exist! mo_multi_param_reg')
     643             :     end select
     644             : 
     645             : 
     646             :     !-------------------------------------------------------------------
     647             :     ! call regionalization of parameters related to LAI
     648             :     ! it is now outside of mHM since LAI is now dynamic variable
     649             :     !-------------------------------------------------------------------
     650             :     call canopy_intercept_param(processMatrix, param(:), &
     651             :          gridded_LAI0, n_subcells1, upper_bound1, lower_bound1, left_bound1, right_bound1, Id0, mask0, &
     652         207 :          nodata_dp, maxInter1(:, :, 1))
     653             : 
     654         414 :   end subroutine mpr
     655             : 
     656             :   ! ----------------------------------------------------------------------------
     657             : 
     658             :   !    NAME
     659             :   !        baseflow_param
     660             : 
     661             :   !    PURPOSE
     662             :   !>       \brief baseflow recession parameter
     663             : 
     664             :   !>       \details This subroutine calculates the baseflow recession parameter
     665             :   !>       based on the geological units at the Level 0 scale. For each level 0
     666             :   !>       cell, it assigns the value specified in the parameter array param for the
     667             :   !>       geological unit in this cell.
     668             :   !>       Global parameters needed (see mhm_parameter.nml):
     669             :   !>       - param(1) = GeoParam(1,:)
     670             :   !>       - param(2) = GeoParam(2,:)
     671             :   !>       - ...
     672             : 
     673             :   !    INTENT(IN)
     674             :   !>       \param[in] "real(dp), dimension(:) :: param"       list of required parameters
     675             :   !>       \param[in] "integer(i4), dimension(:) :: geoUnit0" ids of geological units at L0
     676             : 
     677             :   !    INTENT(OUT)
     678             :   !>       \param[out] "real(dp), dimension(:) :: k2_0" - baseflow recession parameter at Level 0
     679             : 
     680             :   !    HISTORY
     681             :   !>       \authors Stephan Thober, Rohini Kumar
     682             : 
     683             :   !>       \date Dec 2012
     684             : 
     685             :   ! Modifications:
     686             :   ! Stephan Thober  Dec 2013 - changed intent(inout) to intent(out)
     687             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     688             : 
     689         414 :   subroutine baseflow_param(param, geoUnit0, k2_0)
     690             : 
     691         207 :     use mo_mpr_global_variables, only : GeoUnitList
     692             :     !$ use omp_lib
     693             : 
     694             :     implicit none
     695             : 
     696             :     ! list of required parameters
     697             :     real(dp), dimension(:), intent(in) :: param
     698             : 
     699             :     ! ids of geological units at L0
     700             :     integer(i4), dimension(:), intent(in) :: geoUnit0
     701             : 
     702             :     ! - baseflow recession parameter at Level 0
     703             :     real(dp), dimension(:), intent(out) :: k2_0
     704             : 
     705             :     ! loop variable
     706             :     integer(i4) :: ii
     707             : 
     708             :     ! geo unit
     709             :     integer(i4), dimension(1) :: gg
     710             : 
     711             : 
     712         207 :     if (size(param) .ne. size(geoUnitList)) &
     713           0 :           call error_message(' mo_multi_param_reg: baseflow_param: size mismatch, subroutine baseflow parameters ')
     714             : 
     715     9583982 :     k2_0 = nodata_dp
     716             : 
     717             :     !$OMP PARALLEL
     718             :     !$OMP DO PRIVATE(gg) SCHEDULE(STATIC)
     719     9583982 :     do ii = 1, size(k2_0)
     720             :        ! get parameter index in geoUnitList
     721   105421525 :        gg = minloc(abs(geoUnitList - geoUnit0(ii)))
     722     9583982 :        k2_0(ii) = param(gg(1))
     723             :     end do
     724             :     !$OMP END DO
     725             :     !$OMP END PARALLEL
     726             : 
     727         207 :   end subroutine baseflow_param
     728             : 
     729             :   ! ----------------------------------------------------------------------------
     730             : 
     731             :   !    NAME
     732             :   !        snow_acc_melt_param
     733             : 
     734             :   !    PURPOSE
     735             :   !>       \brief Calculates the snow parameters.
     736             : 
     737             :   !>       \details This subroutine calculates the snow parameters
     738             :   !>       threshold temperature (TT), degree-day factor without precipitation (DD)
     739             :   !>       and maximum degree-day factor (DDmax) as well as increase of degree-day
     740             :   !>       factor per mm of increase in precipitation (IDDP).
     741             : 
     742             :   !>       Global parameters needed (see mhm_parameter.nml):
     743             :   !>       - param(1) = snowTreshholdTemperature
     744             :   !>       - param(2) = degreeDayFactor_forest
     745             :   !>       - param(3) = degreeDayFactor_impervious
     746             :   !>       - param(4) = degreeDayFactor_pervious
     747             :   !>       - param(5) = increaseDegreeDayFactorByPrecip
     748             :   !>       - param(6) = maxDegreeDayFactor_forest
     749             :   !>       - param(7) = maxDegreeDayFactor_impervious
     750             :   !>       - param(8) = maxDegreeDayFactor_pervious
     751             :   !>       INTENT(IN)
     752             : 
     753             :   !    INTENT(IN)
     754             :   !>       \param[in] "real(dp), dimension(8) :: param"    eight global parameters
     755             :   !>       \param[in] "real(dp), dimension(:) :: fForest1" [1] fraction of forest cover
     756             :   !>       \param[in] "real(dp), dimension(:) :: fIperm1"  [1] fraction of sealed area
     757             :   !>       \param[in] "real(dp), dimension(:) :: fPerm1"   [1] fraction of permeable area
     758             : 
     759             :   !    INTENT(OUT)
     760             :   !>       \param[out] "real(dp), dimension(:) :: tempThresh1"  [degreeC] threshold temperature for snow rain
     761             :   !>       \param[out] "real(dp), dimension(:) :: degDayNoPre1" [mm-1 degreeC-1] Degree-day factor with
     762             :   !>       \param[out] "real(dp), dimension(:) :: degDayInc1"   [d-1 degreeC-1]  Increase of the Degree-day
     763             :   !>       \param[out] "real(dp), dimension(:) :: degDayMax1"   [mm-1 degreeC-1] Maximum Degree-day factor
     764             : 
     765             :   !    HISTORY
     766             :   !>       \authors Stephan Thober, Rohini Kumar
     767             : 
     768             :   !>       \date Dec 2012
     769             : 
     770             :   ! Modifications:
     771             :   ! Juliane Mai    Oct 2013 - OLD parametrization --> param(1) = snowTreshholdTemperature
     772             :   !                                                --> param(2) = degreeDayFactor_forest
     773             :   !                                                --> param(3) = degreeDayFactor_impervious
     774             :   !                                                --> param(4) = degreeDayFactor_pervious
     775             :   !                                                --> param(5) = increaseDegreeDayFactorByPrecip
     776             :   !                                                --> param(6) = maxDegreeDayFactor_forest
     777             :   !                                                --> param(7) = maxDegreeDayFactor_impervious
     778             :   !                                                --> param(8) = maxDegreeDayFactor_pervious
     779             :   !                                             -------------------------------
     780             :   !                                             degreeDayFactor_impervious    = degreeDayFactor_forest + delta_1 + delta_2
     781             :   !                                             degreeDayFactor_pervious      = degreeDayFactor_forest + delta_1
     782             :   !                                             maxDegreeDayFactor_forest     = degreeDayFactor_forest + delta_3
     783             :   !                                             maxDegreeDayFactor_impervious = degreeDayFactor_impervious + delta_5
     784             :   !                                                                           = degreeDayFactor_forest + delta_1 + delta_2 + delta_5
     785             :   !                                             maxDegreeDayFactor_pervious   = degreeDayFactor_pervious + delta_4
     786             :   !                                                                           = degreeDayFactor_forest + delta_1 + delta_4
     787             :   !                                             -------------------------------
     788             :   !                                             NEW parametrization
     789             :   !                                                --> param(1) = snowTreshholdTemperature
     790             :   !                                                --> param(2) = degreeDayFactor_forest
     791             :   !                                                --> param(3) = delta_2
     792             :   !                                                --> param(4) = delta_1
     793             :   !                                                --> param(5) = increaseDegreeDayFactorByPrecip
     794             :   !                                                --> param(6) = delta_3
     795             :   !                                                --> param(7) = delta_5
     796             :   !                                                --> param(8) = delta_4
     797             :   ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
     798             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     799             : 
     800         828 :   subroutine snow_acc_melt_param(param, fForest1, fIperm1, fPerm1, tempThresh1, degDayNoPre1, degDayInc1, &
     801         414 :        degDayMax1)
     802             :     implicit none
     803             : 
     804             :     ! eight global parameters
     805             :     real(dp), dimension(8), intent(in) :: param
     806             : 
     807             :     ! [1] fraction of forest cover
     808             :     real(dp), dimension(:), intent(in) :: fForest1
     809             : 
     810             :     ! [1] fraction of sealed area
     811             :     real(dp), dimension(:), intent(in) :: fIperm1
     812             : 
     813             :     ! [1] fraction of permeable area
     814             :     real(dp), dimension(:), intent(in) :: fPerm1
     815             : 
     816             :     ! [degreeC] threshold temperature for snow rain
     817             :     real(dp), dimension(:), intent(out) :: tempThresh1
     818             : 
     819             :     ! [mm-1 degreeC-1] Degree-day factor with
     820             :     real(dp), dimension(:), intent(out) :: degDayNoPre1
     821             : 
     822             :     ! [mm-1 degreeC-1] Maximum Degree-day factor
     823             :     real(dp), dimension(:), intent(out) :: degDayMax1
     824             : 
     825             :     ! [d-1 degreeC-1]  Increase of the Degree-day
     826             :     real(dp), dimension(:), intent(out) :: degDayInc1
     827             : 
     828         414 :     real(dp) :: tmp_degreeDayFactor_forest, tmp_degreeDayFactor_impervious, tmp_degreeDayFactor_pervious
     829             : 
     830         414 :     real(dp) :: tmp_maxDegreeDayFactor_forest, tmp_maxDegreeDayFactor_impervious, tmp_maxDegreeDayFactor_pervious
     831             : 
     832             : 
     833         414 :     tmp_degreeDayFactor_forest = param(2)                                    ! OLD: param(2)
     834         414 :     tmp_degreeDayFactor_impervious = param(2) + param(4) + param(3)              ! OLD: param(3)
     835         414 :     tmp_degreeDayFactor_pervious = param(2) + param(4)                         ! OLD: param(4)
     836         414 :     tmp_maxDegreeDayFactor_forest = param(2) + param(6)   ! OLD: param(6)
     837         414 :     tmp_maxDegreeDayFactor_impervious = param(2) + param(4) + param(3) + param(7)   ! OLD: param(7)
     838         414 :     tmp_maxDegreeDayFactor_pervious = param(2) + param(4) + param(8)   ! OLD: param(8)
     839             : 
     840       15630 :     tempThresh1 = param(1)
     841       15630 :     degDayInc1 = param(5)
     842             : 
     843           0 :     degDayNoPre1 = (&
     844             :          tmp_degreeDayFactor_forest * fForest1 + &
     845           0 :          tmp_degreeDayFactor_impervious * fIperm1 + &
     846       16044 :          tmp_degreeDayFactor_pervious * fPerm1)
     847           0 :     degDayMax1 = (&
     848             :          tmp_maxDegreeDayFactor_forest * fForest1 + &
     849             :          tmp_maxDegreeDayFactor_impervious * fIperm1 + &
     850       16044 :          tmp_maxDegreeDayFactor_pervious * fPerm1)
     851             : 
     852         207 :   end subroutine snow_acc_melt_param
     853             : 
     854             :   ! ----------------------------------------------------------------------------
     855             : 
     856             :   !    NAME
     857             :   !        iper_thres_runoff
     858             : 
     859             :   !    PURPOSE
     860             :   !>       \brief sets the impervious layer threshold parameter for runoff generation
     861             : 
     862             :   !>       \details to be done by Kumar
     863             :   !>       ....
     864             : 
     865             :   !>       Global parameters needed (see mhm_parameter.nml):
     866             :   !>       - param(1) = imperviousStorageCapacity
     867             : 
     868             :   !    INTENT(IN)
     869             :   !>       \param[in] "real(dp), dimension(1) :: param" - given threshold parameter
     870             : 
     871             :   !    INTENT(OUT)
     872             :   !>       \param[out] "real(dp), dimension(:, :, :) :: sealedThresh1"
     873             : 
     874             :   !    HISTORY
     875             :   !>       \authors Stephan Thober, Rohini Kumar
     876             : 
     877             :   !>       \date Dec 2012
     878             : 
     879             :   ! Modifications:
     880             :   ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
     881             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     882             : 
     883         207 :   subroutine iper_thres_runoff(param, sealedThresh1)
     884             :     implicit none
     885             : 
     886             :     ! - given threshold parameter
     887             :     real(dp), dimension(1), intent(in) :: param
     888             : 
     889             :     real(dp), dimension(:, :, :), intent(out) :: sealedThresh1
     890             : 
     891             : 
     892        8229 :     sealedThresh1 = param(1)
     893             : 
     894         414 :   end subroutine iper_thres_runoff
     895             : 
     896             :   ! ----------------------------------------------------------------------------
     897             : 
     898             :   !    NAME
     899             :   !        karstic_layer
     900             : 
     901             :   !    PURPOSE
     902             :   !>       \brief calculates the Karstic percolation loss
     903             : 
     904             :   !>       \details This subroutine calls first the karstic_fraction upscaling
     905             :   !>       routine for determine the karstic fraction area for every Level 1
     906             :   !>       cell. Then, the karstic percolation loss is estimated given two
     907             :   !>       shape parameters by
     908             :   !>       \f[ karstLoss1 = 1 + ( fKarArea * param(1)) *( (-1)**INT(param(2),i4) ) \f]
     909             :   !>       where \f$ karstLoss1 \f$ is the karstic percolation loss and \f$ fKarArea \f$
     910             :   !>       is the fraction of karstic area at level 1
     911             :   !>       Global parameters needed (see mhm_parameter.nml):
     912             :   !>       - param(1) = rechargeCoefficient
     913             :   !>       - param(2) = rechargeFactor_karstic
     914             :   !>       - param(3) = gain_loss_GWreservoir_karstic
     915             : 
     916             :   !    INTENT(IN)
     917             :   !>       \param[in] "real(dp), dimension(3) :: param"           parameters
     918             :   !>       \param[in] "integer(i4), dimension(:) :: geoUnit0"     id of the Karstic formation
     919             :   !>       \param[in] "logical, dimension(:, :) :: mask0"         mask at level 0
     920             :   !>       \param[in] "real(dp), dimension(:) :: SMs_FC0"         [-] soil mositure deficit from field
     921             :   !>       capacity w.r.t to saturation
     922             :   !>       \param[in] "real(dp), dimension(:) :: KsVar_V0"        [-] relative variability of saturated
     923             :   !>       \param[in] "integer(i4), dimension(:) :: Id0"          Cell ids of hi res field
     924             :   !>       \param[in] "integer(i4), dimension(:) :: n_subcells1"  number of l0 cells within a l1 cell
     925             :   !>       \param[in] "integer(i4), dimension(:) :: upper_bound1" upper row of a l1 cell in l0 grid
     926             :   !>       \param[in] "integer(i4), dimension(:) :: lower_bound1" lower row of a l1 cell in l0 grid
     927             :   !>       \param[in] "integer(i4), dimension(:) :: left_bound1"  left col of a l1 cell in l0 grid
     928             :   !>       \param[in] "integer(i4), dimension(:) :: right_bound1" right col of a l1 cell in l0 grid
     929             : 
     930             :   !    INTENT(OUT)
     931             :   !>       \param[out] "real(dp), dimension(:) :: karstLoss1" [-]    Karstic percolation loss
     932             :   !>       \param[out] "real(dp), dimension(:) :: L1_Kp"      [d-1] percolation coefficient
     933             : 
     934             :   !    HISTORY
     935             :   !>       \authors Rohini Kumar, Stephan Thober
     936             : 
     937             :   !>       \date Feb 2013
     938             : 
     939             :   ! Modifications:
     940             :   ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
     941             :   ! Stephan Thober Dec 2013 - changed intent(inout) to intent(out)
     942             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
     943             : 
     944        1242 :   subroutine karstic_layer(param, geoUnit0, mask0, SMs_FC0, KsVar_V0, Id0, n_subcells1, upper_bound1, lower_bound1, &
     945         828 :        left_bound1, right_bound1, karstLoss1, L1_Kp)
     946             : 
     947         207 :     use mo_mpr_global_variables, only : GeoUnitList, geoUnitKar
     948             :     use mo_upscaling_operators, only : L0_fractionalCover_in_Lx, upscale_arithmetic_mean
     949             :     !$ use omp_lib
     950             : 
     951             :     implicit none
     952             : 
     953             :     ! parameters
     954             :     real(dp), dimension(3), intent(in) :: param
     955             : 
     956             :     ! id of the Karstic formation
     957             :     integer(i4), dimension(:), intent(in) :: geoUnit0
     958             : 
     959             :     ! mask at level 0
     960             :     logical, dimension(:, :), intent(in) :: mask0
     961             : 
     962             :     ! [-] soil mositure deficit from field
     963             :     ! capacity w.r.t to saturation
     964             :     real(dp), dimension(:), intent(in) :: SMs_FC0
     965             : 
     966             :     ! [-] relative variability of saturated
     967             :     real(dp), dimension(:), intent(in) :: KsVar_V0
     968             : 
     969             :     ! Cell ids of hi res field
     970             :     integer(i4), dimension(:), intent(in) :: Id0
     971             : 
     972             :     ! number of l0 cells within a l1 cell
     973             :     integer(i4), dimension(:), intent(in) :: n_subcells1
     974             : 
     975             :     ! upper row of a l1 cell in l0 grid
     976             :     integer(i4), dimension(:), intent(in) :: upper_bound1
     977             : 
     978             :     ! lower row of a l1 cell in l0 grid
     979             :     integer(i4), dimension(:), intent(in) :: lower_bound1
     980             : 
     981             :     ! left col of a l1 cell in l0 grid
     982             :     integer(i4), dimension(:), intent(in) :: left_bound1
     983             : 
     984             :     ! right col of a l1 cell in l0 grid
     985             :     integer(i4), dimension(:), intent(in) :: right_bound1
     986             : 
     987             :     ! [-]    Karstic percolation loss
     988             :     real(dp), dimension(:), intent(out) :: karstLoss1
     989             : 
     990             :     ! [d-1] percolation coefficient
     991             :     real(dp), dimension(:), intent(out) :: L1_Kp
     992             : 
     993             :     ! fraction of karstic area
     994         414 :     real(dp), dimension(:), allocatable :: fKarArea
     995             : 
     996             :     ! temporal variable
     997    19168378 :     real(dp), dimension(size(SMs_FC0, 1)) :: tmp
     998             : 
     999             :     integer(i4) :: nGeoUnits
    1000             : 
    1001             :     integer(i4) :: i
    1002             : 
    1003             : 
    1004             :     ! ------------------------------------------------------------------
    1005             :     ! PERCOLATION; 1/Kp = f(Ks)
    1006             :     ! ------------------------------------------------------------------
    1007             :     ! Regionalise Kp with variability of last soil layer property
    1008             :     !$OMP PARALLEL
    1009         414 :     tmp = merge(param(1) * (1.0_dp + SMs_FC0) / (1.0_dp + KsVar_V0), &
    1010    19168378 :          nodata_dp, Id0 .ne. nodata_i4)
    1011             :     !$OMP END PARALLEL
    1012             : 
    1013             :     L1_Kp = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
    1014         414 :          left_bound1, right_bound1, Id0, mask0, nodata_dp, tmp)
    1015             : 
    1016             :     ! minimum constrains
    1017       15630 :     L1_Kp = merge(2.0_dp, L1_Kp, L1_Kp .lt. 2.0_dp)
    1018             : 
    1019         414 :     nGeoUnits = size(geoUnitlist, 1)
    1020             : 
    1021             :     ! 1st calculate fraction of Karstic area
    1022        1242 :     allocate(fKarArea(size(karstLoss1, 1)))
    1023       16044 :     fKarArea = 0.0_dp
    1024             : 
    1025        4554 :     do i = 1, nGeoUnits
    1026        4140 :        if(GeoUnitKar(i) .eq. 0) cycle
    1027           0 :        fKarArea(:) = L0_fractionalCover_in_Lx(geoUnit0, geoUnitlist(i), mask0, &
    1028        4554 :             upper_bound1, lower_bound1, left_bound1, right_bound1, n_subcells1)
    1029             :     end do
    1030             : 
    1031             :     ! 2nd calculate karstLoss1
    1032       15630 :     karstLoss1 = 1.0_dp - (fKarArea * param(2))
    1033             : 
    1034         414 :     deallocate(fKarArea)
    1035             : 
    1036         414 :   end subroutine karstic_layer
    1037             : 
    1038             :   ! ----------------------------------------------------------------------------
    1039             : 
    1040             :   !    NAME
    1041             :   !        canopy_intercept_param
    1042             : 
    1043             :   !    PURPOSE
    1044             :   !>       \brief estimate effective maximum interception capacity at L1
    1045             : 
    1046             :   !>       \details estimate effective maximum interception capacity at L1 for a given
    1047             :   !>       Leaf Area Index field.
    1048             :   !>       Global parameters needed (see mhm_parameter.nml):
    1049             :   !>       Process Case 1:
    1050             :   !>       - param(1) = canopyInterceptionFactor
    1051             : 
    1052             :   !    INTENT(IN)
    1053             :   !>       \param[in] "integer(i4), dimension(:, :) :: processMatrix" indicate processes
    1054             :   !>       \param[in] "real(dp), dimension(:) :: param"               array of global parameters
    1055             :   !>       \param[in] "real(dp), dimension(:, :) :: LAI0"             LAI at level-0(nCells0, time)
    1056             :   !>       \param[in] "integer(i4), dimension(:) :: n_subcells1"      Number of L0 cells within a L1 cell
    1057             :   !>       \param[in] "integer(i4), dimension(:) :: upper_bound1"     Upper row of high resolution block
    1058             :   !>       \param[in] "integer(i4), dimension(:) :: lower_bound1"     Lower row of high resolution block
    1059             :   !>       \param[in] "integer(i4), dimension(:) :: left_bound1"      Left column of high resolution block
    1060             :   !>       \param[in] "integer(i4), dimension(:) :: right_bound1"     Right column of high resolution block
    1061             :   !>       \param[in] "integer(i4), dimension(:) :: Id0"              Cell ids at level 0
    1062             :   !>       \param[in] "logical, dimension(:, :) :: mask0"             mask at level 0 field
    1063             :   !>       \param[in] "real(dp) :: nodata"                            - nodata value
    1064             : 
    1065             :   !    INTENT(OUT)
    1066             :   !>       \param[out] "real(dp), dimension(:, :) :: max_intercept1" max interception at level-1(nCells1, time)
    1067             : 
    1068             :   !    HISTORY
    1069             :   !>       \authors Rohini Kumar
    1070             : 
    1071             :   !>       \date Aug. 2013
    1072             : 
    1073             :   ! Modifications:
    1074             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
    1075             : 
    1076         621 :   subroutine canopy_intercept_param(processMatrix, param, LAI0, n_subcells1, upper_bound1, lower_bound1, left_bound1, &
    1077         621 :        right_bound1, Id0, mask0, nodata, max_intercept1)
    1078             : 
    1079         414 :     use mo_string_utils, only : num2str
    1080             :     use mo_upscaling_operators, only : upscale_arithmetic_mean
    1081             : 
    1082             :     implicit none
    1083             : 
    1084             :     ! indicate processes
    1085             :     integer(i4), dimension(:, :), intent(in) :: processMatrix
    1086             : 
    1087             :     ! array of global parameters
    1088             :     real(dp), dimension(:), intent(in) :: param
    1089             : 
    1090             :     ! LAI at level-0(nCells0, time)
    1091             :     real(dp), dimension(:, :), intent(in) :: LAI0
    1092             : 
    1093             :     ! Number of L0 cells within a L1 cell
    1094             :     integer(i4), dimension(:), intent(in) :: n_subcells1
    1095             : 
    1096             :     ! Upper row of high resolution block
    1097             :     integer(i4), dimension(:), intent(in) :: upper_bound1
    1098             : 
    1099             :     ! Lower row of high resolution block
    1100             :     integer(i4), dimension(:), intent(in) :: lower_bound1
    1101             : 
    1102             :     ! Left column of high resolution block
    1103             :     integer(i4), dimension(:), intent(in) :: left_bound1
    1104             : 
    1105             :     ! Right column of high resolution block
    1106             :     integer(i4), dimension(:), intent(in) :: right_bound1
    1107             : 
    1108             :     ! Cell ids at level 0
    1109             :     integer(i4), dimension(:), intent(in) :: Id0
    1110             : 
    1111             :     ! mask at level 0 field
    1112             :     logical, dimension(:, :), intent(in) :: mask0
    1113             : 
    1114             :     ! - nodata value
    1115             :     real(dp), intent(in) :: nodata
    1116             : 
    1117             :     ! max interception at level-1(nCells1, time)
    1118             :     real(dp), dimension(:, :), intent(out) :: max_intercept1
    1119             : 
    1120             :     integer(i4) :: iStart, iEnd, it
    1121             : 
    1122         207 :     real(dp), dimension(:), allocatable :: max_intercept0
    1123             : 
    1124         207 :     real(dp), dimension(:), allocatable :: gamma_intercept
    1125             : 
    1126             : 
    1127             :     ! ------------------------------------------------------------------
    1128             :     ! Maximum interception parameter
    1129             :     ! ------------------------------------------------------------------
    1130         207 :     select case(processMatrix(1, 1))
    1131             :     case(1)
    1132         207 :        iStart = processMatrix(1, 3) - processMatrix(1, 2) + 1
    1133         207 :        iEnd = processMatrix(1, 3)
    1134             : 
    1135             :        ! allocate space
    1136         621 :        allocate(gamma_intercept(iEnd - iStart + 1))
    1137         621 :        allocate(max_intercept0 (size(Id0, 1)))
    1138             : 
    1139             :        ! estimate max. intercept at Level-0
    1140         414 :        gamma_intercept(:) = param(iStart : iEnd)
    1141             : 
    1142        2709 :        do it = 1, size(LAI0, 2)
    1143             :           !$OMP PARALLEL
    1144   115845612 :           max_intercept0(:) = LAI0(:, it) * gamma_intercept(1)
    1145             :           !$OMP END PARALLEL
    1146             : 
    1147             :           ! Upscale by arithmetic mean
    1148             :           max_intercept1(:, it) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, left_bound1, &
    1149        2709 :                right_bound1, Id0, mask0, nodata, max_intercept0(:))
    1150             : 
    1151             :        end do
    1152             : 
    1153         207 :        deallocate(gamma_intercept)
    1154         207 :        deallocate(max_intercept0)
    1155             :     CASE DEFAULT
    1156         207 :        call error_message('mo_multi_param_reg: This processMatrix=', num2str(processMatrix(1, 1)), ' is not implemented!')
    1157             :     end select
    1158             : 
    1159         207 :   end subroutine canopy_intercept_param
    1160             : 
    1161             : 
    1162             :   ! ----------------------------------------------------------------------------
    1163             : 
    1164             :   !    NAME
    1165             :   !        aerodynamical_resistance
    1166             : 
    1167             :   !    PURPOSE
    1168             :   !>       \brief Regionalization of aerodynamic resistance
    1169             : 
    1170             :   !>       \details estimation of aerodynamical resistance
    1171             :   !>       Global parameters needed (see mhm_parameter.nml):
    1172             :   !>       - param(1) = canopyheigth_forest
    1173             :   !>       - param(2) = canopyheigth_impervious
    1174             :   !>       - param(3) = canopyheigth_pervious
    1175             :   !>       - param(4) = displacementheight_coeff
    1176             :   !>       - param(5) = roughnesslength_momentum_coeff
    1177             :   !>       - param(6) = roughnesslength_heat_coeff
    1178             : 
    1179             :   !    INTENT(IN)
    1180             :   !>       \param[in] "real(dp), dimension(:, :) :: LAI0"         LAI at level-0
    1181             :   !>       \param[in] "integer(i4), dimension(:) :: LCover0"      land cover field
    1182             :   !>       \param[in] "real(dp), dimension(6) :: param"           input parameter
    1183             :   !>       \param[in] "logical, dimension(:, :) :: mask0"         mask at level 0
    1184             :   !>       \param[in] "integer(i4), dimension(:) :: Id0"          Cell ids of hi res field
    1185             :   !>       \param[in] "integer(i4), dimension(:) :: n_subcells1"  number of l0 cells within a l1 cell
    1186             :   !>       \param[in] "integer(i4), dimension(:) :: upper_bound1" upper row of a l1 cell in l0 grid
    1187             :   !>       \param[in] "integer(i4), dimension(:) :: lower_bound1" lower row of a l1 cell in l0 grid
    1188             :   !>       \param[in] "integer(i4), dimension(:) :: left_bound1"  left col of a l1 cell in l0 grid
    1189             :   !>       \param[in] "integer(i4), dimension(:) :: right_bound1" right col of a l1 cell in l0 grid
    1190             : 
    1191             :   !    INTENT(OUT)
    1192             :   !>       \param[out] "real(dp), dimension(:, :) :: aerodyn_resistance1" aerodynmaical resistance
    1193             : 
    1194             :   !    HISTORY
    1195             :   !>       \authors Matthias Zink
    1196             : 
    1197             :   !>       \date Apr 2013
    1198             : 
    1199             :   ! Modifications:
    1200             :   ! Matthias Zink Jun 2017 - moved from mo_multi_scale_param_reg.f90 to mo_mpr_pet.f90
    1201             :   ! Robert Schweppe Jun 2018 - refactoring and reformatting
    1202             : 
    1203          36 :   subroutine aerodynamical_resistance(LAI0, LCover0, param, mask0, Id0, n_subcells1, upper_bound1, lower_bound1, &
    1204          24 :        left_bound1, right_bound1, aerodyn_resistance1)
    1205             : 
    1206         207 :     use mo_common_constants, only : eps_dp
    1207             :     use mo_mpr_constants, only : WindMeasHeight, karman
    1208             :     use mo_upscaling_operators, only : upscale_arithmetic_mean
    1209             : 
    1210             :     implicit none
    1211             : 
    1212             :     ! LAI at level-0
    1213             :     real(dp), dimension(:, :), intent(in) :: LAI0
    1214             : 
    1215             :     ! land cover field
    1216             :     integer(i4), dimension(:), intent(in) :: LCover0
    1217             : 
    1218             :     ! input parameter
    1219             :     real(dp), dimension(6), intent(in) :: param
    1220             : 
    1221             :     ! mask at level 0
    1222             :     logical, dimension(:, :), intent(in) :: mask0
    1223             : 
    1224             :     ! Cell ids of hi res field
    1225             :     integer(i4), dimension(:), intent(in) :: Id0
    1226             : 
    1227             :     ! number of l0 cells within a l1 cell
    1228             :     integer(i4), dimension(:), intent(in) :: n_subcells1
    1229             : 
    1230             :     ! upper row of a l1 cell in l0 grid
    1231             :     integer(i4), dimension(:), intent(in) :: upper_bound1
    1232             : 
    1233             :     ! lower row of a l1 cell in l0 grid
    1234             :     integer(i4), dimension(:), intent(in) :: lower_bound1
    1235             : 
    1236             :     ! left col of a l1 cell in l0 grid
    1237             :     integer(i4), dimension(:), intent(in) :: left_bound1
    1238             : 
    1239             :     ! right col of a l1 cell in l0 grid
    1240             :     integer(i4), dimension(:), intent(in) :: right_bound1
    1241             : 
    1242             :     ! aerodynmaical resistance
    1243             :     real(dp), dimension(:, :), intent(out) :: aerodyn_resistance1
    1244             : 
    1245             :     integer(i4) :: tt
    1246             : 
    1247          12 :     real(dp), dimension(:), allocatable :: maxLAI
    1248             : 
    1249          12 :     real(dp), dimension(:), allocatable :: zm
    1250             : 
    1251          12 :     real(dp), dimension(:), allocatable :: canopy_height0
    1252             : 
    1253          12 :     real(dp), dimension(:), allocatable :: zm_zero, zh_zero, displace
    1254             : 
    1255             :     ! dim 1 = number of cells on level 0,
    1256             :     ! dim2=month of year
    1257          12 :     real(dp), dimension(:, :), allocatable :: aerodyn_resistance0
    1258             : 
    1259             : 
    1260             :     ! initialize some things
    1261      558576 :     allocate(zm                  (size(LCover0, dim = 1))) ; zm = nodata_dp
    1262      558564 :     allocate(zm_zero             (size(LCover0, dim = 1))) ; zm_zero = nodata_dp
    1263      558564 :     allocate(zh_zero             (size(LCover0, dim = 1))) ; zh_zero = nodata_dp
    1264      558564 :     allocate(displace            (size(LCover0, dim = 1))) ; displace = nodata_dp
    1265      558564 :     allocate(canopy_height0      (size(LCover0, dim = 1))) ; canopy_height0 = nodata_dp
    1266     6702672 :     allocate(aerodyn_resistance0 (size(LCover0, dim = 1), size(LAI0, 2))) ; aerodyn_resistance0 = nodata_dp
    1267      558564 :     allocate(maxLAI              (size(LCover0, dim = 1))) ; maxLAI = nodata_dp
    1268             : 
    1269             :     ! regionalization of canopy height
    1270             :     ! substitute with canopy height
    1271      558564 :     canopy_height0 = merge(param(1), canopy_height0, LCover0 == 1)  ! forest
    1272      558564 :     canopy_height0 = merge(param(2), canopy_height0, LCover0 == 2)  ! impervious
    1273             : 
    1274             :     ! old implementation used values from LUT statically for all cells (Jan-Dec, 12 values):
    1275             :     ! 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
    1276          12 :     maxLAI = MAXVAL(LAI0, dim=2)
    1277             : 
    1278         156 :     do tt = 1, size(LAI0, 2)
    1279             : 
    1280             :        ! pervious canopy height is scaled with LAI
    1281     6702768 :        canopy_height0 = merge((param(3) * LAI0(:, tt) / maxLAI), canopy_height0, LCover0 == 3)  ! pervious
    1282             : 
    1283             :        ! estimation of the aerodynamic resistance on the lower level
    1284             :        ! see FAO Irrigation and Drainage Paper No. 56 (p. 19 ff) for more information
    1285     6702624 :        zm = WindMeasHeight
    1286             :        ! correction: if wind measurement height is below canopy height loagarithm becomes negative
    1287     6702768 :        zm = merge(canopy_height0 + zm, zm, ((abs(zm - nodata_dp) .GT. eps_dp) .AND. (zm .LT. canopy_height0)))
    1288             : 
    1289             :        ! zh       = zm
    1290     6702768 :        displace = param(4) * canopy_height0
    1291     6702768 :        zm_zero = param(5) * canopy_height0
    1292     6702768 :        zh_zero = param(6) * zm_zero
    1293             :        !
    1294             :        ! calculate aerodynamic resistance (changes monthly)
    1295     6702624 :        aerodyn_resistance0(:, tt) = log((zm - displace) / zm_zero) * log((zm - displace) / zh_zero) / (karman**2.0_dp)
    1296             :        aerodyn_resistance1(:, tt) = upscale_arithmetic_mean(n_subcells1, upper_bound1, lower_bound1, &
    1297         156 :             left_bound1, right_bound1, Id0, mask0, nodata_dp, aerodyn_resistance0(:, tt))
    1298             : 
    1299             :     end do
    1300             : 
    1301          12 :   end subroutine aerodynamical_resistance
    1302             : 
    1303             : END MODULE mo_multi_param_reg

Generated by: LCOV version 1.16