;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  filename = "/home/gaurav/wrf_exp/wrfout_d02_2014-10-26_00:00:00-fnl-myjbmj.nc"
  a        = addfile(filename,"r")
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  it  = 07
  ;u = wrf_user_getvar(a,"U10",it)   ; U component of the wind at 10m
  u = wrf_user_getvar(a,"U",it)      ; U component of the wind at 10m
  printVarSummary(u)
  ;v = wrf_user_getvar(a,"V10",it)   ; V component of the wind at 10m
  v = wrf_user_getvar(a,"V",it)      ; V component of the wind at 10m
  printVarSummary(v)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 ;u10 = u10*1.94386                    
 ;v10 = v10*1.94386
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  u@units       = "m/s"
  v@units       = "m/s"
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  wks = gsn_open_wks("png","wind")
  printVarSummary(wks)
;---Set common resources for all plots
  res                = True
  res@gsnFrame       = False
  res@gsnDraw        = True
  res@gsnLeftString  = ""
  res@gsnRightString = ""
  ;res@vcMonoLineArrowColor = False
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;---Necessary for contours to be overlaid correctly on WRF projection
  res@tfDoNDCOverlay = True
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;   
;---Wind vector plot
    cmap = read_colormap_file("BlAqGrYeOrReVi200")
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  
  vec_res                  = res
  vec_res@vcFillArrowsOn   = True
  vec_res@vcMinDistanceF   = 0.02
  vec_res@vcRefLengthF     = 0.025 
  vec_res@vcRefMagnitudeF  = 15
  vec_res@vcMinFracLengthF = 0.2
  vec_res@vcLevelPalette   = cmap(6:193,:)
  vec_res@vcWindBarbColor  = "blue"
  vec_res@vcRefAnnoOn      = True
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  opts                     = res
  map_res                  = True
  pltres                   = True
  pressure_level          = (/ 850./)   ; pressure levels to plot
  nlevels                  = dimsizes(pressure_level)     ; number of pressure levels
  p                        = wrf_user_getvar(a, "pressure",it) 
  do level                 = 0,nlevels-1
  pressure                 = pressure_level(level)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    if ( pressure .eq. 850 ) then   
          ;contour                         = wrf_contour(a,wks,u(level,:,:),opts)
          contour                          = wrf_vector(a,wks,u,v,vec_res)
          plot                             = wrf_map_overlays(a,wks,(/contour/),pltres,map_res)
          draw(contour)
 ;        frame(wks)
        end if
    end do
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  
;---Map plot
  ;map_res                      = True
  map_res@mpLandFillColor      = "tan"
  map_res@mpOceanFillColor     = "cornflowerblue"
  map_res@gsnFrame             = True
  map_res@gsnDraw              = True
  map_res@tmXBLabelFontHeightF = 0.025                 ; resize tick labels
  map_res@tmYLLabelFontHeightF = 0.025  
  map_res@mpGeophysicalLineThicknessF = 5.
  map_res@mpOutlineOn           = True
;---Set map resources based on projection on WRF output file
  map_res = wrf_map_resources(a,map_res)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  map = gsn_csm_map(wks,map_res)
  vector = gsn_csm_vector(wks,u(:,:,:),v(:,:,:),vec_res)
  overlay(map,vector)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  draw(map)   ; This will draw all overlaid plots and the map
  frame(wks)
end