; These 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"  
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
;----------------------------------------------------------------------
; Read data
;----------------------------------------------------------------------
  fili    = "narr-a_221_20101201_1200_000.sub.grb"
  f       = addfile (fili, "r")
  var     = f->TMP_221_SFC
  lat2d   = f->gridlat_221                 ; (277 x 349)
  lon2d   = f->gridlon_221

;
; Print some information about the data. Note metadta attached to 
; lat2d / lon2d. This will be used for plotting later.
;
  printVarSummary(var)
  printVarSummary(lat2d)
  printVarSummary(lon2d)
  printMinMax(var,0)
  printMinMax(lat2d,0)
  printMinMax(lon2d,0)

;---Not used in script, but could come in handy
  dimx    = dimsizes( var )
  nlat    = dimx(0)
  mlon    = dimx(1)

;----------------------------------------------------------------------
;  Graphics section
;----------------------------------------------------------------------

;---Open PNG file for graphics
  wks  = gsn_open_wks ("png", "narr_plot")
    
;---Set resources for common for both plots.
  res                        = True

  res@gsnMaximize            = True

  res@cnFillOn               = True                 ; turn color fill  
  res@cnFillPalette          = "BlueWhiteOrangeRed" ; change color map
  res@cnLineLabelsOn         = False                ; turn on/off contour labels
; res@cnLinesOn              = False                ; turn on/off contour lines
; res@cnInfoLabelOn          = False                ; turn on/off contour info label
    
  res@mpGridAndLimbOn        = True
  res@mpGridLineDashPattern  = 2                  ; lat/lon lines as dashed
  res@pmTickMarkDisplayMode  = "Always"           ; turn on nicer map tickmarks
  res@pmTitleZone            = 4                  ; helps move the main title down
  res@gsnAddCyclic           = False              ; regional data, don't add lon cyclic point

;----------------------------------------------------------------------
; Draw plot using lat/lon read off the file. We're using a simple
; cylindrical equidistant projection.
;----------------------------------------------------------------------
  var@lat2d = lat2d     ; For plotting purposes ONLY.
  var@lon2d = lon2d     ; Will remove later.

;---Set resources for cylindrical equidistant projection
  res1                        = res
  res1@mpMinLatF              = -10     ; Zoom in on area
  res1@mpMaxLatF              = 90      ; of interest.
  res1@mpMinLonF              = 140
  res1@mpMaxLonF              = 360
  res1@mpCenterLonF           = 250
  res1@tiMainString           = "Using lat/lon arrays"

;---Create and draw the C.E. plot
  plot = gsn_csm_contour_map(wks,var,res1)

;----------------------------------------------------------------------
; Now try to draw plot using native projection, by using map 
; parameters read off the file. The plot doesn't fit in the exact
; square of the plotting area, but the data looks correctly plotted.
;----------------------------------------------------------------------

;---Do not set these for native projection!
  delete(var@lat2d)
  delete(var@lon2d)

;---Set resources for native projection
  res2 = res

  res2@tfDoNDCOverlay         = True                 ; tell NCL this is a native projetion

  res2@mpProjection           = "LambertConformal" 

  res2@mpLimitMode            = "Corners"            ; choose range of map
  res2@mpLeftCornerLatF       = lat2d@corners(1)     ; lat2d(0,0)
  res2@mpLeftCornerLonF       = lon2d@corners(1)     ; lon2d(0,0)
  res2@mpRightCornerLatF      = lat2d@corners(3)     ; lat2d(nlat-1,mlon-1)
  res2@mpRightCornerLonF      = lon2d@corners(3)     ; lon2d(nlat-1,mlon-1)

  res@mpLambertParallel1F     = lat2d@mpLambertParallel1F
  res@mpLambertParallel2F     = lat2d@mpLambertParallel2F
  res@mpLambertMeridianF      = lat2d@mpLambertMeridianF

  res2@tiMainString           = "Native projection"
  res2@lbOrientation          = "Vertical"           ; just for fun 

;---Create and draw the native plot
  plot = gsn_csm_contour_map(wks,var,res2)

end