clamp: turnover/01.biome.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.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 ;************************************************

    33 ; read data: model

    34 ;************************************************

    36  model_grid = "T42"

    38  model_name = "i01.06cn"

    39 ;model_name = "i01.06casa"

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

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

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

    46  xm     = fm->lon

    47  ym     = fm->lat

    49  pool  = fm->WOODC

    50  flux  = fm->WOODC_ALLOC

    52 ;Units for these variables are:

    53 ;pool: g C/m^2

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

    56  nsec_per_year = 60*60*24*365

    58  flux = flux *  nsec_per_year

    60 ;************************************************

    61 ; read data: observed-biome

    62 ;************************************************

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

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

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

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

    71   classob = tofloat(fo->LAND_CLASS)

    73   nclass = 20

    75 ;*******************************************************************

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

    77 ;********************************************************************

    79   nclassn     = nclass + 1

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

    81 ; print (range)

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

    84 ; particular range, and then take an average.

    86   nr           = dimsizes(range)

    87   nx           = nr-1

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

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

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

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

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

    94 ; get data

    96   base_1D  = ndtooned(classob)

    97   data1_1D = ndtooned(pool)

    98   data2_1D = ndtooned(flux)

   100 ; output

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

   104   do nd=0,1

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

   108     base = base_1D

   110     if(nd.eq.0) then

   111       data = data1_1D

   112     else

   113       data = data2_1D

   114     end if

   116 ; Loop through each range, using base.

   118     do i=0,nr-2

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

   120 ;        print("")

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

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

   123       else

   124 ;        print("")

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

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

   127       end if

   129 ;     Calculate average

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

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

   133         count = dimsizes(idx)

   134       else

   135         yvalues(nd,i)    = yvalues@_FillValue

   136         count            = 0

   137       end if

   139 ;#############################################################

   140 ; set the following 4 classes to _FillValue:

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

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

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

   145          yvalues(nd,i)    = yvalues@_FillValue

   146          count            = 0

   147       end if

   148 ;#############################################################

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

   152 ; Clean up for next time in loop.

   154       delete(idx)

   155     end do

   157     delete(data)

   158   end do

   160 ;============================

   161 ;compute turnover time

   162 ;============================

   164  u = yvalues(0,:)

   165  v = yvalues(1,:)

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

   169  uu = u(good)

   170  vv = v(good)

   172  n_biome = dimsizes(uu)

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

   176  t_biome = uu/vv

   178  t_biome_avg = avg(t_biome)

   180 ;print (t_biome)

   181 ;print (t_biome_avg)

   183 ;===========================

   184 ; for html table - biome

   185 ;===========================

   187   output_html = "table_biome.html"

   189 ; column (not including header column)

   191   col_head = (/"WOOD Flux","WOOD Pool","WOOD Turnover Time"/)

   193   ncol = dimsizes(col_head)

   195 ; row (not including header row)

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

   197 ;                    Unclassified, Permanent Snow and Ice

   199   row_head  = (/"Evergreen Needleleaf Forests" \

   200                ,"Evergreen Broadleaf Forests" \

   201                ,"Deciduous Needleleaf Forest" \

   202                ,"Deciduous Broadleaf Forests" \

   203                ,"Mixed Forests" \

   204                ,"Closed Bushlands" \

   205                ,"Open Bushlands" \

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

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

   208                ,"Grasslands" \

   209                ,"Permanent Wetlands" \

   210                ,"Croplands" \

   211                ,"Cropland/Natural Vegetation Mosaic" \

   212                ,"Barren or Sparsely Vegetated" \

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

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

   215                ,"All Biome" \

   216                /)

   217   nrow = dimsizes(row_head)

   219 ; arrays to be passed to table.

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

   222  do i=0,nrow-2

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

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

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

   226  end do

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

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

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

   231 ;**************************************************

   232 ; html table

   233 ;**************************************************

   235   header_text = "<H1>WOOD Turnover Time:  Model "+model_name+"</H1>"

   237   header = (/"<HTML>" \

   238             ,"<HEAD>" \

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

   240             ,"</HEAD>" \

   241             ,header_text \

   242             /)

   243   footer = "</HTML>"

   245   table_header = (/ \

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

   247        ,"<tr>" \

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

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

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

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

   252        ,"</tr>" \

   253        /)

   254   table_footer = "</table>"

   255   row_header = "<tr>"

   256   row_footer = "</tr>"

   258   lines = new(50000,string)

   259   nline = 0

   261   set_line(lines,nline,header)

   262   set_line(lines,nline,table_header)

   263 ;-----------------------------------------------

   264 ;row of table

   266   do n = 0,nrow-1

   267      set_line(lines,nline,row_header)

   269      txt1  = row_head(n)

   270      txt2  = text4(n,0)

   271      txt3  = text4(n,1)

   272      txt4  = text4(n,2)

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

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

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

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

   279      set_line(lines,nline,row_footer)

   280   end do

   281 ;-----------------------------------------------

   282   set_line(lines,nline,table_footer)

   283   set_line(lines,nline,footer)

   285 ; Now write to an HTML file.

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

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

   288     asciiwrite(output_html,lines(idx))

   289   else

   290    print ("error?")

   291   end if

   293 end

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