; velocity potential.ncl
; Program to compute streamfunction, velocity potential, divergence, and vorticity
; sr 11.20.2000

begin

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

; open file to write interpolated variable to (file should not already exist)
out = addfile("GHA.2005.nc","c")

; open data file to get sigma lev data from
data = addfile("ATMdaily.2005.nc", "r")

;set dims of data
latS = 11
latN = 22
lonL = 19
lonR = 79
timeS = 60
timeE = 151

;set missing value
uvmsg = 1e+36

;create output variables with proper dimensions
dvx   = new ( (/timeS,timeE,latS,latN,lonL,lonR /), float, uvmsg )  ; divergence (same as dv)
vrx   = new ( (/timeS,timeE,latS,latN,lonL,lonR /), float, uvmsg )  ; vorticity  (same as vr)
sf    = new ( (/timeS,timeE,latS,latN,lonL,lonR /), float, uvmsg )  ; stream function
vp    = new ( (/timeS,timeE,latS,latN,lonL,lonR /), float, uvmsg )  ; velocity potential

; Run function for all seasons (cycle through time using a do-loop)
; Here we are calculating al vars for the lowest model level (near the ground), lev=11
; Info on these functions can be found at 
; http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/uv2sfvp.html (sf and vp)
; http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/uv2vrdv.html (div and vort)

do i = 0,4
	uv2vrdvg(data->U(i,11,:,:),data->V(i,11,:,:),vrx(i,:,:),dvx(i,:,:))   ; u,v ==> div and vort 
	uv2sfvpg(data->U(i,11,:,:),data->V(i,11,:,:),sf(i,:,:),vp(i,:,:))     ; u,v ==> stream function + velocity pot
end do

; Assign dimensions to output vars so GrADS can read the files
dvx!0 = "time"
dvx&time = data->time
dvx!1 = "lat"
dvx&lat = data->lat
dvx!2 = "lon"
dvx&lon = data->lon

vrx!0 = "time"
vrx&time = data->time
vrx!1 = "lat"
vrx&lat = data->lat
vrx!2 = "lon"
vrx&lon = data->lon

sf!0 = "time"
sf&time = data->time
sf!1 = "lat"
sf&lat = data->lat
sf!2 = "lon"
sf&lon = data->lon

vp!0 = "time"
vp&time = data->time
vp!1 = "lat"
vp&lat = data->lat
vp!2 = "lon"
vp&lon = data->lon

; Output data to file
; variables
out->dvx = dvx
out->vrx = vrx
out->sf = sf
out->vp = vp
; dimensions
out->lat = data->lat
out->lon = data->lon
out->time = data->time



end