subroutine init(old_qpar,old_qperp,old_e_density,time,dt)
c
c Set up arrays and define values of constants.
c
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 'netcdf.inc'
      include 'itg.par'
      include 'itg.cmn'
c arguments:
      complex old_qpar(lz,mz,nz,nspecz),old_qperp(lz,mz,nz,nspecz)
      complex old_e_density(kz,mz,nz)
      real time(mz),dt(mz)

c netcdf variables
      integer status,ncid
      character filename*80

c Local Variables:
      integer nw
      parameter(nw=500) ! # of steps in bounce-averages (must be even)

      integer i,l,m,n,itime(7),lp,lm,idum,ntablex,ifail,j,l0
      real z1,z2,z3,s21bcE,s21bbE,elliptic,kshift
      real theta,dtheta,dx,x,zkapmax,arg,arglim,dpsidrho
      real epsb,wt2(kz,lz),ellk(kz)
c Needed for debug (end of subroutine):
c      complex unity(lz,mz,nz)
c      real dummy(mz,nz),surfarea,dvdrhon
      real wrkinv(lz),invp(lz,lz),ident,temp(lz,lz)
      integer ntheta0,ntheta1,j1,l1,k
      real kapav,bratiop,taubp,kapp,kap2p
      real zeta1,zeta2,dzeta
      real rkymin,rkymax
      integer ntim,ntimr
      real dt1,tim, rmu1_sav
      real a1,b1,c1,d1,theta1,theta2,bmin,bmax
      integer lin_sav,ne_sav,nfreq_sav,nstp_sav,ntotal_sav,ifilter_sav
      character*30 x05abE
      external x05abE
      external x05aaE
      external s21bcE
      external s21bbE
c
      ldb=ld-1
c
c initialize arrays of m and n modes being used:
c
c nd		# of n modes
c md	# of m modes for each n
c
      malias=0
      nalias=0
      if(nffty .ne. 0 .and. nfftx .ne. 0) then
         malias=nffty
         nalias=nfftx
         if(md .eq. 0 .and. nd .eq. 0) then
            md=(nffty-1)/3+1
            nd=2*((nfftx-1)/3)+1
         endif
      endif

      if(nd.eq.-1) then
         nmin=0
         nmax=1
         nddm=1
         nddp=1
         nrr(1)=1
         nd=1
         goto 413
      endif
c
      if(mod(nd,2).eq.0) nd=nd+1
      if(nd.eq.1) then
         nmin=0
         nmax=0
         nddm=1
         nddp=1
         nrr(1)=0
      else
         nddm=nd/2
         nddp=nd/2+1
         nmin=-nd/2
         nmax=nd/2
         do n=1,nddp
            nrr(n)=n-1
         enddo
c     
         do n=1,nddm
            nrr(n+nddp)=n-nddm-1
         enddo
      endif
c
 413  continue
      if(ld .gt. lz) call aborter(6,
     &     'error: ld> lz (see itg.cmn and itg.par)')
      if(nd .gt. nz) call aborter(6,
     &     'error: nd> nz (see itg.cmn and itg.par)')
      if(md .gt. mz) call aborter(6,
     &     'error: md> mz (see itg.cmn and itg.par)')

      meq=-1	! set defaults in case there is no m=n=0 mode

      do n=1,nd
         if(md.gt.0) then
            mmax(n)=md-1
            mmin(n)=0
         else
            md=mhi-mlow+1
            mmin(n)=mlow
            mmax(n)=mhi
         endif
         do m=1,md
            mrr(m,n)=m+mmin(n)-1
c     
            if(mrr(m,n).eq.0.and.nrr(n).eq.0) then
               meq=m
               neq=n
            endif
c     
         enddo
      enddo

c Set up for real-to-complex y FFT's.
c md = # of independent ky modes (with ky>=0):
      mdspec=2*(md-1)           ! = # of ky points (without dealiasing)
      if(malias .eq. 0) then
         malias=4*(md-1)        ! overkill (more efficient to specify nffty)
      endif
      if(malias .gt. mffty) call aborter(6,
     &     'error: nffty or malias > mffty (see itg.in and itg.par)')
      
      madd=malias-mdspec
      rootm=sqrt(float(malias))
      
c     initialize trig tables for y FFT's:
      idum=0
      call csfftm(0,malias,1,1.0,idum,idum,idum,idum,tabley,idum,0)
c     
c     Set up for complex-to-complex x FFT's.
c     nd = # of independent kx modes (usually odd for symmetry around kx=0)
c     nalias = # of x points (with dealiasing):
      if(nd.eq.1) nalias=1
      if(nalias .eq. 0) then
         nalias=2*(nd-1)        ! overkill (more efficient to specify nfftx)
      endif
      if(nalias .gt. mfftx) call aborter(6,
     &     'error: nfftx or nalias > mfftx (see itg.in and itg.par)')
      
      nadd=nalias-nd
      rootn=sqrt(float(nalias))
      
c     initialize trig tables for x FFT's:
      ntablex=100+2*nzz
      call xmcfft(0,nalias,1,1.0,0.0,idum,idum,0.0,idum,idum,
     &     tablex,ntablex,1.0,idum)
c     
c     finish setting default m=0,n=0 mode if not already specified:
      if(meq .eq. -1) then
         neq=1
         meq=md+1
         if(meq .gt. mz) call aborter(6,'no room for m=0 mode')
      endif
c     
c     rkkin is the coefficient for the Hammett-Perkins
c     Landau damping model.
c     
c     rkkin = chi_1 * sqrt(2) * v_ti / c_s
c     = 2 * (2/pi*Ti/Te)**0.5,
c     
c     
      rkkin=2.*(2./pi)**0.5
      rkkperp=(2./pi)**0.5
      rkkperp2=(pi/(2.))**0.5
c     
c     zDq1 and zBq1 are the coefficients for the 4-pole Landau damping model
c     from Greg's PRL:
c     
c     zDq1 = D1 * sqrt(2) * Sqrt(TiovTe)
c     D1 = 2 sqrt(pi) / (3 pi -8)
c     zBq1 = (3 + 2*B1)*TiovTe
c     zEq1 = (3 + 2*B1)
c     B1 = (32 - 9 pi) / (6 pi - 16)
      
      zDq1 = 2.*sqrt(pi)/(3.0*pi-8.0)*sqrt(2.)
      zBq1 = (3.+2.*(32.-9.*pi)/(6.*pi-16.))
      zEq1 = 3.+2.*(32.-9.*pi)/(6.*pi-16.)
c     
      if(lin.eq.1) igradon=0
c     
      avgflag=0.0
      navgflag=0
      ncycle = 0
      ne=0
c     
      lhistory=float(ld)*.6+1
