initial commit of ij fix

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 and lats):
+            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)

scm/etc/scripts/plot_configs/all_vars_test.ini

@@ -1,15 +1,15 @@
-scm_datasets = output_gabls3_SCM_GFS_v15p2/output.nc, output_gabls3_noahmp_SCM_GFS_v15p2_noahmp/output.nc,
-scm_datasets_labels = GFSv15p2, GFSv15p2-noahmp
-plot_dir = plots_all_vars/
+scm_datasets = output_atomic_ERA5_Jan16T22Jan18T06_SCM_GFS_v16_no_nsst/output.nc, output_atomic_ERA5_Jan16T22Jan18T06_dephy_SCM_GFS_v16_no_nsst/output.nc,
+scm_datasets_labels = GFSv16, GFSv16_dephy
+plot_dir = plots_all_vars_no_nsst/
 obs_file = ../data/comparison_data/gabls3_scm_cabauw_obs_v33.nc
-obs_compare = True
+obs_compare = False
 plot_ind_datasets = True
 time_series_resample = False
 
 [time_slices]
   [[total]]
-    start = 2006, 7, 1, 13
-    end = 2006, 7, 2, 12
+    start = 2020, 1, 16, 22
+    end = 2020, 1, 18, 07
 
 [time_snapshots]