;****************************************************
;This script calculates the velocity potential
;Input: U and V wind components
;****************************************************

begin

ufile = addfile("uwind200hPa_1975.nc","r") 
vfile = addfile("vwind200hPa_1975.nc","r") 

uwnd = ufile->uwnd(:,:,::-1,:)     ;sort coordinates usinf ::-1
printVarSummary(uwnd)

vwnd = vfile->vwnd(:,:,::-1,:)     ;sort coordinates using ::-1
printVarSummary(vwnd)

sfvp = uv2sfvpF(uwnd,vwnd)
sfvp = (/sfvp/1e6/)

copy_VarMeta(uwnd,vp)

vp@long_name = "velocity potential"
vp@units     = "m/s" 

ncdf = addfile("test.nc" ,"c")  ; open output netCDF file

fAtt               = True            ; assign file attributes
fAtt@title         = "Velocity Potential"
fAtt@source_file   =  "https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html"
fAtt@Conventions   = "None"
fAtt@creation_date = systemfunc ("date")
fileattdef( ncdf, fAtt )            ; copy file attributes

ncdf->vp = vp
end