subroutine anal(lvl,step)
      implicit none
      include 'pimc.par'
      integer i,k,j,l,jsp,js,is,il,ik,ncount(mslices),lvl,im,nl
      integer intv,nsl,nt,jm,isp,ism,nj,mi,ind,inf,js1,is1,it
      integer ist,ip,kp,step
      real*8 rg(2),xcm(mdim,mparts),etot(3),ekin,epot(5),cm(mslices)
      real*8 d,evir,grsum(lptable),dspsum(ldsptable),uu,dist(mdim)
      real*8 dd,dd2,yself,evirc,nact,evirclr,etephlr
      real*8 coulomb,dphonon3,dvir3,dcoulomb,coulvir,ppbc
      real*8 sksum(mnsh),dxcm(mdim),dump(mdim)
      external coulomb,dphonon3,dvir3,dcoulomb,coulvir
      data ist/0/
      save ist
      include 'pimc.cm'
      ppbc(ax,ay,l)=anint((ax-ay)/ell(l))*ell(l)
c
c  compute block average
c
      intv=2**(lvl-1)
      nsl=nslices/intv
      do nt=1,intv
c
c  observables
c
       ekin=0.0
       epot(1)=0.0
       epot(2)=0.0
       epot(3)=0.0
       epot(4)=0.0
       epot(5)=0.0
       etot(1)=0.0
       etot(2)=0.0
       etot(3)=0.0
       do l=1,ndim
        do nj=1,np
         xcm(l,nj)=0.d0
        enddo
       enddo
       rg(1)=0.d0
       rg(2)=0.d0
       evir=0.d0
       evirc=0.d0
       evirclr=0.d0
       etephlr=0.d0
       do k=1,nsl
        cm(k)=0.d0
        ncount(k)=0.d0
       enddo
       if (np.gt.1) then
        do k=1,lptable
         grsum(k)=0.0d0
        enddo
        do k=1,nlamb
         sksum(k)=0.0d0
        enddo
       endif
       do k=1,ldsptable
        dspsum(k)=0.0d0
       enddo
       
c prepare virial estimator tables
       call cvirtab

c compute the path centroid needed in coulvir
       do nj=1,np
        do jm=1,nsl
         js=iwrap(nt+(jm-1)*intv)
         do l=1,ndim
          xcm(l,nj)=xcm(l,nj)+x(l,nj,js)
         enddo
        enddo
        do l=1,ndim
         xcm(l,nj)=xcm(l,nj)/nsl
        enddo
       enddo

       do nj=1,np
        do jm=1,nsl
         js=iwrap(nt+(jm-1)*intv)
         jsp=iwrap(js+intv)
         js1=iwrap(js-1)
         dd2=0.d0
         do l=1,ndim
          dd2=dd2+(x(l,nj,js1)-x(l,nj,js))**2
         enddo
         ekin=ekin+dd2
         if (charge.ne.0.0) then
          do mi=1,2
           epot(4)=epot(4)-zp*dcoulomb(x(1,nj,js),proton(1,mi)
     &                              ,x(1,nj,js1),proton(1,mi),ndim)
          enddo
         endif
        enddo     ! over js
c intermolecular contribution
        do mi=1,nj-1
         do l=1,ndim
          dump(l)=ppbc(xcm(l,nj),xcm(l,mi),l)
          dxcm(l)=xcm(l,nj)-xcm(l,mi)-dump(l)
         enddo
         do jm=1,nsl
          js=iwrap(nt+(jm-1)*intv)
          js1=iwrap(js-1)
          if (charge.ne.0.0d0) then
           epot(4)=epot(4)+dcoulomb(x(1,nj,js),x(1,mi,js)
     &                             ,x(1,nj,js1),x(1,mi,js1),ndim)
