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

; histogram normalized by rain and compute correleration

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

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

  data_types    = (/ "Obs",        "Model"       /)

  data_names    = (/ "data.81.nc" , "i01.03cn_1545-1569_ANN_climo.nc" /)

; data_names    = (/ "data.933.nc", "i01.03cn_1545-1569_ANN_climo.nc" /)

; data_names    = (/ "data.81.nc" , "i01.04casa_1605-1629_ANN_climo.nc" /)

; data_names    = (/ "data.933.nc", "i01.04casa_1605-1629_ANN_climo.nc" /)

  filevar_names = (/ (/"PREC_ANN","TNPP_C"/), (/"RAIN","NPP"/) /)

  ndata_types   = dimsizes(data_types)

  data_file_obs = addfile(data_names(0),"r")    ; Open obs file

  data_file_mod = addfile(data_names(1),"r")    ; Open model file

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

; read in data: observed

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

 RAIN1 = tofloat(data_file_obs->PREC_ANN) ; for data.81

;RAIN1 = data_file_obs->PREC              ; for data.933

 NPP1  = data_file_obs->TNPP_C

 xo    = data_file_obs->LONG_DD  

 yo    = data_file_obs->LAT_DD   

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

 print (xo)

 nx = dimsizes(xo)

 do i= 0,nx-1

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

        xo(i) = xo(i)+ 360.

    end if

 end do

 print (xo)

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

; read in data: model       

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

 ai    = data_file_mod->RAIN

 bi    = data_file_mod->NPP

 xi    = data_file_mod->lon      

 yi    = data_file_mod->lat      

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

; interpolate from model grid to observed grid

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

 RAIN2 = linint2_points(xi,yi,ai,True,xo,yo,0) 

 NPP2  = linint2_points(xi,yi,bi,True,xo,yo,0) 

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

; convert unit

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

; Units for these four variables are:

;

; RAIN1 : mm/year

; RAIN2 : mm/s

; NPP1  : g C/m^2/year

; NPP2  : g C/m^2/s

;

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

;

  nsec_per_year = 60*60*24*365                  ; # seconds per year

; Do the necessary conversions.

  RAIN1 = RAIN1 / 1000.

  RAIN2 = (RAIN2/ 1000.) * nsec_per_year

  NPP2     = NPP2 * nsec_per_year

; Redo the units.

  RAIN1@units = "m/yr"

  RAIN2@units = "m/yr"

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

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

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

; print min/max and unit

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

  pminmax(RAIN1,"RAIN1")

  pminmax(RAIN2,"RAIN2")

  pminmax(NPP1,"NPP1")

  pminmax(NPP2,"NPP2")

  RAIN1_1D = ndtooned(RAIN1)

  RAIN2_1D = ndtooned(RAIN2)

  NPP1_1D  = ndtooned(NPP1)

  NPP2_1D  = ndtooned(NPP2)

;

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

; ranges.

;

  nbins       = 15     ; Number of bins to use.

  nicevals    = nice_mnmxintvl(min(RAIN2_1D),max(RAIN2_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(RAIN2_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(RAIN2_1D))

  mn_yvalues   = new((/2,nx/),typeof(RAIN2_1D))

  mx_yvalues   = new((/2,nx/),typeof(RAIN2_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

;

; Start the graphics.

;

  wks = gsn_open_wks("png","xy")

  res             = True

  res@gsnMaximize = True

  res@gsnDraw     = False

  res@gsnFrame    = False

  res@xyMarkLineMode = "Markers"

  res@xyMarkerSizeF   = 0.014

  res@xyMarker       = 16

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

  res@trYMinF        = min(mn_yvalues) - 10.

  res@trYMaxF        = max(mx_yvalues) + 10.

; res@tiMainString  = "Observed vs i01.03cn 81 site"

; res@tiMainString  = "Observed vs i01.03cn 933 site"

  res@tiMainString  = "Observed vs i01.04casa 81 site"

; res@tiMainString  = "Observed vs i01.04casa 933 site"

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

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

;

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

; vertical lines going through the markers.

;

  res@pmLegendDisplayMode    = "Always"

; res@pmLegendWidthF         = 0.1

  res@pmLegendWidthF         = 0.08

  res@pmLegendHeightF        = 0.05

  res@pmLegendOrthogonalPosF = -1.17

; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

; res@pmLegendParallelPosF   =  0.18

  res@pmLegendParallelPosF   =  0.23  ;(rightward)

;  res@lgPerimOn             = False

  res@lgLabelFontHeightF     = 0.015

; res@xyExplicitLegendLabels = (/"observed","model_i01.03cn"/)

  res@xyExplicitLegendLabels = (/"observed","model_i01.04casa"/)

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

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

  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

;

; Here's how to do the legend by hand.

;

; mkres               = True    ; Marker resources

; txres               = True    ; Text resources

; mkres@gsMarkerIndex = 16

; mkres@gsMarkerSizeF = 0.02

; txres@txFontHeightF = 0.02

; txres@txJust        = "CenterLeft"

;

; Change these values if you want to move the marker legend location.

; These values are in the same data space as the plot. 

;

; xlg1_cen = 0.2

; ylg1_cen = 900.

; xlg2_cen = 0.2

; ylg2_cen = 760.

; mkres@gsMarkerColor = "brown"

; lnres@gsLineColor   = "brown"

; lg_mark_legend1 = gsn_add_polymarker(wks,xy,xlg1_cen,ylg1_cen,mkres)

; lg_line_legend1 = gsn_add_polyline(wks,xy,(/xlg1_cen,xlg1_cen/), \

;                                           (/ylg1_cen-60,ylg1_cen+60/),lnres)

; lg_cap_legend11  = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \

;                                           (/ylg1_cen-60,ylg1_cen-60/),lnres)

; lg_cap_legend12  = gsn_add_polyline(wks,xy,(/xlg1_cen-0.1,xlg1_cen+0.1/), \

;                                           (/ylg1_cen+60,ylg1_cen+60/),lnres)

; tx_legend1 = gsn_add_text(wks,xy,"observed",xlg1_cen+0.15,ylg1_cen,txres)

; mkres@gsMarkerColor = "blue"

; lnres@gsLineColor   = "blue"

; lg_mark_legend2 = gsn_add_polymarker(wks,xy,xlg2_cen,ylg2_cen,mkres)

; lg_line_legend2 = gsn_add_polyline(wks,xy,(/xlg2_cen,xlg2_cen/), \

;                                           (/ylg2_cen-60,ylg2_cen+60/),lnres)

; lg_cap_legend21  = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \

;                                           (/ylg2_cen-60,ylg2_cen-60/),lnres)

; lg_cap_legend22  = gsn_add_polyline(wks,xy,(/xlg2_cen-0.1,xlg2_cen+0.1/), \

;                                           (/ylg2_cen+60,ylg2_cen+60/),lnres)

; tx_legend2 = gsn_add_text(wks,xy,"model_i01.03cn",xlg2_cen+0.15,ylg2_cen,txres)

  draw(xy)

  frame(wks)

 u = yvalues(0,:)

 v = yvalues(1,:)

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

 uu = u(good)

 vv = v(good)

 nz = dimsizes(uu)

 print (nz)

 ccr = esccr(uu,vv,0)

 print (ccr)

;old eq

;bias = sum(((vv-uu)/uu)^2)

;M    = (1.- sqrt(bias/nz))*5.

;new eq

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

 M    = (1.- (bias/nz))*5.

 print (bias)

 print (M)

end