diff --git a/modules/util/gfs_reader.py b/modules/util/gfs_reader.py
index 0a2902e491720b2e36f6a30f529474d6f8fcada3..a54a6080725d1f95a903640076ef0eca8caf1eab 100644
--- a/modules/util/gfs_reader.py
+++ b/modules/util/gfs_reader.py
@@ -139,6 +139,24 @@ def get_time_interpolated_layer(ds_0, ds_1, time_0, time_1, time, fld_name, pres
     return intrp_lyr
 
 
+def get_time_interpolated_layer_s(ds_0, ds_1, time_0, time_1, time, fld_name_s, press, lon_lo, lon_hi, lat_lo, lat_hi):
+    lyr_s_0 = []
+    lyr_s_1 = []
+    for fld_name in fld_name_s:
+        lyr = get_horz_layer(ds_0, fld_name, press, lon_lo, lon_hi, lat_lo, lat_hi)
+        lyr_s_0.append(lyr)
+
+        lyr = get_horz_layer(ds_1, fld_name, press, lon_lo, lon_hi, lat_lo, lat_hi)
+        lyr_s_1.append(lyr)
+
+    lyr = xr.concat([lyr_s_0, lyr_s_1], 'time')
+    lyr = lyr.assign_coords(time=[time_0, time_1])
+
+    intrp_lyr = lyr.interp(time=time, method='linear')
+
+    return intrp_lyr
+
+
 def get_vert_profile(xr_dataset, fld_name, lons, lats):
     if lon360:  # convert -180/+180 to 0,360
         lons = np.where(lons < 0, lons + 360, lons)