diff -r 000000000000 -r 0c6405ab2ff4 carbon_sink/11.table1.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/carbon_sink/11.table1.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,447 @@
+;********************************************************
+;using model biome vlass
+;
+; 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"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+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"
+
+;************************************************
+; model name and grid       
+;************************************************
+
+ model_grid = "T42"
+
+ model_name  = "cn"
+ model_name1 = "i01.06cn"
+ model_name2 = "i01.10cn"
+
+;---------------------------------------------------
+; get biome data: model
+
+  biome_name_mod = "Model PFT Class"
+
+  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+  film = "class_pft_"+model_grid+".nc"
+
+  fm = addfile(dirm+film,"r")
+ 
+  classmod = fm->CLASS_PFT
+
+  delete (fm)
+
+; model data has 17 land-type classes
+
+  nclass_mod = 17
+
+;--------------------------------------------------
+; get model data: landfrac and area
+
+ dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/" 
+ film_l = "lnd_T42.nc"
+ fm_l   = addfile (dirm_l+film_l,"r")
+  
+ landfrac = fm_l->landfrac
+ area     = fm_l->area
+
+; change area from km**2 to m**2
+  area = area * 1.e6             
+;---------------------------------------------------
+; read data: model, group 1
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film = model_name1 + "_1980-2004_ANN_climo.nc"
+ fm   = addfile (dirm+film,"r")
+ 
+ NPP1   = fm->NPP
+
+ leafc  = fm->LEAFC
+ woodc  = fm->WOODC
+ frootc = fm->FROOTC
+ VegC   = leafc
+ VegC   = leafc + woodc + frootc
+
+ litterc = fm->LITTERC
+ cwdc    = fm->CWDC
+ LiCwC   = litterc
+ LiCwC   = litterc + cwdc
+
+ SoilC   = fm->SOILC
+
+ delete (fm)
+;--------------------------------------------------- 
+; read data: model, group 2
+
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ film = model_name2 + "_1990-2004_ANN_climo.nc"
+ fm   = addfile (dirm+film,"r")
+
+ NPP2   = fm->NPP 
+ NEE2   = fm->NEE 
+
+;---------------------------------------------------
+; Units for these variables are:
+
+;NPP1: g C/m^2/s
+;NPP2: g C/m^2/s
+;NEE2: g C/m^2/s
+
+;VegC:  g C/m^2
+;LiCwC: g C/m^2
+;SoilC: g C/m^2
+
+ nsec_per_year = 60*60*24*365
+
+; change unit to g C/m^2/year
+  
+ NPP1 = NPP1 *  nsec_per_year
+ NPP2 = NPP2 *  nsec_per_year
+ NEE2 = NEE2 *  nsec_per_year 
+
+;---------------------------------------------------
+; take into account landfrac
+
+ area(:,:)    = area(:,:)    * landfrac(:,:)
+ NPP1(0,:,:)  = NPP1(0,:,:)  * landfrac(:,:)
+ VegC(0,:,:)  = VegC(0,:,:)  * landfrac(:,:)
+ LiCwC(0,:,:) = LiCwC(0,:,:) * landfrac(:,:)
+ SoilC(0,:,:) = SoilC(0,:,:) * landfrac(:,:)
+ NPP2(0,:,:)  = NPP2(0,:,:)  * landfrac(:,:)
+ NEE2(0,:,:)  = NEE2(0,:,:)  * landfrac(:,:)
+
+ data_n = 7
+                
+;*******************************************************************
+; Calculate "nice" bins for binning the data in equally spaced ranges
+;********************************************************************
+
+; using model biome class
+  nclass      = nclass_mod
+
+  range       = fspan(0,nclass,nclass+1)
+
+; print (range)
+; Use this range information to grab all the values in a
+; particular range, and then take an average.
+
+  nx           = dimsizes(range) - 1
+
+;==============================
+; put data into bins
+;==============================
+
+; using observed biome class
+; base  = ndtooned(classob)
+; using model biome class
+  base  = ndtooned(classmod)
+
+; output
+
+  yvalues = new((/data_n,nx/),float)
+  count   = new((/data_n,nx/),float)
+
+  do n = 0,data_n-1
+
+    if(n.eq.0) then
+      data = ndtooned(area)
+    end if
+
+    if(n.eq.1) then
+      data = ndtooned(NPP1)
+    end if
+
+    if(n.eq.2) then
+      data = ndtooned(VegC)
+    end if
+
+    if(n.eq.3) then
+      data = ndtooned(LiCwC)
+    end if
+
+    if(n.eq.4) then
+      data = ndtooned(SoilC)
+    end if
+
+    if(n.eq.5) then
+      data = ndtooned(NPP2)
+    end if
+
+    if(n.eq.6) then
+      data = ndtooned(NEE2)
+    end if
+
+; Loop through each range, using base.
+
+    do i=0,nx-1
+      if (i.ne.(nx-1)) 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
+        if (n.eq.0) then 
+           yvalues(n,i) = sum(data(idx))
+        else 
+           yvalues(n,i) = avg(data(idx))
+        end if
+
+        count(n,i)   = dimsizes(idx)
+      else
+        yvalues(n,i) = yvalues@_FillValue
+        count(n,i)   = 0
+      end if
+
+;#############################################################
+; using model biome class:
+;     set the following 4 classes to _FillValue:
+;     (3)Needleleaf Deciduous Boreal Tree,
+;     (8)Broadleaf Deciduous Boreal Tree,
+;     (9)Broadleaf Evergreen Shrub,
+;     (16)Wheat
+
+      if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
+         yvalues(n,i) = yvalues@_FillValue
+         count(n,i)   = 0
+      end if
+;#############################################################  
+
+;     print(n + ": " + count + " points, avg = " + yvalues(n,i))
+
+      delete(idx)
+    end do 
+
+    delete(data)
+  end do
+
+  delete (base)
+  delete (area)
+  delete (NPP1)
+  delete (VegC)
+  delete (LiCwC)
+  delete (SoilC)
+  delete (NPP2)
+  delete (NEE2)
+
+;----------------------------------------------------------------
+; data for table1
+
+ good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))
+;print (good)
+
+ w = yvalues(0,:)
+ area_g = w(good)
+
+ w = yvalues(1,:)
+ NPP1_g = w(good)
+
+ w = yvalues(2,:)
+ VegC_g = w(good)
+
+ w = yvalues(3,:)
+ LiCwC_g = w(good)
+
+ w = yvalues(4,:)
+ SoilC_g = w(good)
+
+ w = yvalues(5,:)
+ NPP2_g = w(good)
+
+ w = yvalues(6,:)
+ NEE2_g = w(good) 
+
+ n_biome = dimsizes(NPP1_g)
+
+ NPP1_t    = new((/n_biome/),float)
+ VegC_t    = new((/n_biome/),float)
+ LiCwC_t   = new((/n_biome/),float)
+ SoilC_t   = new((/n_biome/),float)
+ NEE2_t    = new((/n_biome/),float)
+ NPP_ratio = new((/n_biome/),float)
+
+ NPP_ratio = NPP2_g/NPP1_g
+
+;-----------------------------------------------------------------
+; data for table2
+
+; change unit from g to Pg (Peta gram)
+ factor_unit = 1.e-15
+
+ NPP1_t    = NPP1_g  * area_g * factor_unit
+ VegC_t    = VegC_g  * area_g * factor_unit
+ LiCwC_t   = LiCwC_g * area_g * factor_unit
+ SoilC_t   = SoilC_g * area_g * factor_unit
+ NEE2_t    = NEE2_g  * area_g * factor_unit
+
+ print (NPP1_t)
+ 
+;-------------------------------------------------------------
+; html table1 data
+
+; column (not including header column)
+
+  col_head  = (/"Area (1.e12m2)" \
+               ,"NPP (gC/m2/yr)" \
+               ,"VegC (gC/m2)" \
+               ,"Litter+CWD (gC/m2)" \
+               ,"SoilC (gC/m2)" \
+               ,"NPP ratio" \
+               ,"NEE (gC/m2/yr)" \
+               /)
+
+  ncol = dimsizes(col_head)
+
+; row (not including header row)                   
+
+; using model biome class:  
+  row_head  = (/"Not Vegetated" \
+               ,"Needleleaf Evergreen Temperate Tree" \
+               ,"Needleleaf Evergreen Boreal Tree" \
+;              ,"Needleleaf Deciduous Boreal Tree" \
+               ,"Broadleaf Evergreen Tropical Tree" \
+               ,"Broadleaf Evergreen Temperate Tree" \
+               ,"Broadleaf Deciduous Tropical Tree" \
+               ,"Broadleaf Deciduous Temperate Tree" \
+;              ,"Broadleaf Deciduous Boreal Tree" \
+;              ,"Broadleaf Evergreen Shrub" \
+               ,"Broadleaf Deciduous Temperate Shrub" \
+               ,"Broadleaf Deciduous Boreal Shrub" \
+               ,"C3 Arctic Grass" \
+               ,"C3 Non-Arctic Grass" \
+               ,"C4 Grass" \
+               ,"Corn" \
+;              ,"Wheat" \                      
+               ,"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",area_g(i)*1.e-12)
+  text4(i,1) = sprintf("%.1f",NPP1_g(i))
+  text4(i,2) = sprintf("%.1f",VegC_g(i))
+  text4(i,3) = sprintf("%.1f",LiCwC_g(i))
+  text4(i,4) = sprintf("%.1f",SoilC_g(i))
+  text4(i,5) = sprintf("%.1f",NPP_ratio(i))
+  text4(i,6) = sprintf("%.1f",NEE2_g(i))
+ end do
+  text4(nrow-1,0) = "-"
+  text4(nrow-1,1) = "-"
+  text4(nrow-1,2) = "-"
+
+;-------------------------------------------------------
+; create html table1
+
+  header_text = "<H1>NEE and Carbon Stocks and Fluxes:  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 Type</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
+       ,"   <th bgcolor=DDDDDD >"+col_head(6)+"</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-2
+     set_line(lines,nline,row_header)
+
+     txt0  = row_head(n)
+     txt1  = text4(n,0)
+     txt2  = text4(n,1)
+     txt3  = text4(n,2)
+     txt4  = text4(n,3)
+     txt5  = text4(n,4)
+     txt6  = text4(n,5)
+     txt7  = text4(n,6)
+
+     set_line(lines,nline,"<th>"+txt0+"</th>")
+     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,"<th>"+txt5+"</th>")
+     set_line(lines,nline,"<th>"+txt6+"</th>")
+     set_line(lines,nline,"<th>"+txt7+"</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.
+
+  output_html = "table_carbon_sink1.html"
+
+  idx = ind(.not.ismissing(lines))
+  if(.not.any(ismissing(idx))) then
+    asciiwrite(output_html,lines(idx))
+  else
+   print ("error?")
+  end if
+
+  delete (col_head)
+  delete (row_head)
+  delete (text4)
+
+end
+