;----------------------------------------------------------------------
; 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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"


begin
;---Make a list of all files we are interested in
  DATADir = "~/ncl/data/20160115-18/"
  FILES = systemfunc (" ls -1 " + DATADir + "rap* ")
  numFILES = dimsizes(FILES)
  do ifil = 0,numFILES-1
;---Start the graphics
  wks = gsn_open_wks("png","jan16storm"+ifil)
  a = addfile(FILES(ifil)+".grb2","r") ; Open the next file
  lat2d   = f->gridlat_0
  lon2d   = f->gridlon_0
  rot     = f->gridrot_0         ; (lat,lon)
  refc = a->REFC_P0_L200_GLC0
  u    = a->UGRD_P0_L7_GLC0
  v    = a->VGRD_P0_L7_GLC0
  
  dimuv   = dimsizes(ugrd)
  klev    = dimuv(0)
  nlat    = dimuv(1)
  mlon    = dimuv(2)

  do kl=10,10    ; 0,klev-1
     U = sin(rot)*vgrd(kl,:,:) + cos(rot)*ugrd(kl,:,:)
     V = cos(rot)*vgrd(kl,:,:) - sin(rot)*ugrd(kl,:,:)
  end do

;---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@cnFillPalette     = "gui_default" ; choose colormap
  cnres@cnLinesOn         = False    ; turn off contour lines
  
  cnres@gsnLeftString     = "REFC ("  + refc@units + ")"
  cnres@gsnRightString    = "UV winds ("  + u@units + ")"

  cnres@cnLevelSelectionMode = "ManualLevels"        ; set manual contour levels
  cnres@cnMinLevelValF       =  -8.                 ; set min contour level
  cnres@cnMaxLevelValF       =   36                  ; set max contour level
  cnres@cnLevelSpacingF      =   4                   ; set contour spacing

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

  cnres@mpMinLatF              = min(lat2d)-1
  cnres@mpMaxLatF              = max(lat2d)+1
  cnres@mpMinLonF              = min(lon2d)-1
  cnres@mpMaxLonF              = max(lon2d)+1
  
  
  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 do
end