      subroutine timestep(
     *     density1,u_par1,t_par1,q_par1,t_perp1,q_perp1,
     *     density2,u_par2,t_par2,q_par2,t_perp2,q_perp2,
     *     e_density1,e_p1,e_r1,e_t1,
     *     e_density2,e_p2,e_r2,e_t2,dt)
c
c     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
c     and Gregory W. Hammett. ALL RIGHTS RESERVED.
c
c mbeer -- switched to Cowley coordinates, May 8, 1992
c
c Revision 3.0  1994/02/03  22:53:53  bdorland
c Upgraded to complex notation.
c
c  Advance equations one time step, not including collisional
c  dissipation
c
      implicit none
      include 'itg.par'
      include 'itg.cmn'
      include 'itg.mom'

      complex kp_density2(lz,mz,nz,nspecz)
      complex kp_upar2(lz,mz,nz,nspecz)
      complex kp_tpar2(lz,mz,nz,nspecz)
      complex kp_qpar2(lz,mz,nz,nspecz)
      complex kp_tperp2(lz,mz,nz,nspecz)
      complex kp_qperp2(lz,mz,nz,nspecz)
      complex kp_vphi(lz,mz,nz,nspecz)
c     complex akp_tpar2(lz,mz,nz,nspecz)
      complex akp_qpar2(lz,mz,nz,nspecz)
      complex akp_qperp2(lz,mz,nz,nspecz)
      complex wd_t(kz,lz,nz),awd_t(kz,lz,nz)

      complex density2(lz,mz,nz,nspecz),u_par2(lz,mz,nz,nspecz),
     * t_par2(lz,mz,nz,nspecz),q_par2(lz,mz,nz,nspecz),
     * t_perp2(lz,mz,nz,nspecz),q_perp2(lz,mz,nz,nspecz)
      complex e_density2(kz,mz,nz),e_p2(kz,mz,nz),e_r2(kz,mz,nz)
      complex e_t2(kz,mz,nz)

      real nuur,nuui,vmom,dt(mz),nu_erk,nuky
      real switch(mz)
      integer i,l,m,n
c     
c     turn off omegad dissipative terms for m=0
      if (iflow.eq.1) then
         switch(1)=1.0
      else
         switch(1)=0.0
      endif
      do m=2,md
         switch(m)=1.0
      enddo

      if (nparmom.eq.3) then
         nuur=nu5r
         nuui=nu5i
         vmom=0.
      else
         nuur=0.0
         nuui=0.0
         vmom=1.
      endif   
c
      do n=1,nd 
         do i=1,nspecies
            do m=1,md
c     
c     load kp_upar2,kp_vphi,etc.
c     
               do l=1,ld-1
                  kp_density2(l,m,n,i)=density2(l,m,n,i)*bmaginv(l)
                  kp_upar2(l,m,n,i)=u_par2(l,m,n,i)*bmaginv(l)
                  kp_tpar2(l,m,n,i)=t_par2(l,m,n,i)*bmaginv(l)
                  kp_vphi(l,m,n,i)=phi_u(l,m,n,i)
               enddo
               kp_density2(ld,m,n,i)=0.
               kp_upar2(ld,m,n,i)=0.
               kp_tpar2(ld,m,n,i)=0.
               kp_vphi(ld,m,n,i)=0.
            enddo
         enddo 
      enddo
c
      call kparfft(kp_density2,0)
      call kparfft(kp_upar2,0)
      call kparfft(kp_tpar2,0)
      call kparfft(kp_vphi,0)

      do n=1,nd
       do i=1,nspecies
        do m=1,md
         do l=1,ld
c
c -----------------density equation----------------------
c
      density1(l,m,n,i)=density(l,m,n,i)-dt(m)*(
     .    nl_density1(l,m,n,i)+nl_density2(l,m,n,i)
     .   +ii*rky(m,n)/Ln(i)*phi_n(l,m,n,i)
     .   +ii*kp_upar2(l,m,n,i)*b_mag(l)*vt(i)
     .   -ii*omegad(l,m,n,i)*(phi_nd(l,m,n,i)+2.*density2(l,m,n,i))
     .   -ii*omegad(l,m,n,i)*(t_par2(l,m,n,i)+t_perp2(l,m,n,i)))
