subroutine ecollis(e_density1,e_p1,e_r1,e_t1,dt)
c
c bounce-averaged Lorentz collision operator for trapped electrons
c gwh 9/11/94
c
c Use implicit time-advancement, for stability and robustness,
c employing standard tri-diagonal techniques.
c
c
      implicit none
      include 'itg.par'
      include 'itg.cmn'
      include 'itg.mom'

      real dt(mz)

c Local Variables:
      real dtnu0(mz),epsb
      integer l,m,n

c************************************************************************

c For simplicity for now, we will ignore the energy dependence of the 
c collision frequency so the collision operator is almost diagonal 
c (except for the phi_ba term, which is treated explicitly, which 
c appears to work fine, particularly given the way Mike set up the 
c Poisson equation solver in its depedence on <<phi>_b>_kappa from 
c the past time step).
c
c Use the definitions:
c
c nueeff= nu_e_eff L_n/v_t, effective electron collision frequency
c MB: clarification: nueeff= nu_e_eff L_ne/v_ti
c
c nu_e_eff = nu_e/(r/R)
c          = nustare * (1+1/Zeff) *3 *sqrt(pi*r/R) *Ln/R/q *sqrt(mi/me)
c MB:  slight error in previous line, it's really nueeff, and missing Ti/Te
c nueeff   = nustare * (1+1/Zeff) *3*sqrt(pi*r/R/2)*Lne/R/q *sqrt(mi*Te/me*Ti)
c where nustare is the usual neoclassical/snap collision
c frequency parameter (ignoring e-e collisions),
c 	nustare = nu_ei*Zeff/(r/R)/omega_b
c 	nu_ei = 1/tau_e = Braginskii's definition
c 	omega_b = sqrt(2*r/R*Te/me)/(q R) = bounce frequency
c MB: TRANSP,SNAP use: omega_b = sqrt(r/R*Te/me)/(q R),
c     nueeff above reflects this sqrt(2) change
c
      epsb=epse/(1+epse)
      do m=1,md
         dtnu0(m) = dt(m)*nueeff*epse/(8*epsb**2)*sqrt(2*pi)/16
      enddo
c
c      if (nemom.eq.3) then

      if(mrr(1,1) .ne. 0) call aborter(6,
     .     'ecollis.f works only if the first m mode is m=0')

c solve for all (k_theta .ne. 0) modes:

      do n=1,nd
        do m=2,md
          do l=2,kd-1

c
c -------------electron density equation-----------------
c
             e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte
     .            *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))
     .            -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

c
c -------------electron pressure equation--------------
c
             e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte
     .            *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))
     .            -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))
c
c  -------------electron r equation--------------
c
             e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte
     .            *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))
     .            -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

          enddo

c handle the boundary at l=kd (everything vanishes for l>kd):
          l=kd
          e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte
     .         *jac(l)*(diff(l)*(               -phi_ba(l,m,n))
     .         -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

          e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte
     .         *jac(l)*(diff(l)*(               -phi_ba(l,m,n))
     .         -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

          e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte
     .         *jac(l)*(diff(l)*(               -phi_ba(l,m,n))
     .         -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

c handle the boundary at the origin at l=1 (no flux across kappa=0
c because of symmetry...):
          l=1
          e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte
     .         *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)))

          e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte
     .         *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)))

          e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte
     .         *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)))

        enddo
      enddo

c
c  -------------electron t equation--------------
c
c Collision operator hasn't been added to the t-equation yet...

      if (nemom.eq.4) then
         call aborter(6,"Collision operator not yet valid for nemom=4")
      endif  


c***********************************************************
c solve for all (k_theta = 0) modes:
c (most of the below is a slightly modified copy of the section
c above for k_theta .ne. 0 modes)

      do n=1,nd
         m=1
          do l=2,kdpass-1

c
c -------------electron density equation-----------------
c
             e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte
     .            *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))
     .            -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

c
c -------------electron pressure equation--------------
c
             e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte
     .            *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))
     .            -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))
c
c  -------------electron r equation--------------
c
             e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte
     .            *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n))
     .            -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

          enddo

