subroutine poisson(phi,den,tp)

c     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
c     and Gregory W. Hammett. ALL RIGHTS RESERVED.
c
      implicit none             
      
      include 'itg.par'
      include 'itg.cmn'

c define the arguments:
      complex den(lz,mz,nz,nspecz),phi(lz,mz,nz),tp(lz,mz,nz,nspecz)
      
c declare local workspace
      complex phiav(mz,nz),nbar(lz,mz,nz),tmp,wgt_avg,numer0
      complex e_denko(kz,mz,nz)
c      complex phiav_old(mz,nz)
      real dz(lz),znorm,delta
c      logical first
c      save phiav_old
c      integer mznz
c      parameter (mznz=mz*nz)
c      data phiav_old /mznz*(0.0,0.0)/

c declare other local variables:
      integer i,l,m,n
c     
c     start executable code:
c     
c iphi00=9  = set phi=0 for all Fourier components, not just (m=0,n=0).
c (no electric field, for neoclassical comparisons):
      if(iphi00 .eq. 9) return

c     
c load stuff for field line average     
c
      if (iperiod.ne.2) then
         l_left=1
         l_right=ld
      endif
c
      do l=l_left,l_right
c BD: backwards compatibility issue:
         if((igeo.eq.0) .and. (iphi00.eq.0)) then
            dz(l)=2.*x0/ldb
         else
            dz(l)=2.*x0/ldb*jacobian(l)
         endif
      enddo
c
c Calculate the real-space ion density (summed over species)
c      
      do n=1,nd
         do m=1,md
            do l=1,ld-1
               nbar(l,m,n)=0.
            enddo
            phiav(m,n)=0.
         enddo
      enddo
c
c     Switch outer loop for better multi-tasking:
c
      do n=1,nd
         do i=1,nspecies
            do m=1,md
               do l=1,ld-1
                  nbar(l,m,n)=nbar(l,m,n)
     .                 +(den(l,m,n,i)*nwgt(l,m,n,i)
     .                 +tp(l,m,n,i)*twgt(l,m,n,i))*g0wgtp(l,m,n,i)
               enddo
            enddo
         enddo
      enddo
c
c GWH:
c Parallel filtering to model finite-size particles.  The perp filter
c was included in nwgt and twgt in flrinit.f.  Formally, the parallel
c filtering should go in n_i_bar (as the particle charges are assigned
c to the grid), and not just as a filtering on phi.  In general geometry,
c with the <phi> calculation, there may be some subtle differences...
c However, I'm not sure which order the particle people actually do it
c in.  (If needed, could change this to be a filtering of phi at the end
c of this subroutine...)
c
c      if(s_par .ne. 0.0) then
c         call bcon(nbar)        ! apply boundary conditions 
c         call kparfft(nbar,13)
c      endif

c     
c If iphi00=2, require field line average of phi = 0 to ensure no
c radial electron transport:
c     
c     includes variation of k_x^2+k_y^2 along field line
c     
c
c Substantial simplification of this calculation can be had for 
c nonlinear runs, where we always know ahead of time that 
c     rky2(m,n)=0 iff m=1 and
c     rkx2(1,1,n).ne.0 for n>1.
c (This is true as long as shr.ne.0.)
c However, I am leaving the old logic for linear runs, in which the 
c layout in (m,n) is not known ahead of time. Dorland.
c 
cfpp$ skip R
      if (iphi00.eq.2 .and. (lin.eq.1 .or. shr.eq.0.)) then
         do n=1,nd
            do m=1,md
               wgt_avg=0.
               if ((rky2(m,n).eq.0.0).and.(rkx2(1,m,n).ne.0.0)) then
                  tmp=0.
                  do l=l_left,l_right-1
                     tmp=tmp+(nbar(l,m,n)-e_denk(l,1,n))
     .                    *pfilter(l,m,n)*jacobian(l)
                  enddo
                  wgt_avg=tmp/denom0(m,n)*(1.-stochf)
               endif
               phiav(m,n)=wgt_avg
            enddo
         enddo
      endif
