diff -r 000000000000 -r 0c6405ab2ff4 turnover/01.biome.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/turnover/01.biome.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,294 @@
+;********************************************************
+; required command line input parameters:
+;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
+;
+; 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 set_line(lines:string,nline:integer,newlines:string) 
+begin
+; add line to ascci/html file
+    
+  nnewlines = dimsizes(newlines)
+  if(nline+nnewlines-1.ge.dimsizes(lines))
+    print("set_line: bad index, not setting anything.") 
+    return
+  end if 
+  lines(nline:nline+nnewlines-1) = newlines
+;  print ("lines = " + lines(nline:nline+nnewlines-1))
+  nline = nline + nnewlines
+  return 
+end
+;**************************************************************
+; Main code.
+begin
+ 
+ plot_type     = "ps"
+ plot_type_new = "png"
+
+;************************************************
+; read data: model       
+;************************************************
+
+ model_grid = "T42"
+
+ model_name = "i01.06cn"
+;model_name = "i01.06casa"
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film = model_name + "_1980-2004_ANN_climo.nc"
+
+ fm   = addfile (dirm+film,"r")
+  
+ xm     = fm->lon  
+ ym     = fm->lat
+
+ pool  = fm->WOODC
+ flux  = fm->WOODC_ALLOC
+
+;Units for these variables are:
+;pool: g C/m^2
+;flux: g C/m^2/s
+
+ nsec_per_year = 60*60*24*365
+  
+ flux = flux *  nsec_per_year 
+      
+;************************************************
+; read data: observed-biome
+;************************************************
+
+  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)
+
+  nclass = 20
+                
+;*******************************************************************
+; 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
+
+  base_1D  = ndtooned(classob)
+  data1_1D = ndtooned(pool)
+  data2_1D = ndtooned(flux)
+
+; output
+
+  yvalues      = new((/2,nx/),float)
+
+  do nd=0,1
+
+;   See if we are doing data1 (nd=0) or data2 (nd=1).
+
+    base = base_1D
+
+    if(nd.eq.0) then
+      data = data1_1D
+    else
+      data = data2_1D
+    end if
+
+; Loop through each range, using base.
+
+    do i=0,nr-2
+      if (i.ne.(nr-2)) then
+;        print("")
+;        print("In range ["+range(i)+","+range(i+1)+")")
+         idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
+      else
+;        print("")
+;        print("In range ["+range(i)+",)")
+         idx = ind(base.ge.range(i))
+      end if
+
+;     Calculate average 
+
+      if(.not.any(ismissing(idx))) then
+        yvalues(nd,i)    = avg(data(idx))
+        count = dimsizes(idx)
+      else
+        yvalues(nd,i)    = yvalues@_FillValue
+        count            = 0
+      end if
+
+;#############################################################
+; set the following 4 classes to _FillValue:
+; Water Bodies(0), Urban and Build-Up(13),
+; Permenant Snow and Ice(15), Unclassified(17)
+
+      if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
+         yvalues(nd,i)    = yvalues@_FillValue
+         count            = 0
+      end if
+;############################################################# 
+
+;     print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
+
+; Clean up for next time in loop.
+
+      delete(idx)
+    end do
+
+    delete(data)
+  end do
+
+;============================
+;compute turnover time
+;============================
+
+ u = yvalues(0,:)
+ v = yvalues(1,:)
+
+ good = ind(.not.ismissing(u) .and. .not.ismissing(v))
+
+ uu = u(good)
+ vv = v(good) 
+
+ n_biome = dimsizes(uu)
+
+ t_biome = new((/n_biome/),float)
+
+ t_biome = uu/vv
+
+ t_biome_avg = avg(t_biome)
+  
+;print (t_biome)
+;print (t_biome_avg)
+
+;===========================
+; for html table - biome
+;===========================
+
+  output_html = "table_biome.html"
+
+; column (not including header column)
+
+  col_head = (/"WOOD Flux","WOOD Pool","WOOD Turnover Time"/)
+
+  ncol = dimsizes(col_head)
+
+; row (not including header row)
+; 4 classes removed: Water Bodies, Urban and Build-Up, 
+;                    Unclassified, Permanent Snow and Ice                   
+
+  row_head  = (/"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" \                                           
+               ,"Cropland/Natural Vegetation Mosaic" \             
+               ,"Barren or Sparsely Vegetated" \                             
+               ,"Woody Savannas (N. Hem.)" \           
+               ,"Savannas (N. Hem.)" \
+               ,"All Biome" \                
+               /)  
+  nrow = dimsizes(row_head)                  
+
+; arrays to be passed to table. 
+  text4 = new ((/nrow, ncol/),string )
+ 
+ do i=0,nrow-2
+  text4(i,0) = sprintf("%.1f",vv(i))
+  text4(i,1) = sprintf("%.1f",uu(i))
+  text4(i,2) = sprintf("%.1f",t_biome(i))
+ end do
+  text4(nrow-1,0) = "-"
+  text4(nrow-1,1) = "-"
+  text4(nrow-1,2) = sprintf("%.1f",t_biome_avg)
+
+;**************************************************
+; html table
+;**************************************************
+
+  header_text = "<H1>WOOD Turnover Time:  Model "+model_name+"</H1>" 
+
+  header = (/"<HTML>" \
+            ,"<HEAD>" \
+            ,"<TITLE>CLAMP metrics</TITLE>" \
+            ,"</HEAD>" \
+            ,header_text \
+            /) 
+  footer = "</HTML>"
+
+  table_header = (/ \
+        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
+       ,"<tr>" \
+       ,"   <th bgcolor=DDDDDD >Biome Class</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
+       ,"</tr>" \
+       /)
+  table_footer = "</table>"
+  row_header = "<tr>"
+  row_footer = "</tr>"
+
+  lines = new(50000,string)
+  nline = 0
+
+  set_line(lines,nline,header)
+  set_line(lines,nline,table_header)
+;-----------------------------------------------
+;row of table
+
+  do n = 0,nrow-1
+     set_line(lines,nline,row_header)
+
+     txt1  = row_head(n)
+     txt2  = text4(n,0)
+     txt3  = text4(n,1)
+     txt4  = text4(n,2)
+
+     set_line(lines,nline,"<th>"+txt1+"</th>")
+     set_line(lines,nline,"<th>"+txt2+"</th>")
+     set_line(lines,nline,"<th>"+txt3+"</th>")
+     set_line(lines,nline,"<th>"+txt4+"</th>")
+
+     set_line(lines,nline,row_footer)
+  end do
+;-----------------------------------------------
+  set_line(lines,nline,table_footer)
+  set_line(lines,nline,footer) 
+
+; Now write to an HTML file.
+  idx = ind(.not.ismissing(lines))
+  if(.not.any(ismissing(idx))) then
+    asciiwrite(output_html,lines(idx))
+  else
+   print ("error?")
+  end if
+
+end
+