#include "fpreproc/f77_dcomplx.h"
!------------------------------------------------------------------------
      subroutine rfxqlo_sum(
     >   nshe,numen,numxsi,enerkev,xsi,xsi_trap,btime,
     >   deb,dec,dab,dac,io,ierr)
!
      use rfxqlo_private
!
      IMPLICIT NONE
c
c  evaluate quasilinear coefficients, summing over all nphi values over
c  all antennae...
c
c  see rfxqlo_rfoper
c
c  calling arguments:
c
      integer, intent(in) :: nshe       ! radial shell index -- this refers
                                        ! to the xplasma grid "id_xqlo" as
                                        ! passed to rfxqlo_rfgrids_set in its
                                        ! compound type argument

      integer, intent(in) :: numen      ! number of energy points
      integer, intent(in) :: numxsi     ! number of pitch points

      real*8, intent(in) :: enerkev(0:numen-1) ! energy grid, *** keV ***
      real*8, intent(in) :: xsi(1:numxsi) ! pitch grid (vpll/v, dimensionless)
      real*8, intent(in) :: xsi_trap    ! trapping pitch angle
                                        !  |xsi(j)| < xsi_trap => trapped

      real*8, intent(in) :: btime(1:numxsi) ! normalized bounce time:
c        (Integral dl/abs(v_parrallel))/(2 pi q Rmaj/v)
c
c  see rfxqlo_qlocon compound type for specification of Rmaj in (m)
c
c--------------
c  quasilinear coefficients -- summed over multiple calls to this routine!
c  caller must ZERO the sum prior to entering an rfxqlo_rfoper call loop!
C
      real*8, intent(out) :: deb(0:numen-1,1:numxsi)
      real*8, intent(out) :: dec(0:numen-1,1:numxsi)
      real*8, intent(out) :: dab(0:numen-1,1:numxsi)
      real*8, intent(out) :: dac(0:numen-1,1:numxsi)
c-----------------
c OUTPUT -- space allocated by caller
c           quasilinear diffusion coefficients
c **these are sums over all harmonics, antennae, nphi**

c df/dt = (1./qjac) ( d/dE ( deb df/dE + dec df/d(xsi) )
c               + d/d(xsi) ( dab df/d(xsi) + dac df/dE ) )

c--The discrete version of the general F.P. equation is given
c       in "Particle Conserving B.C.'s for Fpp", Oct. 17, 1983, gwh,
c       which is based on McCoy, Mirin, &Killeen, "FPPAC: A
c       Two-Dimensional Multispecies Nonlinear Fokker-Planck
c       Package", Computer Physics Communications, V.24, p. 37 (1981).
c
c--A few notes are in order here.  In this code, qjac has had 2 pi q R
c       factored out of it, so that qjac = sqrt(enerev) abs(xsi) btime.
c       Because f must be symmetric in the trapped region, we must
c       also have deb, & dab symmetric in the trapped region
c       while dec & dac must be anti-symmetric.
c       *enforcing symmetry is up to the caller*
c       *rfxqlo_rfoper does not know where the trapped-passing bdy is*
c       Outside the trapped region, there is no symmetry requirement.

c  all coefficients are dimensioned (1:numen,1:numxsi) **by caller**
c  real*8, dimension(:,:) :: deb,dec,dab,dac
c
c  additional inputs and controls:
c     Ai, Zi, harmonic...  (nharm component written here...)
c
      type(rfxqlo_qlocon) :: io
c
      integer, intent(out) :: ierr      ! return status, 0=OK
c
c-----------------------------------------
      logical :: irenorm
      integer :: iharm,iant,iphi,iant_save,iphi_save
      REAL*8, PARAMETER :: ZERO = 0.0d0
#ifdef __DEBUG
      real*8 zwk(4*max(numen,numxsi))
      character*10 shlab
      character*50 nflab
#endif
c-----------------------------------------
c
      deb=ZERO
      dec=ZERO
      dab=ZERO
      dac=ZERO
c
      call rfxqlo_all_renorm_check(ierr)
      if(ierr.ne.0) return
c
      iant_save=selected_antenna
      iphi_save=selected_nphi
c
      do iant=1,n_antenna
         selected_antenna=iant

         if(antenna_array(iant)%power.le.ZERO) cycle

         do iharm=antenna_array(iant)%nharm_min,
     >      antenna_array(iant)%nharm_max
            io%nharm=iharm
            do iphi=1,antenna_array(iant)%num_nphi
               selected_nphi=iphi
c
c  dmc -- for fppmod, the spectrum for the QLO is symmetrized, even if
c  the original antenna spectrum is not symmetric
cxxx               if(antenna_array(iant)%sym_nphi) then
c

               antenna_array(iant)%nphi(iphi) =
     >              - antenna_array(iant)%nphi(iphi)

czzz               call rfxqlo_output_debug('rfxqlo_sum(1)')

               call rfxqlo_rfoper(nshe,numen,numxsi,enerkev,
     >              xsi,xsi_trap,btime,
     >              deb,dec,dab,dac,io,ierr)

               antenna_array(iant)%nphi(iphi) =
     >              - antenna_array(iant)%nphi(iphi)

cxxx               endif

czzz               call rfxqlo_output_debug('rfxqlo_sum(2)')

               call rfxqlo_rfoper(nshe,numen,numxsi,enerkev,
     >            xsi,xsi_trap,btime,
     >            deb,dec,dab,dac,io,ierr)

            enddo
         enddo
      enddo
!
      selected_antenna = iant_save
      selected_nphi = iphi_save

#ifdef __DEBUG
      zwk=ZERO
      write(shlab,'(" shell=",i3)') nshe
      nflab='rfxqlo_rfoper debug'
      call r8_grf3f1(enerkev,deb,xsi,zwk,numen,numxsi,numen,
     >   'energy (plasma frame)','DEB coeff.'//shlab,'pitch',
     >   'keV','DEB',' ',nflab,
     >   rfx_sys%run_label,0)
      call r8_grf3f1(enerkev,dab,xsi,zwk,numen,numxsi,numen,
     >   'energy (plasma frame)','DAB coeff.'//shlab,'pitch',
     >   'keV','DAB',' ',nflab,
     >   rfx_sys%run_label,0)
      call r8_grf3f1(enerkev,dec,xsi,zwk,numen,numxsi,numen,
     >   'energy (plasma frame)','DEC coeff.'//shlab,'pitch',
     >   'keV','DEC',' ',nflab,
     >   rfx_sys%run_label,0)
      call r8_grf3f1(enerkev,dac,xsi,zwk,numen,numxsi,numen,
     >   'energy (plasma frame)','DAC coeff.'//shlab,'pitch',
     >   'keV','DAC',' ',nflab,
     >   rfx_sys%run_label,0)
#endif

      return
      end
!
!------------------------------------------------------------------------
      subroutine rfxqlo_rfoper(nshe,numen,numxsi,enerkev,
     >   xsi,xsi_trap,btime,
     >   deb,dec,dab,dac,io,ierr)
