;---------------------------------------------------------------------- ; This script shows how to create a basic plot of temperature (degC) ; calculated from a WRF file. ; ; This script is meant to be compared to one that uses wrf_xxxx ; functions to do the plotting. external INTERP1D "./interp1d.so" external GET_IJ "./get_ij.so" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl" ;---------------------------------------------------------------------- undef("set_res1_list") function set_res1_list() begin res = True res@pmTickMarkDisplayMode = "Always" res@mpFillOn = False ; turn off map fill res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries res@mpLimbLineColor = "Black" res@mpPerimLineColor = "Black" res@mpPerimLineThicknessF = 10 res@mpUSStateLineThicknessF = 10.0 res@mpGeophysicalLineThicknessF = 10.0 res@mpGridAndLimbOn = False res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@lbLabelBarOn = False res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels ; res@cnLevels = ispan(-9,57,3) ; set the contour levels cmap = read_colormap_file("NCV_rainbow2") cmap(0,:) = 1. ; Make first color white res@cnFillPalette = cmap res@cnLevels = (/0.05,0.10,0.25,0.5,0.75,1.0,1.25,1.50,1.75,2.0,\ 2.5,3.0,3.5,4.0,4.5,5.0/) res@cnFillColors = (/0,105,95,85,75,115,122,129,157,169,176,183,190,197,204,211, \ 218,225,239,247,252/) ; res@cnLevels = (/0.05,0.10,0.25,0.5,0.75,1.0,1.25,1.50,1.75,2.0,\ ; 2.5,3.0,3.5,4.0,4.5,5.0/) ; If you want the first color to be transparent rather than white, then use: ; cmap(0,3) = 0,0 ; Make first color transparent res@cnLinesOn = False ;--- These are required for gsn_csm_contour_map plotting ;res = wrf_map_resources(f,res) res@tfDoNDCOverlay = True res@gsnAddCyclic = False return(res) end undef("set_res2_list") function set_res2_list() begin res = True res@pmTickMarkDisplayMode = "Always" res@mpFillOn = False ; turn off map fill res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries res@mpLimbLineColor = "Black" res@mpPerimLineColor = "Black" res@mpPerimLineThicknessF = 10 res@mpUSStateLineThicknessF = 10.0 res@mpGeophysicalLineThicknessF = 10.0 res@mpGridAndLimbOn = False res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@lbLabelBarOn = False res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels cmap = read_colormap_file("NCV_rainbow2") cmap(0,:) = 1. ; Make first color white res@cnFillPalette = cmap ; res@cnFillPalette="NCV_rainbow2" ; res@cnLevels = (/0.10,0.25,0.5,1.0,2.5,\ ; 5.0,10,15,20,25,30,35,40,45,50/) ; res@cnLevels = (/0.10,0.25,0.5,1.0,1.5,\ ; 2.0,2.5,3,3.5,4.,4.5,5,6,7/) ; res@cnFillColors = (/0,90,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204,211, \ ; 218,225,232,239,247,252/) res@cnLevels = (/0.05,0.10,0.25,0.5,1.0,1.5,\ 2.0,2.5,3,3.5,4.,4.5,5,6,7/) res@cnFillColors = (/0,90,80,70,63,108,115,122,129,157,161,169,176,183,190,197,204,211, \ 218,225,232,239,247,252/) return(res) end undef("set_res3_list") function set_res3_list() begin res = True res@pmTickMarkDisplayMode = "Always" res@mpFillOn = False ; turn off map fill res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries res@mpLimbLineColor = "Black" res@mpPerimLineColor = "Black" res@mpPerimLineThicknessF = 10 res@mpUSStateLineThicknessF = 10.0 res@mpGeophysicalLineThicknessF = 10.0 res@mpGridAndLimbOn = False res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@lbLabelBarOn = False res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame ; res@cnFillPalette = "MPL_gist_ncar" cmap = read_colormap_file("MPL_gist_ncar") cmap(0,:) = 1. ; Make first color white res@cnFillPalette = cmap res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels res@cnLevels = (/ 0.70,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.80,0.81,\ 0.82,0.83,0.84,0.85,0.86,0.87,0.88,\ 0.89,0.90,0.91,0.92,0.93,0.94,0.95,0.96,\ 0.97,0.98,0.99,1.0/) res@cnFillColors = (/ 0,17,21,25,29,33,37,41,\ 45,49,53,57,61,65,69,73,75,77,\ 81 ,83 ,85,\ 87,89,91,93,95,126,122,117,115,113,111,110/) res@cnLinesOn = False ;--- These are required for gsn_csm_contour_map plotting ;res = wrf_map_resources(f,res) res@tfDoNDCOverlay = True res@gsnAddCyclic = False return(res) end undef("set_res4_list") function set_res4_list() begin res = True res@pmTickMarkDisplayMode = "Always" res@mpFillOn = False ; turn off map fill res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries res@mpLimbLineColor = "Black" res@mpPerimLineColor = "Black" res@mpPerimLineThicknessF = 10 res@mpUSStateLineThicknessF = 10.0 res@mpGeophysicalLineThicknessF = 10.0 res@mpGridAndLimbOn = False res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@lbLabelBarOn = False res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels cmap = read_colormap_file("NCV_rainbow2") cmap(0,:) = 1. ; Make first color white res@cnFillPalette = cmap ; res@cnLevels = (/ 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.2,1.8/) res@cnLevels = (/-5,-2,-1.0,-0.5,-0.25,0,0.01,0.02,0.06,0.1,0.2,0.3,0.4,0.5,0.6,\ 0.7,0.8,0.9,1.0,2.0,3.0,5.0/) res@cnFillColors = (/0,90,80,70,63,55,108,115,122,129,157,161,169,176,183,190,197,204,211, \ 218,225,232,239,247,251,252,253/) return(res) end begin radar_deg = 0.50 ;radar_deg = 1.50 del_x = 1.3*1000. plots_dbz = new ( 6, graphic ) plots_zdr = new ( 6, graphic ) plots_crs = new ( 6, graphic ) plots_dp = new ( 6, graphic ) plots_ha = new ( 6, graphic ) plot_summary = new ( 8, graphic ) filename = "wrfout_d02_2013-02-08_15:30:00" f = addfile(filename,"r") print(f) dx = f->RDX(0) dx = 1./dx print("dx = " + dx) dx = dx*0.001 ;print(f) times = wrf_user_list_times(f) ; get times in the file lats = (/ 39.00, 42.50 /) lons = (/ -75.0, -71.5/) ; bounding box for February 08 2013 ; radar = wrf_user_getvar(f,"radar",0) ; Temperature (C) xlat2d = wrf_user_getvar(f,"XLAT",0) ; Temperature (C) xlon2d = wrf_user_getvar(f,"XLONG",0) ; Temperature (C) dims2d = dimsizes(xlat2d) ; latitude longitude of Upton, NY lat_beg = 40.8682 lon_beg = 72.8792 i_loc = 1 j_loc = 1 GET_IJ::get_ij(xlat2d,xlon2d,lat_beg,lon_beg,i_loc,j_loc,dims2d(0),dims2d(1)) i_point = i_loc-1 j_point = j_loc-1 itt = 0 wtype = "png" ; wtype = "pdf" wtype@wkWidth = 2500 ; produces a better looking PNG wtype@wkHeight = 2500 ; wks = gsn_open_wks(wtype,out_name) wks = gsn_open_wks(wtype,"wrf_panels_plot_surf_gsn") ; Create various resource lists res1 = set_res1_list() res2 = set_res1_list() res3 = res1 res4 = res1 res5 = res1 ;do it = 17,22,1 do it = 16,21,1 ;do it = 17,17,1 yr = str_get_cols(times, 0, 3) year = stringtoint(yr) mnth = str_get_cols(times, 5, 6) month = stringtoint(mnth) dy = str_get_cols(times, 8, 9) day = stringtoint(dy) hr = str_get_cols(times, 11, 12) hour = stringtoint(hr) mn = str_get_cols(times, 14, 15) minute = stringtoint(mn) new_time_units = "hours since 1800-01-01 00:00" new_time_array = cd_inv_calendar( year(it), month(it), day(it), \ hour(it), minute(it), 0, new_time_units,0) print("new_time_array = " + new_time_array) ; format = "" ;; defaults to "%H%M EST %d %c %Y" format = "" ;; defaults to "%H%M UTC %d %c %Y" format = "%H%M UTC %d %c %Y" ; new_time_array = new_time_array + 2 new_time_array = new_time_array + 0 time_var=cd_string(new_time_array,format) print("time_var = " + time_var) radar = wrf_user_getvar(f,"SNOWH",it) zdr = wrf_user_getvar(f,"RAINNC",it) crs = wrf_user_getvar(f,"SNOWNC",it) kdp = wrf_user_getvar(f,"GRAUPELNC",it) hail = wrf_user_getvar(f,"HAILNC",it) radar_o = wrf_user_getvar(f,"SNOWH",it-1) zdr_o = wrf_user_getvar(f,"RAINNC",it-1) crs_o = wrf_user_getvar(f,"SNOWNC",it-1) kdp_o = wrf_user_getvar(f,"GRAUPELNC",it-1) hail_o = wrf_user_getvar(f,"HAILNC",it-1) radar = 100.*(radar-radar_o) zdr = zdr-zdr_o printVarSummary(zdr) delete( [/ zdr@description, zdr@units/] ) ; if present delete( [/ radar@description, radar@units/] ) ; if present delete( [/ crs@description, crs@units/] ) ; if present delete( [/ kdp@description, kdp@units/] ) ; if present delete( [/ hail@description, hail@units/] ) ; if present crs = crs-crs_o kdp = kdp-kdp_o hail = hail-hail_o ;printVarSummary(radar) nl = 5 ; level index ;plot = gsn_csm_contour_map(wks,radar(n1:,:),res1) dims3d = dimsizes(radar) print(dims3d) ; set up times ;res1@tiMainString = time_var ;res2@tiMainString = time_var ;res3@tiMainString = time_var ;res4@tiMainString = time_var ;res5@tiMainString = time_var res1@gsnCenterString = time_var res2@gsnCenterString = time_var res3@gsnCenterString = time_var res4@gsnCenterString = time_var res5@gsnCenterString = time_var ; 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(f, 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 print("y_start = " + y_start) print("x_start = " + x_start) print("y_end = " + y_end) print("x_end = " + x_end) res1@mpMinLatF = xlat2d(y_start,x_start) res1@mpMaxLatF = xlat2d(y_end,x_end) res1@mpMinLonF = xlon2d(y_start,x_start) res1@mpMaxLonF = xlon2d(y_end,x_end) res2@mpMinLatF = xlat2d(y_start,x_start) res2@mpMaxLatF = xlat2d(y_end,x_end) res2@mpMinLonF = xlon2d(y_start,x_start) res2@mpMaxLonF = xlon2d(y_end,x_end) res3@mpMinLatF = xlat2d(y_start,x_start) res3@mpMaxLatF = xlat2d(y_end,x_end) res3@mpMinLonF = xlon2d(y_start,x_start) res3@mpMaxLonF = xlon2d(y_end,x_end) res4@mpMinLatF = xlat2d(y_start,x_start) res4@mpMaxLatF = xlat2d(y_end,x_end) res4@mpMinLonF = xlon2d(y_start,x_start) res4@mpMaxLonF = xlon2d(y_end,x_end) res5@mpMinLatF = xlat2d(y_start,x_start) res5@mpMaxLatF = xlat2d(y_end,x_end) res5@mpMinLonF = xlon2d(y_start,x_start) res5@mpMaxLonF = xlon2d(y_end,x_end) res6 = res3 res7 = res4 res8 = res5 res9 = res2 res6@gsnTickMarksOn = False res7@gsnTickMarksOn = False res8@gsnTickMarksOn = False res9@gsnTickMarksOn = False res6@gsnCenterString = "Snow" res7@gsnCenterString = "Graupel" res8@gsnCenterString = "Hail" res9@gsnCenterString = "Total" ; radar_xy = radar zdr_xy = zdr crs_xy = crs dp_xy = kdp ha_xy = hail printMinMax(dp_xy,False) ; calculate differential phase from model KDP ; print("calculating differential phase") ; factor = dx*2 ; dp = factor*kdp print("calculated differential phase") ; radar_plane = radar_xy(y_start:y_end,x_start:x_end) zdr_plane = zdr_xy(y_start:y_end,x_start:x_end) crs_plane = crs_xy(y_start:y_end,x_start:x_end) dp_plane = dp_xy(y_start:y_end,x_start:x_end) ha_plane = ha_xy(y_start:y_end,x_start:x_end) xlat2d_sub = xlat2d(y_start:y_end,x_start:x_end) xlon2d_sub = xlon2d(y_start:y_end,x_start:x_end) radar_plane@lat2d = xlat2d_sub ; Needed for plotting radar_plane@lon2d = xlon2d_sub zdr_plane@lat2d = xlat2d_sub ; Needed for plotting zdr_plane@lon2d = xlon2d_sub crs_plane@lat2d = xlat2d_sub ; Needed for plotting crs_plane@lon2d = xlon2d_sub dp_plane@lat2d = xlat2d_sub ; Needed for plotting dp_plane@lon2d = xlon2d_sub ha_plane@lat2d = xlat2d_sub ; Needed for plotting ha_plane@lon2d = xlon2d_sub plot_radar = gsn_csm_contour_map(wks,radar_plane,res1) plots_dbz(itt) = plot_radar plot_zdr = gsn_csm_contour_map(wks,zdr_plane,res2) plots_zdr(itt) = plot_zdr plot_crs = gsn_csm_contour_map(wks,crs_plane,res3) plots_crs(itt) = plot_crs plot_dp = gsn_csm_contour_map(wks,dp_plane,res4) plots_dp(itt) = plot_dp plot_ha = gsn_csm_contour_map(wks,ha_plane,res5) plots_ha(itt) = plot_ha if (itt.eq.0)then plot_crs = gsn_csm_contour_map(wks,zdr_plane,res9) plot_summary(0) = plot_crs plot_crs = gsn_csm_contour_map(wks,crs_plane,res6) plot_summary(1) = plot_crs plot_dp = gsn_csm_contour_map(wks,dp_plane,res7) plot_summary(2) = plot_dp plot_ha = gsn_csm_contour_map(wks,ha_plane,res8) plot_summary(3) = plot_ha end if if (itt.eq.4)then plot_crs = gsn_csm_contour_map(wks,zdr_plane,res9) plot_summary(4) = plot_crs plot_crs = gsn_csm_contour_map(wks,crs_plane,res6) plot_summary(5) = plot_crs plot_dp = gsn_csm_contour_map(wks,dp_plane,res7) plot_summary(6) = plot_dp plot_ha = gsn_csm_contour_map(wks,ha_plane,res8) plot_summary(7) = plot_ha end if itt = itt+1 end do ; Panel the WRF plots. Type = "Snow Accumulation" pnlres = True ; pnlres@txString = t@description + " (" + t@units + ")" pnlres@txString = Type + " (cm)" pnlres@gsnPanelYWhiteSpacePercent = 2 ; 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 = "cm" ; 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) Type = "Accumulated Liq. Equiv. Prec" pnlres@txString = Type + " (mm)" pnlres@gsnPanelYWhiteSpacePercent = 13 ; Add white space b/w plots. ; pnlres@lbTitleString = "dBZ" ; title string pnlres@lbTitleString = "mm" ; title string gsn_panel(wks,(/plots_zdr/),(/2,3/),pnlres) Type = "Accumulated Liq. Equiv. Snow" pnlres@txString = Type + " (mm)" pnlres@gsnPanelYWhiteSpacePercent = 13 ; Add white space b/w plots. ; pnlres@lbTitleString = " " ; title string gsn_panel(wks,(/plots_crs/),(/2,3/),pnlres) Type = "Graupel Liq. Equiv. Accumulation" pnlres@txString = Type + " (mm)" gsn_panel(wks,(/plots_dp/),(/2,3/),pnlres) Type = "Hail Liq. Equiv. Accumulation" pnlres@txString = Type + " (mm)" gsn_panel(wks,(/plots_ha/),(/2,3/),pnlres) pnlres@txString = "30 Minute Accumulations" + " (mm)" gsn_panel(wks,(/plot_summary/),(/2,4/),pnlres) end