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

; 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"

 system("sed s#model_name#"+model_name+"# table.html > table.html.new")

 system("mv -f table.html.new table.html")

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

 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.

  M_lai_mean = sprintf("%.2f", Mmean)

  system("sed s#M_lai_mean#"+M_lai_mean+"# table.html > table.html.new")

  system("mv -f table.html.new table.html")

  print (M_lai_mean)

 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.

  M_lai_max = sprintf("%.2f", Mmax)

  system("sed s#M_lai_max#"+M_lai_max+"# table.html > table.html.new")

  system("mv -f table.html.new table.html")

  print (M_lai_max)

 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.

  M_lai_phase = sprintf("%.2f", Mphase)

  system("sed s#M_lai_phase#"+M_lai_phase+"# table.html > table.html.new")

  system("mv -f table.html.new table.html")

  print (M_lai_phase)

 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)