;********************************************************

; histogram normalized by rain and compute correleration

;********************************************************

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"

procedure pminmax(data:numeric,name:string)

begin

  print ("min/max " + name + " = " + min(data) + "/" + max(data))

  if(isatt(data,"units")) then

    print (name + " units = " + data@units)

  end if

end

;

; Main code.

;

begin

;nclass = 18

 nclass = 20

;************************************************

; read in data: observed

;************************************************

 diri1  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"

;fili1  = "land_class_T42.nc"

 fili1  = "land_class_T42_new.nc"

 fili2  = "LAI_2000-2005_mean_T42.nc"

 data_file_ob1 = addfile(diri1+fili1,"r")

 data_file_ob2 = addfile(diri1+fili2,"r")

 RAIN1 = tofloat(data_file_ob1->LAND_CLASS)               

 NPP1  = data_file_ob2->LAI      

;************************************************

; read in data: model       

;************************************************

 diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

;fili3  = "i01.03cn_1545-1569_ANN_climo.nc"

 fili3  = "i01.04casa_1605-1629_ANN_climo.nc"

 data_file_model = addfile(diri2+fili3,"r")

 NPP2  = data_file_model->TLAI      

;************************************************

; print min/max and unit

;************************************************

  pminmax(RAIN1,"RAIN1")

  pminmax(NPP1,"NPP1")

  pminmax(NPP2,"NPP2")

  RAIN1_1D = ndtooned(RAIN1)

  NPP1_1D  = ndtooned(NPP1)

  NPP2_1D  = ndtooned(NPP2)

;

; Calculate some "nice" bins for binning the data in equally spaced

; ranges.

;

; nbins       = nclass + 1         ; Number of bins to use.

; nicevals    = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,False)

; nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)

; range       = fspan(nicevals(0),nicevals(1),nvals)

  nclassn     = nclass + 1

  range       = fspan(0,nclassn-1,nclassn)

; print (nicevals)

; print (nvals)

  print (range)

; exit

;

; Use this range information to grab all the values in a

; particular range, and then take an average.

;

  nr      = dimsizes(range)

  nx      = nr-1

  xvalues     = new((/2,nx/),typeof(RAIN1_1D))

  xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.

  dx           = xvalues(0,1) - xvalues(0,0)       ; range width

  dx4          = dx/4                              ; 1/4 of the range

  xvalues(1,:) = xvalues(0,:) - dx/5.

  yvalues      = new((/2,nx/),typeof(RAIN1_1D))

  mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

  mx_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

  do nd=0,1

;

; See if we are doing model or observational data.

;

    if(nd.eq.0) then

      data     = RAIN1_1D

      npp_data = NPP1_1D

    else

      data     = RAIN1_1D

      npp_data = NPP2_1D

    end if

;

; Loop through each range and check for values.

;

    do i=0,nr-2

      if (i.ne.(nr-2)) then

         print("")

         print("In range ["+range(i)+","+range(i+1)+")")

        idx = ind((range(i).le.data).and.(data.lt.range(i+1)))

      else

         print("")

         print("In range ["+range(i)+",)")

        idx = ind(range(i).le.data)

      end if

;

; Calculate average, and get min and max.

;

      if(.not.any(ismissing(idx))) then

        yvalues(nd,i)    = avg(npp_data(idx))

        mn_yvalues(nd,i) = min(npp_data(idx))

        mx_yvalues(nd,i) = max(npp_data(idx))

        count = dimsizes(idx)

      else

        count            = 0

        yvalues(nd,i)    = yvalues@_FillValue

        mn_yvalues(nd,i) = yvalues@_FillValue

        mx_yvalues(nd,i) = yvalues@_FillValue

      end if

;

; Print out information.

;

       print(nd + ": " + count + " points, avg = " + yvalues(nd,i))

       print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))

;

; Clean up for next time in loop.

;

      delete(idx)

    end do

    delete(data)

    delete(npp_data)

  end do

;

; Start the graphics.

;

 u = yvalues(0,:)

 v = yvalues(1,:)

 print (u)

 print (v)

 good = ind(.not.ismissing(u) .and. .not.ismissing(v))

 uu = u(good)

 vv = v(good)

 nz = dimsizes(uu)

 print (nz)

 ccr = esccr(uu,vv,0)

 print (ccr)

;new eq

 bias = sum(abs(vv-uu)/(vv+uu))

 M    = (1.- (bias/nz))*5.

 print (bias)

 print (M)

end