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("data/K.nc","r")
K = f->K
f = addfile("data/third0.nc","r")
th = f->third0
;thm = dim_avg_n(th,(/1,2/))
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)
th1 = th(:,r2:nd(1)-1-r2,r2:nd(2)-1-r2)
nd := dimsizes(v1)

J = new(nd(0),integer)
fg2 = new(nd(0),float)

do i = 0,nd(0)-1
    dcp = ndtooned(v1(i,:,:))
    thi = ndtooned(th1(i,:,:))
    nn = dimsizes(dcp)
    do i1 = 0,nn-2
        do i2 = i1+1,nn-1
            if(dcp(i1).lt.dcp(i2)) then
                xx = dcp(i1)
                dcp(i1) = dcp(i2)
                dcp(i2) = xx
                delete(xx)
                xx = thi(i1)
                thi(i1) = thi(i2)
                thi(i2) = xx
                delete(xx)
            end if
        end do    
    end do 
    ;qsort(dcp)
    ;dcp = dcp(::-1)
    if(K(i).eq.0) then
        J(i) = 0
    end if
    if(K(i).eq.N) then
        J(i) = K(i)
    end if
    if(K(i).ne.0.and.K(i).ne.N) then
        do j = 0,K(i)-1
            ;print(dcp(j)*3600.)
            ;print((N-(j+1))*avg(dcp((j+1)::))/(j+1)*3600.)
            ;print(N*thi(j)/(j+1)*3600.)
            if(dcp(j).gt.-1.*((N-(j+1))*avg(dcp((j+1)::))+N*thi(j))/(j+1))
                J(i) = j+1
                fg2(i) = avg(dcp(J(i)::))
            end if
        end do
    end if
    delete(dcp)
end do
J!0 = "time"
fg2!0 = "time"

printVarSummary(J)
printVarSummary(fg2)
printMinMax(J,True)
printMinMax(fg2,True)

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

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