beta/07.landfrac.ncl
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:f560a0e6dba7
       
     1 ;********************************************************
       
     2 ; take into account landfrac
       
     3 ; note: landfrac from lnd_T42.nc
       
     4 ;       <= lnd_diag_4.0 has correct landfrac
       
     5 ;          lnd_diag_3.1 has wrong   landfrac  
       
     6 ;
       
     7 ; using model biome
       
     8 ;
       
     9 ; required command line input parameters:
       
    10 ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
       
    11 ;
       
    12 ;**************************************************************
       
    13 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
       
    14 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
       
    15 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
       
    16 ;**************************************************************
       
    17 procedure set_line(lines:string,nline:integer,newlines:string) 
       
    18 begin
       
    19 ; add line to ascci/html file
       
    20     
       
    21   nnewlines = dimsizes(newlines)
       
    22   if(nline+nnewlines-1.ge.dimsizes(lines))
       
    23     print("set_line: bad index, not setting anything.") 
       
    24     return
       
    25   end if 
       
    26   lines(nline:nline+nnewlines-1) = newlines
       
    27 ;  print ("lines = " + lines(nline:nline+nnewlines-1))
       
    28   nline = nline + nnewlines
       
    29   return 
       
    30 end
       
    31 ;**************************************************************
       
    32 ; Main code.
       
    33 begin
       
    34  
       
    35  plot_type     = "ps"
       
    36  plot_type_new = "png"
       
    37 
       
    38 ;************************************************
       
    39 ; read data: model       
       
    40 ;************************************************
       
    41  co2_i = 283.1878
       
    42  co2_f = 364.1252
       
    43 
       
    44  model_grid = "T42"
       
    45 
       
    46 ;model_name_i = "i01.07cn"
       
    47 ;model_name_f = "i01.10cn"
       
    48 
       
    49  model_name_i = "i01.07casa"
       
    50  model_name_f = "i01.10casa"
       
    51 
       
    52  model_name = model_name_f
       
    53 
       
    54  dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
    55  film_i = model_name_i + "_1990-2004_ANN_climo.nc"
       
    56  film_f = model_name_f + "_1990-2004_ANN_climo.nc"
       
    57 
       
    58  fm_i   = addfile (dirm+film_i,"r")
       
    59  fm_f   = addfile (dirm+film_f,"r")
       
    60   
       
    61  xm     = fm_f->lon  
       
    62  ym     = fm_f->lat
       
    63 
       
    64  npp_i  = fm_i->NPP
       
    65  npp_f  = fm_f->NPP
       
    66 
       
    67  delete (fm_i)
       
    68  delete (fm_f)
       
    69 
       
    70 ;Units for these variables are:
       
    71 ;npp_i: g C/m^2/s
       
    72 
       
    73  nsec_per_year = 60*60*24*365
       
    74   
       
    75  npp_i = npp_i *  nsec_per_year
       
    76  npp_f = npp_f *  nsec_per_year
       
    77 
       
    78 ;--------------------------------------------------
       
    79 ;get landfrac data
       
    80 
       
    81  dirm= "/fis/cgd/cseg/people/jeff/surface_data/" 
       
    82  film_l = "lnd_T42.nc"
       
    83  fm_l   = addfile (dirm+film_l,"r")
       
    84   
       
    85  landfrac = fm_l->landfrac
       
    86 
       
    87 ;npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:)
       
    88 ;npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:)
       
    89 
       
    90  npp_i = npp_i * conform(npp_i, landfrac, (/1,2/))
       
    91  npp_f = npp_f * conform(npp_f, landfrac, (/1,2/))
       
    92     
       
    93 ;===================================================
       
    94 ; read data: observed at stations
       
    95 ;===================================================
       
    96 
       
    97  station = (/"DukeFACE" \
       
    98             ,"AspenFACE" \
       
    99             ,"ORNL-FACE" \
       
   100             ,"POP-EUROFACE" \
       
   101             /)
       
   102 
       
   103  lat_ob = (/ 35.58,  45.40,  35.54, 42.22/)
       
   104  lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
       
   105  lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
       
   106 ;print (lon_ob)
       
   107 
       
   108  n_sta  = dimsizes(station)
       
   109  beta_4_ob = new((/n_sta/),float)
       
   110  beta_4_ob = 0.60
       
   111 
       
   112 ;===================================================
       
   113 ; get model data at station 
       
   114 ;===================================================
       
   115 
       
   116  npp_i_4  =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
       
   117 
       
   118  npp_f_4  =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
       
   119 
       
   120 ;print (npp_i_4)
       
   121 ;print (npp_f_4)
       
   122 
       
   123 ;============================
       
   124 ;compute beta_4
       
   125 ;============================
       
   126 
       
   127  beta_4 = new((/n_sta/),float)
       
   128 
       
   129  beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
       
   130 
       
   131  beta_4_avg = avg(beta_4)
       
   132 
       
   133 ;print (beta_4)
       
   134 ;print (beta_4_avg)
       
   135 
       
   136 ;M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 3.
       
   137 
       
   138  bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob))) 
       
   139  M_beta  = (1. - (bias/n_sta))*3.
       
   140  
       
   141  print (M_beta)
       
   142 
       
   143 ;=========================
       
   144 ; for html table - station
       
   145 ;=========================
       
   146 
       
   147   output_html = "table_station.html"
       
   148 
       
   149 ; column (not including header column)
       
   150 
       
   151   col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
       
   152 
       
   153   ncol = dimsizes(col_head)
       
   154 
       
   155 ; row (not including header row)
       
   156   row_head = (/"DukeFACE" \
       
   157               ,"AspenFACE" \
       
   158               ,"ORNL-FACE" \
       
   159               ,"POP-EUROFACE" \
       
   160               ,"All Station" \                
       
   161               /)  
       
   162   nrow = dimsizes(row_head)                  
       
   163 
       
   164 ; arrays to be passed to table. 
       
   165   text4 = new ((/nrow, ncol/),string )
       
   166 
       
   167  do i=0,nrow-2
       
   168   text4(i,0) = sprintf("%.1f",lat_ob(i))
       
   169   text4(i,1) = sprintf("%.1f",lon_ob(i))
       
   170   text4(i,2) = sprintf("%.1f",co2_i)
       
   171   text4(i,3) = sprintf("%.1f",co2_f)
       
   172   text4(i,4) = sprintf("%.1f",npp_i_4(0,i))
       
   173   text4(i,5) = sprintf("%.1f",npp_f_4(0,i))
       
   174   text4(i,6) = sprintf("%.2f",beta_4(i))
       
   175   text4(i,7) = "-"
       
   176  end do
       
   177   text4(nrow-1,0) = "-"
       
   178   text4(nrow-1,1) = "-"
       
   179   text4(nrow-1,2) = "-"
       
   180   text4(nrow-1,3) = "-"
       
   181   text4(nrow-1,4) = "-"
       
   182   text4(nrow-1,5) = "-"
       
   183   text4(nrow-1,6) = sprintf("%.2f",beta_4_avg)
       
   184   text4(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob))
       
   185 
       
   186 ;-----------
       
   187 ; html table
       
   188 ;-----------
       
   189 
       
   190   header_text = "<H1>Beta Factor: Model "+model_name+"</H1>" 
       
   191 
       
   192   header = (/"<HTML>" \
       
   193             ,"<HEAD>" \
       
   194             ,"<TITLE>CLAMP metrics</TITLE>" \
       
   195             ,"</HEAD>" \
       
   196             ,header_text \
       
   197             /) 
       
   198   footer = "</HTML>"
       
   199 
       
   200   table_header = (/ \
       
   201         "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
       
   202        ,"<tr>" \
       
   203        ,"   <th bgcolor=DDDDDD >Station</th>" \
       
   204        ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
       
   205        ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
       
   206        ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
       
   207        ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
       
   208        ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
       
   209        ,"   <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
       
   210        ,"   <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
       
   211        ,"   <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
       
   212        ,"</tr>" \
       
   213        /)
       
   214   table_footer = "</table>"
       
   215   row_header = "<tr>"
       
   216   row_footer = "</tr>"
       
   217 
       
   218   lines = new(50000,string)
       
   219   nline = 0
       
   220 
       
   221   set_line(lines,nline,header)
       
   222   set_line(lines,nline,table_header)
       
   223 ;-----------------------------------------------
       
   224 ;row of table
       
   225 
       
   226   do n = 0,nrow-1
       
   227      set_line(lines,nline,row_header)
       
   228 
       
   229      txt1  = row_head(n)
       
   230      txt2  = text4(n,0)
       
   231      txt3  = text4(n,1)
       
   232      txt4  = text4(n,2)
       
   233      txt5  = text4(n,3)
       
   234      txt6  = text4(n,4)
       
   235      txt7  = text4(n,5)
       
   236      txt8  = text4(n,6)
       
   237      txt9  = text4(n,7)
       
   238 
       
   239      set_line(lines,nline,"<th>"+txt1+"</th>")
       
   240      set_line(lines,nline,"<th>"+txt2+"</th>")
       
   241      set_line(lines,nline,"<th>"+txt3+"</th>")
       
   242      set_line(lines,nline,"<th>"+txt4+"</th>")
       
   243      set_line(lines,nline,"<th>"+txt5+"</th>")
       
   244      set_line(lines,nline,"<th>"+txt6+"</th>")
       
   245      set_line(lines,nline,"<th>"+txt7+"</th>")
       
   246      set_line(lines,nline,"<th>"+txt8+"</th>")
       
   247      set_line(lines,nline,"<th>"+txt9+"</th>")
       
   248 
       
   249      set_line(lines,nline,row_footer)
       
   250   end do
       
   251 ;-----------------------------------------------
       
   252   set_line(lines,nline,table_footer)
       
   253   set_line(lines,nline,footer) 
       
   254 
       
   255 ; Now write to an HTML file.
       
   256   idx = ind(.not.ismissing(lines))
       
   257   if(.not.any(ismissing(idx))) then
       
   258     asciiwrite(output_html,lines(idx))
       
   259   else
       
   260    print ("error?")
       
   261   end if
       
   262 
       
   263   delete (col_head)
       
   264   delete (row_head)
       
   265   delete (text4)
       
   266   delete (table_header)
       
   267   delete (idx)
       
   268 
       
   269 ;------------------------------------------------
       
   270 ; read biome data: model
       
   271 
       
   272   biome_name_mod = "Model PFT Class"
       
   273 
       
   274   dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
   275   film = "class_pft_"+model_grid+".nc"
       
   276 
       
   277   fm = addfile(dirm+film,"r")
       
   278  
       
   279   classmod = fm->CLASS_PFT               
       
   280 
       
   281   delete (fm)
       
   282 
       
   283 ; model data has 17 land-type classes
       
   284 
       
   285   nclass_mod = 17
       
   286 
       
   287 ;------------------------------------------------
       
   288 ; read biome data: observed
       
   289 
       
   290   biome_name_ob = "MODIS LandCover"
       
   291 
       
   292   diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
       
   293   filo = "land_class_"+model_grid+".nc"
       
   294 
       
   295   fo = addfile(diro+filo,"r")
       
   296  
       
   297   classob = tofloat(fo->LAND_CLASS)               
       
   298 
       
   299   delete (fo)
       
   300 
       
   301 ; input observed data has 20 land-type classes
       
   302 
       
   303   nclass_ob = 20
       
   304                 
       
   305 ;********************************************************************
       
   306 ; use land-type class to bin the data in equally spaced ranges
       
   307 ;********************************************************************
       
   308 
       
   309 ; using observed biome class  
       
   310 ; nclass      = nclass_ob
       
   311 ; using model biome class
       
   312   nclass      = nclass_mod
       
   313 
       
   314   nclassn     = nclass + 1
       
   315   range       = fspan(0,nclassn-1,nclassn)
       
   316 ; print (range)
       
   317 
       
   318 ; Use this range information to grab all the values in a
       
   319 ; particular range, and then take an average.
       
   320 
       
   321   nr           = dimsizes(range)
       
   322   nx           = nr-1
       
   323   xvalues      = new((/2,nx/),float)
       
   324   xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
       
   325   dx           = xvalues(0,1) - xvalues(0,0)       ; range width
       
   326   dx4          = dx/4                              ; 1/4 of the range
       
   327   xvalues(1,:) = xvalues(0,:) - dx/5.
       
   328 
       
   329 ;==============================
       
   330 ; put data into bins
       
   331 ;==============================
       
   332 
       
   333 ; using observed biome class
       
   334 ; base_1D  = ndtooned(classob)
       
   335 ; using model biome class
       
   336   base_1D  = ndtooned(classmod)
       
   337 
       
   338   data1_1D = ndtooned(npp_i)
       
   339   data2_1D = ndtooned(npp_f)
       
   340 
       
   341 ; output
       
   342 
       
   343   yvalues = new((/2,nx/),float)
       
   344   count   = new((/2,nx/),float)
       
   345 
       
   346   do nd=0,1
       
   347 
       
   348 ;   See if we are doing data1 (nd=0) or data2 (nd=1).
       
   349 
       
   350     base = base_1D
       
   351 
       
   352     if(nd.eq.0) then
       
   353       data = data1_1D
       
   354     else
       
   355       data = data2_1D
       
   356     end if
       
   357 
       
   358 ; Loop through each range, using base.
       
   359 
       
   360     do i=0,nr-2
       
   361       if (i.ne.(nr-2)) then
       
   362 ;        print("")
       
   363 ;        print("In range ["+range(i)+","+range(i+1)+")")
       
   364          idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
       
   365       else
       
   366 ;        print("")
       
   367 ;        print("In range ["+range(i)+",)")
       
   368          idx = ind(base.ge.range(i))
       
   369       end if
       
   370 
       
   371 ;     Calculate average 
       
   372 
       
   373       if(.not.any(ismissing(idx))) then
       
   374         yvalues(nd,i) = avg(data(idx))
       
   375         count(nd,i)   = dimsizes(idx)
       
   376       else
       
   377         yvalues(nd,i) = yvalues@_FillValue
       
   378         count(nd,i)   = 0
       
   379       end if
       
   380 
       
   381 ;#############################################################
       
   382 ;using observed biome class:
       
   383 ; 
       
   384 ;     set the following 4 classes to _FillValue:
       
   385 ;     Water Bodies(0), Urban and Build-Up(13),
       
   386 ;     Permenant Snow and Ice(15), Unclassified(17)
       
   387 
       
   388 ;     if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
       
   389 ;        yvalues(nd,i) = yvalues@_FillValue
       
   390 ;        count(nd,i)   = 0
       
   391 ;     end if
       
   392 ;############################################################# 
       
   393 
       
   394 ;#############################################################
       
   395 ;using model biome class:
       
   396 ;
       
   397 ;     set the following 4 classes to _FillValue:
       
   398 ;     (3)Needleleaf Deciduous Boreal Tree,
       
   399 ;     (8)Broadleaf Deciduous Boreal Tree,
       
   400 ;     (9)Broadleaf Evergreen Shrub,
       
   401 ;     (16)Wheat
       
   402 
       
   403       if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
       
   404          yvalues(nd,i) = yvalues@_FillValue
       
   405          count(nd,i)   = 0
       
   406       end if
       
   407 ;############################################################# 
       
   408 
       
   409 ;     print(nd + ": " + count(nd,i) + " points, avg = " + yvalues(nd,i))
       
   410 
       
   411 ; Clean up for next time in loop.
       
   412 
       
   413       delete(idx)
       
   414     end do
       
   415 
       
   416     delete(data)
       
   417   end do
       
   418 
       
   419 ;============================
       
   420 ;compute beta
       
   421 ;============================
       
   422 
       
   423  u       = yvalues(0,:)
       
   424  v       = yvalues(1,:)
       
   425  u_count = count(0,:)
       
   426  v_count = count(1,:)
       
   427 
       
   428  good = ind(.not.ismissing(u) .and. .not.ismissing(v))
       
   429 
       
   430  uu       = u(good)
       
   431  vv       = v(good)
       
   432  uu_count = u_count(good)
       
   433  vv_count = v_count(good) 
       
   434 
       
   435  n_biome = dimsizes(uu)
       
   436  beta_biome = new((/n_biome/),float)
       
   437 
       
   438  beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
       
   439 
       
   440 ;beta_biome_avg = avg(beta_biome)
       
   441  beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)
       
   442   
       
   443 ;print (beta_biome_avg)
       
   444 
       
   445 ;===========================
       
   446 ; for html table - biome
       
   447 ;===========================
       
   448 
       
   449   output_html = "table_biome.html"
       
   450 
       
   451 ; column (not including header column)
       
   452 
       
   453   col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
       
   454 
       
   455   ncol = dimsizes(col_head)
       
   456 
       
   457 ; row (not including header row)
       
   458 
       
   459 ;----------------------------------------------------
       
   460 ; using observed biome class:
       
   461 ;  
       
   462 ; row_head  = (/"Evergreen Needleleaf Forests" \
       
   463 ;              ,"Evergreen Broadleaf Forests" \
       
   464 ;              ,"Deciduous Needleleaf Forest" \
       
   465 ;              ,"Deciduous Broadleaf Forests" \
       
   466 ;              ,"Mixed Forests" \                      
       
   467 ;              ,"Closed Bushlands" \                   
       
   468 ;              ,"Open Bushlands" \                     
       
   469 ;              ,"Woody Savannas (S. Hem.)" \           
       
   470 ;              ,"Savannas (S. Hem.)" \                 
       
   471 ;              ,"Grasslands" \                         
       
   472 ;              ,"Permanent Wetlands" \                 
       
   473 ;              ,"Croplands" \                                           
       
   474 ;              ,"Cropland/Natural Vegetation Mosaic" \             
       
   475 ;              ,"Barren or Sparsely Vegetated" \                             
       
   476 ;              ,"Woody Savannas (N. Hem.)" \           
       
   477 ;              ,"Savannas (N. Hem.)" \
       
   478 ;              ,"All Biome" \                
       
   479 ;              /)
       
   480 
       
   481 ;----------------------------------------------------
       
   482 ; using model biome class:
       
   483 ;  
       
   484   row_head  = (/"Not Vegetated" \
       
   485                ,"Needleleaf Evergreen Temperate Tree" \
       
   486                ,"Needleleaf Evergreen Boreal Tree" \
       
   487 ;              ,"Needleleaf Deciduous Boreal Tree" \
       
   488                ,"Broadleaf Evergreen Tropical Tree" \
       
   489                ,"Broadleaf Evergreen Temperate Tree" \
       
   490                ,"Broadleaf Deciduous Tropical Tree" \
       
   491                ,"Broadleaf Deciduous Temperate Tree" \
       
   492 ;              ,"Broadleaf Deciduous Boreal Tree" \
       
   493 ;              ,"Broadleaf Evergreen Shrub" \
       
   494                ,"Broadleaf Deciduous Temperate Shrub" \
       
   495                ,"Broadleaf Deciduous Boreal Shrub" \
       
   496                ,"C3 Arctic Grass" \
       
   497                ,"C3 Non-Arctic Grass" \
       
   498                ,"C4 Grass" \
       
   499                ,"Corn" \
       
   500 ;              ,"Wheat" \                      
       
   501                ,"All Biome" \                
       
   502                /)  
       
   503 
       
   504   nrow = dimsizes(row_head)                  
       
   505 
       
   506 ; arrays to be passed to table. 
       
   507   text4 = new ((/nrow, ncol/),string )
       
   508  
       
   509  do i=0,nrow-2
       
   510   text4(i,0) = sprintf("%.1f",co2_i)
       
   511   text4(i,1) = sprintf("%.1f",co2_f)
       
   512   text4(i,2) = sprintf("%.1f",uu(i))
       
   513   text4(i,3) = sprintf("%.1f",vv(i))
       
   514   text4(i,4) = sprintf("%.2f",beta_biome(i))
       
   515  end do
       
   516   text4(nrow-1,0) = "-"
       
   517   text4(nrow-1,1) = "-"
       
   518   text4(nrow-1,2) = "-"
       
   519   text4(nrow-1,3) = "-"
       
   520   text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
       
   521 
       
   522 ;**************************************************
       
   523 ; html table
       
   524 ;**************************************************
       
   525 
       
   526   header_text = "<H1>Beta Factor: Model "+model_name+"</H1>" 
       
   527 
       
   528   header = (/"<HTML>" \
       
   529             ,"<HEAD>" \
       
   530             ,"<TITLE>CLAMP metrics</TITLE>" \
       
   531             ,"</HEAD>" \
       
   532             ,header_text \
       
   533             /) 
       
   534   footer = "</HTML>"
       
   535 
       
   536   table_header = (/ \
       
   537         "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
       
   538        ,"<tr>" \
       
   539        ,"   <th bgcolor=DDDDDD >Biome Class</th>" \
       
   540        ,"   <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
       
   541        ,"   <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
       
   542        ,"   <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
       
   543        ,"   <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
       
   544        ,"   <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
       
   545        ,"</tr>" \
       
   546        /)
       
   547   table_footer = "</table>"
       
   548   row_header = "<tr>"
       
   549   row_footer = "</tr>"
       
   550 
       
   551   lines = new(50000,string)
       
   552   nline = 0
       
   553 
       
   554   set_line(lines,nline,header)
       
   555   set_line(lines,nline,table_header)
       
   556 ;-----------------------------------------------
       
   557 ;row of table
       
   558 
       
   559   do n = 0,nrow-1
       
   560      set_line(lines,nline,row_header)
       
   561 
       
   562      txt1  = row_head(n)
       
   563      txt2  = text4(n,0)
       
   564      txt3  = text4(n,1)
       
   565      txt4  = text4(n,2)
       
   566      txt5  = text4(n,3)
       
   567      txt6  = text4(n,4)
       
   568 
       
   569      set_line(lines,nline,"<th>"+txt1+"</th>")
       
   570      set_line(lines,nline,"<th>"+txt2+"</th>")
       
   571      set_line(lines,nline,"<th>"+txt3+"</th>")
       
   572      set_line(lines,nline,"<th>"+txt4+"</th>")
       
   573      set_line(lines,nline,"<th>"+txt5+"</th>")
       
   574      set_line(lines,nline,"<th>"+txt6+"</th>")
       
   575 
       
   576      set_line(lines,nline,row_footer)
       
   577   end do
       
   578 ;-----------------------------------------------
       
   579   set_line(lines,nline,table_footer)
       
   580   set_line(lines,nline,footer) 
       
   581 
       
   582 ; Now write to an HTML file.
       
   583   idx = ind(.not.ismissing(lines))
       
   584   if(.not.any(ismissing(idx))) then
       
   585     asciiwrite(output_html,lines(idx))
       
   586   else
       
   587    print ("error?")
       
   588   end if
       
   589 
       
   590 end
       
   591