;----------------------------------------------------------------------
; wrf_nogsn_6.ncl
;----------------------------------------------------------------------
; Concepts illustrated:
;   - Using wrf_xxxx scripts to plot WRF-ARW data
;   - Overlaying filled contours, and vectors on a map
;   - Setting the correct WRF map projection using wrf_map_resources
;   - Setting lots of resources to customize a WRF plot
;----------------------------------------------------------------------
; This script is meant to show the difference between plotting WRF
; data using wrf_xxx scripts, and using gsn_csm_xxx scripts.
;----------------------------------------------------------------------
; These files are loaded by default in NCL V6.2.0 and newer
; 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
;---Open WRF output file
  filename = "/global/scratch/umtropea/conus/ctrl/2d/2012/wind/wrf2d_d01_CTRL_U_201210-201212.nc"
  a        = addfile(filename,"r")
  filename2 = "/global/scratch/umtropea/conus/ctrl/2d/2012/wind/wrf2d_d01_CTRL_V_201210-201212.nc"
  b        = addfile(filename2,"r")

;---Read several WRF variables at first time step
  u   = wrf_user_getvar(a,"U",0)    ; 2D U at mass points
  v   = wrf_user_getvar(b,"V",0)    ; 2D V at mass points
  
;Destagger
  u_unstag = wrf_user_unstagger(u,u@stagger)
  v_unstag = wrf_user_unstagger(v,v@stagger)  

;---Now get the lowest (bottommost) level?
  udata = u_unstag
  vdata = v_unstag
  udata2 = udata*1.94386                    ; Convert wind into knots
  vdata2 = vdata*1.94386

;---Change the metadata
  udata2@units      	  = "kts"
  vdata2@units            = "kts"

  wks = gsn_open_wks("x11", "NCL_BARB_TEST1")

;---Wind vector plot
  vec_res                       = True
  vec_res@vcLineArrowThicknessF = 3.0
  vec_res@vcGlyphStyle          = "CurlyVector"

  vector = wrf_vector(a,wks,udata2,vdata2,vec_res)
  
;---Overlay plots on map and draw.
  map = wrf_map_overlays(a,wks,(/vector/),True,True)

end