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

  a = addfile("/home/swati/work/analysis/2_cities/Patna/14aug2014/wrf_out_data/3rd_domain/13_aug.nc","r")

; The specific pressure levels that we want the data interpolated to.
  pressure = (/ 850 /)   ; pressure levels to plot

; First get the variables we will need        

    u  = wrf_user_getvar(a,"ua",-1)        ; u averaged to mass points
    v  = wrf_user_getvar(a,"va",-1)        ; v averaged to mass points
    p  = wrf_user_getvar(a,"pressure",-1) ; used to mention vertical coordinate to interplate to

    uu  = wrf_user_intrp3d( u,p,"h",pressure,0.,False)
    vv  = wrf_user_intrp3d( v,p,"h",pressure,0.,False)

    u_plane  = dim_avg_Wrap(uu(south_north|:,west_east|:, Time|:) )
    v_plane  = dim_avg_Wrap(vv(south_north|:,west_east|:, Time|:) )

    spd     = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec
    spd@description = "Wind Speed"
    spd@units = "m/s"
    u_plane@units = "m/s"
    v_plane@units = "m/s"

;************************************************************************************

 u1= wrf_user_getvar(a,"U10",0)

 copy_VarCoords(u1,u_plane)                         ; copy coord vars to speed
 copy_VarCoords(u1,v_plane)                         ; copy coord vars to speed
 copy_VarCoords(u1,spd)                         ; copy coord vars to speed

;************************************************************************************
;printVarSummary(u_plane)
;exit

;---Set some basic plot options

  wks = gsn_open_wks("png","wind_at_850_using_gsn_function")

  res               = True
  res@gsnDraw  = False                          ; don't draw
  res@gsnFrame = False                          ; don't advance frame
  res@cnFillOn      = True    			; color plot desired
  res@cnLinesOn     = False         		; no contour lines
  res@gsnAddCyclic  = False
  res@cnLineLabelsOn   = False
  res@lbLabelBarOn = False 
  res@cnFillPalette =  "rainbow"                    ; "WhiteBlueGreenYellowRed"

  lat = wrf_user_getvar(a,"XLAT",0)   ; latitude/longitude
  lon = wrf_user_getvar(a,"XLONG",0)  ; required for plotting  

 res@sfXArray = lon
  res@sfYArray = lat

  res@mpMinLatF     = min(lat)
  res@mpMaxLatF     = max(lat)
  res@mpMinLonF     = min(lon)
  res@mpMaxLonF     = max(lon)

  res@pmTickMarkDisplayMode = "Always"   ; nicer tickmarks
  res@mpProjection  = "CylindricalEquidistant"    ; The default
  res@mpDataBaseVersion     = "MediumRes"       ; better and more map outlines
  res@mpDataSetName         = "Earth..4"
  res@mpOutlineBoundarySets = "AllBoundaries"
  res@mpOutlineOn           = True 

;---Necessary for contours to be overlaid correctly on WRF projection
;  res@tfDoNDCOverlay = True

  
  res@gsnScalarContour     = True               ; contours desired
  res@cnFillDrawOrder      = "PreDraw"

  res@vcRefMagnitudeF           = 4.0             ; define vector ref mag
  res@vcRefLengthF              = 0.045           ; define length of vec ref
 ; res@vcRefAnnoOrthogonalPosF   = -1.0            ; move ref vector
 ; res@vcRefAnnoArrowLineColor   = "black"         ; change ref vector color
 ; res@vcRefAnnoArrowUseVecColor = False           ; don't use vec color for ref

  res@vcGlyphStyle            = "CurlyVector"     ; turn on curly vectors
  res@vcLineArrowColor        = "black"           ; change vector color
  res@vcLineArrowThicknessF   = 2.0               ; change vector thickness
  res@vcVectorDrawOrder       = "PostDraw"        ; draw vectors last

  plot=gsn_csm_vector_scalar_map_ce(wks,u_plane,v_plane,spd,res)