From d3969b66ae6f35bd5aaba79fb50a4468eb54e9cb Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Tue, 14 May 2024 16:03:53 -0400 Subject: [PATCH] save progress; add vertical forcing option --- scm/etc/scripts/UFS_case_gen.py | 119 +++++++------ .../plot_configs/all_vars_test_UFS.ini | 160 ++++++++++++++++++ scm/etc/scripts/scm_analysis.py | 4 +- scm/etc/scripts/scm_read_obs.py | 44 ++++- 4 files changed, 269 insertions(+), 58 deletions(-) create mode 100644 scm/etc/scripts/plot_configs/all_vars_test_UFS.ini diff --git a/scm/etc/scripts/UFS_case_gen.py b/scm/etc/scripts/UFS_case_gen.py index c421babea..79f46ca9d 100755 --- a/scm/etc/scripts/UFS_case_gen.py +++ b/scm/etc/scripts/UFS_case_gen.py @@ -61,20 +61,21 @@ parser = argparse.ArgumentParser() group1 = parser.add_mutually_exclusive_group(required=True) -group1.add_argument('-l', '--location', help='longitude and latitude in degress E and N, respectively, separated by a space', nargs=2, type=float) -group1.add_argument('-ij', '--index', help='i,j indices within the tile (if known - bypasses search for closest model point to lon/lat location)', nargs=2, type=int) -parser.add_argument('-d', '--date', help='date corresponding to initial conditions in YYYYMMDDHHMMSS format', required=False) -parser.add_argument('-i', '--in_dir', help='path to the directory containing the UFS initial conditions', required=True) -parser.add_argument('-g', '--grid_dir', help='path to the directory containing the UFS supergrid files (AKA "fix" directory)', required=True) -parser.add_argument('-f', '--forcing_dir', help='path to the directory containing the UFS history files', required=True) -parser.add_argument('-n', '--case_name', help='name of case', required=True) -parser.add_argument('-t', '--tile', help='tile of desired point (if known - bypasses tile search if present)', type=int, choices=range(1,8)) -parser.add_argument('-a', '--area', help='area of grid cell in m^2', type=float) -parser.add_argument('-oc', '--old_chgres', help='flag to denote that the initial conditions use an older data format (pre-chgres_cube)', action='store_true') -parser.add_argument('-lam', '--lam', help='flag to signal that the ICs and forcing is from a limited-area model run', action='store_true') -parser.add_argument('-sc', '--save_comp', help='flag to create UFS reference file for comparison', action='store_true') -parser.add_argument('-near','--use_nearest', help='flag to indicate using the nearest UFS history file gridpoint',action='store_true') -parser.add_argument('-fm', '--forcing_method', help='method used to calculate forcing (1=tendencies from dycore, 2=advective terms calculated, 3=total time tendency terms calculated)', type=int, choices=range(1,3), default=2) +group1.add_argument('-l', '--location', help='longitude and latitude in degress E and N, respectively, separated by a space', nargs=2, type=float) +group1.add_argument('-ij', '--index', help='i,j indices within the tile (if known - bypasses search for closest model point to lon/lat location)', nargs=2, type=int) +parser.add_argument('-d', '--date', help='date corresponding to initial conditions in YYYYMMDDHHMMSS format', required=False) +parser.add_argument('-i', '--in_dir', help='path to the directory containing the UFS initial conditions', required=True) +parser.add_argument('-g', '--grid_dir', help='path to the directory containing the UFS supergrid files (AKA "fix" directory)', required=True) +parser.add_argument('-f', '--forcing_dir', help='path to the directory containing the UFS history files', required=True) +parser.add_argument('-n', '--case_name', help='name of case', required=True) +parser.add_argument('-t', '--tile', help='tile of desired point (if known - bypasses tile search if present)', type=int, choices=range(1,8)) +parser.add_argument('-a', '--area', help='area of grid cell in m^2', type=float) +parser.add_argument('-oc', '--old_chgres', help='flag to denote that the initial conditions use an older data format (pre-chgres_cube)', action='store_true') +parser.add_argument('-lam', '--lam', help='flag to signal that the ICs and forcing is from a limited-area model run', action='store_true') +parser.add_argument('-sc', '--save_comp', help='flag to create UFS reference file for comparison', action='store_true') +parser.add_argument('-near','--use_nearest', help='flag to indicate using the nearest UFS history file gridpoint',action='store_true') +parser.add_argument('-fm', '--forcing_method', help='method used to calculate forcing (1=tendencies from dycore, 2=advective terms calculated, 3=total time tendency terms calculated)', type=int, choices=range(1,4), default=2) +parser.add_argument('-vert','--vertical_forcing',help='method used for vertical advection (1=calculated from history files, 2=vertical velocity provided)', type=int, choices=range(1,3), default=1) ######################################################################################## # @@ -96,6 +97,7 @@ def parse_arguments(): save_comp = args.save_comp use_nearest = args.use_nearest forcing_method = args.forcing_method + vertical_forcing = args.vertical_forcing #validate args if not os.path.exists(in_dir): @@ -135,7 +137,7 @@ def parse_arguments(): raise Exception(message) return (location, index, date_dict, in_dir, grid_dir, forcing_dir, tile, \ - area, case_name, old_chgres, lam, save_comp, use_nearest, forcing_method) + area, case_name, old_chgres, lam, save_comp, use_nearest, forcing_method, vertical_forcing) ######################################################################################## # @@ -1103,7 +1105,7 @@ def centered_diff_derivative(y, dx): return deriv -def get_UFS_forcing_data_advective_tendency(dir, i, j, tile, neighbors, dx, dy, nlevs, lam, save_comp_data): +def get_UFS_forcing_data_advective_tendency(dir, i, j, tile, neighbors, dx, dy, nlevs, lam, save_comp_data, vertical_forcing): # Determine UFS history file format (tiled/quilted) if lam: @@ -1270,42 +1272,43 @@ def get_UFS_forcing_data_advective_tendency(dir, i, j, tile, neighbors, dx, dy, # grid(True) # show() - v_advec_u[t,0] = 0.0 - v_advec_u[t,nlevs-1] = 0.0 - v_advec_v[t,0] = 0.0 - v_advec_v[t,nlevs-1] = 0.0 - v_advec_T[t,0] = 0.0 - v_advec_T[t,nlevs-1] = 0.0 - v_advec_qv[t,0] = 0.0 - v_advec_qv[t,nlevs-1] = 0.0 - for k in range(1, nlevs-1): - w_asc = max(smoothed_w[t,k], 0.0) - w_des = min(smoothed_w[t,k], 0.0) - gradient_T_asc = (temp[t,k,center,center] - temp[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) - gradient_T_des = (temp[t,k+1,center,center] - temp[t,k,center,center])/(pres[t,k+1]-pres[t,k]) - #gradient_T_asc = (smoothed_T[t,k,center,center] - smoothed_T[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) - #gradient_T_des = (smoothed_T[t,k+1,center,center] - smoothed_T[t,k,center,center])/(pres[t,k+1]-pres[t,k]) - gradient_qv_asc = (qv[t,k,center,center] - qv[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) - gradient_qv_des = (qv[t,k+1,center,center] - qv[t,k,center,center])/(pres[t,k+1]-pres[t,k]) - #gradient_qv_asc = (smoothed_qv[t,k,center,center] - smoothed_qv[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) - #gradient_qv_des = (smoothed_qv[t,k+1,center,center] - smoothed_qv[t,k,center,center])/(pres[t,k+1]-pres[t,k]) - #gradient_u_asc = (u_wind[t,k,center,center] - u_wind[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) - #gradient_u_des = (u_wind[t,k+1,center,center] - u_wind[t,k,center,center])/(pres[t,k+1]-pres[t,k]) - gradient_u_asc = (smoothed_u[t,k,center,center] - smoothed_u[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) - gradient_u_des = (smoothed_u[t,k+1,center,center] - smoothed_u[t,k,center,center])/(pres[t,k+1]-pres[t,k]) - #gradient_v_asc = (v_wind[t,k,center,center] - v_wind[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) - #gradient_v_des = (v_wind[t,k+1,center,center] - v_wind[t,k,center,center])/(pres[t,k+1]-pres[t,k]) - gradient_v_asc = (smoothed_v[t,k,center,center] - smoothed_v[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) - gradient_v_des = (smoothed_v[t,k+1,center,center] - smoothed_v[t,k,center,center])/(pres[t,k+1]-pres[t,k]) - - rho = pres[t,k]/(rdgas*temp[t,k,center,center]) - adiabatic_exp_comp_term = -(w_asc + w_des)*grav/cp - v_advec_u[t,k] = rho*grav*(w_asc*gradient_u_asc + w_des*gradient_u_des) - v_advec_v[t,k] = rho*grav*(w_asc*gradient_v_asc + w_des*gradient_v_des) - v_advec_T[t,k] = rho*grav*(w_asc*gradient_T_asc + w_des*gradient_T_des) + adiabatic_exp_comp_term - v_advec_qv[t,k] = rho*grav*(w_asc*gradient_qv_asc + w_des*gradient_qv_des) + if (vertical_forcing == 1): + v_advec_u[t,0] = 0.0 + v_advec_u[t,nlevs-1] = 0.0 + v_advec_v[t,0] = 0.0 + v_advec_v[t,nlevs-1] = 0.0 + v_advec_T[t,0] = 0.0 + v_advec_T[t,nlevs-1] = 0.0 + v_advec_qv[t,0] = 0.0 + v_advec_qv[t,nlevs-1] = 0.0 + for k in range(1, nlevs-1): + w_asc = max(smoothed_w[t,k], 0.0) + w_des = min(smoothed_w[t,k], 0.0) + gradient_T_asc = (temp[t,k,center,center] - temp[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) + gradient_T_des = (temp[t,k+1,center,center] - temp[t,k,center,center])/(pres[t,k+1]-pres[t,k]) + #gradient_T_asc = (smoothed_T[t,k,center,center] - smoothed_T[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) + #gradient_T_des = (smoothed_T[t,k+1,center,center] - smoothed_T[t,k,center,center])/(pres[t,k+1]-pres[t,k]) + gradient_qv_asc = (qv[t,k,center,center] - qv[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) + gradient_qv_des = (qv[t,k+1,center,center] - qv[t,k,center,center])/(pres[t,k+1]-pres[t,k]) + #gradient_qv_asc = (smoothed_qv[t,k,center,center] - smoothed_qv[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) + #gradient_qv_des = (smoothed_qv[t,k+1,center,center] - smoothed_qv[t,k,center,center])/(pres[t,k+1]-pres[t,k]) + #gradient_u_asc = (u_wind[t,k,center,center] - u_wind[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) + #gradient_u_des = (u_wind[t,k+1,center,center] - u_wind[t,k,center,center])/(pres[t,k+1]-pres[t,k]) + gradient_u_asc = (smoothed_u[t,k,center,center] - smoothed_u[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) + gradient_u_des = (smoothed_u[t,k+1,center,center] - smoothed_u[t,k,center,center])/(pres[t,k+1]-pres[t,k]) + #gradient_v_asc = (v_wind[t,k,center,center] - v_wind[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) + #gradient_v_des = (v_wind[t,k+1,center,center] - v_wind[t,k,center,center])/(pres[t,k+1]-pres[t,k]) + gradient_v_asc = (smoothed_v[t,k,center,center] - smoothed_v[t,k-1,center,center])/(pres[t,k]-pres[t,k-1]) + gradient_v_des = (smoothed_v[t,k+1,center,center] - smoothed_v[t,k,center,center])/(pres[t,k+1]-pres[t,k]) + + rho = pres[t,k]/(rdgas*temp[t,k,center,center]) + adiabatic_exp_comp_term = -(w_asc + w_des)*grav/cp + v_advec_u[t,k] = rho*grav*(w_asc*gradient_u_asc + w_des*gradient_u_des) + v_advec_v[t,k] = rho*grav*(w_asc*gradient_v_asc + w_des*gradient_v_des) + v_advec_T[t,k] = rho*grav*(w_asc*gradient_T_asc + w_des*gradient_T_des) + adiabatic_exp_comp_term + v_advec_qv[t,k] = rho*grav*(w_asc*gradient_qv_asc + w_des*gradient_qv_des) - if (False): + if (True): plot(h_advec_T[t,:]*86400.0,pres[t,:],"r") # #y_av = movingaverage(y, 10) plot(v_advec_T[t,:]*86400.0,pres[t,:],"b") @@ -1413,6 +1416,7 @@ def get_UFS_forcing_data_advective_tendency(dir, i, j, tile, neighbors, dx, dy, for var2d in vars2d: var2d["values"] = np.asarray(var2d["values"]) sec_in_hr = 3600. + time[0] = 0. forcing = { "time": time*sec_in_hr, "wa": smoothed_w, @@ -1440,7 +1444,7 @@ def get_UFS_forcing_data_advective_tendency(dir, i, j, tile, neighbors, dx, dy, ######################################################################################## # ######################################################################################## -def write_SCM_case_file(state, surface, oro, forcing, case, date): +def write_SCM_case_file(state, surface, oro, forcing, case, date, vertical_forcing): """Write all data to a netCDF file in the DEPHY-SCM format""" # Working types @@ -1507,7 +1511,10 @@ def write_SCM_case_file(state, surface, oro, forcing, case, date): nc_file.adv_qt = forcing_off nc_file.adv_rv = forcing_off nc_file.adv_rt = forcing_off - nc_file.forc_wa = forcing_off + if (vertical_forcing == 2): + nc_file.forc_wa = forcing_on + else: + nc_file.forc_wa = forcing_off nc_file.forc_wap = forcing_off nc_file.forc_geo = forcing_off nc_file.nudging_ua = forcing_off @@ -2050,7 +2057,7 @@ def main(): #read in arguments (location, indices, date, in_dir, grid_dir, forcing_dir, tile, area, case_name, - old_chgres, lam, save_comp, use_nearest, forcing_method) = parse_arguments() + old_chgres, lam, save_comp, use_nearest, forcing_method, vertical_forcing) = parse_arguments() (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) @@ -2083,7 +2090,7 @@ def main(): logging.critical(message) raise Exception(message) elif forcing_method == 2: - (forcing_data, comp_data) = get_UFS_forcing_data_advective_tendency(forcing_dir, hist_i, hist_j, tile, neighbors, dx, dy, state_data["nlevs"], lam, save_comp) + (forcing_data, comp_data) = get_UFS_forcing_data_advective_tendency(forcing_dir, hist_i, hist_j, tile, neighbors, dx, dy, state_data["nlevs"], lam, save_comp, vertical_forcing) elif forcing_method == 3: message = "Chose the total time tendency. This is removed for simplicity in this script." logging.critical(message) @@ -2094,7 +2101,7 @@ def main(): raise Exception(message) # Write SCM case file - fileOUT = write_SCM_case_file(state_data, surface_data, oro_data, forcing_data, case_name, date) + fileOUT = write_SCM_case_file(state_data, surface_data, oro_data, forcing_data, case_name, date, vertical_forcing) # write out to compare SCM output to (atmf for state variables and sfcf for physics # tendencies) diff --git a/scm/etc/scripts/plot_configs/all_vars_test_UFS.ini b/scm/etc/scripts/plot_configs/all_vars_test_UFS.ini new file mode 100644 index 000000000..25e2d3e3d --- /dev/null +++ b/scm/etc/scripts/plot_configs/all_vars_test_UFS.ini @@ -0,0 +1,160 @@ +scm_datasets = output_control_c96_n003_total_SCM_GFS_v16/output.nc, output_control_c96_n003_vert_SCM_GFS_v16/output.nc, +scm_datasets_labels = total, vert +plot_dir = plots_all_vars_UFS/ +obs_file = ../data/comparison_data/control_c192_n000_comp_data.nc +obs_compare = True +plot_ind_datasets = True +time_series_resample = False + +[time_slices] + [[all]] + start = 2021, 3, 22, 6 + end = 2021, 3, 23, 6 + +[time_snapshots] + +[plots] + [[profiles_mean]] + vars = qv, T, u, v, qc, ql, qi, w_ls, u_g, v_g, dT_dt_rad_forc, rad_cloud_fraction + vars_labels = 'specific humidity ($g$ $kg^{-1}$)', 'T (K)', 'u ($m$ $s^{-1}$)', 'v ($m$ $s^{-1}$)', 'total cloud water (res + sgs) ($g$ $kg^{-1}$)', 'res. liquid cloud water ($g$ $kg^{-1}$)', 'res. ice cloud water ($g$ $kg^{-1}$)', 'forcing vertical velocity ($m$ $s^{-1}$)', 'forcing geo. u ($m$ $s^{-1}$)', 'forcing geo. v ($m$ $s^{-1}$)','radiative forcing ($K$ $s^{-1}$)','cloud fraction' + vert_axis = pres_l + vert_axis_label = 'p (Pa)' + y_inverted = True + y_log = False + y_min_option = min #min, max, val (if val, add y_min = float value) + y_max_option = max #min, max, val (if val, add y_max = float value) + conversion_factor = 1000.0, 1.0, 1.0, 1.0, 1000.0, 1000.0, 1000.0, 1.0, 1.0, 1.0, 1.0, 1.0 + + [[profiles_mean_multi]] + [[[thil_adv]]] + vars = h_advec_thil, v_advec_thil + vars_labels = 'horizontal adv. $\theta_{il}$','vertical adv. $\theta_{il}$' + x_label = '$K$ $day^{-1}$' + conversion_factor = 86400.0 + [[[qt_adv]]] + vars = h_advec_qt, v_advec_qt + vars_labels = 'horizontal adv. $q_t$','vertical adv. $q_t$' + x_label = '$g$ $kg^{-1}$ $day^{-1}$' + conversion_factor = 8.64E7 + [[[conv_mass_flux]]] + vars = upd_mf, dwn_mf, det_mf + vars_labels = 'updraft mass flux','downdraft mass flux','detrainment mass flux' + x_label = '$kg$ $m^{-2}$ $s^{-1}$' + [[[T_tendencies]]] + vars = T_force_tend, dT_dt_pbl, dT_dt_deepconv, dT_dt_shalconv, dT_dt_micro, dT_dt_lwrad, dT_dt_swrad, dT_dt_ogwd, dT_dt_cgwd + vars_labels = 'force', 'PBL', 'Deep Con', 'Shal Con', 'MP', 'LW', 'SW', 'OGWD', 'CGWD' + x_label = '$K$ $day^{-1}$' + conversion_factor = 8.64E4 + [[[T_force_phys]]] + vars = T_force_tend, dT_dt_phys + vars_labels = 'forcing', 'physics' + x_label = '$K$ $day^{-1}$' + conversion_factor = 8.64E4 + [[[q_tendencies]]] + vars = qv_force_tend, dq_dt_pbl, dq_dt_deepconv, dq_dt_shalconv, dq_dt_micro + vars_labels = 'force', 'PBL', 'Deep Con', 'Shal Con', 'MP' + x_label = '$g$ $kg^{-1}$ $day^{-1}$' + conversion_factor = 8.64E7 + [[[q_force_phys]]] + vars = qv_force_tend, dq_dt_phys + vars_labels = 'forcing', 'physics' + x_label = '$g$ $kg^{-1}$ $day^{-1}$' + conversion_factor = 8.64E7 + [[[oz_tendencies]]] + vars = doz_dt_pbl, doz_dt_prodloss, doz_dt_oz, doz_dt_T, doz_dt_ovhd + vars_labels = 'PBL', 'prod/loss', 'oz', 'T', 'ovhd' + x_label = '$g$ $kg^{-1}$ $day^{-1}$' + conversion_factor = 8.64E7 + [[[u_tendencies]]] + vars = du_dt_pbl, du_dt_ogwd, du_dt_deepconv, du_dt_shalconv, du_dt_cgwd + vars_labels = 'PBL', 'OGWD', 'Deep Con', 'Shal Con', 'CGWD' + x_label = '$m$ $s^{-1}$ $day^{-1}$' + conversion_factor = 8.64E4 + [[[u_force_phys]]] + vars = u_force_tend, du_dt_phys + vars_labels = 'forcing', 'physics' + x_label = '$m$ $s^{-1}$ $day^{-1}$' + conversion_factor = 8.64E4 + [[[v_tendencies]]] + vars = dv_dt_pbl, dv_dt_ogwd, dv_dt_deepconv, dv_dt_shalconv, dv_dt_cgwd + vars_labels = 'PBL', 'OGWD', 'Deep Con', 'Shal Con', 'CGWD' + x_label = '$m$ $s^{-1}$ $day^{-1}$' + conversion_factor = 8.64E4 + [[[v_force_phys]]] + vars = v_force_tend, dv_dt_phys + vars_labels = 'forcing', 'physics' + x_label = '$m$ $s^{-1}$ $day^{-1}$' + conversion_factor = 8.64E4 + [[[water_path]]] + vars = rad_cloud_lwp, rad_cloud_iwp, rad_cloud_rwp, rad_cloud_swp + vars_labels = 'cloud liquid water path','cloud ice water path','rain water path','snow water path' + x_label = '$g$ $m^{-2}$' + [[[effective_radius]]] + vars = rad_eff_rad_ql, rad_eff_rad_qi, rad_eff_rad_qr, rad_eff_rad_qs + vars_labels = 'cloud liquid effective radius','cloud ice effective radius','rain effective radius','snow effective radius' + x_label = '$\mu m$' + + [[profiles_instant]] + + [[time_series]] + vars = pres_s, lhf, shf , T_s , tprcp_inst, tprcp_rate_inst, tau_u, tau_v, pwat, sfc_net_sw, tprcp_rate_accum + vars_labels = 'surface pressure','sfc latent heat flux ($W$ $m^{-2}$)','sfc sensible heat flux ($W$ $m^{-2}$)','surface forcing temperature (K)','instantaneous LWE total precipitation ($m$)','instantaneous LWE total precipitation rate ($mm$ $h^{-1}$)','sfc u stress ($Pa$)','sfc v stress ($Pa$)','precipitable water ($cm$)', 'sfc net SW (total) ($W$ $m^{-2}$)', 'accumulated total precipitation rate ($mm$ $h^{-1}$)' + conversion_factor = 1.0, 1.0, 1.0, 1.0, 1.0, 3.6E6, 1.0, 1.0, 0.1, 1.0, 3.6E6 + + [[time_series_multi]] + [[[soil_T]]] + vars = soil_T, soil_T, soil_T, soil_T + levels = 1, 2, 3, 4 + vars_labels = '10 cm','40 cm','100 cm','200 cm' + y_label = '$K$' + [[[soil_moisture]]] + vars = soil_moisture, soil_moisture, soil_moisture, soil_moisture + levels = 1, 2, 3, 4 + vars_labels = '10 cm','40 cm','100 cm','200 cm' + y_label = '$m^3$ $m^{-3}$' + [[[soil_moisture_unfrozen]]] + vars = soil_moisture_unfrozen, soil_moisture_unfrozen, soil_moisture_unfrozen, soil_moisture_unfrozen + levels = 1, 2, 3, 4 + vars_labels = '10 cm','40 cm','100 cm','200 cm' + y_label = '$m^3$ $m^{-3}$' + [[[surface_prcp_inst]]] + vars = tprcp_inst, mp_prcp_inst, dcnv_prcp_inst, scnv_prcp_inst + vars_labels = 'total','microphysics','deep conv.','shal conv.' + y_label = '$m$' + [[[surface_prcp_accum]]] + vars = tprcp_accum, ice_accum, snow_accum, graupel_accum, conv_prcp_accum + vars_labels = 'total','ice','snow','graupel','conv' + y_label = '$m$' + [[[surface_prcp_rate_accum]]] + vars = tprcp_rate_accum, ice_rate_accum, snow_rate_accum, graupel_rate_accum, conv_prcp_rate_accum + vars_labels = 'total','ice','snow','graupel','conv' + y_label = '$mm$ $h^{-1}$' + conversion_factor = 3.6E6 + obs_var = tprcp_rate_accum + [[[surface_sw_down]]] + vars = sfc_dwn_sw, sfc_dwn_sw_dir_nir, sfc_dwn_sw_dif_nir, sfc_dwn_sw_dir_vis, sfc_dwn_sw_dif_vis + vars_labels = 'sfc downwelling SW (total)', 'sfc downwelling SW (direct near-infrared)','sfc downwelling SW (diffuse near-infrared)','sfc downwelling SW (direct uv+vis)','sfc downwelling SW (diffuse uv+vis)' + y_label = '($W$ $m^{-2}$)' + [[[surface_sw_up]]] + vars = sfc_up_sw, sfc_up_sw_dir_nir, sfc_up_sw_dif_nir, sfc_up_sw_dir_vis, sfc_up_sw_dif_vis + vars_labels = 'sfc upwelling SW (total)','sfc upwelling SW (direct near-infrared)','sfc upwelling SW (diffuse near-infrared)', 'sfc upwelling SW (direct uv+vis)','sfc upwelling SW (diffuse uv+vis)' + y_label = '($W$ $m^{-2}$)' + [[[surface_lw]]] + vars = sfc_dwn_lw, sfc_up_lw_land, sfc_up_lw_water, sfc_up_lw_ice + vars_labels = 'sfc downwelling LW (total)','sfc upwelling LW (land)','sfc upwelling LW (water)','sfc upwelling LW (ice)' + y_label = '($W$ $m^{-2}$)' + + [[contours]] + vars = qv, T, u, v, qc, ql, qi + vars_labels = 'specific humidity ($g$ $kg^{-1}$)', 'T (K)', 'u ($m$ $s^{-1}$)', 'v ($m$ $s^{-1}$)', 'total cloud water (res + sgs) ($g$ $kg^{-1}$)', 'res. liquid cloud water ($g$ $kg^{-1}$)', 'res. ice cloud water ($g$ $kg^{-1}$)' + vert_axis = pres_l + vert_axis_label = 'p (Pa)' + y_inverted = True + y_log = False + y_min_option = val #min, max, val (if val, add y_min = float value) + y_min = 10000.0 + y_max_option = val #min, max, val (if val, add y_max = float value) + y_max = 100000.0 + x_ticks_num = 10 + y_ticks_num = 10 + conversion_factor = 1000.0, 1.0, 1.0, 1.0, 1000.0, 1000.0, 1000.0 diff --git a/scm/etc/scripts/scm_analysis.py b/scm/etc/scripts/scm_analysis.py index 34c85c529..18c369d21 100755 --- a/scm/etc/scripts/scm_analysis.py +++ b/scm/etc/scripts/scm_analysis.py @@ -856,6 +856,8 @@ def replace_fill_with_nan(nc_ds, var_name, var, group, time_diag, pres_l, datase obs_dict = sro.read_LASSO_obs(obs_file, time_slices, date_inst) elif('gabls3' in case_name.strip()): obs_dict = sro.read_gabls3_obs(obs_file, time_slices, date_inst) + elif('control' in case_name.strip()): + obs_dict = sro.read_UFS_comp(obs_file, time_slices, date_inst) try: os.makedirs(plot_dir) @@ -1124,7 +1126,7 @@ def replace_fill_with_nan(nc_ds, var_name, var, group, time_diag, pres_l, datase data_time_slice_series_rs = data_time_slice_series.resample(resample_string).mean() #print obs_data_time_slice.shape, obs_date_range.shape, data_time_slice_series_rs.shape - + print(obs_date_range, obs_data_time_slice) spr.plot_time_series_multi(obs_date_range, [data_time_slice_series_rs], [scm_datasets_labels[i]], 'date', label, ind_dir + '/time_series_' + plot_name + plot_ext, obs_time = obs_date_range, obs_values = obs_data_time_slice, line_type='color', color_index=i, conversion_factor=conversion_factor) elif(obs_delta_seconds < data_delta_seconds): print('The case where observations are more frequent than model output has not been implmented yet... ') diff --git a/scm/etc/scripts/scm_read_obs.py b/scm/etc/scripts/scm_read_obs.py index 59ac43c34..8d3440914 100644 --- a/scm/etc/scripts/scm_read_obs.py +++ b/scm/etc/scripts/scm_read_obs.py @@ -337,4 +337,46 @@ def read_gabls3_obs(obs_file, time_slices, date): 'sfc_dwn_sw': obs_sw_dn, 'sfc_up_sw': obs_sw_up, 'sfc_rad_net_land': obs_sfc_rad_net, 'gflux': -1*obs_gflux, 't2m':obs_t2m, 'q2m':obs_q2m, 'ustar':obs_ustar,'u10m':obs_u10m, 'v10m':obs_v10m, 'hpbl':obs_hpbl, 'tsfc':obs_tsk} - return return_dict \ No newline at end of file + return return_dict + +def read_UFS_comp(obs_file, time_slices, date): + + obs_fid = Dataset(obs_file, 'r') + obs_fid.set_auto_mask(False) + + obs_time_seconds_elapsed = obs_fid.variables['time'][:] + obs_datetime_string = obs_fid.getncattr('startDate') + + #get initial date from global file attribute + obs_init_datetime = datetime.datetime.strptime(obs_datetime_string, '%Y%m%d%H%M%S') + + obs_date = [] + for i in range(obs_time_seconds_elapsed.size): + obs_date.append(obs_init_datetime + datetime.timedelta(seconds = obs_time_seconds_elapsed[i])) + obs_date = np.array(obs_date) + + obs_time_slice_indices = [] + for time_slice in time_slices: + start_date = datetime.datetime(time_slices[time_slice]['start'][0], time_slices[time_slice]['start'][1],time_slices[time_slice]['start'][2], time_slices[time_slice]['start'][3], time_slices[time_slice]['start'][4]) + end_date = datetime.datetime(time_slices[time_slice]['end'][0], time_slices[time_slice]['end'][1],time_slices[time_slice]['end'][2], time_slices[time_slice]['end'][3], time_slices[time_slice]['end'][4]) + start_date_index = np.where(obs_date >= start_date)[0][0] + end_date_index = np.where(obs_date >= end_date)[0][0] + obs_time_slice_indices.append([start_date_index, end_date_index]) + + + obs_pres_l = obs_fid.variables['levs'][0,:] + obs_T = obs_fid.variables['temp'][:] + obs_qv = obs_fid.variables['qv'][:] + obs_u = obs_fid.variables['u'][:] + obs_v = obs_fid.variables['v'][:] + obs_shf = obs_fid.variables['shtfl'][:] + obs_lhf = obs_fid.variables['lhtfl'][:] + + obs_fid.close() + + return_dict = {'time': obs_time_seconds_elapsed, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, 'time_h': obs_time_seconds_elapsed/3600.0, + 'pres_l': obs_pres_l, 'T': obs_T, 'qv': obs_qv, 'u': obs_u, 'v': obs_v, 'shf': obs_shf, 'lhf': obs_lhf} + + return return_dict + + \ No newline at end of file