;Load vital libraries for NCL-WRF Functions. 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/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Begin Program and extract data to plot with. begin ;a=addfile(filename,"r") a=addfile("/tera9/brianjs/model/WRF_3.5/Thompson_MYJ/O_6_21_2010_Thompson_MYJ/wrfout_d01_2010-06-20_12:00:00.nc","r") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Specify retrieval data for all times and initiate TIME LOOP. times = wrf_user_getvar(a,"times",-1) ntimes = dimsizes(times) times_in_file = a->Times ; The specific pressure levels that we want the data interpolated to. height_levels = (/200,400,600,800,1000,1200,1400,1600,1800,2000/) nlevels = dimsizes(height_levels) ; number of height levels ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; do it = 3,8 title_file = chartostring(times_in_file(it,0:3)) +"_"+ chartostring(times_in_file(it,5:6)) +"_"+ chartostring(times_in_file(it,8:9)) +"_"+ chartostring(times_in_file(it,11:12)) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; First get the variables we will need terrain = a->HGT uvmet = wrf_user_getvar(a,"uvmet",-1) ; get rotated u and v comp of wind for all times z = wrf_user_getvar(a, "z",it) ; grid point height ter = wrf_user_getvar(a,"ter",it) ; model terrain height q_in = a->QVAPOR ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Interpolate/convert height field to AGL instead of SLP for proper plotting of level height field winds nheight = conform(z,ter,(/1,2/)) AGLheight = z - nheight copy_VarCoords(z,AGLheight) lat = a->XLAT(3,:,0) lon = a->XLONG(3,0,:) ;Retrieve u and v component of winds umet = uvmet(0,:,:,:,:) ; extract u component vmet = uvmet(1,:,:,:,:) ; extract v component ;Interpolate/Calculate wind barbs u_plane = wrf_user_intrp3d(umet(it,:,:,:),AGLheight,"h",height_levels,0.,False) v_plane = wrf_user_intrp3d(vmet(it,:,:,:),AGLheight,"h",height_levels,0.,False) q_out = wrf_user_intrp3d(q_in(it,:,:,:),AGLheight,"h",height_levels,0.,False)*1000 ;g/kg u_plane!0 = "height_levels" u_plane&height_levels = height_levels v_plane!0 = "height_levels" v_plane&height_levels = height_levels u_plane!1 = "lat" u_plane&lat = lat v_plane!1 = "lat" v_plane&lat = lat u_plane!2 = "lon" u_plane&lon = lon v_plane!2 = "lon" v_plane&lon = lon ;Apply gaussian filter to remove waves that 13 km RUC cannot resolve (i.e. 2delta-x or 2*13 = 26 km). For lat and lon, 13 weights are applied in each ;direction (for even division of 4 km spacing for 26 km waves). Standard deviation is assumed to be 1. ;weights = filwgts_normal(13, 1, 0) ;wgt = sqrt(weights/sum(weights)) ;u_lat_filter = wgt_runave_n_Wrap(u_plane,wgt,-1,1) ;u_lon_filter = wgt_runave_n_Wrap(u_plane,wgt,-1,2) ;v_lat_filter = wgt_runave_n_Wrap(v_plane,wgt,-1,1) ;v_lon_filter = wgt_runave_n_Wrap(v_plane,wgt,-1,2) ;q_lat_filter = wgt_runave_n_Wrap(q_outs,wgt,-1,1) ;q_lon_filter = wgt_runave_n_Wrap(q_outs,wgt,-1,2) ;Calculate mean (0-2 km) wind field ;u = u_lon_filter ;v = v_lon_filter ;q_out = q_lon_filter u = u_plane v = v_plane u_avg = (u(0,:,:)+u(1,:,:)+u(2,:,:)+u(3,:,:)+u(4,:,:)+u(5,:,:)+u(6,:,:)+u(7,:,:)+u(8,:,:)+u(9,:,:))/10 v_avg = (v(0,:,:)+v(1,:,:)+v(2,:,:)+v(3,:,:)+v(4,:,:)+v(5,:,:)+v(6,:,:)+v(7,:,:)+v(8,:,:)+v(9,:,:))/10 u_avg!0 = "lat" u_avg&lat = lat u_avg!1 = "lon" u_avg&lon = lon v_avg!0 = "lat" v_avg&lat = lat v_avg!1 = "lon" v_avg&lon = lon ;Calculate Convergence in a deep (0-2 km) layer dudx = new((/10,399,399/),float) dudx = 0 dvdy = new((/10,399,399/),float) dvdy = 0 dqdx = new((/10,399,399/),float) dqdx = 0 dqdy = new((/10,399,399/),float) dqdy = 0 do k = 0,9 do i=3,395 do j=3,395 dudx(k,i,j)=(u_plane(k,i,j+1)-u_plane(k,i,j-1))/8000 dvdy(k,i,j)=(v_plane(k,i+1,j)-v_plane(k,i-1,j))/8000 dqdx(k,i,j)=(q_out(k,i,j+1)-q_out(k,i,j-1))/8000 dqdy(k,i,j)=(q_out(k,i+1,j)-q_out(k,i-1,j))/8000 end do end do end do mfc_adv = u_plane mfc_adv = ((u_plane*dqdx) + (v_plane*dqdy))/1e-4 mfc_conv = u_plane mfc_conv = (q_out*(dudx+dvdy))/1e-3 mfc_tot = (mfc_adv+mfc_conv) c = mfc_tot deep_convergence = (c(0,:,:) + c(1,:,:) + c(2,:,:) + c(3,:,:) + c(4,:,:) + c(5,:,:) + c(6,:,:)+ c(7,:,:)+ c(8,:,:)+ c(9,:,:))/10 deep_convergence!0 = "lat" deep_convergence!1 = "lon" deep_convergence&lat = lat deep_convergence&lon = lon ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Establish Plots (type of plot, title, color map for contours, indications for manual changes). wks = gsn_open_wks("png","WRF_Deep_MFC_no_filter"+"_"+title_file) gsn_define_colormap(wks,"MPL_BuGn") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Restricting domain to lat/lon of choice (same as WRF) nlon = dimsizes(lon) nlat = dimsizes(lat) mpres = True mpres@gsnDraw = False mpres@gsnFrame = False mpres@gsnRightString = "" mpres@mpFillOn = False ; turn off map fill mpres@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last mpres@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries mpres@mpLimitMode = "Corners" ; Zoom in on the plot area. mpres@mpLeftCornerLonF = lon(0) mpres@mpRightCornerLonF = lon(nlon-1) mpres@mpLeftCornerLatF = lat(0) mpres@mpRightCornerLatF = lat(nlat-1) mpres@tmXTOn = "False" mpres@tmYROn = "False" mpres@mpGridAndLimbDrawOrder = "Predraw" ; set to predraw phase (mandatory) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Utilize coloring and labeling resources for the winds ; Add colors to windbarbs vcres = True vcres@gsnDraw = False vcres@gsnFrame = False vcres@vcLevelSelectionMode= "ExplicitLevels" ; select winds at certain vertical levels vcres@vcLevels = ispan(5,65,5) ; interval of winds to be plotted and colored vcres@vcRefAnnoOn = False ; turns off the ref vector vcres@vcRefLengthF = 0.045 ; set length of ref vector vcres@vcGlyphStyle = "WindBarb" ; turn on wind barbs vcres@vcMonoWindBarbColor = True ; color barbs by scalar vcres@vcWindBarbLineThicknessF = 2 ; wind barb thickness vcres@vcMinDistanceF = .07 vcres@vcVectorDrawOrder = "Predraw" ; set to predraw phase (mandatory) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;draw contours of convergence cnres = True cnres@gsnDraw = False cnres@gsnFrame = False cnres@cnLevelSelectionMode = "ExplicitLevels" cnres@cnLevels=(/0.2,0.4,0.6,0.8,1,2,4,6,9,10,12,14,16/) cnres@cnFillOn = True cnres@cnLinesOn = False cnres@cnInfoLabelOn = False cnres@cnFillDrawOrder = "Predraw" ; set to predraw phase (mandatory) cnres@lbOrientation = "vertical" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Text for the label bar txres = True txres@txFontHeightF = .017 ; Set the font height txres@txFontThicknessF = 2. label = "gkg~S~-1~N~" gsn_text_ndc(wks,label,.82,.735,txres) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Text for time label t_title = (/"","","","2100 UTC","0000 UTC","0300 UTC","0600 UTC","0900 UTC","1200 UTC"/) txres2 = True txres2@txFontHeightF = .016 ; Set the font height txres2@txFontThicknessF = 2. label = t_title(it) gsn_text_ndc(wks,label,.705,.740,txres2) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Overlay vectors and contours upon the map cnid = gsn_csm_contour(wks,deep_convergence,cnres) vcid = gsn_csm_vector(wks,u_avg,v_avg,vcres) ; Plot the wind barbs mpid = gsn_csm_map(wks,mpres) overlay(mpid,cnid) ; overlay convergence contours over base map overlay(mpid,vcid) ; overlay the wind barbs over base map maximize_output(wks,True) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end do end