cfpp$ noskip R
c
c Nonlinear replacement loop:
c
      if (iphi00.eq.2 .and. lin.eq.0) then
         do n=2,nd
            do l=l_left,l_right-1
               phiav(1,n)=phiav(1,n)+
     .              (nbar(l,1,n)-e_denk(l,1,n))*pfilter(l,1,n)
     .              /denom0(1,n)*(1.-stochf)*jacobian(l)
c
c  In adiabatic electron code, phiav was (T_i/T_e)<Phi>, now it is <Phi>
c
            enddo
         enddo
      endif
c     
c     Solve Poisson's equation
c  
c      do i=1,3  ! for iteration of Poisson's Eq.
c
      do n=1,nd
         do m=1,md
            do l=1,ld-1
               phi(l,m,n)=(nbar(l,m,n)-e_denk(l,m,n)+tiovte*(
     .          +phiav(m,n)+phi_bk(l,m,n)-phiav_old(m,n)))*pfilter(l,m,n)
c     don't use phiav, just phi_bk from previous timestep (doesn't work)
c               phi(l,m,n)=(nbar(l,m,n)-e_denk(l,m,n)+tiovte*(
c     .          +phi_bk(l,m,n)))*pfilter(l,m,n)
c
c     save e_denk
c               e_denko(l,m,n)=e_denk(l,m,n)
            enddo
c            write(*,*) 'phiav ',phiav(1,1)
c            write(*,*) 'phi   ',phi(3,1,1)
c            write(*,*) 'pfilter ',pfilter(3,1,1)
c            write(*,*) 
         enddo
      enddo

c Add external phi for Rosenbluth-Hinton test.
c (GWH 4/25/97):
      if(pmag .eq. 0.0) then
         do n=1,nd
            do l=1,ld
               phi(l,1,n)=phi(l,1,n)-ii/rkx(l,1,n)
            enddo
         enddo
      endif

c      call baphi(phi)
c      call kapav(phi_ba,den)
c      do n=1,nd
c         do m=1,md
c            do l=1,ld-1
c               e_denk(l,m,n)=e_denko(l,m,n)
c            enddo
c         enddo
c      enddo   
c      enddo

      if (epse.eq.0.0) then
         do m=1,md
            do n=1,nd
               phiav_old(m,n)=0.0
            enddo
         enddo
      else
         do m=1,md
            do n=1,nd
               phiav_old(m,n)=phiav(m,n)
            enddo
         enddo
      endif
c
c GWH ?? 9/19/94: The use of phi_bk (and phiav_old) from the previous
c time step appears to be a kind of 1-iteration approximation.  We 
c should eventually test whether additional iterations on each call
c to poisson.f are required to get a converged solution to the
c integral equations implied by the bounce and kappa averages...
c
c I am confused about why the Poisson Eq. is so complicated now.
c
c Also, more thought needs to go into a general geometry formulation.
c This works for now, where  pfilter(l,1,n) is independent of l in 
c the large-aspect ratio limit...
c
c BD 6/22/95:
c Should be right for general geometry now.  pfilter doesn't pick up any
c new theta dependence.
c 
c     
c set Phi(ky=0, kz=0)=0 if iphi00=0     m=1, n.ne.1 modes
c     
      if (iphi00.eq.0 .and. nd.ne.1) then
         do n=2,nd
            wgt_avg=0.0
            znorm=0.
            do l=l_left,l_right-1
               znorm=znorm+dz(l)
               wgt_avg=wgt_avg+phi(l,1,n)*dz(l)
            enddo
            do l=l_left,l_right-1
               phi(l,1,n)=phi(l,1,n)-wgt_avg/znorm
            enddo
         enddo
      endif
c     
c Do boundary conditions
c
      if(s_par .ne. 0.0) then
         call kparfft(phi,13)
      endif

      call bcon(phi)

      return
      end