c
c -------------parallel momentum equation----------------
c
      u_par1(l,m,n,i)=u_par(l,m,n,i)-dt(m)*(
     .    nl_upar1(l,m,n,i)+nl_upar2(l,m,n,i)
     .   +ii*kp_vphi(l,m,n,i)*vt(i)*charge(i)/tau(i)
     .   +ii*(kp_tpar2(l,m,n,i)+kp_density2(l,m,n,i))*b_mag(l)*vt(i)
     .   -ii*omegad(l,m,n,i)*2.*((2.+nuui*switch(m))*u_par2(l,m,n,i)
     .   +0.5*vmom*(q_par2(l,m,n,i)+q_perp2(l,m,n,i)))
     .   +2.*nuur*abs(omegad(l,m,n,i)*switch(m))*u_par2(l,m,n,i)
     .   +(t_perp2(l,m,n,i)+density2(l,m,n,i)+phi_qperp(l,m,n,i)
     .   *charge(i)/tau(i))*bgrad(l)*vt(i))
         enddo
        enddo
       enddo
      enddo
c
c     ---------three pole closure: parallel moments----------
c     
      if(nparmom.eq.3) then
         do n=1,nd
            do i=1,nspecies
               do m=1,md
                  do l=1,ld
                     q_par2(l,m,n,i)=rkkin*t_par2(l,m,n,i)
                  enddo
               enddo
            enddo  
         enddo
      endif
c
c     load kp_qpar2
c
      do n=1,nd
         do i=1,nspecies
            do m=1,md
               do l=1,ld-1
                  kp_qpar2(l,m,n,i)=q_par2(l,m,n,i)*bmaginv(l)
               enddo
               kp_qpar2(ld,m,n,i)=0.0
            enddo 
         enddo
      enddo

      if (nparmom.eq.3) then
         call kparfft(kp_qpar2,2)
      else
         call kparfft(kp_qpar2,0)
      endif
c     
      do n=1,nd
       do i=1,nspecies
        do m=1,md
         do l=1,ld
c
c -------------parallel temperature equation-----------
c
         t_par1(l,m,n,i)=t_par(l,m,n,i)-dt(m)*(
     .     nl_tpar(l,m,n,i)
     .    +2.*ii*kp_upar2(l,m,n,i)*b_mag(l)*vt(i)
     .    +ii*rky(m,n)/Ln(i)*phi_tpar(l,m,n,i)
     .    +ii*kp_qpar2(l,m,n,i)*b_mag(l)*vt(i)
     .    -ii*omegad(l,m,n,i)*2.*
     .      ( 1.5*(density2(l,m,n,i)+t_par2(l,m,n,i))
     .      +(1.5+nu1i*switch(m))*t_par2(l,m,n,i)
     .      +(0.5+nu2i*switch(m))*t_perp2(l,m,n,i)
     .      -0.5*(t_perp2(l,m,n,i)+density2(l,m,n,i))
     .      +0.5*phi_tpard(l,m,n,i) )
     .    +2.*nu1r*abs(omegad(l,m,n,i)*switch(m))*t_par2(l,m,n,i) 
     .    +2.*nu2r*abs(omegad(l,m,n,i)*switch(m))*t_perp2(l,m,n,i)
     .    +2.*(u_par2(l,m,n,i)+q_perp2(l,m,n,i))*bgrad(l)*vt(i) 
     .    +2./3.*nu_ii(i)*(t_par2(l,m,n,i)-t_perp2(l,m,n,i)) )
           enddo
         enddo
       enddo
      enddo
c
      if(nparmom.eq.4) then
c
c      load kp_tpar2,akp_qpar2
c
         do n=1,nd
            do i=1,nspecies
               do m=1,md
                  do l=1,ld-1
                     kp_tpar2(l,m,n,i)=t_par2(l,m,n,i)
c BD: not sure which we should be using; I favor commented line
c                     akp_qpar2(l,m,n,i)=q_par2(l,m,n,i)*bmaginv(l)
                     akp_qpar2(l,m,n,i)=q_par2(l,m,n,i)
                  enddo
                  kp_tpar2(ld,m,n,i)=0.0
                  akp_qpar2(ld,m,n,i)=0.0
               enddo 
            enddo
         enddo
         call kparfft(kp_tpar2,0)
         call kparfft(akp_qpar2,1)
c
         do n=1,nd
          do i=1,nspecies
           do m=1,md
            do l=1,ld
