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

;--- Create desired grid
  latS =    0
  latN =   70
  lonW = -120
  lonE =    0

  dlat =  2.0
  dlon =  2.0

  nlat = toint((latN-latS)/dlat) + 1
  mlon = toint((lonE-lonW)/dlon) + 1

  lat  = fspan(latS, latN, nlat)
  lon  = fspan(lonW, lonE, mlon)

  lat@units = "degrees_north"
  lon@units = "degrees_east"

 ;print(lat)
 ;print(lon)

  count     = new( (/nlat,mlon/), "float", 1e20) 
  count!0   = "lat"
  count!1   = "lon"
  count&lat =  lat
  count&lon =  lon
  printVarSummary(count)

;--- Read data
 ;???
  clat = ...
  clon = ...
  npts = dimsizes(clat)

;--- Bin data

  count = 0

  do n=0,npts-1
     jl = toint((clat(n)-latS)/dlat) 
     il = toint((clon(n)-lonW)/dlon) 
     count(jl,il) = count(jl,il) + 1
  end do
  print("count: min="+min(count)+"   max="+max(count))

  count = where(count.eq.0, count@_FillValue,count)

;************************************************
; create plot
;************************************************
  wks = gsn_open_wks("x11" ,"colinz")                ; open a ps file
  gsn_define_colormap(wks,"BlAqGrYeOrRe")        ; choose colormap

  res                       = True     ; plot mods desired
  res@gsnMaximize           = True
  res@gsnSpreadColors       = True     ; use full range of color map
  res@gsnAddCyclic          = False    

  res@cnFillOn              = True     ; turn on color fill
  res@cnFillMode            = "RasterFill"       ; Raster Mode
  res@cnLinesOn             = False    ; turn of contour lines
  res@cnLevelSpacingF       = 1.0      ; contour spacing
  res@lbLabelAutoStride     = True

  res@mpMinLatF             = latS
  res@mpMaxLatF             = latN
  res@mpMinLonF             = lonW
  res@mpMaxLonF             = lonE
  res@mpCenterLonF          = (lonE+lonW)*0.5
  res@mpGridAndLimbOn       = True  
  res@mpGridLineDashPattern = 2             ; Dashed lines
  res@mpGridLatSpacingF     = 5.0
  res@mpGridLonSpacingF     = 10.0

 ;res@gsnLeftString         = "..."
  res@gsnCenterString       = "Gulf Track Density"
 ;res@gsnRightString        = "..."

  plot = gsn_csm_contour_map_ce(wks,count, res)