; 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" begin ; ; The WRF ARW input file. ; This needs to have a ".nc" appended, so just do it. a = addfile("/public3/jtwiles/Tubbs_Fire_Analysis/TUBBS_FIRE_WRFOUT/wrfout_d01_2017-10-09_04:00:00.nc","r") ; We generate plots, but what kind do we prefer? type = "x11" ; type = "pdf" ; type = "png" ; type = "ncgm" wks = gsn_open_wks(type,"d1_X_sec_theta_Total_wind_Uandw_vector") ; Set some basic resources res = True res@MainTitle = "REAL-TIME WRF" res@vpWidthF = .9 ; overwrite basic plot size res@vpHeightF = 1.0 res@tmXBLabelFontHeightF = 0.05 res@tmYLLabelFontHeightF = 0.06 pltres = True ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FirstTime = 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) ;--------------------------------------------------------------- do it = 0,ntimes-1 ; TIME LOOP print("Working on time: " + times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots th = wrf_user_getvar(a,"th",it) ; theta rh = wrf_user_getvar(a,"rh",it) ; relative humidity tc = wrf_user_getvar(a,"tc",it) ; T in deg c u = wrf_user_getvar(a,"ua",it) v = wrf_user_getvar(a,"va",it) w = wrf_user_getvar(a,"wa",it) z = wrf_user_getvar(a, "z",it) ; grid point height p = wrf_user_getvar(a, "pressure",it) ; grid point height ; w_amplify = w ; w_amplify = 8*w ;--------------------------------------------------------------- ; Cross section opts = True ; specifying start and end points plane = new(4,float) plane = (/ 46,59 , 100,73 /) ;(/ 28,60 , 111,75 /) ; start x;y & end x;y point rh_plane = wrf_user_intrp3d(rh,p,"v",plane,0.,opts) th_plane = wrf_user_intrp3d(th,p,"v",plane,0.,opts) tc_plane = wrf_user_intrp3d(tc,p,"v",plane,0.,opts) u_plane = wrf_user_intrp3d(u,p,"v",plane,0.,opts) v_plane = wrf_user_intrp3d(v,p,"v",plane,0.,opts) w_plane = wrf_user_intrp3d(w,p,"v",plane,0.,opts) ; w_amplify_plane = wrf_user_intrp3d(w_amplify,p,"v",plane,0.,opts) w_amplify = w_plane w_amplify = w_plane*10 spd = (u_plane^2 + w_plane^2)^0.5 ; m/sec spd@description = "Wind Speed" spd@units = "m/s" u_wind = u_plane*0 ; Let's create nice labels - only have to do this once if ( FirstTime ) then zmax = 51. ; Place top at model top or close to zmax zz = wrf_user_intrp3d(p,p,"v",plane,0.,opts) ; z_ind = ind(.not.ismissing(zz(:,0))) ; zmin = zz(z_ind(0),0) ; delete(z_ind) zmin = 1000 nice_levs = floor((zmin-zmax)/50)*50 zmax = zmin - nice_levs dims = dimsizes(zz) zmax_pos = dims(0)-1 do imax = 1,dims(0)-1 if ( .not.ismissing(zz(imax,0)) .and. zz(imax,0) .ge. zmax ) then zmax_pos = imax end if end do zspan = zmax_pos zmax = zz(zmax_pos,0) nz = floattoint((zmin-zmax)/50+1) FirstTime = False end if ; Options for XY Plots opts_xy = res opts_xy@tiYAxisString = "Pressure (hPa)" opts_xy@tiXAxisString = "Distance (km)" opts_xy@AspectRatio = 0.75 opts_xy@cnMissingValPerimOn = True opts_xy@cnMissingValFillColor = 0 opts_xy@cnMissingValFillPattern = 11 opts_xy@tmYLMode = "Explicit" opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks opts_xy@tmYLLabels = sprintf("%.0f",fspan(zmin,zmax,nz)) ; Create labels opts_xy@tiXAxisFontHeightF = 0.020 opts_xy@tiYAxisFontHeightF = 0.010 opts_xy@tmXBMajorLengthF = 0.02 opts_xy@tmYLMajorLengthF = 0.02 opts_xy@tmYLLabelFontHeightF = 0.0001 opts_xy@PlotOrientation = th_plane@Orientation opts_xy@tiYAxisFontHeightF = 0.03 ;opts_xy@lgLabelFontHeightF = 1.5 ; Plotting options for RH opts_rh = opts_xy opts_rh@pmLabelBarOrthogonalPosF = -0.07 opts_rh@ContourParameters = (/ 0., 100., 10. /) opts_rh@cnFillOn = True opts@cnFillColors = (/"White","Chartreuse","Green",\ "Green3","Green4", \ "ForestGreen","PaleGreen4"/) ; Plotting options for Theta opts_th = opts_xy opts_th@cnInfoLabelOrthogonalPosF = 0.00 opts_th@gsnContourLineThicknessesScale = 2.0 opts_th@ContourParameters = (/ 290., 350., 2. /) ; opts_th@FieldTitle = "" opts_th@cnInfoLabelOn = False ; Plotting options for Temperature opts_tc = opts_xy opts_tc@cnLinesOn = False ; turn off contour lines opts_tc@ContourParameters = (/ 0., 30., 2. /) opts_tc@cnFillOn = True ; color plot desired ; Plotting options for Total Wind (U and V) opts_spd = opts_xy opts_spd@cnFillOn = True opts_spd@ContourParameters = (/ 0., 38., 2. /) opts_spd@pmLabelBarOrthogonalPosF = -0.1 opts_spd@cnFillPalette = "NCV_jaisnd" opts_spd@vcLevelPalette = "NCV_jaisnd" ; Plotting options for vertical velocity opts_t = opts_xy opts_t@cnLinesOn = False opts_t@ContourParameters = (/ -2., 2., .4/) opts_t@cnFillOn = True ; plotting option for u_wind = 0 ; u_wind = u ; u_wind = 0 ; u_wind_plane = wrf_user_intrp3d(u_wind,p,"v",plane,0.,opts) opts_u = opts_xy opts_u@cnLineColor = "Red" ; opts_u@cnInfoLabelOrthogonalPosF = 0.00 opts_u@gsnContourLineThicknessesScale = 3.0 ; Plotting options for Wind Vector opts = opts_xy opts@FieldTitle = "Wind vector" ; overwrite Field Title ; opts@vcRefAnnoOrthogonalPosF = -0.12 ; opts@vcRefAnnoJust = "BottomRight" ; opts@vcRefLengthF = 0.03 ; opts@vcWindBarbLineThicknessF = 6.0 opts@NumVectors = 40 ; wind barb density opts@vcRefAnnoOn = True ; define reference vector opts@vcRefMagnitudeF = 20.0 ; define vector ref mag opts@vcMinDistanceF = 0.04 ; Min distance between vectors opts@vcGlyphStyle = "LineArrow" opts@vcLineArrowThicknessF = 1.5 ; Get the contour info for the rh and temp contour_th = wrf_contour(a,wks,th_plane(0:zmax_pos,:),opts_th) contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),opts_rh) contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc) contour_w = wrf_contour(a,wks,w_plane(0:zmax_pos,:),opts_t) contour_spd = wrf_contour(a,wks,spd(0:zmax_pos,:),opts_spd) vector = wrf_vector(a,wks,u_plane,w_amplify,opts) contour_u = wrf_contour(a,wks,u_wind(0:zmax_pos,:),opts_u) ; MAKE PLOTS plot = wrf_overlays(a,wks,(/contour_th,vector,contour_spd/),pltres) ; Delete options and fields, so we don't have carry over delete(opts_th) delete(opts_rh) delete(opts_tc) delete(opts) delete(opts_t) delete(opts_spd) delete(opts_u) delete(th_plane) delete(rh_plane) delete(tc_plane) ; delete(u_plane) ; delete(v_plane) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do ; END OF TIME LOOP end