diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
index 0b5bbbb3e..fde7052ba 100644
--- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
@@ -534,8 +534,10 @@ 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
+
                  do is = isnow,0
-                   tsnoxy(ix,is)  = tgxy(ix)
+                   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