; ===============================================
; plot vectors EP flux
; ===============================================
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"
;---------------------------------
begin

;---read in---
ncdfin=addfile("test_vector_EP_flux.nc","r")
xT20=ncdfin->xT20 
zT20=ncdfin->zT20
delete(ncdfin)

;==================================================
; plot stuff
;-----------------
plots1=new( 1,"graphic")

 wks    = gsn_open_wks ("ps", "test_vector_EP_flux"  )
gsn_define_colormap(wks,"ncl_default")

 res            = True
; res@gsnDraw    = False         ; Don't draw the plot or advance the
; res@gsnFrame   = False         ; Don't advance frame

;--vector stuff---
res@vcMinDistanceF  = 0.03
;res@vcRefMagnitudeF = 1.
res@vcRefMagnitudeF = 8.79354e+07   ; setting this to the upper 1 percentile
res@vcMaxMagnitudeF = res@vcRefMagnitudeF * 12.0   ; this is about the biggest the maxium vector can be if ref length is 0.1
res@vcRefLengthF    = 0.1              ; define length of vec ref
res@vcMinFracLengthF = 0.1  ; still needed to get most of the vectors to have visible length
;res@vcGlyphStyle = "curlyvector"

printMinMax(xT20,0)
printMinMax(zT20,0)
res@tiMainString= "vector EP Flux"
res@vcMapDirection = False
x = xT20(lev|:,lat|:) 
z = zT20(lev|:,lat|:)
;x = x * 1e-6
;z = z * 1e-6
s = sqrt (x * x + z * z)
print (num(s .gt. 4.21171e+07))
print (num(s .le. 4.21171e+07))
opt = True
opt@PrintStat=True
st = stat_dispersion(s,opt)
xt = stat_dispersion(x,opt)
zt = stat_dispersion(z,opt)

printMinMax(s,0)
plots1     = gsn_csm_vector (wks, x,z, res )

getvalues plots1 ;     
  "vcRefLengthF" : ref_len
  "vcRefMagnitudeF" : ref_mag
  "vcMinDistanceF" : ref_dist
end getvalues
print("ref len is " + ref_len)
print("ref mag is " + ref_mag)
print("ref dist is " + ref_dist)

;pres = True
;pres@gsnPaperOrientation = "portrait"
;gsn_panel(wks,plots1,(/1,1/),pres)
end