; ***********************************************

; using gsn_table for all

; combine 19.metric_plot.ncl and 24.lines.ncl

; ***********************************************

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

;load "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"

;************************************************

begin

;************************************************

; read in model data

;************************************************

  diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

; fili2 = "b30.061m_401_425_MONS_climo_atm.nc"

  fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc"

  g     = addfile(diri2+fili2,"r")

  x     = g->CO2

  xi    = g->lon

  yi    = g->lat

  xdim  = dimsizes(x)

  nlev  = xdim(1)

  y     = x(:,0,:,:)

; printVarSummary (y)

; get the co2 at the lowest level

  y     = x(:,nlev-1,:,:)

; change to unit of observed (u mol/mol)

; Model_units [=] kgCO2 / kgDryAir

; 28.966 = molecular weight of dry air

; 44.       = molecular weight of CO2

; u mol = 1e-6 mol

  factor = (28.966/44.) * 1e6

  y      = y * factor

  y@_FillValue = 1.e36

  y@units      = "u mol/mol"

; y = where(y0 .lt. 287.,y@_FillValue,y)

; printVarSummary (y)

; print (min(y)+"/"+max(y))

;************************************************

; read in data: observed

;************************************************

 diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/ob/"

 fili  = "co2_globalView_98.nc"

 g     = addfile (diri+fili,"r")

 val   = g->CO2_SEAS  

 lon   = g->LON 

 lat   = g->LAT

 sta   = chartostring(g->STATION) 

 delete (g)

;print (sta(0))

 ncase = dimsizes(lat)

;print (ncase)

;**************************************************************

; get only the lowest level at each station 

;**************************************************************

 lat_tmp = lat

 lat_tmp@_FillValue = 1.e+36

 do n = 0,ncase-1

    if (.not. ismissing(lat_tmp(n))) then 

       indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)

       if (dimsizes(indexes) .gt. 1) then

          lat_tmp(indexes(1:)) = lat_tmp@_FillValue

       end if

       delete (indexes)

    end if

 end do

 indexes = ind(.not. ismissing(lat_tmp))

;print (dimsizes(indexes))

;print (indexes)

 lat_ob = lat(indexes)

 lon_ob = lon(indexes)

 val_ob = val(indexes,:)

;printVarSummary (val_ob)

;print (lat_ob +"/"+lon_ob)

;************************************************************

; interpolate model data into observed station

; note: model is 0-360E, 90S-90N

;************************************************************

; to be able to handle observation at (-89.98,-24.80)

  print (yi(0))

  yi(0) = -90.

  i = ind(lon_ob .lt. 0.)

  lon_ob(i) = lon_ob(i) + 360.  

  yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)

  val_model = yo(pts|:,time|:)

  val_model_0 = val_model

; printVarSummary (val_model)

; print (min(val_model)+"/"+max(val_model))

; remove annual mean

  val_model = val_model - conform(val_model,dim_avg(val_model),0)

; print (min(val_model)+"/"+max(val_model))

  nzone = 4

;*******************************************************************

; for table -- zone

;*******************************************************************

; table header name

  table_header_name = "Zone" 

; column (not including header column)

  col_header_zone = (/"Stations","Amplitude Ratio", \

                      "Correlation Coef","M score","Combined Score"/)

  ncol_zone       = dimsizes(col_header_zone ) 

; row (not including header row)

  row_header_zone = (/"60N-90N","30N-60N","EQ-30N","90S-EQ","Total"/)  

  nrow_zone       = dimsizes(row_header_zone)                  

; arrays to be passed to table. 

  value_zone = new ((/nrow_zone, ncol_zone/),string ) 

;--------------------------------------------------------------------

  table_length_zone = 0.4 

; Table header

  ncr1  = (/1,1/)               ; 1 row, 1 column

  x1    = (/0.005,0.15/)        ; Start and end X

  y1    = (/0.900,0.995/)       ; Start and end Y

  text1 = table_header_name

  res1               = True

  res1@txFontHeightF = 0.03

  res1@gsFillColor   = "CornFlowerBlue"

; Column header (equally space in x2)

  ncr2  = (/1,ncol_zone/)         ; 1 rows, 5 columns

  x2    = (/x1(1),0.995/)          ; start from end of x1

  y2    = y1                      ; same as y1

  text2 = col_header_zone

  res2               = True

  res2@txFontHeightF = 0.015

  res2@gsFillColor   = "Gray"

