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"

 
  nlat  = 360
  mlon  = 720   
  lat   = latGlobeFo(nlat,"lat","latitude","degrees_north")
; lat=LAT(106:253)
  lon   = lonGlobeFo(mlon,"lon","longitude","degrees_east")
  lon   = (/lon -180./)
  lon&lon=lon
; print(lat)
; print(lon)
;*****************************************************************
; Variables to hold binned quantities
;*****************************************************************
  gbin  = new ( (/nlat,mlon/), float ) 
  gknt  = new ( (/nlat,mlon/), integer) 

  gbin  = 0.0                      ; initialization
  gknt  = 0  

;*****************************************************************
 diri = "./"
 fili = systemfunc("cd "+diri+" ; ls 2A25.20090113*.HDF")
 nfil = dimsizes( fili )
;*****************************************************************
;  Loop over the files
;*****************************************************************
  tStrt = systemfunc("date")         ; time the loop (wall clock)

  vNam = "rain"

  do nf=0,nfil-1
     print(nf+"   "+fili(nf))
     f       = addfile(diri+fili(nf), "r")
                                        ; read data
     xs     := f->$vNam$                ; (nscan,nray,ncell1)=>(0,1,2)
                                        ;   9239  49   80
     x      := where(xs.ge.toshort(0), xs*0.01, -9999.0)
     x@_FillValue = -9999.0
     printMinMax(x,0)
     lat2d  := ( f->Latitude )          ; (nscan,nray)
     lon2d  := ( f->Longitude)          ; (nscan,nray)
     LAT1D  := ndtooned( conform(x,lat2d, (/0,1/)) )
     LON1D  := ndtooned( conform(x,lon2d, (/0,1/)) )

     bin_sum(gbin,gknt,lon,lat,LON1D,LAT1D,ndtooned(x)) 
  end do

  wallClockElapseTime(tStrt, "Main File Loop", 0)
;*****************************************************************
; User must perform averaging
;*****************************************************************
  gknt     = where(gknt.eq.0 , gknt@_FillValue, gknt)
  gbin     = gbin/gknt*24
  gbin!0   = "lat"
  gbin!1   = "lon"
  gbin&lat =  lat
  gbin&lon =  lon

  copy_VarCoords(gbin, gknt)        ; copy coords

  if (isfilevaratt(f, vNam, "long_name")) then
      gbin@long_name = "BINNED: "+vNam
      gknt@long_name = "BINNED COUNT: "+vNam
  end if

  if (isfilevaratt(f, vNam, "units")) then
      gbin@units     = f->$vNam$@units
  end if
 file_out="out.bin"     ;output file
;print(gbin)
   fbinwrite(file_out,gbin)

  PLOT = True
  if (PLOT) then
      plot    = new ( 2, "graphic")

      wks     = gsn_open_wks("ps" ,"2A")
      gsn_define_colormap(wks, "BlAqGrYeOrReVi200")  
      
      res                     = True          ; Use plot options
      res@gsnDraw             = False
      res@gsnFrame            = False

      res@cnFillOn            = True          ; Turn on color fill
      res@cnFillMode          = "RasterFill"  ; Turn on raster color
      res@cnLinesOn           = False         ; Turn off contourn lines

;---This resource not needed in V6.1.0
      res@gsnSpreadColors     = True          ; use full colormap

;---This resource defaults to True in NCL V6.1.0
      res@lbLabelAutoStride   = True          ; Every other label

;      res@mpCenterLonF        = 180           ; set map center at 180

      res@mpMinLonF = 115
      res@mpMaxLonF = 150
      res@mpMinLatF = 5
      res@mpMaxLatF = 20


      res@lbOrientation       = "vertical"

      plot(0) = gsn_csm_contour_map_ce(wks,gbin, res)
      plot(1) = gsn_csm_contour_map_ce(wks,gknt, res)

;*****************************************************************
; Panel   
;*****************************************************************

      resP                     = True
     ;resP@txString            = "TRMM: 2A25: "+yyyy+"/"+month+"/"+day
      res@gsnMaximize          = True          ; max size 
      resP@gsnPaperOrientation = "portrait"
      gsn_panel(wks,plot,(/2,1/),    resP)
  end if