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
N = 384
r = 64
N = N/r
N = N*N
print("N = "+N)

f = addfile("~/WORK4/CASPO/mb_jiwen_no_shallow/m"+r+"/data/dcape.nc","r")
v0 = f->dcape
nd = dimsizes(v0)
r2 = 64/r
v1 = v0(:,r2:nd(1)-1-r2,r2:nd(2)-1-r2)
nd := dimsizes(v1)
v2 = new(nd(0),float)
K  = new(nd(0),integer)
v2 = 0.
do i = 0,nd(0)-1
    aa = 0.
    do i1 = 0,nd(1)-1
        do i2 = 0,nd(2)-1
            if(v1(i,i1,i2).le.0) then
                v2(i) = v2(i)+v1(i,i1,i2)
                aa = aa+1
            end if
        end do
    end do
    if(aa.eq.0) then
        v2(i) = 0.
        K(i)  = N
    else
        if(aa.eq.N) then
            v2(i) = 0.
            K(i)  = 0
        else
            v2(i) = v2(i)/aa
            K(i)  = N-floattointeger(aa)
            ;v2(i) = v2(i)*int2flt((N-K(i))/K(i))
        end if     
    end if    
end do
v2!0 = "time"
K!0  = "time"
printVarSummary(v2)
printVarSummary(K)
;printMinMax(v2,True)
;printMinMax(K,True)

fn = "data/fgn.nc"
system("rm -f "+fn)
filout = addfile(fn,"c")
filout->fgn = v2

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