diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90
index f53ab3928..6c810c622 100644
--- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.fv3.F90
@@ -593,8 +593,24 @@ subroutine GFS_phys_time_vary_init (
 
                  isnow = nint(snowxy(ix))+1 ! snowxy <=0.0, dzsno >= 0.0
 
+! using stc and tgxy to linearly interpolate the snow temp for each layer
+!Calculate the total thickness
+!                total_thickness = sum(dzsno)
+! Calculate the temperature difference
+!                temp_diff=tgxy(ix)-stc(ix,1)
+! Calculate the mid-points and interpolate temperatures for each layer
+!                 accumulated_thickness = 0.0
+!                do is = isnow, 0
+!                  accumulated_thickness = accumulated_thickness + dzsno(is)
+!                  mid_points(is) = accumulated_thickness - dzsno(is) / 2.0
+!                  layer_temp = tgxy(ix) + (mid_points(is) / total_thickness) * temp_diff
+!                  tsnoxy(ix,is) = layer_temp
+!                end do
+
                  do is = isnow,0
-                   tsnoxy(ix,is)  = tgxy(ix)
+!                  tsnoxy(ix,is)  = tgxy(ix)
+!                  tsnoxy(ix,is) =  tgxy(ix) + (( sum(dzsno(isnow:is)) -0.5*dzsno(is) )/snd)*(tgxy(ix)-stc(ix,1))
+                   tsnoxy(ix,is) =  tgxy(ix) + (( sum(dzsno(isnow:is)) -0.5*dzsno(is) )/snd)*(stc(ix,1)-tgxy(ix))
                    snliqxy(ix,is) = zero
                    snicexy(ix,is) = one * dzsno(is) * weasd(ix)/snd
                  enddo
diff --git a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90 b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90
index e519e472c..3451f807e 100644
--- a/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90
+++ b/physics/SFC_Models/Land/Noahmp/module_sf_noahmplsm.F90
@@ -1989,7 +1989,7 @@ subroutine energy (parameters,ice    ,vegtyp ,ist    ,nsnow  ,nsoil  , & !in
 
   real (kind=kind_phys), parameter                   :: mpe    = 1.e-6
   real (kind=kind_phys), parameter                   :: psiwlt = -150.  !metric potential for wilting point (m)
-  real (kind=kind_phys), parameter                   :: z0     = 0.002  ! bare-soil roughness length (m) (i.e., under the canopy)
+  real (kind=kind_phys), parameter                   :: z0     = 0.015  ! bare-soil roughness length (m) (i.e., under the canopy)
 
 ! ---------------------------------------------------------------------------------------------------
 ! initialize fluxes from veg. fraction
@@ -2626,10 +2626,10 @@ subroutine csnow (parameters,isnow   ,nsnow   ,nsoil   ,snice   ,snliq   ,dzsnso
 ! thermal conductivity of snow
 
   do iz = isnow+1, 0
-!     tksno(iz) = 3.2217e-6*bdsnoi(iz)**2.           ! stieglitz(yen,1965)
+     tksno(iz) = 3.2217e-6*bdsnoi(iz)**2.           ! stieglitz(yen,1965)
 !    tksno(iz) = 2e-2+2.5e-6*bdsnoi(iz)*bdsnoi(iz)   ! anderson, 1976
 !    tksno(iz) = 0.35                                ! constant
-    tksno(iz) = 2.576e-6*bdsnoi(iz)**2. + 0.074    ! verseghy (1991)
+!   tksno(iz) = 2.576e-6*bdsnoi(iz)**2. + 0.074    ! verseghy (1991)
 !    tksno(iz) = 2.22*(bdsnoi(iz)/1000.)**1.88      ! douvill(yen, 1981)
   enddo
 
@@ -5817,7 +5817,8 @@ subroutine thermalz0(parameters,    fveg,          z0m, z0mg,       zlvl,
 
       if (opt_trs == z0heqz0m) then
 
-        z0m_out = exp(fveg * log(z0m)      + (1.0 - fveg) * log(z0mg))
+!       z0m_out = exp(fveg * log(z0m)      + (1.0 - fveg) * log(z0mg))
+        z0m_out = fveg * z0m      + (1.0 - fveg) * z0mg
         z0h_out = z0m_out
 
       elseif (opt_trs == chen09) then
@@ -5834,7 +5835,7 @@ subroutine thermalz0(parameters,    fveg,          z0m, z0mg,       zlvl,
         endif
 
         z0h_out = exp( fveg        * log(z0m * exp(-czil*0.4*258.2*sqrt(ustarx*z0m))) + &
-                      (1.0 - fveg) * log(max(z0m/exp(kb_sigma_f0),1.0e-6)) )
+                      (1.0 - fveg) * log(max(z0mg/exp(kb_sigma_f0),1.0e-6)) )
 
       elseif (opt_trs == tessel) then