dirasc = "./"
        filasc = systemfunc("cd "+dirasc+" ; ls test.vanucia.txt")   
        pthasc = dirasc+filasc
        nfil   = dimsizes(filasc)
        print(filasc)
        print("----")

;------ Create grid coordinates [double] 
		
        lonW  = (/-119.995d/)    
        lonE  = (/-24.005d/) 

        latS  = (/-55.995d/)  
        latN  = (/32.995d/)   

       ;dgrad = 0.01d         ; vary small grid spacing
	dgrad = 1.00d   

;----- Read the values in as 1D array of strings 
        data  = asciiread(pthasc,-1,"string")

;----- Open pc variable (for ALL latitude/longitude point)
        lat   = stringtodouble( str_get_field(data(1:), 8, ",") )    
        lon   = stringtodouble( str_get_field(data(1:), 9, ",") )  
        pc    = stringtodouble( str_get_field(data(1:), 10, ",") )
        npc   = num(.not.ismissing(pc))

;----- Select pc variable within the grid area
        ii    = ind(lat.ge.latS .and. lat.le.latN .and. \
                    lon.ge.lonW .and. lon.le.lonE ) 
        lat  := lat(ii)        ; reassign 
        lon  := lon(ii)
	var   = pc(ii)
        N     = dimsizes(var)
        delete(pc)

print("npc="+npc+"   N="+N)   
print("min(lat)="+min(lat)+"   max(lat)="+max(lat))
print("min(lon)="+min(lon)+"   max(lon)="+max(lon))
print("min(var)="+min(var)+"   max(var)="+max(var))
       

	nlat = toint(abs(latN-latS)/dgrad)+1
	mlon = toint(abs(lonE-lonW)/dgrad)+1
print("nlat="+nlat+"   mlon="+mlon)     

	lat2 = fspan(latS,latN,nlat)
	lon2 = fspan(lonW,lonE,mlon)

	lat2@long_name = "latitude" 
	lat2@units     = "degrees_north" 
	lat2!0         = "lat2" 
	lat2&lat2      =  lat2

	lon2@long_name = "longitude" 
	lon2@units     = "degrees_east" 
	lon2!0         = "lon2" 
	lon2&lon2      =  lon2
	
;----- Create a new var to pc's values
	new_data           =  new((/nlat,mlon/),"double", 1d20)
	new_data           =  0     ; start
	new_data!0         = "lat2" 
	new_data!1         = "lon2" 
	new_data&lat2      =  lat2
	new_data&lon2      =  lon2
printVarSummary(new_data)
printMinMax(new_data, 0)
print("---------------")

;------Sum pc values
	do i=0,N-1
           new_data({lat(i)},{lon(i)}) = new_data({lat(i)},{lon(i)})+(/var(i)/) 
	end do
      
	new_data = where(new_data.ge.0, new_data@_FillValue,  new_data) 
printVarSummary(new_data)
printMinMax(new_data, 0)
print("---------------")

;---- Selecting just non-missing values
    
data1d  := ndtooned(new_data)            ; reduce to 1D array
ivalid  := ind(.not.ismissing(data1d))   ; extract only the non-missing values.
xvalid  := data1d(ivalid)                ; to float
printVarSummary(xvalid)
printMinMax(xvalid, 0)
print("---------------")
         
 system("/bin/rm -f xvalid.txt")
 asciiwrite("xvalid.txt", sprintf("%9.2f", xvalid)) 

;------------ Plot 
wks = gsn_open_wks("png","test.vanucia")

;;setvalues NhlGetWorkspaceObjectId()
;;  "wsMaximumSize" : 900000000
;;end setvalues

   res                     = True    
   res@mpFillOn            = False              
   res@gsnAddCyclic        = False           
   res@cnFillOn            = True            
   res@cnLinesOn           = False
   res@cnLineLabelsOn      = False
   res@cnInfoLabelOn       = False
   res@cnFillMode          = "RasterFill"       ; Raster Mode 

   res@mpMinLatF           = latS     ; range to zoom in on
   res@mpMaxLatF           = latN
   res@mpMinLonF           = lonW
   res@mpMaxLonF           = lonE
  
   mnmxint = nice_mnmxintvl( min(new_data), max(new_data), 18, False)
   res@cnLevelSelectionMode = "ManualLevels"
   res@cnMinLevelValF       = mnmxint(0)
   res@cnMaxLevelValF       = mnmxint(1)
   res@cnLevelSpacingF      = mnmxint(2)
print(res)

   plot = gsn_csm_contour_map(wks,new_data,res)