; Example script to produce plots for a WRF real-data run, ; with the ARW coordinate dynamics option. ; Plot data on a cross section ; This script will plot data at a set angle through a specified point ; This script adds lon/lat info along X-axis load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl" external GET_IJ "./get_ij.so" ; external file required to find center location of cross section begin plot_rad1 = new ( 4, graphic ) plot_rad2 = new ( 4, graphic ) ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. a = addfile("wrfout_d02_2013-02-08_15:30:00.nc","r") ; my additional code -- set bounds xlat = wrf_user_getvar(a, "XLAT",0) xlon = wrf_user_getvar(a, "XLONG",0) ter = wrf_user_getvar(a, "HGT",0) dims2d=dimsizes(xlat) ; lats = (/ 39.00, 42.50 /) ; lons = (/ -75.0, -69.87 /) ; lats = (/ 39.00, 42.50 /) ; lons = (/ -76.0, -70.0/) lats = (/ 39.50, 42.00 /) lons = (/ -75.0, -71.0/) loc = wrf_user_ll_to_ij(a, lons, lats, True) ; 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 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("x_start = " + x_start) print("x_end = " + x_end) print("y_start = " + y_start) print("y_end = " + y_end) ; latitude longitude of Upton, NY xlat_zoom = xlat(y_start:y_end,x_start:x_end) xlon_zoom = xlon(y_start:y_end,x_start:x_end) ter_zoom = ter(y_start:y_end,x_start:x_end) dims_zoom = dimsizes(xlat_zoom) lat_beg = 40.8682 lon_beg = 72.8792 i_loc = 1 j_loc = 1 GET_IJ::get_ij(xlat_zoom,xlon_zoom,lat_beg,lon_beg,i_loc,j_loc,dims_zoom(0),dims_zoom(1)) i_point = i_loc-1 j_point = j_loc-1 print("i_point = " + i_point) print("j_point = " + j_point) print(xlat(j_point,i_point)) print(xlon(j_point,i_point)) ; We generate plots, but what kind do we prefer? ; type = "x11" type = "pdf" ; type = "ps" ; type = "ncgm" wks = gsn_open_wks(type,"plt_radar_temp_cross_sections_90_deg") ; wks = gsn_open_wks(type,"plt_radar_temp_cross_sections_00_deg") ; Set some basic resources ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTime = True FirstTimeMap = True times = wrf_user_getvar(a,"times",-1) ; get times in the file ntimes = dimsizes(times) ; number of times in the file 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) print("hour = " + hour) ; not needed with zoom ; mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file ; nd = dimsizes(mdims) ;--------------------------------------------------------------- ; do it = 0,ntimes-1,2 ; TIME LOOP ; do it = 20,20,1 do it = 17,22,1 it_val = it new_time_units = "hours since 1800-01-01 00:00" new_time_units = new_time_units new_time_array = cd_inv_calendar( year(it_val), month(it_val), day(it_val), \ hour(it_val), minute(it_val), 0, new_time_units,0) print("new_time_array = " + new_time_array) ; format = "" ;; defaults to "%H%M IST %d %c %Y" format = "%H%M UTC %d %c %Y" ; format = "%d %c %Y" new_time_array = new_time_array time_var=cd_string(new_time_array,format) print("time_var = " + time_var) print("Working on time: " + times(it) ) ; res@TimeLabel = times(it) ; Set Valid time to use on plots tc = wrf_user_getvar(a,"tc",it) ; T in C ; rh = wrf_user_getvar(a,"rh",it) ; relative humidity fft = wrf_user_getvar(a,"fft_zh",it) ; total reflectivity ff1 = wrf_user_getvar(a,"ff1_zh",it) ; total reflectivity ff2 = wrf_user_getvar(a,"ff2_zh",it) ; total reflectivity ff3 = wrf_user_getvar(a,"ff3_zh",it) ; total reflectivity ff4 = wrf_user_getvar(a,"ff4_zh",it) ; total reflectivity ff5 = wrf_user_getvar(a,"ff5_zh",it) ; total reflectivity ff6 = wrf_user_getvar(a,"ff6_zh",it) ; total reflectivity ff7 = wrf_user_getvar(a,"ff7_zh",it) ; total reflectivity z = wrf_user_getvar(a, "z",it) ; grid point height tc_zoom = tc(:,y_start:y_end,x_start:x_end) fft_zoom = fft(:,y_start:y_end,x_start:x_end) ff1_zoom = ff1(:,y_start:y_end,x_start:x_end) ff2_zoom = ff2(:,y_start:y_end,x_start:x_end) ff3_zoom = ff3(:,y_start:y_end,x_start:x_end) ff4_zoom = ff4(:,y_start:y_end,x_start:x_end) ff5_zoom = ff5(:,y_start:y_end,x_start:x_end) ff6_zoom = ff6(:,y_start:y_end,x_start:x_end) ff7_zoom = ff7(:,y_start:y_end,x_start:x_end) z_zoom = z(:,y_start:y_end,x_start:x_end) if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 6. ; We are only interested in the first 6km nz = floattoint(zmax + 1) end if ;--------------------------------------------------------------- ; do ip = 1, 3 ; we are doing 3 plots do ip = 2, 2 ; we are doing 3 plots ; do ip = 1, 1 ; we are doing 3 plots ; all with the pivot point (plane) in the center of the domain ; at angles 0, 45 and 90 ; ; | ; angle=0 is | ; | ; plane = new(2,float) ; note this part of the code assumes that the arrays have been created (x,y) ; plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is center of domain (x,y) plane = (/ i_point, j_point /) ; pivot point is Upton opts = False if(ip .eq. 1) then angle = 90. X_plane = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) X_desc = "longitude" end if if(ip .eq. 2) then angle = 0. X_plane = wrf_user_intrp2d(xlat_zoom,plane,angle,opts) X_desc = "latitude" end if if(ip .eq. 3) then ; angle = 45. angle = 135. X_plane = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) X_desc = "longitude" end if fft_plane = wrf_user_intrp3d(fft_zoom,z_zoom,"v",plane,angle,opts) ff1_plane = wrf_user_intrp3d(ff1_zoom,z_zoom,"v",plane,angle,opts) ff2_plane = wrf_user_intrp3d(ff2_zoom,z_zoom,"v",plane,angle,opts) ff3_plane = wrf_user_intrp3d(ff3_zoom,z_zoom,"v",plane,angle,opts) ff4_plane = wrf_user_intrp3d(ff4_zoom,z_zoom,"v",plane,angle,opts) ff5_plane = wrf_user_intrp3d(ff5_zoom,z_zoom,"v",plane,angle,opts) ff6_plane = wrf_user_intrp3d(ff6_zoom,z_zoom,"v",plane,angle,opts) ff7_plane = wrf_user_intrp3d(ff7_zoom,z_zoom,"v",plane,angle,opts) tc_plane = wrf_user_intrp3d(tc_zoom,z_zoom,"v",plane,angle,opts) printVarSummary(fft_plane) ; Find the index where 6km is - only need to do this once if ( FirstTime ) then zz = wrf_user_intrp3d(z_zoom,z,"v",plane,angle,opts) b = ind(zz(:,0) .gt. zmax*1000. ) zmax_pos = b(0) - 1 if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then zspan = b(0) - 1 else zspan = b(0) end if delete(zz) delete(b) FirstTime = False end if ; X-axis lables dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 nx = floattoint( (xmax-xmin)/2 + 3) ; ; xmin_new = nx-xmin ; xmax_new = nx+xmin ; print("xmin_new = " + xmin_new) ; print("xmax_new = " + xmax_new) ;--------------------------------------------------------------- ; Options for XY Plots opts_xy = True opts_xy@tiXAxisString = X_desc opts_xy@tiYAxisString = "Height (km)" opts_xy@cnMissingValPerimOn = True opts_xy@cnMissingValFillColor = 0 opts_xy@cnMissingValFillPattern = 11 opts_xy@tmXTOn = False opts_xy@tmYROn = False opts_xy@tmXBMode = "Explicit" opts_xy@tmXBValues = fspan(0,xspan,nx) ; Create tick marks opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,nx)) ; Create labels print("opts_xy@tmXBLabels" + opts_xy@tmXBLabels) opts_xy@tmXBLabelFontHeightF = 0.015 opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_xy@tmYLLabels = sprintf("%.1f",fspan(zmin,zmax,nz)) ; Create labels opts_xy@tiXAxisFontHeightF = 0.020 opts_xy@tiYAxisFontHeightF = 0.020 opts_xy@tmXBMajorLengthF = 0.02 opts_xy@tmYLMajorLengthF = 0.02 opts_xy@tmYLLabelFontHeightF = 0.015 ; opts_xy@PlotOrientation = tc_plane@Orientation ; Plotting options for RH opts_fft = opts_xy ; opts_fft = True opts_fft@cnFillOn = True ; opts_fft@ContourParameters = (/ -10., 60., 5. /) opts_fft@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels opts_fft@cnLevels = ispan(-9,57,3) ; set the contour levels cmap = read_colormap_file("NCV_rainbow2") cmap(0,:) = 1. ; Make first color white opts_fft@cnFillPalette = cmap opts_fft@pmLabelBarOrthogonalPosF = -0.1 ; opts_fft@lbLabelBarOn = False opts_fft@gsnDraw = False ; don't draw opts_fft@gsnFrame = False ; don't advance frame opts_fft@cnLinesOn = False ; this gives a warning (once) that it's not a valid setting, but it is. ; I think it happens because I reuse the opts for different variables. opts_fft@cnLineLabelsOn = False opts_fft@lbLabelBarOn = False opts_fft@cnInfoLabelOn = False opts_ff1 = opts_fft opts_ff2 = opts_fft opts_ff3 = opts_fft opts_ff4 = opts_fft opts_ff5 = opts_fft opts_ff6 = opts_fft opts_ff7 = opts_fft ; this gives an not valid at this time, but it is valid. Strange. opts_fft@gsnLeftString = "Total Radar Reflectivity" opts_fft@gsnRightString = "dBZ" opts_ff1@gsnLeftString = "Radar Reflectivity of Liquid Water" opts_ff2@gsnLeftString = "Radar Reflectivity of Plates" opts_ff3@gsnLeftString = "Radar Reflectivity of Columns" opts_ff4@gsnLeftString = "Radar Reflectivity of Dendrites" opts_ff5@gsnLeftString = "Radar Reflectivity of Snow" opts_ff6@gsnLeftString = "Radar Reflectivity of Graupel" opts_ff7@gsnLeftString = "Radar Reflectivity of Hail" ; opts_fft@cnFillColors = (/"White","White","White", \ ; "White","Chartreuse","Green", \ ; "Green3","Green4", \ ; "ForestGreen","PaleGreen4"/) ; res@lbLabelBarOn = False ; res@gsnDraw = False ; don't draw ; res@gsnFrame = False ; don't advance frame ; Plotting options for Temperature opts_tc = opts_xy opts_tc@cnInfoLabelZone = 1 opts_tc@cnInfoLabelSide = "Top" opts_tc@cnInfoLabelPerimOn = True opts_tc@cnInfoLabelOrthogonalPosF = -0.00005 ; opts_tc@ContourParameters = (/ 5. /) opts_tc@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels opts_tc@cnLevels = ispan(-16,8,2) ; set the contour levels opts_tc@gsnDraw = False ; don't draw opts_tc@gsnFrame = False ; don't advance frame opts_tc@cnInfoLabelOn = False opts_tc@gsnContourZeroLineThicknessF = 3.5 ; Get the contour info for the rh and temp plot_tc1 = gsn_contour(wks,tc_plane(0:zmax_pos,:),opts_tc) plot_tc2 = gsn_contour(wks,tc_plane(0:zmax_pos,:),opts_tc) plot_tc3 = gsn_contour(wks,tc_plane(0:zmax_pos,:),opts_tc) plot_tc4 = gsn_contour(wks,tc_plane(0:zmax_pos,:),opts_tc) plot_rad1(0) = gsn_csm_contour(wks,fft_plane(0:zmax_pos,:),opts_fft) plot_rad2(0) = gsn_csm_contour(wks,ff1_plane(0:zmax_pos,:),opts_ff1) plot_rad1(1) = gsn_csm_contour(wks,ff2_plane(0:zmax_pos,:),opts_ff2) plot_rad1(2) = gsn_csm_contour(wks,ff3_plane(0:zmax_pos,:),opts_ff3) plot_rad2(1) = gsn_csm_contour(wks,ff5_plane(0:zmax_pos,:),opts_ff5) plot_rad2(2)= gsn_csm_contour(wks,ff6_plane(0:zmax_pos,:),opts_ff6) plot_rad1(3) = gsn_csm_contour(wks,ff4_plane(0:zmax_pos,:),opts_ff4) plot_rad2(3) = gsn_csm_contour(wks,ff7_plane(0:zmax_pos,:),opts_ff7) ;--------------------------------------------------------------- ; MAKE PLOTS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False ; Panel the WRF plots. Type = "QC + QR" pnlres = True ; pnlres@txString = t@description + " (" + t@units + ")" ; pnlres@txString = Type + " at " + radar_deg + " deg elevation (g/cm^3)" ; pnlres@txString = "Vertical Cross Sections" pnlres@txString = time_var 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 overlay(plot_rad1(0),plot_tc1) ; NhlRemoveOverlay(plot_rad1(0),plot_tc,False) overlay(plot_rad1(1),plot_tc2) ; NhlRemoveOverlay(plot_rad1(1),plot_tc,False) overlay(plot_rad1(2),plot_tc3) ; NhlRemoveOverlay(plot_rad1(2),plot_tc,False) overlay(plot_rad1(3),plot_tc4) gsn_panel(wks,(/plot_rad1/),(/2,2/),pnlres) NhlRemoveOverlay(plot_rad1(0),plot_tc1,False) NhlRemoveOverlay(plot_rad1(1),plot_tc2,False) NhlRemoveOverlay(plot_rad1(2),plot_tc3,False) NhlRemoveOverlay(plot_rad1(3),plot_tc4,False) overlay(plot_rad2(0),plot_tc1) overlay(plot_rad2(1),plot_tc2) overlay(plot_rad2(2),plot_tc3) overlay(plot_rad2(3),plot_tc4) gsn_panel(wks,(/plot_rad2/),(/2,2/),pnlres) delete(opts_xy) delete(opts_tc) delete(opts_fft) delete(tc_plane) delete(fft_plane) delete(ff1_plane) delete(ff2_plane) delete(ff3_plane) delete(ff4_plane) delete(ff5_plane) delete(ff6_plane) delete(ff7_plane) delete(X_plane) end do ; make next cross section ; Delete options and fields, so we don't have carry over end do ; END OF TIME LOOP end