begin
  stime = get_cpu_time()
  filename = "crain_composite_2015.09.09.00.00_UTC"
  data = asciiread(filename, -1, "float")
  header = readAsciiHead(filename,2)
  time			= str_get_field(header(1),2," ")
  print (time)

; read data
  radar                = data(7::3)
  lat                  = data(8::3)
  lon                  = data(9::3)

  lat@units = "degrees_north"    ; a good idea to set these
  lon@units = "degrees_east"

;---The radar data seems to have two sets of missing values: -999 and -9999
  radar@_FillValue = -9999
  radar@_FillValue = -999

  etime = get_cpu_time()
  print("Reading the data took " + (etime-stime) + " seconds.")

  printMinMax(lat,0)
  printMinMax(lon,0)
  printMinMax(radar,0)
  print (dimsizes(radar))
  print (dimsizes(lat))
  print (dimsizes(lon))

  wks = gsn_open_wks("png","./test_complex")

  nt = 0
  res                      = True
  res@mpMinLatF            = 20.0
  res@mpMaxLatF            = 50.0
  res@mpMinLonF            = 120
  res@mpMaxLonF            = 150
  res@mpFillOn             = True
  res@mpOutlineOn          = True
  res@mpDataBaseVersion    = "MediumRes"
  res@cnFillOn             = True
  res@cnLinesOn            = False
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/0.1, 5, 20, 50, 75, 100, 150, 210/)
  res@cnFillColors         = (/"white","cadetblue","blue3","chartreuse4","gold", \
                               "darkorange","lightpink", "brown2","blueviolet"/)
  res@sfXArray             = lon      ; Associates lat/lon locations
  res@sfYArray             = lat      ; with data values.
; res@gsnAddCyclic         = False    ; Don't add cyclic point

  stime = get_cpu_time()

;---Draw smoothed contours (the default behavior of NCL)
  res@tiMainString = time + " (smooth contours)"
  contour_pre = gsn_csm_contour_map(wks,radar,res)

  etime = get_cpu_time()
  print("smooth contours took " + (etime-stime) + " seconds.")

;---Draw raster contours
  stime = get_cpu_time()

  res@tiMainString = time + " (raster contours)"
  res@cnFillMode = "RasterFill"
  contour_pre = gsn_csm_contour_map(wks,radar,res)

  etime = get_cpu_time()
  print("raster contours took " + (etime-stime) + " seconds.")

;---Draw markers at locations where the data is > 0.1
  mkres = True
  mkres@gsMarkerIndex = 16  ; filled dot
  mkres@gsMarkerColor = "black"
  draw(contour_pre)
  ii = ind(radar.ge.0.1)
  gsn_polymarker(wks,contour_pre,lon(ii),lat(ii),mkres)
  frame(wks)  
end