subroutine clean(nmax,y,yc,m,kmin,nq,d)
      parameter(im=100,ii=100000000,nx=50000,mm=15,small=0.0001) 
      dimension y(nmax),yc(nmax),r(mm),ju(nx),c(mm,mm),cm(mm),
     .  jh(0:im*im),jpntr(nx),nlist(nx)
      dimension zcm(mm,nx)                                                ! \TR

      if(nmax.gt.nx.or.m.gt.mm) pause "Make mm/nx larger."
      sr=2*small+m-2                                        ! ${\rm tr}(1/r)=1$
      do 10 i=1,m
        r(i)=sr
 10     if(i.eq.m.or.i.eq.1) r(i)=sr/small
      do 20 i=1,nmax
 20     yc(i)=y(i)
      do 30 istep=1,2                                                     ! \TR
        eps=d
        iu=nmax-m+1
        do 40 i=1,iu
 40        ju(i)=i+m-1
 1      call base(nmax,y,1,m,jh,jpntr,eps)
        iunp=0
        do 50 nn=1,iu                                         ! find neighbours
          n=ju(nn)
          call neigh(nmax,y,n,nmax,1,m,jh,jpntr,eps,nlist,nfound)
          if(nfound.lt.kmin) then                 ! not enough neighbours found
            iunp=iunp+1                                   ! mark for next sweep
            ju(iunp)=n
          else                                        ! fine: enough neighbours
            do 90 i=1,m                                 ! centre of mass vector
              s=0
              do 100 np=1,nfound
 100            s=s+y(nlist(np)-m+i)
 90           cm(i)=s/nfound
            if(istep.eq.1) then               ! just store centre of mass ! \TR
              do 110 i=1,m                                                ! \TR
 110            zcm(i,n)=cm(i)                                            ! \TR
            else                                                          ! \TR
              do 120 i=1,m              ! corrected centre of mass vector ! \TR
                s=0                                                       ! \TR
                do 130 np=1,nfound                                        ! \TR
 130              s=s+zcm(i,nlist(np))                                    ! \TR
 120            cm(i)=2*cm(i)-s/nfound                                    ! \TR
              do 140 i=1,m                          ! compute covariance matrix
                do 140 j=i,m
                  s=0
                  do 150 np=1,nfound
                    jm=nlist(np)-m
 150                s=s+(y(jm+i)-cm(i))*(y(jm+j)-cm(j))
                  c(i,j)=r(i)*r(j)*s/nfound
 140              c(j,i)=c(i,j)
              call eigen(c,m)                  ! find eigenvectors (increasing)
              do 160 i=1,m
                s=0
                do 170 iq=1,nq
                  do 170 j=1,m
 170                s=s+(y(n-m+j)-cm(j))*c(i,iq)*c(j,iq)*r(j)
 160            yc(n-m+i)=yc(n-m+i)-s/r(i)/r(i)
            endif
          endif
 50       continue
        iu=iunp
        print*, "With ", eps, iunp, " uncorrected" !c
        eps=eps*sqrt(2.)
 30     if(iunp.ne.0) goto 1
      end