;==============================================================
;               SAPHIR Satellite Data (MT1)
;                  Radiometer Swath(s)
;==============================================================
;                          MAIN
;==============================================================
begin
  vNam1  = "TOA"
  netCDF =  True
  filo   = "2012_10_25_21.BINAVG_POIS.nc"     ; output file name

;==============================================================
  NLAT = 1200
  NLON = 2400 
  lt   = latGlobeFo(NLAT,"lat","latitude","degrees_north")

  latS = -30.0
  latN =  30.0
  lat  = lt({latS:latN})  
  printVarSummary(lat)
  printMinMax(lat,0)
  print("-----")
  
  lon  = lonGlobeFo(NLON,"lon","longitude","degrees_east")
  printVarSummary(lon)
  printMinMax(lon,0)
  print("-----")

;===============================================================
; Read HDF5: Use full path 
;
; NOTE: The attributes Scale_Factor, Add_Offset and Fill_Value
;           are type "string". This is very 'unusual.'
;       _FillValue should match variable type with which  
;           it is associated.
;===============================================================

  diri = "./"
  fili = systemfunc("cd "+diri+"; ls MT1*.h5")
  f    = addfile(fili, "r")

  vnlon = "/ScienceData/Longitude_TOA_Pixels"
  ulon  = f->$vnlon$
  ulon@_FillValue = toushort(ulon@Fill_Value)  
  lon2d = ulon*tofloat(ulon@Scale_Factor) + tofloat(ulon@Add_Offset)

  vnlat = "/ScienceData/Colatitude_TOA_Pixels"
  ulat  = f->$vnlat$
  ulat@_FillValue = toushort(ulat@Fill_Value)
  lat2d = ulat*tofloat(ulat@Scale_Factor) + tofloat(ulat@Add_Offset)
  lat2d = 90.- lat2d        ; colatitude ==> latitude
  
  vnrh1 = "/ScienceData/Flux_Longwave_TOA"
  urh1  = f->$vnrh1$
  urh1@_FillValue = toushort(urh1@Fill_Value)
  data  = urh1*tofloat(urh1@Scale_Factor) + tofloat(urh1@Add_Offset)

;==============================================================
; NCL wants CF/COARDS recognized attributes
;==============================================================

 lon2d@units     = "degree_east"    
 lat2d@units     = "degree_north"
 data@long_name  = urh1@Long_Name
 data@units      = urh1@Units

 printVarSummary(lat2d)     ; [2383] x [51]
 printMinMax(lat2d,0)
 print("-----")

 printVarSummary(lon2d)     ; [2383] x [51]
 printMinMax(lon2d,0)
 print("-----")

 printVarSummary(data)      ; [2383] x [51]
 printMinMax(data,0)
 print("-----")

;==============================================================
; Data distribution
;==============================================================
 opt   = True
 opt@PrintStat = True
 stats = stat_dispersion(data, opt )    ; statistical overview

;==============================================================
; Perform bin averaging
;==============================================================
 nlat  = dimsizes(lat)
 GBIN  = new ( (/nlat,NLON/), float ) 
 GKNT  = new ( (/nlat,NLON/), integer ) 
 
 GBIN  = 0.0                      ; initialization
 GKNT  = 0
 bin_sum(GBIN,GKNT,lon,lat,ndtooned(lon2d),ndtooned(lat2d),ndtooned(data))
 GKNT     = where(GKNT.eq.0 , GKNT@_FillValue, GKNT)
 GBIN     = GBIN/GKNT

 GBIN!0   = "lat"
 GBIN!1   = "lon"
 GBIN&lat =  lat
 GBIN&lon =  lon
 copy_VarCoords(GBIN, GKNT)        ; copy coords

 GKNT@long_name = "BINNED COUNT: "+data@long_name
 GBIN@long_name = "BINNED: "+data@long_name
 GBIN@units     = data@units

 printVarSummary(GBIN)
 printMinMax(GBIN,0)
 nmsg = num(ismissing(GBIN))
 print("GBIN: nmsg="+nmsg)
 print("-----")

 printVarSummary(GKNT)
 printMinMax(GKNT,0)
 nmsg = num(ismissing(GKNT))
 print("GKNT: nmsg="+nmsg)
 print("-----")

;------------------------------------------------------------------
; set the poisson_grid_fill arguments
;------------------------------------------------------------------
 GBIN_FILL = GBIN  
 
 nscan = 2000             ; usually *much* fewer
 eps   = 0.1              ; variable dependent
 gtype = True             ; GBIN is cyclic            
 guess = 0              
 relc  = 0.6              ; standard relaxation coef
 begTime = get_cpu_time() 
 poisson_grid_fill(GBIN_FILL, gtype, guess, nscan, eps, relc, 0)
 print("poisson_grid_fill: " + (get_cpu_time() - begTime))

 printVarSummary(GBIN_FILL)
 printMinMax(GBIN_FILL,0)
 nmsg = num(ismissing(GBIN_FILL))
 print("GBIN_FILL: nmsg="+nmsg)
 print(GBIN(:,1200)+"   "+GBIN_FILL(:,1200))
 print("-----")