!
! dmc 22 Jan 1997:  starting updown asymmetry generalization
!  YofBr, ThofBr, omppofr vs. theta(star) functions now defined for
!  both the upper and lower half plane, in asymmetric cases; in
!  symmetric cases only the upper half plane is defined.
! dmc 23 Jan 1997:  OK -- new subroutines rfx_zoneres and rfx_thupdown have
!  been created to evaluate the updown asymmetric sensitive quantities
!  I have convinced myself that the code supplied is asymptotically
!  correct, i.e. asymmetric algorithm with symmetric data => same
!   results as symmetric algorithm with symmetric data
!   ** but ** these changes should be checked by an RF physics expert
!----------------
!
!  evaluates the RF contributions to the F.P. operator.
!
! This routine became good for essentially general geometry sometime
! after the below date.  However, it still does assume B = Btor = K/R.
! Also, some of its algorithms may fail for bean-shaped plasmas.
!
! 20-Sep-1988 GWH:  moving towards non-circular geometry, although still
! assuming B = Btor = K/R.
!
! dmc -- B = g(psi)/R now.
!
! 2 RFOPER calls take about 3.5-4.8 seconds with lookup of resonance
! position in non-circular geometry.
!
! I have worked fairly hard at optimizing this routine.
! At most one could get another 15% speedup if the fintqe1d calls were
! replaced by instantaneous inline code, or another 12% if the acos calls
! could be replaced by single square root calls.  I don't think it is
! worth trying to do anymore optimization.
!
! includes several new effects different from the standard
! version before Sept. 1985:
! 1. generalized to kpar .ne. 0.
! 2. improved handling of the singularities.
! 3. Inclusion of E_minus J_n+1 terms.
!
! 31-Oct-1987: GWH, I increased the speed of the this routine by almost
! a factor of 4 by taking calls to interpolate on a table of
! bessel function values out of the inner loop! RFOPER now takes 2.5 secs
! per call.  Eliminating the calls to bsjns only saves about .1 secs.
!
      use rfxqlo_private
!
      IMPLICIT NONE
c
c  calling arguments:
c
      integer, intent(in) :: nshe       ! radial shell index -- this refers
                                        ! to the xplasma grid "id_xqlo" as
                                        ! passed to rfxqlo_rfgrids_set in its
                                        ! compound type argument

      integer, intent(in) :: numen      ! number of energy points
      integer, intent(in) :: numxsi     ! number of pitch points

      real*8, intent(in) :: enerkev(0:numen-1) ! energy grid, *** KeV ***
      real*8, intent(in) :: xsi(1:numxsi) ! pitch grid (vpll/v, dimensionless)
      real*8, intent(in) :: xsi_trap    ! trapping pitch angle
                                        !  |xsi(j)| < xsi_trap => trapped

      real*8, intent(in) :: btime(1:numxsi) ! normalized bounce time:
c        (Integral dl/abs(v_parrallel))/(2 pi q Rmaj/v)
c
c  see rfxqlo_qlocon compound type for specification of Rmaj in (m)
c
c--------------
c  quasilinear coefficients -- summed over multiple calls to this routine!
c  caller must ZERO the sum prior to entering an rfxqlo_rfoper call loop!
C
      real*8, intent(inout) :: deb(0:numen-1,1:numxsi)
      real*8, intent(inout) :: dec(0:numen-1,1:numxsi)
      real*8, intent(inout) :: dab(0:numen-1,1:numxsi)
      real*8, intent(inout) :: dac(0:numen-1,1:numxsi)
c-----------------
c OUTPUT -- space allocated by caller
c           quasilinear diffusion coefficients
c **these are sums over all harmonics, antennae, nphi**

c df/dt = (1./qjac) ( d/dE ( deb df/dE + dec df/d(xsi) )
c               + d/d(xsi) ( dab df/d(xsi) + dac df/dE ) )

c--The discrete version of the general F.P. equation is given
c       in "Particle Conserving B.C.'s for Fpp", Oct. 17, 1983, gwh,
c       which is based on McCoy, Mirin, &Killeen, "FPPAC: A
c       Two-Dimensional Multispecies Nonlinear Fokker-Planck
c       Package", Computer Physics Communications, V.24, p. 37 (1981).
c
c--A few notes are in order here.  In this code, qjac has had 2 pi q R
c       factored out of it, so that qjac = sqrt(enerev) abs(xsi) btime.
c       Because f must be symmetric in the trapped region, we must
c       also have deb, & dab symmetric in the trapped region
c       while dec & dac must be anti-symmetric.
c       *enforcing symmetry is up to the caller*
c       *rfxqlo_rfoper does not know where the trapped-passing bdy is*
c       Outside the trapped region, there is no symmetry requirement.

c  all coefficients are dimensioned (1:numen,1:numxsi) **by caller**
c  real*8, dimension(:,:) :: deb,dec,dab,dac
c
c  additional inputs and controls:
c     Ai, Zi, harmonic...
c
      type(rfxqlo_qlocon), intent(in) :: io
c
      integer, intent(out) :: ierr      ! return status, 0=OK
c
c  dmc Feb 1993 -- replace IMSL routine BSJNS with netlib routine
c  rjbesl -- tested OK Feb 23 1993 on DECstation;
c  cf TFTR.92 68828Z90 et al

C============
      integer :: maxen,maxxsi,nharm,iant,iphi
      logical :: lsym
      integer :: jtrap1,jtrap2
      real*8 :: rmaj,ab,rzbeam,bpcorr
      real*8 :: enerev(0:numen-1)
      real*8 :: ck1,ck2
C============
C idecl:  explicitize implicit INTEGER declarations:
      REAL*8, PARAMETER :: ZERO = 0.0d0
      REAL*8, PARAMETER :: HALF = 0.5d0
      REAL*8, PARAMETER :: C4TH = 0.25d0
      REAL*8, PARAMETER :: C10TH = 0.1d0
      REAL*8, PARAMETER :: ONE = 1.0d0
      REAL*8, PARAMETER :: TWO = 2.0d0
      REAL*8, PARAMETER :: THREE = 3.0d0
      REAL*8, PARAMETER :: FOUR = 4.0d0
      REAL*8, PARAMETER :: CEPS3 = 1.0d-3
      REAL*8, PARAMETER :: CEPS6 = 1.0d-6
      REAL*8, PARAMETER :: CBIG3 = 1.0d3
      REAL*8, PARAMETER :: CBIG6 = 1.0d6
!
      INTEGER ir,j,i,ingot,ibes,isign
      real*8 zbesh,zbesh2,zwpbes,zibes
!
!  items removed from old fppmodr declared here:
!
      integer njxsist,nticrf
      real*8 xsistq