c
c --------------parallel heat flux equation---------------
c
             q_par1(l,m,n,i)=q_par(l,m,n,i)-dt(m)*(
     .         nl_qpar(l,m,n,i)
c BD: This looks like an error to me: 
     .        +zBq1*ii*kp_tpar2(l,m,n,i)*vt(i)
c BD: I favor the commented lines:
c     .        +(3-zBq1)*t_par2(l,m,n,i)*bgrad(l)*vt(i) 
c     .        +zDq1*akp_qpar2(l,m,n,i)*vt(i)*b_mag(l)
     .        +zDq1*akp_qpar2(l,m,n,i)*vt(i)
     .        -ii*omegad(l,m,n,i)*
     .         ( (nu5i*switch(m)+6.)*u_par2(l,m,n,i)
     .         + (nu6i*switch(m)-3.)*q_par2(l,m,n,i)
     .         + (nu7i*switch(m)-3.)*q_perp2(l,m,n,i) )
     .        +abs(omegad(l,m,n,i)*switch(m))*
     .         ( nu5r*u_par2(l,m,n,i)
     .         + nu6r*q_par2(l,m,n,i)
     .         + nu7r*q_perp2(l,m,n,i) )
     .        +nu_ii(i)*q_par2(l,m,n,i) )
            enddo
           enddo
          enddo
         enddo
c
      endif
c
c ----------one pole closure: perpendicular moments------------
c
      if(nperpmom.eq.1) then
         do n=1,nd
            do i=1,nspecies
               do m=1,md
                  do l=1,ld
                     q_perp2(l,m,n,i)=
     .                    rkkperp*(t_perp2(l,m,n,i)
     .                    +charge(i)/tau(i)*phi_qperp(l,m,n,i))
                  enddo
               enddo
            enddo
         enddo
      endif
c
c      load kp_qperp2
c
      do n=1,nd
         do i=1,nspecies
            do m=1,md
               do l=1,ld-1
                  kp_qperp2(l,m,n,i)=q_perp2(l,m,n,i)*bmaginv(l)**2
               enddo
               kp_qperp2(ld,m,n,i)=0.0
            enddo 
         enddo
      enddo
      if (nperpmom.eq.1) then
         call kparfft(kp_qperp2,2)
      else
         call kparfft(kp_qperp2,0)
      endif
c
      do n=1,nd
       do i=1,nspecies
        do m=1,md
         do l=1,ld
c
c ----------perpendicular temperature equation-------------- 
c
          t_perp1(l,m,n,i)=t_perp(l,m,n,i)-dt(m)*(
     .      nl_tperp1(l,m,n,i)+nl_tperp2(l,m,n,i)+nl_tperp3(l,m,n,i)
     .     +ii*rky(m,n)/Ln(i)*phi_tperp(l,m,n,i)
     .     +ii*kp_qperp2(l,m,n,i)*b_mag(l)**2*vt(i)
     .     -ii*2.*omegad(l,m,n,i)*
     .      ( t_perp2(l,m,n,i)+density2(l,m,n,i)
     .      +(0.5+nu3i*switch(m))*t_par2(l,m,n,i)
     .      +(1.+nu4i*switch(m))*t_perp2(l,m,n,i)
     .      -0.5*(density2(l,m,n,i)+t_par2(l,m,n,i))
     .      +0.5*phi_tperpd(l,m,n,i) )
     .     +2.*nu4r*abs(omegad(l,m,n,i)*switch(m))*t_perp2(l,m,n,i) 
     .     +2.*nu3r*abs(omegad(l,m,n,i)*switch(m))*t_par2(l,m,n,i)
     .     -u_par2(l,m,n,i)*bgrad(l)*vt(i)
     .     +nu_ii(i)/3.*(t_perp2(l,m,n,i)-t_par2(l,m,n,i))  )
          enddo
         enddo
        enddo
      enddo
c
      if(nperpmom.ge.2) then
c
c     load kp_tperp2,akp_qperp2
c     
         do n=1,nd
            do i=1,nspecies
               do m=1,md
                  do l=1,ld-1
