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

 ; combine scatter, histogram, global and zonal plots

 ; compute all correlation coef and M score

 ; add 1-to-1 line in scatter plots

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

 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"

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

 procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \

                  ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \

                  ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \

                  ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \

                  ,dx4[1]:numeric \

                   )

 begin

 ; Calculaee "nice" bins for binning the data in equally spaced ranges.

 ; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D

 ; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4

   nbins       = 15     ; Number of bins to use.

   nicevals    = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)

   nvals       = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)

   range       = fspan(nicevals(0),nicevals(1),nvals)

 ; Use this range information to grab all the values in a

 ; particular range, and then take an average.

   nr           = dimsizes(range)

   nx           = nr-1

 ; print (nx)

 ; xvalues      = new((/2,nx/),typeof(RAIN1_1D))

   xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.

   dx           = xvalues(0,1) - xvalues(0,0)       ; range width

   dx4          = dx/4                              ; 1/4 of the range

   xvalues(1,:) = xvalues(0,:) - dx/5.

 ; yvalues      = new((/2,nx/),typeof(RAIN1_1D))

 ; mn_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

 ; mx_yvalues   = new((/2,nx/),typeof(RAIN1_1D))

   do nd=0,1

 ; See if we are doing model or observational data.

     if(nd.eq.0) then

       data     = RAIN1_1D

       npp_data = NPP1_1D

     else

       data     = RAIN2_1D

       npp_data = NPP2_1D

     end if

 ; Loop through each range and check for values.

     do i=0,nr-2

       if (i.ne.(nr-2)) then

 ;        print("")

 ;        print("In range ["+range(i)+","+range(i+1)+")")

         idx = ind((range(i).le.data).and.(data.lt.range(i+1)))

       else

 ;        print("")

 ;        print("In range ["+range(i)+",)")

         idx = ind(range(i).le.data)

       end if

 ; Calculate average, and get min and max.

       if(.not.any(ismissing(idx))) then

         yvalues(nd,i)    = avg(npp_data(idx))

         mn_yvalues(nd,i) = min(npp_data(idx))

         mx_yvalues(nd,i) = max(npp_data(idx))

         count = dimsizes(idx)

       else

         count            = 0

         yvalues(nd,i)    = yvalues@_FillValue

         mn_yvalues(nd,i) = yvalues@_FillValue

         mx_yvalues(nd,i) = yvalues@_FillValue

       end if

 ; Print out information.

 ;      print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))

 ;      print("Min/Max:  " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))

 ; Clean up for next time in loop.

       delete(idx)

     end do

     delete(data)

     delete(npp_data)

   end do 

 end

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

 ; Main code.

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

 begin

  plot_type     = "ps"

  plot_type_new = "png"

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

 ; read model data

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

 ;film = "i01.06cn_1798-2004_ANN_climo.nc"

 ;model_name = "06cn"

 ;model_grid = "T42"

 ;film = "i01.06casa_1798-2004_ANN_climo.nc"

 ;model_name = "06casa"

 ;model_grid = "T42"

  film = "i01.10casa_1948-2004_ANN_climo.nc"

  model_name = "10casa"

  model_grid = "T42"

 ;film = "i01.10cn_1948-2004_ANN_climo.nc"

 ;model_name = "10cn"

 ;model_grid = "T42"

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

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

  fm   = addfile (dirm+film,"r")

  nppmod0  = fm->NPP

  rainmod0 = fm->RAIN

  xm       = fm->lon  

  ym       = fm->lat

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

 ; read ob data

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

 ;(1) data at 81 sites

  dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

  fil81 = "data.81.nc"

  f81   = addfile (dir81+fil81,"r")

  id81   = f81->SITE_ID  

  npp81  = f81->TNPP_C

  rain81 = tofloat(f81->PREC_ANN)

  x81    = f81->LONG_DD  

  y81    = f81->LAT_DD

  id81@long_name  = "SITE_ID"

 ;change longitude from (-180,180) to (0,360)

 ;for model data interpolation

  do i= 0,dimsizes(x81)-1

     if (x81(i) .lt. 0.) then

         x81(i) = x81(i)+ 360.

     end if

  end do

 ;print (x81)

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

 ;(2) data at 933 sites

  dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

  fil933 = "data.933.nc"

  f933   = addfile (dir933+fil933,"r")

  id933   = f933->SITE_ID  

  npp933  = f933->TNPP_C

  rain933 = f933->PREC

  x933    = f933->LONG_DD  

  y933    = f933->LAT_DD 

  id933@long_name  = "SITE_ID"

 ;change longitude from (-180,180) to (0,360)

 ;for model data interpolation

  do i= 0,dimsizes(x933)-1

     if (x933(i) .lt. 0.) then

         x933(i) = x933(i)+ 360.

     end if

  end do

 ;print (x933)

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

 ;(3) global data, interpolated into model grid

  dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"

  filglobe = "Npp_"+model_grid+"_mean.nc"

  fglobe   = addfile (dirglobe+filglobe,"r")

  nppglobe0 = fglobe->NPP

  nppglobe  = nppglobe0

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

 ; interpolate model data into ob sites

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

  nppmod   = nppmod0(0,:,:)

  rainmod  = rainmod0(0,:,:)

  delete (nppmod0)

  delete (rainmod0)

  nppmod81  =linint2_points(xm,ym,nppmod,True,x81,y81,0)

  nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)

  rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)

  rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)

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

 ; unified units

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

 ; Units for these variables are:

 ;

 ; rain81  : mm/year

 ; rainmod : mm/s

 ; npp81   : g C/m^2/year

 ; nppmod81: g C/m^2/s

 ; nppglobe: g C/m^2/year

 ;

 ; We want to convert these to "m/year" and "g C/m^2/year".

   nsec_per_year = 60*60*24*365                 

   rain81    = rain81 / 1000.

   rainmod81 = (rainmod81/ 1000.) * nsec_per_year

   nppmod81  = nppmod81 * nsec_per_year

   rain933    = rain933 / 1000.

   rainmod933 = (rainmod933/ 1000.) * nsec_per_year

   nppmod933  = nppmod933 * nsec_per_year

   nppmod  = nppmod * nsec_per_year

   npp81@units      = "gC/m^2/yr"

   nppmod81@units   = "gC/m^2/yr"

   npp933@units     = "gC/m^2/yr"

   nppmod933@units  = "gC/m^2/yr"

   nppmod@units     = "gC/m^2/yr"

   nppglobe@units   = "gC/m^2/yr"

   rain81@units     = "m/yr"

   rainmod81@units  = "m/yr"

   rain933@units    = "m/yr"

   rainmod933@units = "m/yr"

   npp81@long_name      = "NPP (gC/m2/year)"

   npp933@long_name     = "NPP (gC/m2/year)"

   nppmod81@long_name   = "NPP (gC/m2/year)"

   nppmod933@long_name  = "NPP (gC/m2/year)"

   nppmod@long_name     = "NPP (gC/m2/year)"

   nppglobe@long_name   = "NPP (gC/m2/year)"

   rain81@long_name     = "PREC (m/year)"

   rain933@long_name    = "PREC (m/year)"

   rainmod81@long_name  = "PREC (m/year)"

   rainmod933@long_name = "PREC (m/year)"

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

 ;(A)-1 table of site81 (1) (the first 41 sites)

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

   table_length = 0.95 

 ; table header name

   table_header_name = "Site ID" 

 ; column (not including header column)

   col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/) 

   ncol = dimsizes(col_header) 

 ; row (not including header row) 

   nrow       = 41

   row_header = new ((/nrow/),string )

   row_header(0:nrow-1) = id81(0:nrow-1) 

 ; Table header

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

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

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

   text1 = table_header_name

   res1               = True

   res1@txFontHeightF = 0.015

   res1@gsFillColor   = "CornFlowerBlue"

 ; Column header (equally space in x2)

   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

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

   y2    = y1                    ; same as y1

   text2 = col_header

   res2               = True

   res2@txFontHeightF = 0.015

   res2@gsFillColor   = "Gray"

 ; Row header (equally space in y2)

   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   x3    = x1                         ; same as x1

   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   text3 = row_header

   res3               = True

   res3@txFontHeightF = 0.010

   res3@gsFillColor   = "Gray"

 ; Main table body

   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   x4    = x2                      ; Start and end x

   y4    = y3                      ; Start and end Y

   text4 = new((/nrow,ncol/),string)

   color_fill4           = new((/nrow,ncol/),string)

   color_fill4           = "white"

 ; color_fill4(:,ncol-1) = "grey"

 ; color_fill4(nrow-1,:) = "green"

   res4               = True       ; Set up resource list

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

   res4@txFontHeightF = 0.015

   res4@gsFillColor   = color_fill4

   delete (color_fill4)

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

 ; for table value

   do n = 0,nrow-1

      text4(n,0) = sprintf("%5.2f", y81(n))  

      text4(n,1) = sprintf("%5.2f", x81(n))    

      text4(n,2) = sprintf("%5.2f", npp81(n))     

      text4(n,3) = sprintf("%5.2f", rain81(n))          

   end do

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

   plot_name = "table_site81_ob1"

   wks = gsn_open_wks (plot_type,plot_name)

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

 ; for table title

   gRes               = True

   gRes@txFontHeightF = 0.02

 ; gRes@txAngleF      = 90

   title_text = "Observation at 81 Sites (1)"

   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

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

   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)

   clear (wks)

   delete (row_header)

   delete (res3)

   delete (ncr3)

   delete (text3)

   delete (res4)

   delete (ncr4)

   delete (text4)

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

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

   system("rm "+plot_name+"-1."+plot_type_new)

   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

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

 ;(A)-2 table of site81 (2) (the last 40 sites)

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

   table_length = 0.95 

 ; table header name

   table_header_name = "Site ID" 

 ; column (not including header column)

   col_header = (/"Latitude","Longitude","NPP (gC/m2/year)","RAIN (m/year)"/) 

   ncol = dimsizes(col_header) 

 ; row (not including header row) 

   nrow       = 40

   row_header = new ((/nrow/),string )

   row_header(0:nrow-1) = id81(41:41+nrow-1) 

 ; Table header

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

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

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

   text1 = table_header_name

   res1               = True

   res1@txFontHeightF = 0.015

   res1@gsFillColor   = "CornFlowerBlue"

 ; Column header (equally space in x2)

   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

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

   y2    = y1                    ; same as y1

   text2 = col_header

   res2               = True

   res2@txFontHeightF = 0.015

   res2@gsFillColor   = "Gray"

 ; Row header (equally space in y2)

   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   x3    = x1                         ; same as x1

   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   text3 = row_header

   res3               = True

   res3@txFontHeightF = 0.010

   res3@gsFillColor   = "Gray"

 ; Main table body

   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   x4    = x2                      ; Start and end x

   y4    = y3                      ; Start and end Y

   text4 = new((/nrow,ncol/),string)

   color_fill4           = new((/nrow,ncol/),string)

   color_fill4           = "white"

 ; color_fill4(:,ncol-1) = "grey"

 ; color_fill4(nrow-1,:) = "green"

   res4               = True       ; Set up resource list

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

   res4@txFontHeightF = 0.015

   res4@gsFillColor   = color_fill4

   delete (color_fill4)

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

 ; for table value

   do n = 0,nrow-1

      text4(n,0) = sprintf("%5.2f", y81(n+41))  

      text4(n,1) = sprintf("%5.2f", x81(n+41))    

      text4(n,2) = sprintf("%5.2f", npp81(n+41))     

      text4(n,3) = sprintf("%5.2f", rain81(n+41))          

   end do

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

   plot_name = "table_site81_ob2"

   wks = gsn_open_wks (plot_type,plot_name)

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

 ; for table title

   gRes               = True

   gRes@txFontHeightF = 0.02

 ; gRes@txAngleF      = 90

   title_text = "Observation at 81 Sites (2)"

   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

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

   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)

   clear (wks)

   delete (row_header)

   delete (res3)

   delete (ncr3)

   delete (text3)

   delete (res4)

   delete (ncr4)

   delete (text4)

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

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

   system("rm "+plot_name+"-1."+plot_type_new)

   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

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

 ;(A)-3 scatter plot, ob 933

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

  plot_name = "scatter_ob_933"

  title     = "Observed 933 sites"

  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

  res                   = True                  ; plot mods desired

  res@tiMainString      = title                 ; add title

  res@xyMarkLineModes   = "Markers"             ; choose which have markers

  res@xyMarkers         =  16                   ; choose type of marker  

  res@xyMarkerColor     = "red"                 ; Marker color

  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

  res@tmLabelAutoStride = True                  ; nice tick mark labels

  plot  = gsn_csm_xy (wks,id933,npp933,res)     ; create plot

  frame(wks)

  clear (wks)

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

  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  system("rm "+plot_name+"-*."+plot_type_new)

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

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

 ;(A)-4 scatter plot, model 933

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

  plot_name = "scatter_model_933"

  title     = "Model "+ model_name +" 933 sites"

  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

  res                   = True                  ; plot mods desired

  res@tiMainString      = title                 ; add title

  res@xyMarkLineModes   = "Markers"             ; choose which have markers

  res@xyMarkers         =  16                   ; choose type of marker  

  res@xyMarkerColor     = "red"                 ; Marker color

  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

  res@tmLabelAutoStride = True                  ; nice tick mark labels

  plot  = gsn_csm_xy (wks,id933,nppmod933,res)    ; create plot

  frame(wks)

  clear (wks)

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

  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  system("rm "+plot_name+"-*."+plot_type_new)

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

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

 ;(A)-5 scatter plot, model vs ob 81

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

  plot_name = "scatter_model_vs_ob_81"

  title     = model_name +" vs ob 81 sites"

  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

  res                   = True                  ; plot mods desired

  res@tiMainString      = title                 ; add title

  res@xyMarkLineModes   = "Markers"             ; choose which have markers

  res@xyMarkers         =  16                   ; choose type of marker  

  res@xyMarkerColor     = "red"                 ; Marker color

  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

  res@tmLabelAutoStride = True                  ; nice tick mark labels

  res@gsnDraw           = False

  res@gsnFrame          = False                 ; don't advance frame yet

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

 ;compute correlation coef. and M

  ccr81 = esccr(nppmod81,npp81,0)

 ;print (ccr81)

 ;new eq

  bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81))) 

  M81  = (1. - (bias/dimsizes(y81)))*5. 

  print (M81)

  tRes  = True

  tRes@txFontHeightF = 0.025

  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"

  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)

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

  plot  = gsn_csm_xy (wks,npp81,nppmod81,res)       ; create plot

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

 ; add polyline

  dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)

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

  draw (plot)

  frame(wks)

  clear (wks)

  delete (dum)

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

  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  system("rm "+plot_name+"-*."+plot_type_new)

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

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

 ;(A)-6 scatter plot, model vs ob 933

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

  plot_name = "scatter_model_vs_ob_933"

  title     = model_name +" vs ob 933 sites"

  wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

  res                   = True                  ; plot mods desired

  res@tiMainString      = title                 ; add title

  res@xyMarkLineModes   = "Markers"             ; choose which have markers

  res@xyMarkers         =  16                   ; choose type of marker  

  res@xyMarkerColor     = "red"                 ; Marker color

  res@xyMarkerSizeF     = 0.01                  ; Marker size (default 0.01)

  res@tmLabelAutoStride = True                  ; nice tick mark labels

  res@gsnDraw           = False

  res@gsnFrame          = False                 ; don't advance frame yet

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

 ;compute correlation coef. and M

  ccr933 = esccr(nppmod933,npp933,0)

 ;print (ccr933)

 ;new eq

  bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))

  M933 = (1. - (bias/dimsizes(y933)))*5.

  print (M933)

  tRes  = True

  tRes@txFontHeightF = 0.025

  correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"

  gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)

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

  plot  = gsn_csm_xy (wks,npp933,nppmod933,res)    ; create plot

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

 ; add polyline

  dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)

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

  draw (plot)

  frame(wks)

  clear (wks)

  delete (dum)

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

  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  system("rm "+plot_name+"-*."+plot_type_new)

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

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

 ;(A)-7 table of site81 (1) (the first 41 sites), model vs ob

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

   table_length = 0.95 

 ; table header name

   table_header_name = "Site ID" 

 ; row (not including header row) 

   nrow       = 41

   row_header = new ((/nrow/),string )

   row_header(0:nrow-1) = id81(0:nrow-1) 

 ; Table header

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

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

   y1    = (/0.85 ,0.95/)        ; Start and end Y

   text1 = table_header_name

   res1               = True

   res1@txFontHeightF = 0.015

   res1@gsFillColor   = "CornFlowerBlue"

 ; Column header1 (equally space in x2)

   delete (col_header)

   delete (ncr2)

   delete (text2)

   col_header = (/"Latitude","Longitude"/) 

   ncol = dimsizes(col_header) 

   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

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

   y2    = y1                    ; same as y1

   text2 = col_header

   res2               = True

   res2@txFontHeightF = 0.015

   res2@gsFillColor   = "Gray"

 ; Column header2 (equally space in x2)

   col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)

   col_header2 = (/"ob",model_name \

                  ,"ob",model_name \

                  /)

   ncol1 = dimsizes(col_header1)

   ncol2 = dimsizes(col_header2)

   ncr21  = (/1,ncol1/)            ; 1 rows, 4 columns

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

   yhalf  = (y1(0)+y1(1))*0.5

   y21    = (/yhalf,y1(1)/)        ; top half of y1

   text21 = col_header1

   res21               = True

   res21@txFontHeightF = 0.015

   res21@gsFillColor   = "Gray"

   ncr22  = (/1,ncol2/)            ; 1 rows, 12 columns

   x22    = x21                    ; start from end of x1

   y22    = (/y1(0),yhalf/)        ; bottom half of y1

   text22 = col_header2

   res22               = True

   res22@txFontHeightF = 0.015

   res22@gsFillColor   = "Gray"

 ; Row header (equally space in y2)

   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   x3    = x1                         ; same as x1

   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   text3 = row_header

   res3               = True

   res3@txFontHeightF = 0.010

   res3@gsFillColor   = "Gray"

 ; Main table body-1

   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   x4    = x2                      ; Start and end x

   y4    = y3                      ; Start and end Y

   text4 = new((/nrow,ncol/),string)

   color_fill4           = new((/nrow,ncol/),string)

   color_fill4           = "white"

 ; color_fill4(:,ncol-1) = "grey"

 ; color_fill4(nrow-1,:) = "green"

   res4               = True       ; Set up resource list

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

   res4@txFontHeightF = 0.015

   res4@gsFillColor   = color_fill4

   delete (color_fill4)

 ; Main table body-2

   ncr42  = (/nrow,ncol2/)           ; nrow rows, ncol columns

   x42    = x21                      ; Start and end x

   y42    = y3                       ; Start and end Y

   text42 = new((/nrow,ncol2/),string)

   color_fill42           = new((/nrow,ncol2/),string)

   color_fill42           = "white"

 ; color_fill42(:,ncol-1) = "grey"

 ; color_fill42(nrow-1,:) = "green"

   res42               = True       ; Set up resource list

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

   res42@txFontHeightF = 0.015

   res42@gsFillColor   = color_fill42

   delete (color_fill42)

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

 ; for table value

   do n = 0,nrow-1

      text4(n,0) = sprintf("%5.2f", y81(n))  

      text4(n,1) = sprintf("%5.2f", x81(n))              

   end do

   do n = 0,nrow-1

      text42(n,0) = sprintf("%5.2f", npp81(n))  

      text42(n,1) = sprintf("%5.2f", nppmod81(n))    

      text42(n,2) = sprintf("%5.2f", rain81(n))     

      text42(n,3) = sprintf("%5.2f", rainmod81(n))          

   end do

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

   plot_name = "table_site81_model_vs_ob1"

   wks = gsn_open_wks (plot_type,plot_name)

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

 ; for table title

   gRes               = True

   gRes@txFontHeightF = 0.02

 ; gRes@txAngleF      = 90

   title_text = "Model vs Observation at 81 Sites (1)"

   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

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

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

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

   gsn_table(wks,ncr21,x21,y21,text21,res21)

   gsn_table(wks,ncr22,x22,y22,text22,res22)

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

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

   gsn_table(wks,ncr42,x42,y42,text42,res42) 

   frame(wks)

   clear (wks)

   delete (col_header)

   delete (row_header)

   delete (res1)

   delete (ncr1)

   delete (text1)

   delete (res2)

   delete (ncr2)

   delete (text2)

   delete (res21)

   delete (ncr21)

   delete (text21)

   delete (res22)

   delete (ncr22)

   delete (text22)

   delete (res3)

   delete (ncr3)

   delete (text3)

   delete (res4)

   delete (ncr4)

   delete (text4)

   delete (res42)

   delete (ncr42)

   delete (text42)

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

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

   system("rm "+plot_name+"-1."+plot_type_new)

   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

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

 ;(A)-8 table of site81 (2) (the last 40 sites), model vs ob

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

   table_length = 0.95 

 ; table header name

   table_header_name = "Site ID" 

 ; row (not including header row)

   nrow       = 40

   row_header = new ((/nrow/),string )

   row_header(0:nrow-1) = id81(41:41+nrow-1)   

 ; Table header

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

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

   y1    = (/0.85 ,0.95/)        ; Start and end Y

   text1 = table_header_name

   res1               = True

   res1@txFontHeightF = 0.015

   res1@gsFillColor   = "CornFlowerBlue"

 ; Column header1 (equally space in x2)

   col_header = (/"Latitude","Longitude"/) 

   ncol = dimsizes(col_header) 

   ncr2  = (/1,ncol/)            ; 1 rows, ncol columns

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

   y2    = y1                    ; same as y1

   text2 = col_header

   res2               = True

   res2@txFontHeightF = 0.015

   res2@gsFillColor   = "Gray"

 ; Column header2 (equally space in x2)

 ; col_header1 = (/"NPP (gC/m2/year)","RAIN (m/year)"/)

   col_header2 = (/"ob",model_name \

                  ,"ob",model_name \

                  /)

   ncol1 = dimsizes(col_header1)

   ncol2 = dimsizes(col_header2)

   ncr21  = (/1,ncol1/)            ; 1 rows, 4 columns

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

   yhalf  = (y1(0)+y1(1))*0.5

   y21    = (/yhalf,y1(1)/)        ; top half of y1

   text21 = col_header1

   res21               = True

   res21@txFontHeightF = 0.015

   res21@gsFillColor   = "Gray"

   ncr22  = (/1,ncol2/)            ; 1 rows, 12 columns

   x22    = x21                    ; start from end of x1

   y22    = (/y1(0),yhalf/)        ; bottom half of y1

   text22 = col_header2

   res22               = True

   res22@txFontHeightF = 0.015

   res22@gsFillColor   = "Gray"

 ; Row header (equally space in y2)

   ncr3  = (/nrow,1/)                 ; nrow rows, 1 columns

   x3    = x1                         ; same as x1

   y3    = (/1.0-table_length,y1(0)/) ; end at start of y1

   text3 = row_header

   res3               = True

   res3@txFontHeightF = 0.010

   res3@gsFillColor   = "Gray"

 ; Main table body-1

   ncr4  = (/nrow,ncol/)           ; nrow rows, ncol columns

   x4    = x2                      ; Start and end x

   y4    = y3                      ; Start and end Y

   text4 = new((/nrow,ncol/),string)

   color_fill4           = new((/nrow,ncol/),string)

   color_fill4           = "white"

 ; color_fill4(:,ncol-1) = "grey"

 ; color_fill4(nrow-1,:) = "green"

   res4               = True       ; Set up resource list

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

   res4@txFontHeightF = 0.015

   res4@gsFillColor   = color_fill4

   delete (color_fill4)

 ; Main table body-2

   ncr42  = (/nrow,ncol2/)           ; nrow rows, ncol columns

   x42    = x21                      ; Start and end x

   y42    = y3                       ; Start and end Y

   text42 = new((/nrow,ncol2/),string)

   color_fill42           = new((/nrow,ncol2/),string)

   color_fill42           = "white"

 ; color_fill42(:,ncol-1) = "grey"

 ; color_fill42(nrow-1,:) = "green"

   res42               = True       ; Set up resource list

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

   res42@txFontHeightF = 0.015

   res42@gsFillColor   = color_fill42

   delete (color_fill42)

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

 ; for table value

   do n = 0,nrow-1

      text4(n,0) = sprintf("%5.2f", y81(n+41))  

      text4(n,1) = sprintf("%5.2f", x81(n+41))

   end do

   do n = 0,nrow-1

      text42(n,0) = sprintf("%5.2f", npp81(n+41))  

      text42(n,1) = sprintf("%5.2f", nppmod81(n+41))    

      text42(n,2) = sprintf("%5.2f", rain81(n+41))     

      text42(n,3) = sprintf("%5.2f", rainmod81(n+41))          

   end do

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

   plot_name = "table_site81_model_vs_ob2"

   wks = gsn_open_wks (plot_type,plot_name)

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

 ; for table title

   gRes               = True

   gRes@txFontHeightF = 0.02

 ; gRes@txAngleF      = 90

   title_text = "Model vs Observation at 81 Sites (2)"

   gsn_text_ndc(wks,title_text,0.50,0.975,gRes)

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

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

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

   gsn_table(wks,ncr21,x21,y21,text21,res21)

   gsn_table(wks,ncr22,x22,y22,text22,res22)

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

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

   gsn_table(wks,ncr42,x42,y42,text42,res42) 

   frame(wks)

   clear (wks)

   delete (col_header)

   delete (row_header)

   delete (res1)

   delete (ncr1)

   delete (text1)

   delete (res2)

   delete (ncr2)

   delete (text2)

   delete (res21)

   delete (ncr21)

   delete (text21)

   delete (res22)

   delete (ncr22)

   delete (text22)

   delete (res3)

   delete (ncr3)

   delete (text3)

   delete (res4)

   delete (ncr4)

   delete (text4)

   delete (res42)

   delete (ncr42)

   delete (text42)

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

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

   system("rm "+plot_name+"-1."+plot_type_new)

   system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

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

 ;(B) histogram 81

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

 ;get data

   RAIN1_1D = ndtooned(rain81)

   RAIN2_1D = ndtooned(rainmod81)

   NPP1_1D  = ndtooned(npp81)

   NPP2_1D  = ndtooned(nppmod81)

 ; number of bin

   nx = 8

   xvalues      = new((/2,nx/),float)

   yvalues      = new((/2,nx/),float)

   mn_yvalues   = new((/2,nx/),float)

   mx_yvalues   = new((/2,nx/),float)

   dx4          = new((/1/),float)

   get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \

          ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)

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

 ;compute correlation coeff and M score

  u = yvalues(0,:)

  v = yvalues(1,:)

  good = ind(.not.ismissing(u) .and. .not.ismissing(v))

  uu = u(good)

  vv = v(good)

  ccr81h = esccr(uu,vv,0)

 ;print (ccr81h)

 ;new eq

  bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))

  M81h = (1.- (bias/dimsizes(uu)))*5.

  print (M81h)

  delete (u)

  delete (v)

  delete (uu)

  delete (vv)

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

 ; histogram res

   resh                = True

   resh@gsnMaximize    = True

   resh@gsnDraw        = False

   resh@gsnFrame       = False

   resh@xyMarkLineMode = "Markers"

   resh@xyMarkerSizeF  = 0.014

   resh@xyMarker       = 16

   resh@xyMarkerColors = (/"Brown","Blue"/)

   resh@trYMinF        = min(mn_yvalues) - 10.

   resh@trYMaxF        = max(mx_yvalues) + 10.

   resh@tiYAxisString  = "NPP (g C/m2/year)"

   resh@tiXAxisString  = "Precipitation (m/year)"

   max_bar = new((/2,nx/),graphic)

   min_bar = new((/2,nx/),graphic)

   max_cap = new((/2,nx/),graphic)

   min_cap = new((/2,nx/),graphic)

   lnres = True

   line_colors = (/"brown","blue"/)

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

 ;(B)-1 histogram plot, ob 81 site 

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

   plot_name = "histogram_ob_81"

   title     = "Observed 81 site"

   resh@tiMainString  = title

   wks   = gsn_open_wks (plot_type,plot_name)    

   xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)

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

 ;Attach the vertical bar and the horizontal cap line 

   delete (x1)

   delete (x2)

   delete (y1)

   delete (y2)

   do nd=0,0

     lnres@gsLineColor = line_colors(nd)

     do i=0,nx-1

       if(.not.ismissing(mn_yvalues(nd,i)).and. \

          .not.ismissing(mx_yvalues(nd,i))) then

 ;

 ; Attach the vertical bar, both above and below the marker.

 ;

         x1 = xvalues(nd,i)

         y1 = yvalues(nd,i)

         y2 = mn_yvalues(nd,i)

         plx = (/x1,x1/)

         ply = (/y1,y2/)

         min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

         y2 = mx_yvalues(nd,i)

         plx = (/x1,x1/)

         ply = (/y1,y2/)

         max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

 ;

 ; Attach the horizontal cap line, both above and below the marker.

 ;

         x1 = xvalues(nd,i) - dx4

         x2 = xvalues(nd,i) + dx4

         y1 = mn_yvalues(nd,i)

         plx = (/x1,x2/)

         ply = (/y1,y1/)

         min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

         y1 = mx_yvalues(nd,i)

         plx = (/x1,x2/)

         ply = (/y1,y1/)

         max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)

       end if

     end do

   end do

   draw(xy)

   frame(wks)

   clear (wks)

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

  system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)

  system("rm "+plot_name+"-*."+plot_type_new)

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

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

 ;(B)-2 histogram plot, model vs ob 81 site 

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

   plot_name = "histogram_model_vs_ob_81"

   title     = model_name+ " vs Observed 81 site"

   resh@tiMainString  = title

   wks   = gsn_open_wks (plot_type,plot_name)    ; open workstation

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

 ; Add a boxed legend using the more simple method, which won't have

 ; vertical lines going through the markers.

   resh@pmLegendDisplayMode    = "Always"

 ; resh@pmLegendWidthF         = 0.1

   resh@pmLegendWidthF         = 0.08

   resh@pmLegendHeightF        = 0.05

   resh@pmLegendOrthogonalPosF = -1.17