; Example of using panels with WRF data load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" load "WRFUserARW.ncl" external GET_IJ "./get_ij.so" begin dx = 1.3 ; set bounding box ; bounding box for January 23 2016 ; lats = (/ 39.00, 43.50 /) ; lons = (/ -77.0, -69.50 /) ; bounding box for January 23 2016 lats = (/ 39.00, 42.50 /) lons = (/ -76.0, -70.0/) ; bounding box for February 08 2013 ; lats = (/ 38.00, 42.50 /) ; lons = (/ -77.0, -68.75 /) ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. ; a = addfile("../wrfout_d01_2000-01-24_12:00:00.nc","r") ;Type = "Snow" ;Type = "Hail" Type = "Total" height = 1000. elev = 0.5 ; height = 5000. ; height = 6000. ; height = 2500. heightkm = height/1000. ; We generate plots, but what kind do we prefer? ; type = "x11" ; type = "pdf" ; type = "ps" ; type = "ncgm" date = "02_08_2013" ; out_name = "wrf_refl_dates_noass" + date + "_" + heightkm +"_km" ; out_name = "wrf_" + Type + "_refl_cont_" + date + "_" + heightkm +"_km" ; out_name = "wrf_" + Type + "_refl_graup_only_" + date + "_" + heightkm +"_km" ; out_name = "wrf_" + Type + "_refl_graupel_hail_" + date + "_" + heightkm +"_km" ; out_name = "wrf_" + Type + "_refl_graupel_hail_43_" + date + "_" + heightkm +"_km" out_name = "wrf_radar_" + date + "_" + elev +"_deg" wtype = "png" ; wtype = "pdf" wtype@wkWidth = 7500 ; produces a better looking PNG wtype@wkHeight = 7500 wks = gsn_open_wks(wtype,out_name) gsn_define_colormap(wks,"rainbow") ; overwrite the .hluresfile color map ; Set some basic resources ;---Set some text resources txres = True txres@txFontHeightF = 0.01 txres@txPosXF = 0.1 res = True res@NoHeaderFooter = True ; Switch headers and footers off pltres = True pltres@PanelPlot = True ; Indicate these plots are to be paneled. mpres = True mpres = True mpres@mpFillOn = False mpres@mpDataBaseVersion = "MediumRes" ; mpres@mpDataBaseVersion = "HighRes" mpres@mpDataSetName = "Earth..2" mpres@mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpGridLineColor = 0 mpres@mpLimbLineColor = "Black" mpres@mpPerimLineColor = "Black" mpres@mpPerimLineThicknessF = 20 mpres@mpUSStateLineThicknessF = 30.0 mpres@mpGeophysicalLineThicknessF = 30.0 ; get xlat,xlon from model for interpolation ; filename="/data4/wrf/JULY_10_SBM/wrfout_d01_2013-07-10_12:00:00_SBM_CONT" filename="wrfout_d02_2013-02-08_15:30:00" a = addfile(filename,"r") print(a) time = 0 xlat2d_m=wrf_user_getvar(a,"XLAT",time) xlon2d_m=wrf_user_getvar(a,"XLONG",time) xland=wrf_user_getvar(a,"XLAND",time) dims2d = dimsizes(xlat2d_m) ; latitude longitude of Upton, NY lat_beg = 40.8682 lon_beg = 72.8792 i_loc = 1 j_loc = 1 GET_IJ::get_ij(xlat2d_m,xlon2d_m,lat_beg,lon_beg,i_loc,j_loc,dims2d(0),dims2d(1)) i_point = i_loc-1 j_point = j_loc-1 ; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y) ; subtract one since we want to use it as an index in NCL loc = wrf_user_ll_to_ij(a, lons, lats, True) x_start = loc(0,0) - 1 x_end = loc(0,1) - 1 y_start = loc(1,0) - 1 y_end = loc(1,1) - 1 mpres@ZoomIn = True mpres@Xstart = x_start mpres@Ystart = y_start mpres@Xend = x_end mpres@Yend = y_end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; times = wrf_user_getvar(a,"times",-1) ; get all times in the file print("times = " + times) ntimes=dimsizes(times) plots_dbz = new ( 6, graphic ) plots_zdr = new ( 6, graphic ) plots_cc = new ( 6, graphic ) plots_dp = new ( 6, graphic ) ;-rw-r--r--. 1 wrf games 27995692 Apr 19 17:29 KOKX_V06_20170314_160117_pnt_out.json ;rw-r--r--. 1 wrf games 26174552 Apr 19 17:29 KOKX_V06_20170314_165819_pnt_out.json ;rw-r--r--. 1 wrf games 22986338 Apr 19 17:29 KOKX_V06_20170314_180259_pnt_out.json ;rw-r--r--. 1 wrf games 18712708 Apr 19 17:29 KOKX_V06_20170314_190018_pnt_out.json ;rw-r--r--. 1 wrf games 16111132 Apr 19 17:29 KOKX_V06_20170314_200357_pnt_out.json ;rw-r--r--. 1 wrf games 13916104 Apr 19 17:29 KOKX_V06_20170314_205442_pnt_out.json ;rw-r--r--. 1 wrf games 29713148 Apr 19 17:34 KBOX_V06_20170314_155812_pnt_out.json ;rw-r--r--. 1 wrf games 29440092 Apr 19 17:34 KBOX_V06_20170314_170158_pnt_out.json ;rw-r--r--. 1 wrf games 28797564 Apr 19 17:35 KBOX_V06_20170314_175956_pnt_out.json ;rw-r--r--. 1 wrf games 26374882 Apr 19 17:35 KBOX_V06_20170314_185755_pnt_out.json ;rw-r--r--. 1 wrf games 23078016 Apr 19 17:35 KBOX_V06_20170314_200140_pnt_out.json ;rw-r--r--. 1 wrf games 19020678 Apr 19 17:35 KBOX_V06_20170314_210051_pnt_out.json itt = 0 ; do itime = 24,ntimes-2,2 do itime = 17,22 print("Working on time = " + times(itime)) it = itime if (itt.eq.0)then date = "0000 UTC 02 09 2013" end if if (itt.eq.1)then date = "0030 UTC 02 09 2013" end if if (itt.eq.2)then date = "0100 UTC 02 09 2013" end if if (itt.eq.3)then date = "0130 UTC 02 09 2013" end if if (itt.eq.4)then date = "0200 UTC 02 09 2013" end if if (itt.eq.5)then date = "0230 UTC 02 09 2013" end if if (itt.eq.6)then date = "2100 UTC 02 09 2013" end if if (itt.eq.7)then date = "2130 UTC 02 09 2013" end if if (itt.eq.8)then date = "2200 UTC 02 09 2013" end if if (itt.eq.9)then date = "2230 UTC 02 09 2013" end if if (itt.eq.10)then date = "2300 UTC 02 09 2013" end if if (itt.eq.11)then date = "2330 UTC 02 09 2013" end if radar = wrf_user_getvar(a,"fft_zh",it) zdr = wrf_user_getvar(a,"fft_zdr",it) crs = wrf_user_getvar(a,"fft_crs",it) kdp = wrf_user_getvar(a,"fft_kdp",it) ; refl_max = dim_max_n(radar,0) z = wrf_user_getvar(a, "z",it) ; grid point height refl_zkm=wrf_user_intrp3d(radar,z,"h",height,0.,False) print("itime = " + itime) ; t = wrf_user_getvar(a,"T2",itime) ; T2 in Kelvin ; slp = wrf_user_getvar(a,"slp",itime) ; slp do j=0,dims2d(0)-1 do i=0,dims2d(1)-1 radius = sqrt((j-j_point)^2 + (i-i_point)^2) if (radius/dx.le.150)then continue else refl_zkm(j,i)= 0 end if end do end do ; These long_name and units attributes are not necessary. I ; did this so I got better output from print statements below. ; ; printVarSummary(obs_o) ; printVarSummary(obs_o) ; printVarSummary(lat_o) ; printVarSummary(lon_o) ; printMinMax(obs_o,0) ; printMinMax(lat_o,0) ; printMinMax(lon_o,0) ;---Get slice of lat/lon data at halfway points. grid=refl_zkm grid@lat2d = xlat2d_m ; this is for plotting grid@lon2d = xlon2d_m ; using gsn_csm routines printVarSummary(grid) printMinMax(grid,0) ; Generate contours. res1 = res res1@cnFillOn = True res1@lbLabelBarOn = False ; Turn off individual label bars so we can res1@cnFillOn = True ; res1@cnFillPalette = "NCV_jet" ; set color map res1@cnLinesOn = False res1@cnFillMode = "RasterFill" ; this mode is fastest res1@cnLineLabelsOn = False ; do not use line labels res1@cnFillOn = True ; color fill res1@cnLinesOn = False ; do not draw contour lines ; cmap = read_colormap_file("BlueYellowRed") cmap = read_colormap_file("NCV_jet") cmap(0,:) = 1. ; Make first color white ; res1@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels ; res1@cnLevels = (/0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1/) ; 10 levels res1@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels res1@cnLevels = ispan(-10,55,5) ; set the contour levels res1@cnFillPalette = cmap ; If you want the first color to be transparent rather than white, then use: ; cmap(0,3) = 0,0 ; Make first color transparent var_zoom = grid(y_start:y_end,x_start:x_end) contour = wrf_contour(a,wks,var_zoom,res1) ; slp_res = res ; slp_res@ContourParameters = (/ 990., 1026., 2. /) ; slp_res@cnInfoLabelFontHeightF = 0.027 ; con_slp = wrf_contour(a,wks,slp,slp_res) ; Overlay contours on a map ;pltres@NoTitles = True pltres@CommonTitle = True ; pltres@PlotTitle = times(itime) pltres@PlotTitle = date ; plots(itime) = wrf_map_overlays(a,wks,(/contour,con_slp/),pltres,mpres) plots_dbz(itt) = wrf_map_overlays(a,wks,(/contour/),pltres,mpres) itt = itt+1 delete(contour) end do ; Panel the WRF plots. pnlres = True ; pnlres@txString = t@description + " (" + t@units + ")" pnlres@txString = "Polar " + Type + " Reflectivity (dbZ)" + " at " + heightkm + " km" pnlres@gsnPanelYWhiteSpacePercent = 13 ; Add white space b/w plots. pnlres@gsnPanelLabelBar = True ; Turn on common labelbar pnlres@lbLabelAutoStride = True ; Spacing of lbar labels. pnlres@lbBoxMinorExtentF = 0.13 pnlres@lbLabelFontHeightF = .015 pnlres@lbPerimOn = False pnlres@lbLabelFont = "Helvetica" ; label font pnlres@lbTitleOn = True ; turn on title pnlres@lbTitleString = "dbZ" ; title string ; pnlres@lbTitlePosition = "Top" ; title position pnlres@lbTitlePosition = "Bottom" ; title position pnlres@lbTitleFontHeightF= .015 ; make title smaller pnlres@lbTitleDirection = "Across" ; title direction pnlres@pmLabelBarWidthF = 0.5 pnlres@pmLabelBarHeightF = 0.15 pnlres@lbTopMarginF = 0.001 pnlres@pmLabelBarOrthogonalPosF = 0.02 ; position wrt plot pnlres@lbTitleOffsetF = -0.4 gsn_panel(wks,(/plots_dbz/),(/2,3/),pnlres) end