; 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" external GET_IJ "./get_ij.so" ; external file required to find center location of cross section begin ; ; 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_maps_w_angle_line") ; 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 ; 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 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 z = wrf_user_getvar(a, "z",it) ; grid point height 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 ; 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 ; 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 + 1) ;--------------------------------------------------------------- ; MAKE PLOTS if (FirstTimeMap) then lat_plane = wrf_user_intrp2d(xlat_zoom,plane,angle,opts) lon_plane = wrf_user_intrp2d(xlon_zoom,plane,angle,opts) mpres = True mpres@ZoomIn = True mpres@Xstart = x_start mpres@Ystart = y_start mpres@Xend = x_end mpres@Yend = y_end pltres = True pltres@FramePlot = False optsM = res optsM@NoHeaderFooter = True optsM@cnFillOn = True optsM@lbTitleOn = False contour = wrf_contour(a,wks,ter_zoom,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 delete(X_plane) ; Delete options and fields, so we don't have carry over end do ; make next cross section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False end do ; END OF TIME LOOP end