;LOAD VITAL LIBRARIES FOR WRF-NCL FUNCTIONS AND LAMONT, OK VERTICAL PROFILER OBSERVATIONAL DATA 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 IMPORT DATA TO PLOT. begin a=addfile("/tera9/brianjs/profiler/lamont/6_15_2012_cat/cat.cdf","r") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; create time information and strip out the day and hour (insert into wind data arrays later) times2 = ispan(0,12,1) time = (/"06/14-12","","06/14-18","","06/15-00","","06/15-06","","06/15-12","","06/15-18","","06/16-00"/) ;print(time) ;printVarSummary(time) ;print(time) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;EXTRACT WIND COMPONENTS NEEDED FOR NCL, EXCLUDE BAD DATA, CONVERT UNITS, AND REARRANGE VARIABLE ARRAYS ;extract height variable to assign to wind arrays and plot on y-axis ;choose lowest 40 levels (0-39) for the low powered setting to evaluate lowest level of atmosphere in detail h = a->height_p(0:39,0) ;setup height string to explicitly label y-axis ht = (/"0","","","","","0.5","","","","","1.0","","","","","1.5","","","","","2.0","","","","","2.5",\ "","","","","3.0"/) ;print(h) ;retrieve u and v components (3-hrly) of wind as needed in NCL and exclude bad data (999999 or -9999) u = a->u_wind(12:48:3,0:39,0) u_new = u u_new = where(u.eq.999999,-9999,u) ;u_new = where(u.gt.80,-9999,u) v = a->v_wind(12:48:3,0:39,0) v_new = v v_new = where(v.eq.999999,-9999,v) ;v_new = where(v.gt.80,-9999,v) ;printVarSummary(u_new) ;printVarSummary(u) ;print(u) ;print(v) ;exit ;newtime = a->time(12:48:3) ;cal = cd_calendar(newtime,-3) ;print(cal) ;convert to knots and assign time and height arrays defined earlier ukt = u_new ukt = u_new*1.94386 ukt!0 = "times2" ukt×2 = times2 ukt!1 = "h" ukt&h = h vkt = v_new vkt = v_new*1.94386 vkt!0 = "times2" vkt×2 = times2 vkt!1 = "h" vkt&h = h ;print(h) ;printVarSummary(ukt) ;printVarSummary(vkt) ;printMinMax(ukt,True) ;printMinMax(vkt,True) ;obtain wind speed (magnitude) and assign time and height arrays defined earlier wnd = ukt(:,0:39) wnd = sqrt((ukt(:,0:39)^2)+(vkt(:,0:39)^2)) wnd_new = wnd wnd_new = where(wnd.gt.75,-9999,wnd) wnd_new!0 ="times2" wnd_new×2 = times2 wnd_new!1 = "h" wnd_new&h = h ;print(h) ;printVarSummary(wnd) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;PLOT OUTPUT AS A POST SCRIPT, NAME THE FILE, AND CHOOSE AN NCL COLOR TABLE TO PLOT THE WINDS wks = gsn_open_wks("x11","lamont__wind_profiler_June_15_2012") gsn_define_colormap(wks,"WhiteBlueGreenYellowRed") ; choose color map ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;DESIGN THE VISUAL OUTPUT FOR THE PLOT USING RESOURCE FUNCTIONS res2D = True ; Set basic resources res2D@gsnDraw = False ; Don't draw individual plot. res2D@gsnFrame = False ; Don't advance frame. res2D@vpXF = 0.15 ; x location res2D@vpYF = 0.90 ; y location res2D@vpWidthF = 0.70 ; plot width res2D@vpHeightF = 0.40 ; plot height res2D@tiMainString = "Lamont, OK Data Wind Profile ~C~ ~Z75~ 2012/06/15 Case" res2D@tiMainFontThicknessF = 4.;2. ; title thickness res2D@tiXAxisFontThicknessF = 3.;2. ; x axis title font thickness res2D@tiYAxisFontThicknessF = 3.;2. ; y axis title font thickness res2D@tiXAxisString = "Month/Day-Hour (UTC)" ; x axis title res2D@tiYAxisString = "Height (km)" ; y axis title res2D@tiXAxisFontHeightF = 0.016 ; x axis title font size res2D@tmXBMode = "Explicit" ; x axis marker specifications res2D@tmXBLabelFontThicknessF = 2. ; x axis label font thickness res2D@tmXBValues = times2 ; x axis label values res2D@tmXBLabels = time ; x axis labels res2D@tmXBLabelJust = "CenterCenter" ; x axis positions res2D@tmXBLabelFontHeightF = .012 ; x axis label fint height res2D@tmYROn = False ; y axis tick marks on RHS turned off res2D@tmXTOn = False ; x axis top of graph turned off res2D@tmYLMode = "Manual" ; y axis marker specifications res2D@tmYLValues = h ; y values are the heights res2D@tmYLLabels = ht res2D@tmYLTickStartF = 0. res2D@tmYLTickEndF = 3. res2D@tmYLTickSpacingF = .5 res2D@tmYMaxF = 3. res2D@tmYLDataBottomF = 0.0 ; y axis lowest value res2D@tmYLDataTopF = 3.0 ; y axis highest value res2D@tmYLLabelFontThicknessF = 2. ; y axis label font thickness res2D@vcLevelSelectionMode = "ExplicitLevels" ; select winds at certain vertical levels res2D@vcLevels = ispan(5,65,5) ; interval of winds to be plotted and colored res2D@vcRefAnnoOn = False ; turns off the ref vector res2D@vcMinDistanceF = .01; 0.025 ; density of wind barbs res2D@vcRefLengthF = 0.040 ; set length of ref vector res2D@vcGlyphStyle = "WindBarb" ; turn on wind barbs res2D@vcMapDirection = False ; render directional coordinate in vertical profile res2D@gsnSpreadColors = True ; use full colormap res2D@gsnSpreadColorStart = 48;-1 ; don't use added gray res2D@gsnSpreadColorEnd = -2 ; don't use added gray res2D@gsnLeftString = "" res2D@gsnRightString = "" ; "knots" res2D@vcGlyphStyle = "WindBarb" ; choose wind barbs res2D@vcMonoWindBarbColor = False ; color barbs by scalar res2D@vcWindBarbLineThicknessF = 4.;3;2. ; wind barb thickness res2D@lbOrientation = "vertical" ; Label Bar orientation res2D@lbLabelFontThicknessF = 2. ; Label bar font thickness res2D@pmLabelBarWidthF = .075 ; Label bar width ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Text for the label bar txres = True txres@txFontHeightF = .015 ;.02 ; Set the font height txres@txFontThicknessF = 2. ;1.5 ;2. label = "knots" gsn_text_ndc(wks,label,.9,.81,txres) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;This trick reverses the var coords so height is first followed by time (to generate correct plot) ukt_ht_time = ukt(h|:,times2|:) vkt_ht_time = vkt(h|:,times2|:) wnd_ht_time = wnd_new(h|:,times2|:) ;printVarSummary(ukt_ht_time) ;printVarSummary(vkt_ht_time) ;printVarSummary(wnd_ht_time) windlayer = gsn_csm_vector_scalar(wks,ukt_ht_time,vkt_ht_time,wnd_ht_time,res2D) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Overlay current map with larger frame (nesting current data in larger set) getvalues windlayer ; get the X/Y axis min/max for use in the loglin plot "trXMinF" : trxmin "trXMaxF" : trxmax "trYMinF" : trymin "trYMaxF" : trymax end getvalues loglin = create "logling" logLinPlotClass wks ; draw a loglin plot, with expanded X/Y axis "trXMinF" : trxmin-1 "trXMaxF" : trxmax+1 "trYMinF" : trymin "trYMaxF" : 3.0 ; the next 4 resource settings may not be needed "vpXF" : .15 ; set the X-axis NDC starting point "vpYF" : .8 ; set the Y-axis NDC starting point "vpWidthF" : .7 ; set the width of the plot in NDC units "vpHeightF" : .45 ; set the height of the plot in NDC units end create overlay(loglin,windlayer) ; overlay plot with the loglin plot draw(loglin) ; draw the plot frame(wks) ; advance the frame end