clamp: turnover/03.biome_ob.ncl@4be95183fbcd

     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"

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

     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 ; casa only

    68 ; all the plant pools (leaf, wood, and fine root) and

    69 ; coarse woody debris (cwd) and litter pools for

    70 ; CASA need to be divided by 1200.  The soil flux

    71 ; and turnover time are fine and do not need to be adjusted.

    73   if (k .ne. n_comp-1) then

    74       flux = flux/1200.

    75   end if

    76 ;************************************************

    77 ; read data: observed-biome

    78 ;************************************************

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

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

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

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

    87   classob = tofloat(fo->LAND_CLASS)

    89   nclass = 20

    91 ;*******************************************************************

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

    93 ;********************************************************************

    95   nclassn     = nclass + 1

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

    97 ; print (range)

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

   100 ; particular range, and then take an average.

   102   nr           = dimsizes(range)

   103   nx           = nr-1

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

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

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

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

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

   110 ; get data

   112   base_1D  = ndtooned(classob)

   113   data1_1D = ndtooned(pool)

   114   data2_1D = ndtooned(flux)

   116 ; output

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

   119   count   = new((/2,nx/),float)

   121   do nd=0,1

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

   125     base = base_1D

   127     if(nd.eq.0) then

   128       data = data1_1D

   129     else

   130       data = data2_1D

   131     end if

   133 ; Loop through each range, using base.

   135     do i=0,nr-2

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

   137 ;        print("")

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

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

   140       else

   141 ;        print("")

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

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

   144       end if

   146 ;     Calculate average

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

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

   150         count(nd,i)   = dimsizes(idx)

   151       else

   152         yvalues(nd,i) = yvalues@_FillValue

   153         count(nd,i)   = 0

   154       end if

   156 ;#############################################################

   157 ; set the following 4 classes to _FillValue:

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

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

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

   162          yvalues(nd,i) = yvalues@_FillValue

   163          count(nd,i)   = 0

   164       end if

   165 ;#############################################################

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

   169 ; Clean up for next time in loop.

   171       delete(idx)

   172     end do

   174     delete(data)

   175   end do

   177 ;============================

   178 ;compute turnover time

   179 ;============================

   181  u       = yvalues(0,:)

   182  v       = yvalues(1,:)

   183  u_count = count(0,:)

   184  v_count = count(1,:)

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

   188  uu       = u(good)

   189  vv       = v(good)

   190  uu_count = u_count(good)

   191  vv_count = v_count(good)

   193  n_biome = dimsizes(uu)

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

   196  t_biome = uu/vv

   198 ;t_biome_avg = avg(t_biome)

   199  t_biome_avg = sum(uu*uu_count)/sum(vv*vv_count)

   201 ;print (t_biome)

   202 ;print (t_biome_avg)

   204 ;===========================

   205 ; for html table - biome

   206 ;===========================

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

   210 ; column (not including header column)

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

   214   ncol = dimsizes(col_head)

   216 ; row (not including header row)

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

   218 ;                    Unclassified, Permanent Snow and Ice

   220   row_head  = (/"Evergreen Needleleaf Forests" \

   221                ,"Evergreen Broadleaf Forests" \

   222                ,"Deciduous Needleleaf Forest" \

   223                ,"Deciduous Broadleaf Forests" \

   224                ,"Mixed Forests" \

   225                ,"Closed Bushlands" \

   226                ,"Open Bushlands" \

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

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

   229                ,"Grasslands" \

   230                ,"Permanent Wetlands" \

   231                ,"Croplands" \

   232                ,"Cropland/Natural Vegetation Mosaic" \

   233                ,"Barren or Sparsely Vegetated" \

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

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

   236                ,"All Biome" \

   237                /)

   238   nrow = dimsizes(row_head)

   240 ; arrays to be passed to table.

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

   243  do i=0,nrow-2

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

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

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

   247  end do

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

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

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

   252 ;**************************************************

   253 ; html table

   254 ;**************************************************

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

   258   header = (/"<HTML>" \

   259             ,"<HEAD>" \

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

   261             ,"</HEAD>" \

   262             ,header_text \

   263             /)

   264   footer = "</HTML>"

   266   table_header = (/ \

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

   268        ,"<tr>" \

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

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

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

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

   273        ,"</tr>" \

   274        /)

   275   table_footer = "</table>"

   276   row_header = "<tr>"

   277   row_footer = "</tr>"

   279   lines = new(50000,string)

   280   nline = 0

   282   set_line(lines,nline,header)

   283   set_line(lines,nline,table_header)

   284 ;-----------------------------------------------

   285 ;row of table

   287   do n = 0,nrow-1

   288      set_line(lines,nline,row_header)

   290      txt1  = row_head(n)

   291      txt2  = text4(n,0)

   292      txt3  = text4(n,1)

   293      txt4  = text4(n,2)

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

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

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

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

   300      set_line(lines,nline,row_footer)

   301   end do

   302 ;-----------------------------------------------

   303   set_line(lines,nline,table_footer)

   304   set_line(lines,nline,footer)

   306 ; Now write to an HTML file.

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

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

   309     asciiwrite(output_html,lines(idx))

   310   else

   311    print ("error?")

   312   end if

   314   delete (idx)

   315 end do

   316 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--