; Compare the observed precipitation for the Feb 21 storm of the third OSE case in 2016 to the GFS runs ; Plotting differences to see comparisons better between CTL, DROP, DROP_noNPP, and noNPP ; ; Uses the Stage 4 analysis based on multi-sensor hourly/6-hourly Stage III analyses ; on local 4km polar-stereographic grids produced by the 12 River Forecast Centers (RFCs) in CONUS ; and one RFC in Alaska. The Alaska dataset is not included in the CONUS so it is a separate dataset ; Stage 4 includes manual QC performed on stage III data. Data has hourly, 6-hourly, and 24 accumlated, but we will focus ; on 6-hourly since that is what GFS has to easily compare to ; AK Precip is Level 3 ; Downloaded from NCAR ; Website: https://data.eol.ucar.edu/dataset/21.093 ; ; Appropriate Citation if used in publication or presentation: ; ; Lin, Y. 2011. GCIP/EOP Surface: Precipitation NCEP/EMC 4KM Gridded Data (GRIB) Stage IV Data, ; Version 1.0. UCAR/NCAR - Earth Observing Laboratory. http://data.eol.ucar.edu/dataset/21.093. Accessed 23 Sep 2016. ; Also: Data provided by NCAR/EOL under sponsorship of the National Science Foundation. http://data.eol.ucar.edu/ ; Andrew Kren ; October 26, 2016 ;;;;;;;;;;;;;;;;;;;;;;;;;;;; begin ; user's settings ;;;;;;;;;;;;;;;;;;;;;;;; start_date = "20160221" ; forecast initiated date utc_start = "18" ; forecast initialization hour interval = 6 ; output every 6 hours end_date = "20160225" ; forecast end date end_utc = "18" ; forecast end hour region = "alaska" ; region to plot precipitation ; for alaska, set to "alaska" ; for conus, set to "conus" ; GFS model run settings gfs_name = "DROP" ; output title for plots, this is your case name ; CTL, noNPP, DROP, DROP_noNPP GFS_DIR = "/scratch4/BMC/shout/ptmp/Andrew.Kren/pre2016c3_corr/" ; or pre2016c3_corr, pre2016ctl, pre2016c3_corr_no_npp, pre2016_no_npp err = NhlGetErrorObjectId() ; these 4 lines suppress grib warnings amongst others setvalues err "errLevel" : "Fatal" ; only report Fatal errors end setvalues ; read in higher resolution grid for lat and lon wfile = "/scratch4/BMC/shout/Andrew.Kren/ECMWF/reanalysis/ecmwf_reanalysis_sfc_highres.nc" b = addfile(wfile,"r") lat_gfs = b->latitude ; 0.25 degree x 0.25 grid, same as the pgbq*gfs output files lon_gfs = b->longitude lat_gfs = lat_gfs(::-1) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (region.eq."conus") then filename_convention = "ST4" ; filename convection of grib or netcdf files else filename_convention = "AK" end if obs_dir = "/scratch4/BMC/shout/Andrew.Kren/obs_precip/" + region + "/6hourly/" if (region.eq."conus") then lat_low = 32 ; coordinates used to plot the verification region for CA lat_high = 46 lon_low = 231 lon_high = 244 else lat_low = 55 ; coordinates used to plot the verification region for AK lat_high = 65 lon_low = 193 lon_high = 220 end if if (region.eq."conus") then ; defining array sizes of each dataset nlat = 881 nlon = 1121 else nlat = 530 nlon = 460 end if ; --------------------------- end of user's settings ----------------------------------------- fcst_hr = ispan(0,96,6) ; forecast hours to read year_st = stringtointeger(str_get_cols(start_date,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+utc_start)) indice_end = ind(dates.eq.(end_date+end_utc)) ; add dates 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 accum_pcp = new((/fc_files,nlat,nlon/),float) ; accumulated precipitation over the total duration selected above accum_pcp!0="time" accum_pcp!1="lat" accum_pcp!2="lon" accum_pcp@_FillValue = 1e20 hourly_pcp = accum_pcp do i=1,fc_files-1 ; loop over all daily files print("Processing GFS Init Run: " + tostring(dates_plot(i-1))) print("Reading file: " + filename_convention + "." + dates_plot(i) + ".06h") w = addfile(obs_dir+filename_convention + "." + dates_plot(i) + ".06h","r") ; read obs data pcp = w->A_PCP_GDS5_SFC_acc6h ; mm units lat = w->g5_lat_0 lon = w->g5_lon_1 delete(w) ; read gfs filename for this particular time period if (fcst_hr(i).lt.10) then filename = GFS_DIR + "pgbq0" + fcst_hr(i) + ".gfs." + start_date + utc_start + ".grb2" else filename = GFS_DIR + "pgbq" + fcst_hr(i) + ".gfs." + start_date + utc_start + ".grb2" end if print("Reading GFS file: " + filename) ; read gfs data w = addfile(filename,"r") if (fcst_hr(i).lt.5) then pcp_gfs = w->A_PCP_GDS0_SFC_acc ; mm units pcp_gfs = pcp_gfs(::-1,:) else pcp_gfs = w->A_PCP_GDS0_SFC_acc6h ; mm units pcp_gfs = pcp_gfs(::-1,:) end if delete(w) print("Regridding GFS to observed grid...") ; regrid GFS data to regional conus or Alaska grid regrid_gfs_pcp = rgrid2rcm_Wrap(lat_gfs,lon_gfs,pcp_gfs,lat,lon,1) ; compute difference in obs and GFS hourly_pcp(i,:,:) = (/regrid_gfs_pcp-pcp/) ; save pcp print(hourly_pcp(i,:,:)) exit ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; beginning of plot section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (region.eq."conus") then ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; outfile = "6hr_obs_RFC_pcp_" + region + "_init_" + start_date + utc_start + "_"+ dates_plot(i) wks = gsn_open_wks("png",outfile) gsn_define_colormap(wks,"WhBlGrYeRe") res = True res@gsnAddCyclic = False res@pmTickMarkDisplayMode = "Always" ; turn on fancy tickmarks res@sfXArray = lon res@sfYArray = lat res@trGridType = "TriangularMesh" res@mpProjection = "LambertConformal" res@mpLimitMode = "Corners" res@mpLeftCornerLatF = 50 res@mpLeftCornerLonF = 225 res@mpRightCornerLatF = 20 res@mpRightCornerLonF = 283 res@mpLambertParallel1F = 25 res@mpLambertParallel2F = 49 res@mpLambertMeridianF = 255 res@cnLinesOn = False res@gsnMaximize = True res@cnFillOn = True ; color plot desired res@cnLineLabelsOn = False ; turn off contour lines res@tiMainString = "6-hr Precipitation (mm) ending on: " + valid_month(i) + "/" + valid_day(i) + "/" + valid_year(i) + " " + valid_hour(i) + "Z" res@mpMinLatF = 24 ; specify min lat res@mpMaxLatF = 52 res@mpMinLonF = 230 res@mpMaxLonF = 295 res@mpCenterLonF = 245 res@gsnRightString = "" res@gsnLeftString = "" res@mpGridAndLimbOn = True res@mpGridLatSpacingF = 10.0 res@mpGridLonSpacingF = 10.0 res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ;res@mpOutlineBoundarySets = "National" ; turn on country boundaries res@mpGeophysicalLineThicknessF = 2.5 ; thickness of boundaries res@mpUSStateLineThicknessF = 2.5 res@gsnFrame = False ; don't advance frame res@gsnDraw = False ; don't draw plot res@cnFillOpacityF = 0.9 res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 2 ; set the minimum contour level res@cnMaxLevelValF = 60 ; set the maximum contour level res@cnLevelSpacingF = 4 ; set the interval between contours res@lbOrientation = "vertical" res@pmLabelBarOrthogonalPosF = -0.01 ; move label bar closer to y-axis of map plot = gsn_csm_contour_map(wks,pcp,res) gres = True gres@gsLineColor = "Red" gres@gsLineThicknessF = 5.0 lat2 = (/ lat_low, lat_low, lat_high, lat_high, lat_low/) lon2 = (/lon_low, lon_high, lon_high, lon_low, lon_low/) ;dum = gsn_add_polyline(wks,plot,lon2,lat2,gres) draw(plot) frame(wks) ; advance frame delete(res) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; else ; region is alaska ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; outfile = "6hr_obs_RFC_pcp_" + region + "_init_" + start_date + utc_start + "_"+ dates_plot(i) wks = gsn_open_wks("png",outfile) gsn_define_colormap(wks,"WhBlGrYeRe") res = True res@gsnAddCyclic = False res@pmTickMarkDisplayMode = "Always" ; turn on fancy tickmarks res@sfXArray = lon res@sfYArray = lat res@trGridType = "TriangularMesh" res@mpProjection = "LambertConformal" res@mpLimitMode = "Corners" res@mpLeftCornerLatF = 80 res@mpLeftCornerLonF = 150 res@mpRightCornerLatF = 40 res@mpRightCornerLonF = 230 res@mpLambertParallel1F = 80 res@mpLambertParallel2F = 40 res@mpLambertMeridianF = 195 res@cnLinesOn = False res@gsnMaximize = True res@cnFillOn = True ; color plot desired res@cnLineLabelsOn = False ; turn off contour lines res@tiMainString = "6-hr Precip. Diff. (mm) ending on: " + valid_month(i) + "/" + valid_day(i) + "/" + valid_year(i) + " " + valid_hour(i) + "Z" res@mpMinLatF = 24 ; specify min lat res@mpMaxLatF = 52 res@mpMinLonF = 230 res@mpMaxLonF = 295 res@mpCenterLonF = 245 res@gsnRightString = "" res@gsnLeftString = "" res@mpGridAndLimbOn = True res@mpGridLatSpacingF = 10.0 res@mpGridLonSpacingF = 10.0 ;res@mpOutlineBoundarySets = "GeophysicalAndUSStates" res@mpOutlineBoundarySets = "National" ; turn on country boundaries res@mpGeophysicalLineThicknessF = 2.5 ; thickness of boundaries res@mpUSStateLineThicknessF = 2.5 res@gsnFrame = False ; don't advance frame res@gsnDraw = False ; don't draw plot res@cnFillOpacityF = 0.9 res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 2 ; set the minimum contour level res@cnMaxLevelValF = 60 ; set the maximum contour level res@cnLevelSpacingF = 4 ; set the interval between contours res@lbOrientation = "vertical" res@pmLabelBarOrthogonalPosF = -0.01 ; move label bar closer to y-axis of map plot = gsn_csm_contour_map(wks,hourly_pcp(i,:,:),res) gres = True gres@gsLineColor = "Red" gres@gsLineThicknessF = 5.0 lat2 = (/ lat_low, lat_low, lat_high, lat_high, lat_low/) lon2 = (/lon_low, lon_high, lon_high, lon_low, lon_low/) dum = gsn_add_polyline(wks,plot,lon2,lat2,gres) draw(plot) frame(wks) ; advance frame delete(res) exit ;;;;;;;;; end if ;;;;;;;;; delete([/pcp/]) end do print("Calculating accumulated precipitation") ntstrt = 0 ; counter variables for accumulated precipitation ntlast = 0 hourly_pcp(0,:,:) = 0.0 ; set to zero for accurate accumulation do t=0,fc_files-1 accum_pcp(t,:,:) = dim_sum_n(hourly_pcp(ntstrt:ntlast,:,:),0) ntlast = ntlast + 1 ; increment for accumulated precipitation end do print("Making accumulated pcp plots...") ;;;;;;;;;;;;;;;;;;; ; Plot accumulated precipitation ;;;;;;;;;;;;;;;;;;; do i=0,fc_files-1 if (region.eq."conus") then outfile = "accum_obs_RFC_pcp_" + region + "_init_" + start_date + utc_start + "_"+ dates_plot(i) wks = gsn_open_wks("png",outfile) setvalues NhlGetWorkspaceObjectId() "wsMaximumSize" : 500000000 end setvalues gsn_define_colormap(wks,"WhBlGrYeRe") res = True res@gsnAddCyclic = False res@pmTickMarkDisplayMode = "Always" ; turn on fancy tickmarks res@sfXArray = lon res@sfYArray = lat res@trGridType = "TriangularMesh" res@mpProjection = "LambertConformal" res@mpLimitMode = "Corners" res@mpLeftCornerLatF = 50 res@mpLeftCornerLonF = 225 res@mpRightCornerLatF = 20 res@mpRightCornerLonF = 283 res@mpLambertParallel1F = 25 res@mpLambertParallel2F = 49 res@mpLambertMeridianF = 255 res@cnLinesOn = False res@gsnMaximize = True res@cnFillOn = True ; color plot desired res@cnLineLabelsOn = False ; turn off contour lines res@tiMainString = "Accum. Precipitation (mm) ending on: " + valid_month(i) + "/" + valid_day(i) + "/" + valid_year(i) + " " + valid_hour(i) + "Z" res@mpMinLatF = 24 ; specify min lat res@mpMaxLatF = 52 res@mpMinLonF = 230 res@mpMaxLonF = 295 res@mpCenterLonF = 245 res@gsnRightString = "" res@gsnLeftString = "" res@mpGridAndLimbOn = True res@mpGridLatSpacingF = 10.0 res@mpGridLonSpacingF = 10.0 res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ;res@mpOutlineBoundarySets = "National" ; turn on country boundaries res@mpGeophysicalLineThicknessF = 2.5 ; thickness of boundaries res@mpUSStateLineThicknessF = 2.5 res@gsnFrame = False ; don't advance frame res@gsnDraw = False ; don't draw plot res@cnFillOpacityF = 0.9 res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 2 ; set the minimum contour level res@cnMaxLevelValF = 120 ; set the maximum contour level res@cnLevelSpacingF = 4 ; set the interval between contours res@lbOrientation = "vertical" res@pmLabelBarOrthogonalPosF = -0.01 ; move label bar closer to y-axis of map plot = gsn_csm_contour_map(wks,accum_pcp(i,:,:),res) gres = True gres@gsLineColor = "Red" gres@gsLineThicknessF = 5.0 lat2 = (/ lat_low, lat_low, lat_high, lat_high, lat_low/) lon2 = (/lon_low, lon_high, lon_high, lon_low, lon_low/) ;dum = gsn_add_polyline(wks,plot,lon2,lat2,gres) draw(plot) frame(wks) ; advance frame delete(res) else ;;;; outfile = "accum_obs_RFC_pcp_" + region + "_init_" + start_date + utc_start + "_"+ dates_plot(i) wks = gsn_open_wks("png",outfile) setvalues NhlGetWorkspaceObjectId() "wsMaximumSize" : 500000000 end setvalues gsn_define_colormap(wks,"WhBlGrYeRe") res = True res@gsnAddCyclic = False res@pmTickMarkDisplayMode = "Always" ; turn on fancy tickmarks res@sfXArray = lon res@sfYArray = lat res@trGridType = "TriangularMesh" res@mpProjection = "LambertConformal" res@mpLimitMode = "Corners" res@mpLeftCornerLatF = 80 res@mpLeftCornerLonF = 150 res@mpRightCornerLatF = 40 res@mpRightCornerLonF = 230 res@mpLambertParallel1F = 80 res@mpLambertParallel2F = 40 res@mpLambertMeridianF = 195 res@cnLinesOn = False res@gsnMaximize = True res@cnFillOn = True ; color plot desired res@cnLineLabelsOn = False ; turn off contour lines res@tiMainString = "Accum. Precipitation (mm) ending on: " + valid_month(i) + "/" + valid_day(i) + "/" + valid_year(i) + " " + valid_hour(i) + "Z" res@mpMinLatF = 24 ; specify min lat res@mpMaxLatF = 52 res@mpMinLonF = 230 res@mpMaxLonF = 295 res@mpCenterLonF = 245 res@gsnRightString = "" res@gsnLeftString = "" res@mpGridAndLimbOn = True res@mpGridLatSpacingF = 10.0 res@mpGridLonSpacingF = 10.0 ;res@mpOutlineBoundarySets = "GeophysicalAndUSStates" res@mpOutlineBoundarySets = "National" ; turn on country boundaries res@mpGeophysicalLineThicknessF = 2.5 ; thickness of boundaries res@mpUSStateLineThicknessF = 2.5 res@gsnFrame = False ; don't advance frame res@gsnDraw = False ; don't draw plot res@cnFillOpacityF = 0.9 res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 2 ; set the minimum contour level res@cnMaxLevelValF = 120 ; set the maximum contour level res@cnLevelSpacingF = 4 ; set the interval between contours res@lbOrientation = "vertical" res@pmLabelBarOrthogonalPosF = -0.01 ; move label bar closer to y-axis of map plot = gsn_csm_contour_map(wks,accum_pcp(i,:,:),res) gres = True gres@gsLineColor = "Red" gres@gsLineThicknessF = 5.0 lat2 = (/ lat_low, lat_low, lat_high, lat_high, lat_low/) lon2 = (/lon_low, lon_high, lon_high, lon_low, lon_low/) dum = gsn_add_polyline(wks,plot,lon2,lat2,gres) draw(plot) frame(wks) ; advance frame delete(res) end if end do print("Program ended successfully") end