MODULE cbl_friction_vel_module !* This MODULE contains the SUBROUTINE [[comp_friction_vel]] and two ! FUNCTIONS ([[psim]] and [[psis]]) needed to ! evaluate the friction velocity over each land point/tile ! given the wind speed and the current estimate of the Monin-Obukhov ! stability parameter \(\xi\). Outputs `friction_vel`, `psim` and ! `psis` are used in [[define_canopy]]. IMPLICIT NONE PUBLIC comp_friction_vel PUBLIC psim PUBLIC psis PRIVATE CONTAINS SUBROUTINE comp_friction_vel(friction_vel, iter, mp, CVONK, CUMIN, CPI_C, & zetar, zref_uv, zref_tq, z0m, ua ) !*## Purpose ! ! This SUBROUTINE evaluates the value of the friction velocity \(u_*\) as used ! during the iteration loop of the Monin-Obukhov (MO) similarity theory in ! [[define_canopy]]. \(u_*\) is used when quantifying the stability ! parameter `canopy%zetar` in [[update_zetar]], and when evaluating the ! resistance network components `rt0` and `rt1` in [[define_canopy]]. ! The friction velocity quantifies the magnitude of turbulent mixing that ! is occurring as well as the flux of momentum from the atmosphere to ! the land. ! ! The output `friction_vel` is known as `canopy%us` elsewhere in the code. ! ! !## Method ! ! The basic formula is derived from the modified log-law for the ! vertical profile of wind speed near to the ground namely ! ! \( u_{*} = U(z_{ref}) / ( \log [z_{ref}/z_{0m}] - ! \psi_m[\xi * z_{ref}/z_{refTq}] ! + \psi_m[\xi * z_{0m}/z_{refTq}] ) \) ! ! with \(\psi_m\) the integrated similarity function given by [[psim]], ! \(U\) the wind speed at height \(z_{ref}\), ! \(\xi\) the current value of the stability parameter, and \(z_{0m}\) the ! roughness length. ! A minimum value is applied to the input wind speed, `ua`= \(U(z_{ref})\), ! to assist with convergence of the MO iteration in light wind conditions. ! Small and large value limits are applied to the evaluated \(u_*\). ! ! **NOTE** Most literature references assume that the reference levels ! for wind and temperature are equal inside the arguments to the ! stability functions \(\psi_m\) and \(\psi_s\) - in CABLE this need ! not be true. ! \(\xi\) is always defined relative to the reference level for temperature, ! in \(\psi_m\) the reference height for wind is needed. ! Consequently the calls to [[psim]] include conversion ! factors to account for the different reference levels. ! `zref_uv` and `zref_tq` passed to comp_friction_vel are the heights above ! the displacement height `rough%disp`. ! !## References !  ! - section 3.1, equations 1-9. IMPLICIT NONE INTEGER, INTENT(IN) :: mp !! size of cable vector of land points (-) REAL, INTENT(OUT) :: friction_vel(mp) !! friction velocity (ms\(^{-1}\)) INTEGER, INTENT(IN) :: iter !! MO iteration counter (-) ! physical constants REAL, INTENT(IN) :: CVONK !! von Karman constant (-) REAL, INTENT(IN):: CUMIN !! minimum value of wind speed (ms\(^{-1}\)) ! maths & other constants REAL, INTENT(IN) :: CPI_C !! \(\pi\) (-) REAL, INTENT(IN) :: zetar(mp,iter) !! stability parameter - see [[update_zetar]] `canopy%zetar` (-) REAL, INTENT(IN) :: zref_uv(mp) !! reference height for wind `rough%zref_uv` (m) REAL, INTENT(IN) :: zref_tq(mp) !! reference height for temperature and humidity `rough%zref_tq` (m) REAL, INTENT(IN) :: z0m(mp) !! roughness length `rough%z0m` (m) REAL, INTENT(IN) :: ua(mp) !! wind speed `met%ua` (ms\(^{-1}\)) !local vars REAL :: lower_limit(mp), rescale(mp) REAL :: psim_1(mp), psim_2(mp), psim_arg(mp) REAL :: z_eff(mp) !INH: Ticket #138 %us is defined based on U(rough%zref_uv) ! but zetar based on rough%zref_tq - changes to ensure consistency !NB no RSL incorporated here psim_1 = psim( zetar(:,iter) * zref_uv/zref_tq, mp, CPI_C ) !rescale = CVONK * MAX( ua, SPREAD(CUMIN,1,mp ) ) rescale = CVONK * MAX( ua, CUMIN ) z_eff = zref_uv / z0m psim_arg = zetar(:,iter) * z0m / zref_tq psim_2 = psim( psim_arg, mp, CPI_C ) lower_limit = rescale / ( LOG(z_eff) - psim_1 + psim_2 ) friction_vel= MIN( MAX(1.e-6, lower_limit), 10.0 ) RETURN END SUBROUTINE comp_friction_vel FUNCTION psim(zeta, mp, CPI_C ) RESULT(r) !*## Purpose ! ! Evaluates the integrated similarity function for momentum transfer, ! \(\psi_m(\xi)\). ! !## Method ! ! Uses the Businger-Dyer form for unstable conditions (\(\xi<0\)) and the ! Beljaars-Holtslag form for stable conditions (\(\xi>0\)). ! ! This function is used in the evaluation of the friction velocity in ! [[comp_friction_vel]]. ! !## References ! - [Beljaars and Holtslag (1991)](https://doi.org/10.1175/1520-0450(1991)030<0327:FPOLSF>2.0.CO;2) ! - [Businger et al. (1971)](https://doi.org/10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2) ! - [Dyer (1974)](https://doi.org/10.1007/BF00240838) ! mrr, 16-sep-92 (from function psi: mrr, edinburgh 1977) ! computes integrated stability function psim(z/l) (z/l=zeta) ! for momentum, using the businger-dyer form for unstable cases ! and the Beljaars and Holtslag (1991) form for stable cases. INTEGER, INTENT(IN) :: mp !! size of cable vector of land points (-) REAL, INTENT(IN), DIMENSION(mp) :: zeta !! current value of stability parameter \(\xi\) (-) REAL, INTENT(IN) :: CPI_C !! \( \pi \) (-) ! function result REAL, DIMENSION(mp) :: r ! result of FUNCTION REAL, PARAMETER :: & gu = 16.0, & ! gs = 5.0, & ! a = 1.0, & ! b = 0.667, & ! xc = 5.0, & ! d = 0.35 ! REAL, DIMENSION(mp) :: & x, & ! z, & ! stable, & ! unstable ! z = 0.5 + SIGN(0.5,zeta) ! z=1 in stable, 0 in unstable ! Beljaars and Holtslag (1991) for stable stable = -a*zeta - b*(zeta - xc/d)*EXP( -d*zeta) - b*xc/d x = (1.0 + gu*ABS(zeta))**0.25 unstable = ALOG((1.0+x*x)*(1.0+x)**2/8) - 2.0*ATAN(x) + CPI_C*0.5 r = z*stable + (1.0-z)*unstable END FUNCTION psim ELEMENTAL FUNCTION psis(zeta) RESULT(r) !*## Purpose ! ! Evaluates the integrated similarity function for turbulent transfer ! of scalars, \(\psi_s(\xi)\). ! !## Method ! ! Uses the Businger-Dyer form for unstable conditions (\(\xi<0\)) and the ! Beljaars-Holtslag form for stable conditions (\(\xi>0\)). ! ! This function is used in the evaluation of the resistance network component ! `rt1usc` in [[define_canopy]]. ! !## References ! - [Beljaars and Holtslag (1991)](https://doi.org/10.1175/1520-0450(1991)030<0327:FPOLSF>2.0.CO;2) ! - [Businger et al. (1971)](https://doi.org/10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2) ! - [Dyer (1974)](https://doi.org/10.1007/BF00240838) ! mrr, 16-sep-92 (from function psi: mrr, edinburgh 1977) ! computes integrated stability function psis(z/l) (z/l=zeta) ! for scalars, using the businger-dyer form for unstable cases ! and the webb form for stable cases. see paulson (1970). REAL, INTENT(IN) :: zeta !! current value of stability parameter \(\xi\) (-) REAL, PARAMETER :: & gu = 16.0, & ! gs = 5.0, & ! a = 1.0, & ! b = 0.667, & ! c = 5.0, & ! d = 0.35 REAL :: & r, & !! OUT \(\psi_s (\xi) \) stzeta, & ! ustzeta, & ! z, & ! y, & ! stable, & ! unstable z = 0.5 + SIGN(0.5,zeta) ! z=1 in stable, 0 in unstable ! Beljaars and Holtslag (1991) for stable stzeta = MAX(0.,zeta) stable = -(1.+2./3.*a*stzeta)**(3./2.) - & b*(stzeta-c/d)*EXP(-d*stzeta) - b*c/d + 1. y = (1.0 + gu*ABS(zeta))**0.5 unstable = 2.0 * alog((1+y)*0.5) r = z*stable + (1.0-z)*unstable END FUNCTION psis END MODULE cbl_friction_vel_module