Subroutine init
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     
      use itg_data
      use io
      use flr_terms
      use in, only: wdat
      use nl, only: init_nl
      implicit none
      include 'netcdf.inc'

c     netcdf variables
      integer status,ncid
      character filename*80

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

      integer i,l,m,n,itime(7),lp,lm,idum,ifail,j,l0,mdspec
      real z1,z2,z3,s21bcE,s21bbE,elliptic,kshift
      real theta,dtheta,dx,x,zkapmax,arg,arglim,dpsidrho
      real epsb
c     real, dimension(:,:), allocatable :: wt2
      real, dimension(:), allocatable :: ellk

c     Needed for debug (end of subroutine):
c     complex unity(lz,mz,nz)
c     real dummy(mz,nz),surfarea,dvdrhon
c     real, dimension(size(potential, 1)) :: wrkinv
c     real, dimension(size(potential, 1), size(potential, 1)) :: temp
      real ident
      integer ntheta0,ntheta1,j1,l1,k
      real kapav,bratiop,taubp,kapp,kap2p
      real zeta1,zeta2,dzeta
      real rkymin,rkymax
      integer ntim,ntimr,ndum1,ndum2
      real dt1,tim, rmu1_sav
      real a1,b1,c1,d1,theta1,theta2,bmin,b_max
      integer lin_sav,ne_sav,nfreq_sav,nstp_sav,ntotal_sav,ifilter_sav
      character*30 x05abE
      external x05abE
      external x05aaE
      external s21bcE
      external s21bbE
     
      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  /=  0 .and. nfftx  /=  0) then
         malias=nffty
         nalias=nfftx
         if(md  ==  0 .and. nd  ==  0) then
            md=(nffty-1)/3+1
            nd=2*((nfftx-1)/3)+1
         endif
      endif

      if(nd == -1) then
         nmin=0
         nmax=1
         nddm=1
         nddp=1
         allocate(nrr(1))
         nrr=1
         nd=1
         goto 413
      endif
     
      if(mod(nd,2) == 0) nd=nd+1
      if(nd == 1) then
         nmin=0
         nmax=0
         nddm=1
         nddp=1
         allocate(nrr(1))
         nrr=0
      else
         allocate(nrr(nd))
         nddm=nd/2
         nddp=nd/2+1
         nmin=-nd/2
         nmax=nd/2
         do n=1,nddp
            nrr(n)=n-1
         enddo
     
         do n=1,nddm
            nrr(n+nddp)=n-nddm-1
         enddo
      endif
     
 413  continue

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

      if(epse > 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
         
         allocate (ellk(kdpass))
         allocate (omegade(kdpass, md, nd))
         allocate (kap(kdpass+1), dkap(kdpass), taub(kdpass), bratio(kdpass))
         allocate (jac(kdpass), diff(0:kdpass))
      endif

      allocate (mrr(md, nd))
      call allocate 
      
      do n=1,nd
         if(md > 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   
            if(mrr(m,n) == 0.and.nrr(n) == 0) then
               meq=m
               neq=n
            endif
         enddo
      enddo

      write(9,wdat)
      write(9,1100) (nrr(n),mmin(n),mmax(n),n=1,nd)
 1100 format(3x,'at n =',i4,5x,'mmin =',i4,3x,'mmax =',i4)


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  ==  0) then
         malias=4*(md-1)        ! overkill (more efficient to specify nffty)
      endif
c     if(malias  >  mffty) call aborter(6,
c     &     'error: nffty or malias > mffty (see itg.in and itg_data.f90)')
      
      madd=malias-mdspec
      rootm=sqrt(float(malias))
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 == 1) nalias=1
      if(nalias  ==  0) then
         nalias=2*(nd-1)        ! overkill (more efficient to specify nfftx)
      endif
c     if(nalias  >  mfftx) call aborter(6,
c     &     'error: nfftx or nalias > mfftx (see itg.in and itg_data.f90)')
      
      nadd=nalias-nd
      rootn=sqrt(float(nalias))

      if(lin == 0) call init_nl
      
c     finish setting default m=0,n=0 mode if not already specified:
      if(meq  ==  -1) then
         neq=1
         meq=md+1
c     if(meq  >  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.)
     
      if(lin == 1) igradon=0
     
      avgflag=0.0
      navgflag=0
      ncycle = 0
      ne=0

      lhistory=float(ld)*.6+1
     
      ntotal = nstp/nfreq
     
c     r(l) is the grid for the radial coordinate:
     
      do l=1,ld
         r(l)=2.*x0*float(l-1)/ldb-x0
      enddo
     
      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) == 0.or.(lin == 1 .and. layout == 'linear')) then
               if (nmax /= 0) then
                  r0(m,n)=z0*float(nrr(n))/nmax
               else
                  r0(m,n)=z0*float(nrr(n))
               endif
            else
               if (nmax /= 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 == -1.) ifilter=0
      if(icrit == 2) then
         ifilter=1
         rmu1=max(rmu1,0.01)
      endif

      if(ifilter > 0) then
         allocate (rdiff(md, nd))
         rdiff=rmu1
      else
         rmu1=float(n_p)        !so it will appear in old output
      endif
