 f1 = addfile("geo_em.d01.nc", "r")
 f2 = addfile("geo_em.d02.nc", "r")
 f3 = addfile("geo_em.d03.nc", "r")

 var1 = wrf_user_getvar(f1,"HGT_M",0)
 lat1 = wrf_user_getvar(f1,"XLAT",0)
 lon1 = wrf_user_getvar(f1,"XLONG",0)
 var2 = wrf_user_getvar(f2,"HGT_M",0)
 lat2 = wrf_user_getvar(f2,"XLAT",0)
 lon2 = wrf_user_getvar(f2,"XLONG",0)
 var3 = wrf_user_getvar(f3,"HGT_M",0)
 lat3 = wrf_user_getvar(f3,"XLAT",0)
 lon3 = wrf_user_getvar(f3,"XLONG",0)

 var1@lat2d = lat1
 var1@lon2d = lon1
 var2@lat2d = lat2
 var2@lon2d = lon2
 var3@lat2d = lat3
 var3@lon2d = lon3

 dom_dims = dimsizes(var1)
 dom_rank = dimsizes(dom_dims)
 nx1 = dom_dims(dom_rank - 1) - 1
 ny1 = dom_dims(dom_rank - 2) - 1
 dom_dims = dimsizes(var2)
 dom_rank = dimsizes(dom_dims)
 nx2 = dom_dims(dom_rank - 1) - 1
 ny2 = dom_dims(dom_rank - 2) - 1  
 dom_dims = dimsizes(var3)
 dom_rank = dimsizes(dom_dims)
 nx3 = dom_dims(dom_rank - 1) - 1
 ny3 = dom_dims(dom_rank - 2) - 1

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

 ; Set some contouring resources.

 cnres                 = True
 cnres@gsnDraw         = False   ; Don't draw plot (will do later)
 cnres@gsnFrame        = False   ; Don't advance framce  (will do later)
 cnres@gsnMaximize     = True    ; Maximize plot in frame
 cnres@cnInfoLabelOn   = False   ; Turn off annotations
 cnres@gsnSpreadColors = True
 cnres@gsnLeftString   = ""
 cnres@gsnCenterString = ""
 cnres@gsnRightString  = ""
 cnres@tiXAxisString   = ""
 cnres@tiYAxisString   = ""

 ; Set contour resources

 cnres@cnLinesOn       = False
 cnres@cnLineLabelsOn  = False
 cnres@cnInfoLabelOn   = False
 cnres@cnFillOn        = True
 cnres@cnFillMode      = "AreaFill"    ; This is the default.
 cnres@lbOrientation   = "Vertical"
;gsn_define_colormap(wks,"MPL_YlOrBr")
 cnres@cnFillPalette = "OceanLakeLandSnow"
 ; Set contour levels.
 cnres@cnLevelSelectionMode = "ExplicitLevels"
 cnres@cnLevels             = (/0, 150, 350, 500, 750, \
                                1000, 1250, 1500, 1750, 2000, 2250, 3000, 4000, 5000, 6000/)
 
 res = cnres   ; Copy contour resources for 1st plot.

 ; Add map resources
  res@mpProjection          = "Mercator"  
  res@mpDataBaseVersion     = "MediumRes"          ; Default is LowRes
  res@mpOutlineDrawOrder    = "PostDraw"           ; Draw map outlines last
  res@mpGridAndLimbOn       = False                ; Turn off lat/lon lines
  res@pmTickMarkDisplayMode = "Always"             ; Turn on map tickmarks
  res@cnSpanFillPalette     = False                ; Use partial color map
  
  res@mpLimitMode        = "Corners"               ; Portion of map to zoom
  res@mpLeftCornerLatF   = lat1(0,0)
  res@mpLeftCornerLonF   = lon1(0,0)
  res@mpRightCornerLatF  = lat1(ny1,nx1)
  res@mpRightCornerLonF  = lon1(ny1,nx1)

  ; Create 1st set of contours over a map
  map = gsn_csm_contour_map(wks, var1, res)
  
  res2              = cnres  ; Copy contour resources.
  res2@lbLabelBarOn = False  ; Labelbar already created in 1st plot  
  res2@gsnMaximize  = False  ; Use maximization from original plot

  ; Create 2nd set of contours 
  plot2 = gsn_csm_contour(wks, var2, res2)

  res3              = cnres  ; Copy contour resources.
  res3@gsnMaximize  = False  ; Use maximization from original plot
  res3@lbLabelBarOn = False  ; Labelbar already created in 1st plot 

  ; Create 3rd set of contours (finest resolution)
  plot3 = gsn_csm_contour(wks, var3, res3) 

  ; Overlay 2nd and 3rd contour plots on original map/contour plot
  overlay(map, plot2)
  overlay(map, plot3)

  ; Draw map and the three contour plots.
  draw(map)

  ; This section is to draw some polylines delineating the
  ; area on the map of the three datasets.

  lnres    = True
  lnres@gsLineThicknessF  = 1.5
  
  xbox = (/lon2(0,0),lon2(0,nx2),lon2(ny2,nx2),lon2(ny2,0),lon2(0,0)/)
  ybox = (/lat2(0,0),lat2(0,nx2),lat2(ny2,nx2),lat2(ny2,0),lat2(0,0)/)

  xbox_out = new(dimsizes(xbox),typeof(xbox))
  ybox_out = new(dimsizes(ybox),typeof(ybox)) 

  datatondc(map, xbox, ybox, xbox_out, ybox_out)
  gsn_polyline_ndc(wks, xbox_out, ybox_out, lnres)

  xbox = (/lon3(0,0),lon3(0,nx3),lon3(ny3,nx3),lon3(ny3,0),lon3(0,0)/)
  ybox = (/lat3(0,0),lat3(0,nx3),lat3(ny3,nx3),lat3(ny3,0),lat3(0,0)/)
  
  datatondc(map, xbox, ybox, xbox_out, ybox_out)
  gsn_polyline_ndc(wks, xbox_out, ybox_out, lnres)
  
  ; Now advance the frame.

   frame(wks)
