diff --git a/gsmphys/module_sf_noahmplsm.f90 b/gsmphys/module_sf_noahmplsm.f90
index 370f98bc..e486d8a6 100644
--- a/gsmphys/module_sf_noahmplsm.f90
+++ b/gsmphys/module_sf_noahmplsm.f90
@@ -172,6 +172,7 @@ module module_sf_noahmplsm
     real (kind=kind_phys) :: den                !tree density (no. of trunks per m2)
     real (kind=kind_phys) :: rc                 !tree crown radius (m)
     real (kind=kind_phys) :: mfsno              !snowmelt m parameter ()
+    real (kind=kind_phys) :: scffac             !snow cover factor (m)
     real (kind=kind_phys) :: saim(12)           !monthly stem area index, one-sided
     real (kind=kind_phys) :: laim(12)           !monthly leaf area index, one-sided
     real (kind=kind_phys) :: sla                !single-side leaf area per kg [m2/kg]
@@ -1657,7 +1658,7 @@ subroutine energy (parameters,ice    ,vegtyp ,ist    ,nsnow  ,nsoil  , & !in
      if(snowh.gt.0.)  then
          bdsno    = sneqv / snowh
          fmelt    = (bdsno/100.)**parameters%mfsno
-         fsno     = tanh( snowh /(2.5* z0 * fmelt))
+         fsno     = tanh( snowh /(parameters%scffac * fmelt))
 !        print*,'bdsno=',bdsno,sneqv,snowh,parameters%mfsno,fmelt,fsno
      endif
 
diff --git a/gsmphys/noahmp_tables.f90 b/gsmphys/noahmp_tables.f90
index f9421c7b..e5baec3e 100644
--- a/gsmphys/noahmp_tables.f90
+++ b/gsmphys/noahmp_tables.f90
@@ -88,11 +88,20 @@ module noahmp_tables
      &                               0.30, 0.30, 0.00, 0.00, 0.00, 0.00,     &
      &                               0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /
 
+    ! C. He 12/17/2020: optimized MFSNO values dependent on land type based on evaluation with SNOTEL SWE and MODIS SCF, surface albedo
     real (kind=kind_phys) :: mfsno_table(mvt)       !snowmelt curve parameter ()
-      data  ( mfsno_table(i),i=1,mvt) /  2.50, 2.50, 2.50, 2.50, 2.50, 2.50, &
-     &                               2.50, 2.50, 2.50, 2.50, 2.50, 2.50,     &
-     &                               2.50, 2.50, 2.50, 2.50, 2.50, 2.50,     &
-     &                               2.50, 2.50, 0.00, 0.00, 0.00, 0.00,     &
+      data  ( mfsno_table(i),i=1,mvt) /  1.00, 1.00, 1.00, 1.00, 1.00, 2.00, &
+     &                               2.00, 2.00, 2.00, 2.00, 3.00, 3.00,     &
+     &                               4.00, 4.00, 2.50, 3.00, 3.00, 3.50,     &
+     &                               3.50, 3.50, 0.00, 0.00, 0.00, 0.00,     &
+     &                               0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /
+
+    ! C. He 12/17/2020: optimized snow cover factor (m) in SCF formulation to replace original constant 2.5*z0,z0=0.002m, based on evaluation with SNOTEL SWE and MODIS SCF, surface albedo
+    real (kind=kind_phys) :: scffac_table(mvt)       !snow cover factor (m)
+      data  ( scffac_table(i),i=1,mvt) /  0.005, 0.005, 0.005, 0.005, 0.005, 0.008, &
+     &                               0.008, 0.010, 0.010, 0.010, 0.010, 0.007,      &
+     &                               0.021, 0.013, 0.015, 0.008, 0.015, 0.015,      &
+     &                               0.015, 0.015, 0.00, 0.00, 0.00, 0.00,          &
      &                               0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /
 
 !
diff --git a/gsmphys/sfc_noahmp_drv.f b/gsmphys/sfc_noahmp_drv.f
index 89f48c3f..f44db654 100644
--- a/gsmphys/sfc_noahmp_drv.f
+++ b/gsmphys/sfc_noahmp_drv.f
@@ -963,6 +963,7 @@ subroutine transfer_mp_parameters (vegtype,soiltype,slopetype,    &
         parameters%den    =    den_table(vegtype)       !tree density (no. of trunks per m2)
         parameters%rc     =     rc_table(vegtype)       !tree crown radius (m)
         parameters%mfsno  =  mfsno_table(vegtype)       !snowmelt m parameter ()
+        parameters%scffac =  scffac_table(vegtype)      !snow cover factor (m)
         parameters%saim   =   saim_table(vegtype,:)     !monthly stem area index, one-sided
         parameters%laim   =   laim_table(vegtype,:)     !monthly leaf area index, one-sided
         parameters%sla    =    sla_table(vegtype)       !single-side leaf area per kg [m2/kg]