diff --git a/modules/aeolus/aeolus_amv.py b/modules/aeolus/aeolus_amv.py
index 0e70081d6597396e3169a6ef554d5513561fec5b..e7f0f8cbb71fea1b5c78a43170ae068f1da49351 100644
--- a/modules/aeolus/aeolus_amv.py
+++ b/modules/aeolus/aeolus_amv.py
@@ -300,7 +300,6 @@ def run_best_fit_driver(amv_dir, source, raob_path, gfs_path, product_dir, produ
             for key_idx, key in enumerate(keys):
                 gfs_at_raob_dct[key] = (gfs_spd[key_idx], gfs_dir[key_idx], gfs_press)
 
-        #bf_dct = run_best_fit(m_d, raob_dct, gfs_file)
         bf_dct = run_best_fit(m_d, raob_dct, gfs_at_raob_dct)
 
         prd_dct = None
@@ -321,6 +320,7 @@ def run_best_fit_driver(amv_dir, source, raob_path, gfs_path, product_dir, produ
         keys = list(ab_dct.keys())
         for key in keys:
             tup = ab_dct.get(key)
+            print(tup[0].shape, tup[1].shape, tup[2].shape, tup[3].shape)
             amvs_list.append(tup[0])
             bfs_list.append(tup[1])
             rb_list.append(tup[2])
@@ -503,31 +503,12 @@ def run_best_fit(raob_to_amv_dct, raob_dct, gfs_at_raob_dct):
     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)
-    #     gfs_press = xr_dataset['pressure levels']
-    #     gfs_press = gfs_press.values
-    #     gfs_press = gfs_press[::-1]
-    #
-    #     lon_idx = LatLonTuple._fields.index('lon')
-    #     lat_idx = LatLonTuple._fields.index('lat')
-    #
-    #     uv_wind = get_vert_profile_s(xr_dataset, ['u-wind', 'v-wind'], locs[:, lon_idx], locs[:, lat_idx], method='nearest')
-    #     uv_wind = uv_wind.values
-    #     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]
-
     for key_idx, key in enumerate(keys):
         raob = raob_dct.get(key)
         raob_prs = raob[:, 0]
         raob_spd = raob[:, 3]
         raob_dir = raob[:, 2]
-        gfs = gfs_at_raob_dct.get(key)
+        gfs_spd, gfs_dir, gfs_press = gfs_at_raob_dct.get(key)
 
         amvs = raob_to_amv_dct.get(key)
         num_amvs = amvs.shape[1]
@@ -556,8 +537,7 @@ def run_best_fit(raob_to_amv_dct, raob_dct, gfs_at_raob_dct):
             raob_match_list.append(tup)
 
             if do_gfs_best_fit:
-                #bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, gfs_spd[key_idx], gfs_dir[key_idx], gfs_press)
-                bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, gfs[0], gfs[1], gfs[2])
+                bf = best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, gfs_spd, gfs_dir, gfs_press)
                 bf_gfs_list.append(bf)
 
         out_dct[key] = (np.transpose(amvs, axes=[1, 0]), np.array(bf_list), np.array(raob_match_list), np.array(bf_gfs_list))