From 81e618f5617c143fbb3ea3707c0357eae9de5c80 Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Wed, 19 Feb 2025 14:26:48 -0700 Subject: [PATCH] add Man's caseGen changes --- scm/etc/scripts/UFS_case_gen.py | 95 ++---------------------- scm/etc/scripts/UFS_case_gen_HYD.py | 66 ++-------------- scm/etc/scripts/scm_analysis.py | 3 + scm/etc/scripts/scm_plotting_routines.py | 16 ++-- scm/etc/scripts/scm_read_obs.py | 64 ++++++++++++++++ test/memo | 3 + test/plot_scm_out.py | 34 +++++---- 7 files changed, 111 insertions(+), 170 deletions(-) create mode 100644 test/memo diff --git a/scm/etc/scripts/UFS_case_gen.py b/scm/etc/scripts/UFS_case_gen.py index 0575ae79c..4074e9953 100755 --- a/scm/etc/scripts/UFS_case_gen.py +++ b/scm/etc/scripts/UFS_case_gen.py @@ -773,10 +773,7 @@ def get_IC_data_from_UFS_history(dir, i, j, lam, tile): icmr = nc_file['icmr'][0,::-1,i,j] pressfc = nc_file['pressfc'][0,i,j] tmp = nc_file['tmp'][0,::-1,i,j] - #mz - w_ls = nc_file['dzdt'][0,::-1,i,j] - omga = nc_file['omga'][0,::-1,i,j] - + delz = np.asarray(delz) zh = np.zeros(nlevs) zh[0] = 0.5*delz[0] @@ -793,8 +790,6 @@ def get_IC_data_from_UFS_history(dir, i, j, lam, tile): "dz": dz, "ua": np.asarray(ugrd), "va": np.asarray(vgrd), - 'wa': np.asarray(w_ls), - 'wap': np.asarray(omga), "qv": np.asarray(spfh), "o3": np.asarray(o3mr), "ql": np.asarray(clwmr), @@ -2400,9 +2395,6 @@ def get_UFS_forcing_data(nlevs, state_IC, location, use_nearest, forcing_dir, gr u_lay = [] v_lay = [] time_hr = [] - #mz - w_lay = [] - omga_lay = [] # Get grid from UFS IC data (ic_grid_lon, ic_grid_lat) = get_initial_lon_lat_grid(grid_dir, tile, lam) @@ -2444,9 +2436,6 @@ def get_UFS_forcing_data(nlevs, state_IC, location, use_nearest, forcing_dir, gr qv_data = regridder(nc_file['spfh'][0,::-1,:,:]) u_data = regridder(nc_file['ugrd'][0,::-1,:,:]) v_data = regridder(nc_file['vgrd'][0,::-1,:,:]) - #mz - w_data = regridder(nc_file['dzdt'][0,::-1,:,:]) - omga_data= regridder(nc_file['omga'][0,::-1,:,:]) i_get = 0 j_get = 0 # Same grids for history file (data_grid) to IC file (ic_grid). @@ -2456,9 +2445,6 @@ def get_UFS_forcing_data(nlevs, state_IC, location, use_nearest, forcing_dir, gr qv_data = nc_file['spfh'][0,::-1,:,:] u_data = nc_file['ugrd'][0,::-1,:,:] v_data = nc_file['vgrd'][0,::-1,:,:] - #mz - w_data = nc_file['dzdt'][0,::-1,:,:] - omga_data= nc_file['omga'][0,::-1,:,:] i_get = i j_get = j # end if @@ -2469,29 +2455,10 @@ def get_UFS_forcing_data(nlevs, state_IC, location, use_nearest, forcing_dir, gr qv_data = nc_file['spfh'][0,::-1,:,:] u_data = nc_file['ugrd'][0,::-1,:,:] v_data = nc_file['vgrd'][0,::-1,:,:] - #mz - w_data = nc_file['dzdt'][0,::-1,:,:] - omga_data= nc_file['omga'][0,::-1,:,:] j_get = tile_jj i_get = tile_ii # end if (use-nearest) - #mz placeholder to smooth - - - def smooth(data, window_size=5): - return np.convolve(data.flatten(), np.ones(window_size)/window_size, mode='same').reshape(data.shape) - - - # Smooth the data - w_data_smooth = smooth(w_data, window_size=5) - omga_data_smooth = smooth(omga_data, window_size=5) - - w_data = w_data_smooth - omga_data = omga_data_smooth - - #mz - # Store surface pressure. ps.append(ps_data[j_get,i_get]) @@ -2521,9 +2488,6 @@ def smooth(data, window_size=5): qv_lay.append(qv_data[:,j_get,i_get]) u_lay.append(u_data[:,j_get,i_get]) v_lay.append(v_data[:,j_get,i_get]) - #mz - w_lay.append(w_data[:,j_get,i_get]) - omga_lay.append(omga_data[:,j_get,i_get]) time_hr.append(nc_file['time'][0]) # Close file @@ -2538,9 +2502,6 @@ def smooth(data, window_size=5): qv_lay = np.asarray(qv_lay) u_lay = np.asarray(u_lay) v_lay = np.asarray(v_lay) - #mz - w_lay = np.asarray(w_lay) - omga_lay = np.asarray(omga_lay) time_hr = np.asarray(time_hr) # Compute virtual temperature. @@ -2767,21 +2728,12 @@ def smooth(data, window_size=5): qv_rev = np.zeros([1,nlevs]) u_rev = np.zeros([1,nlevs]) v_rev = np.zeros([1,nlevs]) - #mz - w_rev = np.zeros([1,nlevs]) - omga_rev = np.zeros([1,nlevs]) - tv_layr = np.zeros([n_files+1,nlevs]) qv_layr = np.zeros([n_files+1,nlevs]) u_layr = np.zeros([n_files+1,nlevs]) v_layr = np.zeros([n_files+1,nlevs]) p_layr = np.zeros([n_files+1,nlevs]) p_levr = np.zeros([n_files+1,nlevs+1]) - #mz - w_layr = np.zeros([n_files+1,nlevs]) - omga_layr = np.zeros([n_files+1,nlevs]) - - # # First timestep... @@ -2816,7 +2768,7 @@ def smooth(data, window_size=5): v_init_rev = state_IC["va"][::-1] v_rev_new = fv3_remap.map1_ppm(nlevs, from_p, v_init_rev[np.newaxis, :], 0.0, \ nlevs, to_p, 0, 0, -1, kord_tm ) - + # Store p_layr[0,:] = p_lay[0,:] p_levr[0,:] = p_lev[0,:] @@ -2824,10 +2776,6 @@ def smooth(data, window_size=5): u_layr[0,:] = u_rev_new[0,::-1] tv_layr[0,:] = tv_rev_new[0,::-1] qv_layr[0,:] = qv_rev_new[0,::-1] - #mz: not in IC - w_layr[0,:] = 0.0 - omga_layr[0,:] = 0.0 - # Subsequent timestep(s). (exact-mode only) if exact_mode: @@ -2855,16 +2803,6 @@ def smooth(data, window_size=5): v_rev[0,:] = v_lay[t,::-1] v_rev_new = fv3_remap.map1_ppm(nlevs, from_p, v_rev, 0.0, nlevs, to_p, \ 0, 0, -1, kord_tm ) - #mz - # vertical wind @ time > 0 - w_rev[0,:] = w_lay[t,::-1] - w_rev_new = fv3_remap.map1_ppm(nlevs, from_p, w_rev, 0.0, nlevs, to_p, \ - 0, 0, -1, kord_tm ) - # Meridional wind @ time > 0 - omga_rev[0,:] = omga_lay[t,::-1] - omga_rev_new = fv3_remap.map1_ppm(nlevs, from_p, omga_rev, 0.0, nlevs, to_p, \ - 0, 0, -1, kord_tm ) - # Store p_layr[t+1,:] = p_lay[t+1,:] p_levr[t+1,:] = p_lev[t+1,:] @@ -2872,9 +2810,6 @@ def smooth(data, window_size=5): qv_layr[t+1,:] = qv_rev_new[0,::-1] u_layr[t+1,:] = u_rev_new[0,::-1] v_layr[t+1,:] = v_rev_new[0,::-1] - #mz - w_layr[t+1,:] = w_rev_new[0,::-1] - omga_layr[t+1,:] = omga_rev_new[0,::-1] # end for # @@ -2884,9 +2819,6 @@ def smooth(data, window_size=5): qv_layr[t+2,:] = qv_layr[t+1,:] u_layr[t+2,:] = u_layr[t+1,:] v_layr[t+2,:] = v_layr[t+1,:] - #mz - w_layr[t+2,:] = w_layr[t+1,:] - omga_layr[t+2,:] = omga_layr[t+1,:] # end if (exact-mode) # Temperature @@ -2896,8 +2828,6 @@ def smooth(data, window_size=5): stateInit = {"p_lay": p_lay[0,:], \ "t_lay": t_lay[0,:], \ "qv_lay": qv_lay[0,:], \ - "w_lay": w_lay[0,:], \ - "omga_lay": omga_lay[0,:], \ "u_lay": u_lay[0,:], \ "v_lay": v_lay[0,:]} @@ -3042,10 +2972,6 @@ def smooth(data, window_size=5): tot_advec_qv = np.zeros((nlevs,ntimes),dtype=float) tot_advec_u = np.zeros((nlevs,ntimes),dtype=float) tot_advec_v = np.zeros((nlevs,ntimes),dtype=float) - #mz - w_ls = np.zeros((nlevs,ntimes),dtype=float) - omga = np.zeros((nlevs,ntimes),dtype=float) - time[0] = 0.0 time[1] = sec_in_hr*time_hr[0] - time_setback #forcing period should extend from beginning of diagnostic period to right BEFORE the next one @@ -3061,11 +2987,6 @@ def smooth(data, window_size=5): tot_advec_u[:,1] = tot_advec_u[:,0] tot_advec_v[:,0] = dvdt_adv[0,:] tot_advec_v[:,1] = tot_advec_v[:,0] - #mz - w_ls[:,0] = w_lay[0,:] - w_ls[:,1] = w_ls[:,0] - omga[:,0] = omga_lay[0,:] - omga[:,1] = omga[:,0] for t in range(1,n_files): time[2*t] = sec_in_hr*time_hr[t-1] @@ -3082,12 +3003,6 @@ def smooth(data, window_size=5): tot_advec_u[:,2*t+1] = tot_advec_u[:,2*t] tot_advec_v[:,2*t] = dvdt_adv[t,:] tot_advec_v[:,2*t+1] = tot_advec_v[:,2*t] - #mz - w_ls[:,2*t] = w_lay[t,:] - w_ls[:,2*t+1] = w_ls[:,2*t] - omga[:,2*t] = omga_lay[t,:] - omga[:,2*t+1] = omga[:,2*t] - # end for elif (time_method == 'gradient'): #this produced wonky results in the SCM; avoid until investigated more print('Forcing can be interpolated in time since the forcing terms are assumed to follow a constant time-gradient.') @@ -3140,14 +3055,14 @@ def smooth(data, window_size=5): # end if # - #mz w_ls = np.zeros((nlevs,ntimes),dtype=float) - #mz omega = np.zeros((nlevs,ntimes),dtype=float) + w_ls = np.zeros((nlevs,ntimes),dtype=float) + omega = np.zeros((nlevs,ntimes),dtype=float) rad_heating = np.zeros((nlevs,ntimes),dtype=float) forcing = { "time": time, "wa": w_ls.swapaxes(0,1), - "wap": omga.swapaxes(0,1), + "wap": omega.swapaxes(0,1), "tnta_rad": rad_heating.swapaxes(0,1), "ps_forc": np.ones(ntimes)*ps[0], "pa_forc": pressure_forc.swapaxes(0,1), diff --git a/scm/etc/scripts/UFS_case_gen_HYD.py b/scm/etc/scripts/UFS_case_gen_HYD.py index 9ab0a347d..8788d016d 100755 --- a/scm/etc/scripts/UFS_case_gen_HYD.py +++ b/scm/etc/scripts/UFS_case_gen_HYD.py @@ -774,9 +774,7 @@ def get_IC_data_from_UFS_history(dir, i, j, lam, tile): icmr = nc_file['icmr'][0,::-1,i,j] pressfc = nc_file['pressfc'][0,i,j] tmp = nc_file['tmp'][0,::-1,i,j] - #mz - omga = nc_file['omga'][0,::-1,i,j] - + #mz: calculate delz in hydrostatic equilibrium delz = np.zeros(nlevs) rho = np.zeros(nlevs) @@ -802,7 +800,6 @@ def get_IC_data_from_UFS_history(dir, i, j, lam, tile): "dz": dz, "ua": np.asarray(ugrd), "va": np.asarray(vgrd), - 'wap': np.asarray(omga), "qv": np.asarray(spfh), "o3": np.asarray(o3mr), "ql": np.asarray(clwmr), @@ -2468,9 +2465,6 @@ def get_UFS_forcing_data(nlevs, state_IC, location, use_nearest, forcing_dir, gr qv_data = nc_file['spfh'][0,::-1,:,:] u_data = nc_file['ugrd'][0,::-1,:,:] v_data = nc_file['vgrd'][0,::-1,:,:] - #mz - w_data = nc_file['dzdt'][0,::-1,:,:] - omga_data= nc_file['omga'][0,::-1,:,:] i_get = i j_get = j # end if @@ -2481,21 +2475,10 @@ def get_UFS_forcing_data(nlevs, state_IC, location, use_nearest, forcing_dir, gr qv_data = nc_file['spfh'][0,::-1,:,:] u_data = nc_file['ugrd'][0,::-1,:,:] v_data = nc_file['vgrd'][0,::-1,:,:] - omga_data= nc_file['omga'][0,::-1,:,:] j_get = tile_jj i_get = tile_ii # end if (use-nearest) - - def smooth(data, window_size=5): - return np.convolve(data.flatten(), np.ones(window_size)/window_size, mode='same').reshape(data.shape) - - - # Smooth the data - omga_data_smooth = smooth(omga_data, window_size=5) - omga_data = omga_data_smooth - - # Store surface pressure. ps.append(ps_data[j_get,i_get]) @@ -2525,7 +2508,6 @@ def smooth(data, window_size=5): qv_lay.append(qv_data[:,j_get,i_get]) u_lay.append(u_data[:,j_get,i_get]) v_lay.append(v_data[:,j_get,i_get]) - omga_lay.append(omga_data[:,j_get,i_get]) time_hr.append(nc_file['time'][0]) # Close file @@ -2540,7 +2522,6 @@ def smooth(data, window_size=5): qv_lay = np.asarray(qv_lay) u_lay = np.asarray(u_lay) v_lay = np.asarray(v_lay) - omga_lay = np.asarray(omga_lay) time_hr = np.asarray(time_hr) # Compute virtual temperature. @@ -2767,19 +2748,12 @@ def smooth(data, window_size=5): qv_rev = np.zeros([1,nlevs]) u_rev = np.zeros([1,nlevs]) v_rev = np.zeros([1,nlevs]) - #mz - w_rev = np.zeros([1,nlevs]) - omga_rev = np.zeros([1,nlevs]) - tv_layr = np.zeros([n_files+1,nlevs]) qv_layr = np.zeros([n_files+1,nlevs]) u_layr = np.zeros([n_files+1,nlevs]) v_layr = np.zeros([n_files+1,nlevs]) p_layr = np.zeros([n_files+1,nlevs]) p_levr = np.zeros([n_files+1,nlevs+1]) - omga_layr = np.zeros([n_files+1,nlevs]) - - # # First timestep... @@ -2822,9 +2796,6 @@ def smooth(data, window_size=5): u_layr[0,:] = u_rev_new[0,::-1] tv_layr[0,:] = tv_rev_new[0,::-1] qv_layr[0,:] = qv_rev_new[0,::-1] - #mz: not in IC - omga_layr[0,:] = 0.0 - # Subsequent timestep(s). (exact-mode only) if exact_mode: @@ -2852,12 +2823,6 @@ def smooth(data, window_size=5): v_rev[0,:] = v_lay[t,::-1] v_rev_new = fv3_remap.map1_ppm(nlevs, from_p, v_rev, 0.0, nlevs, to_p, \ 0, 0, -1, kord_tm ) - #mz - # vertical pressure wind @ time > 0 - omga_rev[0,:] = omga_lay[t,::-1] - omga_rev_new = fv3_remap.map1_ppm(nlevs, from_p, omga_rev, 0.0, nlevs, to_p, \ - 0, 0, -1, kord_tm ) - # Store p_layr[t+1,:] = p_lay[t+1,:] p_levr[t+1,:] = p_lev[t+1,:] @@ -2865,7 +2830,6 @@ def smooth(data, window_size=5): qv_layr[t+1,:] = qv_rev_new[0,::-1] u_layr[t+1,:] = u_rev_new[0,::-1] v_layr[t+1,:] = v_rev_new[0,::-1] - omga_layr[t+1,:] = omga_rev_new[0,::-1] # end for # @@ -2875,8 +2839,6 @@ def smooth(data, window_size=5): qv_layr[t+2,:] = qv_layr[t+1,:] u_layr[t+2,:] = u_layr[t+1,:] v_layr[t+2,:] = v_layr[t+1,:] - #mz - omga_layr[t+2,:] = omga_layr[t+1,:] # end if (exact-mode) # Temperature @@ -2886,7 +2848,6 @@ def smooth(data, window_size=5): stateInit = {"p_lay": p_lay[0,:], \ "t_lay": t_lay[0,:], \ "qv_lay": qv_lay[0,:], \ - "omga_lay": omga_lay[0,:], \ "u_lay": u_lay[0,:], \ "v_lay": v_lay[0,:]} @@ -3031,10 +2992,6 @@ def smooth(data, window_size=5): tot_advec_qv = np.zeros((nlevs,ntimes),dtype=float) tot_advec_u = np.zeros((nlevs,ntimes),dtype=float) tot_advec_v = np.zeros((nlevs,ntimes),dtype=float) - #mz - w_ls = np.zeros((nlevs,ntimes),dtype=float) - omga = np.zeros((nlevs,ntimes),dtype=float) - time[0] = 0.0 time[1] = sec_in_hr*time_hr[0] - time_setback #forcing period should extend from beginning of diagnostic period to right BEFORE the next one @@ -3050,9 +3007,6 @@ def smooth(data, window_size=5): tot_advec_u[:,1] = tot_advec_u[:,0] tot_advec_v[:,0] = dvdt_adv[0,:] tot_advec_v[:,1] = tot_advec_v[:,0] - #mz - omga[:,0] = omga_lay[0,:] - omga[:,1] = omga[:,0] for t in range(1,n_files): time[2*t] = sec_in_hr*time_hr[t-1] @@ -3069,10 +3023,6 @@ def smooth(data, window_size=5): tot_advec_u[:,2*t+1] = tot_advec_u[:,2*t] tot_advec_v[:,2*t] = dvdt_adv[t,:] tot_advec_v[:,2*t+1] = tot_advec_v[:,2*t] - #mz - omga[:,2*t] = omga_lay[t,:] - omga[:,2*t+1] = omga[:,2*t] - # end for elif (time_method == 'gradient'): #this produced wonky results in the SCM; avoid until investigated more print('Forcing can be interpolated in time since the forcing terms are assumed to follow a constant time-gradient.') @@ -3125,14 +3075,14 @@ def smooth(data, window_size=5): # end if # - #mz w_ls = np.zeros((nlevs,ntimes),dtype=float) - #mz omega = np.zeros((nlevs,ntimes),dtype=float) + w_ls = np.zeros((nlevs,ntimes),dtype=float) + omega = np.zeros((nlevs,ntimes),dtype=float) rad_heating = np.zeros((nlevs,ntimes),dtype=float) forcing = { "time": time, - #"wa": w_ls.swapaxes(0,1), - "wap": omga.swapaxes(0,1), + "wa": w_ls.swapaxes(0,1), + "wap": omega.swapaxes(0,1), "tnta_rad": rad_heating.swapaxes(0,1), "ps_forc": np.ones(ntimes)*ps[0], "pa_forc": pressure_forc.swapaxes(0,1), @@ -3229,12 +3179,10 @@ def write_SCM_case_file(state, surface, oro, forcing, init, case, date, forcing_ nc_file.adv_rv = forcing_off nc_file.adv_rt = forcing_off if (vertical_method == 2): - #mz - nc_file.forc_wa = forcing_off - nc_file.forc_wap = forcing_on + nc_file.forc_wa = forcing_on else: nc_file.forc_wa = forcing_off - nc_file.forc_wap = forcing_off + nc_file.forc_wap = forcing_off if (geos_wind_forcing): nc_file.forc_geo = forcing_on else: diff --git a/scm/etc/scripts/scm_analysis.py b/scm/etc/scripts/scm_analysis.py index 9e2b11ada..fa94a4e15 100755 --- a/scm/etc/scripts/scm_analysis.py +++ b/scm/etc/scripts/scm_analysis.py @@ -849,6 +849,9 @@ def replace_fill_with_nan(nc_ds, var_name, var, group, time_diag, pres_l, datase obs_dict = sro.read_gabls3_obs(obs_file, time_slices, date_inst) elif('MOSAiC' in case_name.strip()): obs_dict = sro.read_MOSAiC_obs(obs_file, time_slices, date_inst) + #mz + elif('fm1_sfs' in case_name.strip()): + obs_dict = sro.read_UFS_obs(obs_file, time_slices, date_inst) try: os.makedirs(plot_dir) diff --git a/scm/etc/scripts/scm_plotting_routines.py b/scm/etc/scripts/scm_plotting_routines.py index aa100bfd6..fb5f57094 100755 --- a/scm/etc/scripts/scm_plotting_routines.py +++ b/scm/etc/scripts/scm_plotting_routines.py @@ -61,7 +61,9 @@ def plot_profile_multi(z, values, labels, x_label, y_label, filename, obs_z=None return fig = plt.figure() - colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + #add orange + #colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + colors = ['#e41a1c','#ffa500','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] linestyles = ['-','--','-.',':'] markers = ['+','o','*','s','d','^','v','p','h'] linewidth_val = 2.0 @@ -302,7 +304,8 @@ def plot_time_series_multi(time, values, labels, x_label, y_label, filename, obs fig = plt.figure() lines = [] - colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + #colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + colors = ['#e41a1c','#ffa500','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] linestyles = ['-','--','-.',':'] markers = ['+','o','*','s','d','^','v','p','h'] linewidth_val = 2.0 @@ -507,7 +510,8 @@ def plot_profile_multi_ens(z, values, labels, x_label, y_label, filename, obs_z= fig = plt.figure() - colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + #colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + colors = ['#e41a1c','#ffa500','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] linestyles = ['-','--','-.',':'] markers = ['+','o','*','s','d','^','v','p','h'] hatches = ['/','\\','|','-'] @@ -684,7 +688,8 @@ def plot_time_series_multi_ens(time, values, labels, x_label, y_label, filename, values_ens_processed.append(np.percentile(values[i], [0, 10, 25, 50, 75, 90, 100], axis=0, interpolation='nearest')) lines = [] - colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + #colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + colors = ['#e41a1c','#ffa500','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] linestyles = ['-','--','-.',':'] markers = ['+','o','*','s','d','^','v','p','h'] hatches = ['/','\\','|','-'] @@ -814,7 +819,8 @@ def plot_time_series_multi_ens(time, values, labels, x_label, y_label, filename, def plot_scatter_multi(x, y, x_label, y_label, filename, x_lim=[], y_lim=[], color_index=None): plt.rc('text', usetex=latex_labels) - colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + #colors = ['#e41a1c','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] + colors = ['#e41a1c','#ffa500','#4daf4a','#377eb8','#984ea3','#ff7f00','#a65628','#f781bf','#ffff33'] if color_index is None: color_index = 0 else: diff --git a/scm/etc/scripts/scm_read_obs.py b/scm/etc/scripts/scm_read_obs.py index ba9dc712d..68c9edf2f 100644 --- a/scm/etc/scripts/scm_read_obs.py +++ b/scm/etc/scripts/scm_read_obs.py @@ -6,6 +6,70 @@ import math import forcing_file_common as ffc +#mz +def read_UFS_obs(obs_file, time_slices, date): + + obs_time_slice_indices = [] + + obs_fid = Dataset(obs_file, 'r') + obs_fid.set_auto_mask(False) + + obs_time = 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.size): + obs_date.append(obs_init_datetime + datetime.timedelta(seconds = obs_time[i])) + obs_date = np.array(obs_date) + + 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]) + print (start_date, end_date, start_date_index, end_date_index, obs_date[start_date_index], obs_date[end_date_index]) + + #find the index corresponding to the start of the simulations + obs_start_index = np.where(obs_date == date[0][0])[0] + obs_time = obs_time - obs_time[obs_start_index] + + #pressure stored in kPa (pressure is 0 for initial time) + obs_pres = obs_fid.variables['levs'][1:,:] + + #get average pressure levels + obs_pres_l = np.mean(obs_pres[:,:], (0)) + + #obs_theta = obs_fid.variables['potential_temp'][1:,:] + #print obs_theta[0,:] + #obs_T = (obs_pres/ffc.p0)**(ffc.R_dry/ffc.c_p)*obs_theta + + obs_T = obs_fid.variables['temp'][1:,:] + + #kg/kg + obs_q = obs_fid.variables['qv'][1:,:] + #print obs_q[0,:] + #obs_cld = obs_fid.variables['cloud_fraction'][1:,:] + #obs_qi= obs_fid.variables['qi'][:] + #obs_ql= obs_fid.variables['ql'][:] + obs_u = obs_fid.variables['u'][:] + obs_v = obs_fid.variables['v'][:] + + + #print obs_cld[0,:] + + return_dict = {'time': obs_time, 'date': obs_date, 'time_slice_indices': obs_time_slice_indices, 'pres_l': obs_pres_l, + 'T': obs_T, 'qv': obs_q, 'u': obs_u, 'v': obs_v} + + obs_fid.close() + + return return_dict + +#mz + def read_MOSAiC_obs(obs_file, time_slices, date): obs_time_slice_indices = [] diff --git a/test/memo b/test/memo new file mode 100644 index 000000000..ac5c244ce --- /dev/null +++ b/test/memo @@ -0,0 +1,3 @@ +./cmp_scmout.py -fbl /scratch2/BMC/ufs-phys/Man.Zhang/_case_gen/ccpp-scm_v7_hr_v4/scm/run/fm1_sfs_15c/output_fm1_sfs_NH_noW_15c_n008_SCM_GFS_v17_p8_ugwpv1/output.nc -frt /scratch2/BMC/ufs-phys/Man.Zhang/_case_gen/ccpp-scm_v7_hr_v4/scm/run/fm1_sfs_15c/output_fm1_sfs_HYD_noW_15c_n008_SCM_GFS_v17_p8_ugwpv1/output.nc + +./cmp_scmout.py -fbl /scratch2/BMC/ufs-phys/Man.Zhang/_case_gen/ccpp-scm_v7_hr_v4/scm/run/fm1_sfs_15c_12hly/output_fm1_sfs_NH_noW_15c_n000_SCM_GFS_v17_p8_ugwpv1/output.nc -frt /scratch2/BMC/ufs-phys/Man.Zhang/_case_gen/ccpp-scm_v7_hr_v4/scm/run/fm1_sfs_15c_12hly/output_fm1_sfs_HYD_noW_15c_n000_SCM_GFS_v17_p8_ugwpv1/output.nc diff --git a/test/plot_scm_out.py b/test/plot_scm_out.py index 67b35aa3b..062d66018 100755 --- a/test/plot_scm_out.py +++ b/test/plot_scm_out.py @@ -114,15 +114,15 @@ def plot_results(file_bl, file_rt=None, vars2plt=None): # Baselines and RTs on same plot if plot_diff: plt.subplot(2,1,1) plt.title(SCM_BL[var].description) - plt.plot(x1, y1, color='red') - if plot_diff: plt.plot(x2, y2, color='blue') + plt.plot(x1, y1, color='orange') + if plot_diff: plt.plot(x2, y2, color='red') plt.ylabel('('+SCM_BL[var].units+')') plt.xlabel('(hours)') # Difference (Baseline-MRT) if plot_diff: plt.subplot(2,1,2) - plt.title("Difference (red - blue)") + plt.title("Difference (Orange - red)") plt.plot(x1, y1 - y2, color='black') plt.plot(x1, np.zeros(len(y1)), color='grey',linestyle='dashed') plt.ylabel('('+SCM_BL[var].units+')') @@ -156,17 +156,18 @@ def plot_results(file_bl, file_rt=None, vars2plt=None): # Finally, make figure. if (np.size(x1) > 1): #fig = plt.figure(figsize=(13,10)) - ;z_min = min(np.min(z1), np.min(z2)) + + z_min = min(np.min(z1), np.min(z2)) z_max = max(np.max(z1), np.max(z2)) - z_min = -1*z_max + #z_min = -1*z_max fig = plt.figure(figsize=(5,10)) #mz Compute limits for color bar vmin1, vmax1 = np.min(z1), np.max(z1) if file_rt is not None: plt.subplot(3,1,1) - #plt.contourf(x1, y1, z1, 20, cmap='PuBu', vmin=vmin1, vmax=vmax1) - plt.contourf(x1, y1, z1, levels=np.linspace(z_min, z_max, 20), cmap='gist_ncar', vmin=z_min, vmax=z_max) - plt.ylim(1000,0.01) + plt.contourf(x1, y1, z1, 20, cmap='gist_ncar', vmin=vmin1, vmax=vmax1) + #plt.contourf(x1, y1, z1, levels=np.linspace(z_min, z_max, 20), cmap='gist_ncar', vmin=z_min, vmax=z_max) + plt.ylim(1000,100) plt.xlim(0,np.max(x1)) plt.ylabel('(Pa)') plt.xlabel('(hours)') @@ -174,22 +175,23 @@ def plot_results(file_bl, file_rt=None, vars2plt=None): cbr.set_label('('+SCM_BL[var].units+')') # Set custom y-ticks for the first subplot - y_ticks = [1000, 900, 850, 700, 500, 250, 100, 50, 0.1] - plt.yticks(y_ticks, fontsize=8) + #y_ticks = [1000, 900, 850, 700, 500, 250, 100, 50, 0.1] + y_ticks = [1000, 925, 850, 700, 500, 250, 100] + plt.yticks(y_ticks, fontsize=10) # Add only y-axis grid lines plt.grid(axis='y', linestyle='--', linewidth=0.5, color='gray') if file_rt is not None: # SCM RTs plt.subplot(3,1,2) - #plt.contourf(x2, y2, z2, 20, vmin=vmin1, vmax=vmax1, cmap='gist_ncar') - plt.contourf(x2, y2, z2, levels=np.linspace(z_min, z_max, 20), cmap='gist_ncar', vmin=z_min, vmax=z_max) - plt.ylim(1000,0.01) + plt.contourf(x2, y2, z2, 20, vmin=vmin1, vmax=vmax1, cmap='gist_ncar') + #plt.contourf(x2, y2, z2, levels=np.linspace(z_min, z_max, 20), cmap='gist_ncar', vmin=z_min, vmax=z_max) + plt.ylim(1000,100) plt.xlim(0,np.max(x1)) plt.ylabel('(Pa)') plt.xlabel('(hours)') cbr = plt.colorbar() cbr.set_label('('+SCM_RT[var].units+')') - plt.yticks(y_ticks, fontsize=8) + plt.yticks(y_ticks, fontsize=10) # Add only y-axis grid lines plt.grid(axis='y', linestyle='--', linewidth=0.5, color='gray') # end if @@ -204,12 +206,12 @@ def plot_results(file_bl, file_rt=None, vars2plt=None): c3 = plt.contourf(x2, y2, dz, 20, cmap='bwr', vmin=vmin, vmax=-vmin) # plt.title("Difference (top - middle)", fontsize=8) # plt.contourf(x2, y2, dz, 20, cmap='bwr') - plt.ylim(1000,0.01) + plt.ylim(1000,100) plt.ylabel('(Pa)') plt.xlabel('(hours)') cbr = plt.colorbar(c3) cbr.set_label('('+SCM_RT[var].units+')') - plt.yticks(y_ticks, fontsize=8) + plt.yticks(y_ticks, fontsize=10) # Add only y-axis grid lines plt.grid(axis='y', linestyle='--', linewidth=0.5, color='gray') # end if (no differences exist)