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" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl" begin ;; creates a generic cross section plot. graphics_type = "png" graphics_file = "cross_section_test" wrf_dir = "/media/schaefer_hd/FIRE/" wrf_file = "wrfout_d03_2012-08-29_22:20:00" + ".nc" ; your default wrfout files ; won't have an NC suffix dir = "/media/schaefer_hd/" fire_folder = "FIRE/" sensible_folder = "SENSIBLE/" files = systemfunc("ls " + dir + fire_folder + "wrfout_d03_2012-08-29_22:00:00") nfiles = dimsizes(files) ; open the wrf files sf = addfile(wrf_dir + wrf_file, "r") ; get Time coordinates for WRF files Time_strings = chartostring( sf->Times ) Time = wrf_times_c( sf->Times , 0 ) ; classic netcdf time coordinates ;print(Time) ; get pressure and height (two frequently used vertical layers) wrf_z3d = wrf_user_getvar(sf, "height", -1) wrf_p3d = wrf_user_getvar(sf, "pressure", -1) ; get longitude and latitude wrf_lat2d = wrf_user_getvar(sf, "lat", 0) ; our grid doesn't move... wrf_lon2d = wrf_user_getvar(sf, "lon", 0) wks = gsn_open_wks(graphics_type, graphics_file) do nf = 0 ,nfiles-1 f_1 = addfile(files(nf) + ".nc" ,"r") fname_2 = str_sub_str(files(nf),"FIRE","SENSIBLE") f_2 = addfile(fname_2 + ".nc" ,"r") ; get a few fields that go well in cross sections wrf_theta3d = wrf_user_getvar(f_1, "theta", -1) wrf_rh3d = wrf_user_getvar(f_1, "rh", -1) wrf_u3d = wrf_user_getvar(f_1, "ua", -1) ; unstaggared wrf_v3d = wrf_user_getvar(f_1, "va", -1) ; unstaggared wrf_w3d = wrf_user_getvar(f_1, "wa", -1) ; unstaggared ; now get dimensions ntdsdydx = dimsizes(wrf_z3d) nt = ntdsdydx(0) ; time steps ns = ntdsdydx(1) ; sigma levels ny = ntdsdydx(2) nx = ntdsdydx(3) ;;;;;;; Now this is where you define the crossections! ;;;;;;;;;;;;;;;;;;;;;;;; ; ; select your crossection region... ; in this example we are doing a simple crossection along "x" ; ; if you want to to south->north then change the code accordingly ; i_start = 66 i_end = 167 di = i_end - i_start j_start = 60 j_end = 161 dj = j_end - j_start nx_cross = j_end - j_start + 1 ;;;;;;; And now we do the vertical levels for the crossection ;;;;;;;;;;;;;;;;;; ;;;;;;; (change as you need) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; nz_cross = 1001 ; My terrain avgerages 1400 so I am subtracting that off for ; AGL heights....you could subtract wrf_user_getvar "ter" ; from your heights to get the correct AGL height_levels = fspan(1400,9900,nz_cross) * 1e-3 ; makes a BIG height array height_levels@units = "km" ; from 0->8500m AGL by 8.5m steps height_levels@description = "Height MSL (km)" ; These heights are MSL as that is the waw ; WRF stores them....convert to AGL below ; since height_levels is a true monotonic coordiate, we can and SHOULD name its ; dimension after itself. height_levels!0 = "height_levels" height_levels&height_levels = height_levels ;========================================================================= ; This is setting up the heights to be AGL in the end height_levels2 = fspan(0,8500,nz_cross) * 1e-3 height_levels2@units = "km" height_levels2@description = "Height AGL (km)" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;; Create some arrays for the transects ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; theta_cross = new( (/ nt, nz_cross, nx_cross /), float) theta_cross!0 = "Time" theta_cross&Time = Time theta_cross!1 = "height_levels" theta_cross&height_levels = height_levels ua_cross = theta_cross va_cross = theta_cross wa_cross = theta_cross rh_cross = theta_cross theta_cross@description = "Potential Temperature" ua_cross@description = "U-Wind" va_cross@description = "V-Wind" wa_cross@description = "W-Wind" rh_cross@description = "Relative Humidity" theta_cross@units = "K" ua_cross@units = "m s-1" va_cross@units = "m s-1" wa_cross@units = "m s-1" rh_cross@units = "%" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;; Now make the crossecitons ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; do t = 0, nt-1 do j = 0, nx_cross-1 j_st = j + j_start i_st = j + i_start theta_cross(t,:,j) = (/ wrf_interp_1d(wrf_theta3d(t,:,j_st,i_st), \ wrf_z3d(t,:,j_st,i_st)/1000, \ height_levels ) /) ua_cross(t,:,j) = (/ wrf_interp_1d(wrf_u3d(t,:,j_st,i_st), \ wrf_z3d(t,:,j_st,i_st)/1000, \ height_levels ) /) va_cross(t,:,j) = (/ wrf_interp_1d(wrf_u3d(t,:,j_st,i_st), \ wrf_z3d(t,:,j_st,i_st)/1000, \ height_levels ) /) wa_cross(t,:,j) = (/ wrf_interp_1d(wrf_w3d(t,:,j_st,i_st), \ wrf_z3d(t,:,j_st,i_st)/1000, \ height_levels ) /) rh_cross(t,:,j) = (/ wrf_interp_1d(wrf_rh3d(t,:,j_st,i_st), \ wrf_z3d(t,:,j_st,i_st)/1000, \ height_levels ) /) end do end do ;;;;;;; The interp function messes up the dimensions so we need to fix em ;;;;;; theta_cross!0 = "Time" theta_cross&Time = Time ua_cross!0 = "Time" ua_cross&Time = Time va_cross!0 = "Time" va_cross&Time = Time wa_cross!0 = "Time" wa_cross&Time = Time theta_cross!1 = "height_levels" theta_cross&height_levels = height_levels2 ua_cross!1 = "height_levels" ua_cross&height_levels = height_levels2 va_cross!1 = "height_levels" va_cross&height_levels = height_levels2 wa_cross!1 = "height_levels" wa_cross&height_levels = height_levels2 t = 0 i = 0 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; plot stuff ; wks = gsn_open_wks(graphics_type, graphics_file) ; open a graphics unit gsn_define_colormap(wks,"BlAqGrYeOrRevi200") ; choose color map res = True res@gsnDraw = False ; Don't draw individual plot. res@gsnFrame = False ; Don't advance frame. res@gsnMaximize = True res@vpHeightF = .35 res@vpWidthF = .9 res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks res@cnFillOn = True ; turn on color for contours res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off contour line labels res@gsnSpreadColors = True ; use full color map res@lbLabelAutoStride = True res@cnMinLevelValF = 320 res@cnMaxLevelValF = 346 res@cnLevelSpacingF = 1 res@cnLevelSelectionMode = "ManualLevels" res@sfYArray = height_levels2 ; Assigning AGL values to y-axis resvector = True resvector@gsnDraw = False ; Don't draw individual plot. resvector@gsnFrame = False ; Don't advance frame. resvector@vcRefMagnitudeF = 5.0 ; make vectors larger resvector@vcRefLengthF = 0.050 ; ref vec length resvector@vcGlyphStyle = "CurlyVector" ; turn on curly vectors resvector@vcMinDistanceF = 0.025 ; thin out vectors resvector@vcRefAnnoOn = False resvector@gsnLeftString = "" ; needed for the overlay resvector@gsnRightString = "" ; needed for the overlay res_hum = True res_hum@gsnDraw = False ; Don't draw individual plot. res_hum@gsnFrame = False ; Don't advance frame. ;res_hum@gsnShadeFillType = "pattern" ; pattern fill ;res_hum@gsnShadeMid = 12 ; use pattern 17 (stippling) res_hum@cnLinesOn = True ; turn off contour lines res_hum@cnLineLabelsOn = True ; turn off contour line labels res_hum@cnLevelSelectionMode = "ManualLevels" res_hum@cnMinLevelValF = 5 res_hum@cnMaxLevelValF = 70 res_hum@cnLevelSpacingF = 5 res_hum@gsnLeftString = "" ; needed for the overlay res_hum@gsnRightString = "" ; needed for the overlay res_hum@cnLineThicknessF = 1.5 res_hum@cnLineLabelPerimOn = False res_hum@cnLineLabelInterval = 4 res_hum@cnInfoLabelOn = False res_hum@cnLineLabelBackgroundColor = -1 res_hum@cnLineColor = "white" res_hum@cnConstFLabelOn = False do t = 0, nt-1 res@tiMainString = "Plume Transect" plot = gsn_csm_contour(wks,theta_cross(t,:,:),res) ; contour the variable plotvector = gsn_csm_vector(wks,ua_cross(t,:,:),wa_cross(t,:,:),resvector) overlay(plot,plotvector) ; result is windlayer over the other two ; makes a stipple pattern for rh's between 90% and 100% plotrh = gsn_csm_contour(wks,rh_cross(t,:,:),res_hum) ;plotrh = gsn_contour_shade(plotrh,5.,70.,res_hum) overlay(plot,plotrh) ; draw(plot) ; draw this overlay result frame(wks) ; finally advance the frame once all on page end do end do end