snapshot...

modules/amv/intercompare.py

@@ -923,13 +923,6 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
     imin = imin[0][0]
     imin = kk[0][imin]
 
-    # TDR - Feb 2020
-    if imin == fcst_num_levels - 1:
-        if verbose:
-            print('AMV location lat,lon,prs ({0},{1},{2}) failed to find min vector difference in layers around AMV'.format(amv_lat,amv_lon,amv_prs))
-        return bf_data
-
-
 # Use a parabolic fit to find the best-fit pressure.
 # p2 - Minimized model pressure at level imin (hPa)
 # v2 - Minimized vector difference at level imin (m/s)
@@ -940,11 +933,13 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
 
     p2 = fcst_prs[imin]
     v2 = VecDiff[imin]
+    SatwindBestFitPress = p2
+    SatwindBestFitU = fcst_uwind[imin]
+    SatwindBestFitV = fcst_vwind[imin]
 
 # assumes fcst data level 0 at surface and (fcst_num_levels-1) at model top
 # if bottom model level
     if imin == 0 or imin == fcst_num_levels - 1:
-        SatwindBestFitPress = p2
         if verbose:
             print('AMV location lat,lon,alt ({0},{1},{2}) min vector difference at top/bottom'.format(amv_lat,amv_lon,amv_prs))
     else:
@@ -956,9 +951,7 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
 # if top of allowed region
         if p3 < TopPress:
             SatwindBestFitPress = p2
-
-# check not collinear
-        elif v1 != v2 and v2 != v3:
+        elif v1 != v2 and v2 != v3: # check not collinear
             SatwindBestFitPress = p2 - (0.5 *
             ((((p2 - p1) * (p2 - p1) * (v2 - v3)) - ((p2 - p3) * (p2 - p3) * (v2 - v1))) /
             (((p2 - p1) * (v2 - v3)) - ((p2 - p3) * (v2 - v1)))))
@@ -967,14 +960,9 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
                     print('Best Fit not found between two pressure layers')
                     print('SatwindBestFitPress {0} p1 {1} p2 {2} p3 {3} imin {4}'.format(SatwindBestFitPress,p1,p2,p3,imin))
                 SatwindBestFitPress = p2
-        else:
-            SatwindBestFitPress = p2
 
 # Find best fit U and V by linear interpolation.
-    if p2 == SatwindBestFitPress:
-        SatwindBestFitU = fcst_uwind[imin]
-        SatwindBestFitV = fcst_vwind[imin]
-    else:
+    if p2 != SatwindBestFitPress:
         if p2 < SatwindBestFitPress:
             LevBelow = imin - 1
             LevAbove = imin
@@ -990,6 +978,7 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
 
 # Check to see if the best fit pressure is constrained.
     SatwindGoodConstraint = 1
+    flag = 0
     vdiff = np.sqrt((amv_uwind - SatwindBestFitU) ** 2 + (amv_vwind - SatwindBestFitV) ** 2)
     if vdiff > 4.0:
         SatwindGoodConstraint = 0
@@ -1003,12 +992,10 @@ def best_fit(amv_spd, amv_dir, amv_prs, amv_lat, amv_lon, fcst_spd, fcst_dir, fc
         SatwindGoodConstraint = 0
         flag = 1
 
-
     if SatwindGoodConstraint == 1:
         bfit_prs = SatwindBestFitPress
         bfit_u = SatwindBestFitU
         bfit_v = SatwindBestFitV
-        flag = 0
     else:
         bfit_prs = undef
         bfit_u = undef
@@ -1106,8 +1093,8 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
         a3 = fcst_alt[imin+1]
         v3 = vec_diff[imin+1]
 
-        # check not co-linear
-        if a1 != a2 and a2 != a3:
+        # check not co-linear and below alt_max
+        if a3 < alt_max and a1 != a2 and a2 != a3:
             sat_wind_best_fit_alt = minimize_quadratic(a1, a2, a3, v1, v2, v3)
             if sat_wind_best_fit_alt < a1 or sat_wind_best_fit_alt > a3:
                 sat_wind_best_fit_alt = a2
@@ -1128,6 +1115,7 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
 
     # check contraints
     good_constraint = 1
+    flag = 0
 
     # Check for overall good agreement: want vdiff less than 4 m/s
     vdiff = np.sqrt((amv_uwind - sat_wind_best_fit_u) ** 2 + (amv_vwind - sat_wind_best_fit_v) ** 2)
@@ -1144,7 +1132,7 @@ def best_fit_altitude(amv_spd, amv_dir, amv_alt, amv_lat, amv_lon, fcst_spd, fcs
         flag = 1
 
     if good_constraint == 1:
-        bf_tup = (sat_wind_best_fit_u, sat_wind_best_fit_v, sat_wind_best_fit_alt, 0)
+        bf_tup = (sat_wind_best_fit_u, sat_wind_best_fit_v, sat_wind_best_fit_alt, flag)
     else:
         bf_tup = (np.nan, np.nan, np.nan, flag)