;----------------------------------------------------------------------
; This script overlays vectors over color contours over a map.
;
; It additionally zooms in over US/Canada and adds some shapefile
; outlines. The outlines are drawn in a thick white line you can
; see them. It's not necessarily the best choice in general for
; publication!
;----------------------------------------------------------------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
;---Open file and read data
  filename = "~/ncl/data/20160115-18/rap_130_20160115_0000_000.grb2"
  a    = addfile(filename,"r")
  refc = a->REFC_P0_L200_GLC0
  u    = a->UGRD_P0_L7_GLC0
  v    = a->VGRD_P0_L7_GLC0

  lat = a->gridlat_0
  lon = a->gridlat_0

;---Start the graphics
  wks = gsn_open_wks("x11","mov6")

;---Set some contour resources
  cnres                   = True

  cnres@gsnMaximize       = True     ; maximize plot in frame
  cnres@gsnDraw           = False    ; turn off draw
  cnres@gsnFrame          = False    ; turn off frame

  cnres@cnFillOn          = True     ; turn on color fill
  cnres@cnLinesOn         = False    ; turn off contour lines
  
  cnres@gsnLeftString     = "REFC ("  + refc@units + ")"
  cnres@gsnRightString    = "UV winds ("  + u@units + ")"


  cnres@mpOutlineOn       = False    ; We're going to add 
  cnres@mpFillOn          = False    ; shapefile outlines.

  cnres@mpMinLatF         =   35
  cnres@mpMaxLatF         =   47
  cnres@mpMinLonF         =  -80
  cnres@mpMaxLonF         =  -67

  contour_plot = gsn_csm_contour_map(wks,refc,cnres)

;---Set some vector resources
  vcres                   = True

  vcres@gsnMaximize       = True     ; maximize plot in frame
  vcres@gsnDraw           = False    ; turn off draw
  vcres@gsnFrame          = False    ; turn off frame

  vcres@vcMinDistanceF    = 0.005    ; thin the vectors
  vcres@vcRefLengthF      = 0.05     ; width/num_vectors (~25 vectors)

  vcres@vcLineArrowColor  = "Gray25"

  vcres@gsnRightString    = ""       ; turn off sub-titles
  vcres@gsnLeftString     = ""

  vector_plot = gsn_csm_vector(wks,u,v,vcres)

  overlay(contour_plot,vector_plot)

;---Attach shapefile polylines to contour_plot
  lnres                  = True 
  lnres@gsLineColor      = "White"
  lnres@gsLineThicknessF = 2.5
  id1 = gsn_add_shapefile_polylines(wks,contour_plot,"~/ncl/data/shape/statesp020.shp",lnres)
  
  
;---Drawing the contour plot will draw the two sets of polylines.
  draw(contour_plot)
  frame(wks)
end