ij index fix for issue #531

scm/doc/TechGuide/chap_cases.rst

@@ -460,7 +460,8 @@ generating multiple cases at once.
    UFS_forcing_ensemble_generator.py [-h] -d DIR -n CASE_NAME
    (-lonl LON_1 LON_2 -latl LAT_1 LAT_2 -nens NENSMEMBERS |
    -lons [LON_LIST] -lats [LAT_LIST] |
-   -fxy [LON_LAT_FILE])
+   -fxy [LON_LAT_FILE] |
+   -tile [TILE_NUMBER] -is [I_INDEX_LIST] -js [J_INDEX_LIST])
    [-dt TIMESTEP] [-cres C_RES] [-sdf SUITE] [-sc] [-near] [-fm] [-vm] [-wn] [-geos]
 
 Mandatory arguments:
@@ -477,6 +478,8 @@ Mandatory arguments:
    -  ``--lon_list (-lons)`` AND ``--lat_list (-lats)``: longitude and latitude of cases
 
    -  ``--lonlat_file (fxy)``: file containing longitudes and latitudes
+
+   -  ``--tile (-tile) `` AND ``--i_list (-is)`` AND ``--j_list (-js)``: tile, i- and j-indices of cases
 
 Optional arguments:
 

scm/etc/scripts/UFS_case_gen.py

@@ -583,7 +583,8 @@ def find_loc_indices_UFS_IC(loc, dir, lam, tile, indices):
         logging.debug(message)
 
     #find index of closest point in the tile if indices are not specified
-    theta = 0.0
+    theta = missing_value
+    dist_min = missing_value
     if not indices:
         (tile_j, tile_i, point_lon, point_lat, dist_min, theta) = find_loc_indices(loc, dir, tile, lam)
         message = 'nearest IC file indices in tile {} - (i,j): ({},{})'.format(tile, tile_i, tile_j)
@@ -595,7 +596,7 @@ def find_loc_indices_UFS_IC(loc, dir, lam, tile, indices):
         tile_i = indices[0]
         tile_j = indices[1]
         #still need to grab the lon/lat if the tile and indices are supplied
-        (point_lon, point_lat) = find_lon_lat_of_indices(indices, grid_dir, tile, lam)
+        (point_lon, point_lat) = find_lon_lat_of_indices(indices, dir, tile, lam)
 
         message = 'This index has a central longitude/latitude of [{0},{1}]'.format(point_lon,point_lat)
         logging.debug(message)
@@ -643,16 +644,19 @@ def get_initial_lon_lat_grid(dir, tile, lam):
 
     return (longitude, latitude)    
 
-def compare_hist_and_IC_points(location, hist_lon, hist_lat, IC_lon, IC_lat, hist_dist, angle_to_hist_point, IC_dist, angle_to_IC_point):
+def compare_hist_and_IC_points(location, indices, tile, hist_lon, hist_lat, IC_lon, IC_lat, hist_dist, angle_to_hist_point, IC_dist, angle_to_IC_point):
     #determine distance and angle from IC point to hist point
     dist = haversine_distance(IC_lat,IC_lon,hist_lat,hist_lon)
     theta = math.atan2(math.sin(math.radians(hist_lon - IC_lon))*math.cos(math.radians(hist_lat)), math.cos(math.radians(IC_lat))*math.sin(math.radians(hist_lat)) - math.sin(math.radians(IC_lat))*math.cos(math.radians(hist_lat))*math.cos(math.radians(hist_lon - IC_lon)))
     theta_deg = math.fmod((math.degrees(theta) + 360.0), 360.0)
 
     message = 'Location summary'
-    message += '\n The location as entered is {} deg E {} deg N'.format(location[0],location[1])
+    if indices:
+        message += '\n The location as entered is i={} and j={} in tile {}'.format(indices[0],indices[1],tile)
+    else:
+        message += '\n The location as entered is {} deg E {} deg N'.format(location[0],location[1])
     message += '\n The closest point in the UFS history file is {} deg E {} deg N, or {} km away at a bearing of {} deg'.format(hist_lon, hist_lat, hist_dist, angle_to_hist_point)
-    message += '\n The closest point in the UFS IC files (native grid) is {} deg E {} deg N, or {} km away at a bearing of {} deg'.format(IC_lon, IC_lat, IC_dist, angle_to_IC_point)
+    message += '\n The closest point in the UFS IC files (native grid) is {} deg E {} deg N, or {} km away at a bearing of {} deg'.format(IC_lon, IC_lat, IC_dist if not indices else "0", angle_to_IC_point if not indices else "N/A")
     message += '\n Therefore, the history point (used to define the SCM grid column) is {} km away from the closest UFS native grid poing at a bearing of {} deg'.format(dist, theta_deg)
     logging.info(message)
 
@@ -3685,12 +3689,15 @@ def main():
      old_chgres, lam, save_comp, use_nearest, forcing_method, vertical_method, geos_wind_forcing, wind_nudge) = parse_arguments()
 
     #find indices corresponding to both UFS history files and initial condition (IC) files
