; sigtopres.ncl
; Code originally from http://ccr.aos.wisc.edu/climate_modeling/modeling/data_processing/ncl/sigtopres.foam.web.ncl
;
; sigtopres.foam.web.ncl
; This program converts sigma levs to pressure levs using the NCL function "vinth2p".
; For info on vinth2p go to http://ngwww.ucar.edu/ngdoc/ng/ref/ncl/functions/vinth2p.html
; Note: You should NOT use this function to interpolate temp and geopot below ground (see web page above).

; S.R.  11.20.2000 
; Questions? Email srausche@students.wisc.edu

begin

; open file to write interpolated variable to (file should not already exist)
out = addfile("foam.11k.pressurevars.nc","c")

; open data file to get sigma lev data from
data = addfile("/grove3/iesuw/foam/ha.F1.11k.seae486605.nc", "r")

; choose levs you want to interpolate to (can have as many or as few as you like)
; here we are interpolating to 8 pressure levels 
lev = (/ 1000.,900.,850.,700.,500.,400.,300.,200. /)

; define some attributes of the variable "lev" so GrADS will be able to read the ouput file 
lev@units = "mb"
lev@long_name="pressure level"

; Define the variables and their dimensions.
; Variables have dimensions of (time,level,lat,lon)
; In this example, we are creating two pressure levels for seasonal data,
; so our variables are dimensioned (5(time),2(lev),40(lat),48(lon))  

omega = new((/5,8,40,48/),float)
u = new((/5,8,40,48/),float)
v = new((/5,8,40,48/),float)
zht = new((/5,8,40,48/),float)
q =  new((/5,8,40,48/),float)

; interpolate data using functions
; we use a "do" loop to calculate all seasons 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; arguments for vinth2p function are:

; function vinth2p(
; datai     : numeric, 3 or 4 dimensional data dimensioned time x lev x lat x lon or lev x lat x lon
; hbcofa[*] : numeric, Hybrid coefficients A, a single vector dimensioned the same as the lev dimension of datai,
; hbcofb[*] : numeric, Hybrid coefficients B, a single vector dimensioned the same as the lev dimension of datai
; plevo[*]  : numeric, An array of desired output pressure levels (mb)
; psfc      : numeric, 2 or 3 dimensional surface pressure data (Pa). These two or dimensions should 
;                      agree with the time, lat, and lon dimensions of datai 
; intyp[1]  : integer, Single scalar value to determine interpolation type 1 - LINEAR, 2 - LOG, 3 - LOG LOG 
; p0[1]     : numeric, Single scalar numeric value that is the surface reference pressure (mb)
; ii[1]     : integer, Not used. Set to 1,
; kxtrp[1]  : logical, If False no extrapolation is done when the pressure level is outside of the range of psfc
; )    
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


do i = 0,4
	u(i,:,:,:) = vinth2p(data->U(i,:,:,:), data->hyam, data->hybm, lev, data->PS(i,:,:), 1 ,1000, 1, False)
	omega(i,:,:,:) = vinth2p(data->OMEGA(i,:,:,:),data->hyam, data->hybm, lev, data->PS(i,:,:), 1 ,1000, 1, False)
        v(i,:,:,:) = vinth2p(data->V(i,:,:,:), data->hyam, data->hybm, lev, data->PS(i,:,:), 1 ,1000, 1, False)
        zht(i,:,:,:) = vinth2p(data->HT1(i,:,:,:), data->hyam, data->hybm, lev, data->PS(i,:,:), 1 ,1000, 1, False)
        q(i,:,:,:) = vinth2p(data->Q(i,:,:,:), data->hyam, data->hybm, lev, data->PS(i,:,:), 1 ,1000, 1, False)
end do


;get surface pressure (needed for grads VINT function) and divide by 100 to convert from pascals to mb
ps = data->PS(:,:,:)/100


; assign metadata to variables
; get time, lat, lon from orig data file, but use newly created lev variable
; need to do this so GrADS will read the file

q!0 = "time"
q&time = data->time
q!1 = "lev"
q&lev = lev
q!2 = "lat"
q&lat = data->lat
q!3 = "lon"
q&lon = data->lon

u!0 = "time"
u&time = data->time 
u!1 = "lev"
u&lev = lev 
u!2 = "lat"
u&lat = data->lat
u!3 = "lon"
u&lon = data->lon

v!0 = "time"
v&time = data->time
v!1 = "lev"
v&lev = lev
v!2 = "lat"
v&lat = data->lat
v!3 = "lon"
v&lon = data->lon

zht!0 = "time"
zht&time = data->time
zht!1 = "lev"
zht&lev = lev
zht!2 = "lat"
zht&lat = data->lat
zht!3 = "lon"
zht&lon = data->lon

omega!0 = "time"
omega&time = data->time
omega!1 = "lev"
omega&lev = lev
omega!2 = "lat"
omega&lat = data->lat
omega!3 = "lon"
omega&lon = data->lon

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

; write out data to file
; write variables
out->u = u
out->v = v
out->zht = zht
out->omega = omega
out->ps = ps

; write metadata 
out->lat = data->lat
out->lon = data->lon
out->time = data->time
out->lev = lev

end