; 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
;***************************************************************
; These library files are loaded by default in NCL V6.2.0 and newer
;
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" 
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
;**************************************************************
;
;;===============================================================
; 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")
  nFiles  = dimsizes(iFiles)

  pltDir   = "./"
  pltType  = "pdf" ; png
  pltName  = "smap_single_picture"
  pltTitle =  "all to regional"
  pltPath = pltDir+pltName
  wks = gsn_open_wks(pltType,pltPath)

; ------2. read in regional boundary -------------------------------
  minLatBoundary  = 25
  maxLatBoundary  = 35
  minLonBoundary  = 110
  maxLonBoundary  = 120

;------3. read in the reginal boundary in SMAP ------------------
  nFiles  = dimsizes(iFiles)  
  f       = addfile(iFiles(0), "r")
  ;ncl_filedump iFile(0)|less
  lat2d     = f->cell_lat
  printVarSummary(lat2d)
  ;print(lat2d)

  lon2d     = f->cell_lon
  printVarSummary(lon2d)
  ;print(lon2d)

;---Manually add a _FillValue to lat/lon; not sure why 
;---_FillValue not associated with the variable
  lat2d@_FillValue = -9999.0
  lon2d@_FillValue = -9999.0
 ;printMinMax(lat2d, 0)
 ;printMinMax(lon2d, 0)

; ------4.set the cordination -------------------------------

 nSmapLat = dimsizes(lat2d(:,0))
 nSmapLon = dimsizes(lon2d(0,:))
 minSmapLat = min(lat2d)
 maxSmapLat = max(lat2d)
 minSmapLon = min(lon2d)
 maxSmapLon = max(lon2d)

 lat=fspan(minSmapLat, maxSmapLat, nSmapLat)   
    lat!0="lat"
    lat&lat=lat
    lat@units="degrees_north"
    lat@long_name="Latitude"

 lon=fspan(minSmapLon,maxSmapLon, nSmapLon)
    lon!0="lon"
    lon&lon=lon
    lon@units="degrees_east"
    lon@long_name="Longitude"
   
 lat2d!0="lat"
 lat2d!1="lon"
 lat2d&lon=lon
 lat2d&lat=lat

 lon2d!0="lat"
 lon2d!1="lon"
 lon2d&lon=lon
 lon2d&lat=lat

 res                  = True              ; Plot modes desired.
 res@gsnMaximize      = False;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

;---Resources for plotting original (source) data
 res@sfXArray         = lon2d
 res@sfYArray         = lat2d
 res@trGridType       = "TriangularMesh"  

 res@tiMainString     ="word to regional"
 res@gsnLeftString    = "Soil Moisture root zone"    ; default long_name is too long

;do i =0,nFiles-1
  f= addfile(iFiles(0),"r")
  sm   = f->sm_rootzone
  printVarSummary(sm)
  sm!0="lat"
  sm!1="lon"
  sm&lon=lon
  sm&lat=lat
  printVarSummary(sm)
  
  sm_regional = sm({minLatBoundary:maxLatBoundary},{minLonBoundary:maxLonBoundary})
  printVarSummary(sm_regional)
  nSmLat = dimsizes(sm_regional(:,0)) 
  nSmLon = dimsizes(sm_regional(0,:))
  
  ;print(f->time)
  res@gsnRightString = f->time
  plot_smap = gsn_csm_contour_map(wks,sm_regional,res)
;end do