c handle the boundary at l=kdpass (no flux beyond max kappa):
c (These expressions could be simplified by combining into the above
c do loops if phi_ba(kdpass+1) was defined, since diff(kdpass)=0):

          l=kdpass
          e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte
     .         *jac(l)*(diff(l)*(               -phi_ba(l,m,n))
     .         -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

          e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte
     .         *jac(l)*(diff(l)*(               -phi_ba(l,m,n))
     .         -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

          e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte
     .         *jac(l)*(diff(l)*(               -phi_ba(l,m,n))
     .         -diff(l-1)*(phi_ba(l,m,n)-phi_ba(l-1,m,n)))

c handle the boundary at the origin at l=1 (no flux across kappa=0
c because of symmetry...):
          l=1
          e_density1(l,m,n)=e_density1(l,m,n)-dtnu0(m)*tiovte
     .         *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)))

          e_p1(l,m,n)=e_p1(l,m,n)-dtnu0(m)/2*tiovte
     .         *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)))

          e_r1(l,m,n)=e_r1(l,m,n)-dtnu0(m)*0.3125*tiovte
     .         *jac(l)*(diff(l)*(phi_ba(l+1,m,n)-phi_ba(l,m,n)))

      enddo

c**********************************************************************
c
c Set up coefficients of the tridiagonal equation to be solved,
c of the form:
c
c a(l)*den(l-1) + b(l)*den(l) + c(l)*den(l+1) = den_old(l)
c

c
c For dt(k_theta) this gets a little more complicated:
c
c     density moment
c
      do m=1,md
         do l=1,kdpass
            aec(l,m)=-dtnu0(m)*jac(l)*diff(l-1)
            cec(l,m)=-dtnu0(m)*jac(l)*diff(l)
            bec(l,m)=1-aec(l,m)-cec(l,m)
         enddo
         aec(1,m)=0.0
         cec(kdpass,m)=0.0
      enddo
c set up coefficients of the inversion formula:
      do m=1,md
         zbeta(1,m)=bec(1,m)
         do l=2,kdpass
            zgamma(l,m)=cec(l-1,m)/zbeta(l-1,m)
            zbeta(l,m)=bec(l,m)-aec(l,m)*zgamma(l,m)
         enddo
      enddo
c recursive-solution to tridagonal equations doesn't vectorize so
c put it as the outer loop.  (When optimizing for multitasking,
c may want to multitask over m (and over different fluid moments), 
c and vectorize over n.)
c
c Try to speed up by combining the various loops...

      do m=1,md
         do n=1,nd
            e_density1(1,m,n)= e_density1(1,m,n)/zbeta(1,m)
         enddo
      enddo

c solve for (k_theta .ne. 0) modes:
      do l=2,kd
       do m=2,md
        do n=1,nd
           e_density1(l,m,n)=(e_density1(l,m,n)
     .          -aec(l,m)*e_density1(l-1,m,n))/zbeta(l,m)
        enddo
       enddo
      enddo

      do l=kd-1,1,-1
       do m=2,md
        do n=1,nd
           e_density1(l,m,n)=e_density1(l,m,n)
     .          -zgamma(l+1,m)*e_density1(l+1,m,n)
         enddo
       enddo
      enddo

c solve for (k_theta .eq. 0) modes:
      do l=2,kdpass
       m=1
        do n=1,nd
           e_density1(l,m,n)=(e_density1(l,m,n)
     .          -aec(l,m)*e_density1(l-1,m,n))/zbeta(l,m)
        enddo
      enddo

      do l=kdpass-1,1,-1
       m=1
        do n=1,nd
           e_density1(l,m,n)=e_density1(l,m,n)
     .          -zgamma(l+1,m)*e_density1(l+1,m,n)
         enddo
      enddo
c
c     p moment
c
      do m=1,md
         do l=1,kdpass
            aec(l,m)=-dtnu0(m)/2*jac(l)*diff(l-1)
            cec(l,m)=-dtnu0(m)/2*jac(l)*diff(l)
            bec(l,m)=1-aec(l,m)-cec(l,m)
         enddo
         aec(1,m)=0.0
         cec(kdpass,m)=0.0
      enddo