!
C============
C idecl:  explicitize implicit REAL declarations:
      REAL*8 zpi,omegrf,zrmax,omeg0,ckovm,zra,zbres,zrres,zkparres,
     > zpsi2,rfinteq1d,zpsi1,zeovc,zmaxk,zargmax,zdel,
     > emcoef,zcon1,a,zw,zwup,zvup2,zvmid2,zvup,zvmid,
     > zbup,zrtbup,zbmid,zrtbmid,x,w,zv2,zv,b,zsqrtb,x1,x2,arg1,
     > arg2,singx,sinint1,sinint2,arg,vparres,bres,rresds,omres,
     > zdyavg,singer,vperp2,vperp,zthdsra,zthdsrb,zomppa,zomppb,
     > zompp,znvp2ovo,zomdd1,zomdd2,zomdd,zt2inv,ztest,ztaub2,argum,
     > theta2,zfix,rfdiff,zsym
c
      real*8 zkperpa,zkperpb,zeplusa,zeplusb
      complex*16 zcrfpola,zcrfpolb
C============
      REAL*8 zBrequ ! lookup in ThofBr array to find Theta(zBrequ,r)
c
      real*8 znphi                      ! nphi ~ npll
      real*8 zfreq                      ! antenna frequency (Hz)
      integer :: ishe                   ! shell index
c
      COMPLEX*16 zkpcorr

c  subgrid information...

      integer :: ixmulti,ix,ivmulti,nxsi_sg,ne_sg,ii,jj,iii,jjj
      integer :: maxne_sg
      real*8, dimension(:), allocatable :: xsi_sg,en_sg,zfvol
      real*8 :: zwvv,zwvva,zvtot

#ifdef __DEBUG
      real*8 zwk(4*max(numen,numxsi))
      character*10 shlab
      real*8, dimension(:,:,:), allocatable :: zdeb,zdec,zdab,zdac
      real*8, dimension(:,:,:,:), allocatable :: zepla,zbesa,zkpca,zfixa
      integer lnout
      logical lchar
      character*50 nflab
#endif

      data zeovc /1.6021d-20/
cxx      data Emcoef /ONE/                 !1 = turn on, 0= turn off E_minus

c-----------------------------------------------
      maxen = numen-1
      maxxsi = numxsi
      enerev = CBIG3*enerkev            ! -> eV

      ! velocity space subgridding
      ivmulti=rfx_grid%nvmulti          ! must be an odd integer (1,3,...)
      zwvv = (ONE/ivmulti)/ivmulti      ! 1/(ivmulti**2)
      zwvva = (ONE/ivmulti)/(ivmulti/2 + 1) ! for edge...

      ! real space subgridding
      ixmulti=rfx_grid%nxmulti
      allocate(zfvol(ixmulti))
      zvtot=ZERO
      do ix=1,ixmulti
         ishe = (nshe-1)*ixmulti + ix
         zfvol(ix)=dv_qlo(ishe)
         zvtot=zvtot+zfvol(ix)
      enddo
      zfvol=zfvol/zvtot

      ne_sg=1+ivmulti*maxen
      nxsi_sg=1+ivmulti*(numxsi-1)
      allocate(en_sg(0:ne_sg-1),xsi_sg(nxsi_sg))
      maxne_sg=ne_sg-1

      ii=-1
      do i=0,maxen-1
         do iii=0,ivmulti-1
            ii=ii+1
            en_sg(ii)=enerev(i)+iii*(enerev(i+1)-enerev(i))/ivmulti
         enddo
      enddo
      en_sg(ii+1)=enerev(maxen)

      jj=0
      jtrap1=0
      do j=1,maxxsi-1
         if(jtrap1.eq.0) then
            if(xsi(j).gt.-abs(xsi_trap)) jtrap1=j
         endif
         if(xsi(j).lt.abs(xsi_trap)) jtrap2=j
         do jjj=0,ivmulti-1
            jj=jj+1
            xsi_sg(jj)=xsi(j)+jjj*(xsi(j+1)-xsi(j))/ivmulti
         enddo
      enddo
      xsi_sg(jj+1)=xsi(maxxsi)

      lsym=antenna_array(selected_antenna)%sym_nphi

#ifdef __DEBUG
      zwk=ZERO
      write(shlab,'(" shell=",i3)') nshe
      allocate(zdeb(0:ne_sg-1,nxsi_sg,ixmulti),
     >   zdec(0:ne_sg-1,nxsi_sg,ixmulti))
      zdeb=ZERO; zdec=ZERO
      allocate(zdab(0:ne_sg-1,nxsi_sg,ixmulti),
     >   zdac(0:ne_sg-1,nxsi_sg,ixmulti))
      zdab=ZERO; zdac=ZERO
      allocate(zepla(0:ne_sg-1,nxsi_sg,2,ixmulti),
     >   zbesa(0:ne_sg-1,nxsi_sg,2,ixmulti))
      zepla=ZERO; zbesa=ZERO
      allocate(zkpca(0:ne_sg-1,nxsi_sg,2,ixmulti),
     >   zfixa(0:ne_sg-1,nxsi_sg,2,ixmulti))
      zkpca=ZERO; zfixa=ZERO
#endif

      ck1 = 1.0d0/cmpm                  ! convert k_perp: /m -> /cm
      ck2 = ck1/vpsv                    ! convert E_plus: V/m -> stat-volt/cm

cdbg      call zdebout		! debug output

      ierr = 0

      nticrf = rfx_grid%nchi_in

      rmaj = io%rmaj*100.0d0            ! -> cm
      ab = io%ai
      rzbeam  = io%zi
      nharm = io%nharm
      bpcorr = antenna_array(selected_antenna)%bpcorr

      iant = selected_antenna
      iphi = selected_nphi

      zpi=abs(acos(-ONE))
      zfreq=antenna_array(selected_antenna)%freq
      omegrf=2*zpi*zfreq                !rf frequency (radians/sec)

      !----------------------------------------------------------------
      ! loop over ixmulti subshells...
      do 2000 ix = 1,ixmulti

      ishe = (nshe-1)*ixmulti + ix

      ir=2*ishe-1
      zRmax=RatBmin(ir)
      omeg0=9.579d9*rZbeam*bpcorr*gzc(ishe)/zRmax/Ab  !gyro-freq @ Bmin
      ckovm=9.579e11  !k(ergs/ev) over m(grams)
      xsistq=1-nharm*Omeg0/omegrf     !kpar=0 resonant pitch angle

      zra = rZbeam/Ab

c find position of the kpar=0 resonance layer:
      zBres=omegrf/nharm/rZbeam/9.579d3*Ab
      zRres=bpcorr*CBIG6*gzc(ishe)/zBres ! gzc is R*Bt w/ para/diamag. corr.
c find kparallel at resonance:
      znphi=antenna_array(selected_antenna)%nphi(selected_nphi)
      zkparres=znphi/zRres

      zPsi2=rZbeam*rfinteq1d(xb_qlo(ishe),Psi,ZERO,delxPsi,numPsi)
      zPsi1=rZbeam*rfinteq1d(xb_qlo(ishe-1),Psi,ZERO,delxPsi,numPsi)

      zPsi2=zPsi2/zeovc/rZbeam
      zPsi1=zPsi1/zeovc/rZbeam

