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

 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.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

   data_types    = (/ "Obs",        "Model"       /)

   data_names    = (/ "data.81.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)

 ; data_names    = (/ "data.81.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)

   filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)

   ndata_types   = dimsizes(data_types)

   data_file_obs = addfile(data_names(0),"r")    ; Open obs file

   data_file_mod = addfile(data_names(1),"r")    ; Open model file

 ;

 ; Read four variables from files.

 ;

   PREC_ANN = tofloat(data_file_obs->PREC_ANN)

   TNPP_C   = data_file_obs->TNPP_C

   RAIN     = data_file_mod->RAIN

   NPP    = data_file_mod->NPP

 ;

 ; Units for these four variables are:

 ;

 ; PREC_ANN : mm/year

 ; TNPP_C   : g C/m^2/year

 ; RAIN     : mm/s

 ; NPP    : g C/m^2/s

 ;

 ; We want to convert these to "m/year" and "g C/m^2/year".

 ;

   nsec_per_year = 60*60*24*365                  ; # seconds per year

 ; Do the necessary conversions.

   PREC_ANN = PREC_ANN / 1000.

   RAIN     = (RAIN / 1000.) * nsec_per_year

   NPP    = NPP * nsec_per_year

 ; Redo the units.

   PREC_ANN@units = "m/yr"

   RAIN@units     = "m/yr"

   NPP@units    = "gC/m^2/yr"

   TNPP_C@units   = "gC/m^2/yr"

   pminmax(PREC_ANN,"PREC_ANN")

   pminmax(TNPP_C,"TNPP_C")

   pminmax(RAIN,"RAIN")

   pminmax(NPP,"NPP")

   RAIN_1D     = ndtooned(RAIN)

   NPP_1D    = ndtooned(NPP)

   PREC_ANN_1D = ndtooned(PREC_ANN)

   TNPP_C_1D   = ndtooned(TNPP_C)

 ;

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

 ; ranges.

 ;

   nbins       = 15     ; Number of bins to use.

   nicevals    = nice_mnmxintvl(min(RAIN_1D),max(RAIN_1D),nbins,True)

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

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

 ;

 ; 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(RAIN_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(RAIN_1D))

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

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

   do nd=0,1

 ;

 ; See if we are doing model or observational data.

 ;

     if(nd.eq.0) then

       data     = PREC_ANN_1D

       npp_data = TNPP_C_1D

     else

       data     = RAIN_1D

       npp_data = NPP_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(data_types(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

   xvalues@long_name = "Mean Annual precipitation (m/year)"

   yvalues@long_name = "NPP (g C/m2/year)"                 

 ;

 ; Start the graphics.

 ;

 ; wks = gsn_open_wks("x11","npp")

   wks = gsn_open_wks("png","npp")

   res             = True

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

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

   res@gsnMaximize = False

   res@gsnDraw     = False

   res@gsnFrame    = False

   res@xyMarkLineMode = "Markers"

   res@xyMarkerSizeF   = 0.014

   res@xyMarker       = 16

 ; res@xyMarkerColors = (/"Gray25","Gray50"/)

   res@xyMarkerColors = (/"brown","blue"/)

   res@trYMinF        = min(mn_yvalues) - 10.

   res@trYMaxF        = max(mx_yvalues) + 10.

   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

   draw(xy)

   frame(wks)

 end