      subroutine omove(iflag,p,n)
#define SIMPLE_SPRNG
#include "sprng_f.h"
c driver for one body moves
      include 'upi.cm'
      integer iflag,n,ix,ipack,nsteps,levels,itspan,iconf,l
      integer nmovers,istp,itorg,itar,ii,movers,intread
      real*8 action,t0,t1
      dimension ix(mpack),ipack(mpack4),movers(mnmovers)
      character p(n)*(*)
      save nsteps,levels,itspan,iconf
      include 'cpbc.cm'

! janor add
      integer ijar,jjar,kjar
      
      integer nbinjar(3),nbinr
      real*8 rbinjar
      real*8 xtmp,ytmp,ztmp,rtmp,Vtmp
      real*8 rxtmp,rytmp,rztmp,rrtmp
      logical b1jar,b2jar,b3jar
! end janor add

      if(iflag.eq.0) then
c interpret the input
        if(n.ne.4) then
          write (*,*) ' input data in omove inappropriate'
          stop
        endif

        nsteps=intread(p(3))
        levels=intread(p(4))
        call mcheck(levels,mlevels,'levels','mlevels','omove')
        write (6,5) nsteps,levels
5       format(' number of single particle steps/block ',i5
     +  /' levels ',i5)
        write (6,8) drfmax
8       format(' maximum drift ' ,e12.5)
        write (6,*) ' gamma= ',(gamma(1,l),l=1,ntypes)
        itspan=2**levels
        if(itspan.gt.nslices)then
          write (6,*) ' itspan.gt.nslices in omove',itspan,nslices
          stop
        endif
        iconf=0
      endif
c initialize fermion matrices
      call fdminit(iref,irev,action,iflag)
      call second(t0)
      nmovers=1
      nvstp=0
      do istp=1,nsteps

c select time origin at random
      itorg=1+sprng()*nslices
      itar=iwrap(itorg+itspan)
      call pmove(itar)

c try a move of each particle with hbs2m>0
      do ii=1,nparts
      movers(1)=ii
      nvstp=nvstp+1
      if(sprng().lt. gamma(1,idtype(ii) ) )
     +     call weave(itorg,nmovers,movers,levels)
       enddo


c find averages
      call anal(itorg,itspan,1)

! janor add
      ntjar=ntjar+1

      if(write_r2_file)then
         xtmp=0.d0
         ytmp=0.d0
         do kjar=1,nslices
            do jjar=nfpty(1),nlpty(1)
               xtmp=xtmp+r(1,jjar,kjar)**2
               ytmp=ytmp+r(2,jjar,kjar)**2
            end do
         end do
         write(58,580)xtmp+ytmp
 580     format(1pe12.4)
      end if

      if(ntjar.gt.ntsteadyjar)then
         do kjar=1,nslices
            do jjar=nfpty(1),nlpty(1)
               nbinjar(1)=(r(1,jjar,kjar)+.5*ell(1))/dxdistjar+1.d0
               nbinjar(2)=(r(2,jjar,kjar)+.5*ell(2))/dydistjar+1.d0
               nbinjar(3)=(r(3,jjar,kjar)+.5*ell(3))/dzdistjar+1.d0
               b1jar=nbinjar(1).lt.1 .or. nbinjar(1).gt.nxdistjar
               b2jar=nbinjar(2).lt.1 .or. nbinjar(2).gt.nydistjar
               b3jar=nbinjar(3).lt.1 .or. nbinjar(3).gt.nzdistjar
               if(b1jar.or.b2jar.or.b3jar)then
                  write(*,*)'You did the density binning wrong!'
                  stop
               else
                  densdistxy(nbinjar(1),nbinjar(2))=
     &                densdistxy(nbinjar(1),nbinjar(2))+1.d0
                  densdistxz(nbinjar(1),nbinjar(3))=
     &                densdistxz(nbinjar(1),nbinjar(3))+1.d0
               end if

               rbinjar=sqrt(r(1,jjar,kjar)**2+r(2,jjar,kjar)**2)
               nbinr=rbinjar/drdistjar+1.d0
               if(nbinr.le.nrdistjar)then
                  densdistr(nbinr)=densdistr(nbinr)+1.d0
               end if
            end do
         end do      

         if(mod(ntjar,ntskipjar).eq.0)then
            open(unit=54,file='dens_dist_params.dat',status='unknown')
            open(unit=55,file='dens_distxy.dat',status='unknown')
            open(unit=56,file='dens_distxz.dat',status='unknown')
            open(unit=57,file='dens_distr.dat' ,status='unknown')
            rxtmp=ntjar-ntsteadyjar
            write(54,540)rxtmp,rxtmp,rxtmp,rxtmp
            rxtmp=nxdistjar
            rytmp=nydistjar
            rztmp=nzdistjar
            rrtmp=nrdistjar
            write(54,540)rxtmp,rytmp,rztmp,rrtmp
            write(54,540)dxdistjar,dydistjar,dzdistjar,drdistjar
            write(54,540)ell(1),ell(2),ell(3),ellrjar
540         format(4(1pe16.8))
            do jjar=1,nydistjar
               do ijar=1,nxdistjar
                  xtmp=(ijar-.5)*dxdistjar-.5*ell(1)
                  ytmp=(jjar-.5)*dydistjar-.5*ell(2)
                  write(55,550)xtmp,ytmp,densdistxy(ijar,jjar)
               end do
            end do
 550        format(3(1pe12.4))
            do kjar=1,nzdistjar
               do ijar=1,nxdistjar
                  xtmp=(ijar-.5)*dxdistjar-.5*ell(1)
                  ztmp=(kjar-.5)*dzdistjar-.5*ell(3)
                  write(56,560)xtmp,ztmp,densdistxz(ijar,kjar)
               end do
            end do
 560        format(3(1pe12.4))
            do ijar=1,nrdistjar
               rtmp=(ijar-.5)*drdistjar
               Vtmp=3.141592d0*rtmp*drdistjar*ell(3)
               write(57,570)rtmp,Vtmp,dendistr(ijar)
            end do
 570        format(3(1pe12.4))
            close(unit=54)
            close(unit=55)
            close(unit=56)
            close(unit=57)
         end if
      end if

! end janor add

c reference point move?
!     if(iref.ne.0)call refmove
c write out configuration?
      iconf=iconf+1
      if(mod(iconf,iwrc).eq.0) then
      call wpack(r,nparts,ndim,el2,ix,ipack,nslices,mnp,join
     +,iref,irev,iperm,mcpw)
      endif
c try insertion?
!     if(sprng().lt.gammai)call insert
c try displacement?
      call displace
c check time?
        if(tcheck.gt.0) then
         call second(t1)
         if(t1-t0.gt.tcheck) then
         write (*,*) ' istp = ',istp
         go to 100
         endif
        endif
      enddo

100   continue
      call second(t1)
      anorm(jdcr)=t1-t0
      call pars

      return
      end
