Merge branch 'mpaiao-pr-elm_outvars' into next (PR #7328)

This pull request expands the list of variables in the ELM output, and most changes were taken from CLM5.
This adds aerodynamic resistances, friction velocity, and canopy air properties to the ELM output (temperature, humidity, wind speed) to the history file.
These are helpful for assessing the model representation of turbulent fluxes. All new variables are set as "inactive" by default.

The pull request has two very minor additional edits:

Also borrowed from CLM5, this replaces the indices 1 and 2 for the resistances with local variables "above_canopy" and "below_canopy". This change is simply to improve code readability.
Revised the comments and variable output names for variable obu. I simply replaced "Monin-Obukhov length" with "Obukhov length", which is the correct name of the variable according to the American Meteorological Society.

[BFB]

Author: peterdschwartz

components/elm/src/biogeophys/BareGroundFluxesMod.F90

@@ -87,7 +87,7 @@ subroutine BareGroundFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, &
     real(r8) :: dth(bounds%begp:bounds%endp)     ! diff of virtual temp. between ref. height and surface
     real(r8) :: dthv                             ! diff of vir. poten. temp. between ref. height and surface
     real(r8) :: dqh(bounds%begp:bounds%endp)     ! diff of humidity between ref. height and surface
-    real(r8) :: obu(bounds%begp:bounds%endp)     ! Monin-Obukhov length (m)
+    real(r8) :: obu(bounds%begp:bounds%endp)     ! Obukhov length scale (m)
     real(r8) :: ur(bounds%begp:bounds%endp)      ! wind speed at reference height [m/s]
     real(r8) :: um(bounds%begp:bounds%endp)      ! wind speed including the stablity effect [m/s]
     real(r8) :: temp1(bounds%begp:bounds%endp)   ! relation for potential temperature profile
@@ -250,7 +250,7 @@ subroutine BareGroundFluxes(bounds, num_nolakeurbanp, filter_nolakeurbanp, &
          z0hg_patch(p) = z0hg_col(c)
          z0qg_patch(p) = z0qg_col(c)
 
-         ! Initialize Monin-Obukhov length and wind speed
+         ! Initialize Obukhov length scale and wind speed
 
          call MoninObukIni(ur(p), thv(c), dthv, zldis(p), z0mg_patch(p), um(p), obu(p))
          num_iter(p) = 0._r8

components/elm/src/biogeophys/CanopyFluxesMod.F90

@@ -17,6 +17,8 @@ module CanopyFluxesMod
   use elm_varctl            , only : use_hydrstress
   use elm_varpar            , only : nlevgrnd, nlevsno
   use elm_varcon            , only : namep
+  use elm_varcon            , only : mm_epsilon
+  use elm_varcon            , only : pa_to_kpa
   use pftvarcon             , only : crop, nfixer
   use decompMod             , only : bounds_type
   use PhotosynthesisMod     , only : Photosynthesis, PhotosynthesisTotal, Fractionation, PhotoSynthesisHydraulicStress
@@ -168,21 +170,16 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
     real(r8), parameter :: ria  = 0.5_r8             ! free parameter for stable formulation (currently = 0.5, "gamma" in Sakaguchi&Zeng,2008)
 
     real(r8) :: zldis(bounds%begp:bounds%endp)       ! reference height "minus" zero displacement height [m]
-    real(r8) :: zeta                                 ! dimensionless height used in Monin-Obukhov theory
     real(r8) :: wc                                   ! convective velocity [m/s]
     real(r8) :: ugust_total(bounds%begp:bounds%endp) ! gustiness including convective velocity [m/s]
     real(r8) :: dth(bounds%begp:bounds%endp)         ! diff of virtual temp. between ref. height and surface
     real(r8) :: dthv(bounds%begp:bounds%endp)        ! diff of vir. poten. temp. between ref. height and surface
     real(r8) :: dqh(bounds%begp:bounds%endp)         ! diff of humidity between ref. height and surface
-    real(r8) :: obu(bounds%begp:bounds%endp)         ! Monin-Obukhov length (m)
-    real(r8) :: um(bounds%begp:bounds%endp)          ! wind speed including the stablity effect [m/s]
     real(r8) :: ur(bounds%begp:bounds%endp)          ! wind speed at reference height [m/s]
-    real(r8) :: uaf(bounds%begp:bounds%endp)         ! velocity of air within foliage [m/s]
     real(r8) :: temp1(bounds%begp:bounds%endp)       ! relation for potential temperature profile
     real(r8) :: temp12m(bounds%begp:bounds%endp)     ! relation for potential temperature profile applied at 2-m
     real(r8) :: temp2(bounds%begp:bounds%endp)       ! relation for specific humidity profile
     real(r8) :: temp22m(bounds%begp:bounds%endp)     ! relation for specific humidity profile applied at 2-m
-    real(r8) :: ustar(bounds%begp:bounds%endp)       ! friction velocity [m/s]
     real(r8) :: tstar                                ! temperature scaling parameter
     real(r8) :: qstar                                ! moisture scaling parameter
     real(r8) :: thvstar                              ! virtual potential temperature scaling parameter
@@ -235,7 +232,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
     real(r8) :: efpot                                ! potential latent energy flux [kg/m2/s]
     real(r8) :: efe(bounds%begp:bounds%endp)         ! water flux from leaf [mm/s]
     real(r8) :: efsh                                 ! sensible heat from leaf [mm/s]
-    real(r8) :: obuold(bounds%begp:bounds%endp)      ! monin-obukhov length from previous iteration
+    real(r8) :: obuold(bounds%begp:bounds%endp)      ! Obukhov length scale from previous iteration
     real(r8) :: tlbef(bounds%begp:bounds%endp)       ! leaf temperature from previous iteration [K]
     real(r8) :: ecidif                               ! excess energies [W/m2]
     real(r8) :: err(bounds%begp:bounds%endp)         ! balance error
@@ -319,6 +316,13 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
     real(r8) :: prev_tau_diff(bounds%begp:bounds%endp) ! Previous difference in iteration tau
 
     character(len=64) :: event !! timing event
+
+    ! Indices for raw and rah
+    integer, parameter :: above_canopy = 1         ! Above canopy
+    integer, parameter :: below_canopy = 2         ! Below canopy
+
+    ! Lower bound for VPD (based on CLM)
+    real(r8), parameter :: vpd_min = 50._r8
     !------------------------------------------------------------------------------
 
     associate(                                                               &
@@ -444,6 +448,18 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
          eflx_sh_soil         => veg_ef%eflx_sh_soil        , & ! Output: [real(r8) (:)   ]  sensible heat flux from soil (W/m**2) [+ to atm]
          eflx_sh_veg          => veg_ef%eflx_sh_veg         , & ! Output: [real(r8) (:)   ]  sensible heat flux from leaves (W/m**2) [+ to atm]
          eflx_sh_grnd         => veg_ef%eflx_sh_grnd        , & ! Output: [real(r8) (:)   ]  sensible heat flux from ground (W/m**2) [+ to atm]
+         rah_above            => frictionvel_vars%rah_above_patch , & ! Output: [real(r8) (:)   ]  above-canopy sensible heat flux resistance [s/m]
+         rah_below            => frictionvel_vars%rah_above_patch , & ! Output: [real(r8) (:)   ]  below-canopy sensible heat flux resistance [s/m]
+         raw_above            => frictionvel_vars%raw_below_patch , & ! Output: [real(r8) (:)   ]  above-canopy water vapour flux resistance [s/m]
+         raw_below            => frictionvel_vars%raw_below_patch , & ! Output: [real(r8) (:)   ]  below-canopy water vapour flux resistance [s/m]
+         ustar                => frictionvel_vars%ustar_patch     , & ! Output: [real(r8) (:)   ]  friction velocity [m/s]
+         um                   => frictionvel_vars%um_patch        , & ! Output: [real(r8) (:)   ]  wind speed including the stablity effect [m/s]
+         uaf                  => frictionvel_vars%uaf_patch       , & ! Output: [real(r8) (:)   ]  canopy air wind speed [m/s]
+         taf                  => frictionvel_vars%taf_patch       , & ! Output: [real(r8) (:)   ]  canopy air temperature [K]
+         qaf                  => frictionvel_vars%qaf_patch       , & ! Output: [real(r8) (:)   ]  canopy air specific humidity [kg/kg]
+         obu                  => frictionvel_vars%obu_patch       , & ! Output: [real(r8) (:)   ]  Obukhov length scale [m]
+         zeta                 => frictionvel_vars%zeta_patch      , & ! Output: [real(r8) (:)   ]  dimensionless stability parameter 
+         vpd                  => frictionvel_vars%vpd_patch       , & ! Output: [real(r8) (:)   ]  vapour pressure deficit [kPa]
          begp                 => bounds%begp                               , &
          endp                 => bounds%endp                                 &
          )
@@ -747,7 +763,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
          p = filterp(f)
          c = veg_pp%column(p)
 
-         ! Initialize Monin-Obukhov length and wind speed
+         ! Initialize Obukhov length scale and wind speed
 
          call MoninObukIni(ur(p), thv(c), dthv(p), zldis(p), z0mv(p), um(p), obu(p))
          num_iter(p) = 0._r8
@@ -784,8 +800,8 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
 
             ! Determine aerodynamic resistances
             ram1(p)  = 1._r8/(ustar(p)*ustar(p)/um(p))
-            rah(p,1) = 1._r8/(temp1(p)*ustar(p))
-            raw(p,1) = 1._r8/(temp2(p)*ustar(p))
+            rah(p,above_canopy) = 1._r8/(temp1(p)*ustar(p))
+            raw(p,above_canopy) = 1._r8/(temp2(p)*ustar(p))
 
             ! Forbid removing more than 99% of wind speed in a time step.
             ! This is mainly to avoid convergence issues since this is such a
@@ -842,18 +858,25 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
 
             !! Sakaguchi changes for stability formulation ends here
 
-            rah(p,2) = 1._r8/(csoilcn*uaf(p))
-            raw(p,2) = rah(p,2)
+            rah(p,below_canopy) = 1._r8/(csoilcn*uaf(p))
+            raw(p,below_canopy) = rah(p,below_canopy)
             if (use_lch4) then
-               grnd_ch4_cond(p) = 1._r8/(raw(p,1)+raw(p,2))
+               grnd_ch4_cond(p) = 1._r8/(raw(p,above_canopy)+raw(p,below_canopy))
             end if
 
             ! Stomatal resistances for sunlit and shaded fractions of canopy.
             ! Done each iteration to account for differences in eah, tv.
 
-            svpts(p) = el(p)                         ! pa
-            eah(p) = forc_pbot(t) * qaf(p) / 0.622_r8   ! pa
+            svpts(p) = el(p)                         ! Pa
+            eah(p) = forc_pbot(t) * qaf(p) / mm_epsilon   ! Pa
             rhaf(p) = eah(p)/svpts(p)
+
+            ! variables for history fields
+            rah_above(p)  = rah(p,above_canopy)
+            raw_above(p)  = raw(p,above_canopy)
+            rah_below(p)  = rah(p,below_canopy)
+            raw_below(p)  = raw(p,below_canopy)
+            vpd(p)        = max((svpts(p) - eah(p)), vpd_min) * pa_to_kpa ! kPa
          end do
 
          ! Modification for shrubs proposed by X.D.Z
@@ -941,9 +964,9 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
             ! Sensible heat conductance for air, leaf and ground
             ! Moved the original subroutine in-line...
 
-            wta    = 1._r8/rah(p,1)             ! air
+            wta    = 1._r8/rah(p,above_canopy)  ! air
             wtl    = (elai(p)+esai(p))/rb(p)    ! leaf
-            wtg(p) = 1._r8/rah(p,2)             ! ground
+            wtg(p) = 1._r8/rah(p,below_canopy)  ! ground
             wtshi  = 1._r8/(wta+wtl+wtg(p))
             wtl0(p) = wtl*wtshi         ! leaf
             wtg0    = wtg(p)*wtshi      ! ground
@@ -1005,7 +1028,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
             ! Air has same conductance for both sensible and latent heat.
             ! Moved the original subroutine in-line...
 
-            wtaq    = frac_veg_nosno(p)/raw(p,1)                        ! air
+            wtaq    = frac_veg_nosno(p)/raw(p,above_canopy)             ! air
             wtlq    = frac_veg_nosno(p)*(elai(p)+esai(p))/rb(p) * rpp   ! leaf
 
             !Litter layer resistance. Added by K.Sakaguchi
@@ -1016,10 +1039,10 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
 
             ! add litter resistance and Lee and Pielke 1992 beta
             if (delq(p) < 0._r8) then  !dew. Do not apply beta for negative flux (follow old rsoil)
-               wtgq(p) = frac_veg_nosno(p)/(raw(p,2)+rdl)
+               wtgq(p) = frac_veg_nosno(p)/(raw(p,below_canopy)+rdl)
             else
                if (do_soilevap_beta()) then
-                  wtgq(p) = soilbeta(c)*frac_veg_nosno(p)/(raw(p,2)+rdl)
+                  wtgq(p) = soilbeta(c)*frac_veg_nosno(p)/(raw(p,below_canopy)+rdl)
                endif
             end if
 
@@ -1112,7 +1135,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
             taf(p) = wtg0*t_grnd(c) + wta0(p)*thm(p) + wtl0(p)*t_veg(p)
             qaf(p) = wtlq0(p)*qsatl(p) + wtgq0*qg(c) + forc_q(t)*wtaq0(p)
 
-            ! Update Monin-Obukhov length and wind speed including the
+            ! Update Obukhov length scale and wind speed including the
             ! stability effect
 
             dth(p) = thm(p)-taf(p)
@@ -1123,13 +1146,13 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
             qstar = temp2(p)*dqh(p)
 
             thvstar = tstar*(1._r8+0.61_r8*forc_q(t)) + 0.61_r8*forc_th(t)*qstar
-            zeta = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c))
+            zeta(p) = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c))
 
