subroutine levy
       implicit none
       include 'pimc.par'
       integer il,l,js,jsm,jsp,i,is,nj,mi,ik,jkmin,k
       integer lm,jm,jk,step,iflag,jl,lvl,jdn,jup,jmid,kdiv
       integer mover(mmovers),itspan,itar,iorg,kspan,ndiv,js1,is1
       real*8 eps,sg,xold(mdim,0:mmovers),pold,pnew,xm,ppo,ppn,pwrap
       real*8 deltac,cold,cnew
       real*8 dio,din,ysold,ysnew,pold1,pnew1,ronf,deltaew,ewold,ewnew
       real*8 dbgact,dbgepact,dbgepact1,dbgep12,dbgcpot,epold,epnew
       real*8 gaussn,coulomb
       real*4 etime,tim(2)
       parameter (eps=1.d-32)
       external gaussn,coulomb
       save itspan,iflag,lvl
       character*80 messagge
       include 'pimc.cm'
       data iflag/0/
c do one block of MC dynamics
c Levy path reconstruction by bisection
c set the number of movers to a power of 2
       if (iflag.eq.0) then
        iflag=1
        lvl=dint(dlog(dfloat(nmovers))/dlog(2.d0))
        nmovers=2**lvl-1
        itspan=nmovers+1
        write(*,*) 'effective number of movers=',nmovers
        write(2,*)
        write(2,*) 'effective number of movers=',nmovers
        nrap=nslices/max(1,nmovers)
       endif
       do step=1,nstep
        if (nmovers.gt.0) then
c select the origin at random
         call rmar(ronf)
         iorg=dint(nslices*ronf)+1
         do nj=1,np
c loops recursively over time origins
          do lm=1,nrap
           iorg=iwrap(iorg+nmovers)
           itar=iwrap(iorg+itspan)
           js=iorg
           do l=1,ndim
            xold(l,0)=x(l,nj,js)
           enddo
c loops over levels
           do jm=1,nmovers+1
            mover(jm)=iwrap(js+jm)
            do l=1,ndim
             xold(l,jm)=x(l,nj,mover(jm))
            enddo
c           write (*,*) js,iwrap(js+itspan),mover(jm)
           enddo
           kspan=itspan
c          goto 777
           ndiv=1
           do jl=lvl,1,-1
            jdn=js
            kspan=kspan/2
            sg=dsqrt(dfloat(kspan))*sigma
            do kdiv=1,ndiv
             jmid=iwrap(jdn+kspan)
             jup=iwrap(jmid+kspan)
             do l=1,ndim
c compute the trial free particle term
              xm=(x(l,nj,jdn)+x(l,nj,jup))/2.d0
              x(l,nj,jmid)=xm+sg*gaussn()
             enddo
             jdn=jup
            enddo
            ndiv=2*ndiv
           enddo
c777        continue
           deltac=0.d0
c computes the coulomb trial action
!          if (charge.ne.0.d0) then
            cold=0.d0
            cnew=0.d0
            do jm=1,nmovers+1
             js=mover(jm)
             js1=iwrap(js-1)
             do mi=1,np
              if (mi.ne.nj) then
               cold=cold+coulomb(xold(1,jm),x(1,mi,js),
     &                           xold(1,jm-1),x(1,mi,js1),ndim)
               cnew=cnew+coulomb(x(1,nj,js),x(1,mi,js),
     &                           x(1,nj,js1),x(1,mi,js1),ndim)
              endif
             enddo
             do mi=1,2
              cold=cold-coulomb(xold(1,jm),proton(1,mi)
     &                         ,xold(1,jm-1),proton(1,mi),ndim)
              cnew=cnew-coulomb(x(1,nj,js),proton(1,mi)
     &                         ,x(1,nj,js1),proton(1,mi),ndim)
             enddo
            enddo
            deltac=-cnew+cold
!          endif
           ntry=ntry+1
           call rmar(ronf)
           ronf=dlog(max(eps,ronf))
           if ( deltac.ge.ronf) then
c    ** .. accept **
            nacc = nacc + 1
c update the action
            cpot=cpot-deltac
            action=cpot
           else
c    **  .. reject and replace **
            do jm=1,nmovers
             js=mover(jm)
             do l=1,ndim
              x(l,nj,js) = xold(l,jm)
             enddo
            enddo
           endif  ! accept or reject
          enddo !  over nrap
         enddo  ! over particles
        endif   ! (nmovers>0)
c attempt to displace the entire path
        call rmar(ronf)
        if (np.gt.1.and.ronf.le.gammadsp) call displace
c computing averages
        if (mod(step,nanal).eq.0) then
         tt=etime(tim)
         if (idbg.ge.2) then
          messagge='levy before anal: step'
!         call compactdbg(messagge,step)
         endif
         call anal(1,step)
         ttanal=ttanal+etime(tim)-tt
        endif
        if (mod(step,nspill).eq.0)
     &      write (31) (((x(l,nj,js),l=1,ndim),nj=1,np),js=1,nslices)
     $                ,((proton(l,nj),l=1,ndim),nj=1,2)

       enddo         ! over steps

       return
       end