diff --git a/.github/workflows/ci_run_scm_rts.yml b/.github/workflows/ci_run_scm_rts.yml index 3176274f7..632a44fbb 100644 --- a/.github/workflows/ci_run_scm_rts.yml +++ b/.github/workflows/ci_run_scm_rts.yml @@ -194,7 +194,7 @@ jobs: - name: Download SCM RT baselines run: | cd ${dir_bl} - wget ftp://ftp.rap.ucar.edu:/pub/ccpp-scm/rt-baselines-${{matrix.build-type}}.zip + wget https://dtcenter.ucar.edu/ccpp/users/rt/rt-baselines-${{matrix.build-type}}.zip unzip rt-baselines-${{matrix.build-type}}.zip - name: Compare SCM RT output to baselines diff --git a/test/cmp_rt2bl.py b/test/cmp_rt2bl.py index 91dcb27d5..f29ed5218 100755 --- a/test/cmp_rt2bl.py +++ b/test/cmp_rt2bl.py @@ -14,22 +14,20 @@ # parser = argparse.ArgumentParser() -parser.add_argument('-drt', '--dir_rt', help='Directory containing SCM RT output') -parser.add_argument('-dbl', '--dir_bl', help='Directory containing SCM RT baselines') -parser.add_argument('-plt', '--do_plot', help='If true, create plots of SCM RT output', action='store_true') +parser.add_argument('-drt', '--dir_rt', help='Directory containing SCM RT output', required=True) +parser.add_argument('-dbl', '--dir_bl', help='Directory containing SCM RT baselines', required=True) # def parse_args(): args = parser.parse_args() dir_rt = args.dir_rt dir_bl = args.dir_bl - do_plot = args.do_plot - return (dir_rt, dir_bl, do_plot) + return (dir_rt, dir_bl) # def main(): # - (dir_rt, dir_bl, do_plot) = parse_args() + (dir_rt, dir_bl) = parse_args() # error_count = 0 @@ -48,7 +46,7 @@ def main(): # Create plots between RTs and baselines (only if differences exist) if (result != 0): - plot_files = plot_results(file_rt, file_bl, do_plot) + plot_files = plot_results(file_bl, file_rt) # Setup output directories for plots. result = os.system("mkdir -p scm_rt_out/"+run["case"]+"/"+run["suite"]) diff --git a/test/cmp_scmout.py b/test/cmp_scmout.py new file mode 100755 index 000000000..b2a7886b1 --- /dev/null +++ b/test/cmp_scmout.py @@ -0,0 +1,52 @@ +#!/usr/bin/env python + +############################################################################## +# +# This script compares SCM RT output to baselines. +# +############################################################################## +import os +import sys +from os.path import exists +import argparse +from plot_scm_out import plot_results + +# +parser = argparse.ArgumentParser() +parser.add_argument('-fbl', '--file_bl', help='File containing SCM RT baselines', required=True) +parser.add_argument('-frt', '--file_rt', help='File containing SCM RT output') + +# +def parse_args(): + args = parser.parse_args() + file_rt = args.file_rt + file_bl = args.file_bl + return (file_bl, file_rt) + +# +def main(): + # + (file_bl, file_rt) = parse_args() + + plot_files = plot_results(file_bl, file_rt) + + # Put plots in local directory (scm_plots/) + result = os.system("mkdir -p scm_plots/") + com = "mv" + for plot_file in plot_files: + com = com + " " + plot_file + # end for + result = os.system(com+" scm_plots/") + + # Housekeeping + if file_rt is not None: + result = os.system('tar -cvf scm_out_abs.tar scm_plots/*.png') + result = os.system('tar -cvf /scratch1/data_untrusted/Dustin.Swales/scm_out_abs.tar scm_plots/*.png') + else: + result = os.system('tar -cvf scm_out_diff.tar scm_plots/*.png') + result = os.system('tar -cvf /scratch1/data_untrusted/Dustin.Swales/scm_out_diff.tar scm_plots/*.png') + # end if + result = os.system('rm -rf scm_plots/') +# +if __name__ == '__main__': + main() diff --git a/test/plot_scm_out.py b/test/plot_scm_out.py index 7641e80f9..d8474071e 100755 --- a/test/plot_scm_out.py +++ b/test/plot_scm_out.py @@ -1,56 +1,73 @@ #!/usr/bin/env python -############################################################################## +############################################################################### # -# This script compares SCM RT output to baselines. +# This script compares SCM output from two simulations. +# Defined here as Baseline (BL) and Regression Test (RT). # -############################################################################## +############################################################################### from netCDF4 import Dataset import numpy as np import matplotlib.pyplot as plt # -def plot_results(file_BL, file_RT, plot_RT_only): - # List of SCM output fields to plot (This is everything) - vars2plot = ["time_inst","time_diag","time_swrad","time_lwrad","time_rad","pres", \ - "pres_i","sigma","sigma_i","pres_s","qv","T","u","v","ql","qi","qc", \ - "qv_force_tend","T_force_tend","u_force_tend","v_force_tend","w_ls", \ - "u_g","v_g","dT_dt_rad_forc","h_advec_thil","h_advec_qt", \ - "v_advec_thil","v_advec_qt","T_s","lhf","shf","tprcp_inst", \ - "tprcp_rate_inst","t2m","q2m","ustar","tsfc","tau_u","tau_v","upd_mf", \ - "dwn_mf","det_mf","sfc_up_lw_land","sfc_up_lw_ice","sfc_up_lw_water", \ - "sfc_up_sw_dir_nir","sfc_up_sw_dif_nir","sfc_up_sw_dir_vis", \ - "sfc_up_sw_dif_vis","sfc_dwn_sw_dir_nir","sfc_dwn_sw_dif_nir", \ - "sfc_dwn_sw_dir_vis","sfc_dwn_sw_dif_vis","mp_prcp_inst", \ - "dcnv_prcp_inst","scnv_prcp_inst","pwat","dT_dt_lwrad","dT_dt_swrad", \ - "dT_dt_pbl","dT_dt_deepconv","dT_dt_shalconv","dT_dt_micro", \ - "dT_dt_ogwd","dT_dt_cgwd","dT_dt_phys","dT_dt_nonphys","dq_dt_pbl", \ - "dq_dt_deepconv","dq_dt_shalconv","dq_dt_micro","dq_dt_phys", \ - "dq_dt_nonphys","doz_dt_pbl","doz_dt_prodloss","doz_dt_oz","doz_dt_T", \ - "doz_dt_ovhd","doz_dt_phys","doz_dt_nonphys","du_dt_pbl","du_dt_ogwd", \ - "du_dt_deepconv","du_dt_cgwd","du_dt_shalconv","du_dt_phys", \ - "du_dt_nonphys","dv_dt_pbl","dv_dt_ogwd","dv_dt_deepconv","dv_dt_cgwd", \ - "dv_dt_shalconv","dv_dt_phys","dv_dt_nonphys","sfc_dwn_sw","sfc_up_sw", \ - "sfc_net_sw","sfc_dwn_lw","gflux","u10m","v10m","hpbl","tprcp_accum", \ - "ice_accum","snow_accum","graupel_accum","conv_prcp_accum", \ - "tprcp_rate_accum","ice_rate_accum","snow_rate_accum", \ - "graupel_rate_accum","conv_prcp_rate_accum","max_cloud_fraction", \ - "toa_total_albedo","vert_int_lwp_mp","vert_int_iwp_mp", \ - "vert_int_lwp_cf","vert_int_iwp_cf"] - # Use subset of available SCM output for plots. - vars2plot = ["qv","T","u","v","ql","qi","qc","sfc_dwn_sw","sfc_up_sw", \ - "sfc_net_sw","sfc_dwn_lw", "u10m","v10m","hpbl","gflux", \ - "qv_force_tend","T_force_tend","u_force_tend","v_force_tend","w_ls", \ - "u_g","v_g","dT_dt_rad_forc","h_advec_thil","h_advec_qt", \ - "v_advec_thil","v_advec_qt","T_s","lhf","shf","tprcp_inst", \ - "tprcp_rate_inst","t2m","q2m","ustar","tsfc","tau_u","tau_v"] +def plot_results(file_bl, file_rt=None, vars2plt=None): +#def plot_results(file_bl, file_rt, plot_bl, plot_rt, plot_all, debug): + # List of SCM output fields to plot + vars2plot_ALL = \ + ["pres", "pres_i","sigma","sigma_i","pres_s","qv","T","u","v","ql", \ + "qi","qc","qv_force_tend","T_force_tend","u_force_tend", \ + "v_force_tend","w_ls","u_g","v_g","dT_dt_rad_forc","h_advec_thil", \ + "h_advec_qt", "v_advec_thil","v_advec_qt","T_s","lhf","shf", \ + "tprcp_inst","tprcp_rate_inst","t2m","q2m","ustar","tsfc","tau_u", \ + "tau_v","upd_mf","dwn_mf","det_mf","sfc_up_lw_land","sfc_up_lw_ice", \ + "sfc_up_lw_water","sfc_up_sw_dir_nir","sfc_up_sw_dif_nir", \ + "sfc_up_sw_dir_vis","sfc_up_sw_dif_vis","sfc_dwn_sw_dir_nir", \ + "sfc_dwn_sw_dif_nir","sfc_dwn_sw_dir_vis","sfc_dwn_sw_dif_vis", \ + "mp_prcp_inst","dcnv_prcp_inst","scnv_prcp_inst","pwat", \ + "dT_dt_lwrad","dT_dt_swrad","dT_dt_pbl","dT_dt_deepconv", \ + "dT_dt_shalconv","dT_dt_micro","dT_dt_ogwd","dT_dt_cgwd", \ + "dT_dt_phys","dT_dt_nonphys","dq_dt_pbl","dq_dt_deepconv", \ + "dq_dt_shalconv","dq_dt_micro","dq_dt_phys","dq_dt_nonphys", \ + "doz_dt_pbl","doz_dt_prodloss","doz_dt_oz","doz_dt_T","doz_dt_ovhd", \ + "doz_dt_phys","doz_dt_nonphys","du_dt_pbl","du_dt_ogwd","dv_dt_cgwd",\ + "du_dt_deepconv","du_dt_cgwd","du_dt_shalconv","du_dt_phys", \ + "du_dt_nonphys","dv_dt_pbl","dv_dt_ogwd","dv_dt_deepconv", \ + "dv_dt_shalconv","dv_dt_phys","dv_dt_nonphys","sfc_dwn_sw", \ + "sfc_up_sw","sfc_net_sw","sfc_dwn_lw","gflux","u10m","v10m","hpbl", \ + "tprcp_accum","ice_accum","snow_accum","graupel_accum", \ + "conv_prcp_accum","tprcp_rate_accum","ice_rate_accum", \ + "snow_rate_accum","graupel_rate_accum","conv_prcp_rate_accum", \ + "max_cloud_fraction","toa_total_albedo","vert_int_lwp_mp", \ + "vert_int_iwp_mp","vert_int_lwp_cf","vert_int_iwp_cf"] + + # Smaller subset of SCM outputs to plot. + vars2plot_SUB = \ + ["qv","T","u","v","ql","qi","qc","sfc_dwn_sw","sfc_up_sw", \ + "sfc_net_sw","sfc_dwn_lw", "u10m","v10m","hpbl","gflux", \ + "qv_force_tend","T_force_tend","u_force_tend","v_force_tend","w_ls", \ + "u_g","v_g","dT_dt_rad_forc","h_advec_thil","h_advec_qt", \ + "v_advec_thil","v_advec_qt","T_s","lhf","shf","tprcp_inst", \ + "tprcp_rate_inst","t2m","q2m","ustar","tsfc","tau_u","tau_v"] + vars2plot_DEBUG = \ + ["qv","u10m"] - # Open SCM datasets - SCM_BL = Dataset(file_BL) - SCM_RT = Dataset(file_RT) + # Which fields to plot? (default is subset of full fields) + if vars2plt is None: + vars2plot = vars2plot_SUB + else: + vars2plot = vars2plt + # end if - plot_diff = True - if plot_RT_only: - plot_diff = False + # Only plot differences if two files provided + plot_diff = False + if file_rt is not None: + plot_diff = True + # end if + + # Open SCM datasets + SCM_BL = Dataset(file_bl) + if file_rt is not None: + SCM_RT = Dataset(file_rt) # end if plot_files = [] @@ -62,7 +79,13 @@ def plot_results(file_BL, file_RT, plot_RT_only): timeD = SCM_BL[var].dimensions[0] timeD = timeD[0:len(timeD)-4] x1 = SCM_BL[timeD][:].squeeze()/3600. #seconds -> hours - x2 = SCM_RT[timeD][:].squeeze()/3600. #seconds - >hours + if file_rt is not None: + x2 = SCM_RT[timeD][:].squeeze()/3600. #seconds - >hours + # If temporal dimensions disagree, con't compute deltas from experiments, turn off difference plots. + if (len(x1) != len(x2)): + plot_diff = False + # end if + # end if # Is this a 2D (time, x) variable? (Really 1D since x=1 in SCM) is2D = False if (len(SCM_BL[var].dimensions)==2): @@ -72,90 +95,88 @@ def plot_results(file_BL, file_RT, plot_RT_only): if (len(SCM_BL[var].shape) != 3): if (is2D): y1 = SCM_BL[var][:,0].squeeze() - y2 = SCM_RT[var][:,0].squeeze() + if file_rt is not None: + y2 = SCM_RT[var][:,0].squeeze() + # end if else: y1 = SCM_BL[var][:] - y2 = SCM_RT[var][:] + if file_rt is not None: + y2 = SCM_RT[var][:] + # end if # endif # Make figure if (np.size(x1) > 1): fig = plt.figure(figsize=(13,10)) - # Baselines and SCM RTs on same plot - if plot_diff: - plt.subplot(2,1,1) - # end if + + # Baselines and RTs on same plot + if plot_diff: plt.subplot(2,1,1) plt.title(SCM_BL[var].description) - if plot_RT_only: - plt.plot(x1, y1, color='blue') - else: - plt.plot(x1, y1, color='blue') - plt.plot(x2, y2, color='black') - # end if + plt.plot(x1, y1, color='blue') + if plot_diff: plt.plot(x2, y2, color='black') plt.ylabel('('+SCM_BL[var].units+')') plt.xlabel('(hours)') - # Difference (Baseline-SCMRT) + # Difference (Baseline-MRT) if plot_diff: plt.subplot(2,1,2) plt.title("Difference (blue - black)") plt.plot(x1, y1 - y2, color='red') plt.plot(x1, np.zeros(len(y1)), color='grey',linestyle='dashed') - plt.ylabel('('+SCM_RT[var].units+')') + plt.ylabel('('+SCM_BL[var].units+')') plt.xlabel('(hours)') - # + # Save figure fileOUT = 'scm.' + var +'.png' plt.savefig(fileOUT) - # plot_files.append(fileOUT) # three-dimensional variables elif len(SCM_BL[var].shape) == 3: z1 = np.transpose(SCM_BL[var][:,:,0]).squeeze() - z2 = np.transpose(SCM_RT[var][:,:,0]).squeeze() + if file_rt is not None: + z2 = np.transpose(SCM_RT[var][:,:,0]).squeeze() + # end if # vertical axis y1 = SCM_BL["pres"][0,:].squeeze()*0.01 - y2 = SCM_RT["pres"][0,:].squeeze()*0.01 + if file_rt is not None: + y2 = SCM_RT["pres"][0,:].squeeze()*0.01 + # end if nlev = SCM_BL[var][:,:,0].shape[1] # Layer (nlev) quantities are the default, so check for case where we have an # interface (nlev+1) variable to plot. if (SCM_BL[var][:,:,0].shape[1] > len(y1)): y1 = SCM_BL["pres_i"][0,:].squeeze()*0.01 - y2 = SCM_RT["pres_i"][0,:].squeeze()*0.01 + if file_rt is not None: + y2 = SCM_RT["pres_i"][0,:].squeeze()*0.01 + # end if # endif - # Comppute differences and determine valid plot range(s). - dz = z1-z2 - # Finally, make figure. if (np.size(x1) > 1): fig = plt.figure(figsize=(13,10)) - if plot_RT_only: - plt.contourf(x2, y2, z2, 20, cmap='YlGnBu') - plt.ylim(1000,200) - plt.ylabel('(Pa)') - plt.xlabel('(hours)') - cbr = plt.colorbar() - cbr.set_label('('+SCM_RT[var].units+')') - else: - # Baselines - plt.subplot(3,1,1) - plt.title(SCM_BL[var].description, fontsize=12) - plt.contourf(x1, y1, z1, 20, cmap='YlGnBu') - plt.ylim(1000,200) - plt.ylabel('(Pa)') - cbr = plt.colorbar() - cbr.set_label('('+SCM_RT[var].units+')') + if file_rt is not None: plt.subplot(3,1,1) + plt.contourf(x1, y1, z1, 20, cmap='YlGnBu') + plt.ylim(1000,200) + plt.xlim(0,np.max(x1)) + plt.ylabel('(Pa)') + plt.xlabel('(hours)') + cbr = plt.colorbar() + cbr.set_label('('+SCM_BL[var].units+')') + if file_rt is not None: # SCM RTs plt.subplot(3,1,2) plt.contourf(x2, y2, z2, 20, cmap='YlGnBu') plt.ylim(1000,200) + plt.xlim(0,np.max(x1)) plt.ylabel('(Pa)') plt.xlabel('(hours)') cbr = plt.colorbar() cbr.set_label('('+SCM_RT[var].units+')') - # Only plot differences if they exist. - if np.count_nonzero(dz) > 0: + # end if + # Only plot differences if requested, and only if they are non-zero. + if plot_diff: + dz = z1-z2 + if (np.count_nonzero(dz) > 0): plt.subplot(3,1,3) plt.title("Difference (top - middle)", fontsize=8) plt.contourf(x2, y2, dz, 20, cmap='bwr') @@ -165,9 +186,11 @@ def plot_results(file_BL, file_RT, plot_RT_only): cbr = plt.colorbar() cbr.set_label('('+SCM_RT[var].units+')') # end if (no differences exist) + # end if (plot differences) + # Save figure fileOUT = 'scm.' + var +'.png' + plt.show() plt.savefig(fileOUT) - # plot_files.append(fileOUT) # end if (Have enought pts to plot?) # end if (fields exist?)