; 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

; United States Boundary (rough estimate)                                                         
  latS = 20
  latN = 60
  lonW = -120
  lonE = -75

  ji = region_ind(lat2d, lon2d, latS, latN, lonW, lonE)
  lat_strt = ji(0)
  lat_last = ji(1)
  lon_strt = ji(2)
  lon_last = ji(3)
;
; 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)

  var_subset =   var(lat_strt:lat_last,lon_strt:lon_last)
  lat_subset = lat2d(lat_strt:lat_last,lon_strt:lon_last)
  lon_subset = lon2d(lat_strt:lat_last,lon_strt:lon_last)
  printVarSummary(var_subset)
  printVarSummary(lat_subset)
  printVarSummary(lon_subset)
  printMinMax(var_subset,0)
  printMinMax(lat_subset,0)
  printMinMax(lon_subset,0)

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

;---Open PNG file for graphics
  wks  = gsn_open_wks ("png", "narr_plot_zoom")
    
;---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

  var_subset@lat2d = lat_subset    ; IMPORTANT FOR PLOTTING!
  var_subset@lon2d = lon_subset

;---Set lambert conformal and zoom in on area of interest.
  res@mpProjection                = "LambertConformal" 
  res@mpLimitMode                 = "Corners"                ; choose range of map
  res@mpLeftCornerLatF            = lat2d(lat_strt,lon_strt)
  res@mpLeftCornerLonF            = lon2d(lat_strt,lon_strt)
  res@mpRightCornerLatF           = lat2d(lat_last,lon_last)
  res@mpRightCornerLonF           = lon2d(lat_last,lon_last)

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

  res@mpGeophysicalLineThicknessF = 2.5

  res@tiMainString                = "Zoom in on area of interest"
  plot = gsn_csm_contour_map(wks,var_subset,res)
end