; 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 from a a given point A to point B ; Vertical coordinate is pressure load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;load "WRFUserARW.ncl" begin ;****** ;****** User Specified Pressure levels ;****** You can define whichever levels you want either ;****** using explicit array definition or something like fspan or ispan ;****** The levels you select here will automatically plot and zoom on ;****** the levels of interest ;****** NOTE: make sure levels are increasing (e.g. from bottom to top) ;****** ;;pnew = (/ 1000., 975., 950., 925., 900., 875., 850., 800., 750., 700., \ ;; 650., 600., 500., 400., 300., 200., 100. /) pnew = ispan(1000,800,2) * 1. nlevs = dimsizes(pnew) ;****** ;****** User Specified LAT/LON plane ;****** Can work for spanning a range of latitudes or a range of longitudes ;****** if you specify a range then it will average over the longitude dimension ;****** leaving a latitude line on the x-axis ;****** lats = (/60.5042, 60.5042/) lons = (/-150.500, -149.833/) ; lats = (/39.568611, 39.568611/) ; lons = (/-106.1794, -105.9794/) ;****** ;****** ;****** User Specified LAT/LON plane ;****** ;****** type = "x11" ;*************************************** ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. ;*************************************** a = addfile("../wrfprd/wrfout_d03_2019-11-30_04:30:00.nc","r") ;****** ;****** We generate plots, but what kind do we prefer? ;****** wks = gsn_open_wks(type,"CrossSection_horizontalspeed__W_E_0430UTC_domain3") ;;gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map ;****** ;****** What times and how many time steps are in the data set? ;****** FirstTime = True times = wrf_user_list_times(a) ; get times in the file ntimes = dimsizes(times) ; number of times in the file ;****** ;****** Some static plotting resources ;****** pres = True pres@gsnFrame = False pres@gsnMaximize = True pres@gsnPanelRowSpec = True pres@gsnPanelLabelBar = True pres@lbLabelAutoStride = True ;****** ;****** Begin time loop ;****** do it = 0,ntimes-1 ; TIME LOOP print("Working on time: " + times(it) ) ;----- First get the variables we will need uv = wrf_user_getvar(a,"uvmet",it) ; u wind component th1 = wrf_user_getvar(a,"th",it) ; theta uv@_FillValue = -999. ; w1 = wrf_user_getvar(a,"wa",it) ; w1@_FillValue = -999. p1 = wrf_user_getvar(a, "pressure",it) ; pressure is our vertical coordinate T1 = wrf_user_getvar(a,"th",it) ; theta (T+300) T1@_FillValue = -999. ; Calc speed, convert from m/s to kts ws = sqrt(uv(0,:,:,:)*uv(0,:,:,:) + uv(1,:,:,:)*uv(1,:,:,:)) * 1.9425 ws@description = "Wind Velocity" ws@units = "knots" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;; Select lat/lon plane and create output arrays ;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; plane = new(4,float) xlat = a->XLAT(0,:,:) xlon = a->XLONG(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 plane = (/ x_start,y_start , x_end,y_end /) ; approx. start x;y and end x;y point ndimz = dimsizes( xlat(y_start:y_end, x_start:x_end) ) nlat = ndimz(0) nlon = ndimz(1) if(nlon.eq.1) then w_plane = new( (/nlevs,nlat/), typeof(ws) ) th_plane = new( (/nlevs,nlat/), typeof(th1) ) idim = 1 end if if(nlat.eq.1) then w_plane = new( (/nlevs,nlon/), typeof(ws) ) th_plane = new( (/nlevs,nlon/), typeof(th1) ) idim = 0 end if if(nlon.gt.1 .and. nlat.gt.1) then print("Both lat and lon span a range: Averaging over rightmost dimension") w_plane = new( (/nlevs,nlat/), typeof(ws) ) th_plane = new( (/nlevs,nlat/), typeof(th1) ) idim = 1 end if w = ws (:,y_start:y_end, x_start:x_end) w_plane !0 = "lev" w_plane !1 = "dist" w_plane &lev = pnew w_plane &dist = ispan(0,dimsizes(w_plane(0,:))-1,1) w_plane &lev@units = "hPa" th = th1(:,y_start:y_end, x_start:x_end) T = T1 (:,y_start:y_end, x_start:x_end) th_plane!0 = "lev" th_plane!0 = "lev" th_plane!1 = "dist" th_plane!1 = "dist" th_plane&lev = pnew th_plane&lev = pnew th_plane&dist = ispan(0,dimsizes(th_plane(0,:))-1,1) th_plane&dist = ispan(0,dimsizes(th_plane(0,:))-1,1) th_plane&lev@units = "hPa" th_plane&lev@units = "hPa" p = p1 (:,y_start:y_end, x_start:x_end) ;******* ;******* Loop over cross sections ;******* do ip = 0, 0 ; we are just going to do one crossection for now ; Interpolate data vertically (in z) px = p do ll = 0,nlevs-1 w_plane (ll,:) = dim_avg_n( wrf_interp_3d_z( w ,px, pnew(ll)), idim ) th_plane(ll,:) = dim_avg_n( wrf_interp_3d_z( T ,px, pnew(ll)), idim ) th_plane(ll,:) = dim_avg_n( wrf_interp_3d_z( th,px, pnew(ll)), idim ) end do ;**** ;**** ;**** Options for XY Plots ;**** ;**** opts_xy = True opts_xy@gsnDraw = False opts_xy@gsnFrame = False opts_xy@gsnMaximize = True opts_xy@cnMissingValPerimOn = True opts_xy@cnMissingValFillColor = 0 opts_xy@cnMissingValFillPattern = 0 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@lbAutoManage = True opts_xy@cnInfoLabelOn = False ; Plotting options for Theta opts_th = opts_xy opts_th@cnFillOn = False opts_th@cnLevelSpacingF = 0.25 ; set contour spacing ; Plotting options for Wind Speed opts_ws = opts_xy opts_ws@cnLevelSpacingF = 2 ; set contour spacing opts_ws@cnFillOn = True opts_ws@gsnSpreadColors = True ; Plotting options for Temperature opts_th = opts_xy opts_th@cnFillOn = False ; Color filled plot resources (W) opts_xy@tiXAxisString = "Number of Grid Points" opts_xy@tiYAxisString = "Pressure (mb)" opts_xy@tiMainString = "Horizontal Speed "+times(it) opts_xy@cnFillOn = True opts_xy@cnLevelSelectionMode = "ManualLevels" opts_xy@cnMaxLevelValF = 50. opts_xy@cnMinLevelValF = 10. opts_xy@cnLevelSpacingF = 2 opts_xy@cnLinesOn = False ; no contour lines opts_xy@cnLineLabelsOn = False ; no contour labels opts_xy@gsnSpreadColors = True opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = ispan(0,nlevs-1,1) * 1. opts_xy@tmYLLabels = w_plane&lev+"" ; Can change out variables to plot below ; CAUTION: I set explicit contour levels above for the fill figure ; so if you change out Ri Number then change contour levels above ; Make Plots but don't draw them yet plot_fill = gsn_contour( wks, w_plane , opts_xy ) plot_line = gsn_contour( wks, th_plane, opts_th ) ; Overlay plots overlay( plot_fill, plot_line ) ; Plot and send to output type gsn_panel(wks,(/ plot_fill /), (/1/), pres) frame(wks) ; ; Plot and send to output type ; gsn_panel(wks,(/ plot_fill /), (/1/), pres) ; frame(wks) end do ; make next cross section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do ; END OF TIME LOOP end