c     changed mbeer (old way didn't work for large x0)
c     
      nfreq=(nstp-1)/nez+1
      nfreqcnt=1
      if (nfreq.ne.0) then
         ntotal=nstp/nfreq
      else
         ntotal=1
      endif
c     
c     r(l) is the grid for the radial coordinate:
c     
      do l=1,ld
         r(l)=2.*x0*float(l-1)/ldb-x0
      enddo
c     
      do l=2,ldb
         dr(l)=r(l+1)-r(l-1)
      enddo
      dr(1)=2.*(r(2)-r(1))
      dr(ld)=2.*(r(ld)-r(ldb))

      write(9,*) 'Theta_0:'
      do m=1,md
         do n=1,nd
            if (mrr(m,n).eq.0.or.(lin.eq.1 .and. layout.eq.'linear')) then
               if (nmax.ne.0) then
                  r0(m,n)=z0*float(nrr(n))/nmax
               else
                  r0(m,n)=z0*float(nrr(n))
               endif
            else
               if (nmax.ne.0) then
                  r0(m,n)=z0*float(nrr(n))/nmax/float(mrr(m,n))
               else
                  r0(m,n)=z0*float(nrr(n))/float(mrr(m,n))
               endif
            endif
            write(9,*) m,n,' r0= ',r0(m,n)
         enddo
      enddo
c     
c Set up diffusivities:
c
      if(rmu1.eq.-1.) ifilter=0
      if(icrit.eq.2) then
         ifilter=1
         rmu1=max(rmu1,0.01)
      endif
      if(ifilter.gt.0) then
         do n=1,nd
            do m=1,md
               rdiff(m,n)=rmu1
            enddo
         enddo
      else
         rmu1=float(n_p)     !so it will appear in old output
      endif
c
c Prepare for theta FFT's
c
      if (iperiod.eq.3) then 
         call connect
      else
         ndom=1
         ncon(1)=md*nd
         scon(1)=1
         do m=1,md
            do n=1,nd
               mcon(m,n)=n+(m-1)*nd
               conpos(m,n)=1
               icon(m,n)=1
            enddo
         enddo   
      endif  
c     
c     define ky grid
c     
      do n=1,nd
         do m=1,md
            rky(m,n)=(mrr(m,n)/y0) ! ky
            rky2(m,n)=(mrr(m,n)/y0)**2 ! ky**2
         enddo
      enddo
c     
c     center of theta axis moved to zero, r0=theta0, mbeer
c     
c     b_mag(l) is the magnitude of the magnetic field B normalized to B_0
c     bmaginv(l) is 1/B
c     bgrad(l) is 1/B dB/dl    
c     gpar is the coefficient of d/dz (1/qR if igeo=0)
c
c     All in the circular concentric small aspect ratio limit unless igeo=1.
c     
      if(igeo.eq.1) then
         call geometry(dpsidrho)
         do l=1,ld
            bmaginv(l)=1./b_mag(l)
         enddo
c Use goto statement only because of compiler vectorization problems:
         goto 100
      endif

      dpsidrho=1.0
c
c Backwards compatibility:
c
      do l=1,ld
         gradrho(l)=1.0
      enddo

      gpar=epsn/qsf
      do l=1,ld
         b_mag(l)=1./(1.+eps*cos(r(l)))
         bmaginv(l)=1./b_mag(l)
         bgrad(l)=gpar*eps*sin(r(l))*b_mag(l)
      enddo   
c     
c     define kx grid:
c     
c     The current version (6.0.1.1) of Cray's fpp has a bug in the following
c     nested loops, unless vectorization is turned off:
cfpp$ skip R
      do n=1,nd
         do m=1,md
c     
c     Need switch for m=0 modes.
c     
            do l=1,ld
               rkx(l,m,n)=(mrr(m,n)/y0)*((r(l)-r0(m,n))*shr
     &              -alpha*(sin(r(l))-sin(r0(m,n))))
               if (m.eq.meq .and. md.ne.1) 
     .              rkx(l,m,n)=mrr(m+1,n)/y0*(-r0(m+1,n))*shr
               if (md.eq.1 .and. mrr(m,n).eq.0) rkx(l,m,n)=1./y0
c     
c     local if qsf < 0
c     
               if (qsf.lt.0) rkx(l,m,n)=0.0
c     
c     rkx(l,m,n)=rkx(l,m,n)*bmaginv(l)
c     This will be done properly sometime...
c
               rkx2(l,m,n)=rkx(l,m,n)**2
c               rkperp2(l,m,n)=rky2(m,n)+rkx2(l,m,n)
c
c     should be (k_perp rho)**2
c
               rkperp2(l,m,n)=(rky2(m,n)+rkx2(l,m,n))*bmaginv(l)**2
            enddo
         enddo
      enddo
c     
cfpp$ noskip R
c     
c     omegad(l) is the w_d profile (for toroidal terms) 
c     
cfpp$ skip R
      do i=1,nspecies
         do m=1,md
            do n=1,nd
               do l=1,ld
                  omegad(l,m,n,i)=
     .                 epsn*rky(m,n)*rho(i)*vt(i)*
     .                 ( cos(r(l))+(shr*(r(l)-r0(m,n))
     .                 -alpha*(sin(r(l))-sin(r0(m,n))))*sin(r(l)))
c     
c     special case for k_theta=0, k_r=const.     
c     
                  if ((m.eq.meq).or.((md.eq.1).and.mrr(m,n).eq.0)) then
                     omegad(l,m,n,i)=epsn*
     .                    rho(i)*vt(i)*rkx(l,m,n)*sin(r(l))
	             if(iflow.eq.0) omegad(l,m,n,i)=0.
                  endif
c     
c     local if qsf<0
c     
                  if (qsf.lt.0) omegad(l,m,n,i)=epsn*rky(m,n)
     .                 *rho(i)*vt(i)
               enddo
            enddo
         enddo
      enddo
cfpp$ noskip R

 100  continue
c
c Define the Jacobian of the coordinate transformation:
c (The abs() is probably overkill.)
c
      do l=1,ld
         jacobian(l)=abs(dpsidrho/(gpar*b_mag(l)))
c         write(*,*) r(l),'J_itgc= ',jacobian(l)
      enddo
c     
c Debug:
c
c      i=1
c      do n=1,nd
c         do m=1,md
c            do l=1,ld
c               write(*,2000) l,m,n,i,omegad(l,m,n,i)
c               write(*,2001) l,bgrad(l),
c     .              epsn*eps/qsf*sin(r(l))*b_mag(l),b_mag(l)
c               write(*,2002) l,m,n,i,rkperp2(l,m,n)
c               write(*,2003) l,gpar,gradrho(l),r(l)
c               write(*,*) l,gbd0(l)+cvd0(l),-shr*sin(r(l))
c               write(*,*) l,gbd(l)+cvd(l),cos(r(l))+shr*r(l)*sin(r(l))
c               write(*,*) 
c            enddo
c            write(*,*) 
c         enddo
c      enddo
 2000 format('wd(',i2,',',i2,',',i2,',',i2,')= ',f13.6)
 2001 format('grad ln(B)(',i2,')= ',f13.6,' ',f13.6,' ',f13.6)
 2002 format('kperp**2(',i2,',',i2,',',i2,',',i2,')= ',f13.6)
 2003 format('gradpar(',i2,')= ',f13.6,' |grad.rho|= ',f13.6,' theta= ',
     .     f13.6)

c
c     Find cfl multipliers:
c     
c     cflx=float(ldb)/x0
c
c This line gives trouble if nd=1, but we don't run nonlinearly with nd=1,
c so just replace it with something similar:
c
c      cflx=shr*z0*float(nd)/(float(nd/2)*2*pi*y0)
      if(nd.eq.1) then
         cflx=shr*z0/(2.*pi*y0)
      else
         cflx=shr*z0*float(nd)/(float(nd/2)*2*pi*y0)
      endif
      cfly=float(mdspec)/(2.*pi*y0)
c     
c     kap(l) = kappa pitch-angle grid for bounce-averaged electrons
c     

      if(epse.gt.0.0) then
c This criterion for kd needs to be upgraded for the case x0>pi, as does
c much of the trapped electron calculation:
      kd=ldb/2   ! covers just the trapped particles
      kdpass=ldb ! covers all particles, including passing particles

      if (igeo.eq.0) then
c
c     circular concentric limit
c
      do l=1,kd
         j=ldb/2+l
         kap(l)=sin(r(j)/2.)
         kap(l)=max(0.0,kap(l))
      enddo
c     kap(1)=kap(2)/4  ! put the first grid point slightly off kap=0
         
      do l=1,kd
         ellk(l)=s21bbE(0.,1.-kap(l)**2,1.,ifail)
      enddo
         
      epsb=epse/(1.+epse)

      else
c
c     general geometry
c
      bmin=b_mag(1)
      bmax=b_mag(1)
      do l=1,ld
         if (b_mag(l).gt.bmax) bmax=b_mag(l)
         if (b_mag(l).lt.bmin) bmin=b_mag(l)
      enddo
      write(6,*) 'bmin,bmax: ',bmin,bmax

      epsb=(bmax-bmin)/2/bmax
      do l=1,kd
         j=ldb/2+l
         kap(l)=sqrt((1-bmin/b_mag(j))/2/epsb)
         kap(l)=max(0.0,kap(l))
      enddo

      endif

      if(epse .gt. 0.0) then
         zkapmax=1./sqrt(2.*epsb)
      else
         zkapmax=2.0
      endif

c eventually would like the kap grid to be more closely space
c around the kap=1 point for kap>1 as well as for kap<1.  For now
c just use a uniform grid, but at least put the first point at
c kap>1 equally spaced from the boundary as the kap<1 point, so
c that kap=1 is exactly between two grid points:

      do l=kd+1,kdpass
         kap(l)=1.0+(1.0-kap(kd))+(l-kd-1)*(zkapmax-1.0)/(kdpass-kd)
      enddo
      kshift=(2-kap(kd)-zkapmax*(1-kap(kd)))/kap(kd)
      do l=kd+1,kdpass
         kap(l)=(1-kap(kd-l+1+kd))*(zkapmax-kshift)+kshift
         ellk(l)=s21bbE(0.,1.-1./kap(l)**2,1.,ifail)
      enddo
         
c     define a mirror point on the other side of the maximum kappa
c     to make differencing formulas easier:
      kap(kdpass+1)=zkapmax+(zkapmax-kap(kdpass))
         
c     Define dkap as the width of the zone about each grid point:
c     dkap(1)=kap(2)/2
c     do l=2,kdpass
c        dkap(l)=(kap(l+1)-kap(l-1))/2
c     enddo
      do l=1,kd-1
         dkap(l)=kap(l+1)-kap(l)
      enddo
      dkap(kd)=1.-kap(kd)
      dkap(kd+1)=kap(kd+1)-1.
      do l=kd+2,kdpass
         dkap(l)=kap(l)-kap(l-1)
      enddo   
c
c bounce averaged precession frequency:
c
cfpp$ skip R
      do m=1,md
         do n=1,nd
            do l=1,kd
               if (kap(l).ne.1.0) then
                  elliptic=1.-kap(l)**2/3.*
     .                 s21bcE(0.,1.-kap(l)**2,1.,ifail)
     .                 /s21bbE(0.,1.-kap(l)**2,1.,ifail)
               else
                  elliptic=0.
               endif  
               omegade(l,m,n)=2.*epsn*rky(m,n)*(elliptic-0.5
     .                 +2.*shr*(elliptic+kap(l)**2-1.)
     .      -alpha*2./3.*(1-kap(l)**2-(1-2*kap(l)**2)*elliptic))
c     
c     special case for k_theta=0, k_r=const.     
c     
               if ((m.eq.meq).or.((md.eq.1).and.mrr(m,n).eq.0))
     .              omegade(l,m,n)=0.
c     
c     local if qsf<0
c
c If igeo=1, qsf can be negative, so in the new section, 
c don't keep this option:
c
               if (qsf.lt.0 .and. igeo.eq.0) 
     .              omegade(l,m,n)=2.*epsn*rky(m,n)
            enddo
         enddo
      enddo
cfpp$ noskip R
      endif ! end of if(epse.eq.0.0)

c     
c     Twisted periodicity boundary conditions parameters:
c     
      if (iperiod.eq.2) then
         l_left=float(ld-1)/2.*(1.-xp/x0)+1.+.5
         l_right=float(ld-1)/2.*(1.+xp/x0)+1.+.5
         write(9,*) 'l_left, l_right: ', l_left, l_right
      else
         l_left=1
         l_right=ld
      endif

      do i=1,nspecies
         do n=1,nd
            do m=1,md
               do l=1,ld
                  pfilter2(l,m,n)=0.
                  potential(l,m,n)=0.
                  density(l,m,n,i)=0.
                  u_par(l,m,n,i)=0.
                  t_par(l,m,n,i)=0.
                  q_par(l,m,n,i)=0.
                  t_perp(l,m,n,i)=0.
                  q_perp(l,m,n,i)=0.
c
c The next two lines prevent an uninitialized variable problem
c for a 3+1 model:
c     
                  old_qpar(l,m,n,i)=0.
                  old_qperp(l,m,n,i)=0.
                  nl_density1(l,m,n,i)=0.
                  nl_density2(l,m,n,i)=0.
                  nl_upar1(l,m,n,i)=0.
                  nl_upar2(l,m,n,i)=0.
                  nl_tpar(l,m,n,i)=0.
                  nl_qpar(l,m,n,i)=0.
                  nl_tperp1(l,m,n,i)=0.
                  nl_tperp2(l,m,n,i)=0.
                  nl_tperp3(l,m,n,i)=0.
                  nl_qperp1(l,m,n,i)=0.
                  nl_qperp2(l,m,n,i)=0.
               enddo
            enddo
         enddo
      enddo
      do n=1,nd
         do m=1,md
            do l=1,kdpass
               e_density(l,m,n)=0.
               e_p(l,m,n)=0.
               e_r(l,m,n)=0.
               e_t(l,m,n)=0.
c  fix initialization problem:
               old_e_density(l,m,n)=0.
               nl_edensity(l,m,n)=0.
               nl_ep(l,m,n)=0.
               nl_er(l,m,n)=0.
               nl_et(l,m,n)=0.
               phi_ba(l,m,n)=0.0
            enddo
            do l=1,ld
               phi_bk(l,m,n)=0.
               e_denk(l,m,n)=0.
            enddo   
         enddo
      enddo   
      if(nread.eq.1) then 
c
c Read previous *.ncp results file:
c
         ntim=0
         ntimr=0
         ne_sav=ne           ! save some values
         ntotal_sav=ntotal
         nstp_sav=nstp
         nfreq_sav=nfreq
         lin_sav=lin
         ifilter_sav=ifilter
         rmu1_sav=rmu1

         if(binary) then
            open(unit=7,file=runname(1:lrunname)//'.bresp',
     .           form='unformatted',status='old')
            call bwread(7,ntim,ntimr,tim,dt1,time,dt)
            close(unit=7)
         else
            filename=runname(1:lrunname)//'.ncp'
            status = nf_open(filename, 0, ncid) ! read only
            IF (status .ne. nf_noerr) 
     .          print *, filename, ' not found'
            call wreadnc(ncid,ntim,ntimr,tim,dt1,time,dt)
            status = nf_close(ncid)
         endif
c
c     restore B variation of kperp rho which was removed for post
c
      do l=1,ld
         do m=1,md
            do n=1,nd
               rkperp2(l,m,n)=rkperp2(l,m,n)*bmaginv(l)**2
            enddo
         enddo
      enddo

c Reset some values just read from the previous results file,
c for example, restart time indexing from the beginning:
         ne=ne_sav
         ntotal=ntotal_sav
         nstp=nstp_sav
         nfreq_sav=nfreq
         lin=lin_sav
         ifilter=ifilter_sav
         rmu1=rmu1_sav

      else
         call pertb
c
c If icrit=1,2, advance each k_theta with a different dt.
c     
         if(icrit.ne.0) then
            rkymin=rky(1,1)
            if(rkymin.eq.0.) then
               rkymin=1.
               if(md.gt.1) rkymin=rky(2,1)
            endif
            
            rkymax=rky(md,1)
            if(rkymax.eq.0.) rkymax=1.
            
            n=1
            do m=1,md
               time(m)=0.0
               dt(m)=rkymax/max(rkymin,rky(m,n))*dt0
            enddo
         else
            do m=1,md
               dt(m)=dt0
            enddo
         endif
      endif
c     
c Allow for k_perp dependent etas for linear estimates:
c
      do i=1,nspecies
         do n=1,nd
            do m=1,md
               etak(m,n,i)=eta(i)
               etak_par(m,n,i)=eta_par(i)
            enddo
         enddo
      enddo
      call flrinit(1)
c
c Get the time and date of this run (used to label plots):
c
      call x05aaE(itime)
      if(debug_plotlabel .eq. ' ' ) then
         date=x05abE(itime)
      else
         date=debug_plotlabel
      endif
c
c weights for bounce averages
c
      if (epse.gt.0.0) then
c The algorithm for the integration below could probably be more
c effficient and converge faster (so that nw doesn't need to be
c so large), perhaps using unequally spaced points (like sccf.f 
c in my FPP code) or using a canned routine.  
c But it works for now...
c
      if (igeo.eq.0) then
c
c     circular concentric limit
c
      do l=1,kdpass
         do j=1,ld
            wgtba1(l,j)=0.
         enddo
      enddo

      dtheta=2.*x0/float(ldb)
c
c     trapped
c
      do l=2,kd
         ntheta0=ldb/2+l
         ntheta1=ldb/2+2-l
         j=ntheta1
c use variable transformation: sin(theta/2)=kappa sin(x)
         if (sin(r(j)/2)/kap(l).lt.-1.0) then 
            zeta1=-pi/2
         else
            zeta1=asin(sin(r(j)/2)/kap(l))
         endif   
         zeta2=asin(sin((r(j)+dtheta/2)/2)/kap(l))
c         write(*,101) kap(l),r(j),r(j)+dtheta/2,
c     .        sin(r(j)/2)/kap(l),sin((r(j)+dtheta/2)/2)/kap(l),zeta1,zeta2
 101     format(7f15.9)
         dzeta=zeta2-zeta1
         do i=1,nw
            x=(float(i)-0.5)/float(nw)*dzeta+zeta1
            dx=1./float(nw)*dzeta
            arg=1-kap(l)**2*sin(x)**2
            wgtba1(l,j)=wgtba1(l,j)+
     .           2.*dx*sqrt((1-epse+2*epse*arg)/arg/(1+epse))
         enddo   
         do j=ntheta1+1,ntheta0-1
            zeta1=asin(sin((r(j)-dtheta/2)/2)/kap(l))
            zeta2=asin(sin((r(j)+dtheta/2)/2)/kap(l))
c            write(*,101) kap(l),r(j)-dtheta/2,r(j)+dtheta/2,
c     . sin((r(j)-dtheta/2)/2)/kap(l),sin((r(j)+dtheta/2)/2)/kap(l),zeta1,zeta2
               dzeta=zeta2-zeta1
            do i=1,nw
               x=(float(i-1)+0.5)/float(nw)*dzeta+zeta1
               dx=1./float(nw)*dzeta
               arg=1-kap(l)**2*sin(x)**2
               wgtba1(l,j)=wgtba1(l,j)+
     .              2.*dx*sqrt((1-epse+2*epse*arg)/arg/(1+epse))
            enddo   
         enddo 
         j=ntheta0
         zeta1=asin(sin((r(j)-dtheta/2)/2)/kap(l))
         if (sin(r(j)/2)/kap(l).gt.1.0) then
            zeta2=pi/2
         else   
            zeta2=asin(sin(r(j)/2)/kap(l))
         endif   
         dzeta=zeta2-zeta1
c         write(*,101) kap(l),r(j)-dtheta/2,r(j),
c     .        sin((r(j)-dtheta/2)/2)/kap(l),sin(r(j)/2)/kap(l),zeta1,zeta2
         do i=1,nw
            x=(float(i)-0.5)/float(nw)*dzeta+zeta1
            dx=1./float(nw)*dzeta
            arg=1-kap(l)**2*sin(x)**2
            wgtba1(l,j)=wgtba1(l,j)+
     .           2.*dx*sqrt((1-epse+2*epse*arg)/arg/(1+epse))
         enddo   
      enddo
c
c     passing
c
      do l=kd+1,kdpass
         do j=1,ldb
            wgtba1(l,j)=0.
            do i=1,nw
               x=(float(i-1)/float(nw)+0.5/float(nw))*dtheta-dtheta/2
               dx=1./float(nw)
               theta=r(j)+x
               arg=2*epse/(1+epse*cos(theta))*(kap(l)**2-sin(theta/2)**2)
               wgtba1(l,j)=wgtba1(l,j)+dx/sqrt(arg/2./epsb)
            enddo
            wgtba1(l,j)=wgtba1(l,j)*dtheta
         enddo
      enddo   
c
c ntheta0,ntheta1 are chosen so that:
c For trapped particles, the bounce-average symmetrically involves 
c phi points above and below the midplane.
c For passing particles, the bounce-average involves all phi points
c from l=1 to ldb (the ld point is identical to the l=1 point).
c

c kap=0 limiting value:
      wgtba1(1,ldb/2+1)=2*pi  ! limiting value...

c calculate bounce-time numerically (for passing particles also...):
      do l=1,kd
         ntheta0=ldb/2+l
         ntheta1=ldb/2+2-l
         do j=ntheta1,ntheta0
            taub(l)=taub(l)+wgtba1(l,j)
         enddo   
      enddo
      do l=kd+1,kdpass
         taub(l)=0.
         do j=1,ldb
            taub(l)=taub(l)+wgtba1(l,j)
         enddo   
      enddo

c debug:
      do l=1,kdpass
         write(6,102) l,kap(l),taub(l),4*ellk(l),omegade(l,2,1)/rky(2,1)/epsn
 102     format( 'l,kap(l),taub(l),4*ellk(l) =',i3,4f14.7) 
      enddo

      do l=1,kdpass
         ellk(l)=taub(l)/4
      enddo

c
c ?? I think I've now generalized all of this subroutine to use
c a grid-centered approach, with bounce-averages symmetrically 
c involving phi above and below the midplane...
c
      do j=1,ld
         do l=1,kdpass
            wgtba2(l,j)=0.
         enddo
      enddo

      do j=ldb/2+1,ld
         l0=j-ldb/2
         do l=l0,kdpass
            if (l.le.kd) then
               wgtba2(l,j)=sqrt(2*epse/(1+epse*cos(r(j))))*
     .              (sqrt( (kap(l)+dkap(l))**2 - sin(r(j)/2)**2 )
     .              -sqrt( kap(l)**2 - sin(r(j)/2)**2 ))
            else
               wgtba2(l,j)=sqrt(2*epse/(1+epse*cos(r(j))))*
     .              (sqrt( kap(l)**2 - sin(r(j)/2)**2 )
     .              -sqrt( (kap(l)-dkap(l))**2 - sin(r(j)/2)**2 ))
            endif   
            wgtba2(l,j)=wgtba2(l,j)/dkap(l)
c B(theta)=1/(1+eps*cos(theta)) included to get particle conservation.
c symmetry, for theta<0:
c ?? GWH: also need to put a 1/B into the flux-surface averages...
            wgtba2(l,ldb-j+2)=wgtba2(l,j)
         enddo
      enddo   

c ?? GWH: The bounce-averages (wgtba1), kappa-averages (wgtba2), 
c and the form of the coefficients of the pitch-angle-scattering 
c operator, are related in various ways.  Some of these relations 
c are important for particle conservation, and so it may be possible 
c to derive some of these coefficients from each other in a way
c that preserves conservation properties exactly (or at least to
c within roundoff error).  For now, just stick to the above methods,
c as they should preserve the conservation properties to second
c order in the grid spacings...

c load matrix for poisson inversion
      do j=1,ldb
         l1=1+abs(j-ldb/2-1)
         do l=l1,kdpass
            if ((l.eq.l1).and.(l1.lt.kd)) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)/2
            elseif ((l.gt.l1).and.(l.lt.kd)) then
               wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j))/2
            elseif ((l.eq.kd).and.(l1.ne.kd)) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j)/2
            elseif ((l.eq.kd).and.(l1.eq.kd)) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)
            elseif (l.eq.kd+1) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j)/2
            elseif ((l.gt.kd+1).and.(l.lt.kdpass)) then
               wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j))/2
            elseif (l.eq.kdpass) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)/2
            endif   
         enddo   
         do j1=1,ldb
            pmat(j,j1)=0.
            l0=1+max(abs(j-ldb/2-1),abs(j1-ldb/2-1))
            do l=l0,kdpass
               pmat(j,j1)=pmat(j,j1)+wt2(l,j)*wgtba1(l,j1)/taub(l)
            enddo
            temp(j,j1)=pmat(j,j1)
         enddo
      enddo
