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