diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90
index ff1ce1688c..ef1b30d2d6 100644
--- a/src/biogeophys/CanopyFluxesMod.F90
+++ b/src/biogeophys/CanopyFluxesMod.F90
@@ -228,6 +228,7 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
     use clm_varcon         , only : sb, cpair, hvap, vkc, grav, denice, c_to_b
     use clm_varcon         , only : denh2o, tfrz, tlsai_crit, alpha_aero
     use clm_varcon         , only : c14ratio
+    use clm_varcon         , only : c_water, c_dry_biomass
     use perf_mod           , only : t_startf, t_stopf
     use QSatMod            , only : QSat
     use CLMFatesInterfaceMod, only : hlm_fates_interface_type
@@ -420,7 +421,9 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
     real(r8) :: uuc(bounds%begp:bounds%endp)       ! undercanopy windspeed
     real(r8) :: carea_stem                         !cross-sectional area of stem
 
-    ! biomass parameters
+    ! Biomass heat storage tuning parameters
+    ! These parameters can be used to account for differences
+    ! in vegetation shape.
     real(r8), parameter :: k_vert = 0.1            !vertical distribution of stem
     real(r8), parameter :: k_cyl_vol = 1.0         !departure from cylindrical volume
     real(r8), parameter :: k_cyl_area = 1.0        !departure from cylindrical area
@@ -740,12 +743,13 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
 
             ! calculate specify heat capacity of vegetation
             ! as weighted averaged of dry biomass and water
-            ! lma_dry has units of kg dry mass/m2 here (table 2 of bonan 2017) 
-            ! cdry_biomass = 1400 J/kg/K, cwater = 4188 J/kg/K
-            cp_leaf(p)  = leaf_biomass(p) * (1400._r8*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*4188._r8)
+            ! lma_dry has units of kg dry mass/m2 here
+            ! (Appendix B of Bonan et al., GMD, 2018) 
+
+            cp_leaf(p)  = leaf_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water)
 
             ! cp-stem will have units J/k/ground_area
-            cp_stem(p) = stem_biomass(p) * (1400._r8*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*4188._r8)
+            cp_stem(p) = stem_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water)
             ! adjust for departure from cylindrical stem model
             cp_stem(p) = k_cyl_vol * cp_stem(p)
 
@@ -1471,15 +1475,19 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
 
          ! Downward longwave radiation below the canopy
 
-         dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(c) + &
-              emv(p)*emg(c)*sb*(tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))*(1.-frac_rad_abs_by_stem(p))+ts_ini(p)**3*(ts_ini(p) + 4._r8*dt_stem(p))*frac_rad_abs_by_stem(p))
+         dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(c) &
+              + emv(p)*emg(c)*sb*(tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) &
+              *(1.-frac_rad_abs_by_stem(p))+ts_ini(p)**3*(ts_ini(p) &
+              + 4._r8*dt_stem(p))*frac_rad_abs_by_stem(p))
 
          ! Upward longwave radiation above the canopy
 
          ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(c) &
-              + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb*((1.-frac_rad_abs_by_stem(p)) &
-              *tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))+frac_rad_abs_by_stem(p)*ts_ini(p)**3*(ts_ini(p) &
-              + 4._r8*dt_stem(p))) + emg(c)*(1._r8-emv(p))*sb*lw_grnd)
+              + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb &
+              *((1.-frac_rad_abs_by_stem(p))*tlbef(p)**3 &
+              *(tlbef(p) + 4._r8*dt_veg(p)) &
+              +frac_rad_abs_by_stem(p)*ts_ini(p)**3*(ts_ini(p)+ 4._r8*dt_stem(p))) &
+              + emg(c)*(1._r8-emv(p))*sb*lw_grnd)
 
 
          ! Calculate the skin temperature as a weighted sum of all the ground and vegetated fraction
diff --git a/src/biogeophys/EnergyFluxType.F90 b/src/biogeophys/EnergyFluxType.F90
index 07829cc199..1f1aa64a9e 100644
--- a/src/biogeophys/EnergyFluxType.F90
+++ b/src/biogeophys/EnergyFluxType.F90
@@ -441,12 +441,12 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp)
          avgflag='A', long_name='sensible heat from veg', &
          ptr_patch=this%eflx_sh_veg_patch, set_lake=0._r8, c2l_scale_type='urbanf')
 
-    this%eflx_sh_stem_patch(begp:endp) = spval
-    call hist_addfld1d (fname='FSH_STEM', units='W/m^2',  &
-         avgflag='A', long_name='sensible heat from stem', &
-         ptr_patch=this%eflx_sh_stem_patch, c2l_scale_type='urbanf',default = 'inactive')
-
     if (use_biomass_heat_storage) then
+       this%eflx_sh_stem_patch(begp:endp) = spval
+       call hist_addfld1d (fname='FSH_STEM', units='W/m^2',  &
+            avgflag='A', long_name='sensible heat from stem', &
+            ptr_patch=this%eflx_sh_stem_patch, c2l_scale_type='urbanf',default = 'inactive')
+
        this%dhsdt_canopy_patch(begp:endp) = spval
        call hist_addfld1d (fname='DHSDT_CANOPY', units='W/m^2',  &
             avgflag='A', long_name='change in canopy heat storage', &
diff --git a/src/main/clm_varcon.F90 b/src/main/clm_varcon.F90
index a6044dcb3e..ff9f5e0d2c 100644
--- a/src/main/clm_varcon.F90
+++ b/src/main/clm_varcon.F90
@@ -81,6 +81,8 @@ module clm_varcon
   real(r8), public :: alpha_aero = 1.0_r8                           ! constant for aerodynamic parameter weighting
   real(r8), public :: tlsai_crit = 2.0_r8                           ! critical value of elai+esai for which aerodynamic parameters are maximum
   real(r8), public :: watmin = 0.01_r8                              ! minimum soil moisture (mm)
+  real(r8), public :: c_water = 4188_r8                             ! specific heat of water   [J/kg/K]
+  real(r8), public :: c_dry_biomass  = 1400_r8                      ! specific heat of dry biomass
 
   real(r8), public :: re = SHR_CONST_REARTH*0.001_r8                ! radius of earth (km)