c set up coefficients of the inversion formula:
      do m=1,md
         zbeta(1,m)=bec(1,m)
         do l=2,kdpass
            zgamma(l,m)=cec(l-1,m)/zbeta(l-1,m)
            zbeta(l,m)=bec(l,m)-aec(l,m)*zgamma(l,m)
         enddo
      enddo
c recursive-solution to tridagonal equations doesn't vectorize so
c put it as the outer loop.  (When optimizing for multitasking,
c may want to multitask over m (and over different fluid moments), 
c and vectorize over n.)
c
c Try to speed up by combining the various loops...

      do m=1,md
         do n=1,nd
            e_p1(1,m,n)= e_p1(1,m,n)/zbeta(1,m)
         enddo
      enddo

c solve for (k_theta .ne. 0) modes:
      do l=2,kd
       do m=2,md
        do n=1,nd
           e_p1(l,m,n)=(e_p1(l,m,n)
     .          -aec(l,m)*e_p1(l-1,m,n))/zbeta(l,m)
        enddo
       enddo
      enddo

      do l=kd-1,1,-1
       do m=2,md
        do n=1,nd
           e_p1(l,m,n)=e_p1(l,m,n)
     .          -zgamma(l+1,m)*e_p1(l+1,m,n)
         enddo
       enddo
      enddo

c solve for (k_theta .eq. 0) modes:
      do l=2,kdpass
       m=1
        do n=1,nd
           e_p1(l,m,n)=(e_p1(l,m,n)
     .          -aec(l,m)*e_p1(l-1,m,n))/zbeta(l,m)
        enddo
      enddo

      do l=kdpass-1,1,-1
       m=1
        do n=1,nd
           e_p1(l,m,n)=e_p1(l,m,n)
     .          -zgamma(l+1,m)*e_p1(l+1,m,n)
         enddo
      enddo
c
c     r moment
c
      do m=1,md
         do l=1,kdpass
            aec(l,m)=-dtnu0(m)*0.3125*jac(l)*diff(l-1)
            cec(l,m)=-dtnu0(m)*0.3125*jac(l)*diff(l)
            bec(l,m)=1-aec(l,m)-cec(l,m)
         enddo
         aec(1,m)=0.0
         cec(kdpass,m)=0.0
      enddo
c set up coefficients of the inversion formula:
      do m=1,md
         zbeta(1,m)=bec(1,m)
         do l=2,kdpass
            zgamma(l,m)=cec(l-1,m)/zbeta(l-1,m)
            zbeta(l,m)=bec(l,m)-aec(l,m)*zgamma(l,m)
         enddo
      enddo
c recursive-solution to tridagonal equations doesn't vectorize so
c put it as the outer loop.  (When optimizing for multitasking,
c may want to multitask over m (and over different fluid moments), 
c and vectorize over n.)
c
c Try to speed up by combining the various loops...

      do m=1,md
         do n=1,nd
            e_r1(1,m,n)= e_r1(1,m,n)/zbeta(1,m)
         enddo
      enddo

c solve for (k_theta .ne. 0) modes:
      do l=2,kd
       do m=2,md
        do n=1,nd
           e_r1(l,m,n)=(e_r1(l,m,n)
     .          -aec(l,m)*e_r1(l-1,m,n))/zbeta(l,m)
        enddo
       enddo
      enddo

      do l=kd-1,1,-1
       do m=2,md
        do n=1,nd
           e_r1(l,m,n)=e_r1(l,m,n)
     .          -zgamma(l+1,m)*e_r1(l+1,m,n)
         enddo
       enddo
      enddo

c solve for (k_theta .eq. 0) modes:
      do l=2,kdpass
       m=1
        do n=1,nd
           e_r1(l,m,n)=(e_r1(l,m,n)
     .          -aec(l,m)*e_r1(l-1,m,n))/zbeta(l,m)
        enddo
      enddo

      do l=kdpass-1,1,-1
       m=1
        do n=1,nd
           e_r1(l,m,n)=e_r1(l,m,n)
     .          -zgamma(l+1,m)*e_r1(l+1,m,n)
         enddo
      enddo

      return
      end