c
c     invert
c
      ifail=0
c     call f01aaE(temp,lz,ldb/2+1,pmat,lz,wrkinv,ifail)
c     not used

      else
c
c     general geometry bounce average weights
c   
      do l=1,kdpass
         do j=1,ld
            wgtba1(l,j)=0.
         enddo
      enddo

      dtheta=2.*x0/float(ldb)
c
c     trapped  (assume B and J vary linearly over theta grid pt. interval)
c

c     special case for l=2
c     assume up-down symmetric for now...
c
c     c1<0
c
      l=2
      ntheta0=ldb/2+l
      ntheta1=ldb/2+2-l
      
c gpar is not a function of theta any more...
c      b1=1/gpar(ldb/2+1)
c      a1=(1/gpar(ldb/2)-b1)/dtheta**2
      b1=1/gpar
      a1=0.
      d1=b_mag(ldb/2+1)
      c1=(b_mag(ldb/2)-d1)/dtheta**2
      c1=-c1/bmin*(1-2*epsb*kap(l)**2)
      d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
      j=ntheta1
c      wgtba1(l,j)=a1/c1/2*(r(j)+dtheta/2)*sqrt(c1*(r(j)+dtheta/2)**2+d1)
c     . +(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))
c     .     *asin(sqrt(-c1/d1)*(r(j)+dtheta/2))
c     . -(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))*(-pi/2)
      wgtba1(l,j)=b1/sqrt(-c1)*(asin(sqrt(-c1/d1)*(r(j)+dtheta/2))+pi/2)
      j=ntheta1+1