!set up an array of bessel function values:
!first find the maximum value of kperp, and check that our Bessel
!function table covers necessary size:

      zmaxk=ZERO
      do 11 j=1,nticrf
 11   zmaxk=max(zmaxk,wavefields(iphi,iant)%skperp(ishe,j)*ck1)

c maximum argument of Bess function table (*1.5 for safety):
      zargmax=TWO*sqrt(TWO*ckovm/Ab*enerev(maxen))*zmaxk/omeg0
      if(zargmax .lt. ONE) zargmax=ONE
      call rfx_bestabl(zargmax)

c calculate a few constants outside of the major loop:
      zcon1=sqrt(2*ckovm/ab)
      a=omegrf/nharm/omeg0

c*********************************************
c  Loop over energy and pitch angle grid     *
c*********************************************
      do 1510 ii=0,maxne_sg

c calculate a few constants outside of the pitch angle loop:
c (the vax compiler does not do a good job of automatically
c optimizing these loops, so we need to help it manually.)

         zW=en_sg(ii)
         if(ii .ne. maxne_sg)then
            zWup=(en_sg(ii)+en_sg(ii+1))/TWO
         else
            zWup=zW
         endif

         zvup2=2*ckovm/ab*zWup
         zvmid2=2*ckovm/ab*zW
         zvup=sqrt(zvup2)
         zvmid=sqrt(zvmid2)
         zbup=C4TH*zvup2*(zkparres/nharm/omeg0)**2
         zrtbup=sqrt(zbup)
         zbmid=C4TH*zvmid2*(zkparres/nharm/omeg0)**2
         zrtbmid=sqrt(zbmid)

         i = ii/ivmulti                 ! energy index for deb & dec

         do 1500 jj=1,nxsi_sg

c***************************************
c   evaluate deb and dec	       *
c***************************************
!  energy diffusion terms are evaluated between grid points.
!       deb(i,j) is deb halfway between i and i+1.
!       dec(i,j) is dec halfway between i and i+1.
!	deb(i,1)=deb(i,maxxsi)=dec(i,1)=dec(i,maxxsi)=0.0
!	    because vperp**2=0 there.  This is enforced in this
!	    subroutine to ensure the proper boundary conditions.

            x=xsi_sg(jj)                ! pitch angle
            j = 1 + (jj-1+(ivmulti/2))/ivmulti ! pitch index for deb & dec

            if(i .eq. maxen)go to 1200
            if(j .eq. 1 .or. j .eq. maxxsi) go to 1200 ! see above comment

            W=zWup
            zv2=zvup2
            zv=zvup
            b=zbup
            zsqrtb=zrtbup

!handling the xsi singularity:
! the xsi singularity is integrable, the integral can be found
! in Gradshteyn and Ryzhik.

            if(jj .eq. 1)then
               x1=xsi_sg(1)
            else
               x1=(xsi_sg(jj)+xsi_sg(jj-1))/2
            endif
            if(jj .eq. nxsi_sg)then
               x2=xsi_sg(nxsi_sg)
            else
               x2=(xsi_sg(jj)+xsi_sg(jj+1))/2
            endif

            arg1=1-a*(1-x1**2)+b*(1-x1**2)**2
            arg2=1-a*(1-x2**2)+b*(1-x2**2)**2

c if we are not near the singularity then we can use a simple
c average value:
            if(arg1 .le. ZERO.and. arg2 .le. ZERO) goto 970 !no resonance
            if(arg1 .gt. C10TH.and. arg2 .gt. C10TH.and. abs(arg1-arg2)
     1         .lt. C10TH*max(arg1,arg2))then
               singx=sqrt(HALF*(ONE/arg1+ONE/arg2))/zv/nharm
               goto 942
            endif

!the following formula for the singularity breaks down for kpar .eq. 0
! use an alternate formula in that case:
            if(b/a .lt. CEPS6)go to 930

            if(arg1 .ge. ZERO)then
               sinint1=log(abs(2*zsqrtb*sqrt(arg1)-2*b*(1-x1**2)+a))
     1            /2/zsqrtb
            else
               sinint1=log(sqrt(abs(a**2-4*b)))/2/zsqrtb
            endif

            if(arg2 .ge. ZERO)then
               sinint2=log(abs(2*zsqrtb*sqrt(arg2)-2*b*(1-x2**2)+a))
     1            /2/zsqrtb
            else
               sinint2=log(sqrt(abs(a**2-4*b)))/2/zsqrtb
            endif

            go to 940
!this is the kpar .eq. 0 formula:
 930        arg1=1-a*(1-x1**2)
            sinint1=0
            if(arg1 .ge. ZERO)sinint1=sqrt(arg1)/a
            arg2=1-a*(1-x2**2)
            sinint2=ZERO
            if(arg2 .ge. ZERO)sinint2=sqrt(arg2)/a

 940        singx=abs(sinint1-sinint2)/abs(x)/abs(x2-x1)
            singx=singx/nharm/zv

! singx is the value of 1/|n vpar-kpar mu c/e|_res averaged
! over a xsi zone.
 942        if(singx .le. ZERO)go to 970 !no resonance in this xsi zone.

!loop over two possible roots:
            do 960 isign=-1,1,2
               arg=1-a*(1-x**2)+b*(1-x**2)**2

!allow for the possibility that a tangent resonance occurs somewhere
! in this zone, although the center of the zone may miss the resonance:
               if(arg .lt. ZERO)arg=ZERO
               vparres=zv*(isign*sqrt(arg)
     1            +sign(ONE,zkparres)*zsqrtb*(1-x**2))
! skip if vpar and x aren't of the same sign.  remember to symmetrize
! <Q> for trapped particles later.
               if(lsym.or.(j.lt.jtrap1).or.(j.gt.jtrap2)) then
                  zsym=ONE
                  if( (x .ge. ZERO.and. vparres .lt. ZERO) .or.
     1               (x .lt. ZERO.and. vparres .ge. ZERO)) go to 960
               else
                  zsym=HALF
               endif
               Bres=(omegrf-zkparres*vparres)/
     >            nharm/rZbeam/9.579d3*Ab
c position of the true Doppler-shifted resonance:
               Rresds=bpcorr*CBIG6*gzc(ishe)/Bres
               zBrequ=Bres
! Omega-ci at resonance:
               omres=9.579d3*rZbeam*Bres/Ab

!Calculate d Omega / dl singularity:
! (a factor of Rq has been removed to cancel an Rq in BTIME)

c  (dmc) calculate up-down averaged resonance line/zone intersection
c  length:
               call rfx_zoneres(ishe,zBrequ,zdyavg)
               if(zdyavg.le.ZERO) go to 960 !resonance didn't intersect

               singer=zdyavg*
     1            Bres*Rresds**2/Omres/abs(zPsi2-zPsi1)/
     2            Rmaj/qzc(ishe)

