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
r2 = 64/r

f = addfile("~/WORK3/jiwen/m"+r+"/wf.nc","r")
wf = f->w
printVarSummary(wf)
f = addfile("data/thetavf.nc","r")
thvf = f->thetav

nd = dimsizes(wf)
wf1 = wf(:,:,r2:nd(2)-1-r2,r2:nd(3)-1-r2)
thvf1 = thvf(:,:,r2:nd(2)-1-r2,r2:nd(3)-1-r2)

wm = dim_avg_n(wf1,(/2,3/))
tm = dim_avg_n(thvf1,(/2,3/))

nd := dimsizes(wf1)

var = new((/nd(0),nd(1)/),float)

do i = 0,nd(1)-1
    do k = 0,nd(0)-1
        wp = wf1(k,i,:,:)-wm(k,i)
        tp = thvf1(k,i,:,:)-tm(k,i)

        vp = wp*tp
        delete(wp)
        delete(tp)
        var(k,i) = dim_avg_n(dim_avg_n(vp,0),0)
        delete(vp)
    end do
end do
var!0 = "z(top_bottom)"
var!1 = "time"
printVarSummary(var)
printMinMax(var,True)

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