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