diff --git a/scm/etc/scripts/UFS_case_gen.py b/scm/etc/scripts/UFS_case_gen.py index 4074e9953..e00055e12 100755 --- a/scm/etc/scripts/UFS_case_gen.py +++ b/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) diff --git a/scm/etc/scripts/UFS_forcing_ensemble_generator.py b/scm/etc/scripts/UFS_forcing_ensemble_generator.py index 1cacf9ec6..66e90ab85 100755 --- a/scm/etc/scripts/UFS_forcing_ensemble_generator.py +++ b/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) diff --git a/scm/etc/scripts/plot_configs/all_vars_test.ini b/scm/etc/scripts/plot_configs/all_vars_test.ini index a831d8945..e0da7716b 100644 --- a/scm/etc/scripts/plot_configs/all_vars_test.ini +++ b/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]