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

; combine  scatter, histogram, global and zonal plots

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

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

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

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

 plot_type     = "ps"

 plot_type_new = "png"

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

; read in ob data

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

;(1)

 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)

 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)

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

 filglobe = "Npp_T31_mean.nc"

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

 nppglobe0 = fglobe->NPP

 nppglobe  = nppglobe0

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

; read in model data

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

 model_name = "b30.061n"

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

 film = "b30.061n_1995-2004_ANN_climo_lnd.nc"

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

 nppmod0  = fm->NPP

 rainmod0 = fm->RAIN

 xm       = fm->lon  

 ym       = fm->lat

 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)

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

; 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) plot scatter ob 81

 plot_name = "scatter_ob_81"

 title     = "Observed 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

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

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

 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) plot scatter model 81

 plot_name = "scatter_model_81"

 title     = "Model "+ model_name +" 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

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

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

 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)

;(E) scatter 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

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

;compute correlation coef. and M

 ccr81 = esccr(nppmod81,npp81,0)

 print (ccr81)

;new eq

 bias = sum(abs(nppmod81-npp81)/(nppmod81+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

 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)

;(F) scatter 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

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

;compute correlation coef. and M

 ccr933 = esccr(nppmod933,npp933,0)

 print (ccr933)

;new eq

 bias = sum(abs(nppmod933-npp933)/(nppmod933+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

 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)

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

;(G) histogram 81

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

;get data

  RAIN1_1D = ndtooned(rain81)

  RAIN2_1D = ndtooned(rainmod81)

  NPP1_1D  = ndtooned(npp81)

  NPP2_1D  = ndtooned(nppmod81)

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

  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

  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

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

;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)/(vv+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"/)

;=================================================================

; histogram ob 81 site only

;

  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 

  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)

        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)

;===========================================================================

; histogram model vs ob 81 site 

  plot_name = "histogram_mod-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

; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)

; resh@pmLegendParallelPosF   =  0.18

  resh@pmLegendParallelPosF   =  0.88  ;(rightward)

; resh@lgPerimOn              = False

  resh@lgLabelFontHeightF     = 0.015

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

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

  tRes  = True

  tRes@txFontHeightF = 0.025

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

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

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

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

;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 (RAIN1_1D)

 delete (RAIN2_1D)

 delete (NPP1_1D)

 delete (NPP2_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)   

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

;(H) histogram 933

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

;get data

  RAIN1_1D = ndtooned(rain933)

  RAIN2_1D = ndtooned(rainmod933)

  NPP1_1D  = ndtooned(npp933)

  NPP2_1D  = ndtooned(nppmod933)

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

  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

  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

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

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

 ccr933h = esccr(uu,vv,0)

 print (ccr933h)

;new eq

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

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

 print (M933h)

 delete (u)

 delete (v)

 delete (uu)

 delete (vv)

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

; histogram res

  delete (resh)

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

;=================================================================

; histogram ob 933 site only

  plot_name = "histogram_ob_933"

  title     = "Observed 933 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 

  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)

        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)

  delete (xy)

  clear (wks)

;===========================================================================

; histogram model vs ob 933 site 

  plot_name = "histogram_mod-ob_933"

  title     = model_name+ " vs Observed 933 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

; resh@pmLegendOrthogonalPosF = -1.00  ;(downward)

; resh@pmLegendParallelPosF   =  0.18

  resh@pmLegendParallelPosF   =  0.88  ;(rightward)

; resh@lgPerimOn              = False

  resh@lgLabelFontHeightF     = 0.015

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

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

  tRes  = True

  tRes@txFontHeightF = 0.025

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

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

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

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

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

  delete (xy)

  clear (wks)

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

; global contour 

;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@gsnSpreadColors      = True            ; use full colormap

  resg@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

  resg@cnMinLevelValF       = 0.              ; Min level

  resg@cnMaxLevelValF       = 2200.           ; Max level

  resg@cnLevelSpacingF      = 200.            ; interval

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

;global contour ob

  delta = 0.00000000001

  nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)

  plot_name = "global_ob"

  title     = "Observed MODIS MOD 17"

  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,nppglobe,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_model"

  title     = "Model "+ model_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,nppmod,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 vs ob

  plot_name = "global_model_vs_ob"

  wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

  gsn_define_colormap(wks,"gui_default")     ; choose colormap

  delete (plot)

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

  resg@gsnFrame             = False          ; Do not draw plot 

  resg@gsnDraw              = False          ; Do not advance frame

;(d) compute correlation coef and M score

  uu1 = ndtooned(nppmod)

  vv1 = ndtooned(nppglobe)

  delete (good) 

  good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))

  ug = uu1(good)

  vg = vv1(good)

  ccrG = esccr(ug,vg,0)

  print (ccrG)

  MG   = (ccrG*ccrG)* 5.0

  print (MG)

; plot correlation coef

  gRes  = True

  gRes@txFontHeightF = 0.02

  gRes@txAngleF = 90

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

  gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)

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

;(a) ob

  title     = "Observed MODIS MOD 17"

  resg@tiMainString  = title

  plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)       

;(b) model

  title     = "Model "+ model_name

  resg@tiMainString  = title

  plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg) 

;(c) model-ob

  zz = nppmod

  zz = nppmod - nppglobe

  title = "Model_"+model_name+" - Observed"

  resg@cnMinLevelValF  = -500           ; Min level

  resg@cnMaxLevelValF  =  500.          ; Max level

  resg@cnLevelSpacingF =  50.           ; interval

  resg@tiMainString    = title

  plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 

  pres                            = True        ; panel plot mods desired

  pres@gsnPanelYWhiteSpacePercent = 5           ; increase white space around

                                                ; indiv. plots in panel

  pres@gsnMaximize                = True        ; fill the page

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