clamp: carbon_sink/11.table1.ncl@4be95183fbcd

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

     2 ;using model biome vlass

3 ;

     4 ; required command line input parameters:

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

6 ;

     7 ; histogram normalized by rain and compute correleration

     8 ;**************************************************************

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

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

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

    12 ;**************************************************************

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

    14 begin

    15 ; add line to ascci/html file

    17   nnewlines = dimsizes(newlines)

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

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

    20     return

    21   end if

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

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

    24   nline = nline + nnewlines

    25   return

    26 end

    27 ;**************************************************************

    28 ; Main code.

    29 begin

    31  plot_type     = "ps"

    32  plot_type_new = "png"

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

    35 ; model name and grid

    36 ;************************************************

    38  model_grid = "T42"

    40  model_name  = "cn"

    41  model_name1 = "i01.06cn"

    42  model_name2 = "i01.10cn"

    44 ;---------------------------------------------------

    45 ; get biome data: model

    47   biome_name_mod = "Model PFT Class"

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

    50   film = "class_pft_"+model_grid+".nc"

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

    54   classmod = fm->CLASS_PFT

    56   delete (fm)

    58 ; model data has 17 land-type classes

    60   nclass_mod = 17

    62 ;--------------------------------------------------

    63 ; get model data: landfrac and area

    65  dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/"

    66  film_l = "lnd_T42.nc"

    67  fm_l   = addfile (dirm_l+film_l,"r")

    69  landfrac = fm_l->landfrac

    70  area     = fm_l->area

    72 ; change area from km**2 to m**2

    73   area = area * 1.e6

    74 ;---------------------------------------------------

    75 ; read data: model, group 1

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

    78  film = model_name1 + "_1980-2004_ANN_climo.nc"

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

    81  NPP1   = fm->NPP

    83  leafc  = fm->LEAFC

    84  woodc  = fm->WOODC

    85  frootc = fm->FROOTC

    86  VegC   = leafc

    87  VegC   = leafc + woodc + frootc

    89  litterc = fm->LITTERC

    90  cwdc    = fm->CWDC

    91  LiCwC   = litterc

    92  LiCwC   = litterc + cwdc

    94  SoilC   = fm->SOILC

    96  delete (fm)

    97 ;---------------------------------------------------

    98 ; read data: model, group 2

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

   101  film = model_name2 + "_1990-2004_ANN_climo.nc"

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

   104  NPP2   = fm->NPP

   105  NEE2   = fm->NEE

   107 ;---------------------------------------------------

   108 ; Units for these variables are:

   110 ;NPP1: g C/m^2/s

   111 ;NPP2: g C/m^2/s

   112 ;NEE2: g C/m^2/s

   114 ;VegC:  g C/m^2

   115 ;LiCwC: g C/m^2

   116 ;SoilC: g C/m^2

   118  nsec_per_year = 60*60*24*365

   120 ; change unit to g C/m^2/year

   122  NPP1 = NPP1 *  nsec_per_year

   123  NPP2 = NPP2 *  nsec_per_year

   124  NEE2 = NEE2 *  nsec_per_year

   126 ;---------------------------------------------------

   127 ; take into account landfrac

   129  area(:,:)    = area(:,:)    * landfrac(:,:)

   130  NPP1(0,:,:)  = NPP1(0,:,:)  * landfrac(:,:)

   131  VegC(0,:,:)  = VegC(0,:,:)  * landfrac(:,:)

   132  LiCwC(0,:,:) = LiCwC(0,:,:) * landfrac(:,:)

   133  SoilC(0,:,:) = SoilC(0,:,:) * landfrac(:,:)

   134  NPP2(0,:,:)  = NPP2(0,:,:)  * landfrac(:,:)

   135  NEE2(0,:,:)  = NEE2(0,:,:)  * landfrac(:,:)

   137  data_n = 7

   139 ;*******************************************************************

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

   141 ;********************************************************************

   143 ; using model biome class

   144   nclass      = nclass_mod

   146   range       = fspan(0,nclass,nclass+1)

   148 ; print (range)

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

   150 ; particular range, and then take an average.

   152   nx           = dimsizes(range) - 1

   154 ;==============================

   155 ; put data into bins

   156 ;==============================

   158 ; using observed biome class

   159 ; base  = ndtooned(classob)

   160 ; using model biome class

   161   base  = ndtooned(classmod)

   163 ; output

   165   yvalues = new((/data_n,nx/),float)

   166   count   = new((/data_n,nx/),float)

   168   do n = 0,data_n-1

   170     if(n.eq.0) then

   171       data = ndtooned(area)

   172     end if

   174     if(n.eq.1) then

   175       data = ndtooned(NPP1)

   176     end if

   178     if(n.eq.2) then

   179       data = ndtooned(VegC)

   180     end if

   182     if(n.eq.3) then

   183       data = ndtooned(LiCwC)

   184     end if

   186     if(n.eq.4) then

   187       data = ndtooned(SoilC)

   188     end if

   190     if(n.eq.5) then

   191       data = ndtooned(NPP2)

   192     end if

   194     if(n.eq.6) then

   195       data = ndtooned(NEE2)

   196     end if

   198 ; Loop through each range, using base.

   200     do i=0,nx-1

   201       if (i.ne.(nx-1)) then

   202 ;        print("")

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

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

   205       else

   206 ;        print("")

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

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

   209       end if

   211 ;     Calculate average

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

   214         if (n.eq.0) then

   215            yvalues(n,i) = sum(data(idx))

   216         else

   217            yvalues(n,i) = avg(data(idx))

   218         end if

   220         count(n,i)   = dimsizes(idx)

   221       else

   222         yvalues(n,i) = yvalues@_FillValue

   223         count(n,i)   = 0

   224       end if

   226 ;#############################################################

   227 ; using model biome class:

   228 ;     set the following 4 classes to _FillValue:

   229 ;     (3)Needleleaf Deciduous Boreal Tree,

   230 ;     (8)Broadleaf Deciduous Boreal Tree,

   231 ;     (9)Broadleaf Evergreen Shrub,

   232 ;     (16)Wheat

   234       if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then

   235          yvalues(n,i) = yvalues@_FillValue

   236          count(n,i)   = 0

   237       end if

   238 ;#############################################################

   240 ;     print(n + ": " + count + " points, avg = " + yvalues(n,i))

   242       delete(idx)

   243     end do

   245     delete(data)

   246   end do

   248   delete (base)

   249   delete (area)

   250   delete (NPP1)

   251   delete (VegC)

   252   delete (LiCwC)

   253   delete (SoilC)

   254   delete (NPP2)

   255   delete (NEE2)

   257 ;----------------------------------------------------------------

   258 ; data for table1

   260  good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))

   261 ;print (good)

   263  w = yvalues(0,:)

   264  area_g = w(good)

   266  w = yvalues(1,:)

   267  NPP1_g = w(good)

   269  w = yvalues(2,:)

   270  VegC_g = w(good)

   272  w = yvalues(3,:)

   273  LiCwC_g = w(good)

   275  w = yvalues(4,:)

   276  SoilC_g = w(good)

   278  w = yvalues(5,:)

   279  NPP2_g = w(good)

   281  w = yvalues(6,:)

   282  NEE2_g = w(good)

   284  n_biome = dimsizes(NPP1_g)

   286  NPP1_t    = new((/n_biome/),float)

   287  VegC_t    = new((/n_biome/),float)

   288  LiCwC_t   = new((/n_biome/),float)

   289  SoilC_t   = new((/n_biome/),float)

   290  NEE2_t    = new((/n_biome/),float)

   291  NPP_ratio = new((/n_biome/),float)

   293  NPP_ratio = NPP2_g/NPP1_g

   295 ;-----------------------------------------------------------------

   296 ; data for table2

   298 ; change unit from g to Pg (Peta gram)

   299  factor_unit = 1.e-15

   301  NPP1_t    = NPP1_g  * area_g * factor_unit

   302  VegC_t    = VegC_g  * area_g * factor_unit

   303  LiCwC_t   = LiCwC_g * area_g * factor_unit

   304  SoilC_t   = SoilC_g * area_g * factor_unit

   305  NEE2_t    = NEE2_g  * area_g * factor_unit

   307  print (NPP1_t)

   309 ;-------------------------------------------------------------

   310 ; html table1 data

   312 ; column (not including header column)

   314   col_head  = (/"Area (1.e12m2)" \

   315                ,"NPP (gC/m2/yr)" \

   316                ,"VegC (gC/m2)" \

   317                ,"Litter+CWD (gC/m2)" \

   318                ,"SoilC (gC/m2)" \

   319                ,"NPP ratio" \

   320                ,"NEE (gC/m2/yr)" \

   321                /)

   323   ncol = dimsizes(col_head)

   325 ; row (not including header row)

   327 ; using model biome class:

   328   row_head  = (/"Not Vegetated" \

   329                ,"Needleleaf Evergreen Temperate Tree" \

   330                ,"Needleleaf Evergreen Boreal Tree" \

   331 ;              ,"Needleleaf Deciduous Boreal Tree" \

   332                ,"Broadleaf Evergreen Tropical Tree" \

   333                ,"Broadleaf Evergreen Temperate Tree" \

   334                ,"Broadleaf Deciduous Tropical Tree" \

   335                ,"Broadleaf Deciduous Temperate Tree" \

   336 ;              ,"Broadleaf Deciduous Boreal Tree" \

   337 ;              ,"Broadleaf Evergreen Shrub" \

   338                ,"Broadleaf Deciduous Temperate Shrub" \

   339                ,"Broadleaf Deciduous Boreal Shrub" \

   340                ,"C3 Arctic Grass" \

   341                ,"C3 Non-Arctic Grass" \

   342                ,"C4 Grass" \

   343                ,"Corn" \

   344 ;              ,"Wheat" \

   345                ,"All Biome" \

   346                /)

   347   nrow = dimsizes(row_head)

   349 ; arrays to be passed to table.

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

   352  do i=0,nrow-2

   353   text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12)

   354   text4(i,1) = sprintf("%.1f",NPP1_g(i))

   355   text4(i,2) = sprintf("%.1f",VegC_g(i))

   356   text4(i,3) = sprintf("%.1f",LiCwC_g(i))

   357   text4(i,4) = sprintf("%.1f",SoilC_g(i))

   358   text4(i,5) = sprintf("%.1f",NPP_ratio(i))

   359   text4(i,6) = sprintf("%.1f",NEE2_g(i))

   360  end do

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

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

   363   text4(nrow-1,2) = "-"

   365 ;-------------------------------------------------------

   366 ; create html table1

   368   header_text = "<H1>NEE and Carbon Stocks and Fluxes:  Model "+model_name+"</H1>"

   370   header = (/"<HTML>" \

   371             ,"<HEAD>" \

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

   373             ,"</HEAD>" \

   374             ,header_text \

   375             /)

   376   footer = "</HTML>"

   378   table_header = (/ \

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

   380        ,"<tr>" \

   381        ,"   <th bgcolor=DDDDDD >Biome Type</th>" \

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

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

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

   385        ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \

   386        ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \

   387        ,"   <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \

   388        ,"   <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \

   389        ,"</tr>" \

   390        /)

   391   table_footer = "</table>"

   392   row_header = "<tr>"

   393   row_footer = "</tr>"

   395   lines = new(50000,string)

   396   nline = 0

   398   set_line(lines,nline,header)

   399   set_line(lines,nline,table_header)

   401 ;----------------------------

   402 ;row of table

   404   do n = 0,nrow-2

   405      set_line(lines,nline,row_header)

   407      txt0  = row_head(n)

   408      txt1  = text4(n,0)

   409      txt2  = text4(n,1)

   410      txt3  = text4(n,2)

   411      txt4  = text4(n,3)

   412      txt5  = text4(n,4)

   413      txt6  = text4(n,5)

   414      txt7  = text4(n,6)

   416      set_line(lines,nline,"<th>"+txt0+"</th>")

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

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

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

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

   421      set_line(lines,nline,"<th>"+txt5+"</th>")

   422      set_line(lines,nline,"<th>"+txt6+"</th>")

   423      set_line(lines,nline,"<th>"+txt7+"</th>")

   425      set_line(lines,nline,row_footer)

   426   end do

   427 ;----------------------------

   428   set_line(lines,nline,table_footer)

   429   set_line(lines,nline,footer)

   431 ; Now write to an HTML file.

   433   output_html = "table_carbon_sink1.html"

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

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

   437     asciiwrite(output_html,lines(idx))

   438   else

   439    print ("error?")

   440   end if

   442   delete (col_head)

   443   delete (row_head)

   444   delete (text4)

   446 end

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