! singer is now 1/|d Omega / dl| averaged over a radial zone,
! with Rmaj q divided out.

!vperp at resonance:
               vperp2=zv2-vparres**2
               if(vperp2 .lt. ZERO)vperp2=ZERO
               vperp=sqrt(vperp2)

c theta where Doppler-shifted resonance intersects this radial shell
c (set theta=0 or pi if the resonance intersects this radial zone, but
c not the zone center!):
c NOTE: THIS ALGORITHM WILL NOT WORK PROPERLY IF THE PLASMA IS INDENTED
C OR SUFFICIENTLY HIGH BETA SUCH THAT BMAX IS NOT AT THETA+PI!
c (dmc) -- attempting updown asymmetry generalization to overcome
c  this limitation (23 Jan 1997)

c return (angles btw 0 and 2pi):
c  zthdsra, zomppa: th, d**2 Omega / dl**2 @resonance above Bmin
c  zthdsrb, zomppb: th, d**2 Omega / dl**2 @resonance below Bmin

               call rfx_thupdown(ishe,zra,
     1            zpi,zBrequ,zthdsra,zthdsrb,zomppa,zomppb)

c Defining the correlation time for a single resonance as Tau_c, where
c Tau_c**2 is given in Eq. 2.30 of my thesis:
c
c	2 Tau_c**2 = | \integral d \tau \exp( - i \integral d \tau'
c		    (omega - k_\| v_\| - n \Omega) )  |
c
c The factor of 2 accounts for the fact that there are two resonances
c per bounce.  When the two resonances are well separated so the
c asymptotic form of the Airy function can be used, this reduces to
c Eq. 2.31 of my thesis.
c
c In terms of the variables used in this code, this becomes
c
c 	Tau_c**2 = 2 * Pi *singx * singer * Rmaj * qzc
c
c we now make sure Tau_c is less than the cutoff value given
c by the Airy function.  In my thesis (p. 49), I suggest using the
c maximum value of the Airy function.  Here, I will instead use the value of
c the Airy function at x=0 (where x is the argument of the Airy function.)
c This reduces the correlation time a bit more.

c zompp = d**2 Omega / d l**2:  up-down average
               zompp=HALF*(zomppa+zomppb)
c zomdd = - d**2 (n Omega + kparallel vparallel) / d t**2 at resonance:
cdmc	zomdd=-isign*vparres/singx*zompp+
cdmc	1   (1./singer/rmaj/qzc(ishe))**2 *nharm*vperp2/2./omres
c dmc -- recoded to avoid overflow
               znvp2ovo=nharm*vperp2/TWO/omres
               zomdd1=-isign*vparres/singx*zompp
               zomdd2=(ONE/singer/rmaj/qzc(ishe))**2
               zomdd=zomdd1 + zomdd2*znvp2ovo
c end dmc
               zt2inv=(abs(zomdd))**(TWO/THREE)/7.9d0*2
               ztest=zt2inv*2*zpi*singx*singer*Rmaj*qzc(ishe)
               if(ztest .gt. ONE) then
                  singx=singx/ztest
               endif

c
c we now make sure that (Tau_c)**2 < tau_bounce**2/4.
c We enforce this by reducing singx if necessary.
c
               ztaub2=(btime(j)*2*zpi*qzc(ishe)*Rmaj)**2 /zv2
               ztest=FOUR*2*zpi*singx*singer*Rmaj*qzc(ishe)/ztaub2
               if(ztest .gt. ONE) then
                  singx=singx/ztest
               endif

c dmc:  following index evaluations now correct for updown asymmetry:
c corresponding (nearest!) theta index, above midplane:
               call wfinterp(zthdsra,zkperpa,zeplusa,zcrfpola)
c and below midplane:
               call wfinterp(zthdsrb,zkperpb,zeplusb,zcrfpolb)

c--Kennel-Engelmann quasilinear diffusion coeficient is proportional
c       to | theta | **2.  Here we ignore the Epar contibution
c       and just keep Eplus and Eminus.

c resonance above the midplane:
c--Bessel function argument:
               argum=zkperpa*ck1*vperp/omres
               zibes=argum/delbes
               ibes=zibes
               if(ibes .ge. nmaxbes) then
                  write(lunmsg,*)
     1               'RFOPER warning, argument outside Bessel table!'
                  ibes=nmaxbes-1
                  zwpbes=ONE
               else
                  zwpbes=zibes-ibes
               endif
               zbesh=besa(ibes,nharm)*(ONE-zwpbes)+
     1            besa(ibes+1,nharm)*zwpbes
               zbesh2=besa(ibes,nharm+2)*(ONE-zwpbes)+
     1            besa(ibes+1,nharm+2)*zwpbes

               zkpcorr=zbesh+zcrfpola*zbesh2
               theta2=zeplusa*zeplusa*ck2*ck2*
     1            (AREAL(zkpcorr)*AREAL(zkpcorr)
     1            +AIMAG(zkpcorr)*AIMAG(zkpcorr))
#ifdef __DEBUG
               zbesa(ii,jj,1,ix)=zbesh
               zepla(ii,jj,1,ix)=zeplusa
               zkpca(ii,jj,1,ix)=
     1            (AREAL(zkpcorr)*AREAL(zkpcorr)
     1            +AIMAG(zkpcorr)*AIMAG(zkpcorr))
#endif

c resonance below the midplane:
c--Bessel function argument:
               argum=zkperpb*ck1*vperp/omres
               zibes=argum/delbes
               ibes=zibes
               if(ibes .ge. nmaxbes) then
                  write(lunmsg,*)
     1               'RFOPER warning, argument outside Bessel table!'
                  ibes=nmaxbes-1
                  zwpbes=ONE
               else
                  zwpbes=zibes-ibes
               endif
               zbesh=besa(ibes,nharm)*(ONE-zwpbes)+
     1            besa(ibes+1,nharm)*zwpbes
               zbesh2=besa(ibes,nharm+2)*(ONE-zwpbes)+
     1            besa(ibes+1,nharm+2)*zwpbes

               zkpcorr=zbesh+zcrfpolb*zbesh2
               theta2=theta2+zeplusb*zeplusb*ck2*ck2*
     1            (AREAL(zkpcorr)*AREAL(zkpcorr)
     1            +AIMAG(zkpcorr)*AIMAG(zkpcorr))
#ifdef __DEBUG
               zbesa(ii,jj,2,ix)=zbesh
               zepla(ii,jj,2,ix)=zeplusa
               zkpca(ii,jj,2,ix)=
     1            (AREAL(zkpcorr)*AREAL(zkpcorr)
     1            +AIMAG(zkpcorr)*AIMAG(zkpcorr))
#endif

               zfix=vperp2*singx*singer ! dmc -- avoid overflow, 20 Oct 93
               rfdiff=zcon1*zsym*W*abs(x)/FOUR*8.9873d4
     1            *rZbeam**2
     1            *zfix*theta2          ! dmc -- avoid overflow, 20 Oct 93
