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