co2/43.metric+lines.ncl
changeset 0 0c6405ab2ff4
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/co2/43.metric+lines.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,553 @@
     1.4 +; ***********************************************
     1.5 +; using gsn_table for all
     1.6 +; combine 19.metric_plot.ncl and 24.lines.ncl
     1.7 +; ***********************************************
     1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
     1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
    1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    1.11 +;load "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"
    1.12 +;************************************************
    1.13 +begin
    1.14 +;************************************************
    1.15 +; read in data: observed
    1.16 +;************************************************
    1.17 + diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
    1.18 + fili  = "co2_globalView_98.nc"
    1.19 + g     = addfile (diri+fili,"r")
    1.20 + val   = g->CO2_SEAS  
    1.21 + lon   = g->LON 
    1.22 + lat   = g->LAT
    1.23 + sta   = chartostring(g->STATION) 
    1.24 + delete (g)
    1.25 + 
    1.26 +;print (sta(0))
    1.27 +
    1.28 + ncase = dimsizes(lat)
    1.29 +;print (ncase)
    1.30 +
    1.31 +;**************************************************************
    1.32 +; get only the lowest level at each station 
    1.33 +;**************************************************************
    1.34 + lat_tmp = lat
    1.35 + lat_tmp@_FillValue = 1.e+36
    1.36 + 
    1.37 + do n = 0,ncase-1
    1.38 +    if (.not. ismissing(lat_tmp(n))) then 
    1.39 +       indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
    1.40 +       if (dimsizes(indexes) .gt. 1) then
    1.41 +          lat_tmp(indexes(1:)) = lat_tmp@_FillValue
    1.42 +       end if
    1.43 +       delete (indexes)
    1.44 +    end if
    1.45 + end do
    1.46 +
    1.47 + indexes = ind(.not. ismissing(lat_tmp))
    1.48 +;print (dimsizes(indexes))
    1.49 +;print (indexes)
    1.50 + 
    1.51 + lat_ob = lat(indexes)
    1.52 + lon_ob = lon(indexes)
    1.53 + val_ob = val(indexes,:)
    1.54 +;printVarSummary (val_ob)
    1.55 +;print (lat_ob +"/"+lon_ob)
    1.56 +
    1.57 +;************************************************
    1.58 +; read in model data
    1.59 +;************************************************
    1.60 +  diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
    1.61 +; fili2 = "b30.061m_401_425_MONS_climo_atm.nc"
    1.62 +  fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc"
    1.63 +
    1.64 +  g     = addfile(diri2+fili2,"r")
    1.65 +  x     = g->CO2
    1.66 +  xi    = g->lon
    1.67 +  yi    = g->lat
    1.68 +  xdim  = dimsizes(x)
    1.69 +  nlev  = xdim(1)
    1.70 +  y     = x(:,0,:,:)
    1.71 +; printVarSummary (y)
    1.72 +  
    1.73 +; get the co2 at the lowest level
    1.74 +  y     = x(:,nlev-1,:,:)
    1.75 +
    1.76 +; change to unit of observed (u mol/mol)
    1.77 +; Model_units [=] kgCO2 / kgDryAir
    1.78 +; 28.966 = molecular weight of dry air
    1.79 +; 44.       = molecular weight of CO2
    1.80 +; u mol = 1e-6 mol
    1.81 +
    1.82 +  factor = (28.966/44.) * 1e6
    1.83 +  y      = y * factor
    1.84 +
    1.85 +  y@_FillValue = 1.e36
    1.86 +  y@units      = "u mol/mol"
    1.87 +; y = where(y0 .lt. 287.,y@_FillValue,y)
    1.88 +; printVarSummary (y)
    1.89 +; print (min(y)+"/"+max(y))
    1.90 +
    1.91 +; interpolate into observed station
    1.92 +; note: model is 0-360E, 90S-90N
    1.93 +
    1.94 +; to be able to handle observation at (-89.98,-24.80)
    1.95 +  print (yi(0))
    1.96 +  yi(0) = -90.
    1.97 +
    1.98 +  i = ind(lon_ob .lt. 0.)
    1.99 +  lon_ob(i) = lon_ob(i) + 360.  
   1.100 +
   1.101 +  yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
   1.102 +
   1.103 +  val_model = yo(pts|:,time|:)
   1.104 +  val_model_0 = val_model
   1.105 +; printVarSummary (val_model)
   1.106 +; print (min(val_model)+"/"+max(val_model))
   1.107 +
   1.108 +; remove annual mean
   1.109 +  val_model = val_model - conform(val_model,dim_avg(val_model),0)
   1.110 +; print (min(val_model)+"/"+max(val_model))
   1.111 +
   1.112 +  nzone = 4
   1.113 +
   1.114 +;*******************************************************************
   1.115 +; for table -- zone
   1.116 +;*******************************************************************
   1.117 +; table header name
   1.118 +  table_header_name = "Zone" 
   1.119 +
   1.120 +; column (not including header column)
   1.121 +  col_header_zone = (/"Stations","Amplitude Ratio", \
   1.122 +                      "Correlation Coef","M score","Combined Score"/)
   1.123 +  ncol_zone       = dimsizes(col_header_zone ) 
   1.124 +
   1.125 +; row (not including header row)
   1.126 +  row_header_zone = (/"60N-90N","30N-60N","EQ-30N","90S-EQ","Total"/)  
   1.127 +  nrow_zone       = dimsizes(row_header_zone)                  
   1.128 +
   1.129 +; arrays to be passed to table. 
   1.130 +  value_zone = new ((/nrow_zone, ncol_zone/),string ) 
   1.131 +;--------------------------------------------------------------------
   1.132 +
   1.133 +  table_length_zone = 0.4 
   1.134 +
   1.135 +; Table header
   1.136 +  ncr1  = (/1,1/)               ; 1 row, 1 column
   1.137 +  x1    = (/0.005,0.15/)        ; Start and end X
   1.138 +  y1    = (/0.900,0.995/)       ; Start and end Y
   1.139 +  text1 = table_header_name
   1.140 +  res1               = True
   1.141 +  res1@txFontHeightF = 0.03
   1.142 +  res1@gsFillColor   = "CornFlowerBlue"
   1.143 +
   1.144 +; Column header (equally space in x2)
   1.145 +  ncr2  = (/1,ncol_zone/)         ; 1 rows, 5 columns
   1.146 +  x2    = (/x1(1),0.995/)          ; start from end of x1
   1.147 +  y2    = y1                      ; same as y1
   1.148 +  text2 = col_header_zone
   1.149 +  res2               = True
   1.150 +  res2@txFontHeightF = 0.015
   1.151 +  res2@gsFillColor   = "Gray"
   1.152 +
   1.153 +; Row header (equally space in y2)
   1.154 +  ncr3  = (/nrow_zone,1/)         ; 5 rows, 1 columns
   1.155 +  x3    = x1                      ; same as x1
   1.156 +  y3    = (/1.0-table_length_zone,0.900/) ; end at start of y1
   1.157 +  text3 = row_header_zone
   1.158 +  res3               = True
   1.159 +  res3@txFontHeightF = 0.02
   1.160 +  res3@gsFillColor   = "Gray"
   1.161 +
   1.162 +; Main table body
   1.163 +  ncr4  = (/nrow_zone,ncol_zone/) ; 5 rows, 5 columns
   1.164 +  x4    = x2                      ; Start and end x
   1.165 +  y4    = y3                      ; Start and end Y
   1.166 +  text4 = new((/nrow_zone,ncol_zone/),string)
   1.167 +
   1.168 +  color_fill4      = new((/nrow_zone,ncol_zone/),string)
   1.169 +  color_fill4      = "white"
   1.170 +  color_fill4(:,ncol_zone-1) = "grey"
   1.171 +
   1.172 +  res4               = True       ; Set up resource list
   1.173 +; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.
   1.174 +  res4@txFontHeightF = 0.02
   1.175 +  res4@gsFillColor   = color_fill4
   1.176 +
   1.177 +  delete (color_fill4)
   1.178 +
   1.179 +;*******************************************************************
   1.180 +; for table -- station
   1.181 +;*******************************************************************
   1.182 +; column (not including header column)
   1.183 +  col_header_sta = (/"Latitude","Longitude","Amplitude Ratio", \
   1.184 +                      "Correlation Coef","M score","Combined Score"/)
   1.185 +  ncol_sta       = dimsizes(col_header_sta ) 
   1.186 +;-------------------------------------------------------------------
   1.187 +
   1.188 +; Table header
   1.189 +  ncr5  = (/1,1/)               ; 1 row, 1 column
   1.190 +  x5    = (/0.005,0.15/)        ; Start and end X
   1.191 +  y5    = (/0.900,0.995/)       ; Start and end Y
   1.192 +  res5               = True
   1.193 +  res5@txFontHeightF = 0.02
   1.194 +  res5@gsFillColor   = "CornFlowerBlue"
   1.195 +
   1.196 +; Column header (equally space in x2)
   1.197 +  ncr6  = (/1,ncol_sta/)         ; 1 rows, 5 columns
   1.198 +  x6    = (/x5(1),0.995/)          ; start from end of x1
   1.199 +  y6    = y5                      ; same as y1
   1.200 +  text6 = col_header_sta
   1.201 +  res6               = True
   1.202 +  res6@txFontHeightF = 0.012
   1.203 +  res6@gsFillColor   = "Gray"
   1.204 +
   1.205 +; Row header (equally space in y2)
   1.206 +
   1.207 +  res7               = True
   1.208 +  res7@txFontHeightF = 0.015
   1.209 +  res7@gsFillColor   = "Gray"
   1.210 +;--------------------------------------------------------------
   1.211 +; for station line plot
   1.212 +
   1.213 +; for x-axis in xyplot
   1.214 +  mon = ispan(1,12,1)
   1.215 +  mon@long_name = "month"
   1.216 +
   1.217 +  plot_type = "ps"
   1.218 +  plot_type_new = "png"
   1.219 +
   1.220 +  res                   = True                      ; plot mods desired
   1.221 +  res@xyLineThicknesses = (/2.0,2.0,2.0/)           ; make 2nd lines thicker
   1.222 +  res@xyLineColors      = (/"red","black"/)  ; change line color
   1.223 +
   1.224 +; Add a boxed legend using the more simple method
   1.225 +
   1.226 +  res@pmLegendDisplayMode    = "Always"
   1.227 +; res@pmLegendWidthF         = 0.1
   1.228 +  res@pmLegendWidthF         = 0.08
   1.229 +  res@pmLegendHeightF        = 0.06
   1.230 +; res@pmLegendOrthogonalPosF = -1.17
   1.231 +; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
   1.232 +  res@pmLegendOrthogonalPosF = -0.30  ;(downward)
   1.233 +
   1.234 +; res@pmLegendParallelPosF   =  0.18
   1.235 +  res@pmLegendParallelPosF   =  0.23  ;(rightward)
   1.236 +
   1.237 +; res@lgPerimOn             = False
   1.238 +  res@lgLabelFontHeightF     = 0.015
   1.239 +  res@xyExplicitLegendLabels = (/"b30.061n","observed"/)
   1.240 +;-------------------------------------------------------------------
   1.241 +
   1.242 +  do z = 0,nzone-1
   1.243 +
   1.244 +  if (z .eq. 0) then 
   1.245 +;    maximum score for the zone, 60N-90N
   1.246 +     zone = "60N-90N" 
   1.247 +     score_max = 5.0
   1.248 +;    index of stations in this zone
   1.249 +     ind_z = ind(lat_ob .ge. 60.)
   1.250 +;    print (ind_z)
   1.251 +;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
   1.252 +;    print (val_ob(ind_z,:))
   1.253 +;    print (val_model(ind_z,:))
   1.254 +  end if
   1.255 +
   1.256 +  if (z .eq. 1) then 
   1.257 +;    maximum score for the zone, 30N-60N
   1.258 +     zone = "30N-60N" 
   1.259 +     score_max = 5.0
   1.260 +;    index of stations in this zone
   1.261 +     ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
   1.262 +;    print (ind_z)
   1.263 +;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
   1.264 +;    print (val_ob(ind_z,:))
   1.265 +;    print (val_model(ind_z,:))
   1.266 +  end if
   1.267 +
   1.268 +  if (z .eq. 2) then 
   1.269 +;    maximum score for the zone, EQ-30N 
   1.270 +     zone = "EQ-30N"
   1.271 +     score_max = 5.0
   1.272 +;    index of stations in this zone
   1.273 +     ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
   1.274 +;    print (ind_z)
   1.275 +;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
   1.276 +;    print (val_ob(ind_z,:))
   1.277 +;    print (val_model(ind_z,:))
   1.278 +  end if
   1.279 +
   1.280 +  if (z .eq. 3) then 
   1.281 +;    maximum score for the zone, 90S-EQ
   1.282 +     zone = "90S-EQ" 
   1.283 +     score_max = 5.0
   1.284 +;    index of stations in this zone
   1.285 +     ind_z = ind(lat_ob .lt. 0. )
   1.286 +;    print (ind_z)
   1.287 +;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
   1.288 +;    print (val_ob(ind_z,:))
   1.289 +;    print (val_model(ind_z,:))
   1.290 +  end if
   1.291 +
   1.292 + npts = dimsizes(ind_z)
   1.293 + print (npts)
   1.294 +;-------------------------------------------------------------------------
   1.295 +; for metric table computation
   1.296 +
   1.297 + amp_ob        = new((/npts/),float)
   1.298 + amp_model     = new((/npts/),float)
   1.299 +
   1.300 + amp_ratio_sta = new((/npts/),float)
   1.301 + ccr_sta       = new((/npts/),float)
   1.302 + M_sta         = new((/npts/),float)
   1.303 + score_sta     = new((/npts/),float)
   1.304 +
   1.305 +;----------------------
   1.306 +; for table -- station
   1.307 +;----------------------
   1.308 +
   1.309 +; row (not including header row)  
   1.310 +  nrow_sta       = npts                  
   1.311 +
   1.312 +; Table header
   1.313 +  ncr5  = (/1,1/)               ; 1 row, 1 column
   1.314 +  x5    = (/0.005,0.15/)        ; Start and end X
   1.315 +  y5    = (/0.900,0.995/)       ; Start and end Y
   1.316 +  res5               = True
   1.317 +  res5@txFontHeightF = 0.02
   1.318 +  res5@gsFillColor   = "CornFlowerBlue"
   1.319 +
   1.320 +; Column header (equally space in x2)
   1.321 +  ncr6  = (/1,ncol_sta/)         ; 1 rows, 5 columns
   1.322 +  x6    = (/x1(1),0.995/)          ; start from end of x1
   1.323 +  y6    = y1                      ; same as y1
   1.324 +  text6 = col_header_sta
   1.325 +  res6               = True
   1.326 +  res6@txFontHeightF = 0.012
   1.327 +  res6@gsFillColor   = "Gray"
   1.328 +
   1.329 +; Row header (equally space in y2)
   1.330 +  ncr7  = (/nrow_sta,1/)         ; 5 rows, 1 columns
   1.331 +  x7    = x1                      ; same as x1
   1.332 +; y7    = (/1.0-table_length_sta,0.900/) ; end at start of y1
   1.333 +  text7 = new((/nrow_sta/),string)
   1.334 +  res7               = True
   1.335 +  res7@txFontHeightF = 0.015
   1.336 +  res7@gsFillColor   = "Gray"
   1.337 +;-------------------------------------------------------------------------
   1.338 +; for station line plot
   1.339 +
   1.340 +  npts_str = ""
   1.341 +  npts_str = npts
   1.342 +; print (npts_str)
   1.343 +
   1.344 +  plot_data   = new((/2,12,npts/),float)
   1.345 +  plot_data_0 = new((/12,npts/),float)
   1.346 +
   1.347 +  plot_data!0 = "case"
   1.348 +  plot_data!1 = "month"
   1.349 +  plot_data!2 = "pts"
   1.350 +  plot_data@long_name   = "CO2 Seasonal"
   1.351 +
   1.352 +  plot_data_0!0 = "month"
   1.353 +  plot_data_0!1 = "pts"
   1.354 +  plot_data_0@long_name = "CO2"
   1.355 +;--------------------------------------------------------------------------
   1.356 + do n=0,npts-1
   1.357 +    amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:)) 
   1.358 +    amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
   1.359 +
   1.360 +    amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
   1.361 +    ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
   1.362 +    M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
   1.363 +    score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
   1.364 +
   1.365 +    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))
   1.366 +;----------------------------------------------------------------------
   1.367 +; for station line plot
   1.368 +
   1.369 +    plot_data(0,:,n) = (/val_model(ind_z(n),:)/)
   1.370 +    plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)
   1.371 +
   1.372 +    plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)
   1.373 +   
   1.374 +    plot_name = sta(ind_z(n))    
   1.375 +    title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"
   1.376 +;   print (title)
   1.377 +;   print (plot_name)
   1.378 +
   1.379 +    wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
   1.380 +;------------------------------------------
   1.381 +;   for panel plot
   1.382 +   
   1.383 +    plot=new(2,graphic)                        ; create graphic array
   1.384 +    res@gsnFrame     = False                   ; Do not draw plot 
   1.385 +    res@gsnDraw      = False                   ; Do not advance frame
   1.386 +
   1.387 +    pres                            = True     ; panel plot mods desired
   1.388 +    pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
   1.389 +                                               ; indiv. plots in panel
   1.390 +    pres@gsnMaximize                = True     ; fill the page
   1.391 +;------------------------------------------
   1.392 +    res@tiMainString = title                           ; add title
   1.393 +
   1.394 +    plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res)   ; create plot 1
   1.395 +
   1.396 +    plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2
   1.397 +
   1.398 +    gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot
   1.399 +
   1.400 +    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   1.401 +    system("rm "+plot_name+"."+plot_type)
   1.402 +    clear (wks)
   1.403 +;---------------------------------------------------------------------------  
   1.404 + end do
   1.405 +;-------------------------------------------------------------------------
   1.406 +;   for line plot in a zone
   1.407 +
   1.408 +    plot_name = "All_"+npts_str
   1.409 +    title = plot_name + " in "+ zone
   1.410 +;   print (title)
   1.411 +;   print (plot_name)
   1.412 +
   1.413 +    wks = gsn_open_wks (plot_type,plot_name)        ; open workstation
   1.414 +;-----------------------------------------
   1.415 +;   for panel plot
   1.416 +   
   1.417 +    plot=new(2,graphic)                        ; create graphic array
   1.418 +    res@gsnFrame     = False                   ; Do not draw plot 
   1.419 +    res@gsnDraw      = False                   ; Do not advance frame
   1.420 +
   1.421 +    pres                            = True     ; panel plot mods desired
   1.422 +    pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
   1.423 +                                               ; indiv. plots in panel
   1.424 +    pres@gsnMaximize                = True     ; fill the page
   1.425 +;-----------------------------------------
   1.426 +    res@tiMainString = title                                     ; add title
   1.427 +
   1.428 +    plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res)   ; create plot 1
   1.429 +    
   1.430 +    plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2
   1.431 +
   1.432 +    gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot
   1.433 +
   1.434 +    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   1.435 +    system("rm "+plot_name+"."+plot_type)
   1.436 +
   1.437 +    clear (wks)
   1.438 +;   delete (ind_z)
   1.439 +    delete (plot_data)
   1.440 +    delete (plot_data_0)    
   1.441 +;---------------------------------------------------------------------------
   1.442 +; for zone table value
   1.443 +
   1.444 + amp_ratio_zone = avg(amp_ratio_sta)
   1.445 + ccr_zone       = avg(ccr_sta)
   1.446 + M_zone   = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts) 
   1.447 + score_zone   = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
   1.448 +
   1.449 + print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)
   1.450 +
   1.451 +  text4(z,0) = sprintf("%5.2f", npts)
   1.452 +  text4(z,1) = sprintf("%5.2f", amp_ratio_zone)
   1.453 +  text4(z,2) = sprintf("%5.2f", ccr_zone)
   1.454 +  text4(z,3) = sprintf("%5.2f", M_zone)
   1.455 +  text4(z,4) = sprintf("%5.2f", score_zone)  
   1.456 +;---------------------------------------------------------------------------
   1.457 +; plot station table
   1.458 +;----------------------------------
   1.459 +; header value for station table
   1.460 +  text5  = zone 
   1.461 +
   1.462 +; row value for station table 
   1.463 +  text7  = sta(ind_z)
   1.464 +
   1.465 +  if (z .eq. 0) then
   1.466 +     table_length_sta = 0.5
   1.467 +  end if
   1.468 +  if (z .eq. 1) then
   1.469 +     table_length_sta = 0.995
   1.470 +  end if
   1.471 +  if (z .eq. 2) then
   1.472 +     table_length_sta = 0.8
   1.473 +  end if
   1.474 +  if (z .eq. 3) then
   1.475 +     table_length_sta = 0.8
   1.476 +  end if    
   1.477 +
   1.478 +  x7 = x5   
   1.479 +  y7 = (/1.0-table_length_sta,0.900/) ; end at start of y1
   1.480 +
   1.481 +; Main table body
   1.482 +  ncr8  = (/nrow_sta,ncol_sta/) ; 5 rows, 5 columns
   1.483 +  x8    = x6                      ; Start and end x
   1.484 +  y8    = y7                      ; Start and end Y
   1.485 +  text8 = new((/nrow_sta,ncol_sta/),string)
   1.486 +
   1.487 +  color_fill8      = text8
   1.488 +  color_fill8      = "white"
   1.489 +  color_fill8(:,ncol_sta-1) = "grey"
   1.490 +
   1.491 +  res8               = True       ; Set up resource list
   1.492 +  res8@gsnDebug      = True       ; Useful to print NDC row,col values used.
   1.493 +  res8@txFontHeightF = 0.02
   1.494 +  res8@gsFillColor   = color_fill8
   1.495 +
   1.496 +  delete (color_fill8)
   1.497 +
   1.498 +; table value for station table
   1.499 +  text8(:,0) = sprintf("%5.2f", (/lat(ind_z)/))
   1.500 +  text8(:,1) = sprintf("%5.2f", (/lon(ind_z)/))
   1.501 +  text8(:,2) = sprintf("%5.2f", (/amp_ratio_sta/))
   1.502 +  text8(:,3) = sprintf("%5.2f", (/ccr_sta/))
   1.503 +  text8(:,4) = sprintf("%5.2f", (/M_sta/))
   1.504 +  text8(:,5) = sprintf("%5.2f", (/score_sta/))
   1.505 + 
   1.506 +  plot_name = "table_sta." + zone 
   1.507 +  
   1.508 +  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
   1.509 +
   1.510 +  gsn_table(wks,ncr5,x5,y5,text5,res5)
   1.511 +  gsn_table(wks,ncr6,x6,y6,text6,res6)
   1.512 +  gsn_table(wks,ncr7,x7,y7,text7,res7)
   1.513 +  gsn_table(wks,ncr8,x8,y8,text8,res8) 
   1.514 +
   1.515 +  frame(wks)
   1.516 +
   1.517 +  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   1.518 +  system("rm "+plot_name+"."+plot_type)
   1.519 +
   1.520 + delete (ind_z)
   1.521 + delete (amp_model)
   1.522 + delete (amp_ob)
   1.523 + delete (amp_ratio_sta)
   1.524 + delete (ccr_sta)
   1.525 + delete (M_sta)
   1.526 + delete (score_sta)
   1.527 + delete (text7)
   1.528 + delete (text8)
   1.529 + delete (res7)
   1.530 + delete (res8)
   1.531 + clear (wks)
   1.532 + end do
   1.533 +;**************************************************
   1.534 +; plot zone table
   1.535 +;**************************************************
   1.536 + 
   1.537 +  text4(4,0) = sum(stringtofloat(text4(0:3,0))) 
   1.538 +  text4(4,1) = 0.
   1.539 +  text4(4,2) = 0.
   1.540 +  text4(4,3) = 0.
   1.541 +  text4(4,4) = sum(stringtofloat(text4(0:3,4)))
   1.542 +
   1.543 +  plot_name = "table_zone" 
   1.544 +  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
   1.545 +
   1.546 +  gsn_table(wks,ncr1,x1,y1,text1,res1)
   1.547 +  gsn_table(wks,ncr2,x2,y2,text2,res2)
   1.548 +  gsn_table(wks,ncr3,x3,y3,text3,res3)
   1.549 +  gsn_table(wks,ncr4,x4,y4,text4,res4) 
   1.550 +
   1.551 +  frame(wks)
   1.552 +
   1.553 +  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
   1.554 +  system("rm "+plot_name+"."+plot_type)
   1.555 +
   1.556 +end