! 8.9873e4 = (e/k)**2, e is electron charge (cgs), k is erg/ev

               deb(i,j)=deb(i,j)+rfdiff*zwvv*zfvol(ix)
               dec(i,j)=dec(i,j)-rfdiff*zwvv*zfvol(ix)*
     1            (x**2-xsistq)/2/x/W
#ifdef __DEBUG
               zfixa(ii,jj,1,ix)=zfix
               zdeb(ii,jj,ix)=rfdiff
               zdec(ii,jj,ix)=rfdiff*(x**2-xsistq)/2/x/W
#endif
 960        continue
 970        continue

c***************************************
c   evaluate dab and dac	       *
c***************************************
!  angle diffusion term is evaluated halfway between grid points.
!       dab(i,j) is dab halfway between j and j+1:
!       dac(i,j) is dac halfway between j and j+1:
!  We have found that setting dac(i,maxen)=0.0
!  gives better numerical stability, so we enforce this condition
!  below.

 1200       continue
            i = (ii + (ivmulti/2))/ivmulti ! energy index for dab & dac
            j = 1 + (jj-1)/ivmulti      ! pitch index for dab & dac
            if(j .eq. maxxsi .or.
     1         j .eq. maxxsi/2 .or. i .eq. 0) go to 1500
            W=zW
            zv2=zvmid2
            zv=zvmid
            b=zbmid
            zsqrtb=zrtbmid
            x=HALF*(xsi_sg(jj)+xsi_sg(jj+1))

!handling the xsi singularity:

            x1=xsi_sg(jj)
            x2=xsi_sg(jj+1)

            arg1=1-a*(1-x1**2)+b*(1-x1**2)**2
            arg2=1-a*(1-x2**2)+b*(1-x2**2)**2

c if we are not near the singularity then we can use a simple
c average value:
            if(arg1 .le. ZERO.and. arg2 .le. ZERO) goto 1270 !no resonance
            if(arg1 .gt. C10TH.and. arg2 .gt. C10TH.and. abs(arg1-arg2)
     1         .lt. C10TH*max(arg1,arg2))then
               singx=sqrt(HALF*(ONE/arg1+ONE/arg2))/zv/nharm
               goto 1242
            endif

!the following formula for the singularity breaks down for kpar .eq. 0
! use an alternate formula in that case:
            if(b/a .lt. CEPS6)go to 1230

            if(arg1 .ge. ZERO)then
               sinint1=log(abs(2*zsqrtb*sqrt(arg1)-2*b*(1-x1**2)+a))
     1            /2/zsqrtb
            else
               sinint1=log(sqrt(abs(a**2-4*b)))/2/zsqrtb
            endif

            if(arg2 .ge. ZERO)then
               sinint2=log(abs(2*zsqrtb*sqrt(arg2)-2*b*(1-x2**2)+a))
     1	    /2/zsqrtb
            else
               sinint2=log(sqrt(abs(a**2-4*b)))/2/zsqrtb
            endif

            go to 1240
!this is the kpar .eq. 0 formula:
 1230       arg1=1-a*(1-x1**2)
            sinint1=0
            if(arg1 .ge. ZERO)sinint1=sqrt(arg1)/a
            arg2=1-a*(1-x2**2)
            sinint2=ZERO
            if(arg2 .ge. ZERO)sinint2=sqrt(arg2)/a

 1240       singx=abs(sinint1-sinint2)/abs(x)/abs(x2-x1)
            singx=singx/nharm/zv

! singx is the value of 1/|n vpar-kpar mu c/e|_res averaged
! over a xsi zone.
 1242       if(singx .le. ZERO)go to 1270 !no resonance in this xsi zone.

!loop over two possible roots:
            do 1260 isign=-1,1,2
               arg=1-a*(1-x**2)+b*(1-x**2)**2

!allow for the possibility that a tangent resonance occurs somewhere
! in this zone, although the center of the zone may miss the resonance:
               if(arg .lt. ZERO)arg=ZERO
               vparres=zv*(isign*sqrt(arg)
     1            +sign(ONE,zkparres)*zsqrtb*(1-x**2))
! skip if vpar and x aren't of the same sign.  remember to symmetrize
! <Q> for trapped particles later.
               if(lsym.or.(j.lt.jtrap1).or.(j.gt.jtrap2-1)) then
                  zsym=ONE
                  if( (x .ge. ZERO.and. vparres .lt. ZERO) .or.
     1               (x .lt. ZERO.and. vparres .ge. ZERO)) go to 1260
               else
                  zsym=HALF
               endif
               Bres=(omegrf-zkparres*vparres)/
     >            nharm/rZbeam/9.579d3*Ab
c position of the true Doppler-shifted resonance:
               Rresds=bpcorr*CBIG6*gzc(ishe)/Bres
               zBrequ=Bres
! Omega-ci at resonance:
               omres=9.579d3*rZbeam*Bres/Ab

!Calculate d Omega / dl singularity:
! (a factor of Rq has been removed to cancel an Rq in BTIME)

c  (dmc) calculate up-down averaged resonance line/zone intersection
c  length:
               call rfx_zoneres(ishe,zBrequ,zdyavg)
               if(zdyavg.le.ZERO) go to 1260 !resonance didn't intersect

               singer=zdyavg*Bres*Rresds**2/
     1            Omres/abs(zPsi2-zPsi1)/Rmaj/qzc(ishe)

! singer is now 1/|d Omega / dl| averaged over a radial zone,
! with Rmaj q divided out.

!vperp at resonance:
               vperp2=zv2-vparres**2
               if(vperp2 .lt. ZERO)vperp2=ZERO
               vperp=sqrt(vperp2)

c theta where Doppler-shifted resonance intersects this radial shell
c (set theta=0 or pi if the resonance intersects this radial zone, but
c not the zone center!):
c NOTE: THIS ALGORITHM WILL NOT WORK PROPERLY IF THE PLASMA IS INDENTED
C OR SUFFICIENTLY HIGH BETA SUCH THAT BMAX IS NOT AT THETA+PI!
c (dmc) -- attempting updown asymmetry generalization to overcome
c  this limitation (23 Jan 1997)

c return (angles btw 0 and 2pi):
c  zthdsra, zomppa: th, d**2 Omega / dl**2 @resonance above Bmin
c  zthdsrb, zomppb: th, d**2 Omega / dl**2 @resonance below Bmin

               call rfx_thupdown(ishe,zra,
     1            zpi,zBrequ,zthdsra,zthdsrb,zomppa,zomppb)

