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
r1 = 512/r
g = 9.8
var = new((/160,r1,r1/),float)
var = 0.

f = addfile("data/wpthetavp_m.nc","r")
wth0 = f->wpthetavp
f = addfile("data/thetav.nc","r")
th = f->thetav
;printVarSummary(th)
delete(f)

f = addfile("~/WORK4/CASPO/mb_jiwen_no_shallow/m"+r+"/data/jb.nc","r")
jb = f->jb
f = addfile("~/WORK4/CASPO/mb_jiwen_no_shallow/m"+r+"/data/jt.nc","r")
jt = f->jt
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(z)
do i = 0,nd(1)-1
    do i1 = 0,nd(2)-1
        do i2 = 0,nd(3)-1
            do k = jt(i,i1,i2)-1,jb(i,i1,i2)-1
                wth = -1.*(wth0(k,i)-wth0(k+1,i))
                var(i,i1,i2) = var(i,i1,i2)+g*wth/th(k,i,i1,i2)
                delete(wth)
            end do
        end do
    end do 
end do
var!0 = "time"
var!1 = "y"
var!2 = "x"
printVarSummary(var)
printMinMax(var*3600.,True)

fn = "data/third0.nc"
system("rm -f "+fn)
filout = addfile(fn,"c")
filout->third0 = var
end