-            if (zeta >= 0._r8) then     !stable
-               zeta = min(2._r8,max(zeta,0.01_r8))
+            if (zeta(p) >= 0._r8) then     !stable
+               zeta(p) = min(2._r8,max(zeta(p),0.01_r8))
                um(p) = max(ur(p),0.1_r8)
             else                     !unstable
-               zeta = max(-100._r8,min(zeta,-0.01_r8))
+               zeta(p) = max(-100._r8,min(zeta(p),-0.01_r8))
                if ((.not. atm_gustiness) .or. force_land_gustiness) then
                   wc = beta*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8
                   ugust_total(p) = sqrt(ugust(t)**2 + wc**2)
@@ -1138,7 +1161,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
                   um(p) = max(ur(p),0.1_r8)
                end if
             end if
-            obu(p) = zldis(p)/zeta
+            obu(p) = zldis(p)/zeta(p)
 
             if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1
             if (nmozsgn(p) >= 4) obu(p) = zldis(p)/(-0.01_r8)

components/elm/src/biogeophys/FrictionVelocityMod.F90

@@ -24,7 +24,7 @@ module FrictionVelocityMod
   !
   ! !PUBLIC MEMBER FUNCTIONS:
   public :: FrictionVelocity       ! Calculate friction velocity
-  public :: MoninObukIni           ! Initialization of the Monin-Obukhov length
+  public :: MoninObukIni           ! Initialization of the Obukhov length scale
   !
   ! !PRIVATE MEMBER FUNCTIONS:
   private :: StabilityFunc1        ! Stability function for rib < 0.
