begin

;------ Data 

        dirasc = "./"
        filasc = systemfunc("cd "+dirasc+" ; ls 2015*.txt")   
        pthasc = dirasc+filasc
        nfil   = dimsizes(filasc)
       
;------- Read multiple files
 
    do j  = 0, nfil-1   
    
;----- Read the values in as 1D array of strings 
        datas := asciiread(pthasc(j),-1,"string")
        nline  = dimsizes(datas)
;print("nline="+nline)

; NCL code is confusing the - delimeter with a minus sign
; Change delimeter embedded with yyyy-mm-dd to yyyy mm dd
; There may be better ways but (1) converting string to character
;                              (2) converting - character to : character
;                              (3) change : character to : string

        datac := tochar(datas)      ; (nline,:)
        datac(:,4)  = tochar(":")
        datac(:,7)  = tochar(":")
        datac(:,10) = tochar(":")
        datas       = tostring(datac)    ; overwrite original
;print(datas)

        delim  = ":"
        nfld   = str_fields_count(datas(0), delim)  ; look at 1st line only
;print("nfld="+nfld)

        yyyy  := toint(str_get_field(datas,  1, delim) )
          mm  := toint(str_get_field(datas,  2, delim) )
          dd  := toint(str_get_field(datas,  3, delim) )
          hh  := toint(str_get_field(datas,  4, delim) )
          mn  := toint(str_get_field(datas,  5, delim) )
          sc  := toint(str_get_field(datas,  6, delim) )

        tunits = "hours since 2015-01-01 00:00:00" ; "seconds/hours/days since ...."
        time  := cd_inv_calendar(yyyy,mm,dd,hh,mn,sc,tunits, 0)
        time!0 = "time"
        ntimes = dimsizes(time)
;printVarSummary(time)

      delin    = ";"
        lat   := todouble( str_get_field(datas, 3, delin) )    ; fields count starts at 1
        lon   := todouble( str_get_field(datas, 4, delin) )
        VAL   := todouble( str_get_field(datas,5, delin) )
        
;print("time="+time+"  yyyy="+yyyy+"  mm="+mm+"  dd="+dd+"  hh="+hh \
;      +"  mn="+mn+"  sc="+sc+"  LAT="+LAT+"  LON="+LON+"  VAL="+VAL)

;------ Create grid coordinates to 1 km spatial resolution:
		
        lonW  = (/-60.00d/)    
        lonE  = (/-53.00d/) 
        latS  = (/-23.00d/)  
        latN  = (/-13.00d/)   
   
	dgrad = 0.01d    	
	nlat = toint(abs(latN-latS)/dgrad)+1
	mlon = toint(abs(lonE-lonW)/dgrad)+1

	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 store data	       
       new_data = new((/nlat,mlon/),"double")
       new_data           =  0      ;start   
       new_data!0         = "lat2" 
       new_data!1         = "lon2" 
       new_data&lat2      =  lat2
       new_data&lon2      =  lon2

;----- Select the data information 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  := lat            ; need to be same size from lat or lon
	var   = 1              ; value for reference (true or false - later)
        N     = dimsizes(var)
     
;------ Adding the variable - var - for each lat and lon grid:
	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.gt.0, new_data,new_data@_FillValue) 

printVarSummary(new_data)
printMinMax(new_data,0)

WRITE_NETCDF = True
varnames = (/"counts"/)
if(WRITE_NETCDF) then
  fout = "var_2015_"+j+".nc"    
  system("rm -f " + fout)
  ncdf = addfile(fout, "c")
  ncdf->$varnames(0)$ = new_data 
end if

end do 
end