;***************************************************************
; Concepts illustrated: 
;   - Reading a SMAP HDF5 level 4 file with groups
;   - Use 'direct' syntax to access variable within groups
;   - Manually adding _FillValue to latitude and longitude
;   - Plot
;***************************************************************
; 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. 
;***************************************************************

;---1. Read specified h5 files----------------------------------
 ;iFiles   = systemfunc("ls ./*.h5")
  iDir    = "/Users/shea/Data/HDF5/"
  iFiles  = systemfunc("cd "+iDir+"; ls SMAP_L4_SM_gph_20150815T013000_Vv3030_001.h5")
  nFiles   = dimsizes(iFiles)
  print(iFiles)

;---2. Specify plot 'stuff'   ----------------------------------
  pltDir   = "./"
  pltType  = "png"
  pltName  = "smap_single_picture_YiLu"
  pltTitle = "all to regional"
  pltPath  = pltDir+pltName

; ------3. Specify regional boundary -------------------------------
  minLatBoundary  = 25
  maxLatBoundary  = 35
  minLonBoundary  = 110
  maxLonBoundary  = 120

;------4. read in lat/lon in SMAP ------------------
  f       = addfile(iDir+iFiles(0), "r")
  lat2d     = f->cell_lat
  printVarSummary(lat2d)

  lon2d     = f->cell_lon
  printVarSummary(lon2d)

; ------5. set the cordinates -------------------------------
; Generally, the following is *NOT* appropriate. However,
; the grid was tested to be 'rectilinear'. Hence, the following is ok.

  lat   = lat2d(:,0)       ; generally *NOT* appropriate
  lat@long_name="Latitude"
  lat@units="degrees_north"
  lat!0   = "lat"
  lat&lat =  lat

  lon   = lon2d(0,:)       ; generally *NOT* appropriate
  lon@long_name="Longitude"
  lon@units="degrees_east"
  lon!0  ="lon"
  lon&lon= lon
   
; ------6. Plot -------------------------------

  wks = gsn_open_wks(pltType,pltPath)
  res                  = True              ; Plot modes desired.
  res@gsnMaximize      = True              ; Maximize plot
  res@gsnAddCyclic     = False
 
  res@cnFillOn         = True              ; color plot desired
  res@cnLinesOn        = False             ; turn off contour lines
  res@cnLineLabelsOn   = False             ; turn off contour labels
  res@cnFillMode       = "RasterFill"      ; turn raster on      
  res@cnFillPalette    = "BlAqGrYeOrReVi200"
 
  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       = 0.05               ; set min contour level
  res@cnMaxLevelValF       = 0.95               ; set max contour level
  res@cnLevelSpacingF      = 0.05               ; set contour spacing

  extra_space              = 5                  ; arbitrary ... could be 0
  res@mpMinLatF            = minLatBoundary - extra_space 
  res@mpMaxLatF            = maxLatBoundary + extra_space
  res@mpMinLonF            = minLonBoundary - extra_space
  res@mpMaxLonF            = maxLonBoundary + extra_space
;;res@mpFillOn             = False              ; turn off gray continental background

  res@tiMainFontHeightF    = 0.0150

 do i =0,nFiles-1
  f= addfile(iDir+iFiles(i),"r")
  sm   = f->sm_rootzone
  sm@long_name = "Soil Moisture Root Zone"    ; default long_name is too long
  sm!0="lat"
  sm!1="lon"
  sm&lon=lon
  sm&lat=lat
  printVarSummary(sm)
  
  sm_regional = sm({minLatBoundary:maxLatBoundary},{minLonBoundary:maxLonBoundary})
  printVarSummary(sm_regional)

  TIME = cd_calendar( f->time, 0)
  print(TIME)
  ymdh = cd_calendar( f->time,-3)
  res@gsnRightString = ymdh         ; f->time

  res@tiMainString     = iFiles(i)  ; "world to regional"
  plot_smap = gsn_csm_contour_map(wks,sm_regional,res)
 end do