; Compute RMSE and % improvement of CTL and EXP runs from GFS for SHOUT-HRR, SHOUT-ENRR ; Runs through an array of GFS and EXP runs based on arrays to change below ; Outputs plots of RMSE and % improvement and an average percent improvement ; as a table ; Works for various variables of geopotential height, sea-level pressure, vector winds, relative humidity, and temperature ; Level dependent fields work as well by selecting the "level_selector" variable to an appropriate level ; Experiments are verified against the ECMWF analysis (not-reanalysis) ; Andrew Kren ; November 30, 2016 ; Updated December 7, 2016 to run over various levels and variables simultaneously ; Updated in February 6, 2016 to add in Student's t-test to results over each variable and over all cycles averaged totally ; Updated April 5, 2017 to plot the locations of significance in the RMSE plots, but not the relative RMSE ; Updated April 20, 2017 to make a histogram plot of all RMSE % changes for all variables in addition to table of results ; This is for better synopsis of the results ;;;;;;;;;;;;;;;;;;;;;;;;;;;; begin ; data grid N-S[90 to -90]; WE [0,359] nlat=181 nlon=360 ; forecast start and end dates with interval of 6 hours start_date = (/"20161005","20161005","20161005"/) utc_start = (/"06","12","18"/) interval = 6 end_date = (/"20161009","20161009","20161009"/) end_utc = (/"06","12","18"/) atmos_variable = (/"HGT","U","V","RH","TMP"/) ; settings to compute RMSE and %, SLP always computed so no need to put in here level_selector = (/200,300,500,700,850,925/) ; levels for RMSE and % improvement fcst_hr = ispan(0,96,6) ; forecast hours to read case_name = "drop_ctl_hrr" ; or drop_noNPP for output of filenames and netcdf file labels = (/"CTL","EXP"/) ; labels for correlation plots labels2 = (/"CTL","EXP"/) ; labels for table of statistics, don't change these labels ANL_DIR = "/scratch4/BMC/shout/ptmp/Andrew.Kren/ecmwf_analysis/" ; ECMWF analysis for verification ERA_DIR = "/scratch4/BMC/shout/Andrew.Kren/ECMWF/reanalysis/" ; used just for grid CTL_DIR = "/scratch4/BMC/shout/ptmp/Andrew.Kren/pre2016_hrr/" ; pre2016_hrr, pre2016_hrr_no_npp FL1_DIR = "/scratch4/BMC/shout/ptmp/Andrew.Kren/pre2016_hrr_drop_2/" ; pre2016_hrr_drop_2, pre2016_hrr_drop_fullres text2 = (/"METRIC","SE","CONUS","West US","East US","Global"/) ; labels for Table of statistics region_label = (/"SE"/) ; label for plots of Verification region, SE or AK lat_low = 25 ; coordinates used to plot the verification region for SE region lat_high = 35 lon_low = 277 lon_high = 287 ; lat_low = 55 ; coordinates used to plot the verification region for AK ; lat_high = 65 ; lon_low = 193 ; lon_high = 220 ; lat_low = -80 ; coordinates used to plot the northern hemisphere (20-80N) or southern hemisphere (20-80S) ; lat_high = -20 ; lon_low = 0 ; lon_high = 359 lat_low_2 = 25 ; coordinates used to plot the verification region for CONUS lat_high_2 = 50 lon_low_2 = 235 lon_high_2 = 295 lat_low_3 = 25 ; coordinates used to plot the verification region for WEST US lat_high_3 = 50 lon_low_3 = 235 lon_high_3 = 265 lat_low_4 = 25 ; coordinates used to plot the verification region for EAST US lat_high_4 = 50 lon_low_4 = 265 lon_high_4 = 295 siglvl = 0.05 ; level of statistical significance you want for your t-test, 0.05 would be 95% level, etc. ;;;;;; ; end of user's settings ;;;;;;; system("rm -f *.png") saved_pct_variable_ak = new((/dimsizes(atmos_variable),dimsizes(level_selector)/),"float") saved_pct_variable_conus = saved_pct_variable_ak saved_pct_variable_global = saved_pct_variable_ak saved_pct_variable_west = saved_pct_variable_ak saved_pct_variable_east = saved_pct_variable_ak do a=0,dimsizes(atmos_variable)-1 do b=0,dimsizes(level_selector)-1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; read in ERA reanalysis lat/lon information filename = ERA_DIR + "ecmwf_reanalysis.nc" w = addfile(filename,"r") lat = w->latitude lon = w->longitude delete(w) ; get weights for computations global_wgts = new((/nlat,nlon/),"float") wgts = latRegWgt(lat,"float",0) ; weights based on cosine of latitude, for all latitudes do x=0,nlon-1 global_wgts(:,x) = wgts end do delete(x) wgts_ak = global_wgts({lat_low:lat_high},:) ; weights for se region wgts_conus = global_wgts({lat_low_2:lat_high_2},:) ; weights for conus region wgts_west = global_wgts({lat_low_3:lat_high_3},:) ; weights for west conus region wgts_east = global_wgts({lat_low_4:lat_high_4},:) ; weights for east conus region ; arrays to save variables of RSME and % improvement saved_rmse_slp_ak_ctl = new((/dimsizes(start_date),dimsizes(fcst_hr)/),"float") saved_rmse_slp_conus_ctl = saved_rmse_slp_ak_ctl ; RMSE arrays saved_rmse_slp_global_ctl = saved_rmse_slp_ak_ctl saved_rmse_slp_west_ctl = saved_rmse_slp_ak_ctl saved_rmse_slp_east_ctl = saved_rmse_slp_ak_ctl saved_rmse_slp_west_drop = saved_rmse_slp_ak_ctl saved_rmse_slp_east_drop = saved_rmse_slp_ak_ctl saved_rmse_slp_ak_drop = saved_rmse_slp_ak_ctl saved_rmse_slp_conus_drop = saved_rmse_slp_ak_ctl saved_rmse_slp_global_drop = saved_rmse_slp_ak_ctl saved_rmse_hgt_ak_ctl = saved_rmse_slp_ak_ctl saved_rmse_hgt_conus_ctl = saved_rmse_slp_ak_ctl ; RMSE arrays saved_rmse_hgt_global_ctl = saved_rmse_slp_ak_ctl saved_rmse_hgt_ak_drop = saved_rmse_slp_ak_ctl saved_rmse_hgt_conus_drop = saved_rmse_slp_ak_ctl saved_rmse_hgt_global_drop = saved_rmse_slp_ak_ctl saved_rmse_hgt_west_ctl = saved_rmse_slp_ak_ctl saved_rmse_hgt_east_ctl = saved_rmse_slp_ak_ctl saved_rmse_hgt_west_drop = saved_rmse_slp_ak_ctl saved_rmse_hgt_east_drop = saved_rmse_slp_ak_ctl saved_pct_hgt_ak = saved_rmse_slp_ak_ctl saved_pct_hgt_conus = saved_rmse_slp_ak_ctl ; % improvement arrays saved_pct_hgt_global = saved_rmse_slp_ak_ctl saved_pct_slp_ak = saved_rmse_slp_ak_ctl saved_pct_slp_conus = saved_rmse_slp_ak_ctl ; % improvement arrays saved_pct_slp_global = saved_rmse_slp_ak_ctl saved_pct_slp_west = saved_rmse_slp_ak_ctl saved_pct_slp_east = saved_rmse_slp_ak_ctl saved_pct_hgt_west = saved_rmse_slp_ak_ctl saved_pct_hgt_east = saved_rmse_slp_ak_ctl print("Working on computations for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa") do i=0,dimsizes(start_date)-1 ; loop over all GFS cycles as defined in beginning of program print("Current GFS Initialization: " + start_date(i)+utc_start(i)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; year_st = stringtointeger(str_get_cols(start_date(i),0,3)) year_end = year_st + 1 dates = yyyymmddhh_time(year_st,year_end,interval,"string") ; string dates indice_start = ind(dates.eq.(start_date(i)+utc_start(i))) indice_end = ind(dates.eq.(end_date(i)+end_utc(i))) dates_plot = stringtointeger(dates(indice_start:indice_end)) dates_plot@units = "YYYYMMDDHH" valid_month = str_get_cols(tostring(dates_plot),4,5) valid_day = str_get_cols(tostring(dates_plot),6,7) valid_hour = str_get_cols(tostring(dates_plot),8,9) valid_year = str_get_cols(tostring(dates_plot),2,3) fc_files = dimsizes(dates_plot) ; number of forecast files vector_magnitude = 1.0 vector_length = 0.02 scale = 1.0e5 omega = 7.292e-05 colors = (/"Black","red", "green"/) ; colors for legend ; x and y locations for the legends, in normalized coordinates ; x and y locations for the legends, in normalized coordinates xdot = (/0.2,0.4, 0.6/) ydot = (/0.1,0.1,0.1/) xtxt = (/0.4,0.5,0.6/) ytxt = (/0.12,0.12,0.12/) xtxt2 = (/0.27,0.27/) ytxt2 = (/0.5,0.47/) err = NhlGetErrorObjectId() ; these 4 lines suppress grib warnings amongst others setvalues err "errLevel" : "Fatal" ; only report Fatal errors end setvalues ; --------------------------- looping over files ----------------------------------------- do j=0,fc_files-1 ; loop over gfs forecast times for this Initialization run print("Current analysis date for verification: " + dates_plot(j)) ; filename for this particular time period filename_ANL = ANL_DIR + "pgbanl.ecm." + dates_plot(j) + ".grb2" ; filename for this particular time period if (fcst_hr(j).lt.10) then filename_CTL = CTL_DIR + "pgbf0" + fcst_hr(j) + ".gfs." + start_date(i) + utc_start(i) + ".grb2" filename_FL1 = FL1_DIR + "pgbf0" + fcst_hr(j) + ".gfs." + start_date(i) + utc_start(i) + ".grb2" else filename_CTL = CTL_DIR + "pgbf" + fcst_hr(j) + ".gfs." + start_date(i) + utc_start(i) + ".grb2" filename_FL1 = FL1_DIR + "pgbf" + fcst_hr(j) + ".gfs." + start_date(i) + utc_start(i) + ".grb2" end if ; ecmwf analysis files w = addfile(filename_ANL,"r") mslp = w->PRMSL_GDS0_MSL mslp = mslp / 100.0 ; hPa ; check for variable based on atmos_variable setting if (atmos_variable(a).eq."HGT") then hgt_temp = w->HGT_GDS0_ISBL hgt = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."TMP") then hgt_temp = w->TMP_GDS0_ISBL hgt = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."RH") then hgt_temp = w->R_H_GDS0_ISBL hgt = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."U") then hgt_temp = w->U_GRD_GDS0_ISBL hgt = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."V") then hgt_temp = w->V_GRD_GDS0_ISBL hgt = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if delete(w) ; control case w = addfile(filename_CTL,"r") mslp_ctl = w->PRMSL_3_MSL mslp_ctl = mslp_ctl / 100.0 ; hPa ; check for variable based on atmos_variable setting if (atmos_variable(a).eq."HGT") then hgt_temp = w->HGT_3_ISBL hgt_ctl = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."TMP") then hgt_temp = w->TMP_3_ISBL hgt_ctl = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."RH") then hgt_temp = w->R_H_3_ISBL hgt_ctl = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."U") then hgt_temp = w->U_GRD_3_ISBL hgt_ctl = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."V") then hgt_temp = w->V_GRD_3_ISBL hgt_ctl = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if delete(w) ; experiment case w = addfile(filename_FL1,"r") mslp_drop = w->PRMSL_3_MSL mslp_drop = mslp_drop / 100.0 ; hPa ; check for variable based on atmos_variable setting if (atmos_variable(a).eq."HGT") then hgt_temp = w->HGT_3_ISBL hgt_drop = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."TMP") then hgt_temp = w->TMP_3_ISBL hgt_drop = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."RH") then hgt_temp = w->R_H_3_ISBL hgt_drop = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."U") then hgt_temp = w->U_GRD_3_ISBL hgt_drop = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if if (atmos_variable(a).eq."V") then hgt_temp = w->V_GRD_3_ISBL hgt_drop = hgt_temp({level_selector(b)},:,:) delete(hgt_temp) end if delete(w) ; for RMSE and % improvement calculations mslp_ak = mslp({lat_low:lat_high},{lon_low:lon_high}) mslp_conus = mslp({lat_low_2:lat_high_2},{lon_low_2:lon_high_2}) mslp_west = mslp({lat_low_3:lat_high_3},{lon_low_3:lon_high_3}) mslp_east = mslp({lat_low_4:lat_high_4},{lon_low_4:lon_high_4}) mslp_global = mslp hgt_ak = hgt({lat_low:lat_high},{lon_low:lon_high}) hgt_conus = hgt({lat_low_2:lat_high_2},{lon_low_2:lon_high_2}) hgt_west = hgt({lat_low_3:lat_high_3},{lon_low_3:lon_high_3}) hgt_east = hgt({lat_low_4:lat_high_4},{lon_low_4:lon_high_4}) hgt_global = hgt mslp_ak_ctl = mslp_ctl({lat_low:lat_high},{lon_low:lon_high}) mslp_conus_ctl = mslp_ctl({lat_low_2:lat_high_2},{lon_low_2:lon_high_2}) hgt_ak_ctl = hgt_ctl({lat_low:lat_high},{lon_low:lon_high}) hgt_conus_ctl = hgt_ctl({lat_low_2:lat_high_2},{lon_low_2:lon_high_2}) mslp_west_ctl = mslp_ctl({lat_low_3:lat_high_3},{lon_low_3:lon_high_3}) mslp_east_ctl = mslp_ctl({lat_low_4:lat_high_4},{lon_low_4:lon_high_4}) hgt_west_ctl = hgt_ctl({lat_low_3:lat_high_3},{lon_low_3:lon_high_3}) hgt_east_ctl = hgt_ctl({lat_low_4:lat_high_4},{lon_low_4:lon_high_4}) mslp_ak_drop = mslp_drop({lat_low:lat_high},{lon_low:lon_high}) mslp_conus_drop = mslp_drop({lat_low_2:lat_high_2},{lon_low_2:lon_high_2}) hgt_ak_drop = hgt_drop({lat_low:lat_high},{lon_low:lon_high}) hgt_conus_drop = hgt_drop({lat_low_2:lat_high_2},{lon_low_2:lon_high_2}) mslp_west_drop = mslp_drop({lat_low_3:lat_high_3},{lon_low_3:lon_high_3}) mslp_east_drop = mslp_drop({lat_low_4:lat_high_4},{lon_low_4:lon_high_4}) hgt_west_drop = hgt_drop({lat_low_3:lat_high_3},{lon_low_3:lon_high_3}) hgt_east_drop = hgt_drop({lat_low_4:lat_high_4},{lon_low_4:lon_high_4}) ; compute RMSE at this timestep of certain GFS Init. Run ;print("Computing RMSE and % improvement...") rmse_slp_ak_ctl = wgt_arearmse2(mslp_ak_ctl,mslp_ak,wgts_ak,1) ; 1 scalar value for this forecast hour rmse_slp_conus_ctl = wgt_arearmse2(mslp_conus_ctl,mslp_conus,wgts_conus,1) rmse_slp_west_ctl = wgt_arearmse2(mslp_west_ctl,mslp_west,wgts_west,1) ; 1 scalar value for this forecast hour rmse_slp_east_ctl = wgt_arearmse2(mslp_east_ctl,mslp_east,wgts_east,1) rmse_slp_global_ctl = wgt_arearmse2(mslp_ctl,mslp_global,global_wgts,1) rmse_hgt_ak_ctl = wgt_arearmse2(hgt_ak_ctl,hgt_ak,wgts_ak,1) ; 1 scalar value for this forecast hour rmse_hgt_conus_ctl = wgt_arearmse2(hgt_conus_ctl,hgt_conus,wgts_conus,1) rmse_hgt_global_ctl = wgt_arearmse2(hgt_ctl,hgt_global,global_wgts,1) rmse_hgt_west_ctl = wgt_arearmse2(hgt_west_ctl,hgt_west,wgts_west,1) ; 1 scalar value for this forecast hour rmse_hgt_east_ctl = wgt_arearmse2(hgt_east_ctl,hgt_east,wgts_east,1) rmse_slp_ak_drop = wgt_arearmse2(mslp_ak_drop,mslp_ak,wgts_ak,1) ; 1 scalar value for this forecast hour rmse_slp_conus_drop = wgt_arearmse2(mslp_conus_drop,mslp_conus,wgts_conus,1) rmse_slp_global_drop = wgt_arearmse2(mslp_drop,mslp_global,global_wgts,1) rmse_slp_west_drop = wgt_arearmse2(mslp_west_drop,mslp_west,wgts_west,1) ; 1 scalar value for this forecast hour rmse_slp_east_drop = wgt_arearmse2(mslp_east_drop,mslp_east,wgts_east,1) rmse_hgt_ak_drop = wgt_arearmse2(hgt_ak_drop,hgt_ak,wgts_ak,1) ; 1 scalar value for this forecast hour rmse_hgt_conus_drop = wgt_arearmse2(hgt_conus_drop,hgt_conus,wgts_conus,1) rmse_hgt_global_drop = wgt_arearmse2(hgt_drop,hgt_global,global_wgts,1) rmse_hgt_west_drop = wgt_arearmse2(hgt_west_drop,hgt_west,wgts_west,1) ; 1 scalar value for this forecast hour rmse_hgt_east_drop = wgt_arearmse2(hgt_east_drop,hgt_east,wgts_east,1) ; save variables for later usage saved_rmse_slp_ak_ctl(i,j) = rmse_slp_ak_ctl saved_rmse_slp_conus_ctl(i,j) = rmse_slp_conus_ctl saved_rmse_slp_global_ctl(i,j) = rmse_slp_global_ctl saved_rmse_slp_west_ctl(i,j) = rmse_slp_west_ctl saved_rmse_slp_east_ctl(i,j) = rmse_slp_east_ctl saved_rmse_hgt_ak_ctl(i,j) = rmse_hgt_ak_ctl saved_rmse_hgt_conus_ctl(i,j) = rmse_hgt_conus_ctl saved_rmse_hgt_global_ctl(i,j) = rmse_hgt_global_ctl saved_rmse_hgt_west_ctl(i,j) = rmse_hgt_west_ctl saved_rmse_hgt_east_ctl(i,j) = rmse_hgt_east_ctl saved_rmse_slp_ak_drop(i,j) = rmse_slp_ak_drop saved_rmse_slp_conus_drop(i,j) = rmse_slp_conus_drop saved_rmse_slp_global_drop(i,j) = rmse_slp_global_drop saved_rmse_slp_west_drop(i,j) = rmse_slp_west_drop saved_rmse_slp_east_drop(i,j) = rmse_slp_east_drop saved_rmse_hgt_ak_drop(i,j) = rmse_hgt_ak_drop saved_rmse_hgt_conus_drop(i,j) = rmse_hgt_conus_drop saved_rmse_hgt_global_drop(i,j) = rmse_hgt_global_drop saved_rmse_hgt_west_drop(i,j) = rmse_hgt_west_drop saved_rmse_hgt_east_drop(i,j) = rmse_hgt_east_drop ; calculate % improvement based on rmse and save saved_pct_hgt_ak(i,j) = ((rmse_hgt_ak_drop - rmse_hgt_ak_ctl) / (rmse_hgt_ak_ctl) * 100.0) ; % percent increase/reduction in drop case compared to control saved_pct_hgt_conus(i,j) = ((rmse_hgt_conus_drop - rmse_hgt_conus_ctl) / (rmse_hgt_conus_ctl) * 100.0) ; a postive value indicates increase rsme error, negative is improvement saved_pct_hgt_global(i,j) = ((rmse_hgt_global_drop - rmse_hgt_global_ctl) / (rmse_hgt_global_ctl) * 100.0) saved_pct_hgt_west(i,j) = ((rmse_hgt_west_drop - rmse_hgt_west_ctl) / (rmse_hgt_west_ctl) * 100.0) ; % percent increase/reduction in drop case compared to control saved_pct_hgt_east(i,j) = ((rmse_hgt_east_drop - rmse_hgt_east_ctl) / (rmse_hgt_east_ctl) * 100.0) ; a postive value indicates increase rsme error, negative is improvement saved_pct_slp_ak(i,j) = ((rmse_slp_ak_drop - rmse_slp_ak_ctl) / (rmse_slp_ak_ctl) * 100.0) ; % percent increase/reduction in drop case compared to control saved_pct_slp_conus(i,j) = ((rmse_slp_conus_drop - rmse_slp_conus_ctl) / (rmse_slp_conus_ctl) * 100.0) ; a postive value indicates increase rsme error, negative is improvement saved_pct_slp_global(i,j) = ((rmse_slp_global_drop - rmse_slp_global_ctl) / (rmse_slp_global_ctl) * 100.0) saved_pct_slp_west(i,j) = ((rmse_slp_west_drop - rmse_slp_west_ctl) / (rmse_slp_west_ctl) * 100.0) ; % percent increase/reduction in drop case compared to control saved_pct_slp_east(i,j) = ((rmse_slp_east_drop - rmse_slp_east_ctl) / (rmse_slp_east_ctl) * 100.0) ; a postive value indicates increase rsme error, negative is improvement end do ; end loop over forecast hours ; make plots for this GFS run of RMSE and % improvement/degradation for all 3 domains ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; plot section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; max1 = max(saved_rmse_slp_ak_ctl(i,:)) max2 = max(saved_rmse_slp_ak_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_" + region_label + "_slp_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE "+ region_label + " Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_ak_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_ak_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(saved_rmse_hgt_ak_ctl(i,:)) max2 = max(saved_rmse_hgt_ak_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_" + region_label + "_" + atmos_variable(a) + "_" + level_selector(b) + "_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE " + region_label + " Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_ak_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_ak_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;;;; ;; conus ;; max1 = max(saved_rmse_slp_conus_ctl(i,:)) max2 = max(saved_rmse_slp_conus_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_conus_slp_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE CONUS Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_conus_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_conus_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(saved_rmse_hgt_conus_ctl(i,:)) max2 = max(saved_rmse_hgt_conus_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_conus_" + atmos_variable(a) + "_" + level_selector(b) + "_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE CONUS Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_conus_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_conus_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;;;; ;; global ;;;; max1 = max(saved_rmse_slp_global_ctl(i,:)) max2 = max(saved_rmse_slp_global_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_global_slp_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE Global Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_global_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_global_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(saved_rmse_hgt_global_ctl(i,:)) max2 = max(saved_rmse_hgt_global_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_global_" + atmos_variable(a) + "_" + level_selector(b) + "_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE Global Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_global_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_global_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;;; ; west us ;;; max1 = max(saved_rmse_slp_west_ctl(i,:)) max2 = max(saved_rmse_slp_west_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if wks = gsn_open_wks ("png","rmse_west_slp_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE WEST US Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_west_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_west_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(saved_rmse_hgt_west_ctl(i,:)) max2 = max(saved_rmse_hgt_west_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_west_" + atmos_variable(a) + "_" + level_selector(b) + "_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE WEST US Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_west_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_west_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;;; ; east us ;;; max1 = max(saved_rmse_slp_east_ctl(i,:)) max2 = max(saved_rmse_slp_east_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_east_slp_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE EAST US Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_east_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_slp_east_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(saved_rmse_hgt_east_ctl(i,:)) max2 = max(saved_rmse_hgt_east_drop(i,:)) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_east_" + atmos_variable(a) + "_" + level_selector(b) + "_" + case_name + "_" + start_date(i) + utc_start(i)) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE EAST US Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_east_ctl(i,:),res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,saved_rmse_hgt_east_drop(i,:),res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do draw(plot0) frame(wks) delete(res) ;; ; end plot section for loop ;;; end do ; end loop over individual GFS runs ; average RMSE and % change values avg_rmse_slp_ak_ctl = dim_avg_n(saved_rmse_slp_ak_ctl,0) avg_rmse_slp_ak_drop = dim_avg_n(saved_rmse_slp_ak_drop,0) avg_rmse_hgt_ak_ctl = dim_avg_n(saved_rmse_hgt_ak_ctl,0) avg_rmse_hgt_ak_drop = dim_avg_n(saved_rmse_hgt_ak_drop,0) avg_rmse_slp_conus_ctl = dim_avg_n(saved_rmse_slp_conus_ctl,0) avg_rmse_slp_conus_drop = dim_avg_n(saved_rmse_slp_conus_drop,0) avg_rmse_hgt_conus_ctl = dim_avg_n(saved_rmse_hgt_conus_ctl,0) avg_rmse_hgt_conus_drop = dim_avg_n(saved_rmse_hgt_conus_drop,0) avg_rmse_slp_global_ctl = dim_avg_n(saved_rmse_slp_global_ctl,0) avg_rmse_slp_global_drop = dim_avg_n(saved_rmse_slp_global_drop,0) avg_rmse_hgt_global_ctl = dim_avg_n(saved_rmse_hgt_global_ctl,0) avg_rmse_hgt_global_drop = dim_avg_n(saved_rmse_hgt_global_drop,0) avg_rmse_slp_west_ctl = dim_avg_n(saved_rmse_slp_west_ctl,0) avg_rmse_slp_west_drop = dim_avg_n(saved_rmse_slp_west_drop,0) avg_rmse_hgt_west_ctl = dim_avg_n(saved_rmse_hgt_west_ctl,0) avg_rmse_hgt_west_drop = dim_avg_n(saved_rmse_hgt_west_drop,0) avg_rmse_slp_east_ctl = dim_avg_n(saved_rmse_slp_east_ctl,0) avg_rmse_slp_east_drop = dim_avg_n(saved_rmse_slp_east_drop,0) avg_rmse_hgt_east_ctl = dim_avg_n(saved_rmse_hgt_east_ctl,0) avg_rmse_hgt_east_drop = dim_avg_n(saved_rmse_hgt_east_drop,0) avg_pct_slp_ak = dim_avg_n(saved_pct_slp_ak,0) avg_pct_slp_conus = dim_avg_n(saved_pct_slp_conus,0) avg_pct_slp_global = dim_avg_n(saved_pct_slp_global,0) avg_pct_slp_west = dim_avg_n(saved_pct_slp_west,0) avg_pct_slp_east = dim_avg_n(saved_pct_slp_east,0) avg_pct_hgt_ak = dim_avg_n(saved_pct_hgt_ak,0) avg_pct_hgt_conus = dim_avg_n(saved_pct_hgt_conus,0) avg_pct_hgt_global = dim_avg_n(saved_pct_hgt_global,0) avg_pct_hgt_west = dim_avg_n(saved_pct_hgt_west,0) avg_pct_hgt_east = dim_avg_n(saved_pct_hgt_east,0) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; STATISTICAL SIGNIFICANCE TESTING ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; compute statistical significance over all assimilation cycles for slp ; Verification region ave_ctl_verif_region1 = avg_rmse_slp_ak_ctl ave_exp_verif_region1 = avg_rmse_slp_ak_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_slp_ak_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_slp_ak_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference over Verif. Region between CTL and EXP for variable: SLP ARE statistically different over all cycles with probability: " + prob2) else print("Difference over Verif. Region between CTL and EXP for variable: SLP NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_verif_slp = prob1 ; CONUS ave_ctl_verif_region1 = avg_rmse_slp_conus_ctl ave_exp_verif_region1 = avg_rmse_slp_conus_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_slp_conus_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_slp_conus_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference CONUS between CTL and EXP for variable: SLP ARE statistically different over all cycles with probability: " + prob2) else print("Difference CONUS between CTL and EXP for variable: SLP NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_conus_slp = prob1 ; West US ave_ctl_verif_region1 = avg_rmse_slp_west_ctl ave_exp_verif_region1 = avg_rmse_slp_west_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_slp_west_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_slp_west_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference WEST US between CTL and EXP for variable: SLP ARE statistically different over all cycles with probability: " + prob2) else print("Difference WEST US between CTL and EXP for variable: SLP NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_west_slp = prob1 ; East US ave_ctl_verif_region1 = avg_rmse_slp_east_ctl ave_exp_verif_region1 = avg_rmse_slp_east_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_slp_east_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_slp_east_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference EAST US between CTL and EXP for variable: SLP ARE statistically different over all cycles with probability: " + prob2) else print("Difference EAST US between CTL and EXP for variable: SLP NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_east_slp = prob1 ; Globally ave_ctl_verif_region1 = avg_rmse_slp_global_ctl ave_exp_verif_region1 = avg_rmse_slp_global_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_slp_global_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_slp_global_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference Globally between CTL and EXP for variable SLP ARE statistically different over all cycles with probability: " + prob2) else print("Difference Globally between CTL and EXP for variable: SLP NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_glob_slp = prob1 ;;;;;;;;;;;;;;;;;;;;;;;;;; ; need to add in level dependent fields sigificance here ; Verification region ave_ctl_verif_region1 = avg_rmse_hgt_ak_ctl ave_exp_verif_region1 = avg_rmse_hgt_ak_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_hgt_ak_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_hgt_ak_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference over Verif. Region between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa ARE statistically different over all cycles with probability: " + prob2) else print("Difference over Verif. Region between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_verif_hgt = prob1 ; CONUS ave_ctl_verif_region1 = avg_rmse_hgt_conus_ctl ave_exp_verif_region1 = avg_rmse_hgt_conus_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_hgt_conus_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_hgt_conus_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference CONUS between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa ARE statistically different over all cycles with probability: " + prob2) else print("Difference CONUS between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_conus_hgt = prob1 ; West US ave_ctl_verif_region1 = avg_rmse_hgt_west_ctl ave_exp_verif_region1 = avg_rmse_hgt_west_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_hgt_west_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_hgt_west_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference WEST US between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa ARE statistically different over all cycles with probability: " + prob2) else print("Difference WEST US between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_west_hgt = prob1 ; East US ave_ctl_verif_region1 = avg_rmse_hgt_east_ctl ave_exp_verif_region1 = avg_rmse_hgt_east_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_hgt_east_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_hgt_east_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference EAST US between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa ARE statistically different over all cycles with probability: " + prob2) else print("Difference EAST US between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_east_hgt = prob1 ; Globally ave_ctl_verif_region1 = avg_rmse_hgt_global_ctl ave_exp_verif_region1 = avg_rmse_hgt_global_drop var_ctl_verif_region1 = dim_variance_n(saved_rmse_hgt_global_ctl,0) var_exp_verif_region1 = dim_variance_n(saved_rmse_hgt_global_drop,0) ss = dimsizes(dates_plot) ; number of values in the arrays iflag = True ; population variance similar tval_opt = False ; we only want p-value prob1 = ttest(ave_ctl_verif_region1,var_ctl_verif_region1,ss,ave_exp_verif_region1,var_exp_verif_region1,ss,iflag,tval_opt) prob2 = avg(prob1) if (prob2.le.siglvl) then print("Difference Globally between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa ARE statistically different over all cycles with probability: " + prob2) else print("Difference Globally between CTL and EXP for variable: " + atmos_variable(a) + " at level: " + level_selector(b) + " hPa NOT statistically different over all cycles with probability: " + prob2) end if print("Probability for each 6hr lead time: " + prob1) pvalue_glob_hgt = prob1 ;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;; ; END OF STATISTICAL TESTING AREA ;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; avg_rmse_all_ak_slp_ctl = decimalPlaces(avg(avg_rmse_slp_ak_ctl),2,True) avg_rmse_all_ak_hgt_ctl = decimalPlaces(avg(avg_rmse_hgt_ak_ctl),2,True) avg_rmse_all_conus_slp_ctl = decimalPlaces(avg(avg_rmse_slp_conus_ctl),2,True) avg_rmse_all_conus_hgt_ctl = decimalPlaces(avg(avg_rmse_hgt_conus_ctl),2,True) avg_rmse_all_global_slp_ctl = decimalPlaces(avg(avg_rmse_slp_global_ctl),2,True) avg_rmse_all_global_hgt_ctl = decimalPlaces(avg(avg_rmse_hgt_global_ctl),2,True) avg_rmse_all_west_slp_ctl = decimalPlaces(avg(avg_rmse_slp_west_ctl),2,True) avg_rmse_all_west_hgt_ctl = decimalPlaces(avg(avg_rmse_hgt_west_ctl),2,True) avg_rmse_all_east_slp_ctl = decimalPlaces(avg(avg_rmse_slp_east_ctl),2,True) avg_rmse_all_east_hgt_ctl = decimalPlaces(avg(avg_rmse_hgt_east_ctl),2,True) avg_rmse_all_ak_slp_drop = decimalPlaces(avg(avg_rmse_slp_ak_drop),2,True) avg_rmse_all_ak_hgt_drop = decimalPlaces(avg(avg_rmse_hgt_ak_drop),2,True) avg_rmse_all_conus_slp_drop = decimalPlaces(avg(avg_rmse_slp_conus_drop),2,True) avg_rmse_all_conus_hgt_drop = decimalPlaces(avg(avg_rmse_hgt_conus_drop),2,True) avg_rmse_all_global_slp_drop = decimalPlaces(avg(avg_rmse_slp_global_drop),2,True) avg_rmse_all_global_hgt_drop = decimalPlaces(avg(avg_rmse_hgt_global_drop),2,True) avg_rmse_all_west_slp_drop = decimalPlaces(avg(avg_rmse_slp_west_drop),2,True) avg_rmse_all_west_hgt_drop = decimalPlaces(avg(avg_rmse_hgt_west_drop),2,True) avg_rmse_all_east_slp_drop = decimalPlaces(avg(avg_rmse_slp_east_drop),2,True) avg_rmse_all_east_hgt_drop = decimalPlaces(avg(avg_rmse_hgt_east_drop),2,True) avg_pct_all_ak_slp = decimalPlaces(((avg_rmse_all_ak_slp_drop - avg_rmse_all_ak_slp_ctl) / (avg_rmse_all_ak_slp_ctl) * 100.0),2,True) ;% avg_pct_all_ak_hgt = decimalPlaces(((avg_rmse_all_ak_hgt_drop - avg_rmse_all_ak_hgt_ctl) / (avg_rmse_all_ak_hgt_ctl) * 100.0),2,True) ;% avg_pct_all_conus_slp = decimalPlaces(((avg_rmse_all_conus_slp_drop - avg_rmse_all_conus_slp_ctl) / (avg_rmse_all_conus_slp_ctl) * 100.0),2,True) ;% avg_pct_all_conus_hgt = decimalPlaces(((avg_rmse_all_conus_hgt_drop - avg_rmse_all_conus_hgt_ctl) / (avg_rmse_all_conus_hgt_ctl) * 100.0),2,True) ;% avg_pct_all_global_slp = decimalPlaces(((avg_rmse_all_global_slp_drop - avg_rmse_all_global_slp_ctl) / (avg_rmse_all_global_slp_ctl) * 100.0),2,True) ;% avg_pct_all_global_hgt = decimalPlaces(((avg_rmse_all_global_hgt_drop - avg_rmse_all_global_hgt_ctl) / (avg_rmse_all_global_hgt_ctl) * 100.0),2,True) ;% avg_pct_all_west_slp = decimalPlaces(((avg_rmse_all_west_slp_drop - avg_rmse_all_west_slp_ctl) / (avg_rmse_all_west_slp_ctl) * 100.0),2,True) ;% avg_pct_all_west_hgt = decimalPlaces(((avg_rmse_all_west_hgt_drop - avg_rmse_all_west_hgt_ctl) / (avg_rmse_all_west_hgt_ctl) * 100.0),2,True) ;% avg_pct_all_east_slp = decimalPlaces(((avg_rmse_all_east_slp_drop - avg_rmse_all_east_slp_ctl) / (avg_rmse_all_east_slp_ctl) * 100.0),2,True) ;% avg_pct_all_east_hgt = decimalPlaces(((avg_rmse_all_east_hgt_drop - avg_rmse_all_east_hgt_ctl) / (avg_rmse_all_east_hgt_ctl) * 100.0),2,True) ;% ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; save variables for end of program to plot all variables on a single histogram ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; print("Saving atmospheric variables...") saved_pct_variable_ak(a,b) = avg_pct_all_ak_hgt saved_pct_variable_conus(a,b) = avg_pct_all_conus_hgt saved_pct_variable_global(a,b) = avg_pct_all_global_hgt saved_pct_variable_west(a,b) = avg_pct_all_west_hgt saved_pct_variable_east(a,b) = avg_pct_all_east_hgt ; saved slp for each region is avg_pct_all_region_slp... ;;;;;;;;; ; plot avg RMSE and % improvement over all cycles of GFS runs ;;;;;;;;; print("Plotting avg RMSE and %...") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; plot section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; max1 = max(avg_rmse_slp_ak_ctl) max2 = max(avg_rmse_slp_ak_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_" + region_label + "_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE " + region_label + " Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_ak_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_ak_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_verif_slp if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_verif_slp)-1 if (pvalue_verif_slp(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_verif_slp.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(avg_rmse_hgt_ak_ctl) max2 = max(avg_rmse_hgt_ak_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_" + region_label + "_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE " + region_label + " Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa "+ atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_ak_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_ak_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_verif_hgt if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_verif_hgt)-1 if (pvalue_verif_hgt(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_verif_hgt.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;;;; ;; conus ;; max1 = max(avg_rmse_slp_conus_ctl) max2 = max(avg_rmse_slp_conus_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_conus_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE CONUS Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_conus_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_conus_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_conus_slp if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_conus_slp)-1 if (pvalue_conus_slp(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_conus_slp.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(avg_rmse_hgt_conus_ctl) max2 = max(avg_rmse_hgt_conus_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_conus_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE CONUS Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_conus_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_conus_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_conus_hgt if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_conus_hgt)-1 if (pvalue_conus_hgt(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_conus_hgt.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;;;; ;; global ;;;; max1 = max(avg_rmse_slp_global_ctl) max2 = max(avg_rmse_slp_global_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_global_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE Global Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_global_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_global_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_glob_slp if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_glob_slp)-1 if (pvalue_glob_slp(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_glob_slp.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(avg_rmse_hgt_global_ctl) max2 = max(avg_rmse_hgt_global_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_global_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE Global Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_global_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_global_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_glob_hgt if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_glob_hgt)-1 if (pvalue_glob_hgt(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_glob_hgt.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;;;; ;; west us ;; max1 = max(avg_rmse_slp_west_ctl) max2 = max(avg_rmse_slp_west_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if wks = gsn_open_wks ("png","rmse_west_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE WEST US Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_west_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_west_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_west_slp if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_west_slp)-1 if (pvalue_west_slp(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_west_slp.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(avg_rmse_hgt_west_ctl) max2 = max(avg_rmse_hgt_west_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if wks = gsn_open_wks ("png","rmse_west_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE WEST US Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_west_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_west_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_west_hgt if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_west_hgt)-1 if (pvalue_west_hgt(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_west_hgt.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;;; ; east us ;;; max1 = max(avg_rmse_slp_east_ctl) max2 = max(avg_rmse_slp_east_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_east_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE EAST US Domain (SLP)" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = "SLP RMSE (hPa)" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_east_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_slp_east_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_east_slp if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_east_slp)-1 if (pvalue_east_slp(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_east_slp.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;; hgt ;;; max1 = max(avg_rmse_hgt_east_ctl) max2 = max(avg_rmse_hgt_east_drop) if (max1.gt.max2) then max_plot = max1 + 2.0 else max_plot = max2 + 2.0 end if ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","rmse_east_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "RMSE EAST US Domain (" + level_selector(b) + " hPa " + atmos_variable(a) + ")" res@trYMaxF = max_plot res@trYMinF = 0 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour ; res@tmXBValues = (/0,2,4,6,8,10,12/) ; res@tmXBLabels = (/"0","12","24","36","48","60","72"/) ; res@gsnXRefLine = (/4,12/) res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@tiYAxisString = level_selector(b) + " hPa " + atmos_variable(a) + " RMSE" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(0) plot0 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_east_ctl,res) res@xyLineColor = colors_plot(1) res@xyLineThicknessF = 4.0 plot1 = gsn_csm_xy(wks,x_axis,avg_rmse_hgt_east_drop,res) overlay(plot0,plot1) ; Attach a legend txres = True txres@txFontHeightF = 0.016 do k=0,dimsizes(labels)-1 txres@txFontColor = colors(k) gsn_text_ndc(wks,labels(k),xtxt2(k),ytxt2(k),txres) end do ; attach significance ;;;;;; ; determine colors for various sections of line for significance colors_signif = pvalue_east_hgt if (any(colors_signif.le.siglvl)) then colors_signif(ind(colors_signif.le.siglvl)) = 2.0 ; number of polymarkers for significance value total_good = num(colors_signif.eq.2) new_array = new((/total_good/),"float") x_new = new_array y_new = new_array z=0 do k=0,dimsizes(pvalue_east_hgt)-1 if (pvalue_east_hgt(k).le.siglvl) then x_new(z) = tofloat(x_axis(k)) y_new(z) = 0.3 z = z + 1 end if end do end if ; add polymarkers to define significant regions mres = True mres@gsMarkerIndex = 3 mres@gsMarkerSizeF = 15.0 mres@gsMarkerThicknessF = 3.0 mres@gsMarkerColor = "blue3"; blue3 if (any(pvalue_east_hgt.le.siglvl)) then gsn_polymarker(wks,plot0,x_new,y_new,mres) delete(x_new) delete(y_new) delete(new_array) delete(total_good) end if delete(colors_signif) ;;;;;;;;; draw(plot0) frame(wks) delete(res) ;;;; ;;;; ; plot % for 3 domains, hgt and slp ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; max_1 = max(avg_pct_slp_ak) + 2.0 min_1 = min(avg_pct_slp_ak) - 2.0 ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","percent_" + region_label + "_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = region_label + " Domain (SLP)" res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_slp_ak,res) draw(plot0) frame(wks) delete(res) ;; hgt ;;; max_1 = max(avg_pct_hgt_ak) + 2.0 min_1 = min(avg_pct_hgt_ak) - 2.0 ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","percent_" + region_label + "_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = region_label + " Domain " + level_selector(b) + " hPa " + atmos_variable(a) res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_hgt_ak,res) draw(plot0) frame(wks) delete(res) ;; ; conus ;; max_1 = max(avg_pct_slp_conus) + 2.0 min_1 = min(avg_pct_slp_conus) - 2.0 ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","percent_conus_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "CONUS Domain (SLP)" res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_slp_conus,res) draw(plot0) frame(wks) delete(res) ;; hgt ;;; max_1 = max(avg_pct_hgt_conus) + 2.0 min_1 = min(avg_pct_hgt_conus) - 2.0 ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","percent_conus_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "CONUS Domain " + level_selector(b) + " hPa " + atmos_variable(a) res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_hgt_conus,res) draw(plot0) frame(wks) delete(res) ;; ; global ;; max_1 = max(avg_pct_slp_global) + 2.0 min_1 = min(avg_pct_slp_global) - 2.0 ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","percent_global_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "Global Domain (SLP)" res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_slp_global,res) draw(plot0) frame(wks) delete(res) ;; hgt ;;; max_1 = max(avg_pct_hgt_global) + 2.0 min_1 = min(avg_pct_hgt_global) - 2.0 ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","percent_global_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "Global Domain " + level_selector(b) + " hPa " + atmos_variable(a) res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_hgt_global,res) draw(plot0) frame(wks) delete(res) ;;; ; west us ;;; max_1 = max(avg_pct_slp_west) + 2.0 min_1 = min(avg_pct_slp_west) - 2.0 wks = gsn_open_wks ("png","percent_west_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "WEST US Domain (SLP)" res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_slp_west,res) draw(plot0) frame(wks) delete(res) ;; hgt ;;; max_1 = max(avg_pct_hgt_west) + 2.0 min_1 = min(avg_pct_hgt_west) - 2.0 ; do one for ak, global, and conus domains wks = gsn_open_wks ("png","percent_west_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "WEST US Domain " + level_selector(b) + " hPa " + atmos_variable(a) res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_hgt_west,res) draw(plot0) frame(wks) delete(res) ;;; ; east us ;;; max_1 = max(avg_pct_slp_east) + 2.0 min_1 = min(avg_pct_slp_east) - 2.0 wks = gsn_open_wks ("png","percent_east_slp_avg_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "EAST US Domain (SLP)" res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_slp_east,res) draw(plot0) frame(wks) delete(res) ;; hgt ;;; max_1 = max(avg_pct_hgt_east) + 2.0 min_1 = min(avg_pct_hgt_east) - 2.0 wks = gsn_open_wks ("png","percent_east_" + atmos_variable(a) + "_avg_" + level_selector(b) + "_" + case_name) ; send graphics to PNG file res = True ; plot mods desired res@tiMainString = "EAST US Domain " + level_selector(b) + " hPa " + atmos_variable(a) res@trYMaxF = max_1 res@trYMinF = min_1 res@gsnFrame = False ; don't advance frame yet colors_plot = (/"Black","Red"/) ; change line color res@gsnMaximize = True res@gsnDraw = False res@trXMinF = 0 res@trXMaxF = dimsizes(dates_plot)-1 x_axis = ispan(0,dimsizes(dates_plot)-1,1) res@tmXBMode = "Explicit" ; Force labels where we want them, later as function of forecast hour res@tmXBValues = (/0,2,4,6,8,10,12,14,16/) res@tmXBLabels = (/"0","12","24","36","48","60","72","84","96"/) res@gsnXRefLine = (/4,12/) res@gsnYRefLine = 0.0 res@tiYAxisString = "% Improvement/degradation" res@tiXAxisString = "Forecast lead time (hour)" res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmXBLabelFontHeightF = 0.015 res@tmYLLabelFontHeightF = 0.015 res@xyLineThicknessF = 4.0 res@xyLineColor = colors(1) plot0 = gsn_csm_xy(wks,x_axis,avg_pct_hgt_east,res) draw(plot0) frame(wks) delete(res) ;;; ; end of plot section ;;; print("Creating table of all results...") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; create table of all results ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Header ncr2 = (/1,6/) ; 1 row, 6 columns x2 = (/0.000,1.0/) ; Start and end X y2 = (/0.920,1.0/) ; Start and end Y ; Main table body ncr3 = (/6,6/) ; 6 rows, 6 columns x3 = (/0.000,1.0/) ; Start and end X y3 = (/0.5,0.920/) ; Start and end Y text3 = (/ (/"% SLP",tostring_with_format(avg_pct_all_ak_slp,"%6.2f"),tostring_with_format(avg_pct_all_conus_slp,"%6.2f"),tostring_with_format(avg_pct_all_west_slp,"%6.2f"),tostring_with_format(avg_pct_all_east_slp,"%6.2f"),tostring_with_format(avg_pct_all_global_slp,"%6.2f")/), \ (/"% " + level_selector(b) + " " + atmos_variable(a),tostring_with_format(avg_pct_all_ak_hgt,"%6.2f"),tostring_with_format(avg_pct_all_conus_hgt,"%6.2f"),tostring_with_format(avg_pct_all_west_hgt,"%6.2f"),tostring_with_format(avg_pct_all_east_hgt,"%6.2f"),tostring_with_format(avg_pct_all_global_hgt,"%6.2f")/), \ (/"RMS SLP " + labels2(0),tostring_with_format(avg_rmse_all_ak_slp_ctl,"%6.2f"),tostring_with_format(avg_rmse_all_conus_slp_ctl,"%6.2f"),tostring_with_format(avg_rmse_all_west_slp_ctl,"%6.2f"),tostring_with_format(avg_rmse_all_east_slp_ctl,"%6.2f"),tostring_with_format(avg_rmse_all_global_slp_ctl,"%6.2f")/), \ (/"RMS SLP " + labels2(1),tostring_with_format(avg_rmse_all_ak_slp_drop,"%6.2f"),tostring_with_format(avg_rmse_all_conus_slp_drop,"%6.2f"),tostring_with_format(avg_rmse_all_west_slp_drop,"%6.2f"),tostring_with_format(avg_rmse_all_east_slp_drop,"%6.2f"),tostring_with_format(avg_rmse_all_global_slp_drop,"%6.2f")/), \ (/"RMS " + level_selector(b) + " " + atmos_variable(a) + " " + labels2(0),tostring_with_format(avg_rmse_all_ak_hgt_ctl,"%6.2f"),tostring_with_format(avg_rmse_all_conus_hgt_ctl,"%6.2f"),tostring_with_format(avg_rmse_all_west_hgt_ctl,"%6.2f"),tostring_with_format(avg_rmse_all_east_hgt_ctl,"%6.2f"),tostring_with_format(avg_rmse_all_global_hgt_ctl,"%6.2f")/), \ (/"RMS " + level_selector(b) + " " + atmos_variable(a) + " " + labels2(1),tostring_with_format(avg_rmse_all_ak_hgt_drop,"%6.2f"),tostring_with_format(avg_rmse_all_conus_hgt_drop,"%6.2f"),tostring_with_format(avg_rmse_all_west_hgt_drop,"%6.2f"),tostring_with_format(avg_rmse_all_east_hgt_drop,"%6.2f"),tostring_with_format(avg_rmse_all_global_hgt_drop,"%6.2f")/)/) table_file = "table_RMSE_percent_" + atmos_variable(a) + "_" + level_selector(b) + "_" + case_name wks = gsn_open_wks("png",table_file) ; send graphics to PNG file ; ; Header ; res2 = True res2@txFontHeightF = 0.015 res2@txFont = "helvetica-bold" res2@gsFillColor = (/"cornflowerblue","tomato3","tomato3","tomato3","tomato3","tomato3"/) gsn_table(wks,ncr2,x2,y2,text2,res2) ; ; Main body of table. ; res3 = True ; Set up resource list res3@txFontHeightF = 0.014 gsn_table(wks,ncr3,x3,y3,text3,res3) frame(wks) ; Advance the frame. delete(res2) delete(res3) system("convert -trim " + table_file + ".png " + table_file + ".png") end do ; end loop over pressure levels end do ; end loop over atmos_variable ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; make histograms of all results ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; print("Making histogram plots of all results...") y_values = ispan(2,62,2) hist_names = (/"Verification-Region","CONUS","Global","WEST-US","EAST-US"/) variable_names = (/"Sea Level Pressure","200 mb Height","300 mb Height","500 mb Height", "700 mb Height", "850 mb Height", "925 mb Height", \ "200 mb U Wind","300 mb U Wind","500 mb U Wind","700 mb U Wind","850 mb U Wind","925 mb U Wind", \ "200 mb V Wind","300 mb V Wind","500 mb V Wind","700 mb V Wind","850 mb V Wind","925 mb V Wind", \ "200 mb RH","300 mb RH","500 mb RH","700 mb RH","850 mb RH","925 mb RH", \ "200 mb Temp","300 mb Temp","500 mb Temp","700 mb Temp","850 mb Temp","925 mb Temp"/) ; variable names for histogram table color_names = (/"antiquewhite4","aquamarine3","azure3","bisque3", "blue", "blueviolet", "brown3", \ "burlywood3","cadetblue","olivedrab","orange2","orangered2","orchid3", \ "palegoldenrod","palegreen3","paleturquoise4","palevioletred4","peachpuff3","peru", \ "chartreuse3","chocolate","coral3","cornflowerblue","cornsilk4","cyan4", \ "darkgoldenrod1","darkolivegreen","steelblue3","darksalmon","darkviolet","firebrick3"/) ; variable colors for bars in plot, separate color for each variable ; combine variables into one array for min/max determination all_pct_variables = (/saved_pct_variable_ak,saved_pct_variable_conus,saved_pct_variable_global,saved_pct_variable_west,saved_pct_variable_east/) all_pct_slps = (/avg_pct_all_ak_slp,avg_pct_all_conus_slp,avg_pct_all_global_slp,avg_pct_all_west_slp,avg_pct_all_east_slp/) do i=0,4 ; loop over regions ; determine range of max x values for plot max_rmse_value_1 = max(all_pct_variables(i,:,:)) max_rmse_value_2 = max(all_pct_slps(i)) min_rmse_value_1 = min(all_pct_variables(i,:,:)) min_rmse_value_2 = min(all_pct_slps(i)) if (max_rmse_value_1.gt.max_rmse_value_2) then max_x_value = max_rmse_value_1 else max_x_value = max_rmse_value_2 end if if (min_rmse_value_1.gt.min_rmse_value_2) then min_x_value = min_rmse_value_1 else min_x_value = min_rmse_value_2 end if ;;;;;; ; main plot resources ;;;;;;; wks = gsn_open_wks("png",hist_names(i) + "_histogram") ; send graphics to PNG file gsn_define_colormap(wks,"temp1") res = True res@gsnMaximize = True res@gsnFrame = False res@gsnDraw = False res@tiMainString = "% RMSE " + hist_names(i) res@tiMainFontHeightF = 0.03 ;res@vpWidthF = 0.3 ;res@vpHeightF = 0.9 ;res@vpXF = 0.01 ;res@vpYF = 0.99 res@trXMinF = 2.3*min_x_value ; bound the axes by max/min res@trXMaxF = 2.3*max_x_value res@trYMinF = 0 res@trYMaxF = 64 res@gsnTickMarksPointOutward = True ; turn off Y axis tickmarks res@tmYLMajorLengthF = 0 res@tmYLMajorOutwardLengthF = 0.01 ;res@tmXBMajorLenghF = 0.01 ;res@tmXBMajorOutwardLengthF = 0.01 ; create blank plot without grid lines plot_without_ygrid = gsn_blank_plot(wks,res) ; create blank plot with X grid lines res@tmYMmajorGrid = False res@tmYLOn = False res@tmYMajorGridLineDashPattern = 2 ; dashed lines res@tmYMajorGridThicknessF = 1.0 ; default is 2 res@tiXAxisString = "% RMSE Change" res@gsnXRefLine = 0.0 ; reference line to denote pos/neg change in RMSE plot_with_ygrid = gsn_blank_plot(wks,res) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; combine all atmospheric variables into one-dimensional to loop over values an easier way temp_variable_1 = all_pct_variables(i,:,:) temp_variable_2 = ndtooned(temp_variable_1) ; one dimensonal array of all values, should be 31 values temp_variable_3 = array_append_record(all_pct_slps(i),temp_variable_2,0) ; contains all 31 values for this verification region, one dimension temp_variable_3 = temp_variable_3(::-1) ; reverse order of arrays variable_names = variable_names(::-1) do a=0,dimsizes(temp_variable_3)-1 ; loop over all values for this region xbar = new(5,float) ybar = new(5,float) dum2 = new(dimsizes(temp_variable_3),graphic) ;---Set some resources for the bars. bres = True bres@gsEdgesOn = True ; Outline the polygons (bars) ; ; Loop through each city and attach bar to plot that ; doesn't have grid lines. ;---Do longer bar first. bres@gsFillColor = color_names(a) xbar = (/0,temp_variable_3(a),temp_variable_3(a),0,0/) ; values of the % RMSE in x-direction ybar = (/y_values(a)-0.25,y_values(a)-0.25,y_values(a)+0.25,y_values(a)+0.25,y_values(a)-0.25/) ; width of the bars in y-direction dum2(a) = gsn_add_polygon(wks,plot_without_ygrid,xbar,ybar,bres) ; text labels next to bars txres = True ; text mods desired txres@txFontHeightF = 0.0112 ; default size is HUGE! txres@txFont= "helvetica-bold" txres@txAngleF = 0 ; text angle txres@txJust = "CenterLeft" ; puts text on top of bars txres@txFontColor = "black" if (temp_variable_3(a).lt.0) then gsn_text(wks,plot_without_ygrid,variable_names(a),temp_variable_3(a)-0.015*(2.3*min_x_value),y_values(a),txres) else gsn_text(wks,plot_without_ygrid,variable_names(a),temp_variable_3(a)+0.015*(2.3*max_x_value),y_values(a),txres) end if end do ; end loop over variables to plot in histogram draw(plot_with_ygrid) draw(plot_without_ygrid) frame(wks) delete(wks) delete(res) delete(txres) delete(temp_variable_1) delete(temp_variable_2) delete(temp_variable_3) exit end do ; end loop over regions print("Program ended successfully") end