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
r1 = 64
f = addfile("~/WORK4/CASPO/dcape2_jiwen/m"+r1+"/dcape2.nc","r")
dcape = f->dcape
f = addfile("~/WORK4/CASPO/mb_jiwen_no_shallow/m"+r1+"/data/dadt.nc","r")
dadt = f->dadt
f = addfile("~/WORK4/CASPO/mb_jiwen_no_shallow/m"+r1+"/data/jb.nc","r")
jb = f->jb
f = addfile("~/WORK3/jiwen/m"+r1+"/p.nc","r")
p = f->p

mb = dcape
nd = dimsizes(mb)
do i1 = 0,nd(0)-1
    do i2 = 0,nd(1)-1
        do i3 = 0,nd(2)-1
            if(dcape(i1,i2,i3)>0..and.dadt(i1,i2,i3)<0.and.p(jb(i1,i2,i3)-1,i1,i2,i3).ge.600.) then
                mb(i1,i2,i3) = -dcape(i1,i2,i3)*100./dadt(i1,i2,i3)/9.80616
                if(mb(i1,i2,i3).lt.0) then
                    mb(i1,i2,i3) = 0.
                end if    
            else
                mb(i1,i2,i3) = 0.
            end if
        end do
    end do
end do    
mb@description = "mb2"
printVarSummary(mb)
printMinMax(mb,True)

fn = "~/WORK4/CASPO/mb2_jiwen/m"+r1+"/mb2.nc"
system("rm -f "+fn)
filout = addfile(fn,"c")
filout->mb = mb
end