c     
c     Prepare for theta FFT's
c     
      if (iperiod == 3) then 
         call connect
      else
         allocate (mcon(md, nd), ncon(1), scon(1))
         allocate (conpos(md, nd), icon(md, nd))
         ndomz=1
         nconz=1
         ndom=1
         ncon=md*nd
         scon=1
         conpos=1
         icon=1
         do m=1,md
            do n=1,nd
               mcon(m,n)=n+(m-1)*nd
            enddo
         enddo   
      endif  
      allocate (kpar(ld*nconz, ndomz))
      allocate (tablekp(ndomz, 101+2*ld*(nconz+1)))

c     
c     define ky grid
c     
      rky = mrr/y0
      rky2 = (mrr/y0)**2
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 == 1) then
         call geo(dpsidrho)
         bmaginv = 1./b_mag
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
      b_mag   = 1./(1.+eps*cos(r))
      bmaginv = 1./b_mag
      bgrad   = gpar*eps*sin(r)*b_mag
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 == meq .and. md /= 1) 
     .              rkx(l,m,n)=mrr(m+1,n)/y0*(-r0(m+1,n))*shr
               if (md == 1 .and. mrr(m,n) == 0) rkx(l,m,n)=1./y0
c     
c     local if qsf < 0
c     
               if (qsf < 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 == meq).or.((md == 1).and.mrr(m,n) == 0)) then
                     omegad(l,m,n,i)=epsn*
     .                    rho(i)*vt(i)*rkx(l,m,n)*sin(r(l))
	             if(iflow == 0) omegad(l,m,n,i)=0.
                  endif
c     
c     local if qsf<0
c     
                  if (qsf < 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     
      jacobian = abs(dpsidrho/(gpar*b_mag))

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 == 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 > 0.) then
         if (igeo == 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)
            b_max=b_mag(1)
            do l=1,ld
               if (b_mag(l) > b_max) b_max=b_mag(l)
               if (b_mag(l) < bmin) bmin=b_mag(l)
            enddo
            write(6,*) 'bmin,bmax: ',bmin,b_max
            
            epsb=(b_max-bmin)/2/b_max
            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  >  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) /= 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 == meq).or.((md == 1).and.mrr(m,n) == 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 < 0 .and. igeo == 0) 
     .                 omegade(l,m,n)=2.*epsn*rky(m,n)
               enddo
            enddo
         enddo
cfpp$ noskip R
      endif                     ! end of if(epse == 0.0)

c     
c     Twisted periodicity boundary conditions parameters:
c     
      if (iperiod == 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

! Initialize arrays with zeroes

      potential = 0.  ;  density = 0.  ;  u_par  = 0.  ;  t_par = 0.
      q_par     = 0.  ;  t_perp  = 0.  ;  q_perp = 0.

      nl_density1 = 0.  ;  nl_density2 = 0.
      nl_upar1    = 0.  ;  nl_upar2    = 0.  ;  nl_tpar     = 0.
      nl_qpar     = 0.  ;  nl_qperp1   = 0.  ;  nl_qperp2   = 0.
      nl_tperp1   = 0.  ;  nl_tperp2   = 0.  ;  nl_tperp3   = 0.  

      e_denk    = 0.   ;  phi_bk = 0.

      if(epse > 0.) then
         e_density = 0.  ;  e_p    = 0.  ;   e_r       = 0.
         e_t       = 0.  ;  phi_ba = 0.  
         nl_edensity = 0.  ;  nl_ep = 0.  ;  nl_er = 0.  ;  nl_et = 0.
      endif

      call allocate_nez(1+nstp/nprnt, nstp/nfreq)
      if(nread == 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

         ndum1 = 1  ! These arguments are used only by postc.
         ndum2 = 1  ! These arguments are used only by postc.

         if(binary) then
            open(unit=7,file=runname(1:lrunname)//'.bresp',
     .           form='unformatted',status='old')
            call bwread(7,ntim,ntimr,tim,dt1,ndum1,ndum2)
            close(unit=7)
         else
            filename=runname(1:lrunname)//'.ncp'
            status = nf_open(filename, 0, ncid) ! read only
            IF (status  /=  nf_noerr) 
     .           print *, filename, ' not found'
            call wreadnc(ncid,ntim,ntimr,tim,dt1,ndum1,ndum2)
            status = nf_close(ncid)
         endif
c     
c     restore B variation of kperp rho which was removed for post
c     
         do l=1,ld
            rkperp2(l,:,:)=rkperp2(l,:,:)*bmaginv(l)**2
         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 /= 0) then
            rkymin=rky(1,1)
            if(rkymin == 0.) then
               rkymin=1.
               if(md > 1) rkymin=rky(2,1)
            endif
            
            rkymax=rky(md,1)
            if(rkymax == 0.) rkymax=1.
            
            n=1
            do m=1,md
               time(m)=0.0
               dt(m)=rkymax/max(rkymin,rky(m,n))*dt0
            enddo
         else
            dt = dt0
         endif
      endif
c     
c     Allow for k_perp dependent etas for linear estimates:
c     
      allocate (etak(md, nd, nspecies), etak_par(md, nd, nspecies))
      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  ==  ' ' ) then
         date=x05abE(itime)
      else
         date=debug_plotlabel
      endif