; Row header (equally space in y2)

  ncr3  = (/nrow_zone,1/)         ; 5 rows, 1 columns

  x3    = x1                      ; same as x1

  y3    = (/1.0-table_length_zone,0.900/) ; end at start of y1

  text3 = row_header_zone

  res3               = True

  res3@txFontHeightF = 0.02

  res3@gsFillColor   = "Gray"

; Main table body

  ncr4  = (/nrow_zone,ncol_zone/) ; 5 rows, 5 columns

  x4    = x2                      ; Start and end x

  y4    = y3                      ; Start and end Y

  text4 = new((/nrow_zone,ncol_zone/),string)

  color_fill4      = new((/nrow_zone,ncol_zone/),string)

  color_fill4      = "white"

  color_fill4(:,ncol_zone-1) = "grey"

  res4               = True       ; Set up resource list

; res4@gsnDebug      = True       ; Useful to print NDC row,col values used.

  res4@txFontHeightF = 0.02

  res4@gsFillColor   = color_fill4

  delete (color_fill4)

;*******************************************************************

; for table -- station

;*******************************************************************

; column (not including header column)

  col_header_sta = (/"Latitude","Longitude","Amplitude Ratio", \

                      "Correlation Coef","M score","Combined Score"/)

  ncol_sta       = dimsizes(col_header_sta ) 

;-------------------------------------------------------------------

; Table header

  ncr5  = (/1,1/)               ; 1 row, 1 column

  x5    = (/0.005,0.15/)        ; Start and end X

  y5    = (/0.900,0.995/)       ; Start and end Y

  res5               = True

  res5@txFontHeightF = 0.02

  res5@gsFillColor   = "CornFlowerBlue"

; Column header (equally space in x2)

  ncr6  = (/1,ncol_sta/)         ; 1 rows, 5 columns

  x6    = (/x5(1),0.995/)          ; start from end of x1

  y6    = y5                      ; same as y1

  text6 = col_header_sta

  res6               = True

  res6@txFontHeightF = 0.012

  res6@gsFillColor   = "Gray"

; Row header (equally space in y2)

  res7               = True

  res7@txFontHeightF = 0.015

  res7@gsFillColor   = "Gray"

;--------------------------------------------------------------

; for station line plot

; for x-axis in xyplot

  mon = ispan(1,12,1)

  mon@long_name = "month"

  plot_type = "ps"

  plot_type_new = "png"

  res                   = True                      ; plot mods desired

  res@xyLineThicknesses = (/2.0,2.0,2.0/)           ; make 2nd lines thicker

  res@xyLineColors      = (/"red","black"/)  ; change line color

; Add a boxed legend using the more simple method

  res@pmLegendDisplayMode    = "Always"

; res@pmLegendWidthF         = 0.1

  res@pmLegendWidthF         = 0.08

  res@pmLegendHeightF        = 0.06

; res@pmLegendOrthogonalPosF = -1.17

; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

  res@pmLegendOrthogonalPosF = -0.30  ;(downward)

; res@pmLegendParallelPosF   =  0.18

  res@pmLegendParallelPosF   =  0.23  ;(rightward)

; res@lgPerimOn             = False

  res@lgLabelFontHeightF     = 0.015

  res@xyExplicitLegendLabels = (/"b30.061n","observed"/)

;-------------------------------------------------------------------

  do z = 0,nzone-1

  if (z .eq. 0) then 

;    maximum score for the zone, 60N-90N

     zone = "60N-90N" 

     score_max = 5.0

;    index of stations in this zone

     ind_z = ind(lat_ob .ge. 60.)

;    print (ind_z)

;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

;    print (val_ob(ind_z,:))

;    print (val_model(ind_z,:))

  end if

  if (z .eq. 1) then 

;    maximum score for the zone, 30N-60N

     zone = "30N-60N" 

     score_max = 5.0

;    index of stations in this zone

     ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)

;    print (ind_z)

;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

;    print (val_ob(ind_z,:))

;    print (val_model(ind_z,:))

  end if

  if (z .eq. 2) then 

;    maximum score for the zone, EQ-30N 

     zone = "EQ-30N"

     score_max = 5.0

;    index of stations in this zone

     ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)

;    print (ind_z)

;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

