; 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 "/usr/lib/ncarg/nclscripts/csm/gsn_code.ncl" ;load "/usr/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "WRFUserARW.ncl" ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. dir = "./" fili = "wrfout_d03_2015-09-11_223000.nc" filin= dir+fili a = addfile(filin,"r") ; We generate plots, but what kind do we prefer? ;type = "x11" ; type = "pdf" ; type = "ps" ; type = "ncgm" type = "png" pltres = True ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file nd = dimsizes(mdims) xlat = wrf_user_getvar(a, "XLAT",0) xlon = wrf_user_getvar(a, "XLONG",0) ter = wrf_user_getvar(a, "HGT",0) ;--------------------------------------------------------------- ;;do it = 0,ntimes-1,1 ; TIME LOOP it=0 print("Working on time: " + times(it) ) wks = gsn_open_wks(type,"nice"+times(it)) rh = wrf_user_getvar(a,"rh",it) ; relative humidity z = wrf_user_getvar(a, "z",it) ; grid point height ;-------------------------------------------- ; Generating the cross section Vertical coordinate values ; (This is the same as WRFUserARW.ncl) z_max = max(z) z_min = 0. dz = 0.0025 * z_max nlevels = tointeger( z_max/dz ) z_var2d = new( (/nlevels/), typeof(z)) z_var2d(0) = z_min do i=1, nlevels-1 z_var2d(i) = z_var2d(0)+i*dz end do ;----------------------------------------------------------------- if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 10. ; We are only interested in the first 6km nz = floattoint(zmax + 1) FirstTime = False end if ;--------------------------------------------------------------- minlat = 6.502222 maxlat = 6.433889 minlon = -75.826667 maxlon = -75.650278 opt = True locA = wrf_user_ll_to_ij(a,minlon,minlat,opt) locB = wrf_user_ll_to_ij(a,maxlon,maxlat,opt) lonA = locA(0) - 1 latA = locA(1) - 1 lonB = locB(0) - 1 latB = locB(1) - 1 plane = new(4,float) angle = 0. plane = (/ lonA,latA , lonB,latB /) ;--- start x;y & end x;y point opts = opt X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "Longitude" ter_plane = wrf_user_intrp2d(ter,plane,angle,opts) print("Max terrain height in plot " + max(ter_plane)) rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts) rh_plane&Vertical = z_var2d ;printVarSummary(rh_plane) rh_plane2 = rh_plane rh_plane2&Vertical = z_var2d ;printVarSummary(rh_plane2) cross_dims = dimsizes(rh_plane2) rank = dimsizes(cross_dims) ;printVarSummary(rh_plane2) iz_start = cross_dims(0) / 2 print(iz_start) do iz = iz_start,0,1 do ix = 0,cross_dims(rank-1)-1 if ( ismissing(rh_plane2(iz,ix)) ) then rh_plane2(iz,ix) = rh_plane2(iz+1,ix) end if end do end do ; X-axis labels dimsX = dimsizes(X_plane) xmin = X_plane(0) xmax = X_plane(dimsX(0)-1) xspan = dimsX(0)-1 print(xspan) nx=4 ;nx = floattoint( (xmax-xmin)/2 + 1) ;--------------------------------------------------------------- ; Options for XY Plots opts_xy = True opts_xy@gsnFrame = False opts_xy@gsnDraw = False opts_xy@tiXAxisString = X_desc opts_xy@tiYAxisString = "Height (km)" opts_xy@cnMissingValPerimOn = True 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 opts_xy@tmXBLabelFontHeightF = 0.015 opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = fspan(zmin*1000,zmax*1000,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@trYMaxF = zmax*1000 opts_xy@trYMinF = 0 ; Plotting options for RH opts_rh = opts_xy opts_rh@cnLevelSelectionMode = "ManualLevels" opts_rh@cnMinLevelValF = 70. ;opts_rh@ContourParameters = (/ 70., 98., 2. /) opts_rh@cnMaxLevelValF = 98. opts_rh@cnLevelSpacingF = 2. opts_rh@cnFillOn = True opts_rh@cnFillPalette = "WhiteGreen" ; opts_rh@pmLabelBarOrthogonalPosF = 0.1 opts_rh@pmLabelBarHeightF = .075 opts_rh@lbTitlePosition = "Top" opts_rh@lbTitleString = "Relative Humidity (%)" opts_rh@lbTitleFontHeightF = 0.01 opts_rh@lbTitleOffsetF = -.03 ; MAKE PLOTS if (FirstTimeMap) then lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts) lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts) mpres = True pltres = True pltres@FramePlot = False optsM = True optsM@NoHeaderFooter = True optsM@cnFillOn = True optsM@cnFillPalette = "OceanLakeLandSnow" optsM@lbTitleOn = False mpres@mpNationalLineColor = "black" mpres@mpUSStateLineColor = "black" mpres@mpUSStateLineThicknessF = 3.0 mpres@mpNationalLineThicknessF = 3.0 mpres@mpGeophysicalLineThicknessF = 3.0 contour = wrf_contour(a,wks,ter,optsM) plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres) lnres = True lnres@gsLineThicknessF = 5.0 lnres@gsLineColor = "Red" do ii = 0,dimsX(0)-2 gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) end do frame(wks) delete(lon_plane) delete(lat_plane) pltres@FramePlot = True end if ;Contour terrain cross section opts_ter = True opts_ter@gsnYRefLine = 0. opts_ter@gsnAboveYRefLineColor = "black" opts_ter@gsnDraw = False opts_ter@gsnFrame = False opts_ter@tmYLOn = False opts_ter@tmYROn = False opts_ter@tiYAxisOn = False opts_ter@trYMaxF = zmax*1000 opts_ter@trYMinF = 0 ; Make the terrain contours contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter) ; Get the contour info for the rh and temp contour_rh = gsn_csm_contour(wks,rh_plane,opts_rh) rh_plane2&Vertical = z_var2d ;printVarSummary(rh_plane) ;printVarSummary(rh_plane2) contour_rh2 = gsn_csm_contour(wks,rh_plane2,opts_rh) ;--------------------------------------------------------------- ; Choose contour_rh if you don't mind the white at the boundary ; with the mountains. Choose contour_rh2 if you want the white ; filled in. ;draw(contour_rh) draw(contour_rh2) draw(contour_ter) frame(wks) ; Delete options and fields, so we don't have carry over delete(opts_xy) delete(opts_rh) delete(rh_plane) delete(rh_plane2) delete(X_plane) delete(ter_plane) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False ;end do ; END OF TIME LOOP