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
  if(.not.isvar("ndc_overlay")) then
    ndc_overlay = True
  end if

  if(.not.isvar("map_direction")) then
    map_direction = True
  end if

;----------------------------------------------------------------------
; Get the data. You can set the filename on the command line, if 
; desired:
;  ncl 'filename="wrfout_d01_2003-07-13_12:00:00"' wrf_vectors.ncl
;----------------------------------------------------------------------
  if(.not.isvar("filename")) then
    filename = "wrfout_d01_2011-08-08_12:00:00"
  end if
  print("==================================================")
  print("file: " + filename)
  f = addfile(filename+".nc","r")

  u     = wrf_user_getvar(f,"ua",0) 
  v     = wrf_user_getvar(f,"va",0) 
  pf    = wrf_user_getvar(f,"pressure",0) ; Full model pressure [hPa]
  u_950 = wrf_interp_3d_z(u,pf,950.)
  v_950 = wrf_interp_3d_z(v,pf,950.)
  lat2d = f->XLAT(0,:,:)
  lon2d = f->XLONG(0,:,:)

  nlat = dimsizes(lat2d(:,0))
  nlon = dimsizes(lat2d(0,:))
  min_lat_diff = abs(min(lat2d(1,:)-lat2d(0,:)))
  min_lon_diff = abs(min(lon2d(:,1)-lon2d(0,0)))
  max_lat_diff = min_lat_diff
  max_lon_diff = min_lon_diff
  do n=2,nlat-1
    min_lat_diff = abs(min((/min_lat_diff,min(lat2d(n,:)-lat2d(n-1,:))/)))
    max_lat_diff = abs(max((/max_lat_diff,max(lat2d(n,:)-lat2d(n-1,:))/)))
  end do
  do n=2,nlon-1
    min_lon_diff = abs(min((/min_lon_diff,min(lon2d(:,n)-lon2d(:,n-1))/)))
    max_lon_diff = abs(max((/max_lon_diff,max(lon2d(:,n)-lon2d(:,n-1))/)))
  end do
  print("min lat diff = " + min_lat_diff)
  print("max lat diff = " + max_lat_diff)
  print("min lon diff = " + min_lon_diff)
  print("max lon diff = " + max_lon_diff)
;  printMinMax(u_950,0)
;  printMinMax(v_950,0)

;----------------------------------------------------------------------
; Set the plotting options
;----------------------------------------------------------------------
  res                         = True
  res@gsnMaximize             = True
  res@gsnAddCyclic            = False          ; don't add a longitude cyclic point

  res@tfDoNDCOverlay          = ndc_overlay    ; True / False 

  res@pmTickMarkDisplayMode   = "Always"       ; turn on nice tickmarks

  res@vcMapDirection          = map_direction  ; True / False 
  res@vcRefMagnitudeF         = 8.             ; define vector ref mag
  res@vcRefLengthF            = 0.025          ; define length of vec ref
  res@vcMinDistanceF          = 0.02           ; larger means sparser
  res@vcGlyphStyle            = "FillArrow"
  res@vcFillArrowEdgeColor    = "black"
  res@vcFillArrowEdgeThicknessF = 1.0

  res = wrf_map_resources(f, res)

  if(.not.ndc_overlay) then
    u_950@lat2d = lat2d
    u_950@lon2d = lon2d
    v_950@lat2d = lat2d
    v_950@lon2d = lon2d
  end if

  res@tiMainString   = filename
  res@gsnLeftString  = "tfDoNDCOverlay=" + ndc_overlay
  res@gsnRightString = "vcMapDirection=" + map_direction

  figurename = str_sub_str(filename,":","_")

;----------------------------------------------------------------------
; Create the plot
;----------------------------------------------------------------------
  wks = gsn_open_wks("png",figurename + "_ndc_" + ndc_overlay + "_mapdir_" + map_direction)
  vector_plot  = gsn_csm_vector_map(wks,u_950,v_950,res)
end