diff -r 000000000000 -r 0c6405ab2ff4 lai/41.table_mean.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lai/41.table_mean.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,166 @@
+;********************************************************
+; 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
+