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"


;********************************************************************
; Utility functions to get (1) date information (2) facilitate netCDF 
;********************************************************************
function parseFileName ( fNam:string )
; parse the file names and extract data information
;           1         2a
; 0123456789012345678901234
; 2A25.20010101.17823.7.HDF
local onNam, n, output, tmp_c
begin
   nNam = dimsizes( fNam )

       output    = new( 4,integer,"No_FillValue")
       tmp_c     = stringtochar(fNam)
       
       output(0) = stringtointeger((/tmp_c(5:8)/))   ; YYYY 
       output(1) = stringtointeger((/tmp_c(9:10)/))   ; MM
       output(2) = stringtointeger((/tmp_c(11:12)/))   ; DD      
   return (output)   
end

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

function create3d (x[*][*]:numeric, time[*]:numeric)
; convert a 2-dimensional to a 3-dimensional variable
; with time as a coordinate variable.
local dimx, ntim, nlat, mlon, x3d
begin
  dimx    = dimsizes(x)

  ntim    = dimsizes(time)
  nlat    = dimx(0)
  mlon    = dimx(1)
          ; perform expansion
  x3d     = conform_dims( (/ntim,nlat,mlon/), x, (/1,2/))
  copy_VarAtts(x, x3d)
  x3d!0   = "time"
  x3d&time=  time
  copy_VarCoords(x, x3d(0,:,:) )
  return(x3d)
end


;**************************
;Main
;**************************

begin
  
  PLOT  = True
  netCDF= True

  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.20010101*.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
;********************************
  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
;*****************************************************************
; time/date
;*****************************************************************
 ydhm = parseFileName(fili(0))
 yyyy    = ydhm(0)
 month   = ydhm(1)
 day     = ydhm(2)
 hh = 0 
 mn = 0
;*****************************************************************
; netCDF  
;*****************************************************************
  if (netCDF) then
      sec     = 0.0
                                       ; COARDS/CF time
      tunits  = "hours since 1990-1-1 00:00:0.0"
      time    = cd_inv_calendar(yyyy,month,day,hh,mn,sec,tunits, 0)
      time!0  = "time"
      time@units = tunits 
                                       ; gregorian date
      date    = yyyy + month + day
      date!0  = "time"
      date@units = "yyyymmdd"

      diro    = "./"                   ; output directory
      filo    = "TRMM2A25Bin_"+yyyy+day+"_bin_sum"    ; output file name
      fout    = diro+filo+".nc"

      system("/bin/rm -f "+fout)       ; remove any pre-existing file
      ncdf    = addfile(fout ,"c")     ; open output netCDF file

    ;================================================================
    ; create global attributes of the file (not required)
    ;================================================================
      fAtt               = True        ; assign file attributes
      fAtt@title         = "2A25 Swath Binned: " + yyyy+day
      fAtt@source_file   = "TRMM2A25."+yyyy+day+".nc"
     ;fAtt@Conventions   = "CF 1.0"   
      fAtt@creation_date = systemfunc ("date")        
      fileattdef( ncdf, fAtt )         ; copy file attributes    

                                       ; recommended (generally) 
      filedimdef(ncdf,"time",-1,True)  ; make time UNLIMITED dimension 
      ncdf->time     = time
      ncdf->date     = date

      ncdf->$vNam$   = create3d( gbin, time)               
      kNam = vNam+"_knt"
      ncdf->$kNam$   = create3d( gknt, time)               
  end if
end