all/03.co2.ncl
changeset 0 0c6405ab2ff4
child 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:115c14a99b5c
       
     1 ; ***********************************************
       
     2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
       
     3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
       
     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
       
     5 load dirscript + "taylor_diagram.ncl"
       
     6 ;************************************************
       
     7 procedure set_line(lines:string,nline:integer,newlines:string) 
       
     8 begin
       
     9 ; add line to ascci/html file
       
    10     
       
    11   nnewlines = dimsizes(newlines)
       
    12   if(nline+nnewlines-1.ge.dimsizes(lines))
       
    13     print("set_line: bad index, not setting anything.") 
       
    14     return
       
    15   end if 
       
    16   lines(nline:nline+nnewlines-1) = newlines
       
    17 ;  print ("lines = " + lines(nline:nline+nnewlines-1))
       
    18   nline = nline + nnewlines
       
    19   return 
       
    20 end
       
    21 ;****************************************************************************
       
    22 
       
    23 begin
       
    24 
       
    25   plot_type = "ps"
       
    26   plot_type_new = "png"
       
    27 
       
    28 ;-----------------------------------------------------
       
    29 ; edit table.html of current model for movel1_vs_model2
       
    30 
       
    31   if (isvar("compare")) then
       
    32      html_name2 = compare+"/table.html"  
       
    33      html_new2  = html_name2 +".new"
       
    34   end if
       
    35 
       
    36 ;------------------------------------------------------
       
    37 ; edit table.html for current model
       
    38 
       
    39   html_name = model_name+"/table.html"  
       
    40   html_new  = html_name +".new"
       
    41 
       
    42 ;------------------------------------------------------
       
    43 ; read model data
       
    44 
       
    45   fm    = addfile(dirm+film3,"r")
       
    46 
       
    47   x     = fm->CO2
       
    48   xi    = fm->lon
       
    49   yi    = fm->lat
       
    50 
       
    51   delete (fm)
       
    52 
       
    53   xdim  = dimsizes(x)
       
    54   nlev  = xdim(1)
       
    55   y     = x(:,0,:,:)
       
    56   
       
    57 ; get co2 at the lowest level
       
    58   y     = x(:,nlev-1,:,:)
       
    59 
       
    60 ; change to unit of observed (u mol/mol)
       
    61 ; Model_units [=] kgCO2 / kgDryAir
       
    62 ; 28.966 = molecular weight of dry air
       
    63 ; 44.       = molecular weight of CO2
       
    64 ; u mol = 1e-6 mol
       
    65 
       
    66   factor = (28.966/44.) * 1e6
       
    67   y      = y * factor
       
    68 
       
    69   y@_FillValue = 1.e36
       
    70   y@units      = "u mol/mol"
       
    71 
       
    72 ;************************************************
       
    73 ; read data: observed
       
    74 ;************************************************
       
    75   diri  = diro + "co2/"
       
    76   fili  = "co2_globalView_98.nc"
       
    77   g     = addfile (diri+fili,"r")
       
    78 
       
    79   val   = g->CO2_SEAS  
       
    80   lon   = g->LON 
       
    81   lat   = g->LAT
       
    82   sta   = chartostring(g->STATION)
       
    83  
       
    84   delete (g)
       
    85 
       
    86   ncase = dimsizes(lat)
       
    87 ;**************************************************************
       
    88 ; get only the lowest level at each station 
       
    89 ;**************************************************************
       
    90   lat_tmp = lat
       
    91   lat_tmp@_FillValue = 1.e+36
       
    92  
       
    93   do n = 0,ncase-1
       
    94      if (.not. ismissing(lat_tmp(n))) then 
       
    95         indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
       
    96         if (dimsizes(indexes) .gt. 1) then
       
    97            lat_tmp(indexes(1:)) = lat_tmp@_FillValue
       
    98         end if
       
    99         delete (indexes)
       
   100      end if
       
   101   end do
       
   102 
       
   103   indexes = ind(.not. ismissing(lat_tmp))
       
   104  
       
   105   lat_ob = lat(indexes)
       
   106   lon_ob = lon(indexes)
       
   107   val_ob = val(indexes,:)
       
   108 ;************************************************************
       
   109 ; interpolate model data into observed station
       
   110 ; note: model is 0-360E, 90S-90N
       
   111 ;************************************************************
       
   112 ; to be able to handle observation at (-89.98,-24.80)
       
   113   yi(0) = -90.
       
   114 
       
   115   i = ind(lon_ob .lt. 0.)
       
   116   lon_ob(i) = lon_ob(i) + 360.  
       
   117 
       
   118   yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
       
   119 
       
   120   val_model = yo(pts|:,time|:)
       
   121   val_model_0 = val_model
       
   122 ;************************************************************
       
   123 ; remove annual mean
       
   124 ;************************************************************
       
   125   val_model = val_model - conform(val_model,dim_avg(val_model),0)
       
   126 
       
   127 ;*******************************************************************
       
   128 ; res for station line plot
       
   129 ;*******************************************************************
       
   130 ; for x-axis in xyplot
       
   131   mon = ispan(1,12,1)
       
   132   mon@long_name = "month"
       
   133 
       
   134   res                   = True                      ; plot mods desired
       
   135   res@xyLineThicknesses = (/2.0,2.0,2.0/)           ; make 2nd lines thicker
       
   136   res@xyLineColors      = (/"red","blue"/)  ; change line color
       
   137 
       
   138 ; Add a boxed legend using the more simple method
       
   139   res@pmLegendDisplayMode    = "Always"
       
   140 ; res@pmLegendWidthF         = 0.1
       
   141   res@pmLegendWidthF         = 0.08
       
   142   res@pmLegendHeightF        = 0.06
       
   143 ; res@pmLegendOrthogonalPosF = -1.17
       
   144 ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
       
   145   res@pmLegendOrthogonalPosF = -0.30  ;(downward)
       
   146 
       
   147 ; res@pmLegendParallelPosF   =  0.18
       
   148   res@pmLegendParallelPosF   =  0.23  ;(rightward)
       
   149 
       
   150 ; res@lgPerimOn             = False
       
   151   res@lgLabelFontHeightF     = 0.015
       
   152   res@xyExplicitLegendLabels = (/model_name,"observed"/)
       
   153 ;************************************************************
       
   154 ; number of latitude zone
       
   155 ;************************************************************
       
   156   nzone = 4
       
   157 
       
   158 ; saving data for zone
       
   159 ; number of rows for zone table (with data)
       
   160   nrow_zone = nzone 
       
   161 
       
   162 ; number of columns for zone table
       
   163   ncol_zone = 7
       
   164 
       
   165   text = new((/nrow_zone,ncol_zone/),string)
       
   166 
       
   167 do z = 0,nzone-1
       
   168 
       
   169   if (z .eq. 0) then 
       
   170      zone = "60N-90N" 
       
   171      score_max = 5.0
       
   172      ind_z = ind(lat_ob .ge. 60.)
       
   173   end if
       
   174 
       
   175   if (z .eq. 1) then 
       
   176      zone = "30N-60N" 
       
   177      score_max = 5.0
       
   178      ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
       
   179   end if
       
   180 
       
   181   if (z .eq. 2) then 
       
   182      zone = "EQ-30N"
       
   183      score_max = 5.0
       
   184      ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
       
   185   end if
       
   186 
       
   187   if (z .eq. 3) then 
       
   188      zone = "90S-EQ" 
       
   189      score_max = 5.0
       
   190      ind_z = ind(lat_ob .lt. 0. )
       
   191   end if
       
   192 
       
   193   npts = dimsizes(ind_z)
       
   194 
       
   195 ;------------------------------------------------------
       
   196 ; for metric table computation
       
   197   amp_ob        = new((/npts/),float)
       
   198   amp_model     = new((/npts/),float)
       
   199 
       
   200   amp_ratio_sta = new((/npts/),float)
       
   201   ccr_sta       = new((/npts/),float)
       
   202   M_sta         = new((/npts/),float)
       
   203   score_sta     = new((/npts/),float)
       
   204 
       
   205   var_sta       = new((/npts/),float)
       
   206 ;-----------------------------------------------------
       
   207 ; for station line plot
       
   208   npts_str = ""
       
   209   npts_str = npts
       
   210 
       
   211   plot_data   = new((/2,12,npts/),float)
       
   212   plot_data_0 = new((/12,npts/),float)
       
   213 
       
   214   plot_data!0 = "case"
       
   215   plot_data!1 = "month"
       
   216   plot_data!2 = "pts"
       
   217   plot_data@long_name   = "CO2 Annual Cycle (ppm)"
       
   218 
       
   219   plot_data_0!0 = "month"
       
   220   plot_data_0!1 = "pts"
       
   221   plot_data_0@long_name = "CO2 (ppm)"
       
   222 ;--------------------------------------------------------------------------
       
   223   do n=0,npts-1
       
   224 
       
   225      amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:)) 
       
   226      amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
       
   227 
       
   228      amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
       
   229      ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
       
   230      M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
       
   231      score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
       
   232  
       
   233      var_model = stddev(val_model(ind_z(n),:))
       
   234      var_ob    = stddev(val_ob(ind_z(n),:))
       
   235      var_sta(n)= var_model/var_ob 
       
   236 ;----------------------------------------------------------------------
       
   237 ; for station line plot
       
   238 
       
   239      plot_data(0,:,n) = (/val_model(ind_z(n),:)/)
       
   240      plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)
       
   241 
       
   242      plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)
       
   243    
       
   244      plot_name = sta(ind_z(n))    
       
   245      title = plot_name+" ("+lat(ind_z(n))+","+lon(ind_z(n))+")"
       
   246 
       
   247      wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
       
   248 ;------------------------------------------
       
   249 ;    for panel plot
       
   250    
       
   251      plot=new(2,graphic)                        ; create graphic array
       
   252      res@gsnFrame     = False                   ; Do not draw plot 
       
   253      res@gsnDraw      = False                   ; Do not advance frame
       
   254 
       
   255      pres                            = True     ; panel plot mods desired
       
   256      pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
       
   257                                                ; indiv. plots in panel
       
   258      pres@gsnMaximize                = True     ; fill the page
       
   259 ;------------------------------------------
       
   260      res@tiMainString = title                           ; add title
       
   261 
       
   262      plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res)   ; create plot 1
       
   263 
       
   264      plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2
       
   265 
       
   266      gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot
       
   267 
       
   268      delete (wks)
       
   269      delete (plot)
       
   270 
       
   271      system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
       
   272             "rm "+plot_name+"."+plot_type)
       
   273 
       
   274 ;---------------------------------------------------------------------------  
       
   275   end do
       
   276 
       
   277 ;-------------------------------------------------------------------------
       
   278 ; for Taylor plot in a zone
       
   279 
       
   280 ; Cases [Model]
       
   281   case      = (/ model_name /) 
       
   282   nCase     = dimsizes(case )                 ; # of Cases [Cases]
       
   283 
       
   284 ; station compared
       
   285   var       = sta(ind_z) 
       
   286   nVar      = dimsizes(var)                   ; # of stations
       
   287 
       
   288 ; arrays to be passed to taylor plot 
       
   289   ratio      = new ((/nCase, nVar/),typeof(var_sta) )  
       
   290   cc         = new ((/nCase, nVar/),typeof(ccr_sta) ) 
       
   291 
       
   292   ratio(0,:) = var_sta 
       
   293   cc(0,:)    = ccr_sta
       
   294 
       
   295 ;---------------------------------
       
   296 ; create plot
       
   297 
       
   298   rest   = True                           ; default taylor diagram
       
   299         
       
   300   rest@Markers      = (/16/)               ; make all solid fill
       
   301   rest@Colors       = (/"red" /)          
       
   302 ; rest@varLabels    = var
       
   303   rest@caseLabels   = case
       
   304 
       
   305   rOpts  = True
       
   306   rOpts@caseLabels = True
       
   307 ; rOpts@caseLabelsFontHeightF = 0.05
       
   308   rOpts@caseLabelsFontHeightF = 0.15
       
   309 
       
   310   plot_name = "taylor_diagram_"+ zone
       
   311   title = "CO2 Annual Cycle:  "+ zone
       
   312   rest@tiMainString = title
       
   313   
       
   314   wks  = gsn_open_wks (plot_type,plot_name)        ; open workstation  
       
   315   plot = taylor_diagram(wks,ratio,cc,rest)
       
   316 
       
   317   delete (wks)
       
   318   delete (plot)
       
   319 
       
   320   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
       
   321          "rm "+plot_name+"."+plot_type)
       
   322 
       
   323   delete (ratio)
       
   324   delete (cc)
       
   325   delete (var_sta)
       
   326   delete (var)
       
   327 
       
   328 ;-------------------------------------------------------------------------
       
   329 ; for line plot in a zone
       
   330 
       
   331   plot_name = "All_"+npts_str
       
   332   title = plot_name + " in "+ zone
       
   333 
       
   334   wks = gsn_open_wks (plot_type,plot_name)        ; open workstation
       
   335 ;-----------------------------------------
       
   336 ; for panel plot
       
   337    
       
   338   plot=new(2,graphic)                        ; create graphic array
       
   339   res@gsnFrame     = False                   ; Do not draw plot 
       
   340   res@gsnDraw      = False                   ; Do not advance frame
       
   341 
       
   342   pres                            = True     ; panel plot mods desired
       
   343   pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
       
   344                                                ; indiv. plots in panel
       
   345   pres@gsnMaximize                = True     ; fill the page
       
   346 ;-----------------------------------------
       
   347   res@tiMainString = title                                     ; add title
       
   348 
       
   349   plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res)   ; create plot 1
       
   350     
       
   351   plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2
       
   352 
       
   353   gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot
       
   354 
       
   355   delete (wks)
       
   356   delete (plot)
       
   357 
       
   358   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
       
   359          "rm "+plot_name+"."+plot_type)
       
   360 
       
   361   delete (plot_data)
       
   362   delete (plot_data_0)    
       
   363 ;---------------------------------------------------------------------------
       
   364 ; values saved for zone table 
       
   365 
       
   366   amp_ratio_zone = avg(amp_ratio_sta)
       
   367   ccr_zone       = avg(ccr_sta)
       
   368   M_zone         = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts) 
       
   369   score_zone     = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
       
   370 
       
   371   text(z,0) = zone
       
   372   text(z,1) = sprintf("%.0f", npts)
       
   373   text(z,2) = sprintf("%.2f", amp_ratio_zone)
       
   374   text(z,3) = sprintf("%.2f", ccr_zone)
       
   375   text(z,4) = sprintf("%.2f", M_zone)
       
   376   text(z,5) = sprintf("%.2f", score_zone)
       
   377   text(z,6) = zone  
       
   378 
       
   379 ;*******************************************************************
       
   380 ; html table -- station
       
   381 ;*******************************************************************
       
   382   output_html = "score+line_"+zone+".html"
       
   383 
       
   384   header = (/"<HTML>" \
       
   385             ,"<HEAD>" \
       
   386             ,"<TITLE>CLAMP metrics</TITLE>" \
       
   387             ,"</HEAD>" \
       
   388             ,"<H1>Latitude Zone "+zone+": Model "+model_name+"</H1>" \
       
   389             /) 
       
   390   footer = "</HTML>"
       
   391 
       
   392   table_header = (/ \
       
   393         "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
       
   394        ,"<tr>" \
       
   395        ,"   <th bgcolor=DDDDDD >Site Name</th>" \
       
   396        ,"   <th bgcolor=DDDDDD >Latitude</th>" \
       
   397        ,"   <th bgcolor=DDDDDD >Longitude</th>" \
       
   398        ,"   <th bgcolor=DDDDDD >model vs obs.<br>amplitude ratio</th>" \
       
   399        ,"   <th bgcolor=DDDDDD >Correlation Coef.</th>" \
       
   400        ,"   <th bgcolor=DDDDDD >M Score</th>" \
       
   401        ,"   <th bgcolor=DDDDDD >Combined Score</th>" \
       
   402        ,"</tr>" \
       
   403        /)
       
   404   table_footer = "</table>"
       
   405   row_header = "<tr>"
       
   406   row_footer = "</tr>"
       
   407 
       
   408   lines = new(50000,string)
       
   409   nline = 0
       
   410 
       
   411   set_line(lines,nline,header)
       
   412   set_line(lines,nline,table_header)
       
   413 ;-----------------------------------------------
       
   414 ; row of table
       
   415   
       
   416   do n = 0,npts-1
       
   417      set_line(lines,nline,row_header)
       
   418 
       
   419      txt0 = sta(ind_z(n))
       
   420      txt1 = sprintf("%5.2f", (/lat(ind_z(n))/))
       
   421      txt2 = sprintf("%5.2f", (/lon(ind_z(n))/))
       
   422      txt3 = sprintf("%5.2f", (/amp_ratio_sta(n)/))
       
   423      txt4 = sprintf("%5.2f", (/ccr_sta(n)/))
       
   424      txt5 = sprintf("%5.2f", (/M_sta(n)/))
       
   425      txt6 = sprintf("%5.2f", (/score_sta(n)/))
       
   426 
       
   427      set_line(lines,nline,"<th><a href="+txt0+".png>"+txt0+"</a></th>")
       
   428      set_line(lines,nline,"<th>"+txt1+"</th>")
       
   429      set_line(lines,nline,"<th>"+txt2+"</th>")
       
   430      set_line(lines,nline,"<th>"+txt3+"</th>")
       
   431      set_line(lines,nline,"<th>"+txt4+"</th>")
       
   432      set_line(lines,nline,"<th>"+txt5+"</th>")
       
   433      set_line(lines,nline,"<th>"+txt6+"</th>")
       
   434 
       
   435      set_line(lines,nline,row_footer)
       
   436   end do
       
   437 
       
   438 ; last row, summary
       
   439   set_line(lines,nline,row_header)
       
   440 
       
   441   txt0 = "All_"+sprintf("%.0f", (/npts/))
       
   442   txt1 = "-"
       
   443   txt2 = "-"
       
   444   txt3 = sprintf("%5.2f", (/amp_ratio_zone/))
       
   445   txt4 = sprintf("%5.2f", (/ccr_zone/))
       
   446   txt5 = sprintf("%5.2f", (/M_zone/))
       
   447   txt6 = sprintf("%5.2f", (/score_zone/))
       
   448 
       
   449   set_line(lines,nline,"<th><a href="+txt0+".png>"+txt0+"</a></th>")
       
   450   set_line(lines,nline,"<th>"+txt1+"</th>")
       
   451   set_line(lines,nline,"<th>"+txt2+"</th>")
       
   452   set_line(lines,nline,"<th>"+txt3+"</th>")
       
   453   set_line(lines,nline,"<th>"+txt4+"</th>")
       
   454   set_line(lines,nline,"<th>"+txt5+"</th>")
       
   455   set_line(lines,nline,"<th>"+txt6+"</th>")
       
   456 
       
   457   set_line(lines,nline,row_footer)
       
   458 ;-----------------------------------------------
       
   459   set_line(lines,nline,table_footer)
       
   460   set_line(lines,nline,footer) 
       
   461 
       
   462 ; Now write to an HTML file.
       
   463   idx = ind(.not.ismissing(lines))
       
   464   if(.not.any(ismissing(idx))) then
       
   465     asciiwrite(output_html,lines(idx))
       
   466   else
       
   467    print ("error?")
       
   468   end if
       
   469   delete (idx)
       
   470 ;-----------------------------------------------------------------
       
   471 
       
   472   delete (ind_z)
       
   473   delete (amp_model)
       
   474   delete (amp_ob)
       
   475   delete (amp_ratio_sta)
       
   476   delete (ccr_sta)
       
   477   delete (M_sta)
       
   478   delete (score_sta)
       
   479 end do
       
   480 
       
   481 ;*******************************************************************
       
   482 ; html table -- zone
       
   483 ;*******************************************************************
       
   484   output_html = "score+line_vs_ob.html"
       
   485 
       
   486   header = (/"<HTML>" \
       
   487             ,"<HEAD>" \
       
   488             ,"<TITLE>CLAMP metrics</TITLE>" \
       
   489             ,"</HEAD>" \
       
   490             ,"<H1>CO2 Seasonal Cycle Comparisons by Latitude Zone: Model "+model_name+"</H1>" \
       
   491             /) 
       
   492   footer = "</HTML>"
       
   493 
       
   494   delete (table_header)
       
   495   table_header = (/ \
       
   496         "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
       
   497        ,"<tr>" \
       
   498        ,"   <th bgcolor=DDDDDD >Zone</th>" \
       
   499        ,"   <th bgcolor=DDDDDD >Number of Site</th>" \
       
   500        ,"   <th bgcolor=DDDDDD >model vs obs.<br>amplitide ratio</th>" \
       
   501        ,"   <th bgcolor=DDDDDD >Correlation Coef.</th>" \
       
   502        ,"   <th bgcolor=DDDDDD >M Score</th>" \
       
   503        ,"   <th bgcolor=DDDDDD >Combined Score</th>" \
       
   504        ,"   <th bgcolor=DDDDDD >Taylor diagram</th>" \
       
   505        ,"</tr>" \
       
   506        /)
       
   507   table_footer = "</table>"
       
   508   row_header = "<tr>"
       
   509   row_footer = "</tr>"
       
   510 
       
   511   lines = new(50000,string)
       
   512   nline = 0
       
   513 
       
   514   set_line(lines,nline,header)
       
   515   set_line(lines,nline,table_header)
       
   516 ;-----------------------------------------------
       
   517 ;row of table
       
   518  
       
   519   do n = 0,nrow_zone-1
       
   520      set_line(lines,nline,row_header)
       
   521 
       
   522      set_line(lines,nline,"<th><a href=score+line_"+text(n,0)+".html>"+text(n,0)+"</th>")
       
   523      set_line(lines,nline,"<th>"+text(n,1)+"</th>")
       
   524      set_line(lines,nline,"<th>"+text(n,2)+"</th>")
       
   525      set_line(lines,nline,"<th>"+text(n,3)+"</th>")
       
   526      set_line(lines,nline,"<th>"+text(n,4)+"</th>")
       
   527      set_line(lines,nline,"<th>"+text(n,5)+"</th>")
       
   528      set_line(lines,nline,"<th><a href=taylor_diagram_"+text(n,6)+".png>Taylor_diagram</th>")
       
   529 
       
   530      set_line(lines,nline,row_footer)
       
   531   end do
       
   532 
       
   533 ; for the last row
       
   534 
       
   535      txt0 = "All"
       
   536      txt1 = sum(stringtofloat(text(0:3,1))) 
       
   537      txt2 = "-"
       
   538      txt3 = "-"
       
   539      txt4 = "-"
       
   540      txt5 = sum(stringtofloat(text(0:3,5)))
       
   541      txt6 = "-"
       
   542 
       
   543      set_line(lines,nline,row_header)
       
   544 
       
   545      set_line(lines,nline,"<th>"+txt0+"</th>")
       
   546      set_line(lines,nline,"<th>"+txt1+"</th>")
       
   547      set_line(lines,nline,"<th>"+txt2+"</th>")
       
   548      set_line(lines,nline,"<th>"+txt3+"</th>")
       
   549      set_line(lines,nline,"<th>"+txt4+"</th>")
       
   550      set_line(lines,nline,"<th>"+txt5+"</th>")
       
   551      set_line(lines,nline,"<th>"+txt6+"</th>")
       
   552 
       
   553      set_line(lines,nline,row_footer)
       
   554 ;-----------------------------------------------
       
   555   set_line(lines,nline,table_footer)
       
   556   set_line(lines,nline,footer) 
       
   557 
       
   558 ; Now write to an HTML file.
       
   559   idx = ind(.not.ismissing(lines))
       
   560   if(.not.any(ismissing(idx))) then
       
   561     asciiwrite(output_html,lines(idx))
       
   562   else
       
   563    print ("error?")
       
   564   end if
       
   565   delete (idx)
       
   566 ;--------------------------------------------------------------------------
       
   567  
       
   568   M_co2 = txt5
       
   569 
       
   570   if (isvar("compare")) then
       
   571 
       
   572      system("sed -e '1,/M_co2/s/M_co2/"+M_co2+"/' "+html_name2+" > "+html_new2+";"+ \
       
   573             "mv -f "+html_new2+" "+html_name2)
       
   574   end if
       
   575 
       
   576   system("sed s#M_co2#"+M_co2+"# "+html_name+" > "+html_new+";"+ \
       
   577          "mv -f "+html_new+" "+html_name)
       
   578 
       
   579 ;***************************************************************************
       
   580 ; add total score and write to file
       
   581 ;***************************************************************************
       
   582   M_total = M_co2
       
   583 
       
   584   asciiwrite("M_save.co2", M_total)
       
   585  
       
   586 ;***************************************************************************
       
   587 ; output plots
       
   588 ;***************************************************************************
       
   589   output_dir = model_name+"/co2"
       
   590 
       
   591   system("mv *.png *.html " + output_dir)
       
   592 ;*************************************************************************** 
       
   593 end