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