c BD: This looks more "correct", but doesn't work as well: 
c                     akp_qperp2(l,m,n,i)=
c     >                 q_perp2(l,m,n,i)*rkkperp2*bmaginv(l)
                     akp_qperp2(l,m,n,i)=
     >                 q_perp2(l,m,n,i)*rkkperp2
                     kp_tperp2(l,m,n,i)=t_perp2(l,m,n,i)
     >                 +charge(i)/tau(i)*phi_qperp(l,m,n,i)
                  enddo
                  akp_qperp2(ld,m,n,i)=0.0
                  kp_tperp2(ld,m,n,i)=0.0
               enddo
            enddo  
         enddo
         call kparfft(akp_qperp2,1)
         call kparfft(kp_tperp2,0)
c
         do n=1,nd
          do i=1,nspecies
           do m=1,md
            do l=1,ld
c
c -----------perpendicular heat flux equation----------------
c
             q_perp1(l,m,n,i)=q_perp(l,m,n,i)-dt(m)*(
     .         nl_qperp1(l,m,n,i)+nl_qperp2(l,m,n,i)
     .        +nl_qperp3(l,m,n,i)
     .        +ii*kp_tperp2(l,m,n,i)*vt(i)
c BD: This looks more "correct", but doesn't work as well: 
c     .        +akp_qperp2(l,m,n,i)*vt(i)*b_mag(l)
     .        +akp_qperp2(l,m,n,i)*vt(i)
     .        -ii*omegad(l,m,n,i)*
     .         ( (nu8i*switch(m)+1.)*u_par2(l,m,n,i)
     .         + (nu9i*switch(m)-1.)*q_par2(l,m,n,i)
     .         + (nu10i*switch(m)-1.)*q_perp2(l,m,n,i) )
     .        +abs(omegad(l,m,n,i)*switch(m))*
     .         ( nu8r*u_par2(l,m,n,i)
     .         + nu9r*q_par2(l,m,n,i)
     .         + nu10r*q_perp2(l,m,n,i) )
     .         +(t_perp2(l,m,n,i)-t_par2(l,m,n,i)+phi_qperpb(l,m,n,i)
     .         *charge(i)/tau(i))*bgrad(l)*vt(i)
     .         +nu_ii(i)*q_perp2(l,m,n,i) )
            enddo
           enddo
          enddo
         enddo
      endif
c
c bounce averaged electrons:
c
      if (epse.gt.0.0) then
      if (nemom.eq.3) then
        do n=1,nd
          do m=1,md
            do l=1,kd
              wd_t(l,m,n)=nuai*e_density2(l,m,n)+nubi*e_p2(l,m,n)
     .            +nuci*e_r2(l,m,n)
              awd_t(l,m,n)=nuar*e_density2(l,m,n)+nubr*e_p2(l,m,n)
     .            +nucr*e_r2(l,m,n)

            enddo
          enddo
        enddo  
      else
        do n=1,nd
          do m=1,md
            do l=1,kd
              wd_t(l,m,n)=e_t2(l,m,n)
              awd_t(l,m,n)=0.
            enddo
          enddo
        enddo  
      endif  

c     nuky=rmu1 for hyperviscosity in electron eqns
      nuky=0
c      if (ncycle.ge.4000) nuky=nuky/10.
         
      do n=1,nd
        do m=1,md
          do l=1,kd