c Defining the correlation time for a single resonance as Tau_c, where
c Tau_c**2 is given in Eq. 2.30 of my thesis:
c
c	2 Tau_c**2 = | \integral d \tau \exp( - i \integral d \tau'
c		    (omega - k_\| v_\| - n \Omega) )  |
c
c The factor of 2 accounts for the fact that there are two resonances
c per bounce.  When the two resonances are well separated so the
c asymptotic form of the Airy function can be used, this reduces to
c Eq. 2.31 of my thesis.
c
c In terms of the variables used in this code, this becomes
c
c 	Tau_c**2 = 2 * Pi *singx * singer * Rmaj * qzc
c
c we now make sure Tau_c is less than the cutoff value given
c by the Airy function.  In my thesis (p. 49), I suggest using the
c maximum value of the Airy function.  Here, I will instead use the value of
c the Airy function at x=0 (where x is the argument of the Airy function.)
c This reduces the correlation time a bit more.

c zompp = d**2 Omega / d l**2: up-down average
               zompp=HALF*(zomppa+zomppb)
c zomdd = - d**2 (n Omega + kparallel vparallel) / d t**2 at resonance:
cdmc	zomdd=-isign*vparres/singx*zompp+
cdmc	1   (1./singer/rmaj/qzc(ishe))**2 *nharm*vperp2/2./omres
c dmc -- recoded to avoid overflow
               znvp2ovo=nharm*vperp2/TWO/omres
               zomdd1=-isign*vparres/singx*zompp
               zomdd2=(ONE/singer/rmaj/qzc(ishe))**2
               zomdd=zomdd1 + zomdd2*znvp2ovo
c end dmc
               zt2inv=(abs(zomdd))**(TWO/THREE)/7.9d0*TWO
               ztest=zt2inv*2*zpi*singx*singer*Rmaj*qzc(ishe)
               if(ztest .gt. ONE) then
                  singx=singx/ztest
               endif

c
c we now make sure that (Tau_c)**2 < tau_bounce**2/4.
c We enforce this by reducing singx if necessary.
c
               ztaub2=(btime(j)*2*zpi*qzc(ishe)*Rmaj)**2 /zv2
               ztest=FOUR*2*zpi*singx*singer*Rmaj*qzc(ishe)/ztaub2
               if(ztest .gt. ONE) then
                  singx=singx/ztest
               endif

c dmc:  following index evaluations now correct for updown asymmetry:
c corresponding (nearest!) theta index, above midplane:
               call wfinterp(zthdsra,zkperpa,zeplusa,zcrfpola)
c and below midplane:
               call wfinterp(zthdsrb,zkperpb,zeplusb,zcrfpolb)

c--Kennel-Engelmann quasilinear diffusion coeficient is proportional
c       to | theta | **2.  Here we ignore the Epar contibution
c       and just keep Eplus and Eminus.

c resonance above the midplane:
c--Bessel function argument:
               argum=zkperpa*ck1*vperp/omres
               zibes=argum/delbes
               ibes=zibes
               if(ibes .ge. nmaxbes) then
                  write(lunmsg,*)
     1               'RFOPER warning, argument outside Bessel table!'
                  ibes=nmaxbes-1
                  zwpbes=ONE
               else
                  zwpbes=zibes-ibes
               endif
               zbesh=besa(ibes,nharm)*(ONE-zwpbes)+
     1            besa(ibes+1,nharm)*zwpbes
               zbesh2=besa(ibes,nharm+2)*(ONE-zwpbes)+
     1            besa(ibes+1,nharm+2)*zwpbes

               zkpcorr=zbesh+zcrfpola*zbesh2
               theta2=zeplusa*zeplusa*ck2*ck2*
     1            (AREAL(zkpcorr)*AREAL(zkpcorr)
     1            +AIMAG(zkpcorr)*AIMAG(zkpcorr))

c resonance below the midplane:
c--Bessel function argument:
               argum=zkperpb*ck1*vperp/omres
               zibes=argum/delbes
               ibes=zibes
               if(ibes .ge. nmaxbes) then
                  write(lunmsg,*)
     1               'RFOPER warning, argument outside Bessel table!'
                  ibes=nmaxbes-1
                  zwpbes=ONE
               else
                  zwpbes=zibes-ibes
               endif
               zbesh=besa(ibes,nharm)*(ONE-zwpbes)+
     1            besa(ibes+1,nharm)*zwpbes
               zbesh2=besa(ibes,nharm+2)*(ONE-zwpbes)+
     1            besa(ibes+1,nharm+2)*zwpbes

               zkpcorr=zbesh+zcrfpolb*zbesh2
               theta2=theta2+zeplusb*zeplusb*ck2*ck2*
     1            (AREAL(zkpcorr)*AREAL(zkpcorr)
     1            +AIMAG(zkpcorr)*AIMAG(zkpcorr))

               zfix=vperp2*singx*singer ! dmc -- avoid overflow, 20 Oct 93
               rfdiff=zcon1*zsym*W*abs(x)/FOUR*8.9873d4
     1            *rZbeam**2
     1            *zfix*theta2          ! dmc -- avoid overflow, 20 Oct 93
! 8.9873e4 = (e/k)**2, e is electron charge (cgs), k is erg/ev

               if(i .lt. maxen) then
                  dac(i,j)=
     1               dac(i,j)-rfdiff*zwvv*zfvol(ix)*(x**2-xsistq)/2/x/w
                  dab(i,j)=
     1               dab(i,j)+(x**2-xsistq)**2*rfdiff*zwvv*zfvol(ix)/4/
     1               x**2/w**2
               else
                  dab(i,j)=
     1               dab(i,j)+(x**2-xsistq)**2*rfdiff*zwvva*zfvol(ix)/4/
     1               x**2/w**2
               endif
#ifdef __DEBUG
               zfixa(ii,jj,2,ix)=zfix
               zdab(ii,jj,ix)=(x**2-xsistq)**2*rfdiff/4/x**2/w**2
               if(i.lt.maxen) then
                  zdac(ii,jj,ix)=-rfdiff*(x**2-xsistq)/2/x/w
               endif
#endif
 1260       continue
 1270       continue

 1500    continue
 1510 continue

 2000 continue  !  ** end of subshell (ixmulti) loop

