;create slp
;03/05/02

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/csm/shea_util.ncl"

begin

;open file to write variable to
out = addfile("ha.F1_5.preindctl.30yrens477506.psl.nc","c")

;open data file
f = addfile("ha.F1_5.preindctl.30yrens477506.nc","r")

;set some constants
r = 287.04
g = 9.80616
gmoist = 6.5e-3
econst = g/(r*gmoist)
pconst = gmoist/g

;get necessary variables
hyam = f->hyam
hybm = f->hybm
PSFC = f->PS
PHIS = f->PHIS
T    = f->T

;get dimension sizes
ntimes = dimsizes(f->time)
nlevs = dimsizes(hyam) - 1
nlons = dimsizes(f->lon)
nlats = dimsizes(f->lat)

;create new vars
presfac = new((/nlons/),float)
tsurf = new((/nlons/),float)
psl = new((/ntimes,nlats,nlons/),float)

;enter do loops to calc psl
do lon  = 0, nlons - 1
do lat  = 0, nlats - 1
do itim = 0, ntimes - 1

        pnlev = hyam(nlevs)*1000.*100. + hybm(nlevs)*PSFC(itim,lat,lon)
        pnlevm = hyam(nlevs-1)*1000.*100. + hybm(nlevs-1)*PSFC(itim,lat,lon)

        presfac(lon) = (log(PSFC(itim,lat,lon)/pnlev)) / (( log(pnlev/pnlevm) ))
        tsurf(lon) = T(itim,nlevs,lat,lon) + presfac(lon)*(T(itim,nlevs-1,lat,lon) - T(itim,nlevs,lat,lon) )
        psl(itim,lat,lon) = (0.01 * PSFC(itim,lat,lon)) * ((1.0+pconst*PHIS(itim,lat,lon)/tsurf(lon))^econst)

end do
end do
end do
;end do loops

copy_VarCoords(f->T(:,0,:,:),psl)

out->lat = f->lat
out->lon = f->lon
out->time = f->time
out->psl = psl

end