;    print (val_ob(ind_z,:))

;    print (val_model(ind_z,:))

  end if

  if (z .eq. 3) then 

;    maximum score for the zone, 90S-EQ

     zone = "90S-EQ" 

     score_max = 5.0

;    index of stations in this zone

     ind_z = ind(lat_ob .lt. 0. )

;    print (ind_z)

;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

;    print (val_ob(ind_z,:))

;    print (val_model(ind_z,:))

  end if

 npts = dimsizes(ind_z)

 print (npts)

;-------------------------------------------------------------------------

; for metric table computation

 amp_ob        = new((/npts/),float)

 amp_model     = new((/npts/),float)

 amp_ratio_sta = new((/npts/),float)

 ccr_sta       = new((/npts/),float)

 M_sta         = new((/npts/),float)

 score_sta     = new((/npts/),float)

;----------------------

; for table -- station

;----------------------

; row (not including header row)  

  nrow_sta       = npts                  

; Table header

  ncr5  = (/1,1/)               ; 1 row, 1 column

  x5    = (/0.005,0.15/)        ; Start and end X

  y5    = (/0.900,0.995/)       ; Start and end Y

  res5               = True

  res5@txFontHeightF = 0.02

  res5@gsFillColor   = "CornFlowerBlue"

; Column header (equally space in x2)

  ncr6  = (/1,ncol_sta/)         ; 1 rows, 5 columns

  x6    = (/x1(1),0.995/)          ; start from end of x1

  y6    = y1                      ; same as y1

  text6 = col_header_sta

  res6               = True

  res6@txFontHeightF = 0.012

  res6@gsFillColor   = "Gray"

; Row header (equally space in y2)

  ncr7  = (/nrow_sta,1/)         ; 5 rows, 1 columns

  x7    = x1                      ; same as x1

; y7    = (/1.0-table_length_sta,0.900/) ; end at start of y1

  text7 = new((/nrow_sta/),string)

  res7               = True

  res7@txFontHeightF = 0.015

  res7@gsFillColor   = "Gray"

;-------------------------------------------------------------------------

; for station line plot

  npts_str = ""

  npts_str = npts

; print (npts_str)

  plot_data   = new((/2,12,npts/),float)

  plot_data_0 = new((/12,npts/),float)

  plot_data!0 = "case"

  plot_data!1 = "month"

  plot_data!2 = "pts"

  plot_data@long_name   = "CO2 Seasonal"

  plot_data_0!0 = "month"

  plot_data_0!1 = "pts"

  plot_data_0@long_name = "CO2"

;--------------------------------------------------------------------------

 do n=0,npts-1

    amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:)) 

    amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))

    amp_ratio_sta(n) = amp_model(n)/amp_ob(n)

    ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)

    M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)

    score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max

    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))

;----------------------------------------------------------------------

; for station line plot

    plot_data(0,:,n) = (/val_model(ind_z(n),:)/)

    plot_data(1,:,n) = (/val_ob(ind_z(n),:)/)

    plot_data_0(:,n) = (/val_model_0(ind_z(n),:)/)

    plot_name = sta(ind_z(n))    

    title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"

;   print (title)

;   print (plot_name)

    wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

;------------------------------------------

;   for panel plot

    plot=new(2,graphic)                        ; create graphic array

    res@gsnFrame     = False                   ; Do not draw plot 

    res@gsnDraw      = False                   ; Do not advance frame

    pres                            = True     ; panel plot mods desired

    pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

                                               ; indiv. plots in panel

    pres@gsnMaximize                = True     ; fill the page

;------------------------------------------

    res@tiMainString = title                           ; add title

    plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res)   ; create plot 1

    plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,n),res) ; create plot 2

    gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot

    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

    system("rm "+plot_name+"."+plot_type)

    clear (wks)

;---------------------------------------------------------------------------  

 end do

;-------------------------------------------------------------------------

;   for line plot in a zone

    plot_name = "All_"+npts_str

    title = plot_name + " in "+ zone

;   print (title)

;   print (plot_name)

    wks = gsn_open_wks (plot_type,plot_name)        ; open workstation

;-----------------------------------------

;   for panel plot

    plot=new(2,graphic)                        ; create graphic array

    res@gsnFrame     = False                   ; Do not draw plot 

    res@gsnDraw      = False                   ; Do not advance frame

    pres                            = True     ; panel plot mods desired

    pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

                                               ; indiv. plots in panel

    pres@gsnMaximize                = True     ; fill the page