;------------------------------------------------------------------
; Manually eliminate extra boundary values
;------------------------------------------------------------------
 wcStrt  = systemfunc("date")
 begTime = get_cpu_time() 
 do nl=0,NLON-1
    j := ind(.not.ismissing(GBIN(:,nl)))
    if (.not.ismissing(j(0))) then
        nj    = dimsizes(j)
        jStrt = j(0)
        jLast = j(nj-1)
        if (jStrt.ne.0) then
            GBIN_FILL(0:jStrt-1,nl) = GBIN_FILL@_FillValue
        end if
        if (jLast.ne.(nlat-1)) then
            GBIN_FILL(jLast+1:,nl)  = GBIN_FILL@_FillValue
        end if
    end if
 end do
 print("GBIN_FILL loop: " + (get_cpu_time() - begTime))
 wallClockElapseTime(wcStrt, "GBIN_FILL loop:", 0)

 printVarSummary(GBIN_FILL)
 printMinMax(GBIN_FILL,0)
 nmsg = num(ismissing(GBIN_FILL))
 print("GBIN_FILL: nmsg="+nmsg)
 print("-----")

;--------------------------------------------------
; associate lat2d/lon2d for plotting swath
;---------------------------------------------------
  data@lat2d = lat2d
  data@lon2d = lon2d
;--------------------------------------------------
; plots
;---------------------------------------------------
  plt_typ = "png"

  plt_dir = "./"
  plt_nam = "SAPHIR_"+vNam1+"_BINAVG_POIS"
  plt_pth = plt_dir+plt_nam

  wks     = gsn_open_wks(plt_typ, plt_pth)   
  
  plot= new(3,"graphic")  ; space for 3 graphical objects
  
  res=True 
  res@gsnDraw    = False  ; do not immediately draw
  res@gsnFrame   = False  ; do not immediately flush
  res@gsnAddCyclic=False  ; satellite data not cyclic grid
  
  res@cnFillOn  = True             ; enable contour fill
  res@cnLinesOn = False            ; turn off contour lines
  res@cnLineLabelsOn = False       ; turn off contour line labels
  res@cnFillMode     ="CellFill"   ; "CellFill" may be faster
  res@cnFillPalette  = "BlAqGrYeOrReVi200"  
  res@lbLabelBarOn   = False       ; turn off individual lb's
  res@mpFillOn       = False       ; no gray for land
  
  res@gsnCenterString= "source: (" + str_join(tostring(dimsizes(data))," x ") + ")"
  plot(0)=gsn_csm_contour_map(wks,data,res)  

  res@gsnAddCyclic  = True         ; cyclic grid
  res@gsnCenterString= "("+ str_join(tostring(dimsizes(GBIN))," x ") + ")"
  plot(1)=gsn_csm_contour_map(wks,GBIN,res)

  res@gsnCenterString= "GBIN_FILL: (" + str_join(tostring(dimsizes(GBIN_FILL))," x ") + ")"
  plot(2)=gsn_csm_contour_map(wks,GBIN_FILL,res)
  
;---Panel
  resP= True                                  ; panel mods
  resP@txString= "SAPHIR: MT: "+vNam1         ; 6.3.0 and earlier
 ;resP@gsnPanelMainString  = SAPHIR: MT: "+vNam1 ; 6.4.0 and later
  resP@gsnMaximize=True    ; make plot large
 ;resP@gsnPaperOrientation = "portrait"       ; force portrait orientation
  resP@gsnPanelLabelBar    = True             ; add common colorbar
  gsn_panel(wks,plot,(/3,1/),resP)          ; now draw as one plot

;*****************************************************************
; netCDF  
;*****************************************************************
  if (netCDF) then
      diro    = "./"                   ; output directory
      pout    = diro+filo
      system("/bin/rm -f "+pout)       ; remove any pre-existing file
      ncdf    = addfile(pout ,"c")     ; open output netCDF file
;================================================================
; create global attributes of the file (not required)
;================================================================
      fAtt               = True        ; assign file attributes
      fAtt@title         = "SAPHIR Relative Humidity Swath Binned + Poisson " 
      fAtt@creation_date = systemfunc ("date")        
      fileattdef( ncdf, fAtt )         ; copy file attributes    
                                       ; recommended (generally) 
      ncdf->$vNam1$ = GBIN 
  end if
end