forrest@0: ; ***********************************************
forrest@0: ; using gsn_table for all
forrest@0: ; combine 19.metric_plot.ncl and 24.lines.ncl
forrest@0: ; ***********************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ;load "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"
forrest@0: ;************************************************
forrest@0: begin
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; read model data
forrest@0: ;************************************************
forrest@0: 
forrest@0: ; from command line inputs
forrest@0: 
forrest@0: ;--------------------------------------------------
forrest@0: ; edit table.html of current model for movel1_vs_model2
forrest@0: 
forrest@0:  if (isvar("compare")) then
forrest@0:     html_name2 = compare+"/table.html"  
forrest@0:     html_new2  = html_name2 +".new"
forrest@0: 
forrest@0: ; the following only needs to execute once
forrest@0: ;   system("sed 1,/model_nameA/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0: ;          "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: ;          "sed 1,/model_nameB/s//"+model_name+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0: ;          "mv -f "+html_new2+" "+html_name2+";"+ \
forrest@0: ;          "sed s#"+modeln+"#"+model_name+"# "+html_name2+" > "+html_new2+";"+ \
forrest@0: ;          "mv -f "+html_new2+" "+html_name2)
forrest@0:  end if
forrest@0: 
forrest@0: ;-------------------------------------
forrest@0: ; edit table.html for current model
forrest@0: 
forrest@0:  html_name = model_name+"/table.html"  
forrest@0:  html_new  = html_name +".new"
forrest@0: 
forrest@0: ; the following only needs to execute once
forrest@0: ; system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
forrest@0: ;        "mv -f "+html_new+" "+html_name)
forrest@0: ;------------------------------------------------
forrest@0: 
forrest@0:   fm    = addfile(dirm+film3,"r")
forrest@0: 
forrest@0:   x     = fm->CO2
forrest@0:   xi    = fm->lon
forrest@0:   yi    = fm->lat
forrest@0: 
forrest@0:   xdim  = dimsizes(x)
forrest@0:   nlev  = xdim(1)
forrest@0:   y     = x(:,0,:,:)
forrest@0:   
forrest@0: ; get co2 at the lowest level
forrest@0:   y     = x(:,nlev-1,:,:)
forrest@0: 
forrest@0: ; change to unit of observed (u mol/mol)
forrest@0: ; Model_units [=] kgCO2 / kgDryAir
forrest@0: ; 28.966 = molecular weight of dry air
forrest@0: ; 44.       = molecular weight of CO2
forrest@0: ; u mol = 1e-6 mol
forrest@0: 
forrest@0:   factor = (28.966/44.) * 1e6
forrest@0:   y      = y * factor
forrest@0: 
forrest@0:   y@_FillValue = 1.e36
forrest@0:   y@units      = "u mol/mol"
forrest@0: ;************************************************
forrest@0: ; read data: observed
forrest@0: ;************************************************
forrest@0:  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/ob/"
forrest@0:  fili  = "co2_globalView_98.nc"
forrest@0:  g     = addfile (diri+fili,"r")
forrest@0:  val   = g->CO2_SEAS  
forrest@0:  lon   = g->LON 
forrest@0:  lat   = g->LAT
forrest@0:  sta   = chartostring(g->STATION) 
forrest@0:  delete (g)
forrest@0: 
forrest@0:  ncase = dimsizes(lat)
forrest@0: ;**************************************************************
forrest@0: ; get only the lowest level at each station 
forrest@0: ;**************************************************************
forrest@0:  lat_tmp = lat
forrest@0:  lat_tmp@_FillValue = 1.e+36
forrest@0:  
forrest@0:  do n = 0,ncase-1
forrest@0:     if (.not. ismissing(lat_tmp(n))) then 
forrest@0:        indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
forrest@0:        if (dimsizes(indexes) .gt. 1) then
forrest@0:           lat_tmp(indexes(1:)) = lat_tmp@_FillValue
forrest@0:        end if
forrest@0:        delete (indexes)
forrest@0:     end if
forrest@0:  end do
forrest@0: 
forrest@0:  indexes = ind(.not. ismissing(lat_tmp))
forrest@0:  
forrest@0:  lat_ob = lat(indexes)
forrest@0:  lon_ob = lon(indexes)
forrest@0:  val_ob = val(indexes,:)
forrest@0: ;************************************************************
forrest@0: ; interpolate model data into observed station
forrest@0: ; note: model is 0-360E, 90S-90N
forrest@0: ;************************************************************
forrest@0: ; to be able to handle observation at (-89.98,-24.80)
forrest@0:   yi(0) = -90.
forrest@0: 
forrest@0:   i = ind(lon_ob .lt. 0.)
forrest@0:   lon_ob(i) = lon_ob(i) + 360.  
forrest@0: 
forrest@0:   yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
forrest@0: 
forrest@0:   val_model = yo(pts|:,time|:)
forrest@0:   val_model_0 = val_model
forrest@0: 
forrest@0: ; remove annual mean
forrest@0:   val_model = val_model - conform(val_model,dim_avg(val_model),0)
forrest@0: 
forrest@0:   nzone = 4
forrest@0: 
forrest@0: ;*******************************************************************
forrest@0: ; for table -- zone
forrest@0: ;*******************************************************************
forrest@0: ; table header name
forrest@0:   table_header_name = "Zone" 
forrest@0: 
forrest@0: ; column (not including header column)
forrest@0:   col_header_zone = (/"Stations","Amplitude Ratio", \
forrest@0:                       "Correlation Coef","M score","Combined Score"/)
forrest@0:   ncol_zone       = dimsizes(col_header_zone ) 
forrest@0: 
forrest@0: ; row (not including header row)
forrest@0:   row_header_zone = (/"60N-90N","30N-60N","EQ-30N","90S-EQ","Total"/)  
forrest@0:   nrow_zone       = dimsizes(row_header_zone)                  
forrest@0: 
forrest@0: ; arrays to be passed to table. 
forrest@0:   value_zone = new ((/nrow_zone, ncol_zone/),string ) 
forrest@0: ;--------------------------------------------------------------------
forrest@0: 
forrest@0:   table_length_zone = 0.4 
forrest@0: 
forrest@0: ; Table header
forrest@0:   ncr1  = (/1,1/)               ; 1 row, 1 column
forrest@0:   x1    = (/0.005,0.15/)        ; Start and end X
forrest@0:   y1    = (/0.900,0.995/)       ; Start and end Y
forrest@0:   text1 = table_header_name
forrest@0:   res1               = True
forrest@0:   res1@txFontHeightF = 0.03
forrest@0:   res1@gsFillColor   = "CornFlowerBlue"
forrest@0: 
forrest@0: ; Column header (equally space in x2)
forrest@0:   ncr2  = (/1,ncol_zone/)         ; 1 rows, 5 columns
forrest@0:   x2    = (/x1(1),0.995/)          ; start from end of x1
forrest@0:   y2    = y1                      ; same as y1
forrest@0:   text2 = col_header_zone
forrest@0:   res2               = True
forrest@0:   res2@txFontHeightF = 0.015
forrest@0:   res2@gsFillColor   = "Gray"
forrest@0: 
forrest@0: ; Row header (equally space in y2)
forrest@0:   ncr3  = (/nrow_zone,1/)         ; 5 rows, 1 columns
forrest@0:   x3    = x1                      ; same as x1
forrest@0:   y3    = (/1.0-table_length_zone,0.900/) ; end at start of y1
forrest@0:   text3 = row_header_zone
forrest@0:   res3               = True
forrest@0:   res3@txFontHeightF = 0.02
forrest@0:   res3@gsFillColor   = "Gray"
forrest@0: 
forrest@0: ; Main table body
forrest@0:   ncr4  = (/nrow_zone,ncol_zone/) ; 5 rows, 5 columns
forrest@0:   x4    = x2                      ; Start and end x
forrest@0:   y4    = y3                      ; Start and end Y
forrest@0:   text4 = new((/nrow_zone,ncol_zone/),string)
forrest@0: 
forrest@0:   color_fill4      = new((/nrow_zone,ncol_zone/),string)
forrest@0:   color_fill4      = "white"
forrest@0:   color_fill4(:,ncol_zone-1) = "grey"
forrest@0: 
forrest@0:   res4               = True       ; Set up resource list
forrest@0: ; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.
forrest@0:   res4@txFontHeightF = 0.02
forrest@0:   res4@gsFillColor   = color_fill4
forrest@0: 
forrest@0:   delete (color_fill4)
forrest@0: 
forrest@0: ;*******************************************************************
forrest@0: ; for table -- station
forrest@0: ;*******************************************************************
forrest@0: ; column (not including header column)
forrest@0:   col_header_sta = (/"Latitude","Longitude","Amplitude Ratio", \
forrest@0:                       "Correlation Coef","M score","Combined Score"/)
forrest@0:   ncol_sta       = dimsizes(col_header_sta ) 
forrest@0: ;-------------------------------------------------------------------
forrest@0: 
forrest@0: ; Table header
forrest@0:   ncr5  = (/1,1/)               ; 1 row, 1 column
forrest@0:   x5    = (/0.005,0.15/)        ; Start and end X
forrest@0:   y5    = (/0.900,0.995/)       ; Start and end Y
forrest@0:   res5               = True
forrest@0:   res5@txFontHeightF = 0.02
forrest@0:   res5@gsFillColor   = "CornFlowerBlue"
forrest@0: 
forrest@0: ; Column header (equally space in x2)
forrest@0:   ncr6  = (/1,ncol_sta/)         ; 1 rows, 5 columns
forrest@0:   x6    = (/x5(1),0.995/)          ; start from end of x1
forrest@0:   y6    = y5                      ; same as y1
forrest@0:   text6 = col_header_sta
forrest@0:   res6               = True
forrest@0:   res6@txFontHeightF = 0.012
forrest@0:   res6@gsFillColor   = "Gray"
forrest@0: 
forrest@0: ; Row header (equally space in y2)
forrest@0: 
forrest@0:   res7               = True
forrest@0:   res7@txFontHeightF = 0.015
forrest@0:   res7@gsFillColor   = "Gray"
forrest@0: ;--------------------------------------------------------------
forrest@0: ; for station line plot
forrest@0: 
forrest@0: ; for x-axis in xyplot
forrest@0:   mon = ispan(1,12,1)
forrest@0:   mon@long_name = "month"
forrest@0: 
forrest@0:   plot_type = "ps"
forrest@0:   plot_type_new = "png"
forrest@0: 
forrest@0:   res                   = True                      ; plot mods desired
forrest@0:   res@xyLineThicknesses = (/2.0,2.0,2.0/)           ; make 2nd lines thicker
forrest@0:   res@xyLineColors      = (/"red","black"/)  ; change line color
forrest@0: 
forrest@0: ; Add a boxed legend using the more simple method
forrest@0: 
forrest@0:   res@pmLegendDisplayMode    = "Always"
forrest@0: ; res@pmLegendWidthF         = 0.1
forrest@0:   res@pmLegendWidthF         = 0.08
forrest@0:   res@pmLegendHeightF        = 0.06
forrest@0: ; res@pmLegendOrthogonalPosF = -1.17
forrest@0: ; res@pmLegendOrthogonalPosF = -1.00  ;(downward)
forrest@0:   res@pmLegendOrthogonalPosF = -0.30  ;(downward)
forrest@0: 
forrest@0: ; res@pmLegendParallelPosF   =  0.18
forrest@0:   res@pmLegendParallelPosF   =  0.23  ;(rightward)
forrest@0: 
forrest@0: ; res@lgPerimOn             = False
forrest@0:   res@lgLabelFontHeightF     = 0.015
forrest@0:   res@xyExplicitLegendLabels = (/"b30.061n","observed"/)
forrest@0: ;-------------------------------------------------------------------
forrest@0: 
forrest@0:   do z = 0,nzone-1
forrest@0: 
forrest@0:   if (z .eq. 0) then 
forrest@0: ;    maximum score for the zone, 60N-90N
forrest@0:      zone = "60N-90N" 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 60.)
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 1) then 
forrest@0: ;    maximum score for the zone, 30N-60N
forrest@0:      zone = "30N-60N" 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 2) then 
forrest@0: ;    maximum score for the zone, EQ-30N 
forrest@0:      zone = "EQ-30N"
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 3) then 
forrest@0: ;    maximum score for the zone, 90S-EQ
forrest@0:      zone = "90S-EQ" 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .lt. 0. )
forrest@0:   end if
forrest@0: 
forrest@0:  npts = dimsizes(ind_z)
forrest@0: 
forrest@0: ;-------------------------------------------------------------------------
forrest@0: ; for metric table computation
forrest@0: 
forrest@0:  amp_ob        = new((/npts/),float)
forrest@0:  amp_model     = new((/npts/),float)
forrest@0: 
forrest@0:  amp_ratio_sta = new((/npts/),float)
forrest@0:  ccr_sta       = new((/npts/),float)
forrest@0:  M_sta         = new((/npts/),float)
forrest@0:  score_sta     = new((/npts/),float)
forrest@0: 
forrest@0: ;----------------------
forrest@0: ; for table -- station
forrest@0: ;----------------------
forrest@0: 
forrest@0: ; row (not including header row)  
forrest@0:   nrow_sta       = npts                  
forrest@0: 
forrest@0: ; Table header
forrest@0:   ncr5  = (/1,1/)               ; 1 row, 1 column
forrest@0:   x5    = (/0.005,0.15/)        ; Start and end X
forrest@0:   y5    = (/0.900,0.995/)       ; Start and end Y
forrest@0:   res5               = True
forrest@0:   res5@txFontHeightF = 0.02
forrest@0:   res5@gsFillColor   = "CornFlowerBlue"
forrest@0: 
forrest@0: ; Column header (equally space in x2)
forrest@0:   ncr6  = (/1,ncol_sta/)         ; 1 rows, 5 columns
forrest@0:   x6    = (/x1(1),0.995/)          ; start from end of x1
forrest@0:   y6    = y1                      ; same as y1
forrest@0:   text6 = col_header_sta
forrest@0:   res6               = True
forrest@0:   res6@txFontHeightF = 0.012
forrest@0:   res6@gsFillColor   = "Gray"
forrest@0: 
forrest@0: ; Row header (equally space in y2)
forrest@0:   ncr7  = (/nrow_sta,1/)         ; 5 rows, 1 columns
forrest@0:   x7    = x1                      ; same as x1
forrest@0: ; y7    = (/1.0-table_length_sta,0.900/) ; end at start of y1
forrest@0:   text7 = new((/nrow_sta/),string)
forrest@0:   res7               = True
forrest@0:   res7@txFontHeightF = 0.015
forrest@0:   res7@gsFillColor   = "Gray"
forrest@0: ;-------------------------------------------------------------------------
forrest@0: ; for station line plot
forrest@0: 
forrest@0:   npts_str = ""
forrest@0:   npts_str = npts
forrest@0: 
forrest@0:   plot_data   = new((/2,12,npts/),float)
forrest@0:   plot_data_0 = new((/12,npts/),float)
forrest@0: 
forrest@0:   plot_data!0 = "case"
forrest@0:   plot_data!1 = "month"
forrest@0:   plot_data!2 = "pts"
forrest@0:   plot_data@long_name   = "CO2 Seasonal"
forrest@0: 
forrest@0:   plot_data_0!0 = "month"
forrest@0:   plot_data_0!1 = "pts"
forrest@0:   plot_data_0@long_name = "CO2"
forrest@0: ;--------------------------------------------------------------------------
forrest@0:  do n=0,npts-1
forrest@0:     amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:)) 
forrest@0:     amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
forrest@0: 
forrest@0:     amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
forrest@0:     ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
forrest@0:     M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
forrest@0:     score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
forrest@0: 
forrest@0: ;----------------------------------------------------------------------
forrest@0: ; for station line plot
forrest@0: 
forrest@0:     plot_data(0,:,n) = (/val_model(ind_z(n),:)/)
forrest@0:     plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)
forrest@0: 
forrest@0:     plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)
forrest@0:    
forrest@0:     plot_name = sta(ind_z(n))    
forrest@0:     title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"
forrest@0: 
forrest@0:     wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0: ;------------------------------------------
forrest@0: ;   for panel plot
forrest@0:    
forrest@0:     plot=new(2,graphic)                        ; create graphic array
forrest@0:     res@gsnFrame     = False                   ; Do not draw plot 
forrest@0:     res@gsnDraw      = False                   ; Do not advance frame
forrest@0: 
forrest@0:     pres                            = True     ; panel plot mods desired
forrest@0:     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
forrest@0:                                                ; indiv. plots in panel
forrest@0:     pres@gsnMaximize                = True     ; fill the page
forrest@0: ;------------------------------------------
forrest@0:     res@tiMainString = title                           ; add title
forrest@0: 
forrest@0:     plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res)   ; create plot 1
forrest@0: 
forrest@0:     plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2
forrest@0: 
forrest@0:     gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot
forrest@0: 
forrest@0:     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:     system("rm "+plot_name+"."+plot_type)
forrest@0: 
forrest@0:     clear (wks)
forrest@0: ;---------------------------------------------------------------------------  
forrest@0:  end do
forrest@0: ;-------------------------------------------------------------------------
forrest@0: ;   for line plot in a zone
forrest@0: 
forrest@0:     plot_name = "All_"+npts_str
forrest@0:     title = plot_name + " in "+ zone
forrest@0: 
forrest@0:     wks = gsn_open_wks (plot_type,plot_name)        ; open workstation
forrest@0: ;-----------------------------------------
forrest@0: ;   for panel plot
forrest@0:    
forrest@0:     plot=new(2,graphic)                        ; create graphic array
forrest@0:     res@gsnFrame     = False                   ; Do not draw plot 
forrest@0:     res@gsnDraw      = False                   ; Do not advance frame
forrest@0: 
forrest@0:     pres                            = True     ; panel plot mods desired
forrest@0:     pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around
forrest@0:                                                ; indiv. plots in panel
forrest@0:     pres@gsnMaximize                = True     ; fill the page
forrest@0: ;-----------------------------------------
forrest@0:     res@tiMainString = title                                     ; add title
forrest@0: 
forrest@0:     plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res)   ; create plot 1
forrest@0:     
forrest@0:     plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2
forrest@0: 
forrest@0:     gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot
forrest@0: 
forrest@0:     system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:     system("rm "+plot_name+"."+plot_type)
forrest@0: 
forrest@0:     clear (wks)
forrest@0: ;   delete (ind_z)
forrest@0:     delete (plot_data)
forrest@0:     delete (plot_data_0)    
forrest@0: ;---------------------------------------------------------------------------
forrest@0: ; for zone table value
forrest@0: 
forrest@0:  amp_ratio_zone = avg(amp_ratio_sta)
forrest@0:  ccr_zone       = avg(ccr_sta)
forrest@0:  M_zone   = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts) 
forrest@0:  score_zone   = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
forrest@0: 
forrest@0:   text4(z,0) = sprintf("%5.2f", npts)
forrest@0:   text4(z,1) = sprintf("%5.2f", amp_ratio_zone)
forrest@0:   text4(z,2) = sprintf("%5.2f", ccr_zone)
forrest@0:   text4(z,3) = sprintf("%5.2f", M_zone)
forrest@0:   text4(z,4) = sprintf("%5.2f", score_zone)  
forrest@0: ;---------------------------------------------------------------------------
forrest@0: ; plot station table
forrest@0: ;----------------------------------
forrest@0: ; header value for station table
forrest@0:   text5  = zone 
forrest@0: 
forrest@0: ; row value for station table 
forrest@0:   text7  = sta(ind_z)
forrest@0: 
forrest@0:   if (z .eq. 0) then
forrest@0:      table_length_sta = 0.5
forrest@0:   end if
forrest@0:   if (z .eq. 1) then
forrest@0:      table_length_sta = 0.995
forrest@0:   end if
forrest@0:   if (z .eq. 2) then
forrest@0:      table_length_sta = 0.8
forrest@0:   end if
forrest@0:   if (z .eq. 3) then
forrest@0:      table_length_sta = 0.8
forrest@0:   end if    
forrest@0: 
forrest@0:   x7 = x5   
forrest@0:   y7 = (/1.0-table_length_sta,0.900/) ; end at start of y1
forrest@0: 
forrest@0: ; Main table body
forrest@0:   ncr8  = (/nrow_sta,ncol_sta/) ; 5 rows, 5 columns
forrest@0:   x8    = x6                      ; Start and end x
forrest@0:   y8    = y7                      ; Start and end Y
forrest@0:   text8 = new((/nrow_sta,ncol_sta/),string)
forrest@0: 
forrest@0:   color_fill8      = text8
forrest@0:   color_fill8      = "white"
forrest@0:   color_fill8(:,ncol_sta-1) = "grey"
forrest@0: 
forrest@0:   res8               = True       ; Set up resource list
forrest@0: ; res8@gsnDebug      = True       ; Useful to print NDC row,col values used.
forrest@0:   res8@txFontHeightF = 0.02
forrest@0:   res8@gsFillColor   = color_fill8
forrest@0: 
forrest@0:   delete (color_fill8)
forrest@0: 
forrest@0: ; table value for station table
forrest@0:   text8(:,0) = sprintf("%5.2f", (/lat(ind_z)/))
forrest@0:   text8(:,1) = sprintf("%5.2f", (/lon(ind_z)/))
forrest@0:   text8(:,2) = sprintf("%5.2f", (/amp_ratio_sta/))
forrest@0:   text8(:,3) = sprintf("%5.2f", (/ccr_sta/))
forrest@0:   text8(:,4) = sprintf("%5.2f", (/M_sta/))
forrest@0:   text8(:,5) = sprintf("%5.2f", (/score_sta/))
forrest@0:  
forrest@0:   plot_name = "table_sta." + zone 
forrest@0:   
forrest@0:   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0: 
forrest@0:   gsn_table(wks,ncr5,x5,y5,text5,res5)
forrest@0:   gsn_table(wks,ncr6,x6,y6,text6,res6)
forrest@0:   gsn_table(wks,ncr7,x7,y7,text7,res7)
forrest@0:   gsn_table(wks,ncr8,x8,y8,text8,res8) 
forrest@0: 
forrest@0:   frame(wks)
forrest@0: 
forrest@0:   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:   system("rm "+plot_name+"."+plot_type)
forrest@0: 
forrest@0:  delete (ind_z)
forrest@0:  delete (amp_model)
forrest@0:  delete (amp_ob)
forrest@0:  delete (amp_ratio_sta)
forrest@0:  delete (ccr_sta)
forrest@0:  delete (M_sta)
forrest@0:  delete (score_sta)
forrest@0:  delete (text7)
forrest@0:  delete (text8)
forrest@0:  delete (res7)
forrest@0:  delete (res8)
forrest@0:  clear (wks)
forrest@0:  end do
forrest@0: ;**************************************************
forrest@0: ; plot zone table
forrest@0: ;**************************************************
forrest@0:  
forrest@0:   text4(4,0) = sum(stringtofloat(text4(0:3,0))) 
forrest@0:   text4(4,1) = 0.
forrest@0:   text4(4,2) = 0.
forrest@0:   text4(4,3) = 0.
forrest@0:   text4(4,4) = sum(stringtofloat(text4(0:3,4)))
forrest@0: 
forrest@0:   plot_name = "table_zone" 
forrest@0:   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
forrest@0: 
forrest@0:   gsn_table(wks,ncr1,x1,y1,text1,res1)
forrest@0:   gsn_table(wks,ncr2,x2,y2,text2,res2)
forrest@0:   gsn_table(wks,ncr3,x3,y3,text3,res3)
forrest@0:   gsn_table(wks,ncr4,x4,y4,text4,res4) 
forrest@0: 
forrest@0:   frame(wks)
forrest@0: 
forrest@0:   system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
forrest@0:   system("rm "+plot_name+"."+plot_type)
forrest@0: 
forrest@0:   M_co2 = text4(4,4)
forrest@0: 
forrest@0:   if (isvar("compare")) then
forrest@0:      system("sed 1,/M_co2/s//"+M_co2+"/ "+html_name2+" > "+html_new2+";"+ \
forrest@0:             "mv -f "+html_new2+" "+html_name2)
forrest@0:   end if
forrest@0: 
forrest@0:   system("sed s#M_co2#"+M_co2+"# "+html_name+" > "+html_new+";"+ \
forrest@0:          "mv -f "+html_new+" "+html_name) 
forrest@0: ;***************************************************************************
forrest@0: ; output plots
forrest@0: ;***************************************************************************
forrest@0:   output_dir = model_name+"/co2"
forrest@0: 
forrest@0:   system("mv *.png " + output_dir) 
forrest@0: ;*************************************************************************** 
forrest@0: end