module flr_terms

    use itg_data
    implicit none

!
! The following quantities include FLR corrections to Phi:
! phi_n is the FLR-corrected w_* term for the density equation
! phi_u is the FLR-corrected w_* term for the parallel velocity equation
! phi_tperp is the FLR-corrected w_* term for the perp. temperature equation
! phi_qperp is the FLR-corrected w_* term  for the par-perp heat equation
! phi_nd is the FLR-corrected w_d term for the density equation
! phi_tpard is the FLR-corrected w_d term for the perp. temperature equation
! phi_tperpd is the FLR-corrected w_d term  for the par temp. equation
! etc.
!
  complex, dimension(:,:,:,:), allocatable :: phi_n, phi_u, phi_tpar, phi_tperp, &
       phi_qperp, phi_rperp, phi_flr1, phi_flr2, phi_flr3, phi_nd, phi_tpard, &
       phi_tperpd, phi_qperpb

  real, dimension(:,:,:), allocatable :: pfilter2

contains

  subroutine flrinit(init_crit)
!     
!     initializes some arrays used in flr, poisson, timestep, and nlpsc
!       and kparfft
!     
!     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
!     and Gregory W. Hammett. ALL RIGHTS RESERVED.
!

    real rkp2, rgammas, tmu, wgt, zfilter
    real s18ceE, s18cfE
    integer i, j, l, m, n, ifail, lmax, init_crit

    if(.not.allocated(phi_n)) then
       allocate(phi_n(ld,md,nd,nspecies), phi_u(ld,md,nd,nspecies), &
            phi_tpar(ld,md,nd,nspecies), phi_tperp(ld,md,nd,nspecies), &
            phi_qperp(ld,md,nd,nspecies), phi_rperp(ld,md,nd,nspecies), &
            phi_flr1(ld,md,nd,nspecies), phi_flr2(ld,md,nd,nspecies), &
!            phi_flr3(ld,md,nd,nspecies), &
            phi_nd(ld,md,nd,nspecies), &
            phi_tpard(ld,md,nd,nspecies), phi_tperpd(ld,md,nd,nspecies), &
            phi_qperpb(ld,md,nd,nspecies))
       allocate (pfilter2(ld, md, nd))
    endif
!     
!     beginning of flr and poisson initializations
!     flr
!********************************************************************
!
    if(init_crit == 0) then
       do i=1,nspecies
          do n=1,nd
             do m=1,md
                do l=1,ld-1
                   gnwgt(l,m,n,i)=(g0wgt(l,m,n,i)+etak(m,n,i)*flrwgt(l,m,n,i))
                   gtpwgt(l,m,n,i)=etak_par(m,n,i)*g0wgt(l,m,n,i)
                   gtwgt(l,m,n,i)=etak(m,n,i)*(g0wgt(l,m,n,i)+flrwgt2(l,m,n,i))+flrwgt(l,m,n,i)  
!                   grwgt(l,m,n,i)=flrwgt(l,m,n,i)+etak(m,n,i)*(g0wgt(l,m,n,i) &
!                        +flrwgt2(l,m,n,i)+flrwgt(l,m,n,i))
                enddo
             enddo
          enddo
       enddo
       return
    endif

!********************************************************************
!     initialize some arrays only once per run:
!     
    ifail=0
    
    pfilter2=0.
!     
! BD: could reverse i,n loops for better multitasking performance
!
    do i=1,nspecies
       do n=1,nd
          do m=1,md
             do l=1,ld-1
                rkp2=rkperp2(l,m,n)*rho(i)**2
                rgammas=s18cfE(rkp2,ifail)/s18ceE(rkp2,ifail)
!
!     Use "new" Pade approximation
!
                if (iflr == 1) then
                   g0wgt(l,m,n,i)=1./(1.+rkp2/2.)
                   flrwgt(l,m,n,i)=-rkp2/(2.+rkp2)
                   flrwgt2(l,m,n,i)=-rkp2/(1.+rkp2/2.)**2
!     
!     Use "new" FLR model
!     
                else if(iflr == 2.or.iflr == 3.or.iflr == 8.or.iflr == 9) then                 
                   g0wgt(l,m,n,i)=sqrt(s18ceE(rkp2,ifail))
                   flrwgt(l,m,n,i)=-rkp2/2.*(1.-rgammas)
                   flrwgt2(l,m,n,i)=-rkp2/2.*(1.+(1.-rgammas)*(1.-rkp2/2.*(3.+rgammas)))