c      wgtba1(l,j)=a1/c1/2*(r(j)+dtheta/2)*sqrt(c1*(r(j)+dtheta/2)**2+d1)
c     . +(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))
c     .     *asin(sqrt(-c1/d1)*(r(j)+dtheta/2))
c     . -a1/c1/2*(r(j)-dtheta/2)*sqrt(c1*(r(j)-dtheta/2)**2+d1)
c     . -(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))
c     .     *asin(sqrt(-c1/d1)*(r(j)-dtheta/2))
      wgtba1(l,j)=b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)+dtheta/2))
     .           -b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)-dtheta/2))
      j=ntheta0
c      wgtba1(l,j)=
c     . +(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))*(pi/2)
c     . -a1/c1/2*(r(j)-dtheta/2)*sqrt(c1*(r(j)-dtheta/2)**2+d1)
c     . -(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))
c     .     *asin(sqrt(-c1/d1)*(r(j)-dtheta/2))
      wgtba1(l,j)=b1/sqrt(-c1)*(pi/2-asin(sqrt(-c1/d1)*(r(j)-dtheta/2)))

      do l=3,kd
         ntheta0=ldb/2+l
         ntheta1=ldb/2+2-l
         j=ntheta1
c gpar not a function of theta
c         a1=(1/gpar(j+1)-1/gpar(j))/dtheta
c         b1=1/gpar(j)-r(j)/dtheta*(1/gpar(j+1)-1/gpar(j))
         a1=0.
         b1=1/gpar
         c1=(b_mag(j+1)-b_mag(j))/dtheta
         d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j))
         c1=-c1/bmin*(1-2*epsb*kap(l)**2)
         d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c         wgtba1(l,j)=2/c1**2*(a1*c1*(r(j)+dtheta/2)/3-2*a1*d1/3+b1*c1)
