dicussion needed on calculation of aerodynamic resistance

If PET is calculated using Penman-Monteith, aerodynamic resistance needs to be calculated. This is happening in subroutine aerodynamic_resistance in mo_multi_param_reg.f90. It was implemented by Matthias Zink and cites a FAO publication as a source. Now here is the problem: The FAO paper uses a calculation based on a paper of Jacobs(2002) that seems to be valid for the grass reference surface. What I found in a brief search in the internet is that it is most commonly used for canopy_heights for crops of up to 2m. However, we also use that function for forest canopies. The canopy height for forests in mHM is controlled by a parameter (param(1)) and can range from 15 (default) to 40m. The measurement height of wind measurement (zm) is a parameter and set to 10m. To make that function work, the canopy_height0 is added to zm , if the canopy_height0 is higher than zm. I do not know, if that is really the way the function is supposed to be used.

Here is the relevant code:

    ! 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

So my questions are (before I adapt this to the new MPR) @lese, @rakovec, @rkumar, @thober:

• are we okay with using this function for calculation of aerodynamic resistance for forests?
• if yes, is the correction, like Matthias did it, okay, or do we come up with a smoother solution, because now: canopy_height: 9.99 m -> zm = 10.0 m, canopy_height: 10.0 m -> zm = 20.0 m
• if no, what other transfer function do we want?