diff --git a/cime_config/tests.py b/cime_config/tests.py
index b5a48bd2cfb6..5d4e3d247e23 100644
--- a/cime_config/tests.py
+++ b/cime_config/tests.py
@@ -50,7 +50,7 @@
             "ERS.f19_g16.I20TRGSWCNPRDCTCBC.elm-ctc_f19_g16_I20TRGSWCNPRDCTCBC",
             "SMS_Ly2_P1x1.1x1_smallvilleIA.I20TRGSWCNPCROP.elm-lulcc_sville",
             "ERS.r05_r05.ICNPRDCTCBC.elm-cbudget",
-            "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-snowveg_arctic",
+            "ERS_Ld150.ELM_USRDAT.I1850CNPRDCTCBC.elm-snowveg_arctic",
             "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_default_I1850CNPRDCTCBC",
             "ERS.ELM_USRDAT.I1850CNPRDCTCBC.elm-usrpft_codetest_I1850CNPRDCTCBC",
             "ERS.1x1_icycape.I1850GSWCNPRDCTCBC.elm-polygonal_tundra"
diff --git a/components/elm/src/biogeochem/VegStructUpdateMod.F90 b/components/elm/src/biogeochem/VegStructUpdateMod.F90
index 74e135f6f0ea..181fa23da08c 100644
--- a/components/elm/src/biogeochem/VegStructUpdateMod.F90
+++ b/components/elm/src/biogeochem/VegStructUpdateMod.F90
@@ -236,16 +236,19 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, &
          ! Wang and Zeng et al 2007.
          ! Taller vegetation (trees and shrubs) have been updated to use formulation similar to
          ! Sturm et al. 2005; Liston and Hiemstra, 2011; and Belke-Brea et al. 2020
-         if ( woody(ivt(p)) >= 1.0_r8 ) then
+         if ( woody(ivt(p)) == 1.0_r8 ) then
             ol = min( max(snow_depth(c)-hbot(p), 0._r8), htop(p)-hbot(p))
-            fb = 1._r8 - (ol / max(1.e-06_r8, bendresist(ivt(p)) * (htop(p)-hbot(p)))) ** vegshape(ivt(p))
+            fb = (ol / max(1.e-06_r8, (htop(p)-hbot(p)))) 
+         else if (woody(ivt(p)) == 2.0_r8) then
+            ol = min( max(snow_depth(c)-bendresist(ivt(p))*hbot(p), 0._r8), (bendresist(ivt(p))*(htop(p)-hbot(p)))**vegshape(ivt(p)))
+            fb = (ol / max(1.e-06_r8, bendresist(ivt(p))*(htop(p)-hbot(p))**vegshape(ivt(p))))
          else
-            fb = 1._r8 - max(min(snow_depth(c),0.2_r8),0._r8)/0.2_r8   ! 0.2m is assumed
+            fb = max(min(snow_depth(c),0.2_r8),0._r8)/0.2_r8   ! 0.2m is assumed
             !depth of snow required for complete burial of grasses
          endif
 
-         elai(p) = max(tlai(p)*fb, 0.0_r8)
-         esai(p) = max(tsai(p)*fb, 0.0_r8)
+         elai(p) = max(tlai(p)*(1._r8 - fb), 0.0_r8)
+         esai(p) = max(tsai(p)*(1._r8 - fb), 0.0_r8)
 
          ! Fraction of vegetation free of snow
          if ((elai(p) + esai(p)) > 0._r8) then
diff --git a/components/elm/src/main/pftvarcon.F90 b/components/elm/src/main/pftvarcon.F90
index 75cac61b47f8..822c40a2b760 100644
--- a/components/elm/src/main/pftvarcon.F90
+++ b/components/elm/src/main/pftvarcon.F90
@@ -315,7 +315,7 @@ module pftvarcon
   ! NGEE Arctic snow-vegetation interactions
   real(r8), allocatable :: bendresist(:)       ! vegetation resistance to bending under snow loading, 0 to 1 (e.g., Liston and Hiemstra 2011; Sturm et al. 2005)
   real(r8), allocatable :: vegshape(:)         ! shape parameter to modify shrub burial by snow (1 = parabolic, 2 = hemispheric)
-  real(r8), allocatable :: stocking(:)         ! stocking density for pft (stems / hectare)
+  real(r8), allocatable :: stocking(:)         ! stocking density for pft (stems / hectare assumed on param file)
   real(r8), allocatable :: taper(:)            ! ratio of height:radius_breast_height (woody vegetation allometry)
   logical               :: taper_defaults      ! set flag to use taper defaults if not on params file (necessary as import and set values are in different places)
 
@@ -1094,18 +1094,27 @@ subroutine pftconrd
     if (.not. readv ) bendresist(:) = 1._r8
     call ncd_io('vegshape', vegshape, 'read', ncid, readvar=readv, posNOTonfile=.true.)
     if (.not. readv ) vegshape(:) = 1._r8
-   ! check validity
+    ! check validity
     do i = 0, npft-1
       if (bendresist(i) .gt. 1.0_r8 .or. bendresist(i) .le. 0._r8) then
          call endrun(msg="Non-physical selection of bendresist parameter, set between 0 and 1"//errMsg(__FILE__, __LINE__))
       end if
-      if (vegshape(i) .gt. 2.0_r8 .or. vegshape(i) .le. 0._r8) then
-         call endrun(msg="Non-physical selection of vegshape parameter, set between 0 and 2"//errMsg(__FILE__, __LINE__))
+      if (vegshape(i) .gt. 2.0_r8 .or. vegshape(i) .lt. 1._r8) then
+         call endrun(msg="Non-physical selection of vegshape parameter, set as 1 or 2"//errMsg(__FILE__, __LINE__))
       end if
-   end do
+    end do
+
     call ncd_io('stocking', stocking, 'read', ncid, readvar=readv, posNOTonfile=.true.)
-    if (.not. readv ) stocking(:) = 0.1_r8 ! convert previous default of 1000 stems/ha to stems/m2 as had been done in VegStructUpdateMod.F90
-   taper_defaults = .false.
+    if (.not. readv ) then
+      stocking(:) = 0.1_r8 ! convert previous default of 1000 stems/ha to stems/m2 as had been done in VegStructUpdateMod.F90
+    else
+      stocking(:) = stocking(:)/10000._r8 ! conversion from stems/ha to stems/m2 done here now
+    endif
+    if (any(stocking(:) .ge. 1000._r8 )) then
+      call endrun(msg="Non-physical selection of stocking density parameter, implies >1000 plants/m2"//errMsg(__FILE__, __LINE__))
+    endif
+
+    taper_defaults = .false.
     call ncd_io('taper', taper, 'read', ncid, readvar=readv, posNOTonfile=.true.)
     if (.not. readv ) then
       taper(:) = 200._r8 ! pftnames not set to integers yet, so reassign further down.