c     .        *sqrt(c1*(r(j)+dtheta/2)+d1)
cc     .        -2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
cc     .        *sqrt(c1*r(j)+d1)
         wgtba1(l,j)=2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1)
c     .        -2/c1*b1
c     .        *sqrt(c1*r(j)+d1)
         do j=ntheta1+1,ldb/2-1
c gpar not a function of theta
c            a1=(1/gpar(j)-1/gpar(j-1))/dtheta
c            b1=1/gpar(j-1)-r(j-1)/dtheta*(1/gpar(j)-1/gpar(j-1))
            a1=0.
            b1=1/gpar
            c1=(b_mag(j)-b_mag(j-1))/dtheta
            d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1))
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*r(j)+d1)
c     .           -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*(r(j)-dtheta/2)+d1)
            wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j)          +d1)
     .                 -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1)
c gpar not a function of theta
c            a1=(1/gpar(j+1)-1/gpar(j))/dtheta
c            b1=1/gpar(j)-r(j)/dtheta*(1/gpar(j+1)-1/gpar(j))
            a1=0.
            b1=1/gpar
            c1=(b_mag(j+1)-b_mag(j))/dtheta
            d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j))
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=wgtba1(l,j)+2/c1**2*(a1*c1*(r(j)+dtheta/2)/3
c     .           -2*a1*d1/3+b1*c1)*sqrt(c1*(r(j)+dtheta/2)+d1)
c     .           -2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*r(j)+d1)
            wgtba1(l,j)=wgtba1(l,j)+2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1)
     .                             -2/c1*b1*sqrt(c1* r(j)          +d1)
         enddo 
c     
c     near Bmin, use quadratic interpolation
c
         j=ldb/2
c gpar not a function of theta
c            a1=(1/gpar(j)-1/gpar(j-1))/dtheta
c            b1=1/gpar(j-1)-r(j-1)/dtheta*(1/gpar(j)-1/gpar(j-1))
            a1=0.
            b1=1/gpar
            c1=(b_mag(j)-b_mag(j-1))/dtheta
            d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1))
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*r(j)+d1)
c     .           -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*(r(j)-dtheta/2)+d1)
            wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j)          +d1)
     .                 -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1)
c gpar not a function of theta
c            b1=1/gpar(ldb/2+1)
c            a1=(1/gpar(ldb/2)-b1)/dtheta**2
            b1=1/gpar
            a1=0.
            d1=b_mag(ldb/2+1)
            c1=(b_mag(ldb/2)-d1)/dtheta**2
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=wgtba1(l,j)+
c     .         a1/c1/2*(r(j)+dtheta/2)*sqrt(c1*(r(j)+dtheta/2)**2+d1)
c     .         +(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))
c     .           *asin(sqrt(-c1/d1)*(r(j)+dtheta/2))
c     .         -a1/c1/2*r(j)*sqrt(c1*r(j)**2+d1)
c     .         -(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))
c     .           *asin(sqrt(-c1/d1)*r(j))
            wgtba1(l,j)=wgtba1(l,j)
     .         +b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)+dtheta/2))
     .         -b1/sqrt(-c1)*asin(sqrt(-c1/d1)* r(j))
         j=ldb/2+1
c gpar not a function of theta
c            b1=1/gpar(ldb/2+1)
c            a1=(1/gpar(ldb/2)-b1)/dtheta**2
            b1=1/gpar
            a1=0.
            d1=b_mag(ldb/2+1)
            c1=(b_mag(ldb/2)-d1)/dtheta**2
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
            wgtba1(l,j)=
     .          b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)+dtheta/2))
     .         -b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)-dtheta/2))
         j=ldb/2+2
c gpar not a function of theta
c            b1=1/gpar(ldb/2+1)
c            a1=(1/gpar(ldb/2)-b1)/dtheta**2
            b1=1/gpar
            a1=0.
            d1=b_mag(ldb/2+1)
            c1=(b_mag(ldb/2)-d1)/dtheta**2
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=a1/c1/2*r(j)*sqrt(c1*r(j)**2+d1)
c     .         +(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))
c     .           *asin(sqrt(-c1/d1)*r(j))
c     .         -a1/c1/2*(r(j)-dtheta/2)*sqrt(c1*(r(j)-dtheta/2)**2+d1)
c     .         -(b1/sqrt(-c1)-a1*d1/2/c1/sqrt(-c1))
c     .           *asin(sqrt(-c1/d1)*(r(j)-dtheta/2))
            wgtba1(l,j)=
     .          b1/sqrt(-c1)*asin(sqrt(-c1/d1)* r(j)          )
     .         -b1/sqrt(-c1)*asin(sqrt(-c1/d1)*(r(j)-dtheta/2))