c     
c     weights for bounce averages
c     
      if (epse > 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     
         allocate (wgtba1(kdpass, ld), wgtba2(kdpass, ld))
             
         if (igeo == 0) then
c     
c     circular concentric limit
c     
            wgtba1 = 0.

            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) < -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) > 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

            ellk = taub/4

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     
            wgtba2 = 0.

            do j=ldb/2+1,ld
               l0=j-ldb/2
               do l=l0,kdpass
                  if (l <= 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     BD: This is not used, so I commented it out.
c     load matrix for poisson inversion
c     do j=1,ldb
c     l1=1+abs(j-ldb/2-1)
c     do l=l1,kdpass
c     if ((l == l1).and.(l1 < kd)) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)/2
c     elseif ((l > l1).and.(l < kd)) then
c     wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j))/2
c     elseif ((l == kd).and.(l1 /= kd)) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j)/2
c     elseif ((l == kd).and.(l1 == kd)) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)
c     elseif (l == kd+1) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j)/2
c     elseif ((l > kd+1).and.(l < kdpass)) then
c     wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j))/2
c     elseif (l == kdpass) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)/2
c     endif   
c     enddo   
c     do j1=1,ldb
c     pmat(j,j1)=0.
c     l0=1+max(abs(j-ldb/2-1),abs(j1-ldb/2-1))
c     do l=l0,kdpass
c     pmat(j,j1)=pmat(j,j1)+wt2(l,j)*wgtba1(l,j1)/taub(l)
c     enddo
c     temp(j,j1)=pmat(j,j1)
c     enddo
c     enddo
c     
c     invert
c     
            ifail=0
c     call f01aaE(temp,ld,ldb/2+1,pmat,ld,wrkinv,ifail)
c     not used

         else
c     
c     general geometry bounce average weights
c     
            wgtba1 = 0.

            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)
c     c     .        -2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     c     .        *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)
c     c         wgtba1(l,j)=2/c1**2*(a1*c1*r(j)/3-2*a1*d1/3+b1*c1)
c     c     .           *sqrt(c1*r(j)+d1)
c     c     .           -2/c1**2*(a1*c1*(r(j)-dtheta/2)/3-2*a1*d1/3+b1*c1)
c     c     .           *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...):
c This compiler directive is probably ignored by f90??
cffp$ skip R
            omegade=0.
            do m=1,md
               do n=1,nd
                  do l=1,kd
                     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 == meq).or.((md == 1).and.(mrr(m,n) == 0))) 
     .                    omegade(l,m,n)=0.
                     omegade(l,m,n)=omegade(l,m,n)/taub(l)
                  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

            ellk = taub/4
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     
            wgtba2 = 0.

            do j=ldb/2+1,ld
               l0=j-ldb/2
               do l=l0,kdpass
                  if (l <= 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 < 0) then
                        write(6,*) 'error in wgtba2 j,l,a1: ',j,l,a1
                        a1=0
                     endif
                     if (b1 < 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 < 0) then
                        write(6,*) 'error in wgtba2 j,l,a1: ',j,l,a1
                        a1=0
                     endif
                     if (b1 < 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     BD: Since this is not used, I commented it out.
c     load matrix for poisson inversion
c     do j=1,ldb
c     l1=1+abs(j-ldb/2-1)
c     do l=l1,kdpass
c     if ((l == l1).and.(l1 < kd)) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)/2
c     elseif ((l > l1).and.(l < kd)) then
c     wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j))/2
c     elseif ((l == kd).and.(l1 /= kd)) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l-1)*wgtba2(l-1,j)/2
c     elseif ((l == kd).and.(l1 == kd)) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)
c     elseif (l == kd+1) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j)/2
c     elseif ((l > kd+1).and.(l < kdpass)) then
c     wt2(l,j)=(dkap(l)*wgtba2(l,j)+dkap(l+1)*wgtba2(l+1,j))/2
c     elseif (l == kdpass) then
c     wt2(l,j)=dkap(l)*wgtba2(l,j)/2
c     endif   
c     enddo   
c     do j1=1,ldb
c     pmat(j,j1)=0.
c     l0=1+max(abs(j-ldb/2-1),abs(j1-ldb/2-1))
c     do l=l0,kdpass
c     pmat(j,j1)=pmat(j,j1)+wt2(l,j)*wgtba1(l,j1)/taub(l)
c     enddo
c     temp(j,j1)=pmat(j,j1)
c     enddo
c     enddo
c     c
c     invert
c     
            ifail=0
c     call f01aaE(temp,ld,ldb/2+1,pmat,ld,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 == 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
c     c
c     call volflux(unity,unity,surfarea,dummy)
c     if(iperiod /= 2) then
c     l_left=1
c     l_right=ld-1
c     endif
c     c
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

      if(epse > 0.0) deallocate (ellk)
      ntimf = nfreq
      nfreqcnt = 0
      return
      end