load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
begin

; Read file. 
dir_name  = "/Users/shea/Data/HDF4/"
file_name = "2A25.20090515.65505.7.HDF"
hdf4_file=addfile(dir_name+file_name, "r") 

;--- time information; use the first element for identifier

yyyy     = hdf4_file->Year(0) 
mm       = hdf4_file->Month(0)
dd       = hdf4_file->DayOfMonth(0)
hh       = hdf4_file->Hour(0)

time     = yyyy*1000000 + mm*10000 + dd*100 + hh
time!0   = "time"
time@units = "yyyymmddhh"
time@info  = "Granule First Scan UTC Date"
time&time = time    
         print(time)
         print("===")

; Read lat/lon data.
longitude = hdf4_file->Longitude                   ; [nscan | 9247] x [nray | 49] 
latitude  = hdf4_file->Latitude
;        printMinMax(longitude,True)
;        printMinMax(latitude,True)
;        print("==="))
;        print(dimsizes(longitude))
;        print(dimsizes(latitude))
;        print("==="))

; Correct units to follow the CF conventions.
; In the HDF4 file, the attribute values are "degrees". 

longitude@units = "degrees_east"
latitude@units  = "degrees_north"

; Read data for plot.
   RF_short =  hdf4_file->rain                   ; [nscan | 9247] x [nray | 49] x [ncell1 | 80]
                     printVarSummary(RF_short)
                     printMinMax(RF_short,0)      
                     print("===")
;                    asciiwrite("RF_short.txt",RF_short)

   RF_dim   = dimsizes(RF_short)
   print(RF_dim)
   nscan    = RF_dim(0)                            ; 9247
   nray     = RF_dim(1)                            ; 49
   ncell    = RF_dim(2)                            ; 80

   RF_flt   = RF_short/tofloat(RF_short@scale_factor)   ; not sure about this. non-standard
   RF_flt@_FillValue = -999.0                           ; type float
   RF_flt@units      = RF_short@units
   RF_flt@long_name  = RF_short@hdf_name
   copy_VarCoords(RF_short, RF_flt)
                     printVarSummary(RF_flt)
                     printMinMax(RF_flt, True)  
                     print("===")

   RF_flt   = where (RF_flt.lt.0, RF_flt@_FillValue, RF_flt)
                    printVarSummary(RF_flt)
                    printMinMax(RF_flt,True)
  
  ;asciiwrite("RF_with_missing.txt",RF_with_missing)

;--- create height coordinate + meta data
   
   klev       = ncell   ; 80
   lev        = fspan(20, 0.25, klev)
   lev@long_name = "level"
   lev@units  = "km"
   lev!0      = "lev"
   lev&lev    =  lev 
   printVarSummary(lev)
   print("===")

   RF_flt!2   = "lev"
   RF_flt&lev =  lev 
   printVarSummary(RF_flt)   ; [nscan | 9247] x [nray | 49] x [lev | 80]
   print("===")

   do ii= 0,klev-1   ; ncell-1   ; at each level ('cell')
      print("ii="+sprinti("%0.2i", ii)+"  RF_flt: min="+min(RF_flt(:,:,ii)) \
                                               +" max="+max(RF_flt(:,:,ii)) )
   end do
;                   asciiwrite("RF_flt.txt",RF_flt)

;--- new grid (arbitrary)

   nlat  = 1401                      ; new grid coordinates
   mlon  = 2*nlat-2    
   lat   = latGlobeF(nlat,"lat","latitude","degrees_north")
   lon   = lonGlobeF(mlon,"lon","longitude","degrees_east")

;--- Variables to hold binned quantities

   ntim  = 1   ; will force a time coordinate (not necessary)
   GBIN  = new ( (/ntim,klev,nlat,mlon/), float ) 
   GKNT  = new ( (/ntim,klev,nlat,mlon/), integer ) 
   nt    = ntim-1

;--- Binning

   GBIN  = 0.0                      ; initialization
   GKNT  = 0

   do kl=0,klev-1                   ; each level
       bin_sum(GBIN(nt,kl,:,:),GKNT(nt,kl,:,:),lon,lat   \
              ,ndtooned(longitude), ndtooned(latitude),ndtooned(RF_flt(:,:,kl)) )
   end do

;---Avoid division by 0.0

   GKNT = where(GKNT.eq.0 , GKNT@_FillValue, GKNT)
   GBIN = GBIN/GKNT

;--- Add coordinate information and meta data
 
   GBIN!0    = "time"
   GBIN!1    = "lev"
   GBIN!2    = "lat"
   GBIN!3    = "lon"
   GBIN&time =  time
   GBIN&lev  =  lev
   GBIN&lat  =  lat
   GBIN&lon  =  lon
   GBIN@long_name = RF_flt@long_name
   GBIN@units     = RF_flt@units

   copy_VarCoords(GBIN, GKNT)        ; copy coords

;--- Look at interpolated (binned) data

   printVarSummary(GBIN)
   printMinMax(GBIN,0)
   print("====")
   printVarSummary(GKNT)
   printMinMax(GKNT,0)
   print("====")
   print(GBIN&time)
   print("====")

;--- Use stat_dispersion to examine the data

   opt_stat = True
   opt_stat@PrintStat = True
   stat = stat_dispersion(GBIN, opt_stat )
   print("====")

  LAT_SELECT  = 12.0       ; user specified
  lonL        = 75.0
  lonR        = 80.0 

  XWKS   =   gsn_open_wks("png","Vert_Bin_Plot")
             gsn_define_colormap(XWKS,"amwg256")

  res               =   True ; plot mods desired
  res@cnFillOn      =  True ; enable contour fill
  res@cnFillMode    = "RasterFill"    
  res@cnLinesOn     =  False             ; Turn off contour lines
 ;res@cnLevelSelectionMode =  "ExplicitLevels"     
  res@tiMainString         = file_name
  res@tmYLFormat    = "f"                 ; No unnecessary zeroes
  res@tiYAxisString = "level (km)"        ; if not specified, long_name will be used
  res@tmXBTickSpacingF = 1.0
  res@gsnCenterString  = "lat="+LAT_SELECT

  work = GBIN(nt,:,{LAT_SELECT},{lonL:lonR})        ; explicitly extract for convenience

  thresh = 0.4
  work   = where(work.lt.thresh, work@_FillValue, work)  ; don't plot low values  

  printVarSummary(work)     
  printMinMax(work,0)     
  plot = gsn_csm_contour(XWKS,work, res)

end

;; nLabels                   = 8               ; arbitrary  
;; resx@tmXBLabels           = new(nLabels,"string")
;; resx@tmXBMode             = "Explicit"

;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; Draw Full S-N cross section at a specified longitude grid line
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;;  ml = 200
;;  resx@tmXBValues         := toint( fspan(0,nlat-1,nLabels) )
;;  do i=0,nLabels-1
;;     x = lon2d(resx@tmXBValues(i),ml)
;;     y = lat2d(resx@tmXBValues(i),ml)
;;     resx@tmXBLabels(i) = sprintf("%5.1f", y)+"~C~"+sprintf("%5.1f", x)
;;  end do