;*************************************************
;This script computes for the daily vorticity from
;a global data set
;Created: April 26,2017
;Based on NCL example vort_1.ncl
;*************************************************

begin
;************************************************
; variable and file handling
;************************************************
  ufile  = addfile("jan16_uwnd_2017.nc","r")       ; open netcdf file
  u   = short2flt(ufile->uwnd(:,{925},:,:))                             ; pull u off file
  lat = ufile->lat
  lon = ufile->lon
  vfile  = addfile("jan16_vwnd_2017.nc","r")
  v   = short2flt(vfile->vwnd(:,{925},:,:))                             ; pull v off file
   
  u = u(:,::-1,:)
  printVarSummary(u)
  v = v(:,::-1,:)
  printVarSummary(v)

  printVarSummary(lat)
  printVarSummary(lon)

;************************************************
; calculate vorticity on a Fixed Grid
;************************************************         
  vr = uv2vr_cfd (u,v,lat,lon, 2)
  copy_VarMeta(u,vr)
  vr@long_name = "vorticity"
  vr@units     = "per second"


;***********************************************
;Perform centered finite difference

dUdY = center_finite_diff_n(u,lat, False, 0, 1)
dVdX = center_finite_diff_n(v,lon, False, 0, 1)

;***********************************************
; save to a netcdf file
;***********************************************
  diro = "./"
  filo = "vort_jan2017_925hPa.nc"
  ptho = diro+filo
  system("/bin/rm -f "+ptho)
  ncdf = addfile(ptho,"c")

  fAtt = True
  fAtt@title         = "Vorticity "
  fAtt@source_file   = "daily winds from NCEP-NCAR Reanaysis"
  fAtt@Conventions   = "None"
  fAtt@creation_date = systemfunc("date")
  fileattdef(ncdf,fAtt)            ; copy file attributes
  filedimdef(ncdf,"time",-1,True)
  ncdf->vr = vr
  ncdf->dUdY = dUdY
  ncdf->dVdX = dVdX
end