diff --git a/modules/util/gfs_reader.py b/modules/util/gfs_reader.py
index b3cc75b25767ff60125276fe2884deaf99a12ef5..098a071d46848dc289129edfc279358598f34080 100644
--- a/modules/util/gfs_reader.py
+++ b/modules/util/gfs_reader.py
@@ -224,7 +224,7 @@ def get_vert_profile_s(xr_dataset, fld_name_s, lons, lats, method='linear'):
     return intrp_fld
 
 
-def get_point(xr_dataset, fld_name, lons, lats):
+def get_point(xr_dataset, fld_name, lons, lats, method='nearest'):
     if lon360:  # convert -180/+180 to 0,360
         lons = np.where(lons < 0, lons + 360, lons) # convert -180,180 to 0,360
 
@@ -237,7 +237,29 @@ def get_point(xr_dataset, fld_name, lons, lats):
     dim1 = xr.DataArray(lats, dims='k')
     dim2 = xr.DataArray(lons, dims='k')
 
-    intrp_fld = fld.interp(fakeDim1=dim1, fakeDim2=dim2, method='linear')
+    intrp_fld = fld.interp(fakeDim1=dim1, fakeDim2=dim2, method=method)
+
+    return intrp_fld
+
+
+def get_point_s(xr_dataset, fld_name_s, lons, lats, method='nearest'):
+    if lon360:  # convert -180/+180 to 0,360
+        lons = np.where(lons < 0, lons + 360, lons) # convert -180,180 to 0,360
+
+    lat_coords = np.linspace(-90, 90, xr_dataset.fakeDim1.size)
+    lon_coords = np.linspace(0, 359.5, xr_dataset.fakeDim2.size)
+
+    fld_s = []
+    for fld_name in fld_name_s:
+        fld = xr_dataset[fld_name]
+        fld = fld.assign_coords(fakeDim2=lon_coords, fakeDim1=lat_coords)
+        fld_s.append(fld)
+    fld = xr.concat(fld_s, 'fld_dim')
+
+    dim1 = xr.DataArray(lats, dims='k')
+    dim2 = xr.DataArray(lons, dims='k')
+
+    intrp_fld = fld.interp(fakeDim1=dim1, fakeDim2=dim2, method=method)
 
     return intrp_fld