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
r = 64
f = addfile("~/WORK3/jiwen/m"+r+"/tk.nc","r")
T = f->tk
f = addfile("~/WORK3/jiwen/m"+r+"/qv.nc","r")
q = f->qv
T = T*(1.+0.608*q)
f = addfile("~/WORK3/jiwen/m"+r+"/z.nc","r")
z = f->z
f = addfile("~/WORK3/jiwen/m"+r+"/zf.nc","r")
zf = f->z

nd = dimsizes(zf)
Tf = new((/nd(0),nd(1),nd(2),nd(3)/),float)
do i1 = 0,nd(2)-1
    do i2 = 0,nd(3)-1
        do k = 0,nd(1)-1
            Tf(:,k,i1,i2) = wrf_interp_1d(T(:,k,i1,i2),z(:,k,i1,i2),zf(:,k,i1,i2))
        end do    
    end do
end do
Tf(0,:,:,:) = 2.*Tf(1,:,:,:)-Tf(2,:,:,:)
Tf(nd(0)-1,:,:,:) = 2.*Tf(nd(0)-2,:,:,:)-Tf(nd(0)-3,:,:,:)
    
f = addfile("~/WORK3/jiwen/m"+r+"/pf.nc","r")
pf = f->p
var = pot_temp(pf,Tf,-1,False)
var!1 = "time"
printVarSummary(var)
printMinMax(var,True)
    
fn = "data/thetavf.nc"
system("rm -f "+fn)
filout = addfile(fn,"c")
filout->thetav = var
end
