subroutine mpenrgy(p1,p2,wp)
c
c Volume integral of p1*p2
c
      implicit none
      include 'itg.par'
      real p1(lz,mz,nz,nspecz),p2(lz,mz,nz,nspecz),
     *  wpp(mz,nz),wp,coef,zero(mz),dz
      integer l,m,n
      include 'itg.cmn'
c
c     note dr(l)=2.0*dx                   !  Bill's
c     normalize by 0.5/x0 in the x-direction
c           coef=0.5/(x0-2.*xp)
c
c     dz normal
c     normalize by 2*x0 or 2*xp
c
      wp=0.0
c     
      if(iperiod.ne.2) then
         dz=2.*x0/ldb
         coef=1./(2.*x0)
         l_left=1
         l_right=ld
      else
         dz=2.*xp/ldb
         coef=1./(2.*xp)
      endif
c
c The m=0 mode only contributes half as much as the other modes to the 
c total volume integral.  This is clumsy looking, but it works:
c
      zero(1)=1.0
      do m=2,md
         zero(m)=2.0
      enddo

      do n=1,nd
         do m=1,md
            wpp(m,n)=0.0
            do l=l_left,l_right
               wpp(m,n)=wpp(m,n)+p1(l,m,n,1)*p2(l,m,n,1)*dz
            enddo
            wpp(m,n)=coef*wpp(m,n)*zero(m)
            wp=wp+wpp(m,n)
         enddo
      enddo
c
      return
      end