!     
!     Use exp(-b/2) model
!     
                else if(iflr == 4) then
                   g0wgt(l,m,n,i)=exp(-rkp2/2.)
                   flrwgt(l,m,n,i)=-rkp2/2.
                   flrwgt2(l,m,n,i)=-rkp2*(1.-rkp2/4.)
!
!     Try Taylor series
!     
                else if(iflr == 5) then
                   flrwgt(l,m,n,i)=-rkp2/2.
                   flrwgt2(l,m,n,i)=-rkp2/2.
                   g0wgt(l,m,n,i)=1.0
!
!     Like Ron, but Linsker treatment different
!     
                else if(iflr == 6) then
                   g0wgt(l,m,n,i)=s18ceE(rkp2,ifail)
                   flrwgt(l,m,n,i)=-rkp2*(1.-rgammas)
                   flrwgt2(l,m,n,i)=-rkp2*(2.-rgammas)+2.*rkp2**2*(1-rgammas)
!     
!     Full FLR, but no Linsker effect
!     
                else if(iflr == 7) then
                   g0wgt(l,m,n,i)=s18ceE(rkp2,ifail)
                   flrwgt(l,m,n,i)=-rkp2*(1.-rgammas)
                   flrwgt2(l,m,n,i)=-rkp2*(1.-rgammas)
                endif

                flrwgt3(l,m,n,i)=flrwgt(l,m,n,i)**2
                gnwgt(l,m,n,i)=(1.+etak(m,n,i)*flrwgt(l,m,n,i))*g0wgt(l,m,n,i)
                gtpwgt(l,m,n,i)=etak_par(m,n,i)*g0wgt(l,m,n,i)
                gtwgt(l,m,n,i)=etak(m,n,i)*(1.+flrwgt2(l,m,n,i))+flrwgt(l,m,n,i)
                gtwgt(l,m,n,i)=gtwgt(l,m,n,i)*g0wgt(l,m,n,i)
!                grwgt(l,m,n,i)=(flrwgt(l,m,n,i)+etak(m,n,i)*(1.+flrwgt2(l,m,n,i) &
!                     +flrwgt(l,m,n,i)))*g0wgt(l,m,n,i)
!
!     toroidal FLR terms, keep old forms for backward compatability
!
!     Pade in Poisson with new toroidal FLR terms,
!     suggested forms in M.A.Beer Ph.D. thesis
!
                if (iflr == 8) then
                   gnwgtd(l,m,n,i)=(2.+flrwgt(l,m,n,i))*g0wgt(l,m,n,i)/tau(i)
                   gtpwgtd(l,m,n,i)=2.*g0wgt(l,m,n,i)/tau(i)
                   gtwgtd(l,m,n,i)=(1.+2*flrwgt(l,m,n,i) &
                        +flrwgt2(l,m,n,i))*g0wgt(l,m,n,i)/tau(i)
                   gqwgtb(l,m,n,i)=(flrwgt2(l,m,n,i)-flrwgt(l,m,n,i)) &
                        *g0wgt(l,m,n,i)
!
!     Non-O(b) accurate closure in thesis
!
                else if(iflr == 9) then   
                   gnwgtd(l,m,n,i)=(2.+3.*flrwgt(l,m,n,i) &
                        -flrwgt2(l,m,n,i))*g0wgt(l,m,n,i)/tau(i)
                   gtpwgtd(l,m,n,i)=2.*g0wgt(l,m,n,i)/tau(i)
                   gtwgtd(l,m,n,i)=(1.+flrwgt(l,m,n,i) &
                        +4.*flrwgt2(l,m,n,i))*g0wgt(l,m,n,i)/tau(i)  
                   gqwgtb(l,m,n,i)=0.0    ! only tested this model with gpar B=0
!
!     Old forms of toroidal terms, for backward comp. only
!
                else
                   gnwgtd(l,m,n,i)=2.*g0wgt(l,m,n,i)/tau(i)
                   gtpwgtd(l,m,n,i)=2.*g0wgt(l,m,n,i)/tau(i)
                   gtwgtd(l,m,n,i)=(1.+flrwgt(l,m,n,i)*3.)*g0wgt(l,m,n,i)/tau(i)
                   gqwgtb(l,m,n,i)=0.0
                endif
                gnwgtd(l,m,n,i)=gnwgtd(l,m,n,i)*charge(i)
                gtpwgtd(l,m,n,i)=gtpwgtd(l,m,n,i)*charge(i)
                gtwgtd(l,m,n,i)=gtwgtd(l,m,n,i)*charge(i)
                flrwgt(l,m,n,i)=flrwgt(l,m,n,i)*g0wgt(l,m,n,i)
                flrwgt2(l,m,n,i)=flrwgt2(l,m,n,i)*g0wgt(l,m,n,i)
                flrwgt3(l,m,n,i)=flrwgt3(l,m,n,i)*g0wgt(l,m,n,i)
             enddo
          enddo
       enddo
    enddo