!          evirc=evirc+coulvir(x(1,nj,js),x(1,mi,js)
!    &                        ,x(1,nj,js1),x(1,mi,js1)
!    &                        ,dxcm,dump,ndim)
!    &                        ,x(1,nj,1),x(1,mi,1),ndim)
          endif
         enddo  ! over jm
        enddo   ! over mi
        
        if (charge.ne.0.0d0) then
         do js=1,nslices
          js1=iwrap(js+1)
          do is=1,nslices
           is1=iwrap(is+1)
           do l=1,ndim
            evirc=evirc
     &           +(x(l,nj,js)-x(l,nj,is))*
     &             ( ( aja(l,nj,js) + aja(l,nj,js1)
     &                -aja(l,nj,is) - aja(l,nj,is1) 
     &               )/2.d0 
     &               +( bja(l,nj,js) + bja(l,nj,is)
     &                 -bja(l,nj,js1)- bja(l,nj,is1) ) 
     &             )
           enddo
          enddo
         enddo
        endif
       enddo    ! over nj

       ekin=(dble(ndim)/2.d0-cke*ekin/dfloat(nslices*np))/(tau*intv)
       epot(2)=epot(2)*(intv**2)/dfloat(np)
       epot(3)=(epot(3)*(intv**2)+etephlr)/dfloat(np)
       epot(4)=(cpot/beta+epot(4)*dble(intv)/dble(nslices))/dble(np)
       epot(1)=a11*epot(2)+a12*epot(3)+epot(4)
       epot(5)=a11*(epact+epact1)+a12*ep12+cpot
c      write (*,*) 'anal:',epact,epact1,ep12,cpot
       evir=dble(ndim)/2.d0/beta
     &     +a11*evir/2.d0/dfloat(nsl*np)/(tau*intv)
     &     +evirc*charge/4.d0/dble(np*nslices*nslices)+evirclr
        etot(1)=ekin+epot(1)
        etot(3)=evir+epot(1)
        if (np.eq.1) then
         etot(1)=etot(1)-dble(ndim)/2.d0/beta
         etot(3)=etot(3)-dble(ndim)/2.d0/beta
        endif
       etot(2)=epot(5) !                 exp(-epot(5)/12.)
       do nj=1,np
        do js=nt,nslices-1+nt,intv
         do l=1,ndim
          rg(1)=rg(1)+(pbc(x(l,nj,js),xcm(l,nj),l))**2
         enddo
        enddo
       enddo
       if (np.eq.2) then
        ip=1
        do l=1,ndim
         rg(2)=rg(2)+(pbc(xcm(l,ip),xcm(l,ip+1),l))**2
        enddo
       endif
       rg(1)=rg(1)/dfloat(nsl*np)
c effective mass
       do nj=1,np
        js=nt
        do jm=1,nsl-1
         do ik=1,jm-1
          is=iwrap(js-ik*intv)
          k=ik
          if (k.gt.nsl/2) k=nsl/2-(k-nsl/2)
c         write (*,*) js,is,k
          do l=1,ndim
           cm(k)=cm(k)+(pbc(x(l,nj,js),x(l,nj,is),l))**2
          enddo
          ncount(k)=ncount(k)+1
         enddo
         js=js+intv
        enddo 
       enddo
       do k=1,nsl/2
        cm(k)=cm(k)/(2.d0*ndim*alamb*tau*intv*k)
     &              /dfloat(max(ncount(k),1))
       enddo
c g(r) and S(k)
!      if (np.gt.1) then
        do jm=1,nsl
         js=iwrap(nt+(jm-1)*intv)
         do i=1,np
          do j=1,1
           d=0.0
           do l=1,ndim
            d=d+(pbc(x(l,i,js),proton(l,j),l))**2
           enddo
           if (dsqrt(d).le.cutr) then
            ind=1+int(dsqrt(d)*dpi)
            if (ind.gt.lptable) ind=lptable
            grsum(ind)=grsum(ind)+1.d0
           endif
          enddo
         enddo
        enddo
!      endif
         
c update local averages
       call cumul1(ekin,avp(ikin),anormp(ikin),1)
       call cumul1(evir,avp(iev),anormp(iev),1)
       call cumul1(epot,avp(iepot),anormp(iepot),5)
       call cumul1(etot,avp(ietot),anormp(ietot),3)
       call cumul1(rg,avp(irg),anormp(irg),2)
       call cumul1(cm,avp(ieff),anormp(ieff),nsl/2)
       call cumul1(grsum,avp(igr),anormp(igr),lptable)
 
      enddo     ! loop over nt

      return
      end