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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin

do i=1,1  

dummy = asciiread("all_sum_ens1.asc",-1,"float")
  ncol = 6
  npts = dimsizes(dummy)/ncol

  data = onedtond(dummy,(/npts,ncol/))  ; npts x ncol

  data@_FillValue    = -9999999 
  lon_1              = data(:,1)
  lat                = data(:,0)
  slp                = data(:,2)
  hz                 = data(:,3)
  hz1                = data(:,4)
  hz2                = data(:,5)
  lon_1@_FillValue   = -9999999
  lat@_FillValue     = -99999 
  lon_1@_FillValue   = -99999
  lon = where(lon_1.lt.0 .and. lon_1.ge.(-180),lon_1+360,lon_1)

  printMinMax(lon,True)
  printMinMax(lat,True)
  printMinMax(slp,True)
 
 arr = new((/100,110/),integer,-999)
  arr!0 = "lat"
  arr&lat = fspan(53.82,90,100)
  arr&lat@units = "degrees_north"
  arr!1 = "lon"
  arr&lon = fspan(0,358.9,110)
  arr&lon@units = "degrees_east"  
 
  printVarSummary(arr)
  
  arr = 0    ; set array to all 0's
  do gg = 0,dimsizes(slp)-1
     if (.not.ismissing(lat(gg)).and..not.ismissing(lon(gg))) then   ; check to make sure lat/lons are present
        arr({lat(gg)},{lon(gg)}) = (/ arr({lat(gg)},{lon(gg)})+1 /)
     end if
 end do 
  
  wks = gsn_open_wks("x11","density_sum_ens_"+i)
  gsn_define_colormap(wks,"wh-bl-gr-ye-re")    ; not technically necessary
  
  res = True
           res@gsnPolar   = "NH"                         ; specify the hemisphere 
         res@mpMinLatF  = 53 
   ;     res@mpCenterLatF                      = 90
        res@mpCenterLonF                      = 180 
;  res@mpCenterLonF = (res@mpMinLonF+res@mpMaxLonF)/2.
  res@mpFillOn = False
  res@cnFillOn = True
  res@cnLinesOn = False
  res@gsnSpreadColors = True
  res@cnFillMode = "CellFill"      ; fill each grid box
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels = ispan(1,30,1)
;  res@cnFillColors = (/2,27,37,47,57,67,77,87,97,107,117,127,137,147,157,167,177,187,197,202,207,212,217,225,227,231,237,247,255/)
 printVarSummary(arr) 
   res@tiMainString             = "JAS number of storms, ensemble "+i
  map = gsn_csm_contour_map_polar(wks,arr,res)
  delete(data)
  delete(slp)
  delete(lat)
  delete(lon_1)
  delete(hz)
  delete(hz1)
  delete(hz2)
  delete(dummy)
  delete(lon)
  delete(arr)
end do
end