c gpar not a function of theta
c            a1=(1/gpar(j+1)-1/gpar(j))/dtheta
c            b1=1/gpar(j)-r(j)/dtheta*(1/gpar(j+1)-1/gpar(j))
            a1=0.
            b1=1/gpar
            c1=(b_mag(j+1)-b_mag(j))/dtheta
            d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j))
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=wgtba1(l,j)+2/c1**2*(a1*c1*(r(j)+dtheta/2)/3
c     .           -2*a1*d1/3+b1*c1)*sqrt(c1*(r(j)+dtheta/2)+d1)
c     .           -2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*r(j)+d1)
            wgtba1(l,j)=wgtba1(l,j)
     .           +2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1)
     .           -2/c1*b1*sqrt(c1* r(j)          +d1)
         do j=ldb/2+3,ntheta0-1
c            a1=(1/gpar(j)-1/gpar(j-1))/dtheta
c            b1=1/gpar(j-1)-r(j-1)/dtheta*(1/gpar(j)-1/gpar(j-1))
            a1=0.
            b1=1/gpar
            c1=(b_mag(j)-b_mag(j-1))/dtheta
            d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1))
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*r(j)+d1)
c     .           -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*(r(j)-dtheta/2)+d1)
            wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j)          +d1)
     .                 -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1)
c            a1=(1/gpar(j+1)-1/gpar(j))/dtheta
c            b1=1/gpar(j)-r(j)/dtheta*(1/gpar(j+1)-1/gpar(j))
            a1=0.
            b1=1/gpar
            c1=(b_mag(j+1)-b_mag(j))/dtheta
            d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j))
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=wgtba1(l,j)+2/c1**2*(a1*c1*(r(j)+dtheta/2)/3
c     .           -2*a1*d1/3+b1*c1)*sqrt(c1*(r(j)+dtheta/2)+d1)
c     .           -2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*r(j)+d1)
            wgtba1(l,j)=wgtba1(l,j)+2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1)
     .                             -2/c1*b1*sqrt(c1* r(j)          +d1)
         enddo 
         j=ntheta0
c         a1=(1/gpar(j)-1/gpar(j-1))/dtheta
c         b1=1/gpar(j-1)-r(j-1)/dtheta*(1/gpar(j)-1/gpar(j-1))
         a1=0.
         b1=1/gpar
         c1=(b_mag(j)-b_mag(j-1))/dtheta
         d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1))
         c1=-c1/bmin*(1-2*epsb*kap(l)**2)
         d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
cc         wgtba1(l,j)=2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
cc     .           *sqrt(c1*r(j)+d1)
cc     .           -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
cc     .           *sqrt(c1*(r(j)-dtheta/2)+d1)
c         wgtba1(l,j)=
c     .           -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*(r(j)-dtheta/2)+d1)

c         wgtba1(l,j)=2/c1*b1
c     .           *sqrt(c1*r(j)+d1)
c     .           -2/c1*b1
c     .           *sqrt(c1*(r(j)-dtheta/2)+d1)
         wgtba1(l,j)=-2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1)

      enddo
c
c     passing
c
      do l=kd+1,kdpass-1
         j=1
c gpar not a function of theta
c         a1=(1/gpar(j)-1/gpar(ldb))/dtheta
c         b1=1/gpar(ldb)-r(ldb)/dtheta*(1/gpar(j)-1/gpar(ldb))
         a1=0.
         b1=1/gpar
         c1=(b_mag(j)-b_mag(ldb))/dtheta
         d1=b_mag(ldb)-r(ldb)/dtheta*(b_mag(j)-b_mag(ldb))
         c1=-c1/bmin*(1-2*epsb*kap(l)**2)
         d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c         wgtba1(l,j)=2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .        *sqrt(c1*r(j)+d1)
c     .        -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
c     .        *sqrt(c1*(r(j)-dtheta/2)+d1)
         wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j)          +d1)
     .              -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1)
c         a1=(1/gpar(j+1)-1/gpar(j))/dtheta
c         b1=1/gpar(j)-r(j)/dtheta*(1/gpar(j+1)-1/gpar(j))
         a1=0.
         b1=1/gpar
         c1=(b_mag(j+1)-b_mag(j))/dtheta
         d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j))
         c1=-c1/bmin*(1-2*epsb*kap(l)**2)
         d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c         wgtba1(l,j)=wgtba1(l,j)+2/c1**2*(a1*c1*(r(j)+dtheta/2)/3
c     .        -2*a1*d1/3+b1*c1)*sqrt(c1*(r(j)+dtheta/2)+d1)
c     .        -2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .        *sqrt(c1*r(j)+d1)
         wgtba1(l,j)=wgtba1(l,j)
     .        +2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1)
     .        -2/c1*b1*sqrt(c1* r(j)          +d1)
         do j=2,ldb-1
c            a1=(1/gpar(j)-1/gpar(j-1))/dtheta
c            b1=1/gpar(j-1)-r(j-1)/dtheta*(1/gpar(j)-1/gpar(j-1))
            a1=0.
            b1=1/gpar
            c1=(b_mag(j)-b_mag(j-1))/dtheta
            d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1))
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*r(j)+d1)
c     .           -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*(r(j)-dtheta/2)+d1)
            wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j)          +d1)
     .                 -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1)
c            a1=(1/gpar(j+1)-1/gpar(j))/dtheta
c            b1=1/gpar(j)-r(j)/dtheta*(1/gpar(j+1)-1/gpar(j))
            a1=0.
            b1=1/gpar
            c1=(b_mag(j+1)-b_mag(j))/dtheta
            d1=b_mag(j)-r(j)/dtheta*(b_mag(j+1)-b_mag(j))
            c1=-c1/bmin*(1-2*epsb*kap(l)**2)
            d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c            wgtba1(l,j)=wgtba1(l,j)+2/c1**2*(a1*c1*(r(j)+dtheta/2)/3
c     .           -2*a1*d1/3+b1*c1)*sqrt(c1*(r(j)+dtheta/2)+d1)
c     .           -2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .           *sqrt(c1*r(j)+d1)
            wgtba1(l,j)=wgtba1(l,j)
     .           +2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1)
     .           -2/c1*b1*sqrt(c1* r(j)          +d1)
         enddo
         j=ldb   
c         a1=(1/gpar(j)-1/gpar(j-1))/dtheta
c         b1=1/gpar(j-1)-r(j-1)/dtheta*(1/gpar(j)-1/gpar(j-1))
         a1=0.
         b1=1/gpar
         c1=(b_mag(j)-b_mag(j-1))/dtheta
         d1=b_mag(j-1)-r(j-1)/dtheta*(b_mag(j)-b_mag(j-1))
         c1=-c1/bmin*(1-2*epsb*kap(l)**2)
         d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c         wgtba1(l,j)=2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .        *sqrt(c1*r(j)+d1)
c     .        -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
c     .        *sqrt(c1*(r(j)-dtheta/2)+d1)
         wgtba1(l,j)=2/c1*b1*sqrt(c1* r(j)          +d1)
     .              -2/c1*b1*sqrt(c1*(r(j)-dtheta/2)+d1)
c         a1=(1/gpar(1)-1/gpar(j))/dtheta
c         b1=1/gpar(j)-r(j)/dtheta*(1/gpar(1)-1/gpar(j))
         a1=0.
         b1=1/gpar
         c1=(b_mag(1)-b_mag(j))/dtheta
         d1=b_mag(j)-r(j)/dtheta*(b_mag(1)-b_mag(j))
         c1=-c1/bmin*(1-2*epsb*kap(l)**2)
         d1=1-d1/bmin*(1-2*epsb*kap(l)**2)
