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