;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; sub.bin_sum_other.ncl
; Carl Schreck (cjschrec@ncsu.edu)
; January 2015
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Copied from script by Mary Haley that fixes a bug in bin_sum
; Schreck modified it to make the interface more like the original bin_sum
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;----------------------------------------------------------------------
; This is a different binning function, which
; uses the actual bounds of the lat/lon data for
; the binning.
;----------------------------------------------------------------------
procedure bin_sum_other(gbin, gknt, glon, glat, x, y, z )
local dx, dy, \
      dlat, dlon, glatbnd, blonbnd, k, ilt, iln, delta 
begin

  nx = dimsizes( glon )
  ny = dimsizes( glat )

  dx = (glon(1)-glon(0))/2.
  dy = (glat(1)-glat(0))/2.

  dlat = abs(glat(1)-glat(0))
  dlon = abs(glon(1)-glon(0))

  glatbnd = glat(0) - dlat/2.
  glonbnd = glon(0) - dlon/2.

  delta = (glat(1)-glat(0))/2.

; gknt = 0
; gbin = 0
  do k = 0,dimsizes(z)-1
    ilt = ind((glat-delta).le.y(k).and.y(k).lt.(glat+delta))
    iln = ind((glon-delta).le.x(k).and.x(k).lt.(glon+delta))
    if (.not.ismissing(ilt).and..not.ismissing(iln).and.\
        ilt.ge.0 .and. ilt.lt.ny.and.iln.ge.0 .and.iln.lt.nx) then
;;      print("other bin method: point " + y(k) + "/" + x(k) + " " + \
;;            "got binned into lat: " + (glat(ilt)-delta) + " to " + \
;;                                      (glat(ilt)+delta) + ", lon: " + \
;;                                      (glon(iln)-delta) + " to " + \
;;                                      (glon(iln)+delta) + \
;;                                     " (" + ilt + "," + iln + ")")
      gknt(ilt,iln) = gknt(ilt,iln)+1
      gbin(ilt,iln) = gbin(ilt,iln)+z(k)
    else
;---Make sure we get any data equal to upper bounds of lat/lon
      ilt = ind(y(k).eq.(glat+delta))
      iln = ind(x(k).eq.(glon+delta))
      if (.not.ismissing(ilt).and..not.ismissing(iln).and.\
          ilt.ge.0 .and. ilt.lt.ny.and.iln.ge.0 .and.iln.lt.nx) then
;;        print("other bin method: point " + y(k) + "/" + x(k) + " " + \
;;              "got binned into lat: " + (glat(ilt)-delta) + " to " + \
;;                                        (glat(ilt)+delta) + ", lon: " + \
;;                                        (glon(iln)-delta) + " to " + \
;;                                        (glon(iln)+delta))
        gknt(ilt,iln) = gknt(ilt,iln)+1
        gbin(ilt,iln) = gbin(ilt,iln)+z(k)
      end if
    end if
  end do

end ; bin_sum_other