diff -r 000000000000 -r 0c6405ab2ff4 beta/04.biome.ncl
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/beta/04.biome.ncl Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,290 @@
+;********************************************************
+; 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
+
+ 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"
+
+ model_name = model_name_f
+
+ 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)
+
+;*******************************************************************
+; 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
+;============================
+
+ nsec_per_year = 60*60*24*365
+
+ u = yvalues(0,:)
+ v = yvalues(1,:)
+
+ good = ind(.not.ismissing(u) .and. .not.ismissing(v))
+ uu = u(good)* nsec_per_year
+ vv = v(good)* nsec_per_year
+
+ 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 (beta_biome_avg)
+;*******************************************************************
+; for html table
+;*******************************************************************
+
+; column (not including header column)
+
+ col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta"/)
+
+ ncol = dimsizes(col_head)
+
+; row (not including header row)
+ row_head = (/"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" \
+ ,"Cropland/Natural Vegetation Mosaic" \
+ ,"Permanent Snow and Ice" \
+ ,"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("%.2f",co2_i)
+ text4(i,1) = sprintf("%.2f",co2_f)
+ text4(i,2) = sprintf("%.2f",uu(i))
+ text4(i,3) = sprintf("%.2f",vv(i))
+ text4(i,4) = sprintf("%.2f",beta_biome(i))
+ end do
+ text4(nrow-1,0) = "-"
+ text4(nrow-1,1) = "-"
+ text4(nrow-1,2) = "-"
+ text4(nrow-1,3) = "-"
+ text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
+
+;**************************************************
+; html table
+;**************************************************
+ output_html = "table_biome.html"
+
+ header_text = "Beta Factor: Model "+model_name+"
"
+
+ header = (/"" \
+ ,"" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,header_text \
+ /)
+ footer = ""
+
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," | Biome Class | " \
+ ," CO2_i | " \
+ ," CO2_f | " \
+ ," NPP_i | " \
+ ," NPP_f | " \
+ ," Beta | " \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ 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)
+ txt5 = text4(n,3)
+ txt6 = text4(n,4)
+
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+ set_line(lines,nline,""+txt6+" | ")
+
+ 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
+