-    (hist_i, hist_j, hist_lon, hist_lat, hist_dist_min, angle_to_hist_point, neighbors, dx, dy) = find_loc_indices_UFS_history(location, forcing_dir, lam)
-
     (IC_i, IC_j, tile, IC_lon, IC_lat, IC_dist_min, angle_to_IC_point) = find_loc_indices_UFS_IC(location, grid_dir, lam, tile, indices)
+
+    if indices:
+        location = (IC_lon, IC_lat)
+
+    (hist_i, hist_j, hist_lon, hist_lat, hist_dist_min, angle_to_hist_point, neighbors, dx, dy) = find_loc_indices_UFS_history(location, forcing_dir, lam)
 
     #compare the locations of the found history file point and UFS IC point
-    compare_hist_and_IC_points(location, hist_lon, hist_lat, IC_lon, IC_lat, hist_dist_min, angle_to_hist_point, IC_dist_min, angle_to_IC_point)
+    compare_hist_and_IC_points(location, indices, tile, hist_lon, hist_lat, IC_lon, IC_lat, hist_dist_min, angle_to_hist_point, IC_dist_min, angle_to_IC_point)
 
     #read in surface data for the intial conditions at the IC point
     surface_data = get_UFS_surface_data(in_dir, tile, IC_i, IC_j, old_chgres, lam)

scm/etc/scripts/UFS_forcing_ensemble_generator.py

@@ -22,6 +22,9 @@
 parser.add_argument('-lats', '--lat_list',       help='latitudes, separated by a space',       nargs='*', type=float, required=False)
 parser.add_argument('-fxy',  '--lonlat_file',    help='file containing longitudes and latitude',nargs=1,              required=False)
 parser.add_argument('-nens', '--nensmembers',    help='number of SCM UFS ensemble memebers to create',                 type=int,   required=False)
+parser.add_argument('-tile', '--tile',           help='FV3 native grid tile',                                          type=int,   required=False)
+parser.add_argument('-is',   '--i_list',         help='list of i-indices for the given tile of the FV3 native grid, separated by a space', nargs='*', type=int, required=False)
+parser.add_argument('-js',   '--j_list',         help='list of j-indices for the given tile of the FV3 native grid, separated by a space', nargs='*', type=int, required=False)
 parser.add_argument('-dt',   '--timestep',       help='SCM timestep, in seconds',                                      type=int,   default = 3600)
 parser.add_argument('-cres', '--C_RES',          help='UFS spatial resolution',                                        type=int,   default = 96)
 parser.add_argument('-sdf',  '--suite',          help='CCPP suite definition file to use for ensemble',                            default = 'SCM_GFS_v16')
@@ -70,6 +73,15 @@ def main():
         else:
             print("ERROR: Number of longitude/latitudes are inconsistent")
             exit()
+    elif (args.i_list and args.j_list):
+        if (len(args.i_list) == len(args.j_list)):
+            npts = len(args.i_list)
+        else:
+            print("ERROR: Number of i/j indices are inconsistent")
+            exit()
+        if not args.tile:
+            print("ERROR: Must provide a tile number for the given i/j indices")
+            exit()
         # end if
     # end if
 
@@ -121,11 +133,15 @@ def main():
         lat_list = lines[1]
         lats = eval(lat_list)
         npts = len(lons)
+    elif (args.i_list and args.j_list):
+        i_indices = np.asarray(args.i_list)
+        j_indices = np.asarray(args.j_list)
     else:
         print("ERROR: Must provide input points in one of the following formats:")
         print("  Using -nens [] -lonl [] -latl []  (e.g. -nens 20 -lonl 30 40 -latl 30 35)")
         print("  Using -lons [] -lats []           (e.g. -lons 203 204 205 -lats 30 30 30)")
         print("  Using -fxy                        (e.g. -fxy lonlat.txt w/ -lons [] -lats [])")
+        print("  Using -tile [] -is [] -js []      (e.g. -tile 5 -is 5 6 7 -js 40 40 40)")
         exit()
     # end if
 
@@ -158,7 +174,11 @@ def main():
         # Call UFS_case_gen.py
         case_name     = args.case_name +"_n" + str(pt).zfill(3)
         file_scminput = "../../data/processed_case_input/"+case_name+"_SCM_driver.nc"
-        com = "./UFS_case_gen.py -l " +str(lons[pt]) + " " + str(lats[pt]) + \
+        if('lons' in locals() and 'lats' in locals()):
+            loc_string = "-l " +str(lons[pt]) + " " + str(lats[pt])
+        else:
+            loc_string = "-ij " + str(i_indices[pt]) + " " + str(j_indices[pt]) + " -t " + str(args.tile)
+        com = "./UFS_case_gen.py " + loc_string + \
               " -i " + args.dir_ic + " -g " + args.dir_grid + " -f " + args.dir_forcing + " -n " + case_name + com_config
         print(com)
         os.system(com)