load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
  PANEL_ONLY = True

;---Generate some lat/lon arrays
  lat  = ispan(34,52,1)
  lon  = ispan(-84,-61,1)
  nlat = dimsizes(lat)
  nlon = dimsizes(lon)

;---Turn these into coordinate arrays
  lon!0     = "lon"
  lat!0     = "lat"
  lat@units = "degrees_north"
  lon@units = "degrees_east"

;---Create a dummy array. It will be filled with missing values.
  data = new((/nlat,nlon/),float)

;---Fill dummy array, except for last two rows and columns.
  data(:nlat-3,:nlon-3) = generate_2d_array(10, 12, -20., 17., 0, (/nlat-2,nlon-2/))

;---Assign coordinate arrays
  data!0   = "lat"
  data!1   = "lon"
  data&lat = lat
  data&lon = lon

;---Start the plotting section
  wks = gsn_open_wks("png","exterp")

  res                     = True

  res@gsnMaximize         = True        ; Maximize size of plot
  res@gsnSpreadColors     = True        ; Span full color map
  res@cnFillOn            = True        ; Turn on contour fill
  res@cnLinesOn           = False       ; Turn off contour lines
  res@cnLineLabelsOn      = False       ; turn off line labels
  res@gsnAddCyclic        = False       ; Don't try to add cyclic point

  res@mpMinLatF           = min(lat)    ; Zoom in on map area of interest.
  res@mpMaxLatF           = max(lat)  
  res@mpMinLonF           = min(lon)
  res@mpMaxLonF           = max(lon)

  if(PANEL_ONLY) then
    res@lbLabelBarOn = False
    res@gsnDraw      = False
    res@gsnFrame     = False
  end if

;---Create plot of original data
  res@tiMainString      = "Original data with missing values"
  plot_msg = gsn_csm_contour_map(wks,data,res)  ; Create filled contours

;---Fill in the missing values.
  guess     = 1                ; use zonal means
  is_cyclic = False            ; don't try to add cyclic point
  nscan     = 1500             ; usually much less than this
  eps       = 1.e-2            ; variable dependent
  relc      = 0.6              ; relaxation coefficient
  opt       = 0                ; not used
  poisson_grid_fill( data, is_cyclic, guess, nscan, eps, relc, opt)

;---Draw plot with missing values filled in.
  res@tiMainString      = "Original data w/missing values filled in"
  plot_nomsg = gsn_csm_contour_map(wks,data,res)  ; Create filled contours

;---Panel the plots for easy comparison. 
  pres             = True
  pres@gsnMaximize = True
  gsn_panel(wks,(/plot_msg,plot_nomsg/),(/2,1/),pres)

end