program sft2
      parameter (np0=2000)
      character*80 s80
      dimension xx(np0),yy(np0),si(np0),sm(np0),sm2(np0),
     1sm3(np0),nspike(np0),sum1(np0),sum2(np0)
c     nspike   . number of spikes for given frequency

      ia=iargc()
      if(ia.gt.0)then
       call getarg(1,s80)
       read(s80,*)alim
      else
       alim=3.
      endif

      ns=0
      open(9,file='spec.lst',status='old')
1     read(9,900,end=99,err=99)s80
900   format(a80)
      ns=ns+1
      goto 1
99    close(9)
      nst=ns
      write(6,*)ns,' spectra'

      open(9,file='spec.lst')
      do 10 ns=1,nst
      read(9,900)s80
      open(91,file=s80)
      np=0
2     read(91,*,end=88,err=88)x,y
      np=np+1
      if(np.gt.np0)call report('np>np0')
      xx(np)=x
      yy(np)=y
      goto 2
88    close(91)
      if(ns.eq.1)then
       write(6,*)np,' points'
       np1=np
       do 3 i=1,np
       nspike(i)=0
       sum1(i)=0
       sum2(i)=0
       si(i)=yy(i)
       sm(i)=yy(i)
       sm2(i)=yy(i)
3      sm3(i)=yy(i)
      else
       if(np.ne.np1)call report('different number of points')
      endif
      do 4 i=1,np
      ave=(si(i)+sm(i)+sm2(i)+sm3(i))/4.0
      ader=abs(yy(i)-si(i))
      if(ave.ne.0.)ader=ader/abs(ave)
      if(xx(i).gt.2301.0.and.xx(i).lt.2302.0)
     1write(6,*)ns,ave,yy(i),ader,alim
      if(ader.gt.alim)then
       nspike(i)=nspike(i)+1
       sum1(i)=sum1(i)+yy(i)*real(ns)
       sum2(i)=sum2(i)+real(ns)
      endif
4     continue
      do 6 i=1,np
      sm3(i)=sm2(i)
      sm2(i)=sm(i)
      sm(i)=si(i)
6     si(i)=yy(i)
10    continue
      close(9)

      open(40,file='total.prn')
      open(41,file='spike.prn')
      do 7 i=1,np
c     write(6,*)xx(i),nspike(i),(nsspike(i,j),j=1,nspike(i))
      write(41,401)xx(i),nspike(i)
401   format(f18.4,20i6)
7     write(40,400)xx(i),yy(i)-sum1(i)/real(nst)
400   format(2f18.4)
      close(40)
      close(41)
      stop
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end