subroutine pcons
c
c     checks particle conservation 
c
      include 'itg.par'
      include 'itg.cmn'
      complex den(kz,mz,nz),dent(kz,mz,nz),denp(kz,mz,nz)
      complex ntot(nz),nt(nz),np(nz)
      real dtheta
      integer n,l,m
c
c     load trapped, passing, and total density 
c
      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)
c
c     flux surface average total electron density
c
      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
c   
      call kapav(phi_ba,dent)
c
c     flux surface average trapped electron density
c
      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   
c
      call kapav(phi_ba,denp)
c
c     flux surface average passing electron density
c
      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,'; '))

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

      return
      end