c Try Rewoldt's version of a Krook model collision operator, as given
c in Eq. 4 and 6 of Kotschenreuther, Rewoldt, Tang, "Comparison of Initial
c Value and Eigenvalue Codes for Kinetic Toroidal Plasma Instabilities",
c PPPL-2986 (1994) (to be publ. in Comp. Phys. Comm.?).
c
c One of the interesting features of this model is that the loss rate 
c goes to infinity near the trapped-passing boundary, so we will treat 
c this term implicitly below.  For simplicity for now, we will ignore
c the energy dependence of the collision frequency so the collision
c operator is almost diagonal (except for the phi_ba term...).
c (Also, take the delta=0 limit of Rewoldt's full formula.)
c
c nueeff = nu_e_effective = nu_e / (r/R)
c
c             nu_erk=nueeff/4.0*1.31/(1-kap(l)**2)**2
c For now, turn off the Krook model, and try out the full Lorentz
c collisions:
c
c add damping proportional to ky^4 to electron moments
             nu_erk=nuky*rky2(m,n)**2
c
c GWH 9/12/94, it would seem to me that a better formula would be:
c     nu_erk=nu_e/2.0/(1-kap(l))**2
c based on Dimensional arguments and D=dx**2/dt.  This would also
c give more naturally a sqrt(nu) dependence at low collisionality.
c
c
c -------------electron density equation-----------------
c
             e_density1(l,m,n)=e_density(l,m,n)-dt(m)*(
     .            nl_edensity(l,m,n)+
     .            ii*omegade(l,m,n)*1.5*
     .            (e_p2(l,m,n)/tiovte-phi_ba(l,m,n))
     .            +ii*rky(m,n)*phi_ba(l,m,n)
     .            -nu_erk*phi_ba(l,m,n))
             e_density1(l,m,n)=e_density1(l,m,n)/(1+dt(m)*nu_erk)
c
c -------------electron pressure equation--------------
c
             e_p1(l,m,n)=e_p(l,m,n)-dt(m)*(
     .            nl_ep(l,m,n)+
     .            ii*omegade(l,m,n)*2.5*
     .            (e_r2(l,m,n)/tiovte-phi_ba(l,m,n))
     .            +ii*rky(m,n)*phi_ba(l,m,n)*(1.+etae)
     .            -nu_erk*phi_ba(l,m,n))
             e_p1(l,m,n)=e_p1(l,m,n)/(1+dt(m)*nu_erk)
c
c  -------------electron r equation--------------
c
             e_r1(l,m,n)=e_r(l,m,n)-dt(m)*(
     .            nl_er(l,m,n)+
     .            ii*omegade(l,m,n)*3.5*
     .            (wd_t(l,m,n)/tiovte-phi_ba(l,m,n))
     .            +ii*rky(m,n)*phi_ba(l,m,n)*(1.+2.*etae)
     .            +3.5*abs(omegade(l,m,n))*awd_t(l,m,n)/tiovte 
     .            -nu_erk*phi_ba(l,m,n))
             e_r1(l,m,n)=e_r1(l,m,n)/(1+dt(m)*nu_erk)

          enddo
        enddo
      enddo
c
c  -------------electron t equation--------------
c
      if (nemom.eq.4) then
        do n=1,nd
          do m=1,md
            do l=1,kd 
              e_t1(l,m,n)=e_t(l,m,n)-dt(m)*(
     .            nl_et(l,m,n)+
     .            ii*omegade(l,m,n)*4.5*
     .            ((nuai*e_density2(l,m,n)+nubi*e_p2(l,m,n)
     .            +nuci*e_r2(l,m,n)+nudi*e_t2(l,m,n))/tiovte
     .            -phi_ba(l,m,n))
     .            +ii*rky(m,n)*phi_ba(l,m,n)*(1.+3.*etae)
     .            +4.5*abs(omegade(l,m,n))*
     .            (nuar*e_density2(l,m,n)+nubr*e_p2(l,m,n)
     .            +nucr*e_r2(l,m,n)+nudr*e_t2(l,m,n))/tiovte)
              e_t1(l,m,n)=e_t1(l,m,n)/(1+dt(m)*nu_erk)
            enddo
          enddo
        enddo
      endif  

c Evaluate Collisions on trapped electrons:
      if(nueeff .ne. 0.0) call ecollis(e_density1,e_p1,e_r1,e_t1,dt)
      endif

c     
c     Make like ron's wedge
c     
      goto 33
      do n=1,nd
         do i=1,nspecies
            do m=1,md
               if ((abs(r0(m,n)).gt.pi+.001).or.(m.eq.1)) then
                  do l=1,ld
                     density1(l,m,n,i)=cmplx(1.e-15,1.e-15)
                     u_par1(l,m,n,i)=cmplx(1.e-15,1.e-15)
                     t_par1(l,m,n,i)=cmplx(1.e-15,1.e-15)
                     q_par1(l,m,n,i)=cmplx(1.e-15,1.e-15)
                     t_perp1(l,m,n,i)=cmplx(1.e-15,1.e-15)
                     q_perp1(l,m,n,i)=cmplx(1.e-15,1.e-15)
                  enddo
               endif
            enddo
         enddo   
      enddo
 33   continue

      do i=1,nspecies
         call bcon(density1(1,1,1,i))
         call bcon(u_par1(1,1,1,i))
         call bcon(t_par1(1,1,1,i))
         if(nparmom.ge.4) call bcon(q_par1(1,1,1,i))
         call bcon(t_perp1(1,1,1,i))
         if(nperpmom.ge.2) call bcon(q_perp1(1,1,1,i))
      enddo

      return
      end

