; 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 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. ; a = addfile("../wrfout_d01_2000-01-24_12:00:00.nc","r") a = addfile("/gpfs/summit/scratch/nilu4970/WRFTEST_4/WRF_500m_valley_FC_1/test/em_fire/wrfout_d02_0001-01-01_00:00: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_winds_theta_contour_300_2km") ; Set some basic resources res = True res@MainTitle = "WRF 500 m Valley" res@Footer = False pltres = True ter_res = True opts_ter = res opts_ter@gsnYRefLine = 0.0 opts_ter@gsnAboveYRefLineColor = "black" opts_ter@gsnDraw = False opts_ter@gsnFrame = False opts_ter@FieldTitle = "" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 printVarSummary(mdims) nd = dimsizes(mdims) printVarSummary(nd) ;--------------------------------------------------------------- ter = wrf_user_getvar(a,"HGT",0) xlon = wrf_user_getvar(a, "XLONG",0) printVarSummary(xlon) do it = 100,ntimes-1 ; TIME LOOP print("Working on time: " + times(it) ) res@TimeLabel = times(it) ; Set Valid time to use on plots u = wrf_user_getvar(a,"ua",it) ; ua is u averaged to mass points v = wrf_user_getvar(a,"va",it) ; va is v averaged to mass points tc = wrf_user_getvar(a,"th",it) tc_con = u(:10,:,:) ; T in C rh = wrf_user_getvar(a,"rh",it) ; relative humidity t_con = tc(:10,:,:) rh_con = v(:10,:,:) printVarSummary(rh) z = wrf_user_getvar(a,"z",it) ; grid point height printVarSummary(z) z_con = z(:10,:,:) printVarSummary(z_con) if ( FirstTime ) then ; get height info for labels zmin = 0. zmax = max(z_con)/1000. nz = floattoint(zmax/2 + 10) ;nz = floattoint(zmax/2) FirstTime = False end if print(zmax) z_values = fspan(zmin,zmax,nz) printVarSummary(z_values) dimzz = dimsizes(z_values) print(dimzz) opts_ter@trYMaxF = z_values(dimzz-2)*1000 ;--------------------------------------------------------------- do ip = 4, 4 ; 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 | ; | ; if(ip .eq. 1) then angle = 90. end if if(ip .eq. 2) then angle = 0. end if if(ip .eq. 3) then angle = 45. end if if(ip .eq. 4) then angle = 300. end if plane = new(2,float) plane = (/ mdims(nd-2)/2, mdims(nd-2)/2 + 20 /) ; pivot point is center of domain (x,y) opts = False rh_plane = wrf_user_intrp3d(rh_con,z_con,"v",plane,angle,opts) tc_plane = wrf_user_intrp3d(tc_con,z_con,"v",plane,angle,opts) V = sqrt(rh_plane^2 + tc_plane^2) t_plane = wrf_user_intrp3d(t_con,z_con,"v",plane,angle,opts) printVarSummary(V) printVarSummary(rh_plane) ter_plane = wrf_user_intrp2d(ter,plane,angle,opts) X_plane = wrf_user_intrp2d(xlon,plane,angle,opts) printVarSummary(X_plane) rh_plane2 = V tc_plane2 = t_plane cross_dims = dimsizes(rh_plane2) rank = dimsizes(cross_dims) ;printVarSummary(rh_plane2) iz_do = 50 do iz = 0,1 iz_do = iz_do-1 do ix = 0,cross_dims(rank-1)-1 if ( ismissing(rh_plane2(iz_do,ix)) ) then rh_plane2(iz_do,ix) = rh_plane2(iz_do+1,ix) end if if ( ismissing(tc_plane2(iz_do,ix)) ) then tc_plane2(iz_do,ix) = tc_plane2(iz_do+1,ix) end if end do end do dim = dimsizes(rh_plane) ; Find the data span - for use in labels zspan = dim(0) ; Options for XY Plots opts_xy = res opts_xy@tiYAxisString = "Height (km)" 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("%.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@FieldTitle = "Winds (m/s)" opts_rh@description = "Winds (m/s)" opts_rh@ContourParameters = (/ 0., 18., .1 /) opts_rh@pmLabelBarOrthogonalPosF = -0.07 opts_rh@cnFillOn = True ;opts_rh@cnFillColors = (/"White","White","White", \ ; "White","Chartreuse","Green", \ ; "Green3","Green4", \ ; "ForestGreen","PaleGreen4"/) ; opts_rh@cnFillPalette = True opts_rh@cnFillPalette = "MPL_blues" ; Plotting options for Temperature opts_tc = opts_xy ; opts_tc@cnFillOn = True opts_tc@pmLabelBarOrthogonalPosF = -0.07 opts_tc@ContourParameters = (/ 290., 320., 1 /) ; opts_tc@cnFillPalette = "MPL_coolwarm" ; Get the contour info for the rh and temp contour_tc = wrf_contour(a,wks,tc_plane,opts_tc) contour_rh = wrf_contour(a,wks,rh_plane,opts_rh) contour_rh2 = wrf_contour(a,wks,rh_plane2,opts_rh) contour_tc2 = wrf_contour(a,wks,tc_plane2,opts_tc) ;Contour terrain cross section contour_ter = gsn_csm_xy(wks,X_plane,ter_plane,opts_ter) ; MAKE PLOTS ;plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres) ;plot = wrf_overlays(a,wks,(/contour_rh,contour_tc,contour_ter/),pltres) ;plot = wrf_overlays(a,wks,(/contour_rh2,contour_tc2/),pltres) plot = wrf_overlays(a,wks,(/contour_rh2,contour_ter,contour_tc2/),pltres) ; Delete options and fields, so we don't have carry over delete(opts_tc) delete(opts_rh) delete(tc_plane) delete(tc_plane2) delete(rh_plane) delete(rh_plane2) delete(X_plane) delete(ter_plane) end do ; make next cross section ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do ; END OF TIME LOOP end