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

 ; histogram normalized by rain and compute correleration

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

 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 pminmax(data:numeric,name:string)

 begin

   print ("min/max " + name + " = " + min(data) + "/" + max(data))

   if(isatt(data,"units")) then

     print (name + " units = " + data@units)

   end if

 end

 ; Main code.

 begin

  nclass = 20

  plot_type     = "ps"

  plot_type_new = "png"

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

 ; read in data: model       

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

 ;film  = "b30.061n_1995-2004_MONS_climo_lnd.nc"

 ;model_name = "b30.061n"

 ;model_grid = "T31"

 ;film  = "newcn05_ncep_1i_MONS_climo_lnd.nc"

 ;model_name = "newcn"

 ;model_grid = "1.9"

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

 ;model_name = "06cn"

 ;model_grid = "T42"

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

 ;model_name = "06casa"

 ;model_grid = "T42"

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

 ;model_name = "10cn"

 ;model_grid = "T42"

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

  model_name = "10casa"

  model_grid = "T42"

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

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

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

  laimod  = fm->TLAI

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

 ; read in data: observed

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

  ob_name = "MODIS MOD 15A2 2000-2005"

  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"

  filo1  = "land_class_"+model_grid+".nc"

  filo2  = "LAI_2000-2005_MONS_"+model_grid+".nc"

  fo1 = addfile(diro+filo1,"r")

  fo2 = addfile(diro+filo2,"r")

  classob    = tofloat(fo1->LAND_CLASS)               

  laiob      = fo2->LAI

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

 ; for plotting table

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

 ; table header name

   table_header_name = "LAI" 

 ; column (not including header column)

   col_header1 = (/"Mean","Max","Phase","Growth"/)

   col_header2 = (/"ob","model","M" \

                  ,"ob","model","M" \

                  ,"ob","model","M" \

                  ,"ob","model","M" \

                  /)

   ncol1 = dimsizes(col_header1)

   ncol2 = dimsizes(col_header2)

   ncol  = ncol2 

 ; row (not including header row)

   row_header = (/"Water Bodies" \

                 ,"Evergreen Needleleaf Forests" \

                 ,"Evergreen Broadleaf Forests" \

                 ,"Deciduous Needleleaf Forest" \

                 ,"Deciduous Broadleaf Forests" \

                 ,"Mixed Forests" \                      

                 ,"Closed Bushlands" \                   

                 ,"Open Bushlands" \                     

                 ,"Woody Savannas (S. Hem.)" \           

                 ,"Savannas (S. Hem.)" \                 

                 ,"Grasslands" \                         

                 ,"Permanent Wetlands" \                 

                 ,"Croplands" \                          

                 ,"Urban and Built-Up" \                 

                 ,"Cropland/Natural Vegetation Mosaic" \ 

                 ,"Permanent Snow and Ice" \             

                 ,"Barren or Sparsely Vegetated" \       

                 ,"Unclassified" \                       

                 ,"Woody Savannas (N. Hem.)" \           

                 ,"Savannas (N. Hem.)" \

                 ,"All biome average" \                

                 /)  

   nrow = dimsizes(row_header)                  

 ; arrays to be passed to table. 

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

   table_length = 0.995 

 ; Table header

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

   xx1    = (/0.005,0.25/)        ; Start and end X

   yy1    = (/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)

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

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

   yy21    = (/0.9475,0.995/)        ; half of y1

   text21 = col_header1

   res21               = True

   res21@txFontHeightF = 0.015

   res21@gsFillColor   = "Gray"

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

   xx22    = xx21                    ; start from end of x1

   yy22    = (/0.900,0.9475/)       ; half of y1

   text22 = col_header2

   res22               = True

   res22@txFontHeightF = 0.015

   res22@gsFillColor   = "Gray"

 ; Row header (equally space in y2)

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

   xx3    = xx1                      ; same as x1

   yy3    = (/1.0-table_length,0.900/) ; end at start of y1

   text3 = row_header

   res3               = True

   res3@txFontHeightF = 0.01

   res3@gsFillColor   = "Gray"

 ; Main table body

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

   xx4    = xx21                      ; Start and end x

   yy4    = yy3                      ; Start and end Y

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

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

   color_fill4      = "white"

   color_fill4(nrow-1,:) = "CornFlowerBlue"

   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)

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

 ; plot global land class: observed

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

 ;global res

   resg                     = True             ; Use plot options

   resg@cnFillOn            = True             ; Turn on color fill

   resg@gsnSpreadColors     = True             ; use full colormap

 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

 ; resg@lbLabelAutoStride   = True

   resg@cnLinesOn           = False            ; Turn off contourn lines

   resg@mpFillOn            = False            ; Turn off map fill

   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   resg@cnMinLevelValF       = 1.              ; Min level

   resg@cnMaxLevelValF       = 19.             ; Max level

   resg@cnLevelSpacingF      = 1.              ; interval

 ;global contour ob

   classob@_FillValue = 1.e+36

   classob = where(classob.eq.0,classob@_FillValue,classob)

   plot_name = "global_class_ob"

   title     = ob_name

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   plot = gsn_csm_contour_map_ce(wks,classob,resg)   

   frame(wks)

   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)

   clear (wks)

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

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

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

   nclassn     = nclass + 1

   range       = fspan(0,nclassn-1,nclassn)

 ; print (range)

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

 ; particular range, and then take an average.

   nr           = dimsizes(range)

   nx           = nr-1

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

   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.

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

 ;(A) mean

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

 ; get data

   laiob_mean  = dim_avg_Wrap(laiob(lat|:,lon|:,time|:))

   laimod_mean = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))

   DATA11_1D = ndtooned(classob)

   DATA12_1D = ndtooned(laiob_mean)

   DATA22_1D = ndtooned(laimod_mean)

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

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

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

   do nd=0,1

 ; See if we are doing model or observational data.

     if(nd.eq.0) then

       data_ob  = DATA11_1D

       data_mod = DATA12_1D

     else

       data_ob  = DATA11_1D

       data_mod = DATA22_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((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))

       else

 ;        print("")

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

          idx = ind(data_ob.ge.range(i))

       end if

 ; Calculate average, and get min and max.

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

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

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

         mx_yvalues(nd,i) = max(data_mod(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(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_ob)

     delete(data_mod)

   end do

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

 ; compute correlation coef and M score

   u = yvalues(0,:)

   v = yvalues(1,:)

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

   uu = u(good)

   vv = v(good)

   ccrMean = esccr(uu,vv,0)

 ; print (ccrMean)

 ; new eq

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

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

   print (Mmean)

  do i=0,nrow-2

   text4(i,0) = sprintf("%5.2f",u(i))

   text4(i,1) = sprintf("%5.2f",v(i))

   text4(i,2) = "-"

  end do

   text4(nrow-1,0) = sprintf("%5.2f",avg(u))

   text4(nrow-1,1) = sprintf("%5.2f",avg(v))

   text4(nrow-1,2) = sprintf("%5.2f",Mmean)

  delete (u)

  delete (v)

  delete (uu)

  delete (vv)

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

 ; histogram res

   resm                = True

   resm@gsnMaximize    = True

   resm@gsnDraw        = False

   resm@gsnFrame       = False

   resm@xyMarkLineMode = "Markers"

   resm@xyMarkerSizeF  = 0.014

   resm@xyMarker       = 16

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

 ; resm@trYMinF        = min(mn_yvalues) - 10.

 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   resm@trYMinF        = min(mn_yvalues) - 2

   resm@trYMaxF        = max(mx_yvalues) + 4

   resm@tiYAxisString  = "Mean LAI (Leaf Area Index)"

   resm@tiXAxisString  = "Land Cover Type"

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

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

 ; Start the graphics.

   plot_name = "histogram_mean"

   title     = model_name + " vs Observed"

   resm@tiMainString  = title

   wks   = gsn_open_wks (plot_type,plot_name)

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

 ; Add a boxed legend using the more simple method

   resm@pmLegendDisplayMode    = "Always"

 ; resm@pmLegendWidthF         = 0.1

   resm@pmLegendWidthF         = 0.08

   resm@pmLegendHeightF        = 0.05

   resm@pmLegendOrthogonalPosF = -1.17

 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

 ; resm@pmLegendParallelPosF   =  0.18

   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

 ; resm@lgPerimOn              = False

   resm@lgLabelFontHeightF     = 0.015

   resm@xyExplicitLegendLabels = (/"observed",model_name/)

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

   tRes  = True

   tRes@txFontHeightF = 0.025

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

   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

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

 ;Attach the vertical bar and the horizontal cap line 

   do nd=0,1

     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)

         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

         y2 = mx_yvalues(nd,i)

         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),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)

         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

         y1 = mx_yvalues(nd,i)

         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

       end if

     end do

   end do

   draw(xy)

   frame(wks)

   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)

   clear (wks)

  delete (DATA11_1D)

  delete (DATA12_1D)

  delete (DATA22_1D)

 ;delete (range)

 ;delete (xvalues) 

  delete (yvalues)

  delete (mn_yvalues)

  delete (mx_yvalues)

  delete (good)

  delete (max_bar)

  delete (min_bar)

  delete (max_cap)

  delete (min_cap)

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

 ;global res

   resg                     = True             ; Use plot options

   resg@cnFillOn            = True             ; Turn on color fill

   resg@gsnSpreadColors     = True             ; use full colormap

 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

 ; resg@lbLabelAutoStride   = True

   resg@cnLinesOn           = False            ; Turn off contourn lines

   resg@mpFillOn            = False            ; Turn off map fill

   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   resg@cnMinLevelValF       = 0.              ; Min level

   resg@cnMaxLevelValF       = 10.             ; Max level

   resg@cnLevelSpacingF      = 1.              ; interval

 ;global contour ob

   delta = 0.00001

   laiob_mean = where(ismissing(laiob_mean).and.(ismissing(laimod_mean).or.(laimod_mean.lt.delta)),0.,laiob_mean)

   plot_name = "global_mean_ob"

   title     = ob_name

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg)   

   frame(wks)

   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)

   clear (wks)

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

 ;global contour model

   plot_name = "global_mean_model"

   title     = "Model " + model_name

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg)   

   frame(wks)

   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)

   clear (wks)

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

 ;(B) max

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

 ; get data

 ; observed

   laiob_max = laiob(0,:,:)

   s         = laiob(:,0,0)

   laiob_max@long_name = "Leaf Area Index Max"

   dsizes_z = dimsizes(laiob)

   nlat     = dsizes_z(1)

   nlon     = dsizes_z(2)

   do j = 0,nlat-1

   do i = 0,nlon-1

      s = laiob(:,j,i) 

      laiob_max(j,i) = max(s)

   end do

   end do

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

   delete (s)

   delete (dsizes_z)          

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

 ; model

   laimod_max = laimod(0,:,:)

   s          = laimod(:,0,0)

   laimod_max@long_name = "Leaf Area Index Max"

   dsizes_z = dimsizes(laimod)

   nlat     = dsizes_z(1)

   nlon     = dsizes_z(2)

   do j = 0,nlat-1

   do i = 0,nlon-1

      s = laimod(:,j,i) 

      laimod_max(j,i) = max(s)

   end do

   end do

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

   delete (s)

   delete (dsizes_z)          

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

   DATA11_1D = ndtooned(classob)

   DATA12_1D = ndtooned(laiob_max)

   DATA22_1D = ndtooned(laimod_max)

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

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

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

   do nd=0,1

 ; See if we are doing model or observational data.

     if(nd.eq.0) then

       data_ob  = DATA11_1D

       data_mod = DATA12_1D

     else

       data_ob  = DATA11_1D

       data_mod = DATA22_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_ob).and.(data_ob.lt.range(i+1)))

       else

 ;        print("")

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

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

       end if

 ; Calculate average, and get min and max.

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

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

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

         mx_yvalues(nd,i) = max(data_mod(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(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_ob)

     delete(data_mod)

   end do

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

 ; compute correlation coef and M score

   u = yvalues(0,:)

   v = yvalues(1,:)

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

   uu = u(good)

   vv = v(good)

   ccrMax = esccr(uu,vv,0)

 ; print (ccrMax)

 ; new eq

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

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

   print (Mmax)

  do i=0,nrow-2

   text4(i,3) = sprintf("%5.2f",u(i))

   text4(i,4) = sprintf("%5.2f",v(i))

   text4(i,5) = "-"

  end do

   text4(nrow-1,3) = sprintf("%5.2f",avg(u))

   text4(nrow-1,4) = sprintf("%5.2f",avg(v))

   text4(nrow-1,5) = sprintf("%5.2f",Mmax)

  delete (u)

  delete (v)

  delete (uu)

  delete (vv)

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

 ; histogram res

   resm                = True

   resm@gsnMaximize    = True

   resm@gsnDraw        = False

   resm@gsnFrame       = False

   resm@xyMarkLineMode = "Markers"

   resm@xyMarkerSizeF  = 0.014

   resm@xyMarker       = 16

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

 ; resm@trYMinF        = min(mn_yvalues) - 10.

 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   resm@trYMinF        = min(mn_yvalues) - 2

   resm@trYMaxF        = max(mx_yvalues) + 4

   resm@tiYAxisString  = "Max LAI (Leaf Area Index)"

   resm@tiXAxisString  = "Land Cover Type"

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

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

 ; Start the graphics.

   plot_name = "histogram_max"

   title     = model_name + " vs Observed"

   resm@tiMainString  = title

   wks   = gsn_open_wks (plot_type,plot_name)

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

 ; Add a boxed legend using the more simple method

   resm@pmLegendDisplayMode    = "Always"

 ; resm@pmLegendWidthF         = 0.1

   resm@pmLegendWidthF         = 0.08

   resm@pmLegendHeightF        = 0.05

   resm@pmLegendOrthogonalPosF = -1.17

 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

 ; resm@pmLegendParallelPosF   =  0.18

   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

 ; resm@lgPerimOn              = False

   resm@lgLabelFontHeightF     = 0.015

   resm@xyExplicitLegendLabels = (/"observed",model_name/)

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

   tRes  = True

   tRes@txFontHeightF = 0.025

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

   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

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

 ;Attach the vertical bar and the horizontal cap line 

   do nd=0,1

     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)

         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

         y2 = mx_yvalues(nd,i)

         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),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)

         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

         y1 = mx_yvalues(nd,i)

         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

       end if

     end do

   end do

   draw(xy)

   frame(wks)

   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)

   clear (wks)

  delete (DATA11_1D)

  delete (DATA12_1D)

  delete (DATA22_1D)

 ;delete (range)

 ;delete (xvalues) 

  delete (yvalues)

  delete (mn_yvalues)

  delete (mx_yvalues)

  delete (good)

  delete (max_bar)

  delete (min_bar)

  delete (max_cap)

  delete (min_cap)

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

 ;global res

   resg                     = True             ; Use plot options

   resg@cnFillOn            = True             ; Turn on color fill

   resg@gsnSpreadColors     = True             ; use full colormap

 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

 ; resg@lbLabelAutoStride   = True

   resg@cnLinesOn           = False            ; Turn off contourn lines

   resg@mpFillOn            = False            ; Turn off map fill

   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   resg@cnMinLevelValF       = 0.              ; Min level

   resg@cnMaxLevelValF       = 10.             ; Max level

   resg@cnLevelSpacingF      = 1.              ; interval

 ;global contour ob

   delta = 0.00001

   laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max)

   plot_name = "global_max_ob"

   title     = ob_name

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg)   

   frame(wks)

   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)

   clear (wks)

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

 ;global contour model

   plot_name = "global_max_model"

   title     = "Model " + model_name

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg)   

   frame(wks)

   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)

   clear (wks)

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

 ;(C) phase

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

 ; get data

 ; observed

   laiob_phase = laiob(0,:,:)

   s           = laiob(:,0,0)

   laiob_phase@long_name = "Leaf Area Index Max Month"

   dsizes_z = dimsizes(laiob)

   nlat     = dsizes_z(1)

   nlon     = dsizes_z(2)

   do j = 0,nlat-1

   do i = 0,nlon-1

      s = laiob(:,j,i) 

      laiob_phase(j,i) = maxind(s) + 1

   end do

   end do

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

   delete (s)

   delete (dsizes_z)          

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

 ; model

   laimod_phase = laimod(0,:,:)

   s            = laimod(:,0,0)

   laimod_phase@long_name = "Leaf Area Index Max Month"

   dsizes_z = dimsizes(laimod)

   nlat     = dsizes_z(1)

   nlon     = dsizes_z(2)

   do j = 0,nlat-1

   do i = 0,nlon-1

      s = laimod(:,j,i) 

      laimod_phase(j,i) = maxind(s) + 1

   end do

   end do

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

   delete (s)

   delete (dsizes_z)          

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

   DATA11_1D = ndtooned(classob)

   DATA12_1D = ndtooned(laiob_phase)

   DATA22_1D = ndtooned(laimod_phase)

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

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

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

   do nd=0,1

 ; See if we are doing model or observational data.

     if(nd.eq.0) then

       data_ob  = DATA11_1D

       data_mod = DATA12_1D

     else

       data_ob  = DATA11_1D

       data_mod = DATA22_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_ob).and.(data_ob.lt.range(i+1)))

       else

 ;        print("")

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

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

       end if

 ; Calculate average, and get min and max.

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

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

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

         mx_yvalues(nd,i) = max(data_mod(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(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_ob)

     delete(data_mod)

   end do

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

 ; compute correlation coef and M score

   u = yvalues(0,:)

   v = yvalues(1,:)

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

   uu = u(good)

   vv = v(good)

   ccrPhase = esccr(uu,vv,0)

 ; print (ccrPhase)

 ; old eq

 ; bias   = abs(avg(vv)-avg(uu))

 ; new eq

   bias   = avg(abs(vv-uu))

   bias   = where((bias.gt. 6.),12.-bias,bias)

   Mphase = ((6. - bias)/6.)*5.

   print (Mphase)

  do i=0,nrow-2

   text4(i,6) = sprintf("%5.2f",u(i))

   text4(i,7) = sprintf("%5.2f",v(i))

   text4(i,8) = "-"

  end do

   text4(nrow-1,6) = sprintf("%5.2f",avg(u))

   text4(nrow-1,7) = sprintf("%5.2f",avg(v))

   text4(nrow-1,8) = sprintf("%5.2f",Mphase)

  delete (u)

  delete (v)

  delete (uu)

  delete (vv)

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

 ; histogram res

   resm                = True

   resm@gsnMaximize    = True

   resm@gsnDraw        = False

   resm@gsnFrame       = False

   resm@xyMarkLineMode = "Markers"

   resm@xyMarkerSizeF  = 0.014

   resm@xyMarker       = 16

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

 ; resm@trYMinF        = min(mn_yvalues) - 10.

 ; resm@trYMaxF        = max(mx_yvalues) + 10.

   resm@trYMinF        = min(mn_yvalues) - 2

   resm@trYMaxF        = max(mx_yvalues) + 4

   resm@tiYAxisString  = "Max LAI (Leaf Area Index) Month"

   resm@tiXAxisString  = "Land Cover Type"

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

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

 ; Start the graphics.

   plot_name = "histogram_phase"

   title     = model_name + " vs Observed"

   resm@tiMainString  = title

   wks   = gsn_open_wks (plot_type,plot_name)

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

 ; Add a boxed legend using the more simple method

   resm@pmLegendDisplayMode    = "Always"

 ; resm@pmLegendWidthF         = 0.1

   resm@pmLegendWidthF         = 0.08

   resm@pmLegendHeightF        = 0.05

   resm@pmLegendOrthogonalPosF = -1.17

 ; resm@pmLegendOrthogonalPosF = -1.00  ;(downward)

 ; resm@pmLegendParallelPosF   =  0.18

   resm@pmLegendParallelPosF   =  0.88  ;(rightward)

 ; resm@lgPerimOn              = False

   resm@lgLabelFontHeightF     = 0.015

   resm@xyExplicitLegendLabels = (/"observed",model_name/)

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

   tRes  = True

   tRes@txFontHeightF = 0.025

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

   gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)

   xy = gsn_csm_xy(wks,xvalues,yvalues,resm)

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

 ;Attach the vertical bar and the horizontal cap line 

   do nd=0,1

     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)

         min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)

         y2 = mx_yvalues(nd,i)

         max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),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)

         min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

         y1 = mx_yvalues(nd,i)

         max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)

       end if

     end do

   end do

   draw(xy)

   frame(wks)

   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)

   clear (wks)

  delete (DATA11_1D)

  delete (DATA12_1D)

  delete (DATA22_1D)

 ;delete (range)

 ;delete (xvalues) 

  delete (yvalues)

  delete (mn_yvalues)

  delete (mx_yvalues)

  delete (good)

  delete (max_bar)

  delete (min_bar)

  delete (max_cap)

  delete (min_cap)

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

 ;global res

   resg                     = True             ; Use plot options

   resg@cnFillOn            = True             ; Turn on color fill

   resg@gsnSpreadColors     = True             ; use full colormap

 ; resg@cnFillMode          = "RasterFill"     ; Turn on raster color

 ; resg@lbLabelAutoStride   = True

   resg@cnLinesOn           = False            ; Turn off contourn lines

   resg@mpFillOn            = False            ; Turn off map fill

   resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

   resg@cnMinLevelValF       = 1.              ; Min level

   resg@cnMaxLevelValF       = 12.             ; Max level

   resg@cnLevelSpacingF      = 1.              ; interval

 ;global contour ob

   delta = 0.00001

   laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase)

   plot_name = "global_phase_ob"

   title     = ob_name

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg)   

   frame(wks)

   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)

   clear (wks)

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

 ;global contour model

   plot_name = "global_phase_model"

   title     = "Model " + model_name

   resg@tiMainString  = title

   wks = gsn_open_wks (plot_type,plot_name)

   gsn_define_colormap(wks,"gui_default")     ; choose colormap

   delete (plot)

   plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg)   

   frame(wks)

   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)

   clear (wks)

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

 ;(D) grow day

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

 ; get data

   day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)

 ; observed

   laiob_grow = laiob(0,:,:)

   laiob_grow@long_name = "Days of Growing Season"

   dsizes_z = dimsizes(laiob)

   ntime    = dsizes_z(0)

   nlat     = dsizes_z(1)

   nlon     = dsizes_z(2)

   do j = 0,nlat-1

   do i = 0,nlon-1

      nday = 0.

      do k = 0,ntime-1

         if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then

            nday = nday + day_of_data(k)

         end if

      end do

      laiob_grow(j,i) = nday

   end do

   end do