subroutine volsqk(a,b,axb,axb_by_mode)
c
c Volume integral of q*p, where a and b are pitch angle moments 
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 ??BD 9/22/95:
c Not yet generalized for igeo=1.  (I haven't checked the kappa dependence.)

      implicit none
      include 'itg.par'
      include 'itg.cmn'
      complex a(kz,mz,nz),b(kz,mz,nz)
      real axb,zero(mz),c(kz,mz,nz),d(lz,mz,nz)
      real axb_by_mode(mz,nz),dz(lz),vnorm
      integer l,m,n,j,l0
c
      if(iperiod.ne.2) then
         l_left=1
         l_right=ld-1
      endif

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

      zero(1)=1.0
      do m=2,md
c         zero(m)=2.0
         zero(m)=0.5
      enddo
c
c     first multiply pitch angle functions together
c
      do n=1,nd
         do m=2,md
            do l=1,kd
               c(l,m,n)=real(a(l,m,n))*real(b(l,m,n))
     .                 +aimag(a(l,m,n))*aimag(b(l,m,n))
            enddo
         enddo
         m=1
            do l=1,kdpass
               c(l,m,n)=real(a(l,m,n))*real(b(l,m,n))
     .                 +aimag(a(l,m,n))*aimag(b(l,m,n))
            enddo
      enddo
c
c     now average in kappa
c
      do m=2,md
         do n=1,nd
            do j=1,ldb
               l0=1+abs(j-ldb/2-1)
               d(j,m,n)=0.
               do l=l0,kd
                  d(j,m,n)=d(j,m,n)+dkap(l)*c(l,m,n)*wgtba2(l,j)
               enddo
            enddo
         enddo
      enddo

      m=1
         do n=1,nd
            do j=1,ldb
               l0=1+abs(j-ldb/2-1)
               d(j,m,n)=0.
               do l=l0,kd-1
                  d(j,m,n)=d(j,m,n)+dkap(l)*(c(l,m,n)
     .                 +c(l+1,m,n))/2*wgtba2(l,j)
               enddo
               do l=max(kd,l0),kd+1
                  d(j,m,n)=d(j,m,n)+dkap(l)*c(l,m,n)*wgtba2(l,j)
               enddo
               do l=kd+2,kdpass
                  d(j,m,n)=d(j,m,n)+dkap(l)*(c(l,m,n)
     .                 +c(l-1,m,n))/2*wgtba2(l,j)
               enddo
            enddo
         enddo
c
c     finally volume average
c
      axb=0.0
      do n=1,nd
         do m=1,md
            axb_by_mode(m,n)=0.
            vnorm=0.
            do l=l_left,l_right
               axb_by_mode(m,n)=axb_by_mode(m,n)+d(l,m,n)*dz(l)*zero(m)
               vnorm=vnorm+dz(l)
            enddo
            axb_by_mode(m,n)=axb_by_mode(m,n)/vnorm
            axb=axb+axb_by_mode(m,n)
         enddo
      enddo
c
      return
      end