!
! end of initialization section
!********************************************************************
!     
! poisson
!********************************************************************
! initialize some arrays only once per run:
!
! precalculating the g0wgt and g2wgt arrays in FLR, and the
! pfilter array in POISSON, reduces the CPU time spent in these
! subroutines by about 25%.
!     
    do i=1,nspecies
       do n=1,nd
          do m=1,md
             do l=1,ld-1
                rkp2=rkperp2(l,m,n)*rho(i)**2
!     
!     Use "new" Pade approximation if iflr=1:
!     
                if (iflr == 1) then
                   wgt=(1.-1./(1.+rkp2))
                   g0wgtp(l,m,n,i)=2./(2.+rkp2)
                   gewgt(l,m,n,i)=2./(2.+rkp2)
                   twgt(l,m,n,i)=-rkp2/(2.+rkp2)
                   nwgt(l,m,n,i)=1.
!     
!     or use "new" FLR model for iflr=2 or 3:
!     
                else if(iflr == 2.or.iflr == 3) then
                   wgt=(1.-s18ceE(rkp2,ifail))
                   rgammas=s18cfE(rkp2,ifail)/s18ceE(rkp2,ifail)     
                   tmu=1.-rkp2*(1.-rkp2/2.)+rkp2/2.*rgammas*(1.-rkp2*rgammas)
                   twgt(l,m,n,i)=-rkp2/2.*(1.-rgammas)
                   nwgt(l,m,n,i)=1.-rkp2*(1.-rkp2/4.)+rkp2*rgammas/2.*(1.+rkp2*(1.-3.*rgammas/2.))
                   gewgt(l,m,n,i)=sqrt(s18ceE(rkp2,ifail))
                   g0wgtp(l,m,n,i)=gewgt(l,m,n,i)/tmu
!     
!     <J_0> = exp(-b/2) for iflr=4:
!     
                else if(iflr == 4) then
                   wgt=(1.-s18ceE(rkp2,ifail))
                   g0wgtp(l,m,n,i)=exp(-rkp2/2.)
                   gewgt(l,m,n,i)=exp(-rkp2/2.)
                   twgt(l,m,n,i)=-rkp2/2.
                   nwgt(l,m,n,i)=1.
!     
!     or use Taylor series:
!     
                else if(iflr == 5) then
                   wgt=rkp2
                   g0wgtp(l,m,n,i)=1.
                   gewgt(l,m,n,i)=1.
                   twgt(l,m,n,i)=-rkp2/2.
                   nwgt(l,m,n,i)=1.
!     
!     Like Ron, but Linsker treatment different
!     
                else if(iflr == 6) then
                   wgt=(1.-s18ceE(rkp2,ifail))
                   twgt(l,m,n,i)=0.0
                   nwgt(l,m,n,i)=1.
                   g0wgtp(l,m,n,i)=1.
!
!      Full FLR, but no Linsker effect
!     
                else if(iflr == 7) then
                   wgt=(1.-s18ceE(rkp2,ifail))
                   twgt(l,m,n,i)=0.0
                   nwgt(l,m,n,i)=1.
                   g0wgtp(l,m,n,i)=1.
!
!     Pade in Poisson with new toroidal FLR terms
!
                else if(iflr == 8.or.iflr == 9) then
                   wgt=(1.-s18ceE(rkp2,ifail))
                   g0wgtp(l,m,n,i)=1.
                   gewgt(l,m,n,i)=1.
                   twgt(l,m,n,i)=-rkp2/2./(1.+rkp2/2.)**2
                   nwgt(l,m,n,i)=1./(1.+rkp2/2.)
                endif
                
                g0wgtp(l,m,n,i)=g0wgtp(l,m,n,i)*charge(i)*n_I(i)
                pfilter2(l,m,n)=pfilter2(l,m,n)+wgt*n_I(i)*charge(i)**2/tau(i)
                pfilter3(l,m,n,i)=wgt*n_I(i)*charge(i)**2/tau(i)
