;*************************************************
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/wrf/WRFUserARW.ncl"
;*************************************************
; add file example: a = addfile("/home/user/wrf/wrfout_d02_2013-02-12_12:00:00.nc","w") 
; DO NOT FORGET ".nc" !!!!!!!
;*************************************************
  a = addfile("wrfout_d02_2013-02-12_12:05:00.nc","w") 

  velocity_u = a->U           ; get u component wind units  m/s
  velocity_v = a->V           ; get v component wind units  m/s
  rain = a->QVAPOR            ; get the mass point dimensions
  mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file
;****************************************************
;convert wind components to mass points from staggered points
;****************************************************
  v=new((/mdims(0),mdims(1),mdims(2),mdims(3)/),"float")
  do j=0,mdims(2)-1
   v(:,:,j,:)=0.5*(velocity_v(:,:,j,:)+velocity_v(:,:,j+1,:))
  end do

  u=new((/mdims(0),mdims(1),mdims(2),mdims(3)/),"float")
  do i=0,mdims(3)-1
   u(:,:,:,i)=0.5*(velocity_u(:,:,:,i)+velocity_u(:,:,:,i+1))
  end do

;*************************************************
; calculate vorticity: Wrap version maintains meta data
;*************************************************
  vrt = uv2vrG_Wrap(u,v)
;*************************************************
; calculate rotational (non-divergent) wind components 
;*************************************************    
  ur    = new ( dimsizes(u), typeof(u), u@_FillValue )
  vr    = new ( dimsizes(v), typeof(v), v@_FillValue )

  vr2uvg(vrt, ur,vr)

;*************************************************
; assign attributes of the new variables
;*************************************************    
  ur@long_name  = "Non-divergent zonal wind"
  ur@short_name = "u wind"
  ur@units      = "m/s"

  vr@long_name  = "Non-divergent meridional wind"
  vr@short_name = "v wind"
  vr@units      = "m/s"

;*************************************************
;Write new variables
;*************************************************    
  filevardef(a, "UND", typeof(u), getvardims(rain))
  filevardef(a, "VND", typeof(u), getvardims(rain))
  a->UND = ur 
  a->VND = vr