diff -r 000000000000 -r 0c6405ab2ff4 npp/22.histogram.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/npp/22.histogram.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,238 @@
+ 
+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"/) /)
+  filevar_names = (/ (/"PREC_ANN","ANPP_C"/), (/"RAIN","AGNPP"/) /)
+  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 in data: observed
+;************************************************
+ PREC_ANN = tofloat(data_file_obs->PREC_ANN)
+ TNPP_C   = data_file_obs->TNPP_C
+ xo    = data_file_obs->LONG_DD  
+ yo    = data_file_obs->LAT_DD   
+
+; change longitude from (-180,180) to (0,360)
+;print (xo)
+ nx = dimsizes(xo)
+ do i= 0,nx-1
+    if (xo(i) .lt. 0.) then
+        xo(i) = xo(i)+ 360.
+    end if
+ end do
+;print (xo)
+ 
+;************************************************
+; read in data: model       
+;************************************************
+ ai    = data_file_mod->RAIN
+ bi    = data_file_mod->NPP
+ xi    = data_file_mod->lon      
+ yi    = data_file_mod->lat      
+
+;************************************************
+; interpolate from model grid to observed grid
+;************************************************
+ RAIN = linint2_points(xi,yi,ai,True,xo,yo,0) 
+ NPP  = linint2_points(xi,yi,bi,True,xo,yo,0) 
+ 
+;************************************************
+; convert unit
+;************************************************
+; 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("png","npp")
+
+  res             = True 
+  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.
+; res@tiMainString   = "Observed vs i01.03cn_81site"
+  res@tiMainString   = "Observed vs i01.04casa_81site"
+
+  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
+