! BD: Shouldn't these rky2 and rkx2 factors have an additional factor of 
! 1/B**2???
!
! GWH:
! Handle finite-size particle shapes/filtering. [Should check exactly 
! at what stage the particle simulations put in their phi filters, i.e.,
! before or after calculating <phi>?  Also, do they filter in real space
! or in field-line coordinates?  (Doesn't matter in the shearless 
! limit). Do different species have different sizes?]
!
! Use Dimit's style of filter (the parallel filter is elsewhere):
                zfilter=exp(-s_perp**4*rky2(m,n)**2)/(1+s_perp**4*rkx2(l,m,n)**2)
!
! BD: Question: shouldn't this rkperp2 definition have a factor of 1/B**2?
! ??
! Standard filter style:
!               rkperp2=rky2(m,n)+rkx2(l,m,n)         
!               zfilter=exp(-s_perp**2*rkperp2)
                twgt(l,m,n,i)=twgt(l,m,n,i)*zfilter
                nwgt(l,m,n,i)=nwgt(l,m,n,i)*zfilter

! Put particle shapes/filtering in the coefficients used to find n_i_bar:

             enddo
          enddo
       enddo
    enddo

    pfilter = 1.0/(tiovte+pfilter2)
    
!cfpp$ skip R  
    if(iphi00 == 2 .and. (lin == 1 .or. shr == 0.)) then
       do n=1,nd
          do m=1,md
             denom0(m,n)=0.               
             if ((rky2(m,n) == 0.0).and.(rkx2(1,m,n) /= 0.0)) then
                do l=l_left,l_right-1
                   denom0(m,n)=denom0(m,n)+pfilter(l,m,n)*pfilter2(l,m,n)*jacobian(l)
!   added dl/B factor for flux surface average
                enddo
             endif
          enddo
       enddo
    else
       do n=2,nd
          denom0(1,n)=0.
          do l=l_left,l_right-1
             denom0(1,n)=denom0(1,n)+pfilter(l,1,n)*pfilter2(l,1,n)*jacobian(l)
          enddo
       enddo
    endif
!cfpp$ noskip R
      
      
!     end of initialization section
!*********************************************************************
!     
!     for kparfft
!     
    do j=1,ndom
       lmax=ldb*scon(j)
       do l=1,lmax
          if (l <= lmax/2) then
             kpar(l,j)=(l-1)*2*pi/(x0*scon(j)*2.) 
          else
             kpar(l,j)=(l-1-lmax)*2*pi/(x0*scon(j)*2.) 
          endif
!     
!     local if qsf<0
!     
          if (qsf < 0 .and. igeo == 0) kpar(l,j)=shr
       enddo
    enddo

  end subroutine flrinit

  subroutine flr(phi)
!
! calculate the proper FLR-averaged potentials:
!
!     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
!     and Gregory W. Hammett. ALL RIGHTS RESERVED.
!
! Written by Greg Hammett 17-Dec-1990
! Modified by Bill Dorland 4-Mar-1991 
! and again by Bill 15-Sep-1991
! Now pseudo-spectral in radial direction.
! and again by Bill in Jan-1992
! added adiabatic electrons.
! and again by Bill in Mar-1992
! simplified burgeoning code and added nlpm terms.
! mbeer -- changed to Cowley coordinates, kx = s z ky  May 8, 1992
! and many changes since...
!

! define the arguments:
    complex, dimension(:,:,:) :: phi

! declare other local variables:
    integer :: m

    m = size(gnwgt, 4)
! start executable code:
  
    phi_n      = spread(phi,4,m) * gnwgt
    phi_u      = spread(phi,4,m) * g0wgt
    phi_tpar   = spread(phi,4,m) * gtpwgt
    phi_tperp  = spread(phi,4,m) * gtwgt
    phi_qperp  = spread(phi,4,m) * flrwgt
    phi_flr1   = spread(phi,4,m) * flrwgt
    phi_flr2   = spread(phi,4,m) * flrwgt2
    phi_nd     = spread(phi,4,m) * gnwgtd
    phi_tpard  = spread(phi,4,m) * gtpwgtd
    phi_tperpd = spread(phi,4,m) * gtwgtd
    phi_qperpb = spread(phi,4,m) * gqwgtb
!!    if(nperpmom == 3) phi_flr3 = spread(phi,4,m) * grwgt
  
  end subroutine flr

end module flr_terms