;********************************************
; vector_5.ncl
;********************************************
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"
begin

; file handling
;************************************************
   fn  = "atmos.nc" ; define filename
   in  = addfile(fn,"r")                         ; open netcdf file
;************************************************
; read needed variables from file
;************************************************
   T = in->T                                    ; select variable to ave
   W = in->OMEGA
   V = in->V
   P0mb = 1000.
   hyam = in->hyam                              ; get a coefficiants
   hybm = in->hybm                              ; get b coefficiants
   PS   = in->PS                                ; get pressure
;************************************************
; define other arguments required by vinth2p
;************************************************
   interp = 2 
   tlev   =  2
   pnew   = (/ 900,850,800,750,700,650,600,550,500,450,400,350,300,250,200/)
   pnew@units = "mb"          
;************************************************
; interpolate to pressure levels on pressure levels
;************************************************
  t = vinth2p(T,hyam,hybm,pnew,PS,interp,P0mb,tlev,False)
  v = vinth2p(V,hyam,hybm,pnew,PS,interp,P0mb,tlev,False)
  w = vinth2p(W,hyam,hybm,pnew,PS,interp,P0mb,tlev,False)
;************************************************
; Omega is significantly smaller than v, so we will
; scale it so that some vertical motion is visible
;************************************************
 wAve   = avg(w(0,:,:,{170}))           ; used for scaling
 vAve   = avg(v(0,:,:,{170}))
 scale  = fabs(vAve/wAve)
 wscale = w*scale                       ; now scale

 copy_VarCoords(w, wscale)              ; copy coordinate variables 
;***********************************************
; create plot
;***********************************************
 wks   = gsn_open_wks ("ps", "vector" )       ; open workstation
 gsn_define_colormap(wks,"BlAqGrYeOrRevi200")  ; choose color map

 res                 = True                     ; plot mods desired
 res@tiMainString    = "Pressure/Height Vector" ; title

 res@cnLineLabelsOn  = False              ; turn off line labels
 res@cnFillOn        = True               ; turn on color fill
 res@lbLabelStride   = 2                  ; every other color
 
 res@gsnSpreadColors = True               ; use full range of color map

 res@vcRefMagnitudeF = 3.0                ; define vector ref mag
 res@vcRefLengthF    = 0.045              ; define length of vec ref
 res@vcGlyphStyle    = "CurlyVector"      ; turn on curley vectors
 res@vcMinDistanceF  = 0.01               ; thin out vectors
;*****************************************************
; draw plot from pole to pole at 170E
;*****************************************************
 plot  = gsn_csm_pres_hgt_vector(wks,t(0,:,:,{170}),v(0,:,:,{170}),\
                                wscale(0,:,:,{170}),res )  

end