;---File information

 ;dirg  = "/Volumes/MyPassport/harmonie/Ophelia/output/"
  dirg  = "./"                   
  filg  = "fc2017100612+006.grb"
  pthg  = dirg+filg
  fg    = addfile(pthg,"r")
 ;print(fg)    ; same as:   ncl_filedump fc2017100612+006.grb 

;---Extract desired variable

  varName = "SHTFL_GDS3_HTGL_acc6h"
  x     = fg->$varName$
  printVarSummary(x)
  printMinMax(x,0)
  print("-----")

;;   =================> Look at Variable <=================
;;   float SHTFL_GDS3_HTGL_acc6h ( g3_x_0, g3_y_1 )    <=== (439,565)
;;       center :       Dublin
;;       long_name :    Sensible heat flux             <===
;;       units :        W/m^2
;;       _FillValue :   1e+20
;;       coordinates :  g3_lat_0 g3_lon_1              <=== coordinate names
;;       level_indicator :      105
;;       gds_grid_type :        3
;;       parameter_table_version :      253
;;       parameter_number :     122
;;       level :        0
;;       forecast_time :        6
;;       forecast_time_units :  hours
;;       initial_time : 10/06/2017 (12:00)

;---Read appropriate coordinates: Look at Variable

  lat2d = fg->g3_lat_0                      ; <== 2-dimensional: (g3_x_0, g3_y_1)
                                            ; Corners: (24.478, 24.47797, 48.16449, 48.16454)

  lon2d = fg->g3_lon_1                      ; <== 2-dimensional: (g3_x_0, g3_y_1)
                                            ; Corners: (-44.916, -11.0838, -4.569059, -51.43066)

  printVarSummary(lat2d)
  printMinMax(lat2d,0)
  print("-----")
  printVarSummary(lon2d)
  printMinMax(lon2d,0)
  print("-----")

;---Associate the **two-dimensional** lat and lon coordinates
;   Required for plotting only

  x@lat2d = lat2d
  x@lon2d = lon2d

;---Plot

  wks = gsn_open_wks("png" ,"laraQ")      ; send graphics to PNG file

 res                       = True
 res@gsnMaximize           = True
 res@gsnAddCyclic          = False           ; Regional grid

; map plot resources
;res@mpFillOn              = False
 res@mpLimitMode           = "LatLon"
;res@mpDataBaseVersion     = "MediumRes"
;res@mpOutlineBoundarySets = "AllBoundaries"
 space                     = 1.0  ; arbitrary
 res@mpMinLatF             = min(lat2d) - space
 res@mpMaxLatF             = max(lat2d) + space
 res@mpMinLonF             = min(lon2d) - space
 res@mpMaxLonF             = max(lon2d) + space
 
 res@cnFillOn              = True
 res@cnLinesOn             = False
 res@cnFillMode            = "RasterFill"
 res@cnFillPalette         = "amwg"
;res@cnLevelSelectionMode  = "ExplicitLevels"
;res@cnLevels              = (/ 5000., 6000., 7000., 8000., 8500., 9000., \
;                               9500.,10000.,10500.,11000.,11500.,12000., \
;                              12500.,13000.,13500./)

; tickmark resources
 res@pmTickMarkDisplayMode  = "Always"
 res@tmXBLabelFontHeightF   = 0.010

; labelbar resources
;res@pmLabelBarWidthF         = 0.60
;res@txFontHeightF            = 0.012
;res@lbTitleString            = "elevation above mean sea level (feet)"
;res@lbTitleFontHeightF       = 0.012
;res@lbLabelFontHeightF       = 0.008
;res@lbTitleOffsetF           = -0.27
;res@lbBoxMinorExtentF        = 0.15     ; make the labelbar thinner

; title resources
 res@tiMainString      = filg
;res@tiMainOffsetYF    = -0.02    ; Move title down to

 plot = gsn_csm_contour_map(wks,x, res)