1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/beta/07.landfrac.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,591 @@
1.4 +;********************************************************
1.5 +; take into account landfrac
1.6 +; note: landfrac from lnd_T42.nc
1.7 +; <= lnd_diag_4.0 has correct landfrac
1.8 +; lnd_diag_3.1 has wrong landfrac
1.9 +;
1.10 +; using model biome
1.11 +;
1.12 +; required command line input parameters:
1.13 +; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
1.14 +;
1.15 +;**************************************************************
1.16 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.17 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.18 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.19 +;**************************************************************
1.20 +procedure set_line(lines:string,nline:integer,newlines:string)
1.21 +begin
1.22 +; add line to ascci/html file
1.23 +
1.24 + nnewlines = dimsizes(newlines)
1.25 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.26 + print("set_line: bad index, not setting anything.")
1.27 + return
1.28 + end if
1.29 + lines(nline:nline+nnewlines-1) = newlines
1.30 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.31 + nline = nline + nnewlines
1.32 + return
1.33 +end
1.34 +;**************************************************************
1.35 +; Main code.
1.36 +begin
1.37 +
1.38 + plot_type = "ps"
1.39 + plot_type_new = "png"
1.40 +
1.41 +;************************************************
1.42 +; read data: model
1.43 +;************************************************
1.44 + co2_i = 283.1878
1.45 + co2_f = 364.1252
1.46 +
1.47 + model_grid = "T42"
1.48 +
1.49 +;model_name_i = "i01.07cn"
1.50 +;model_name_f = "i01.10cn"
1.51 +
1.52 + model_name_i = "i01.07casa"
1.53 + model_name_f = "i01.10casa"
1.54 +
1.55 + model_name = model_name_f
1.56 +
1.57 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.58 + film_i = model_name_i + "_1990-2004_ANN_climo.nc"
1.59 + film_f = model_name_f + "_1990-2004_ANN_climo.nc"
1.60 +
1.61 + fm_i = addfile (dirm+film_i,"r")
1.62 + fm_f = addfile (dirm+film_f,"r")
1.63 +
1.64 + xm = fm_f->lon
1.65 + ym = fm_f->lat
1.66 +
1.67 + npp_i = fm_i->NPP
1.68 + npp_f = fm_f->NPP
1.69 +
1.70 + delete (fm_i)
1.71 + delete (fm_f)
1.72 +
1.73 +;Units for these variables are:
1.74 +;npp_i: g C/m^2/s
1.75 +
1.76 + nsec_per_year = 60*60*24*365
1.77 +
1.78 + npp_i = npp_i * nsec_per_year
1.79 + npp_f = npp_f * nsec_per_year
1.80 +
1.81 +;--------------------------------------------------
1.82 +;get landfrac data
1.83 +
1.84 + dirm= "/fis/cgd/cseg/people/jeff/surface_data/"
1.85 + film_l = "lnd_T42.nc"
1.86 + fm_l = addfile (dirm+film_l,"r")
1.87 +
1.88 + landfrac = fm_l->landfrac
1.89 +
1.90 +;npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:)
1.91 +;npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:)
1.92 +
1.93 + npp_i = npp_i * conform(npp_i, landfrac, (/1,2/))
1.94 + npp_f = npp_f * conform(npp_f, landfrac, (/1,2/))
1.95 +
1.96 +;===================================================
1.97 +; read data: observed at stations
1.98 +;===================================================
1.99 +
1.100 + station = (/"DukeFACE" \
1.101 + ,"AspenFACE" \
1.102 + ,"ORNL-FACE" \
1.103 + ,"POP-EUROFACE" \
1.104 + /)
1.105 +
1.106 + lat_ob = (/ 35.58, 45.40, 35.54, 42.22/)
1.107 + lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
1.108 + lon_ob = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
1.109 +;print (lon_ob)
1.110 +
1.111 + n_sta = dimsizes(station)
1.112 + beta_4_ob = new((/n_sta/),float)
1.113 + beta_4_ob = 0.60
1.114 +
1.115 +;===================================================
1.116 +; get model data at station
1.117 +;===================================================
1.118 +
1.119 + npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_ob,lat_ob,0)
1.120 +
1.121 + npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_ob,lat_ob,0)
1.122 +
1.123 +;print (npp_i_4)
1.124 +;print (npp_f_4)
1.125 +
1.126 +;============================
1.127 +;compute beta_4
1.128 +;============================
1.129 +
1.130 + beta_4 = new((/n_sta/),float)
1.131 +
1.132 + beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
1.133 +
1.134 + beta_4_avg = avg(beta_4)
1.135 +
1.136 +;print (beta_4)
1.137 +;print (beta_4_avg)
1.138 +
1.139 +;M_beta = abs((beta_4_avg/beta_4_ob) - 1.)* 3.
1.140 +
1.141 + bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))
1.142 + M_beta = (1. - (bias/n_sta))*3.
1.143 +
1.144 + print (M_beta)
1.145 +
1.146 +;=========================
1.147 +; for html table - station
1.148 +;=========================
1.149 +
1.150 + output_html = "table_station.html"
1.151 +
1.152 +; column (not including header column)
1.153 +
1.154 + col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
1.155 +
1.156 + ncol = dimsizes(col_head)
1.157 +
1.158 +; row (not including header row)
1.159 + row_head = (/"DukeFACE" \
1.160 + ,"AspenFACE" \
1.161 + ,"ORNL-FACE" \
1.162 + ,"POP-EUROFACE" \
1.163 + ,"All Station" \
1.164 + /)
1.165 + nrow = dimsizes(row_head)
1.166 +
1.167 +; arrays to be passed to table.
1.168 + text4 = new ((/nrow, ncol/),string )
1.169 +
1.170 + do i=0,nrow-2
1.171 + text4(i,0) = sprintf("%.1f",lat_ob(i))
1.172 + text4(i,1) = sprintf("%.1f",lon_ob(i))
1.173 + text4(i,2) = sprintf("%.1f",co2_i)
1.174 + text4(i,3) = sprintf("%.1f",co2_f)
1.175 + text4(i,4) = sprintf("%.1f",npp_i_4(0,i))
1.176 + text4(i,5) = sprintf("%.1f",npp_f_4(0,i))
1.177 + text4(i,6) = sprintf("%.2f",beta_4(i))
1.178 + text4(i,7) = "-"
1.179 + end do
1.180 + text4(nrow-1,0) = "-"
1.181 + text4(nrow-1,1) = "-"
1.182 + text4(nrow-1,2) = "-"
1.183 + text4(nrow-1,3) = "-"
1.184 + text4(nrow-1,4) = "-"
1.185 + text4(nrow-1,5) = "-"
1.186 + text4(nrow-1,6) = sprintf("%.2f",beta_4_avg)
1.187 + text4(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob))
1.188 +
1.189 +;-----------
1.190 +; html table
1.191 +;-----------
1.192 +
1.193 + header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
1.194 +
1.195 + header = (/"<HTML>" \
1.196 + ,"<HEAD>" \
1.197 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.198 + ,"</HEAD>" \
1.199 + ,header_text \
1.200 + /)
1.201 + footer = "</HTML>"
1.202 +
1.203 + table_header = (/ \
1.204 + "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
1.205 + ,"<tr>" \
1.206 + ," <th bgcolor=DDDDDD >Station</th>" \
1.207 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.208 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.209 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.210 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.211 + ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
1.212 + ," <th bgcolor=DDDDDD >"+col_head(5)+"</th>" \
1.213 + ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
1.214 + ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
1.215 + ,"</tr>" \
1.216 + /)
1.217 + table_footer = "</table>"
1.218 + row_header = "<tr>"
1.219 + row_footer = "</tr>"
1.220 +
1.221 + lines = new(50000,string)
1.222 + nline = 0
1.223 +
1.224 + set_line(lines,nline,header)
1.225 + set_line(lines,nline,table_header)
1.226 +;-----------------------------------------------
1.227 +;row of table
1.228 +
1.229 + do n = 0,nrow-1
1.230 + set_line(lines,nline,row_header)
1.231 +
1.232 + txt1 = row_head(n)
1.233 + txt2 = text4(n,0)
1.234 + txt3 = text4(n,1)
1.235 + txt4 = text4(n,2)
1.236 + txt5 = text4(n,3)
1.237 + txt6 = text4(n,4)
1.238 + txt7 = text4(n,5)
1.239 + txt8 = text4(n,6)
1.240 + txt9 = text4(n,7)
1.241 +
1.242 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.243 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.244 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.245 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.246 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.247 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.248 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.249 + set_line(lines,nline,"<th>"+txt8+"</th>")
1.250 + set_line(lines,nline,"<th>"+txt9+"</th>")
1.251 +
1.252 + set_line(lines,nline,row_footer)
1.253 + end do
1.254 +;-----------------------------------------------
1.255 + set_line(lines,nline,table_footer)
1.256 + set_line(lines,nline,footer)
1.257 +
1.258 +; Now write to an HTML file.
1.259 + idx = ind(.not.ismissing(lines))
1.260 + if(.not.any(ismissing(idx))) then
1.261 + asciiwrite(output_html,lines(idx))
1.262 + else
1.263 + print ("error?")
1.264 + end if
1.265 +
1.266 + delete (col_head)
1.267 + delete (row_head)
1.268 + delete (text4)
1.269 + delete (table_header)
1.270 + delete (idx)
1.271 +
1.272 +;------------------------------------------------
1.273 +; read biome data: model
1.274 +
1.275 + biome_name_mod = "Model PFT Class"
1.276 +
1.277 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.278 + film = "class_pft_"+model_grid+".nc"
1.279 +
1.280 + fm = addfile(dirm+film,"r")
1.281 +
1.282 + classmod = fm->CLASS_PFT
1.283 +
1.284 + delete (fm)
1.285 +
1.286 +; model data has 17 land-type classes
1.287 +
1.288 + nclass_mod = 17
1.289 +
1.290 +;------------------------------------------------
1.291 +; read biome data: observed
1.292 +
1.293 + biome_name_ob = "MODIS LandCover"
1.294 +
1.295 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
1.296 + filo = "land_class_"+model_grid+".nc"
1.297 +
1.298 + fo = addfile(diro+filo,"r")
1.299 +
1.300 + classob = tofloat(fo->LAND_CLASS)
1.301 +
1.302 + delete (fo)
1.303 +
1.304 +; input observed data has 20 land-type classes
1.305 +
1.306 + nclass_ob = 20
1.307 +
1.308 +;********************************************************************
1.309 +; use land-type class to bin the data in equally spaced ranges
1.310 +;********************************************************************
1.311 +
1.312 +; using observed biome class
1.313 +; nclass = nclass_ob
1.314 +; using model biome class
1.315 + nclass = nclass_mod
1.316 +
1.317 + nclassn = nclass + 1
1.318 + range = fspan(0,nclassn-1,nclassn)
1.319 +; print (range)
1.320 +
1.321 +; Use this range information to grab all the values in a
1.322 +; particular range, and then take an average.
1.323 +
1.324 + nr = dimsizes(range)
1.325 + nx = nr-1
1.326 + xvalues = new((/2,nx/),float)
1.327 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.328 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.329 + dx4 = dx/4 ; 1/4 of the range
1.330 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.331 +
1.332 +;==============================
1.333 +; put data into bins
1.334 +;==============================
1.335 +
1.336 +; using observed biome class
1.337 +; base_1D = ndtooned(classob)
1.338 +; using model biome class
1.339 + base_1D = ndtooned(classmod)
1.340 +
1.341 + data1_1D = ndtooned(npp_i)
1.342 + data2_1D = ndtooned(npp_f)
1.343 +
1.344 +; output
1.345 +
1.346 + yvalues = new((/2,nx/),float)
1.347 + count = new((/2,nx/),float)
1.348 +
1.349 + do nd=0,1
1.350 +
1.351 +; See if we are doing data1 (nd=0) or data2 (nd=1).
1.352 +
1.353 + base = base_1D
1.354 +
1.355 + if(nd.eq.0) then
1.356 + data = data1_1D
1.357 + else
1.358 + data = data2_1D
1.359 + end if
1.360 +
1.361 +; Loop through each range, using base.
1.362 +
1.363 + do i=0,nr-2
1.364 + if (i.ne.(nr-2)) then
1.365 +; print("")
1.366 +; print("In range ["+range(i)+","+range(i+1)+")")
1.367 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.368 + else
1.369 +; print("")
1.370 +; print("In range ["+range(i)+",)")
1.371 + idx = ind(base.ge.range(i))
1.372 + end if
1.373 +
1.374 +; Calculate average
1.375 +
1.376 + if(.not.any(ismissing(idx))) then
1.377 + yvalues(nd,i) = avg(data(idx))
1.378 + count(nd,i) = dimsizes(idx)
1.379 + else
1.380 + yvalues(nd,i) = yvalues@_FillValue
1.381 + count(nd,i) = 0
1.382 + end if
1.383 +
1.384 +;#############################################################
1.385 +;using observed biome class:
1.386 +;
1.387 +; set the following 4 classes to _FillValue:
1.388 +; Water Bodies(0), Urban and Build-Up(13),
1.389 +; Permenant Snow and Ice(15), Unclassified(17)
1.390 +
1.391 +; if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
1.392 +; yvalues(nd,i) = yvalues@_FillValue
1.393 +; count(nd,i) = 0
1.394 +; end if
1.395 +;#############################################################
1.396 +
1.397 +;#############################################################
1.398 +;using model biome class:
1.399 +;
1.400 +; set the following 4 classes to _FillValue:
1.401 +; (3)Needleleaf Deciduous Boreal Tree,
1.402 +; (8)Broadleaf Deciduous Boreal Tree,
1.403 +; (9)Broadleaf Evergreen Shrub,
1.404 +; (16)Wheat
1.405 +
1.406 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.407 + yvalues(nd,i) = yvalues@_FillValue
1.408 + count(nd,i) = 0
1.409 + end if
1.410 +;#############################################################
1.411 +
1.412 +; print(nd + ": " + count(nd,i) + " points, avg = " + yvalues(nd,i))
1.413 +
1.414 +; Clean up for next time in loop.
1.415 +
1.416 + delete(idx)
1.417 + end do
1.418 +
1.419 + delete(data)
1.420 + end do
1.421 +
1.422 +;============================
1.423 +;compute beta
1.424 +;============================
1.425 +
1.426 + u = yvalues(0,:)
1.427 + v = yvalues(1,:)
1.428 + u_count = count(0,:)
1.429 + v_count = count(1,:)
1.430 +
1.431 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.432 +
1.433 + uu = u(good)
1.434 + vv = v(good)
1.435 + uu_count = u_count(good)
1.436 + vv_count = v_count(good)
1.437 +
1.438 + n_biome = dimsizes(uu)
1.439 + beta_biome = new((/n_biome/),float)
1.440 +
1.441 + beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
1.442 +
1.443 +;beta_biome_avg = avg(beta_biome)
1.444 + beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)
1.445 +
1.446 +;print (beta_biome_avg)
1.447 +
1.448 +;===========================
1.449 +; for html table - biome
1.450 +;===========================
1.451 +
1.452 + output_html = "table_biome.html"
1.453 +
1.454 +; column (not including header column)
1.455 +
1.456 + col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
1.457 +
1.458 + ncol = dimsizes(col_head)
1.459 +
1.460 +; row (not including header row)
1.461 +
1.462 +;----------------------------------------------------
1.463 +; using observed biome class:
1.464 +;
1.465 +; row_head = (/"Evergreen Needleleaf Forests" \
1.466 +; ,"Evergreen Broadleaf Forests" \
1.467 +; ,"Deciduous Needleleaf Forest" \
1.468 +; ,"Deciduous Broadleaf Forests" \
1.469 +; ,"Mixed Forests" \
1.470 +; ,"Closed Bushlands" \
1.471 +; ,"Open Bushlands" \
1.472 +; ,"Woody Savannas (S. Hem.)" \
1.473 +; ,"Savannas (S. Hem.)" \
1.474 +; ,"Grasslands" \
1.475 +; ,"Permanent Wetlands" \
1.476 +; ,"Croplands" \
1.477 +; ,"Cropland/Natural Vegetation Mosaic" \
1.478 +; ,"Barren or Sparsely Vegetated" \
1.479 +; ,"Woody Savannas (N. Hem.)" \
1.480 +; ,"Savannas (N. Hem.)" \
1.481 +; ,"All Biome" \
1.482 +; /)
1.483 +
1.484 +;----------------------------------------------------
1.485 +; using model biome class:
1.486 +;
1.487 + row_head = (/"Not Vegetated" \
1.488 + ,"Needleleaf Evergreen Temperate Tree" \
1.489 + ,"Needleleaf Evergreen Boreal Tree" \
1.490 +; ,"Needleleaf Deciduous Boreal Tree" \
1.491 + ,"Broadleaf Evergreen Tropical Tree" \
1.492 + ,"Broadleaf Evergreen Temperate Tree" \
1.493 + ,"Broadleaf Deciduous Tropical Tree" \
1.494 + ,"Broadleaf Deciduous Temperate Tree" \
1.495 +; ,"Broadleaf Deciduous Boreal Tree" \
1.496 +; ,"Broadleaf Evergreen Shrub" \
1.497 + ,"Broadleaf Deciduous Temperate Shrub" \
1.498 + ,"Broadleaf Deciduous Boreal Shrub" \
1.499 + ,"C3 Arctic Grass" \
1.500 + ,"C3 Non-Arctic Grass" \
1.501 + ,"C4 Grass" \
1.502 + ,"Corn" \
1.503 +; ,"Wheat" \
1.504 + ,"All Biome" \
1.505 + /)
1.506 +
1.507 + nrow = dimsizes(row_head)
1.508 +
1.509 +; arrays to be passed to table.
1.510 + text4 = new ((/nrow, ncol/),string )
1.511 +
1.512 + do i=0,nrow-2
1.513 + text4(i,0) = sprintf("%.1f",co2_i)
1.514 + text4(i,1) = sprintf("%.1f",co2_f)
1.515 + text4(i,2) = sprintf("%.1f",uu(i))
1.516 + text4(i,3) = sprintf("%.1f",vv(i))
1.517 + text4(i,4) = sprintf("%.2f",beta_biome(i))
1.518 + end do
1.519 + text4(nrow-1,0) = "-"
1.520 + text4(nrow-1,1) = "-"
1.521 + text4(nrow-1,2) = "-"
1.522 + text4(nrow-1,3) = "-"
1.523 + text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
1.524 +
1.525 +;**************************************************
1.526 +; html table
1.527 +;**************************************************
1.528 +
1.529 + header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
1.530 +
1.531 + header = (/"<HTML>" \
1.532 + ,"<HEAD>" \
1.533 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.534 + ,"</HEAD>" \
1.535 + ,header_text \
1.536 + /)
1.537 + footer = "</HTML>"
1.538 +
1.539 + table_header = (/ \
1.540 + "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
1.541 + ,"<tr>" \
1.542 + ," <th bgcolor=DDDDDD >Biome Class</th>" \
1.543 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.544 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.545 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.546 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.547 + ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
1.548 + ,"</tr>" \
1.549 + /)
1.550 + table_footer = "</table>"
1.551 + row_header = "<tr>"
1.552 + row_footer = "</tr>"
1.553 +
1.554 + lines = new(50000,string)
1.555 + nline = 0
1.556 +
1.557 + set_line(lines,nline,header)
1.558 + set_line(lines,nline,table_header)
1.559 +;-----------------------------------------------
1.560 +;row of table
1.561 +
1.562 + do n = 0,nrow-1
1.563 + set_line(lines,nline,row_header)
1.564 +
1.565 + txt1 = row_head(n)
1.566 + txt2 = text4(n,0)
1.567 + txt3 = text4(n,1)
1.568 + txt4 = text4(n,2)
1.569 + txt5 = text4(n,3)
1.570 + txt6 = text4(n,4)
1.571 +
1.572 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.573 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.574 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.575 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.576 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.577 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.578 +
1.579 + set_line(lines,nline,row_footer)
1.580 + end do
1.581 +;-----------------------------------------------
1.582 + set_line(lines,nline,table_footer)
1.583 + set_line(lines,nline,footer)
1.584 +
1.585 +; Now write to an HTML file.
1.586 + idx = ind(.not.ismissing(lines))
1.587 + if(.not.any(ismissing(idx))) then
1.588 + asciiwrite(output_html,lines(idx))
1.589 + else
1.590 + print ("error?")
1.591 + end if
1.592 +
1.593 +end
1.594 +