!============================================================================== ! This source code is part of the ! Australian Community Atmosphere Biosphere Land Exchange (CABLE) model. ! This work is licensed under the CSIRO Open Source Software License ! Agreement (variation of the BSD / MIT License). ! ! You may not use this file except in compliance with this License. ! A copy of the License (CSIRO_BSD_MIT_License_v2.0_CABLE.txt) is located ! in each directory containing CABLE code. ! ! ============================================================================== ! Purpose: Calculate roughness lengths as a function of soil and canopy ! parameters ! ! Contact: Eva.Kowalczyk@csiro.au ! ! History: No significant changes since v1.4b except change to cope with ! split timestep in ACCESS (zref_uv, zref_tq) ! ! ! ============================================================================== MODULE cable_roughness_module USE cable_surface_types_mod, ONLY: ICE_SurfaceType => ICE_cable USE cable_phys_constants_mod, ONLY: CCSD => CSD USE cable_phys_constants_mod, ONLY: CCRD => CRD USE cable_phys_constants_mod, ONLY: CCCD => CCD USE cable_phys_constants_mod, ONLY: CCCW_C => CCW_C USE cable_phys_constants_mod, ONLY: CUSUHM => USUHM USE cable_phys_constants_mod, ONLY: CVONK => VONK USE cable_phys_constants_mod, ONLY: CA33 => A33 USE cable_phys_constants_mod, ONLY: CCTL => CTL USE cable_phys_constants_mod, ONLY: CZDLIN => ZDLIN USE cable_phys_constants_mod, ONLY: CCSW => CSW USE cable_phys_constants_mod, ONLY: CGRAV => GRAV USE cable_other_constants_mod, ONLY : CLAI_THRESH => LAI_THRESH !*# Overview ! ! The procedures contained in this module calculate the roughness parameters ! and the aerodynamic contribution to the resistances controlling the fluxes ! of momentum, heat and water vapour between the land and atmosphere for each ! land point. ! ! The formulations take into account vegetation and snow cover. The dependence ! on atmospheric conditions (i.e. the surface heat fluxes) is incorporated ! later within the [[define_canopy]] subroutine. IMPLICIT NONE real, parameter :: z0soilsn_min = 1.e-7 !! Minimum value for the roughness length for bare soil, \(z_{0,min}\) (m) real, parameter :: z0soilsn_min_PF = 1.e-4 !! Minimum value for the roughness length for permanent ice on land, !! \(z_{0,minPF}\) (m) PRIVATE PUBLIC ruff_resist CONTAINS SUBROUTINE ruff_resist(veg, rough, ssnow, canopy, LAI_pft, HGT_pft, reducedLAIdue2snow ) !* Calculates the roughness parameters and the aerodynamic contribution to the ! resistances controlling the fluxes between the land and atmosphere for each ! land point. ! !## Scientific description ! ! The scientific basis for the formulae is a combination of Localised Near ! Field theory for aerodynamic transfer within canopies (for the heat and ! water vapour fluxes) and bulk formulae for the roughness length and ! displacement height accounting for roughness sublayer effects (for the ! momentum flux). ! ! The calculations take into account whether: ! ! 1. the soil model used is SLI or the default (soilsnow) ! 2. the land point is vegetated (LAI > LAI_THRESH) or not ! ! The primary references for these formulations are ! ! * [Raupach MR (1989)](https://doi.org/10.1002/qj.49711548710) ! * [Raupach MR (1992)](https://doi.org/10.1007/BF00155203) ! * [Raupach MR (1994)](https://doi.org/10.1007/BF00709229) ! * [Kowalczyk et al. (2006)](http://www.cmar.csiro.au/e-print/open/kowalczykea_2006a.pdf) ! !## Outputs ! ! The principal outputs from the MODULE are ! ! * `rough%zref_tq`: reference height above the displacement height for the ! air temperature and humidity, where these forcing observations are deemed ! to have been collected (m) ! * `rough%zref_uv`: height above the displacement height for the ! wind speeds, where the forcing observations are deemed to have been ! collected (m) ! * `rough%hruff`: the canopy height accounting for the presence of snow (m) ! * `canopy%vlaiw`: the leaf area index accounting for the presence of snow ! (m\(^2\) m\(^{-2}\)) ! * `rough%z0soil`: the aerodynamic roughness length for soil (m) ! * `rough%z0ssoilsn`: the aerodynamic roughness length for snow (m) ! * `rough%z0m`: the aerodynamic roughness length for the surface ! (canopy+soil+snow) (m) ! * `rough%disp`: the displacement height of the surface ! (=0.0 if not vegetated) (m) ! * `rough%zruffs`: the depth of the roughness sublayer over vegetated ! surfaces (m) ! * `rough%coexp`: the extinction coefficient for the wind speed profile ! within the canopy, \(c_{0}\) (m\(^{-1}\)). ! ! \(c_{0}\) is defined to be the coefficient within an expontial wind ! profile, i.e. the wind speed at height \(z\) above the ground within ! a canopy of height \(h_c\) is given by ! \( U(z) = U_{h} \exp\{ c_{0} (z-h_c) \} \) where \(U_h\) is the wind ! speed at canopy top. ! ! * `rough%usuh`: the ratio of the friction velocity, \(u_*\), to the wind ! speed at canopy top (-) ! * `rough%rt0us`: *normalised* aerodynamic resistance for the turbulent ! transfer from the soil/snow surface to the displacement height (-) ! * `rough%rt1us`: *normalised* aerodynamic resistance for the turbulent ! transfer from the displacement height to the reference level (-) ! ! `rough%rt1us` is evaluated in three subparts (`rough%rt1usa`, ! `rough%rt1usb`, and `rough%rt1usc`). One of the resistance terms ! (`rough%rt1usc`) is evaluated in subroutine [[define_canopy]]. ! ! Each of the normalized resistances are given by the theoretical formulae ! given by the references. The aerodynamic resistances for the current ! time step are evaluated later by dividing the *normalized resistances* by ! the current time step's friction velocity `canopy%us`. USE cable_common_module, ONLY : cable_user, cable_runtime USE cable_def_types_mod, ONLY : veg_parameter_type, roughness_type, & soil_snow_type, canopy_type, mp !subrs USE hruff_eff_LAI_mod_cbl, ONLY : HgtAboveSnow USE hruff_eff_LAI_mod_cbl, ONLY : LAI_eff !data USE cable_other_constants_mod, ONLY : z0surf_min implicit none !!## Structure !result returned from called subr. Avail. in cross_*paths_module - but unsure !yet real :: HeightAboveSnow(mp) real :: reducedLAIdue2snow(mp) TYPE(roughness_type), INTENT(INOUT) :: rough TYPE (canopy_type), INTENT(INOUT) :: canopy TYPE(soil_snow_type), INTENT(IN) :: ssnow TYPE (veg_parameter_type), INTENT(INOUT) :: veg real :: LAI_pft(mp) real :: HGT_pft(mp) REAL, DIMENSION(mp) :: & xx, & ! =CCCD*LAI; working variable dh ! d/h where d is zero-plane displacement integer :: i ! Set canopy height above snow level: call HgtAboveSnow( HeightAboveSnow, mp, z0soilsn_min, veg%hc, ssnow%snowd, & ssnow%ssdnn ) !* * evaluates the canopy height and leaf area given the presence of snow ! (or not) using [[HgtAboveSnow]] and [[LAI_eff]] rough%hruff = HeightAboveSnow !why are we using veg% for LAI and height-although this is synced now IF(cable_runtime%um_radiation .OR. cable_runtime%esm15 & .OR. cable_runtime%offline ) THEN ! LAI decreases due to snow: formerly canopy%vlaiw CALL LAI_eff( mp, veg%vlai, veg%hc, HeightAboveSnow, & reducedLAIdue2snow ) canopy%vlaiw = reducedLAIdue2snow ENDIF canopy%rghlai = canopy%vlaiw !* * sets the value of soil and snow roughness lengths ! (depends on the configuration of CABLE) IF (cable_user%soil_struc=='default') THEN ! Roughness length of bare soil (m): new formulation- E.Kowalczyk 2014 IF (.NOT.cable_user%l_new_roughness_soil .AND. (.NOT.cable_user%or_evap)) THEN rough%z0soil = 0.0009*min(1.0,canopy%vlaiw) + 1.e-4 rough%z0soilsn = rough%z0soil ELSE rough%z0soil = 0.01*min(1.0,canopy%vlaiw) + 0.02*min(canopy%us**2/CGRAV,1.0) rough%z0soilsn = max(1.e-7,rough%z0soil) ENDIF WHERE( ssnow%snowd .GT. 0.01 ) & rough%z0soilsn = MAX(z0soilsn_min, & rough%z0soil - rough%z0soil*MIN(ssnow%snowd,10.)/10.) WHERE( ssnow%snowd .GT. 0.01 .AND. veg%iveg == ICE_SurfaceType ) & rough%z0soilsn = MAX(rough%z0soilsn, z0soilsn_min_PF ) ELSEIF (cable_user%soil_struc=='sli') THEN rough%z0soil = 0.01*MIN(1.0,canopy%vlaiw) + 0.02*MIN(canopy%us**2/CGRAV,1.0) rough%z0soilsn = MAX(1.e-2,rough%z0soil) ! (1e-2: Mori et al., J Ag Met, 2010) WHERE( ssnow%snowd .GT. 0.01 ) & rough%z0soilsn = MAX( 1.e-2, & rough%z0soil - rough%z0soil*MIN(ssnow%snowd,10.)/10.) ENDIF !| * evaluates the remaining output variables, depending on whether the land ! point is vegetated or not. do i=1,mp if( canopy%vlaiw(i) .LE. CLAI_THRESH .OR. & rough%hruff(i) .LT. rough%z0soilsn(i) ) then ! BARE SOIL SURFACE rough%z0m(i) = rough%z0soilsn(i) rough%rt0us(i) = 0.0 rough%disp(i) = 0.0 ! Reference height zref is height above the displacement height ! Ticket #148: Reference height is height above the displacement height ! noting that this must be above the roughness length and rough%hruff-rough%disp ! (though second case is unlikely to be attained) rough%zref_uv(i) = MAX( 3.5 + rough%z0m(i), rough%za_uv(i) ) rough%zref_tq(i) = MAX( 3.5 + rough%z0m(i), rough%za_tq(i) ) rough%zref_uv(i) = MAX( rough%zref_uv(i), rough%hruff(i)-rough%disp(i) ) rough%zref_tq(i) = MAX( rough%zref_tq(i), rough%hruff(i)-rough%disp(i) ) rough%zruffs(i) = 0.0 rough%rt1usa(i) = 0.0 rough%rt1usb(i) = 0.0 ! Friction velocity/windspeed at canopy height ! eq. 7 Raupach 1994, BLM, vol 71, p211-216 ! (CUSUHM set in physical_constants module): rough%usuh(i) = MIN( SQRT( CCSD + CCRD * ( canopy%vlaiw(i) * 0.5 ) ), CUSUHM ) xx(i) = SQRT( CCCD * MAX( ( canopy%vlaiw(i) * 0.5 ), 0.0005 ) ) ! Displacement height/canopy height: ! eq.8 Raupach 1994, BLM, vol 71, p211-216 dh(i) = 1.0 - ( 1.0 - EXP( -xx(i) ) ) / xx(i) ! Extinction coefficient for wind profile in canopy: ! eq. 3.14, SCAM manual (CSIRO tech report 132) rough%coexp(i) = rough%usuh(i) / ( CVONK * CCCW_C * ( 1.0 - dh(i) ) ) ELSE ! VEGETATED SURFACE ! Friction velocity/windspeed at canopy height ! eq. 7 Raupach 1994, BLM, vol 71, p211-216 ! (CUSUHM set in physical_constants module): rough%usuh(i) = MIN( SQRT( CCSD + CCRD * ( canopy%rghlai(i) * 0.5 ) ), & CUSUHM ) xx(i) = SQRT( CCCD * MAX( ( canopy%rghlai(i) * 0.5 ), 0.0005 ) ) ! eq.8 Raupach 1994, BLM, vol 71, p211-216: dh(i) = 1.0 - ( 1.0 - EXP( -xx(i) ) ) / xx(i) ! Calculate zero-plane displacement: rough%disp(i) = dh(i)* rough%hruff(i) ! Calculate roughness length: rough%z0m(i) = ( (1.0 - dh(i)) * EXP( LOG( CCCW_C ) - 1. + 1. / CCCW_C & - CVONK / rough%usuh(i) ) ) * rough%hruff(i) ! Reference height zref is height above the displacement height ! Reference height zref is height above the displacement height ! Ticket #148: Reference height is height above the displacement height ! noting that this must be above the roughness length and rough%hruff-rough%disp rough%zref_uv(i) = MAX( 3.5 + rough%z0m(i), rough%za_uv(i) ) rough%zref_tq(i) = MAX( 3.5 + rough%z0m(i), rough%za_tq(i) ) rough%zref_uv(i) = MAX( rough%zref_uv(i), rough%hruff(i)-rough%disp(i) ) rough%zref_tq(i) = MAX( rough%zref_tq(i), rough%hruff(i)-rough%disp(i) ) ! find coexp: see notes "simplified wind model ..." eq 34a ! Extinction coefficient for wind profile in canopy: ! eq. 3.14, SCAM manual (CSIRO tech report 132) rough%coexp(i) = rough%usuh(i) / ( CVONK * CCCW_C * ( 1.0 - dh(i) ) ) rough%term2(i) = EXP( 2 * CCSW * canopy%rghlai(i) * & ( 1 - rough%disp(i) / rough%hruff(i) ) ) rough%term3(i) = CA33**2 * CCTL * 2 * CCSW * canopy%rghlai(i) rough%term5(i) = MAX( ( 2. / 3. ) * rough%hruff(i) / rough%disp(i), 1.0 ) rough%term6(i) = EXP( 3. * rough%coexp(i) * ( rough%disp(i) / rough%hruff(i) -1. ) ) rough%term6a(i) = EXP(rough%coexp(i) * ( 0.1 * rough%hruff(i) / rough%hruff(i) -1. )) ! eq. 3.54, SCAM manual (CSIRO tech report 132) rough%rt0us(i) = rough%term5(i) * ( CZDLIN * LOG( & CZDLIN * rough%disp(i) / rough%z0soilsn(i) ) + & ( 1 - CZDLIN ) ) & * ( EXP( 2 * CCSW * canopy%rghlai(i) ) - rough%term2(i) ) & / rough%term3(i) ! See CSIRO SCAM, Raupach et al 1997, eq. 3.49: rough%zruffs(i) = rough%disp(i) + rough%hruff(i) * CA33**2 * CCTL / CVONK / & rough%term5(i) ! See CSIRO SCAM, Raupach et al 1997, eq. 3.51: rough%rt1usa(i) = rough%term5(i) * ( rough%term2(i) - 1.0 ) / rough%term3(i) rough%rt1usb(i) = rough%term5(i) * ( MIN( rough%zref_tq(i) + rough%disp(i), & rough%zruffs(i) ) - rough%hruff(i) ) / & ( CA33**2 * CCTL * rough%hruff(i) ) rough%rt1usb(i) = MAX( rough%rt1usb(i), 0.0 ) ! in case zrufs < rough%hruff END IF END DO !> * if the SLI soil model is used - updates the evaluated `rough%rt0us` IF (cable_user%soil_struc.EQ.'sli') THEN WHERE( canopy%vlaiw .GE. CLAI_THRESH .AND. & rough%hruff .GE. rough%z0soilsn ) ! VEGETATED SURFACE ! vh ! Haverd et al., Biogeosciences 10, 2011-2040, 2013 rough%rt0us = LOG( rough%disp / (0.1 * rough%hruff) ) * & EXP( 2. * CCSW * canopy%rghlai ) * rough%disp & / rough%hruff / (Ca33 ** 2 * Cctl) ENDWHERE ENDIF END SUBROUTINE ruff_resist END MODULE cable_roughness_module