f    = addfile("2015-April-02_OPT.2015040100.nc","r")
;;
  vName= "aod"

  x    = f->$vName$
  printVarSummary(x)   ; (time, iy, jx)
  printMinMax(x,0)
  print("-----")

  xlat = f->xlat       ; (iy, jx)   ==> Cartesian
  xlon = f->xlon

  dimx = dimsizes(x)
  ntim = dimx(0)
  nlat = dimx(1)
  mlon = dimx(2)

;---Extract dates: yyyymmddhh

  ymdh = cd_calendar(x&time, -3)
  print(ymdh)

; Explore data: look at distribution

  opt_sd = True         
  opt_sd@PrintStat = True
  stat_x = stat_dispersion(x, opt_sd )
  print("-----")

;---PLOT: Use "Cartesian" coordinates 'iy', 'ix'
;         Use contour levels based upon data distribution. 
;         This may require iteration to find what is best.

  YMDH = 2015040200                  ; desired date
  nt   = ind(YMDH.eq.ymdh)             ; time index corresponding to the requested data

  pltTyp = "png"                       ; png, x11, pdf, ...
  pltDir = "./"
  pltFil = "rcm_"+YMDH
  pltFil = "rcm"
  pltPth = pltDir+pltFil

; Delete cartesian coordinates associated with the variable
; For this application, NCL does not use this information

  delete( [/ x&jx, x&iy /] )

   wks  = gsn_open_wks(pltTyp,"rcm_"+YMDH)
  ;wks  = gsn_open_wks(pltTyp,"rcm")
  
  res                      = True      ; plot mods desired
  res@gsnMaximize          = True      ; maximize plot size
  res@gsnAddCyclic         = False     ; data are regional   
  res@cnFillOn             = False     ; turn on color
  res@cnFillOn             = True      ; turn on color fill
  res@cnLinesOn            = False     ; turn off contour lines
  res@cnLineLabelsOn       = False     ; turn off contour line labels
 ;res@cnFillPalette        = "BlAqGrYeOrRe"
  res@cnFillPalette        = "WhiteBlueGreenYellowRed"
  res@cnFillMode           = "RasterFill"       ; Raster Mode
 ;res@lbOrientation        = "Vertical"
  res@lbBoxLinesOn         = False     ; Turn off labelbar box lines 

  res@mpFillOn             = False
  res@mpLimitMode          = "Corners"          ; choose region of map
  res@mpLeftCornerLatF     = xlat(0,0)
  res@mpLeftCornerLonF     = xlon(0,0)
  res@mpRightCornerLatF    = xlat(nlat-1,mlon-1)
  res@mpRightCornerLonF    = xlon(nlat-1,mlon-1)
  res@mpDataBaseVersion    = "MediumRes"        ; better map outlines
  res@pmTickMarkDisplayMode= "Always"           ; turn on tickmarks
  
; Map information parameters are file attributes [also associated with 'crs' variable

  res@mpProjection          = "LambertConformal"
  res@mpLambertParallel1F   = f@standard_parallel(0)
  res@mpLambertParallel2F   = f@standard_parallel(1)
  res@mpLambertMeridianF    = f@longitude_of_projection_origin

  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/ 0.25, 0.50, 0.75, 1, 1.25, 1.5, 1.75, 2/)               

  res@sfXArray = xlon                              ; equivalently x@lon2d = xlon
  res@sfYArray = xlat                              ; equivalently x@lat2d = xlat

  res@tiMainString         = "ITCP: RCM: "+YMDH

  plot    = gsn_csm_contour_map(wks,x(nt,:,:),res) ; contour the variable

;
; Draw the lat/lon locations as grid lines. Since the
; lat/lon arrays are not associated with "x", we
; have to pass them explicitly via the gsnCoordsLat/Lon
; attributes.
;
  res@gsnDraw = False
  res@gsnFrame= False

  plot    = gsn_csm_contour_map(wks,x(nt,:,:),res) ; contour the variable
  lnres = True
  lnres@gsnCoordsLat = xlat
  lnres@gsnCoordsLon = xlon
  lnres@gsnCoordsAsLines = True
  ;gsn_coordinates(wks,plot,x(nt,:,:),lnres)