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

 ; 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_ensemble_T42.nc"

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

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

  RAIN1 = tofloat(data_file_ob1->LAND_CLASS)

   z     = data_file_ob2->LAI

   y     = z(0,:,:)

   y@long_name = "Leaf Area Index Max"

   s     = z(:,0,0)

   dsizes_z = dimsizes(z)

   ntime    = dsizes_z(0)

   nlat     = dsizes_z(1)

   nlon     = dsizes_z(2)

   do j = 0,nlat-1

   do i = 0,nlon-1

      s      = z(:,j,i) 

      y(j,i) = max(s)

   end do

   end do

   print (min(y)+"/"+max(y))

   NPP1  = y

   delete (z)

   delete (s)

   delete (y)      

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

 ; read in data: model       

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

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

  fili3  = "i01.03cn_1545-1569_MONS_climo.nc"

 ;fili3  = "i01.04casa_1605-1629_MONS_climo.nc"

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

   z     = data_file_model->TLAI

   y     = z(0,:,:)

   y@long_name = "Leaf Area Index Max"

   s     = z(:,0,0)

   dsizes_z = dimsizes(z)

   ntime    = dsizes_z(0)

   nlat     = dsizes_z(1)

   nlon     = dsizes_z(2)

   do j = 0,nlat-1

   do i = 0,nlon-1

      s      = z(:,j,i) 

      y(j,i) = max(s)

   end do

   end do

   print (min(y)+"/"+max(y))

   NPP2  = y

   delete (z)

   delete (s)

   delete (y)            

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

 ; 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)

   nclass      = nclass + 1

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

 ; 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.

 ;

   wks = gsn_open_wks("ps","xy")

   res             = True

   res@gsnMaximize = True

   res@gsnDraw     = False

   res@gsnFrame    = False

   res@xyMarkLineMode = "Markers"

   res@xyMarkerSizeF   = 0.014

   res@xyMarker       = 16

   res@xyMarkerColors = (/"Brown","Blue"/)

 ; res@trYMinF        = min(mn_yvalues) - 10.

 ; res@trYMaxF        = max(mx_yvalues) + 10.

   res@trYMinF        = min(mn_yvalues) - 2

   res@trYMaxF        = max(mx_yvalues) + 4

   res@tiMainString  = "Observed vs i01.03cn"

 ; res@tiMainString  = "Observed vs i01.04casa"

   res@tiYAxisString = "Max LAI (Leaf Area Index)"

   res@tiXAxisString = "Land Cover Type"

 ;

 ; Add a boxed legend using the more simple method, which won't have

 ; vertical lines going through the markers.

 ;

   res@pmLegendDisplayMode    = "Always"

 ; res@pmLegendWidthF         = 0.1

   res@pmLegendWidthF         = 0.08

   res@pmLegendHeightF        = 0.05

   res@pmLegendOrthogonalPosF = -1.17

 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

 ; res@pmLegendParallelPosF   =  0.18

   res@pmLegendParallelPosF   =  0.23  ;(rightward)

 ;  res@lgPerimOn             = False

   res@lgLabelFontHeightF     = 0.015

   res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/)

 ; res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/)

   xy = gsn_csm_xy(wks,xvalues,yvalues,res)

   max_bar = new((/2,nx/),graphic)

   min_bar = new((/2,nx/),graphic)

   max_cap = new((/2,nx/),graphic)

   min_cap = new((/2,nx/),graphic)

   lnres = True

   line_colors = (/"brown","blue"/)

   do nd=0,1

     lnres@gsLineColor = line_colors(nd)

     do i=0,nx-1

       if(.not.ismissing(mn_yvalues(nd,i)).and. \

          .not.ismissing(mx_yvalues(nd,i))) then

 ;

 ; Attach the vertical bar, both above and below the marker.

 ;

         x1 = xvalues(nd,i)

         y1 = yvalues(nd,i)

         y2 = mn_yvalues(nd,i)

         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

         y2 = mx_yvalues(nd,i)

         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

 ;

 ; Attach the horizontal cap line, both above and below the marker.

 ;

         x1 = xvalues(nd,i) - dx4

         x2 = xvalues(nd,i) + dx4

         y1 = mn_yvalues(nd,i)

         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

         y1 = mx_yvalues(nd,i)

         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

       end if

     end do

   end do

 ;

 ; Here's how to do the legend by hand.

 ;

 ; mkres               = True    ; Marker resources

 ; txres               = True    ; Text resources

 ; mkres@gsMarkerIndex = 16

 ; mkres@gsMarkerSizeF = 0.02

 ; txres@txFontHeightF = 0.02

 ; txres@txJust        = "CenterLeft"

 ;

 ; Change these values if you want to move the marker legend location.

 ; These values are in the same data space as the plot. 

 ;

 ; xlg1_cen = 0.2

 ; ylg1_cen = 900.

 ; xlg2_cen = 0.2

 ; ylg2_cen = 760.

 ; mkres@gsMarkerColor = "brown"

 ; lnres@gsLineColor   = "brown"

 ; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres)

 ; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \

 ;                                           (/ylg1_cen-60,ylg1_cen+60/),lnres)

 ; lg_cap_legend11  = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \

 ;                                           (/ylg1_cen-60,ylg1_cen-60/),lnres)

 ; lg_cap_legend12  = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \

 ;                                           (/ylg1_cen+60,ylg1_cen+60/),lnres)

 ; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres)

 ; mkres@gsMarkerColor = "blue"

 ; lnres@gsLineColor   = "blue"

 ; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres)

 ; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \

 ;                                           (/ylg2_cen-60,ylg2_cen+60/),lnres)

 ; lg_cap_legend21  = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \

 ;                                           (/ylg2_cen-60,ylg2_cen-60/),lnres)

 ; lg_cap_legend22  = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \

 ;                                           (/ylg2_cen+60,ylg2_cen+60/),lnres)

 ; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres)

   draw(xy)

   frame(wks)

   system("convert xy.ps xy.png")

  u = yvalues(0,:)

  v = yvalues(1,:)

  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)

 ;old eq

 ;bias = sum(((vv-uu)/uu)^2)

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

 ;new eq

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

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

  print (bias)

  print (M)

 end