;-----------------------------------------

    res@tiMainString = title                                     ; add title

    plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res)   ; create plot 1

    plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2

    gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot

    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

    system("rm "+plot_name+"."+plot_type)

    clear (wks)

;   delete (ind_z)

    delete (plot_data)

    delete (plot_data_0)    

;---------------------------------------------------------------------------

; for zone table value

 amp_ratio_zone = avg(amp_ratio_sta)

 ccr_zone       = avg(ccr_sta)

 M_zone   = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts) 

 score_zone   = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max

 print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)

  text4(z,0) = sprintf("%5.2f", npts)

  text4(z,1) = sprintf("%5.2f", amp_ratio_zone)

  text4(z,2) = sprintf("%5.2f", ccr_zone)

  text4(z,3) = sprintf("%5.2f", M_zone)

  text4(z,4) = sprintf("%5.2f", score_zone)  

;---------------------------------------------------------------------------

; plot station table

;----------------------------------

; header value for station table

  text5  = zone 

; row value for station table 

  text7  = sta(ind_z)

  if (z .eq. 0) then

     table_length_sta = 0.5

  end if

  if (z .eq. 1) then

     table_length_sta = 0.995

  end if

  if (z .eq. 2) then

     table_length_sta = 0.8

  end if

  if (z .eq. 3) then

     table_length_sta = 0.8

  end if    

  x7 = x5   

  y7 = (/1.0-table_length_sta,0.900/) ; end at start of y1

; Main table body

  ncr8  = (/nrow_sta,ncol_sta/) ; 5 rows, 5 columns

  x8    = x6                      ; Start and end x

  y8    = y7                      ; Start and end Y

  text8 = new((/nrow_sta,ncol_sta/),string)

  color_fill8      = text8

  color_fill8      = "white"

  color_fill8(:,ncol_sta-1) = "grey"

  res8               = True       ; Set up resource list

  res8@gsnDebug      = True       ; Useful to print NDC row,col values used.

  res8@txFontHeightF = 0.02

  res8@gsFillColor   = color_fill8

  delete (color_fill8)

; table value for station table

  text8(:,0) = sprintf("%5.2f", (/lat(ind_z)/))

  text8(:,1) = sprintf("%5.2f", (/lon(ind_z)/))

  text8(:,2) = sprintf("%5.2f", (/amp_ratio_sta/))

  text8(:,3) = sprintf("%5.2f", (/ccr_sta/))

  text8(:,4) = sprintf("%5.2f", (/M_sta/))

  text8(:,5) = sprintf("%5.2f", (/score_sta/))

  plot_name = "table_sta." + zone 

  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  gsn_table(wks,ncr5,x5,y5,text5,res5)

  gsn_table(wks,ncr6,x6,y6,text6,res6)

  gsn_table(wks,ncr7,x7,y7,text7,res7)

  gsn_table(wks,ncr8,x8,y8,text8,res8) 

  frame(wks)

  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  system("rm "+plot_name+"."+plot_type)

 delete (ind_z)

 delete (amp_model)

 delete (amp_ob)

 delete (amp_ratio_sta)

 delete (ccr_sta)

 delete (M_sta)

 delete (score_sta)

 delete (text7)

 delete (text8)

 delete (res7)

 delete (res8)

 clear (wks)

 end do

;**************************************************

; plot zone table

;**************************************************

  text4(4,0) = sum(stringtofloat(text4(0:3,0))) 

  text4(4,1) = 0.

  text4(4,2) = 0.

  text4(4,3) = 0.

  text4(4,4) = sum(stringtofloat(text4(0:3,4)))

  plot_name = "table_zone" 

  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  gsn_table(wks,ncr1,x1,y1,text1,res1)

  gsn_table(wks,ncr2,x2,y2,text2,res2)

  gsn_table(wks,ncr3,x3,y3,text3,res3)

  gsn_table(wks,ncr4,x4,y4,text4,res4) 

  frame(wks)

  system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

  system("rm "+plot_name+"."+plot_type)

;-------------------------------------------------------------------

  system("mkdir temp")

  system("mv *.png temp")

  system("tar cf temp.tar temp")

;------------------------------------------------------------------- 

end