carbon_sink/12.table1+2.ncl
changeset 0 0c6405ab2ff4
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/carbon_sink/12.table1+2.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,593 @@
     1.4 +;********************************************************
     1.5 +;using model biome vlass
     1.6 +;
     1.7 +; required command line input parameters:
     1.8 +;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
     1.9 +;
    1.10 +; histogram normalized by rain and compute correleration
    1.11 +;**************************************************************
    1.12 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
    1.13 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
    1.14 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    1.15 +;**************************************************************
    1.16 +procedure set_line(lines:string,nline:integer,newlines:string) 
    1.17 +begin
    1.18 +; add line to ascci/html file
    1.19 +    
    1.20 +  nnewlines = dimsizes(newlines)
    1.21 +  if(nline+nnewlines-1.ge.dimsizes(lines))
    1.22 +    print("set_line: bad index, not setting anything.") 
    1.23 +    return
    1.24 +  end if 
    1.25 +  lines(nline:nline+nnewlines-1) = newlines
    1.26 +;  print ("lines = " + lines(nline:nline+nnewlines-1))
    1.27 +  nline = nline + nnewlines
    1.28 +  return 
    1.29 +end
    1.30 +;**************************************************************
    1.31 +; Main code.
    1.32 +begin
    1.33 + 
    1.34 + plot_type     = "ps"
    1.35 + plot_type_new = "png"
    1.36 +
    1.37 +;************************************************
    1.38 +; model name and grid       
    1.39 +;************************************************
    1.40 +
    1.41 + model_grid = "T42"
    1.42 +
    1.43 + model_name  = "cn"
    1.44 + model_name1 = "i01.06cn"
    1.45 + model_name2 = "i01.10cn"
    1.46 +
    1.47 +;---------------------------------------------------
    1.48 +; get biome data: model
    1.49 +
    1.50 +  biome_name_mod = "Model PFT Class"
    1.51 +
    1.52 +  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
    1.53 +  film = "class_pft_"+model_grid+".nc"
    1.54 +
    1.55 +  fm = addfile(dirm+film,"r")
    1.56 + 
    1.57 +  classmod = fm->CLASS_PFT
    1.58 +
    1.59 +  delete (fm)
    1.60 +
    1.61 +; model data has 17 land-type classes
    1.62 +
    1.63 +  nclass_mod = 17
    1.64 +
    1.65 +;--------------------------------------------------
    1.66 +; get model data: landfrac and area
    1.67 +
    1.68 + dirm_l= "/fis/cgd/cseg/people/jeff/surface_data/" 
    1.69 + film_l = "lnd_T42.nc"
    1.70 + fm_l   = addfile (dirm_l+film_l,"r")
    1.71 +  
    1.72 + landfrac = fm_l->landfrac
    1.73 + area     = fm_l->area
    1.74 +
    1.75 +; change area from km**2 to m**2
    1.76 +  area = area * 1.e6             
    1.77 +;---------------------------------------------------
    1.78 +; read data: model, group 1
    1.79 +
    1.80 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
    1.81 + film = model_name1 + "_1980-2004_ANN_climo.nc"
    1.82 + fm   = addfile (dirm+film,"r")
    1.83 + 
    1.84 + NPP1   = fm->NPP
    1.85 +
    1.86 + leafc  = fm->LEAFC
    1.87 + woodc  = fm->WOODC
    1.88 + frootc = fm->FROOTC
    1.89 + VegC   = leafc
    1.90 + VegC   = leafc + woodc + frootc
    1.91 +
    1.92 + litterc = fm->LITTERC
    1.93 + cwdc    = fm->CWDC
    1.94 + LiCwC   = litterc
    1.95 + LiCwC   = litterc + cwdc
    1.96 +
    1.97 + SoilC   = fm->SOILC
    1.98 +
    1.99 + delete (fm)
   1.100 +;--------------------------------------------------- 
   1.101 +; read data: model, group 2
   1.102 +
   1.103 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
   1.104 + film = model_name2 + "_1990-2004_ANN_climo.nc"
   1.105 + fm   = addfile (dirm+film,"r")
   1.106 +
   1.107 + NPP2   = fm->NPP 
   1.108 + NEE2   = fm->NEE 
   1.109 +
   1.110 +;---------------------------------------------------
   1.111 +; Units for these variables are:
   1.112 +
   1.113 +;NPP1: g C/m^2/s
   1.114 +;NPP2: g C/m^2/s
   1.115 +;NEE2: g C/m^2/s
   1.116 +
   1.117 +;VegC:  g C/m^2
   1.118 +;LiCwC: g C/m^2
   1.119 +;SoilC: g C/m^2
   1.120 +
   1.121 + nsec_per_year = 60*60*24*365
   1.122 +
   1.123 +; change unit to g C/m^2/year
   1.124 +  
   1.125 + NPP1 = NPP1 *  nsec_per_year
   1.126 + NPP2 = NPP2 *  nsec_per_year
   1.127 + NEE2 = NEE2 *  nsec_per_year 
   1.128 +
   1.129 +;---------------------------------------------------
   1.130 +; take into account landfrac
   1.131 +
   1.132 + area(:,:)    = area(:,:)    * landfrac(:,:)
   1.133 + NPP1(0,:,:)  = NPP1(0,:,:)  * landfrac(:,:)
   1.134 + VegC(0,:,:)  = VegC(0,:,:)  * landfrac(:,:)
   1.135 + LiCwC(0,:,:) = LiCwC(0,:,:) * landfrac(:,:)
   1.136 + SoilC(0,:,:) = SoilC(0,:,:) * landfrac(:,:)
   1.137 + NPP2(0,:,:)  = NPP2(0,:,:)  * landfrac(:,:)
   1.138 + NEE2(0,:,:)  = NEE2(0,:,:)  * landfrac(:,:)
   1.139 +
   1.140 + data_n = 7
   1.141 +                
   1.142 +;*******************************************************************
   1.143 +; Calculate "nice" bins for binning the data in equally spaced ranges
   1.144 +;********************************************************************
   1.145 +
   1.146 +; using model biome class
   1.147 +  nclass      = nclass_mod
   1.148 +
   1.149 +  range       = fspan(0,nclass,nclass+1)
   1.150 +
   1.151 +; print (range)
   1.152 +; Use this range information to grab all the values in a
   1.153 +; particular range, and then take an average.
   1.154 +
   1.155 +  nx           = dimsizes(range) - 1
   1.156 +
   1.157 +;==============================
   1.158 +; put data into bins
   1.159 +;==============================
   1.160 +
   1.161 +; using observed biome class
   1.162 +; base  = ndtooned(classob)
   1.163 +; using model biome class
   1.164 +  base  = ndtooned(classmod)
   1.165 +
   1.166 +; output
   1.167 +
   1.168 +  yvalues = new((/data_n,nx/),float)
   1.169 +  count   = new((/data_n,nx/),float)
   1.170 +
   1.171 +  do n = 0,data_n-1
   1.172 +
   1.173 +    if(n.eq.0) then
   1.174 +      data = ndtooned(area)
   1.175 +    end if
   1.176 +
   1.177 +    if(n.eq.1) then
   1.178 +      data = ndtooned(NPP1)
   1.179 +    end if
   1.180 +
   1.181 +    if(n.eq.2) then
   1.182 +      data = ndtooned(VegC)
   1.183 +    end if
   1.184 +
   1.185 +    if(n.eq.3) then
   1.186 +      data = ndtooned(LiCwC)
   1.187 +    end if
   1.188 +
   1.189 +    if(n.eq.4) then
   1.190 +      data = ndtooned(SoilC)
   1.191 +    end if
   1.192 +
   1.193 +    if(n.eq.5) then
   1.194 +      data = ndtooned(NPP2)
   1.195 +    end if
   1.196 +
   1.197 +    if(n.eq.6) then
   1.198 +      data = ndtooned(NEE2)
   1.199 +    end if
   1.200 +
   1.201 +; Loop through each range, using base.
   1.202 +
   1.203 +    do i=0,nx-1
   1.204 +      if (i.ne.(nx-1)) then
   1.205 +;        print("")
   1.206 +;        print("In range ["+range(i)+","+range(i+1)+")")
   1.207 +         idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
   1.208 +      else
   1.209 +;        print("")
   1.210 +;        print("In range ["+range(i)+",)")
   1.211 +         idx = ind(base.ge.range(i))
   1.212 +      end if
   1.213 +
   1.214 +;     Calculate average
   1.215 + 
   1.216 +      if(.not.any(ismissing(idx))) then
   1.217 +        if (n.eq.0) then 
   1.218 +           yvalues(n,i) = sum(data(idx))
   1.219 +        else 
   1.220 +           yvalues(n,i) = avg(data(idx))
   1.221 +        end if
   1.222 +
   1.223 +        count(n,i)   = dimsizes(idx)
   1.224 +      else
   1.225 +        yvalues(n,i) = yvalues@_FillValue
   1.226 +        count(n,i)   = 0
   1.227 +      end if
   1.228 +
   1.229 +;#############################################################
   1.230 +; using model biome class:
   1.231 +;     set the following 4 classes to _FillValue:
   1.232 +;     (3)Needleleaf Deciduous Boreal Tree,
   1.233 +;     (8)Broadleaf Deciduous Boreal Tree,
   1.234 +;     (9)Broadleaf Evergreen Shrub,
   1.235 +;     (16)Wheat
   1.236 +
   1.237 +      if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
   1.238 +         yvalues(n,i) = yvalues@_FillValue
   1.239 +         count(n,i)   = 0
   1.240 +      end if
   1.241 +;#############################################################  
   1.242 +
   1.243 +;     print(n + ": " + count + " points, avg = " + yvalues(n,i))
   1.244 +
   1.245 +      delete(idx)
   1.246 +    end do 
   1.247 +
   1.248 +    delete(data)
   1.249 +  end do
   1.250 +
   1.251 +  delete (base)
   1.252 +  delete (area)
   1.253 +  delete (NPP1)
   1.254 +  delete (VegC)
   1.255 +  delete (LiCwC)
   1.256 +  delete (SoilC)
   1.257 +  delete (NPP2)
   1.258 +  delete (NEE2)
   1.259 +
   1.260 +;----------------------------------------------------------------
   1.261 +; data for table1
   1.262 +
   1.263 + good = ind(.not.ismissing(yvalues(5,:)) .and. .not.ismissing(yvalues(1,:)))
   1.264 +;print (good)
   1.265 +
   1.266 + w = yvalues(0,:)
   1.267 + area_g = w(good)
   1.268 +
   1.269 + w = yvalues(1,:)
   1.270 + NPP1_g = w(good)
   1.271 +
   1.272 + w = yvalues(2,:)
   1.273 + VegC_g = w(good)
   1.274 +
   1.275 + w = yvalues(3,:)
   1.276 + LiCwC_g = w(good)
   1.277 +
   1.278 + w = yvalues(4,:)
   1.279 + SoilC_g = w(good)
   1.280 +
   1.281 + w = yvalues(5,:)
   1.282 + NPP2_g = w(good)
   1.283 +
   1.284 + w = yvalues(6,:)
   1.285 + NEE2_g = w(good) 
   1.286 +
   1.287 + n_biome = dimsizes(NPP1_g)
   1.288 +
   1.289 + NPP1_t    = new((/n_biome/),float)
   1.290 + VegC_t    = new((/n_biome/),float)
   1.291 + LiCwC_t   = new((/n_biome/),float)
   1.292 + SoilC_t   = new((/n_biome/),float)
   1.293 + NEE2_t    = new((/n_biome/),float)
   1.294 + NPP_ratio = new((/n_biome/),float)
   1.295 +
   1.296 + NPP_ratio = NPP2_g/NPP1_g
   1.297 +
   1.298 +;-----------------------------------------------------------------
   1.299 +; data for table2
   1.300 +
   1.301 +; change unit from g to Pg (Peta gram)
   1.302 + factor_unit = 1.e-15
   1.303 +
   1.304 + NPP1_t    = NPP1_g  * area_g * factor_unit
   1.305 + VegC_t    = VegC_g  * area_g * factor_unit
   1.306 + LiCwC_t   = LiCwC_g * area_g * factor_unit
   1.307 + SoilC_t   = SoilC_g * area_g * factor_unit
   1.308 + NEE2_t    = NEE2_g  * area_g * factor_unit
   1.309 +
   1.310 + print (NPP1_t)
   1.311 + 
   1.312 +;-------------------------------------------------------------
   1.313 +; html table1 data
   1.314 +
   1.315 +; column (not including header column)
   1.316 +
   1.317 +  col_head  = (/"Area (1.e12m2)" \
   1.318 +               ,"NPP (gC/m2/yr)" \
   1.319 +               ,"VegC (gC/m2)" \
   1.320 +               ,"Litter+CWD (gC/m2)" \
   1.321 +               ,"SoilC (gC/m2)" \
   1.322 +               ,"NPP_ratio" \
   1.323 +               ,"NEE (gC/m2/yr)" \
   1.324 +               /)
   1.325 +
   1.326 +  ncol = dimsizes(col_head)
   1.327 +
   1.328 +; row (not including header row)                   
   1.329 +
   1.330 +; using model biome class:  
   1.331 +  row_head  = (/"Not Vegetated" \
   1.332 +               ,"Needleleaf Evergreen Temperate Tree" \
   1.333 +               ,"Needleleaf Evergreen Boreal Tree" \
   1.334 +;              ,"Needleleaf Deciduous Boreal Tree" \
   1.335 +               ,"Broadleaf Evergreen Tropical Tree" \
   1.336 +               ,"Broadleaf Evergreen Temperate Tree" \
   1.337 +               ,"Broadleaf Deciduous Tropical Tree" \
   1.338 +               ,"Broadleaf Deciduous Temperate Tree" \
   1.339 +;              ,"Broadleaf Deciduous Boreal Tree" \
   1.340 +;              ,"Broadleaf Evergreen Shrub" \
   1.341 +               ,"Broadleaf Deciduous Temperate Shrub" \
   1.342 +               ,"Broadleaf Deciduous Boreal Shrub" \
   1.343 +               ,"C3 Arctic Grass" \
   1.344 +               ,"C3 Non-Arctic Grass" \
   1.345 +               ,"C4 Grass" \
   1.346 +               ,"Corn" \
   1.347 +;              ,"Wheat" \                      
   1.348 +               ,"All Biome" \                
   1.349 +               /)  
   1.350 +  nrow = dimsizes(row_head)                  
   1.351 +
   1.352 +; arrays to be passed to table. 
   1.353 +  text4 = new ((/nrow, ncol/),string )
   1.354 + 
   1.355 + do i=0,nrow-2
   1.356 +  text4(i,0) = sprintf("%.1f",area_g(i)*1.e-12)
   1.357 +  text4(i,1) = sprintf("%.1f",NPP1_g(i))
   1.358 +  text4(i,2) = sprintf("%.1f",VegC_g(i))
   1.359 +  text4(i,3) = sprintf("%.1f",LiCwC_g(i))
   1.360 +  text4(i,4) = sprintf("%.1f",SoilC_g(i))
   1.361 +  text4(i,5) = sprintf("%.1f",NPP_ratio(i))
   1.362 +  text4(i,6) = sprintf("%.1f",NEE2_g(i))
   1.363 + end do
   1.364 +  text4(nrow-1,0) = "-"
   1.365 +  text4(nrow-1,1) = "-"
   1.366 +  text4(nrow-1,2) = "-"
   1.367 +
   1.368 +;-------------------------------------------------------
   1.369 +; create html table1
   1.370 +
   1.371 +  header_text = "<H1>NEE and Carbon Stocks and Fluxes:  Model "+model_name+"</H1>" 
   1.372 +
   1.373 +  header = (/"<HTML>" \
   1.374 +            ,"<HEAD>" \
   1.375 +            ,"<TITLE>CLAMP metrics</TITLE>" \
   1.376 +            ,"</HEAD>" \
   1.377 +            ,header_text \
   1.378 +            /) 
   1.379 +  footer = "</HTML>"
   1.380 +
   1.381 +  table_header = (/ \
   1.382 +        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
   1.383 +       ,"<tr>" \
   1.384 +       ,"   <th bgcolor=DDDDDD >Biome Type</th>" \
   1.385 +       ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
   1.386 +       ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
   1.387 +       ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
   1.388 +       ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
   1.389 +       ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
   1.390 +       ,"   <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
   1.391 +       ,"   <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
   1.392 +       ,"</tr>" \
   1.393 +       /)
   1.394 +  table_footer = "</table>"
   1.395 +  row_header = "<tr>"
   1.396 +  row_footer = "</tr>"
   1.397 +
   1.398 +  lines = new(50000,string)
   1.399 +  nline = 0
   1.400 +
   1.401 +  set_line(lines,nline,header)
   1.402 +  set_line(lines,nline,table_header)
   1.403 +
   1.404 +;----------------------------
   1.405 +;row of table
   1.406 +
   1.407 +  do n = 0,nrow-2
   1.408 +     set_line(lines,nline,row_header)
   1.409 +
   1.410 +     txt0  = row_head(n)
   1.411 +     txt1  = text4(n,0)
   1.412 +     txt2  = text4(n,1)
   1.413 +     txt3  = text4(n,2)
   1.414 +     txt4  = text4(n,3)
   1.415 +     txt5  = text4(n,4)
   1.416 +     txt6  = text4(n,5)
   1.417 +     txt7  = text4(n,6)
   1.418 +
   1.419 +     set_line(lines,nline,"<th>"+txt0+"</th>")
   1.420 +     set_line(lines,nline,"<th>"+txt1+"</th>")
   1.421 +     set_line(lines,nline,"<th>"+txt2+"</th>")
   1.422 +     set_line(lines,nline,"<th>"+txt3+"</th>")
   1.423 +     set_line(lines,nline,"<th>"+txt4+"</th>")
   1.424 +     set_line(lines,nline,"<th>"+txt5+"</th>")
   1.425 +     set_line(lines,nline,"<th>"+txt6+"</th>")
   1.426 +     set_line(lines,nline,"<th>"+txt7+"</th>")
   1.427 +
   1.428 +     set_line(lines,nline,row_footer)
   1.429 +  end do
   1.430 +;----------------------------
   1.431 +  set_line(lines,nline,table_footer)
   1.432 +  set_line(lines,nline,footer) 
   1.433 +
   1.434 +; Now write to an HTML file.
   1.435 +
   1.436 +  output_html = "table_carbon_sink1.html"
   1.437 +
   1.438 +  idx = ind(.not.ismissing(lines))
   1.439 +  if(.not.any(ismissing(idx))) then
   1.440 +    asciiwrite(output_html,lines(idx))
   1.441 +  else
   1.442 +   print ("error?")
   1.443 +  end if
   1.444 +
   1.445 +  delete (idx)
   1.446 +
   1.447 +  delete (col_head)
   1.448 +  delete (row_head)
   1.449 +  delete (text4)
   1.450 +  delete (table_header)
   1.451 +
   1.452 +;-----------------------------------------------------------------
   1.453 +; html table2 data
   1.454 +
   1.455 +; column (not including header column)
   1.456 +
   1.457 +  col_head  = (/"NPP (PgC/yr)" \
   1.458 +               ,"VegC (PgC)" \
   1.459 +               ,"Litter+CWD (PgC)" \
   1.460 +               ,"SoilC (PgC)" \
   1.461 +               ,"NEE (PgC/yr)" \
   1.462 +               ,"NPP timeseries" \
   1.463 +               ,"NEE timeseries" \
   1.464 +               ,"Fire timeseries" \
   1.465 +               /)
   1.466 +
   1.467 +  ncol = dimsizes(col_head)
   1.468 +
   1.469 +; row (not including header row)                   
   1.470 +
   1.471 +; using model biome class:  
   1.472 +  row_head  = (/"Not Vegetated" \
   1.473 +               ,"Needleleaf Evergreen Temperate Tree" \
   1.474 +               ,"Needleleaf Evergreen Boreal Tree" \
   1.475 +;              ,"Needleleaf Deciduous Boreal Tree" \
   1.476 +               ,"Broadleaf Evergreen Tropical Tree" \
   1.477 +               ,"Broadleaf Evergreen Temperate Tree" \
   1.478 +               ,"Broadleaf Deciduous Tropical Tree" \
   1.479 +               ,"Broadleaf Deciduous Temperate Tree" \
   1.480 +;              ,"Broadleaf Deciduous Boreal Tree" \
   1.481 +;              ,"Broadleaf Evergreen Shrub" \
   1.482 +               ,"Broadleaf Deciduous Temperate Shrub" \
   1.483 +               ,"Broadleaf Deciduous Boreal Shrub" \
   1.484 +               ,"C3 Arctic Grass" \
   1.485 +               ,"C3 Non-Arctic Grass" \
   1.486 +               ,"C4 Grass" \
   1.487 +               ,"Corn" \
   1.488 +;              ,"Wheat" \                      
   1.489 +               ,"All Biome" \                
   1.490 +               /)  
   1.491 +  nrow = dimsizes(row_head)                  
   1.492 +
   1.493 +; arrays to be passed to table. 
   1.494 +  text4 = new ((/nrow, ncol/),string )
   1.495 + 
   1.496 + do i=0,nrow-2
   1.497 +  text4(i,0) = sprintf("%.1f",NPP1_t(i))
   1.498 +  text4(i,1) = sprintf("%.1f",VegC_t(i))
   1.499 +  text4(i,2) = sprintf("%.1f",LiCwC_t(i))
   1.500 +  text4(i,3) = sprintf("%.1f",SoilC_t(i))
   1.501 +  text4(i,4) = sprintf("%.1f",NEE2_t(i))
   1.502 +  text4(i,5) = "-"
   1.503 +  text4(i,6) = "-"
   1.504 +  text4(i,7) = "-"
   1.505 + end do
   1.506 +  text4(nrow-1,0) = sprintf("%.1f",sum(NPP1_t))
   1.507 +  text4(nrow-1,1) = sprintf("%.1f",sum(VegC_t))
   1.508 +  text4(nrow-1,2) = sprintf("%.1f",sum(LiCwC_t))
   1.509 +  text4(nrow-1,3) = sprintf("%.1f",sum(SoilC_t))
   1.510 +  text4(nrow-1,4) = sprintf("%.1f",sum(NEE2_t))
   1.511 +  text4(nrow-1,5) = "-"
   1.512 +  text4(nrow-1,6) = "-"
   1.513 +  text4(nrow-1,7) = "-"
   1.514 +
   1.515 +;**************************************************
   1.516 +; create html table1
   1.517 +;**************************************************
   1.518 +
   1.519 +  header_text = "<H1>NEE and Carbon Stocks and Fluxes (per biome):  Model "+model_name+"</H1>" 
   1.520 +
   1.521 +  header = (/"<HTML>" \
   1.522 +            ,"<HEAD>" \
   1.523 +            ,"<TITLE>CLAMP metrics</TITLE>" \
   1.524 +            ,"</HEAD>" \
   1.525 +            ,header_text \
   1.526 +            /) 
   1.527 +  footer = "</HTML>"
   1.528 +
   1.529 +  table_header = (/ \
   1.530 +        "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
   1.531 +       ,"<tr>" \
   1.532 +       ,"   <th bgcolor=DDDDDD >Biome Type</th>" \
   1.533 +       ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
   1.534 +       ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
   1.535 +       ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
   1.536 +       ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
   1.537 +       ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
   1.538 +       ,"   <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
   1.539 +       ,"   <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
   1.540 +       ,"   <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
   1.541 +       ,"</tr>" \
   1.542 +       /)
   1.543 +  table_footer = "</table>"
   1.544 +  row_header = "<tr>"
   1.545 +  row_footer = "</tr>"
   1.546 +
   1.547 +  lines = new(50000,string)
   1.548 +  nline = 0
   1.549 +
   1.550 +  set_line(lines,nline,header)
   1.551 +  set_line(lines,nline,table_header)
   1.552 +;-----------------------------------------------
   1.553 +;row of table
   1.554 +
   1.555 +  do n = 0,nrow-1
   1.556 +     set_line(lines,nline,row_header)
   1.557 +
   1.558 +     txt0  = row_head(n)
   1.559 +     txt1  = text4(n,0)
   1.560 +     txt2  = text4(n,1)
   1.561 +     txt3  = text4(n,2)
   1.562 +     txt4  = text4(n,3)
   1.563 +     txt5  = text4(n,4)
   1.564 +     txt6  = text4(n,5)
   1.565 +     txt7  = text4(n,6)
   1.566 +     txt8  = text4(n,7)
   1.567 +
   1.568 +     set_line(lines,nline,"<th>"+txt0+"</th>")
   1.569 +     set_line(lines,nline,"<th>"+txt1+"</th>")
   1.570 +     set_line(lines,nline,"<th>"+txt2+"</th>")
   1.571 +     set_line(lines,nline,"<th>"+txt3+"</th>")
   1.572 +     set_line(lines,nline,"<th>"+txt4+"</th>")
   1.573 +     set_line(lines,nline,"<th>"+txt5+"</th>")
   1.574 +     set_line(lines,nline,"<th>"+txt6+"</th>")
   1.575 +     set_line(lines,nline,"<th>"+txt7+"</th>")
   1.576 +     set_line(lines,nline,"<th>"+txt8+"</th>")
   1.577 +
   1.578 +     set_line(lines,nline,row_footer)
   1.579 +  end do
   1.580 +;-----------------------------------------------
   1.581 +  set_line(lines,nline,table_footer)
   1.582 +  set_line(lines,nline,footer) 
   1.583 +
   1.584 +; Now write to an HTML file.
   1.585 +
   1.586 +  output_html = "table_carbon_sink2.html"
   1.587 +
   1.588 +  idx = ind(.not.ismissing(lines))
   1.589 +  if(.not.any(ismissing(idx))) then
   1.590 +    asciiwrite(output_html,lines(idx))
   1.591 +  else
   1.592 +   print ("error?")
   1.593 +  end if
   1.594 +
   1.595 +end
   1.596 +