diff --git a/modules/util/gfs_reader.py b/modules/util/gfs_reader.py
index 972ec182f5f414e757abed0b5baac4a5ce740ed4..0a2902e491720b2e36f6a30f529474d6f8fcada3 100644
--- a/modules/util/gfs_reader.py
+++ b/modules/util/gfs_reader.py
@@ -125,6 +125,20 @@ def get_horz_layer_s(xr_dataset, fld_names, press, lon_lo=None, lon_hi=None, lat
     return sub_fld
 
 
+def get_time_interpolated_layer(ds_0, ds_1, time_0, time_1, time, fld_name, press, lon_lo, lon_hi, lat_lo, lat_hi):
+    lyr_0 = get_horz_layer(ds_0, fld_name, press, lon_lo, lon_hi, lat_lo, lat_hi)
+    lyr_1 = get_horz_layer(ds_1, fld_name, press, lon_lo, lon_hi, lat_lo, lat_hi)
+
+    lyr = xr.concat([lyr_0, lyr_1], 'time')
+    lyr = lyr.assign_coords(time=[time_0, time_1])
+
+    intrp_lyr = lyr.interp(time=time, method='linear')
+    intrp_lyr = intrp_lyr.expand_dims('channel')
+    intrp_lyr = intrp_lyr.assign_coords(channel=[fld_name])
+
+    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)