subroutine pcons
!
!     checks particle conservation 
!
  use itg_data
  implicit none

  complex, dimension(:,:,:), allocatable :: den, dent, denp
  complex, dimension(:), allocatable :: ntot, nt, np
  
  interface
     subroutine kapav(a, b)
       complex, dimension(:,:,:) :: a, b
     end subroutine kapav
  end interface
  
  real dtheta
  integer n,l,m

  allocate (den(kdpass, md, nd), dent(kdpass, md, nd), denp(kdpass, md, nd))
  allocate (ntot(nd), nt(nd), np(nd))

!
!     load trapped, passing, and total density 
!
  do n=1,nd
     do l=1,kd
        dent(l,1,n)=e_density(l,1,n)-phi_ba(l,1,n)
        denp(l,1,n)=0.
        den(l,1,n)=e_density(l,1,n)-phi_ba(l,1,n)
        do m=2,md
           dent(l,m,n)=0.
           denp(l,m,n)=0.
           den(l,m,n)=0.
        enddo
     enddo
     do l=kd+1,kdpass
        dent(l,1,n)=0.
        denp(l,1,n)=e_density(l,1,n)-phi_ba(l,1,n)
        den(l,1,n)=e_density(l,1,n)-phi_ba(l,1,n)
        do m=2,md
           dent(l,m,n)=0.
           denp(l,m,n)=0.
           den(l,m,n)=0.
        enddo
     enddo
  enddo
  
  call kapav(phi_ba,den)
!
!     flux surface average total electron density
!
  dtheta=2*x0/float(ldb)
  do n=1,nd
     ntot(n)=0.
     do l=1,ldb
        ntot(n)=ntot(n)+(potential(l,1,n)+e_denk(l,1,n))*(1+epse*cos(r(l)))*dtheta
     enddo
  enddo

  call kapav(phi_ba,dent)
!
!     flux surface average trapped electron density
!
  dtheta=2*x0/float(ldb)
  do n=1,nd
     nt(n)=0.
     do l=1,ldb
        nt(n)=nt(n)+(potential(l,1,n)*sqrt(epse*(1+cos(r(l)))/ &
             (1+epse*cos(r(l))))+e_denk(l,1,n))*(1+epse*cos(r(l)))*dtheta
     enddo
  enddo

  call kapav(phi_ba,denp)
!
!     flux surface average passing electron density
!
  dtheta=2*x0/float(ldb)
  do n=1,nd
     np(n)=0.
     do l=1,ldb
        np(n)=np(n)+(potential(l,1,n)*(1.-sqrt(epse*(1+cos(r(l)))/ &
             (1+epse*cos(r(l)))))+e_denk(l,1,n))*(1+epse*cos(r(l)))*dtheta
     enddo
  enddo
         
  do n=1,nd
     write(*,101) n,ntot(n),nt(n),np(n)
  enddo
101 format(' n,tot,trap,pass: ',i3,3(1pe14.6,1pe14.6,'; '))

!  check that flux surface average of <1/taub>_kappa = (1+epse)/2 
!      do l=1,kdpass
!         denp(l,1,1)=sqrt(2*epse/(1+epse))/taub(l)
!      enddo
!      call kapav(phi_ba,denp)
!      np(1)=0
!      do l=1,ldb
!         np(1)=np(1)+e_denk(l,1,1)*(1+epse*cos(r(l)))*dtheta
!      enddo
!      write(*,*) '<<1/taub>_kappa>_theta, (1+epse)/2: ',
!     .      np(1),(1+epse)/2
  

  deallocate (den, dent, denp, ntot, nt, np)

end subroutine pcons