subroutine d1(nmax,y,id,m,ncmin,pln,eln,nmin,kmax)
      parameter(im=100,ii=100000000,nx=50000,tiny=1e-20) 
      dimension y(nmax),jh(0:im*im),ju(nx),d(nx),jpntr(nx),nlist(nx)

      if(nmax.gt.nx) pause "Make nx larger."
      ncomp=nmax-(m-1)*id
      k=int(exp(pln)*(ncomp-2*nmin-1))+1
      if(k.gt.kmax) then
        ncomp=real(ncomp-2*nmin-1)*real(kmax)/k+2*nmin+1
        k=kmax
      endif         
      pln=psi(k)-log(real(ncomp-2*nmin-1))
      print*, 'Mass ', exp(pln),': k=', k, ', N=', ncomp     !c
      eps=exp(pln/m)*sigma(nmax,y)
      iu=ncmin-(m-1)*id
      do 10 i=1,iu
 10     ju(i)=i+(m-1)*id
      eln=0
 1    call base(ncomp+(m-1)*id,y,id,m,jh,jpntr,eps)
      iunp=0
      do 20 nn=1,iu                                           ! find neighbours
        n=ju(nn)
        call neigh(nmax,y,n,nmax,id,m,jh,jpntr,eps,nlist,nfound)
        nf=0
        do 30 ip=1,nfound
          np=nlist(ip)
          nmd=mod(abs(np-n),ncomp)
          if(nmd.le.nmin.or.nmd.ge.ncomp-nmin) goto 30    ! temporal neighbours
          nf=nf+1
          dis=0
          do 40 i=0,m-1
 40         dis=max(dis,abs(y(n-i*id)-y(np-i*id)))
          d(nf)=dis
 30       continue
        if(nf.lt.k) then
          iunp=iunp+1                                     ! mark for next sweep
          ju(iunp)=n
        else
          e=which(k,nf,d)
          eln=eln+log(max(e,tiny))
        endif
 20     continue
      iu=iunp
      print*, eps, iunp  !c
      eps=eps*sqrt(2.)
      if(iunp.ne.0) goto 1
      eln=eln/(ncmin-(m-1)*id)
      end