subroutine volflux(a,b,axb,axb_by_mode)
c
c Volume integral of a*b*grad(rho)/volume integral of grad(rho)
c     
c     For a vector flux F, the quantity that is needed for transport 
c     calculations is dV/drho* < F.grad(rho) >, where <...> signifies a 
c     flux surface average:
c
c     < G > == Int(dtheta dalpha drho J G) / Int(J dtheta dalpha drho)
c
c                                           /
c           == Int(dtheta dalpha drho J G) / dV/drho * delta_rho
c                                         /
c
c     Taking delta_rho -> 0 is the definition of the flux surface average.  
c     Our simulation domain is small enough that averaging over the entire
c     radial domain is ok.
c     
c     Thus, dV/drho = 2.*pi*Int(J dtheta) in our coordinates, since J is
c     only a function of theta.
c
c     This routine actually calculates 
c                                       
c     < F.grad rho > / < |grad rho| > = 
c
c     Int(dtheta dalpha drho J |grad rho| F_rho) / Int(J dtheta dalpha drho)
c     -----------------------------------------------------------------------
c     Int(dtheta dalpha drho J |grad rho| ) / Int(J dtheta dalpha drho)
c
c        F.A delta(rho)
c     =  --------------
c         A  delta(rho)
c
c     so that, for example, the heat flux that is reported is the heat
c     flux per unit area. 
c
c     This quantity is calculated assuming that F is an ExB type flux, and 
c     that a and b are the parts whose product is the generalized radial 
c     (rho direction) component of F, denoted F_rho.
c
c     Note that < |grad rho| > * dV/drho is the area of the flux surface.
c
c     dV/drho == 2.*pi*Int(dtheta J(theta)), where -N*pi <= theta <= N*pi.
c     dV/drho == 2.*pi*vnorm here (the 2*pi is from Int(dalpha))
c
c     This modification for igeo=1 is useful because it reduces to the 
c     previous definition used for the volume-averaged flux in the 
c     circular flux surface limit.  Also, normalization issues associated 
c     with our box actually extending multiples of 2*pi in the theta 
c     direction are cleanly handled with this choice, as follows:
c
c     Our reported heat flux fits into the the fundamental transport 
c     equation 
c
c     (3/2) d nT/dt + 1/V' d/drho V' < Q.grad rho> + ... = 0
c
c     like this: 
c
c     (3/2) d nT/ dt + 1/V' d/drho A Q_itgc + ... = 0
c     
c     where Q_itgc is the heat flux calculated with this subroutine, 
c     and A is the area of the flux surface.  In this equation, 
c     V' should be the physical V' (that is, using pi < theta < pi).  
c
c     This transport equation can be rewritten (using definitions above):
c     
c     (3/2) d nT/ dt + (< |grad rho| >/A) d/drho A Q_itgc + ... = 0.
c
c     In this form, there is no ambiguity about which V' definition to use; 
c     one just calculates < |grad rho| > and the surface area, plugs in
c     the heat flux calculated from itgc, and that's that.  Note that 
c     there will be an additional factor of <|grad rho|> (typically 
c     stabilizing, I believe) that this definition of Q fails to capture.
c
c     The GA people tend (as I understand it) to lump the factor of 
c     <|grad rho|> out front together with the factor in the Q definition,
c     and thus talk about a <|grad rho|**2> effect on the anomalous transport
c     coefficients. 
c
c     Note that the definitions of the gradient scale lengths as 
c     
c     1/L_T == - d ln(T) / d ln(rho)
c
c     become functions of the choice of rho, the flux-surface label.  
c     If one is calculating everything within this ptor-itgc-eik code, 
c     there is no complication arising from this observation.  One may use
c     any definition available.  However, otherwise one must be careful to 
c     keep everything consistent.   Factors of drho_1/drho_2 appear...
c
c     Mike K. and I suggest that since we are performing local simulations
c     (in radius), the default definition of rho that we choose should be
c     a locally defined quantity.  Furthermore, it should be very fast to 
c     calculate (to keep the time spent in the eikcoefs subroutine to a 
c     minimum).  The definition which best seems to fit these two criteria 
c     
c     rho_d = d/D
c
c     where d == the diameter of the flux surface, and D == the diameter of
c     the last closed flux surface.
c
c     Any interpolation formulas or databases, etc., that we develop will 
c     probably use this definition of rho unless someone thinks of a good
c     reason to go with another choice.  
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,dz(lz)
      integer l,m,n
c
c This logic is wrong for input md=-1.     
c
      zero(1)=1.0
      do m=2,md
         zero(m)=0.5
      enddo

      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*gradrho(l)
         else
            dz(l)=2.*x0/ldb*jacobian(l)*gradrho(l)
         endif
      enddo

      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
               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
      return
      end