;;===============================================================
; SMAP values are provided on the global cylindrical EASE-Grid 2.0.
; Each grid cell has a nominal area of approximately 36 x 36 km2 
; regardless of longitude and latitude. Using this projection, all 
; global data arrays have dimensions of 406 rows and 964 columns. 
;;===============================================================
undef("isrectilinear")
function isrectilinear(x:numeric, y:numeric)
;
; Unsupported and not documented on the WWW
;
; Is a grid rectilinear? 
; The 2D handles the case wher 1D arrays were replicated (say) via 'conform'
;
; x - cartesian or longitude: [*] or [*][*]
; y - cartesian or latitude: [*] or [*][*]
;
; For 2D: [*][*]  
; +  the left  dimension is assumed to be 'y/latitude'
; +  the right dimension is assumed to be 'x/longitude'
;
; reclin = isrectilinear(lon2d, lat2d)   ; True or False
;
local xdim, ydim, xrank, yrank, x1d, y1d, xx, yy, nmsg
begin
  xdim = dimsizes(x)
  xrank= dimsizes(xdim)
  ydim = dimsizes(y)
  yrank= dimsizes(ydim)

  if (xrank.ne.yrank) then
      return(False)           ; ; NOT rectilinear
  end if

  if (isatt(x, "_FillValue")) then
      nmsg = num(ismissing(x))
      if (nmsg.ne.0) then
          print("isrectilinear: x has _FillValue: nmsg="+nmsg)
          return(False)       ; no  _FillValue for coordinate
                              ; NOT rectilinear
      end if
  end if

  if (isatt(y, "_FillValue")) then
      nmsg = num(ismissing(y))
      if (nmsg.ne.0) then
          print("isrectilinear: y has _FillValue: nmsg="+nmsg)
          return(False)       ; NOT rectilinear               
                              ; no  _FillValue for coordinate
      end if
  end if

  if (xrank.eq.1 .and.  abs(isMonotonic(x)).eq.1 \
                 .and.  abs(isMonotonic(y)).eq.1 ) then
      return(True)                                 ; rectilinear
  end if

  if (xrank.eq.2) then
      x1d = x(0,:)            
      y1d = y(:,0)
      if (abs(isMonotonic(x1d)).eq.1 .and. abs(isMonotonic(y1d)).eq.1 ) then
          xx = conform(x, x1d, 1)
          yy = conform(y, y1d, 0)
          if (all(xx.eq.x) .and. all(yy.eq.y)) then
              return(True)                         ; rectilinear
          end if
      else
          return(False)                            ; NOT rectilinear
      end if
  end if
end
;-----------------------------------------------------------------------
;                     MAIN
;-----------------------------------------------------------------------

  iDir    = "/Users/shea/Data/HDF5/"
  iFile   = "SMAP_L4_SM_gph_20150815T013000_Vv3030_001.h5"
  iPath   = iDir+iFile
  f       = addfile(iPath, "r")

  lat2d     = f->cell_lat
  printVarSummary(lat2d)
  printMinMax(lat2d,0)
  print("---------------")

  lon2d     = f->cell_lon
  printVarSummary(lon2d)
  printMinMax(lon2d,0)
  print("---------------")
  print("lon2d(0,:)="+sprintf("%9.6f", lon2d(0,:)))
  print("lat2d(:,0)="+sprintf("%9.6f", lat2d(:,0)))

 ;print(sprintf("%9.6f", lon2d)+"  "+sprintf("%9.6f", lat2d))
  print("---------------")

  gridType = isrectilinear(lon2d, lat2d)
  print("gridType="+gridType)