;*************************************************
; Modified version of wind_3.ncl from NCL webpage
;*************************************************
;This script computes the velocity potential and
;divergent winds on a fixed grid using
;spherical harmonics
;*************************************************
begin
;*************************************************
; open file and read in data: data are on a fixed grid
;*************************************************
  ufile = addfile("../uwnd/uwind_anom_ydaymean_1994-2008_JJA_200mb.nc","r")
  vfile = addfile("../vwnd/vwind_anom_ydaymean_1994-2008_JJA_200mb.nc","r")

  u = ufile->uwnd_anom(:,::-1,:)    ;sort coordinates using ::-1
  printVarSummary(u)
  v = vfile->vwnd_anom(:,::-1,:)
  printVarSummary(v)

;*************************************************
; calculate divergence on a fixed grid
;*************************************************
   div = uv2dvF_Wrap(u,v)
;*************************************************
; calculate velocity potential on a fixed grid
;*************************************************
  chi = uv2sfvpF(u,v)
  chi := chi(1,:,:,:)                   ; 1 for velocity potential
  chi = (/chi/1e6/)                     ; arbitrary scale
  chi@long_name = "velocity potential"
  chi@units     = "m/s" 
  printVarSummary(chi)
;*************************************************
; calculate divergent wind component 
; --
; note: the calculation uses a procedure, so memory 
; must be preallocated.
;*************************************************    
  ud    = new ( dimsizes(u), typeof(u), u@_FillValue )
  vd    = new ( dimsizes(v), typeof(v), v@_FillValue )

  dv2uvf(div,ud,vd) ; div  ==> divergent  wind components

  copy_VarCoords(u, ud ) 
  copy_VarCoords(u, vd ) 
  ud@long_name  = "Zonal Divergent Wind"
  ud@units      = u@units
  vd@long_name  = "Meridional Divergent Wind"
  vd@units      = v@units
;*************************************************
; plot results
;*************************************************    
  wks  = gsn_open_wks("png","wind_anom_ydaymean_1994-2008_JJA_200mb")             ; send graphics to PNG file
  cmap = read_colormap_file("nrl_sirkes")
  nc   = dimsizes(cmap(:,0))

  res                      = True

  res@cnFillOn             = True               ; color on
  res@cnLinesOn            = False              ; turn off contour lines

  res@gsnScalarContour     = True               ; vectors over contours
  res@cnFillPalette        = cmap(:nc-4,:)
  res@gsnAddCyclic         = False
  res@vcRefMagnitudeF      = 1.                 ; make vectors larger
  res@vcRefLengthF         = 0.050              ; ref vector length
  res@vcGlyphStyle         = "CurlyVector"      ; turn on curly vectors
  res@vcMinDistanceF       = 0.012              ; thin the vectors
  res@vcRefAnnoOrthogonalPosF  = -0.13         ; Move ref anno into plot

  res@mpLandFillColor       = "gray"            ; change continent color 

  res@tiMainString         = "Velocity Potential via Spherical Harmonics"
  res@gsnCenterString      = "Chi scaled by 1e6"
  res@gsnLeftString        = "Divergent Wind"
  
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF = -2.0
  res@cnMaxLevelValF = 2.0
  res@cnLevelSpacingF = 0.2

  res@gsnMajorLatSpacing = 5.0  
  res@gsnMajorLonSpacing = 10.0 
  res@mpMinLatF = -10.0
  res@mpMaxLatF = 60.0
  res@mpMinLonF = 60.0
  res@mpMaxLonF = 180.0
 
  plot=gsn_csm_vector_scalar_map(wks,ud(0,:,:),vd(0,:,:),chi(0,:,:),res)
end