; 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" 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") ; We generate plots, but what kind do we prefer? ; type = "x11" type = "pdf" ; type = "ps" ; type = "ncgm" wks = gsn_open_wks(type,"plt_zh_CrossSection") ; 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) dims2d=dimsizes(xlat) lats = (/ 39.00, 42.00 /) lons = (/ -75.0, -69.87 /) 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 ; latitude longitude of Upton, NY lat_beg = 40.8682 lon_beg = 72.8792 i_loc = 1 j_loc = 1 GET_IJ::get_ij(xlat,xlon,lat_beg,lon_beg,i_loc,j_loc,dims2d(0),dims2d(1)) i_point = i_loc-1 j_point = j_loc-1 print("i_point = " + i_point) print("j_point = " + j_point) ;--------------------------------------------------------------- ; do it = 0,ntimes-1,2 ; TIME LOOP do it = ntimes/2,ntimes/2,1 ; TIME LOOP 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 zh = wrf_user_getvar(a,"fft_zh",it) ; relative humidity z = wrf_user_getvar(a, "z",it) ; grid point height height = 400. ; 1km zh_plane_400 = wrf_user_intrp3d(zh,z,"h",height,0.,False) height = 600. ; 1km zh_plane_600 = wrf_user_intrp3d(zh,z,"h",height,0.,False) height = 800. ; 1km zh_plane_800 = wrf_user_intrp3d(zh,z,"h",height,0.,False) height = 1000. ; 1km zh_plane_1000 = wrf_user_intrp3d(zh,z,"h",height,0.,False) height = 1250. ; 1km zh_plane_1250 = wrf_user_intrp3d(zh,z,"h",height,0.,False) height = 1500. ; 1km zh_plane_1500 = wrf_user_intrp3d(zh,z,"h",height,0.,False) if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = 10. ; 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_point, j_point /) ; pivot point is center of domain (x,y) 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 = 135. X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) X_desc = "longitude" end if rh_plane = wrf_user_intrp3d(zh,z,"v",plane,angle,opts) tc_plane = wrf_user_intrp3d(tc,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) 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) ;--------------------------------------------------------------- ; 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,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(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_rh = opts_xy opts_rh@cnFillOn = True opts_rh@ContourParameters = (/ -5., 60., 5 /) opts_rh@cnFillColors = (/"DarkOrchid4","DarkOrchid3","DarkOrchid1",\ "DeepSkyBlue4","DeepSkyBlue3","DeepSkyBlue2","DeepSkyBlue1", \ "PaleGreen1","PaleGreen2", \ "LightYellow2","LightYellow1", \ "Gold1","Gold2","Gold3","Gold4", \ "Red4","Red3","Red2","Red1", \ "HotPink4","HotPink3","HotPink2","HotPink1","Maroon4"/) ; 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. /) ; Get the contour info for the rh and temp printVarSummary(tc_plane) var_zoom = tc_plane(:,y_start:y_end) contour_tc = wrf_contour(a,wks,var_zoom(0:zmax_pos,:),opts_tc) var_zoom = rh_plane(:,y_start:y_end) contour_rh = wrf_contour(a,wks,var_zoom(0:zmax_pos,:),opts_rh) ;--------------------------------------------------------------- ; 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 mpres = True mpres = True mpres@mpFillOn = False mpres@mpDataBaseVersion = "MediumRes" ; mpres@mpDataBaseVersion = "HighRes" mpres@mpDataSetName = "Earth..2" mpres@mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres@mpGeophysicalLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpUSStateLineColor = "Black" mpres@mpNationalLineColor = "Black" mpres@mpGridLineColor = 0 mpres@mpLimbLineColor = "Black" mpres@mpPerimLineColor = "Black" mpres@mpUSStateLineThicknessF = 1.0 mpres@mpGeophysicalLineThicknessF = 1.0 mpres@mpGridAndLimbOn = False 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 optsM@ContourParameters = (/ 00., 60., 5. /) optsM@ContourParameters = (/ -5., 60., 5 /) optsM@cnFillOn = True optsM@cnFillColors = (/"DarkOrchid4","DarkOrchid3","DarkOrchid1",\ "DeepSkyBlue4","DeepSkyBlue3","DeepSkyBlue2","DeepSkyBlue1", \ "PaleGreen1","PaleGreen2", \ "LightYellow2","LightYellow1", \ "Gold1","Gold2","Gold3","Gold4", \ "Red4","Red3","Red2","Red1", \ "HotPink4","HotPink3","HotPink2","HotPink1","Maroon4"/) printVarSummary(zh) delete(var_zoom) var_zoom_400 = zh_plane_400(y_start:y_end,x_start:x_end) contour_400 = wrf_contour(a,wks,var_zoom_400,optsM) var_zoom_600 = zh_plane_600(y_start:y_end,x_start:x_end) contour_600 = wrf_contour(a,wks,var_zoom_600,optsM) var_zoom_800 = zh_plane_800(y_start:y_end,x_start:x_end) contour_800 = wrf_contour(a,wks,var_zoom_800,optsM) var_zoom_1000 = zh_plane_1000(y_start:y_end,x_start:x_end) contour_1000 = wrf_contour(a,wks,var_zoom_1000,optsM) var_zoom_1250 = zh_plane_1250(y_start:y_end,x_start:x_end) contour_1250 = wrf_contour(a,wks,var_zoom_1250,optsM) var_zoom_1500 = zh_plane_1500(y_start:y_end,x_start:x_end) contour_1500 = wrf_contour(a,wks,var_zoom_1500,optsM) delete (var_zoom_400) delete (var_zoom_600) delete (var_zoom_800) delete (var_zoom_1000) delete (var_zoom_1250) delete (var_zoom_1500) lnres = True lnres@gsLineThicknessF = 3.0 lnres@gsLineColor = "Red" if (ip.eq.1)then plot_400 = wrf_map_overlays(a,wks,(/contour_400/),pltres,mpres) do ii = 0,dimsX(0)-2 gsn_polyline(wks,plot_400,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) end do frame(wks) plot_600 = wrf_map_overlays(a,wks,(/contour_600/),pltres,mpres) do ii = 0,dimsX(0)-2 gsn_polyline(wks,plot_600,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) end do frame(wks) end if plot_800 = wrf_map_overlays(a,wks,(/contour_800/),pltres,mpres) do ii = 0,dimsX(0)-2 gsn_polyline(wks,plot_800,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) end do frame(wks) if (ip.eq.1)then plot_1000 = wrf_map_overlays(a,wks,(/contour_1000/),pltres,mpres) do ii = 0,dimsX(0)-2 gsn_polyline(wks,plot_1000,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) end do frame(wks) plot_1250 = wrf_map_overlays(a,wks,(/contour_1250/),pltres,mpres) do ii = 0,dimsX(0)-2 gsn_polyline(wks,plot_1250,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) end do frame(wks) plot_1500 = wrf_map_overlays(a,wks,(/contour_1500/),pltres,mpres) do ii = 0,dimsX(0)-2 gsn_polyline(wks,plot_1250,(/lon_plane(ii),lon_plane(ii+1)/),(/lat_plane(ii),lat_plane(ii+1)/),lnres) end do frame(wks) end if delete(lon_plane) delete(lat_plane) pltres@FramePlot = True end if plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres) ; plot x-section ; Delete options and fields, so we don't have carry over delete(opts_xy) delete(opts_tc) delete(opts_rh) delete(tc_plane) delete(rh_plane) delete(X_plane) end do ; make next cross section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTimeMap = False end do ; END OF TIME LOOP end