subroutine volsq(a,b,axb,axb_by_mode)
c
c Volume integral of q*p / Volume
c     
c     wp is volume integral over 
c            (-zlen/2,zlen/2)(-ylen/2,ylen/2)(-xlen/2,xlen/2)
c     wpp is z integral for each m,n
c
c     where zlen is 2.*xp if iperiod=2 and 2.*x0 otherwise.
c     
c BD 6/22/95:
c     dVol picks up a factor of jacobian in general geometry.
c     This would change our previous answers slightly, since we 
c     ignored 1/B(theta) in the previous volume averages.  To recover 
c     the earlier results, the variable "jacobian" should be independent 
c     of theta for igeo=0.  However, this doesn't work because the 
c     jacobian factor in poisson retained its theta dependence.  Thus, 
c     the logic here may look a little odd: ignore the jacobian factor 
c     if igeo=0, (to compare with old results) but not for igeo=1, even 
c     if itor=0.  
c
c     This means that igeo=1, itor=0 nonlinear runs will not give same
c     result as igeo=0 runs!
c     
      implicit none
      include 'itg.par'
      include 'itg.cmn'
      complex a(lz,mz,nz),b(lz,mz,nz)
      real axb_by_mode(mz,nz),axb,zero(mz),vnorm,dz0,dz(lz)
      integer l,m,n,init
      data init/1/
      save init
c
      dz0=2.*x0/ldb
      if(iperiod.ne.2) then
         l_left=1
         l_right=ld-1
      endif

      axb=0.0
c      ylen=2.*pi*y0
c      xlen=2.*pi*y0*(float(nd-1)/2.)/(shr*z0)
c
c This logic is wrong for input md=-1.     
c
      zero(1)=1.0
      do m=2,md
         zero(m)=0.5
      enddo

      do l=l_left,l_right
         if(igeo.eq.0) then
            dz(l)=dz0
         else
            dz(l)=dz0*jacobian(l)
         endif
      enddo

      do n=1,nd
         do m=1,md
            axb_by_mode(m,n)=0.
            vnorm=0.
            do l=l_left,l_right
               vnorm=vnorm+dz(l)
               axb_by_mode(m,n)=axb_by_mode(m,n)
     .              +(real(a(l,m,n))*real(b(l,m,n))
     .              +aimag(a(l,m,n))*aimag(b(l,m,n)))*dz(l)
            enddo
            axb_by_mode(m,n)=axb_by_mode(m,n)/vnorm*zero(m)
            axb=axb+axb_by_mode(m,n)
         enddo
      enddo
c
c Debug:
c      if(init.eq.1) write(*,*) 'int(J dtheta)= ',vnorm
c      if(init.eq.1) init=0
c
      return
      end