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

 ; histogram normalized by rain and compute correleration

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

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"

 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 = 20

  plot_type     = "ps"

  plot_type_new = "png"

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

 ; read data: model       

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

  co2_i = 283.1878

  co2_f = 364.1252

  model_grid = "T42"

 ;model_name_i = "i01.07cn"

 ;model_name_f = "i01.10cn"

  model_name_i = "i01.07casa"

  model_name_f = "i01.10casa"

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

  film_i = model_name_i + "_1990-2004_ANN_climo.nc"

  film_f = model_name_f + "_1990-2004_ANN_climo.nc"

  fm_i   = addfile (dirm+film_i,"r")

  fm_f   = addfile (dirm+film_f,"r")

  npp_i  = fm_i->NPP

  npp_f  = fm_f->NPP

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

 ; read data: observed

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

  ob_name = "MODIS MOD 15A2 2000-2005"

  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"

  filo  = "land_class_"+model_grid+".nc"

  fo = addfile(diro+filo,"r")

  classob    = tofloat(fo->LAND_CLASS)

  class_name = (/"Water Bodies" \

                ,"Evergreen Needleleaf Forests" \

                ,"Evergreen Broadleaf Forests" \

                ,"Deciduous Needleleaf Forest" \

                ,"Deciduous Broadleaf Forests" \

                ,"Mixed Forests" \                      

                ,"Closed Bushlands" \                   

                ,"Open Bushlands" \                     

                ,"Woody Savannas (S. Hem.)" \           

                ,"Savannas (S. Hem.)" \                 

                ,"Grasslands" \                         

                ,"Permanent Wetlands" \                 

                ,"Croplands" \                         

                ,"Urban and Built-Up" \                 

                ,"Cropland/Natural Vegetation Mosaic" \ 

                ,"Permanent Snow and Ice" \             

                ,"Barren or Sparsely Vegetated" \       

                ,"Unclassified" \                       

                ,"Woody Savannas (N. Hem.)" \           

                ,"Savannas (N. Hem.)" \                

                /)  

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

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

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

   nclassn     = nclass + 1

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

 ; print (range)

 ; 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/),float)

   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.

 ; get data

   DATA11_1D = ndtooned(classob)

   DATA12_1D = ndtooned(npp_i)

   DATA22_1D = ndtooned(npp_f)

   yvalues      = new((/2,nx/),float)

   mn_yvalues   = new((/2,nx/),float)

   mx_yvalues   = new((/2,nx/),float)

   do nd=0,1

 ; See if we are doing model or observational data.

     if(nd.eq.0) then

       data_ob  = DATA11_1D

       data_mod = DATA12_1D

     else

       data_ob  = DATA11_1D

       data_mod = DATA22_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((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))

       else

 ;        print("")

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

          idx = ind(data_ob.ge.range(i))

       end if

 ; Calculate average, and get min and max.

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

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

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

         mx_yvalues(nd,i) = max(data_mod(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(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_ob)

     delete(data_mod)

   end do

 ;============================

 ;compute beta

 ;============================

  u = yvalues(0,:)

  v = yvalues(1,:)

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

  uu = u(good)

  vv = v(good)

  ww = class_name(good)

  n_biome = dimsizes(uu)

  beta_biome = new((/n_biome/),float)

  beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)

  beta_biome_avg = avg(beta_biome)

  print("class/beta:  " + ww + "/" + beta_biome)

  print (beta_biome_avg)

 end