diff --git a/modules/util/gfs_reader.py b/modules/util/gfs_reader.py
index e7cea8f3c99e7266429f01735e066b09077b19d7..eb73761d465bb509164232aa9a46dd2f0d81018a 100644
--- a/modules/util/gfs_reader.py
+++ b/modules/util/gfs_reader.py
@@ -187,7 +187,7 @@ def get_time_interpolated_layer_s(xr_dataset_s, time_s, time, fld_name_s, press,
     return intrp_lyr
 
 
-def get_vert_profile(xr_dataset, fld_name, lons, lats):
+def get_vert_profile(xr_dataset, fld_name, lons, lats, method='linear'):
     if lon360:  # convert -180/+180 to 0,360
         lons = np.where(lons < 0, lons + 360, lons)
 
@@ -198,12 +198,12 @@ def get_vert_profile(xr_dataset, fld_name, lons, lats):
     dim2 = xr.DataArray(lons, dims='k')
     dim1 = xr.DataArray(lats, dims='k')
 
-    intrp_fld = fld.interp(fakeDim1=dim1, fakeDim2=dim2, fakeDim0=plevs, method='linear')
+    intrp_fld = fld.interp(fakeDim1=dim1, fakeDim2=dim2, fakeDim0=plevs, method=method)
 
     return intrp_fld
 
 
-def get_vert_profile_s(xr_dataset, fld_name_s, lons, lats):
+def get_vert_profile_s(xr_dataset, fld_name_s, lons, lats, method='linear'):
     if lon360:  # convert -180,+180 to 0,360
         lons = np.where(lons < 0, lons + 360, lons)
 
@@ -219,7 +219,7 @@ def get_vert_profile_s(xr_dataset, fld_name_s, lons, lats):
     dim2 = xr.DataArray(lons, dims='k')
     dim1 = xr.DataArray(lats, dims='k')
 
-    intrp_fld = fld.interp(fakeDim1=dim1, fakeDim2=dim2, fakeDim0=plevs, method='linear')
+    intrp_fld = fld.interp(fakeDim1=dim1, fakeDim2=dim2, fakeDim0=plevs, method=method)
 
     return intrp_fld