c         wgtba1(l,j)=wgtba1(l,j)+2/c1**2*(a1*c1*(r(j)+dtheta/2)/3
c     .        -2*a1*d1/3+b1*c1)*sqrt(c1*(r(j)+dtheta/2)+d1)
c     .        -2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     .        *sqrt(c1*r(j)+d1)
         wgtba1(l,j)=wgtba1(l,j)
     .        +2/c1*b1*sqrt(c1*(r(j)+dtheta/2)+d1)
     .        -2/c1*b1*sqrt(c1* r(j)          +d1)
      enddo   
c
c     need special case for l=kdpass, where denominator is unity
c
      l=kdpass
         j=1
c         a1=(1/gpar(j)-1/gpar(ldb))/dtheta
c         b1=1/gpar(ldb)-r(ldb)/dtheta*(1/gpar(j)-1/gpar(ldb))
         a1=0.
         b1=1/gpar
c         wgtba1(l,j)=a1*r(j)**2/2+b1*r(j)
c     .        -a1*(r(j)-dtheta/2)**2/2-b1*(r(j)-dtheta/2)
         wgtba1(l,j)=b1*dtheta/2
c         a1=(1/gpar(j+1)-1/gpar(j))/dtheta
c         b1=1/gpar(j)-r(j)/dtheta*(1/gpar(j+1)-1/gpar(j))
         a1=0.
         b1=1/gpar
c         wgtba1(l,j)=wgtba1(l,j)+a1*(r(j)+dtheta/2)**2/2
c     .        +b1*(r(j)+dtheta/2)
c     .        -a1*r(j)**2/2-b1*r(j)
         wgtba1(l,j)=wgtba1(l,j)+b1*dtheta/2

         do j=2,ldb-1
c            a1=(1/gpar(j)-1/gpar(j-1))/dtheta
c            b1=1/gpar(j-1)-r(j-1)/dtheta*(1/gpar(j)-1/gpar(j-1))
            a1=0.
            b1=1/gpar
c            wgtba1(l,j)=a1*r(j)**2/2+b1*r(j)
c     .           -a1*(r(j)-dtheta/2)**2/2-b1*(r(j)-dtheta/2)
            wgtba1(l,j)=b1*dtheta/2

c            a1=(1/gpar(j+1)-1/gpar(j))/dtheta
c            b1=1/gpar(j)-r(j)/dtheta*(1/gpar(j+1)-1/gpar(j))
            a1=0.
            b1=1/gpar
c            wgtba1(l,j)=wgtba1(l,j)+a1*(r(j)+dtheta/2)**2/2
c     .           +b1*(r(j)+dtheta/2)
c     .           -a1*r(j)**2/2-b1*r(j)
            wgtba1(l,j)=wgtba1(l,j)+b1*dtheta/2
         enddo
         j=ldb   
c         a1=(1/gpar(j)-1/gpar(j-1))/dtheta
c         b1=1/gpar(j-1)-r(j-1)/dtheta*(1/gpar(j)-1/gpar(j-1))
         a1=0.
         b1=1/gpar
c         wgtba1(l,j)=a1*r(j)**2/2+b1*r(j)
c     .        -a1*(r(j)-dtheta/2)**2/2-b1*(r(j)-dtheta/2)
         wgtba1(l,j)=b1*dtheta/2

c         a1=(1/gpar(1)-1/gpar(j))/dtheta
c         b1=1/gpar(j)-r(j)/dtheta*(1/gpar(1)-1/gpar(j))
         a1=0.
         b1=1/gpar
c         wgtba1(l,j)=wgtba1(l,j)+a1*(r(j)+dtheta/2)**2/2
c     .        +b1*(r(j)+dtheta/2)
c     .        -a1*r(j)**2/2-b1*r(j)
         wgtba1(l,j)=wgtba1(l,j)+b1*dtheta/2

c kap=0 limiting value:
      wgtba1(1,ldb/2+1)=pi/gpar
     .        *sqrt(bmin*dtheta**2/(b_mag(ldb/2)-bmin))

c
c     normalize to epsn/sqrt(2*epsb) to match igeo=0 results
c     taub is normalized to q R_0/(v sqrt(2*epsb))
c
      do l=1,kd
         ntheta0=ldb/2+l
         ntheta1=ldb/2+2-l
         do j=ntheta1,ntheta0
            wgtba1(l,j)=wgtba1(l,j)*epsn*sqrt(2*epsb)
         enddo
      enddo
         
      do l=kd+1,kdpass
         do j=1,ldb
            wgtba1(l,j)=wgtba1(l,j)*epsn*sqrt(2*epsb)
         enddo
      enddo

c
c ntheta0,ntheta1 are chosen so that:
c For trapped particles, the bounce-average symmetrically involves 
c phi points above and below the midplane.
c For passing particles, the bounce-average involves all phi points
c from l=1 to ldb (the ld point is identical to the l=1 point).
c

c calculate bounce-time numerically (for passing particles also...):
      do l=1,kd
         ntheta0=ldb/2+l
         ntheta1=ldb/2+2-l
         do j=ntheta1,ntheta0
            taub(l)=taub(l)+wgtba1(l,j)
         enddo   
      enddo
      do l=kd+1,kdpass
         taub(l)=0.
         do j=1,ldb
            taub(l)=taub(l)+wgtba1(l,j)
         enddo   
      enddo

c calculate omegade numerically (for passing particles also...):

cffp$ skip R
      i=1
      do m=1,md
         do n=1,nd
            do l=1,kd
               omegade(l,m,n)=0.
               ntheta0=ldb/2+l
               ntheta1=ldb/2+2-l
               do j=ntheta1,ntheta0
                  omegade(l,m,n)=omegade(l,m,n)+wgtba1(l,j)
     .                 *omegad(j,m,n,i)/rho(i)/vt(i)
               enddo   
               if ((m.eq.meq).or.((md.eq.1).and.(mrr(m,n).eq.0))) 
     .              omegade(l,m,n)=0.
               omegade(l,m,n)=omegade(l,m,n)/taub(l)
            enddo
            do l=kd+1,kdpass
               omegade(l,m,n)=0.
            enddo
         enddo
      enddo
cffp$ skip R 

c debug:
      do l=1,kdpass
         write(6,112) l,kap(l),taub(l),omegade(l,2,1)/epsn/rky(2,1)
 112     format( 'l,kap(l),taub(l),omegade(l,2,1)/rky(l,2,1)/epsn=',i3,3f14.7) 
      enddo

      do l=1,kdpass
         ellk(l)=taub(l)/4
      enddo

c
c ?? I think I've now generalized all of this subroutine to use
c a grid-centered approach, with bounce-averages symmetrically 
c involving phi above and below the midplane...
c
      do j=1,ld
         do l=1,kdpass
            wgtba2(l,j)=0.
         enddo
      enddo

      do j=ldb/2+1,ld
         l0=j-ldb/2
         do l=l0,kdpass
            if (l.le.kd) then
               a1=1-b_mag(j)/bmin*(1-2*epsb*(kap(l)+dkap(l))**2)
               b1=1-b_mag(j)/bmin*(1-2*epsb*kap(l)**2)
               if (a1.lt.0) then
                  write(6,*) 'error in wgtba2 j,l,a1: ',j,l,a1
                  a1=0
               endif
               if (b1.lt.0) then
                  write(6,*) 'error in wgtba2 j,l,b1: ',j,l,b1
                  b1=0
               endif
               wgtba2(l,j)=sqrt(a1)-sqrt(b1)
            else
               a1=1-b_mag(j)/bmin*(1-2*epsb*kap(l)**2)
               b1=1-b_mag(j)/bmin*(1-2*epsb*(kap(l)-dkap(l))**2)
               if (a1.lt.0) then
                  write(6,*) 'error in wgtba2 j,l,a1: ',j,l,a1
                  a1=0
               endif
               if (b1.lt.0) then
                  write(6,*) 'error in wgtba2 j,l,b1: ',j,l,b1
                  b1=0
               endif
               wgtba2(l,j)=sqrt(a1)-sqrt(b1)
            endif   
            wgtba2(l,j)=wgtba2(l,j)/dkap(l)
