clamp: turnover/02.biome.ncl@0c6405ab2ff4

     1 ;********************************************************

     2 ; required command line input parameters:

     3 ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl

4 ;

     5 ; histogram normalized by rain and compute correleration

     6 ;**************************************************************

     7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"

     8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"

     9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

    10 ;**************************************************************

    11 procedure set_line(lines:string,nline:integer,newlines:string)

    12 begin

    13 ; add line to ascci/html file

    15   nnewlines = dimsizes(newlines)

    16   if(nline+nnewlines-1.ge.dimsizes(lines))

    17     print("set_line: bad index, not setting anything.")

    18     return

    19   end if

    20   lines(nline:nline+nnewlines-1) = newlines

    21 ;  print ("lines = " + lines(nline:nline+nnewlines-1))

    22   nline = nline + nnewlines

    23   return

    24 end

    25 ;**************************************************************

    26 ; Main code.

    27 begin

    29  plot_type     = "ps"

    30  plot_type_new = "png"

    32 ;components

    34  component = (/"Leaf","Wood","Fine_Root","Litter","Coarse_Woody_Debris","Soil"/)

    35  n_comp = dimsizes(component)

    37  field_pool = (/"LEAFC","WOODC","FROOTC","LITTERC","CWDC","SOILC"/)

    38  field_flux = (/"LEAFC_ALLOC","WOODC_ALLOC","FROOTC_ALLOC","LITTERC_LOSS","CWDC_LOSS","SOILC_HR"/)

    40 ;************************************************

    41 ; read data: model

    42 ;************************************************

    44  model_grid = "T42"

    46 ;model_name = "i01.06cn"

    47  model_name = "i01.06casa"

    49  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    50  film = model_name + "_1980-2004_ANN_climo.nc"

    52  fm   = addfile (dirm+film,"r")

    54  do k = 0,n_comp-1

    56  pool  = fm->$field_pool(k)$

    57  flux  = fm->$field_flux(k)$

    59 ;Units for these variables are:

    60 ;pool: g C/m^2

    61 ;flux: g C/m^2/s

    63  nsec_per_year = 60*60*24*365

    65  flux = flux *  nsec_per_year

    67 ;************************************************

    68 ; read data: observed-biome

    69 ;************************************************

    71   ob_name = "MODIS MOD 15A2 2000-2005"

    73   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"

    74   filo  = "land_class_"+model_grid+".nc"

    76   fo = addfile(diro+filo,"r")

    78   classob = tofloat(fo->LAND_CLASS)

    80   nclass = 20

    82 ;*******************************************************************

    83 ; Calculate "nice" bins for binning the data in equally spaced ranges

    84 ;********************************************************************

    86   nclassn     = nclass + 1

    87   range       = fspan(0,nclassn-1,nclassn)

    88 ; print (range)

    90 ; Use this range information to grab all the values in a

    91 ; particular range, and then take an average.

    93   nr           = dimsizes(range)

    94   nx           = nr-1

    95   xvalues      = new((/2,nx/),float)

    96   xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.

    97   dx           = xvalues(0,1) - xvalues(0,0)       ; range width

    98   dx4          = dx/4                              ; 1/4 of the range

    99   xvalues(1,:) = xvalues(0,:) - dx/5.

   101 ; get data

   103   base_1D  = ndtooned(classob)

   104   data1_1D = ndtooned(pool)

   105   data2_1D = ndtooned(flux)

   107 ; output

   109   yvalues      = new((/2,nx/),float)

   111   do nd=0,1

   113 ;   See if we are doing data1 (nd=0) or data2 (nd=1).

   115     base = base_1D

   117     if(nd.eq.0) then

   118       data = data1_1D

   119     else

   120       data = data2_1D

   121     end if

   123 ; Loop through each range, using base.

   125     do i=0,nr-2

   126       if (i.ne.(nr-2)) then

   127 ;        print("")

   128 ;        print("In range ["+range(i)+","+range(i+1)+")")

   129          idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))

   130       else

   131 ;        print("")

   132 ;        print("In range ["+range(i)+",)")

   133          idx = ind(base.ge.range(i))

   134       end if

   136 ;     Calculate average

   138       if(.not.any(ismissing(idx))) then

   139         yvalues(nd,i)    = avg(data(idx))

   140         count = dimsizes(idx)

   141       else

   142         yvalues(nd,i)    = yvalues@_FillValue

   143         count            = 0

   144       end if

   146 ;#############################################################

   147 ; set the following 4 classes to _FillValue:

   148 ; Water Bodies(0), Urban and Build-Up(13),

   149 ; Permenant Snow and Ice(15), Unclassified(17)

   151       if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then

   152          yvalues(nd,i)    = yvalues@_FillValue

   153          count            = 0

   154       end if

   155 ;#############################################################

   157 ;     print(nd + ": " + count + " points, avg = " + yvalues(nd,i))

   159 ; Clean up for next time in loop.

   161       delete(idx)

   162     end do

   164     delete(data)

   165   end do

   167 ;============================

   168 ;compute turnover time

   169 ;============================

   171  u = yvalues(0,:)

   172  v = yvalues(1,:)

   174  good = ind(.not.ismissing(u) .and. .not.ismissing(v))

   176  uu = u(good)

   177  vv = v(good)

   179  n_biome = dimsizes(uu)

   181  t_biome = new((/n_biome/),float)

   183  t_biome = uu/vv

   185  t_biome_avg = avg(t_biome)

   187 ;print (t_biome)

   188 ;print (t_biome_avg)

   190 ;===========================

   191 ; for html table - biome

   192 ;===========================

   194   output_html = "table_"+component(k)+".html"

   196 ; column (not including header column)

   198   col_head = (/component(k)+" Flux",component(k)+" Pool",component(k)+" Turnover Time"/)

   200   ncol = dimsizes(col_head)

   202 ; row (not including header row)

   203 ; 4 classes removed: Water Bodies, Urban and Build-Up,

   204 ;                    Unclassified, Permanent Snow and Ice

   206   row_head  = (/"Evergreen Needleleaf Forests" \

   207                ,"Evergreen Broadleaf Forests" \

   208                ,"Deciduous Needleleaf Forest" \

   209                ,"Deciduous Broadleaf Forests" \

   210                ,"Mixed Forests" \

   211                ,"Closed Bushlands" \

   212                ,"Open Bushlands" \

   213                ,"Woody Savannas (S. Hem.)" \

   214                ,"Savannas (S. Hem.)" \

   215                ,"Grasslands" \

   216                ,"Permanent Wetlands" \

   217                ,"Croplands" \

   218                ,"Cropland/Natural Vegetation Mosaic" \

   219                ,"Barren or Sparsely Vegetated" \

   220                ,"Woody Savannas (N. Hem.)" \

   221                ,"Savannas (N. Hem.)" \

   222                ,"All Biome" \

   223                /)

   224   nrow = dimsizes(row_head)

   226 ; arrays to be passed to table.

   227   text4 = new ((/nrow, ncol/),string )

   229  do i=0,nrow-2

   230   text4(i,0) = sprintf("%.1f",vv(i))

   231   text4(i,1) = sprintf("%.1f",uu(i))

   232   text4(i,2) = sprintf("%.1f",t_biome(i))

   233  end do

   234   text4(nrow-1,0) = "-"

   235   text4(nrow-1,1) = "-"

   236   text4(nrow-1,2) = sprintf("%.1f",t_biome_avg)

   238 ;**************************************************

   239 ; html table

   240 ;**************************************************

   242   header_text = "<H1>"+component(k)+" Turnover Time:  Model "+model_name+"</H1>"

   244   header = (/"<HTML>" \

   245             ,"<HEAD>" \

   246             ,"<TITLE>CLAMP metrics</TITLE>" \

   247             ,"</HEAD>" \

   248             ,header_text \

   249             /)

   250   footer = "</HTML>"

   252   table_header = (/ \

   253         "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \

   254        ,"<tr>" \

   255        ,"   <th bgcolor=DDDDDD >Biome Class</th>" \

   256        ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \

   257        ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \

   258        ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \

   259        ,"</tr>" \

   260        /)

   261   table_footer = "</table>"

   262   row_header = "<tr>"

   263   row_footer = "</tr>"

   265   lines = new(50000,string)

   266   nline = 0

   268   set_line(lines,nline,header)

   269   set_line(lines,nline,table_header)

   270 ;-----------------------------------------------

   271 ;row of table

   273   do n = 0,nrow-1

   274      set_line(lines,nline,row_header)

   276      txt1  = row_head(n)

   277      txt2  = text4(n,0)

   278      txt3  = text4(n,1)

   279      txt4  = text4(n,2)

   281      set_line(lines,nline,"<th>"+txt1+"</th>")

   282      set_line(lines,nline,"<th>"+txt2+"</th>")

   283      set_line(lines,nline,"<th>"+txt3+"</th>")

   284      set_line(lines,nline,"<th>"+txt4+"</th>")

   286      set_line(lines,nline,row_footer)

   287   end do

   288 ;-----------------------------------------------

   289   set_line(lines,nline,table_footer)

   290   set_line(lines,nline,footer)

   292 ; Now write to an HTML file.

   293   idx = ind(.not.ismissing(lines))

   294   if(.not.any(ismissing(idx))) then

   295     asciiwrite(output_html,lines(idx))

   296   else

   297    print ("error?")

   298   end if

   300   delete (idx)

   301 end do

   302 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 26 Jan 2009 22:08:20 -0500
changeset 0	0c6405ab2ff4
permissions	-rw-r--r--