; 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/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" begin ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. dFPath = "/project/k14/basit/wrf-sims/wrf_170/" dFName = "auxhist7_d02_2013-06-04_00:00:00.nc" cFName = "auxhist7_d02_2013-06-04_00:00:00.nc" a = addfile (dFPath+dFName,"r") lat = 21.70 lon = 40.25 stHr = 37 endHr = 37 ; We generate plots, but what kind do we prefer? ; type = "x11" type = "png" wks = gsn_open_wks(type,"vcs_test") gsn_define_colormap(wks,"testcmap") ; Set some basic resources res = True res@MainTitle = "REAL-TIME WRF" res@Footer = False 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) ;--------------------------------------------------------------- opt_coord = True opt_coord@returnInt = False loc_vcs = wrf_user_ll_to_ij(a, lon, lat,opt_coord) i_coord = loc_vcs(0) j_coord = loc_vcs(1) do it = 37,37,1 ;ntimes-1,1 ; TIME LOOP print("Working on time: " + times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots ua = wrf_user_getvar(a,"ua",it) ; T in C th = wrf_user_getvar(a,"theta",it) ; relative humidity p = wrf_user_getvar(a,"pressure",it) z = wrf_user_getvar(a,"z",it) ; grid point height if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 1. ; We are only interested in the first 6km nz = floattoint(zmax + 1) end if ;--------------------------------------------------------------- do ip = 1, 3 ; 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) ; plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is center of domain (x,y) plane = (/i_coord, j_coord/) opts = False if(ip .eq. 1) then angle = 90. X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "longitude" end if if(ip .eq. 2) then angle = 0. X_plane = wrf_user_intrp2d(xlat,plane,angle,opts) X_desc = "latitude" end if if(ip .eq. 3) then angle = 45. X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "longitude" end if ua_plane = wrf_user_intrp3d(ua,z,"v",plane,angle,opts) th_plane = wrf_user_intrp3d(th,z,"v",plane,angle,opts) p_plane = wrf_user_intrp3d(p,z,"v",plane,angle,opts) ; Find the index where 6km is - only need to do this once if ( FirstTime ) then zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts) printVarSummary(zz) print(zz(:,0)) b = ind(zz(:,0) .gt. zmax*1000. ) zmax_pos = b(0) - 1 print("zmax_pos"+zmax_pos) 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 + 1) print("xspan : "+xspan) print("nx : "+nx) print("zspan : "+zspan) print("nz : "+nz) ;--------------------------------------------------------------- ; Options for XY Plots opts_xy = res 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,10) ; Create tick marks opts_xy@tmXBLabels = sprintf("%.1f",fspan(xmin,xmax,10)) ; Create labels 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 = ua_plane@Orientation ; Plotting options for ua opts_ua = opts_xy opts_ua@ContourParameters = (/-20., 20., 1. /) opts_ua@lbTitleOn = True opts_ua@lbTitleString = "" ;"W-E wind component (~N~m s~S-1~N~)" opts_ua@pmLabelBarOrthogonalPosF = -0.1 opts_ua@cnFillOn = True ; Plotting options for theta opts_th = opts_xy opts_th@cnInfoLabelZone = 1 opts_th@cnInfoLabelSide = "Top" opts_th@cnInfoLabelPerimOn = True opts_th@cnInfoLabelOrthogonalPosF = -0.00005 opts_th@cnLevelSpacingF = 1. opts_th@cnLineColor ="grey29" opts_th@cnLineThicknessF = 2.0 ; Plotting options for p opts_p = opts_xy opts_p@cnInfoLabelZone = 1 opts_p@cnInfoLabelSide = "Top" opts_p@cnInfoLabelPerimOn = True opts_p@cnInfoLabelOrthogonalPosF = -0.00005 opts_p@cnLevelSpacingF = 10. opts_p@cnLineColor ="brown" opts_p@cnLineThicknessF = 2.0 ; Get the contour info for the rh and temp contour_th = wrf_contour(a,wks,th_plane(0:zmax_pos,:),opts_th) contour_ua = wrf_contour(a,wks,ua_plane(0:zmax_pos,:),opts_ua) contour_p = wrf_contour(a,wks,p_plane(0:zmax_pos,:),opts_p) ;--------------------------------------------------------------- ; 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 = res optsM@NoHeaderFooter = True optsM@cnFillOn = True optsM@lbTitleOn = False contour = wrf_contour(a,wks,ter,optsM) plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres) lnres = True lnres@gsLineThicknessF = 3.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 plot = wrf_overlays(a,wks,(/contour_ua,contour_th,contour_p/),pltres) ; plot x-section ; Delete options and fields, so we don't have carry over delete(opts_xy) delete(opts_th) delete(opts_ua) delete(opts_p) delete(th_plane) delete(ua_plane) delete(p_plane) delete(X_plane) end do ; make next cross section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False end do ; END OF TIME LOOP end