c B(theta)=1/(1+eps*cos(theta)) included to get particle conservation.
c symmetry, for theta<0:
c ?? GWH: also need to put a 1/B into the flux-surface averages...
            wgtba2(l,ldb-j+2)=wgtba2(l,j)
         enddo
      enddo   

c ?? GWH: The bounce-averages (wgtba1), kappa-averages (wgtba2), 
c and the form of the coefficients of the pitch-angle-scattering 
c operator, are related in various ways.  Some of these relations 
c are important for particle conservation, and so it may be possible 
c to derive some of these coefficients from each other in a way
c that preserves conservation properties exactly (or at least to
c within roundoff error).  For now, just stick to the above methods,
c as they should preserve the conservation properties to second
c order in the grid spacings...

c load matrix for poisson inversion
      do j=1,ldb
         l1=1+abs(j-ldb/2-1)
         do l=l1,kdpass
            if ((l.eq.l1).and.(l1.lt.kd)) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)/2
            elseif ((l.gt.l1).and.(l.lt.kd)) then
               wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j))/2
            elseif ((l.eq.kd).and.(l1.ne.kd)) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j)/2
            elseif ((l.eq.kd).and.(l1.eq.kd)) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)
            elseif (l.eq.kd+1) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j)/2
            elseif ((l.gt.kd+1).and.(l.lt.kdpass)) then
               wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j))/2
            elseif (l.eq.kdpass) then
               wt2(l,j)=dkap(l)*wgtba2(l,j)/2
            endif   
         enddo   
         do j1=1,ldb
            pmat(j,j1)=0.
            l0=1+max(abs(j-ldb/2-1),abs(j1-ldb/2-1))
            do l=l0,kdpass
               pmat(j,j1)=pmat(j,j1)+wt2(l,j)*wgtba1(l,j1)/taub(l)
            enddo
            temp(j,j1)=pmat(j,j1)
         enddo
      enddo
c
c     invert
c
      ifail=0
c      call f01aaE(temp,lz,ldb/2+1,pmat,lz,wrkinv,ifail)
c      not used


      endif   !(igeo if)

c initialization for electron collisions (ecollis)

c need to calculate a few bounce-averages:
c
c bratio = <(B-Bmin)/B>  approx.= <epse*(1-cos(theta))>
c taub = bounce-time*v (only need the relative kappa-dependence)
c
      if (igeo.eq.0) then

      do l=1,kd
         bratio(l)=0.0
         ntheta0=ldb/2+l
         ntheta1=ldb/2+2-l
         do j=ntheta1,ntheta0
            bratio(l)=bratio(l)+cos(r(j))*wgtba1(l,j)
         enddo
         bratio(l)=bratio(l)/taub(l)
      enddo
      bratio(1)=0.0
      do l=kd+1,kdpass
         bratio(l)=0.0
         do j=1,ldb
            bratio(l)=bratio(l)+wgtba1(l,j)*cos(r(j))
         enddo   
         bratio(l)=bratio(l)/taub(l)
      enddo
      do l=2,kdpass
         bratio(l)=epsb*(1.-bratio(l))
      enddo   

      else

      do l=1,kd
         bratio(l)=0.0
         ntheta0=ldb/2+l
         ntheta1=ldb/2+2-l
         do j=ntheta1,ntheta0
            bratio(l)=bratio(l)+(b_mag(j)-bmin)/b_mag(j)*wgtba1(l,j)
         enddo
         bratio(l)=bratio(l)/taub(l)
      enddo
      bratio(1)=0.0
      do l=kd+1,kdpass
         bratio(l)=0.0
         do j=1,ldb
            bratio(l)=bratio(l)+(b_mag(j)-bmin)/b_mag(j)*wgtba1(l,j)
         enddo   
         bratio(l)=bratio(l)/taub(l)
      enddo

      endif


c debug:
         do l=1,kdpass
            write(6,103) l,bratio(l),taub(l)
 103        format('l, bratio(l), taub(l), =',i3,2f14.7) 
         enddo

c
c The generic diffusion operator can be written in the form:
c
c     df/dt = J(x) d/dx ( D(x) df/dx)
c
c Where J(x)=1/V'(x) is a Jacobian which represents conservation properties.
c In our case, x is the the pitch-angle kappa.  The bounce-averaged
c collision operator we use was calculated by J.G. Cordey, Nucl. Fus. 16, 
c 499 (976).  The notation I use is described on p.69-70 of my thesis,
c "Fast Ion Studies of Ion Cyclotron Heating in the PLT Tokamak", 
c G.W. Hammett, Princeton University Ph.D. Disseration, 1986.
c
c Use conservative differencing formulas, including the volume elements
c in the coefficients:
c
         do l=2,kdpass-1
            jac(l)=1.0/kap(l)/taub(l)/ ((kap(l+1)-kap(l-1))/2)
         enddo
         jac(kdpass)=1.0/kap(l)/taub(l)/ ((kap(kdpass)-kap(kdpass-1)))

c deal with kap(1)=0.0, by using kap evaluated half-a-grid away:
         kapav=kap(2)/4.0
         jac(1)=1./kapav/taub(1)/ (kap(2)/2)

         do l=1,kdpass-1
            kapp=(kap(l)+kap(l+1))/2
c            kap2p=(kap(l)**2+kap(l+1)**2)/2
            kap2p=kapp**2
            taubp=(taub(l)+taub(l+1))/2
c            bratiop=(bratio(l)+bratio(l+1))/2
            bratiop=(bratio(l)**0.5+bratio(l+1)**0.5)/2
            bratiop=bratiop**2

c the reason we interpolate sqrt(bratio) instead of bratio, is that bratio
c goes like kap**2 near kap=0.  There are near-cancellations between
c the kap**2 and bratio terms, and one must be careful to interpolate
c both of them in the same way or the diffusion coefficient could go
c negative near kap=0...

            diff(l)=(1-2*epsb*kap2p)*taubp/kapp*
     .           (-bratiop+2*epsb*kap2p)/(kap(l+1)-kap(l))
         enddo
         diff(0)=0.0  ! guard cell value

c evaluate at kap=0 point (taub logarithmically --> infinity at kap=1):
c         diff(kd)=(1-2*epse*kap(kd)**2)*taub(kd)/kap(kd)*
c     .           (-bratio(kd)+2*epse*kap(kd)**2)/(1.0-kap(kd))

c no diffusion of particles beyond maximum kappa grid point:
         diff(kdpass)=0.0

c debug:
         write(6,*) '0, diff(0+0.5),jac =',0,diff(0)
         do l=1,kdpass
            write(6,*) 'l, diff(l+0.5),jac =',l,diff(l),jac(l)
         enddo

      endif   !(epse if)

c end of initialization (executed only on the first call)
c**********************************************************************
      close(unit=12)  ! close the namelist input file (after call geometry)
c
c Debug: Check the surface area calculation:
c     
c      do l=1,ld
c         unity(l,1,1)=1.0
c      enddo
c      do m=2,md
c         do l=1,ld
c            unity(l,m,1)=0.
c         enddo
c      enddo
cc
c      call volflux(unity,unity,surfarea,dummy)
c      if(iperiod.ne.2) then
c         l_left=1
c         l_right=ld-1
c      endif
cc
c      dvdrhon=0.
c      do l=l_left,l_right
c         dvdrhon=dvdrhon+jacobian(l)*2.*x0/float(ldb)*gradrho(l)
c      enddo
c      surfarea=surfarea*dvdrhon*2.*pi
c      write(*,*) 'surface area (itgc)= ',surfarea,' L_n**2'
c      stop
      return
      end