diff --git a/modules/aeolus/aeolus_amv.py b/modules/aeolus/aeolus_amv.py
index b97220b8a2f82db98bfe7bb89ef722f2a71d4bdd..68099c36527eaa7e83bd048cc793d36bc533e749 100644
--- a/modules/aeolus/aeolus_amv.py
+++ b/modules/aeolus/aeolus_amv.py
@@ -10,6 +10,7 @@ from amv.intercompare import best_fit, bin_data_by, get_press_bin_ranges, spd_di
     direction_difference, run_best_fit_gfs2
 import math
 from metpy.units import units
+from util.gfs_reader import get_vert_profile_s
 
 
 amv_file_duration = 60  # minutes
@@ -371,20 +372,31 @@ def analyze2(raob_to_amv_dct, raob_dct, gfs_filename=None):
     keys = list(raob_to_amv_dct.keys())
 
     do_gfs_best_fit = False
+    gfs_press = None
+    gfs_spd = None
+    gfs_dir = None
     if gfs_filename is not None:
+        locs = np.array(keys)
         do_gfs_best_fit = True
         xr_dataset = xr.open_dataset(gfs_filename)
         bf_gfs_list = []
         gfs_press = xr_dataset['pressure levels']
         gfs_press = gfs_press.values
+        gfs_press = gfs_press[::-1]
+
+        uv_wind = get_vert_profile_s(xr_dataset, ['u-wind', 'v-wind'], locs[:, 1], locs[:, 0], method='nearest')
+        uv_wind = uv_wind.data
+        wspd, wdir = spd_dir_from_uv(uv_wind[0, :, :], uv_wind[1, :, :])
+        wspd = wspd.magnitude
+        wdir = wdir.magnitude
+        gfs_spd = wspd[:, ::-1]
+        gfs_dir = wdir[:, ::-1]
+
 
     amvs_list = []
     bf_list = []
     raob_match_list = []
-    for key in keys:
-        rlat = key[0]
-        rlon = key[1]
-
+    for key_idx, key in enumerate(keys):
         raob = raob_dct.get(key)
         nlevs = raob.shape[0]
         raob_prs = raob[:, 0]
@@ -402,9 +414,6 @@ def analyze2(raob_to_amv_dct, raob_dct, gfs_filename=None):
 
             bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, raob_spd, raob_dir, raob_prs)
             bf_list.append(bf)
-            if bf[3] == 0:
-                bspd, bdir = spd_dir_from_uv(bf[0], bf[1])
-                #print(amv_spd, bspd, amv_dir, bdir)
 
             pdiff = amv_prs - raob_prs
             lev_idx = np.argmin(np.abs(pdiff))
@@ -415,7 +424,7 @@ def analyze2(raob_to_amv_dct, raob_dct, gfs_filename=None):
             raob_match_list.append(tup)
 
             if do_gfs_best_fit:
-                bf = run_best_fit_gfs2(xr_dataset, gfs_press, amv_lat, amv_lon, amv_prs, amv_spd, amv_dir)
+                bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, gfs_spd[key_idx], gfs_dir[key_idx], gfs_press)
                 bf_gfs_list.append(bf)
 
     amvs = np.concatenate(amvs_list, axis=1)
diff --git a/modules/amv/intercompare.py b/modules/amv/intercompare.py
index bad32b7a560ce3f45e1ff4399251b9f93483881c..1fafc0ee170e850d285d834c5cd250368cf93520 100644
--- a/modules/amv/intercompare.py
+++ b/modules/amv/intercompare.py
@@ -612,25 +612,6 @@ def run_best_fit_gfs(amvs, gfs_filename, amv_lat_idx=amv_lat_idx, amv_lon_idx=am
 
     return bestfits
 
-def run_best_fit_gfs2(xr_dataset, gfs_press, amv_lat, amv_lon, amv_prs, amv_spd, amv_dir):
-    uv_wind = get_vert_profile_s(xr_dataset, ['u-wind', 'v-wind'], np.array([amv_lon]), np.array([amv_lat]))
-    uv_wind = uv_wind.data
-    wspd, wdir = spd_dir_from_uv(uv_wind[0,:,:], uv_wind[1,:,:])
-
-    rspd = wspd.magnitude
-    rdir = wdir.magnitude
-
-    rspd = rspd.flatten()
-    rdir = rdir.flatten()
-
-    gfs_press = gfs_press[::-1]
-    rspd = rspd[::-1]
-    rdir = rdir[::-1]
-
-    bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, rspd, rdir, gfs_press)
-    bf = (bf[0], bf[1], bf[2], bf[3])
-
-    return bf
 
 def raob_gfs_diff(raob_dct, gfs_filename):
     r_press = 900.0