@@ -62,7 +62,7 @@ subroutine FrictionVelocity(lbn, ubn, fn, filtern, &
     real(r8) , intent(in)    :: z0m     ( lbn: )         ! roughness length over vegetation, momentum [m] [lbn:ubn]
     real(r8) , intent(in)    :: z0h     ( lbn: )         ! roughness length over vegetation, sensible heat [m] [lbn:ubn]
     real(r8) , intent(in)    :: z0q     ( lbn: )         ! roughness length over vegetation, latent heat [m] [lbn:ubn]
-    real(r8) , intent(in)    :: obu     ( lbn: )         ! monin-obukhov length (m) [lbn:ubn]
+    real(r8) , intent(in)    :: obu     ( lbn: )         ! Obukhov length scale (m) [lbn:ubn]
     integer  , intent(in)    :: iter                     ! iteration number
     real(r8) , intent(in)    :: ur      ( lbn: )         ! wind speed at reference height [m/s] [lbn:ubn]
     real(r8) , intent(in)    :: um      ( lbn: )         ! wind speed including the stablity effect [m/s] [lbn:ubn]
@@ -445,7 +445,7 @@ end function StabilityFunc2
   subroutine MoninObukIni (ur, thv, dthv, zldis, z0m, um, obu)
     !$acc routine seq
     ! !DESCRIPTION:
-    ! Initialization of the Monin-Obukhov length.
+    ! Initialization of the Obukhov length scale.
     ! The scheme is based on the work of Zeng et al. (1998):
     ! Intercomparison of bulk aerodynamic algorithms for the computation
     ! of sea surface fluxes using TOGA CORE and TAO data. J. Climate,
@@ -462,7 +462,7 @@ subroutine MoninObukIni (ur, thv, dthv, zldis, z0m, um, obu)
     real(r8), intent(in)  :: zldis ! reference height "minus" zero displacement heght [m]
     real(r8), intent(in)  :: z0m   ! roughness length, momentum [m]
     real(r8), intent(out) :: um    ! wind speed including the stability effect [m/s]
-    real(r8), intent(out) :: obu   ! monin-obukhov length (m)
+    real(r8), intent(out) :: obu   ! Obukhov length scale (m)
     !
     ! !LOCAL VARIABLES:
     real(r8) :: wc    ! convective velocity [m/s]