; 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" begin ; ; 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") ; We generate plots, but what kind do we prefer? ; type = "x11" type = "pdf" ; type = "ps" ; type = "ncgm" date = "01232016" out_name = "obs_refl_" + date wks = gsn_open_wks(type,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@mpUSStateLineThicknessF = 1.0 mpres@mpGeophysicalLineThicknessF = 1.0 ; get xlat,xlon from model for interpolation filename = "/data1/wrf/01_23_2016/GEFS/GEFS_00/wrfout_d04_2016-01-23_00:00:00" a = addfile(filename,"r") 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) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; times = wrf_user_getvar(a,"times",-1) ; get all times in the file plots = new ( 4, graphic ) do itime = 0,3 ; t = wrf_user_getvar(a,"T2",itime) ; T2 in Kelvin ; slp = wrf_user_getvar(a,"slp",itime) ; slp if (itime.eq.0)then date = "1259 UTC 23 12 2016" asc_file = "/data1/wrf/KOKX/KOKX_V06_20160123_125950_out.txt" ; asc_file = "jan26-28_2015.txt" end if if (itime.eq.1)then date = "1331 UTC 23 12 2016" asc_file = "/data1/wrf/KOKX/KOKX_V06_20160123_133101_out.txt" ; asc_file = "jan26-28_2015.txt" end if if (itime.eq.2)then date = "1356 UTC 23 12 2016" asc_file = "/data1/wrf/KOKX/KOKX_V06_20160123_135656_out.txt" ; asc_file = "jan26-28_2015.txt" end if if (itime.eq.3)then date = "1428 UTC 23 12 2016" asc_file = "/data1/wrf/KOKX/KOKX_V06_20160123_142808_out.txt" ; asc_file = "jan26-28_2015.txt" end if num_col = 3 num_obs = numAsciiRow(asc_file) z1 = asciiread(asc_file,(/num_obs,num_col/),"float") lat_o = z1(:,0) lon_o = z1(:,1) obs_o = z1(:,2) ; ; These long_name and units attributes are not necessary. I ; did this so I got better output from print statements below. ; obs_o@long_name = "observations of some sort" obs_o@units = "not sure" lat_o@long_name = "observed latitudes" lat_o@units = "degrees_north" lon_o@long_name = "observed longitudes" lon_o@units = "degrees_east" 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. nlat = dims2d(0) nlon = dims2d(1) xlat_m = xlat2d_m(:,nlon/2) xlon_m = xlon2d_m(nlat/2,:) printVarSummary(xlat_m) printVarSummary(xlon_m) printMinMax(xlat_m,0) printMinMax(xlon_m,0) ; rscan = (/1,0.5,0.25/) ; rscan = (/0.8,0.4,0.20/) ; rscan = (/0.25,0.125,0.05/) rscan = (/0.05/) grid = obj_anal_ic(lon_o,lat_o,obs_o,xlon_m,xlat_m,rscan,False) printVarSummary(grid) printMinMax(grid,0) ; Generate contours. res1 = res res1@cnFillOn = True ; res1@ContourParameters = (/ 0., 60., 5. /) res1@lbLabelBarOn = False ; Turn off individual label bars so we can res1@cnLevelSelectionMode = "ExplicitLevels" res1@cnFillOn = True ; res1@cnFillMode = "RasterFill" ; handles memory better, but can produce blocky contours res1@cnLinesOn = False res1@cnLevels = (/ -5, 0,5,10., \ 15, 20.,25.,30,35,40,45,50/) res1@cnFillColors = (/"White","White","LightSkyBlue1", \ "LightSkyBlue2", \ "LightSkyBlue3","SpringGreen1","SpringGreen4", \ "Yellow","Orange1","Orange3","Red","HotPink3","HotPink1","Violet"/) ; res1@lbAutoManage = False ; control label bar ; res1@pmLabelBarDisplayMode = "Always" ; turns on label bar contour = wrf_contour(a,wks,grid,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(itime) = wrf_map_overlays(a,wks,(/contour/),pltres,mpres) delete(contour) end do ; Panel the WRF plots. pnlres = True ; pnlres@txString = t@description + " (" + t@units + ")" pnlres@txString = "Observed Reflectivity (dbZ)" 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 gsn_panel(wks,(/plots/),(/2,2/),pnlres) end