subroutine poisson(phi,den,tp)

!     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
!     and Gregory W. Hammett. ALL RIGHTS RESERVED.
!
  use itg_data
  implicit none             
      
! define the arguments:
  complex, dimension(:,:,:,:) :: den, tp
  complex, dimension(:,:,:) :: phi
      
  interface 
     subroutine bcon(a)
       complex, dimension(:,:,:) :: a
     end subroutine bcon
  end interface
  interface
     subroutine kparfft(a, i)
       complex, dimension(:,:,:,:) :: a
       integer :: i
     end subroutine kparfft
  end interface

! declare local workspace
  complex, dimension(size(phi, 2), size(phi, 3)) :: phiav
  complex, dimension(size(den, 1), size(den, 2), size(den, 3)) :: nbar
  complex tmp,wgt_avg
  real, dimension(:), allocatable :: dz
  real znorm

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

!     
! load stuff for field line average     
!
  if (iperiod /= 2) then
     l_left=1
     l_right=ld
  endif

  allocate (dz(l_left:l_right))
  do l=l_left,l_right
! BD: backwards compatibility issue:
     if((igeo == 0) .and. (iphi00 == 0)) then
        dz(l)=2.*x0/ldb
     else
        dz(l)=2.*x0/ldb*jacobian(l)
     endif
  enddo
!
! Calculate the real-space ion density (summed over species)
!      
  nbar = 0.
  phiav = 0.

  do i=1,nspecies
     do l=1,ld-1
        nbar(l,:,:)=nbar(l,:,:)+(den(l,:,:,i)*nwgt(l,:,:,i) &
             +tp(l,:,:,i)*twgt(l,:,:,i))*g0wgtp(l,:,:,i)
     enddo
  enddo
 
! GWH:
! Parallel filtering to model finite-size particles.  The perp filter
! was included in nwgt and twgt in flrinit.f.  Formally, the parallel
! filtering should go in n_i_bar (as the particle charges are assigned
! to the grid), and not just as a filtering on phi.  In general geometry,
! with the <phi> calculation, there may be some subtle differences...
! However, I'm not sure which order the particle people actually do it
! in.  (If needed, could change this to be a filtering of phi at the end
! of this subroutine...)
!
!      if(s_par  /=  0.0) then
!         call bcon(nbar)        ! apply boundary conditions 
!         call kparfft(nbar,13)
!      endif

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


!      call baphi(phi)
!      call kapav(phi_ba,den)

  if (epse > 0.0) then
     phiav_old = phiav
  else
     phiav_old = 0.0
  endif
!
! GWH ?? 9/19/94: The use of phi_bk (and phiav_old) from the previous
! time step appears to be a kind of 1-iteration approximation.  We 
! should eventually test whether additional iterations on each call
! to poisson.f are required to get a converged solution to the
! integral equations implied by the bounce and kappa averages...
!
! I am confused about why the Poisson Eq. is so complicated now.
!
! Also, more thought needs to go into a general geometry formulation.
! This works for now, where  pfilter(l,1,n) is independent of l in 
! the large-aspect ratio limit...
!
! BD 6/22/95:
! Should be right for general geometry now.  pfilter doesn't pick up any
! new theta dependence.
! 
!     
! set Phi(ky=0, kz=0)=0 if iphi00=0     m=1, n /= 1 modes
!     
  
  if (iphi00 == 0 .and. nd /= 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
  deallocate (dz)
!     
! Do boundary conditions
!
  if(s_par  /=  0.0) then
! BD
! This needs to be fixed.  It was a bug, but probably didn't show up for nspecies = 1
!     call kparfft(phi,13)
     write(*,*) 'Error in bcon.  Check it out!'
  endif
  
  call bcon(phi)
  
  
end subroutine poisson