diff --git a/modules/aeolus/aeolus_amv.py b/modules/aeolus/aeolus_amv.py
index e5ebe51846d6e5d9d9a0b2b465f58b7ab6aadc75..7710b3f62486308ac31ce46d3f7c54b360c4cb86 100644
--- a/modules/aeolus/aeolus_amv.py
+++ b/modules/aeolus/aeolus_amv.py
@@ -415,7 +415,10 @@ def analyze2(raob_to_amv_dct, raob_dct):
     print('Hgh Spd min/max/mean: {0:.2f}  {1:.2f}  {2:.2f}'.format(good_amvs[hgh, sidx].min(), good_amvs[hgh, sidx].max(), good_amvs[hgh, sidx].mean()))
     print('Hgh Press min/max/mean: {0:.2f}  {1:.2f}  {2:.2f}'.format(good_amvs[hgh, pidx].min(), good_amvs[hgh, pidx].max(), good_amvs[hgh, pidx].mean()))
 
+    # Comparison to Level of Best Fit (LBF) ------------------------------------------------------------------
     bin_size = 200.0
+    bin_ranges = get_press_bin_ranges(50, 1050, bin_size=bin_size)
+
     vld_bf = bfs[:, 3] == 0
     keep_idxs = vld_bf
 
@@ -433,7 +436,6 @@ def analyze2(raob_to_amv_dct, raob_dct):
     print('press bias/rms: {0:.2f}  {1:.2f} '.format(pd_mean, np.sqrt(pd_mean**2 + pd_std**2)))
     print('------------------------------------------')
 
-    bin_ranges = get_press_bin_ranges(50, 1050, bin_size=bin_size)
     bin_pres = bin_data_by(diff, amv_p, bin_ranges)
 
     amv_spd = good_amvs[keep_idxs, sidx]
@@ -484,6 +486,75 @@ def analyze2(raob_to_amv_dct, raob_dct):
               .format(int(x_values[i]), num_pres[i], np.average(np.abs(bin_pres[i])), np.average(bin_pres[i]),
                       np.average(np.abs(bin_spd[i])), np.average(bin_spd[i]), np.average(np.abs(bin_dir[i])), np.average(bin_dir[i])))
 
+    # Direct comparison to RAOB profile ---------------------------------------------------------------
+    vld = raob_match[:, 3] == 0
+    keep_idxs = vld
+
+    amv_p = good_amvs[keep_idxs, pidx]
+    bf_p = raob_match[keep_idxs, 2]
+    diff = amv_p - bf_p
+    mad = np.average(np.abs(diff))
+    bias = np.average(diff)
+    print('********************************************************')
+    print('Number of good best fits: ', bf_p.shape[0])
+    print('press, MAD: {0:.2f}'.format(mad))
+    print('press, bias: {0:.2f}'.format(bias))
+    pd_std = np.std(diff)
+    pd_mean = np.mean(diff)
+    print('press bias/rms: {0:.2f}  {1:.2f} '.format(pd_mean, np.sqrt(pd_mean**2 + pd_std**2)))
+    print('------------------------------------------')
+    bin_pres = bin_data_by(diff, amv_p, bin_ranges)
+
+
+    amv_spd = good_amvs[keep_idxs, sidx]
+    amv_dir = good_amvs[keep_idxs, didx]
+    bf_spd, bf_dir = raob_match[keep_idxs, 0], raob_match[keep_idxs, 1]
+
+    #diff = amv_spd * units('m/s') - bf_spd
+    #diff = diff.magnitude
+    diff = amv_spd - bf_spd
+    spd_mad = np.average(np.abs(diff))
+    spd_bias = np.average(diff)
+    print('spd, MAD: {0:.2f}'.format(spd_mad))
+    print('spd, bias: {0:.2f}'.format(spd_bias))
+    spd_mean = np.mean(diff)
+    spd_std = np.std(diff)
+    print('spd MAD/bias/rms: {0:.2f}  {1:.2f}  {2:.2f}'.format(np.average(np.abs(diff)), spd_mean, np.sqrt(spd_mean**2 + spd_std**2)))
+    print('-----------------')
+    bin_spd = bin_data_by(diff, amv_p, bin_ranges)
+
+    dir_diff = direction_difference(amv_dir, bf_dir)
+    print('dir, MAD: {0:.2f}'.format(np.average(np.abs(dir_diff))))
+    print('dir bias: {0:.2f}'.format(np.average(dir_diff)))
+    print('-------------------------------------')
+    bin_dir = bin_data_by(dir_diff, amv_p, bin_ranges)
+
+    amv_u, amv_v = uv_from_spd_dir(good_amvs[keep_idxs, sidx], good_amvs[keep_idxs, didx])
+    bf_u, bf_v = uv_from_spd_dir(raob_match[keep_idxs, 0], raob_match[keep_idxs, 1])
+    u_diffs = amv_u - bf_u
+    v_diffs = amv_v - bf_v
+
+    vd = np.sqrt(u_diffs**2 + v_diffs**2)
+    vd_mean = np.mean(vd)
+    vd_std = np.std(vd)
+    print('VD bias/rms: {0:.2f}  {1:.2f}'.format(vd_mean, np.sqrt(vd_mean**2 + vd_std**2)))
+    print('******************************************************')
+
+    x_values = []
+    num_pres = []
+    num_spd = []
+    num_dir = []
+    print('level   num cases   hgt MAD/bias    spd MAD/bias     dir MAD/bias')
+    print('-------------------------------------------------------------------')
+    for i in range(len(bin_ranges)):
+        x_values.append(np.average(bin_ranges[i]))
+        num_pres.append(bin_pres[i].shape[0])
+        num_spd.append(bin_spd[i].shape[0])
+        num_dir.append(bin_dir[i].shape[0])
+
+        print('{0:d}        {1:d}       {2:.2f}/{3:.2f}       {4:.2f}/{5:.2f}        {6:.2f}/{7:.2f}'
+              .format(int(x_values[i]), num_pres[i], np.average(np.abs(bin_pres[i])), np.average(bin_pres[i]),
+                      np.average(np.abs(bin_spd[i])), np.average(bin_spd[i]), np.average(np.abs(bin_dir[i])), np.average(bin_dir[i])))
 
 # imports the S4 NOAA output
 # filename: full path as a string, '/home/user/filename'