#ifdef __DEBUG
      ix=1
      nflab='rfxqlo_rfoper debug'
      call fpgrlab_nf(nflab)
      do
         write(lnout(0),*) trim(nflab)
         write(lnout(0),9901) ix,ixmulti
 9901    format(/' debug plot options:'/
     >      ' 1: deb contribution'/' 2: dec contribution'/
     >      ' 3: dab contribution'/' 4: dac contribution'/
     >      ' 5: Bessel fcn eval at resonance'/
     >      ' 6: mod(Bessel + crfpol corr).'/
     >      ' 7: Eplus  at resonance'/
     >      ' 8: singularity correction at resonance'/
     >      ' X: reset subshell index: ',i3,' (btw 1 and ',i3,')'/
     >      ' Q: quit')
         call uread('rfxqlo_rfoper: debug plot option:$')
         if(lchar('Q')) then
            exit
         else if(lchar('X')) then
            call ureadi('rfxqlo_rfoper: debug: new subshell index:$',ix)
            ix=max(1,min(ixmulti,ix))
         else if(lchar('1')) then
            call r8_grf3f1(en_sg*CEPS3,zdeb(0,1,ix),xsi_sg,zwk,
     >         ne_sg,nxsi_sg,ne_sg,
     >         'energy (plasma frame)','DEB contrib.'//shlab,'pitch',
     >         'keV','DEB',' ',nflab,
     >         rfx_sys%run_label,0)
         else if(lchar('2')) then
            call r8_grf3f1(en_sg*CEPS3,zdec(0,1,ix),xsi_sg,zwk,
     >         ne_sg,nxsi_sg,ne_sg,
     >         'energy (plasma frame)','DEC contrib.'//shlab,'pitch',
     >         'keV','DEC',' ',nflab,
     >         rfx_sys%run_label,0)
         else if(lchar('3')) then
            call r8_grf3f1(en_sg*CEPS3,zdab(0,1,ix),xsi_sg,zwk,
     >         ne_sg,nxsi_sg,ne_sg,
     >         'energy (plasma frame)','DAB contrib.'//shlab,'pitch',
     >         'keV','DAB',' ',nflab,
     >         rfx_sys%run_label,0)
         else if(lchar('4')) then
            call r8_grf3f1(en_sg*CEPS3,zdac(0,1,ix),xsi_sg,zwk,
     >         ne_sg,nxsi_sg,ne_sg,
     >         'energy (plasma frame)','DAC contrib.'//shlab,'pitch',
     >         'keV','DAC',' ',nflab,
     >         rfx_sys%run_label,0)
         else if(lchar('5')) then
            call r8_grf3f2(en_sg*CEPS3,zbesa(0,1,1,ix),zbesa(0,1,2,ix),
     >         xsi_sg,zwk,ne_sg,nxsi_sg,ne_sg,
     >         'energy (plasma frame)','Bessel @reson. (up)'//shlab,
     >         'Bessel @reson. (down)'//shlab,'pitch',
     >         'keV',' ',' ',nflab,
     >         rfx_sys%run_label,0)
         else if(lchar('6')) then
            call r8_grf3f2(en_sg*CEPS3,zkpca(0,1,1,ix),zkpca(0,1,2,ix),
     >         xsi_sg,zwk,ne_sg,nxsi_sg,ne_sg,
     >         'energy (plasma frame)','Bes + crfpol (up)'//shlab,
     >         'Bes + crfpol (down)'//shlab,'pitch',
     >         'keV',' ',' ',nflab,
     >         rfx_sys%run_label,0)
         else if(lchar('7')) then
            call r8_grf3f2(en_sg*CEPS3,zepla(0,1,1,ix),zepla(0,1,2,ix),
     >         xsi_sg,zwk,ne_sg,nxsi_sg,ne_sg,
     >         'energy (plasma frame)','Eplus @reson. (up)'//shlab,
     >         'Eplus @reson. (down)'//shlab,'pitch',
     >         'keV',' ',' ',nflab,
     >         rfx_sys%run_label,0)
         else if(lchar('8')) then
            call r8_grf3f2(en_sg*CEPS3,zfixa(0,1,1,ix),zfixa(0,1,2,ix),
     >         xsi_sg,zwk,ne_sg,nxsi_sg,ne_sg,
     >         'energy (plasma frame)','sing fix (deb/dec)'//shlab,
     >         'sing fix (dab/dac)'//shlab,'pitch',
     >         'keV',' ',' ',nflab,
     >         rfx_sys%run_label,0)
         else
            call uermsg('?rfxqlo_rfoper: unrecognized option.',' ',0)
         endif
      enddo
#endif

c--------------------------------------------
c  exit

 9999 continue
#ifdef __DEBUG
      deallocate(zdeb,zdec,zdab,zdac,zepla,zbesa,zkpca,zfixa)
#endif
      deallocate(xsi_sg,en_sg,zfvol)
      return
c--------------------------------------------
      contains
         subroutine wfinterp(zth,zkperp,zeplus,zcrfpol)
         !  linearly interpolate waveform data in poloidal direction
         real*8, intent(in) :: zth      ! poloidal angle coordinate
         real*8, intent(out) :: zkperp  ! kperp
         real*8, intent(out) :: zeplus  ! |E_plus|
         complex*16, intent(out) :: zcrfpol ! polarization

         real*8 zindx,zincr
         integer iindx

         zindx = ONE+zth*nticrf/TWO/zpi
         iindx = min(nticrf,int(zindx))
         zincr = zindx-iindx

         zkperp = (ONE-zincr)*wavefields(iphi,iant)%skperp(ishe,iindx)+
     >      zincr*wavefields(iphi,iant)%skperp(ishe,iindx+1)
         zeplus = (ONE-zincr)*wavefields(iphi,iant)%eplus(ishe,iindx)+
     >      zincr*wavefields(iphi,iant)%eplus(ishe,iindx+1)
         zcrfpol = (ONE-zincr)*wavefields(iphi,iant)%crfpol(ishe,iindx)+
     >      zincr*wavefields(iphi,iant)%crfpol(ishe,iindx+1)

         end subroutine wfinterp
      end
c-----------------------------------------------------------------
      subroutine rfx_bestabl(zargmax)

c  check/update Bessel function table
c  dmc -- do not remake this on every call to rfxqlo_rfoper...

      use rfxqlo_private
      implicit none
      intrinsic bessel_jn

      real*8, intent(in) :: zargmax

c-------------------------------
      real*8 :: zdel,zddum
      real*8, dimension(:), allocatable :: besn
      logical :: new_table
      integer :: i, ingot
      REAL*8, PARAMETER :: ZERO = 0.0d0
c-------------------------------

      new_table=.FALSE.
      if(.not.allocated(besa)) then
         new_table=.TRUE.
      else if(size(besa,2).lt.nharm_max+2) then
         new_table=.TRUE.
      else if(zargmax.gt.argmax_save) then
         new_table=.TRUE.
      endif

      if(new_table) then
         if(allocated(besa)) then
            deallocate(besa,arga)
         endif
      else
         return
      endif

c  allow margin, factor of 2, for future calls...
      argmax_save = max(argmax_save,2*zargmax)

      zdel=0.025d0
      nmaxbes=argmax_save/zdel + 0.5d0
      nmaxbes=max(100,nmaxbes)

      zdel=argmax_save/nmaxbes
      allocate(besa(0:nmaxbes,nharm_max+2),arga(0:nmaxbes))
      delbes=zdel

      allocate(besn(nharm_max+2))
      arga(0)=0
      besa(0,1)=1
      besa(0,2:nharm_max+2)=0
      do i=1,nmaxbes
         arga(i)=i*delbes
         zddum=arga(i)
         Besn(1:nharm_max+2) = BESSEL_JN(0, nharm_max+1, zddum)
         besa(i,1:nharm_max+2)=besn
      enddo

      deallocate(besn)

      return
      end
