subroutine orbrz_set_rf_info2

  use ezspline, only : ezspline_free,ezspline_init,ezspline_error,&
       ezspline_setup,ezspline_interp2,ezspline_gradient2,ezlinear_init
  use orbrzv
  use nbi_com,only: twopi,orbrzv_toric_prftot,orbrzv_toric_prf,nsbeam,lcentr,&
       &ledge,mj,epslon,orbrzv_nr,orbrzv_nz,xe,abeams,xzbeams,bmaxm,&
       &orbrzv_toric_frnm,aplasm,backz,orbrzv_toric_prftot0,orbrzv_prftot,&
       &orbrzv_prf,ebdmaxa,ng,bminm,nonlin,vmins,zel,zmp,zc
  use nbi_mpi_com
  use logmod
  !use rfxqlo_interface
  !use xplasma_obj_instance

  implicit none

  integer, parameter :: lup_max = 35

  !type(rfxqlo_ion) :: zions
  !type(rfxqlo_dims) :: zdims
  !type(rfxqlo_antenna) :: zant

  real(fp), allocatable, dimension(:) :: zx, zth 
  real(fp), allocatable, dimension(:,:) :: zkp1, zep1
  complex(c_double_complex), allocatable, dimension(:,:) :: zcrfpol1
  real(fp) :: orbrzv_tmp(16)
  complex(c_double_complex) :: zctmp

  integer ierr, iwarn, mpierr, inx, inth, ii, jj, bc_x(2), bc_th(2), &
        & ksym, iant, inphi, imode, inant, imaxnphi, iwfdim, ioffset
  integer :: idx_nx, idx_nth, idx_nant, idx_nnphi
  real(fp) :: zms, zqs, zvmax, zomg_max, zomg_min, zkvmax, zdx, zdth, &
            & zfreq, zfreq_min, zkl, zbdbq_m, zxe2cgs, zrho, ztheta, &
            & zr, zz, zpsi, zksym, zkvmin

  real(fp) :: zdpsi(1:2), zerf, zderf(1:2), zlrf_inv
  real(fp), allocatable, dimension(:,:) :: zlrfi

  real(fp) :: zrr, zdtdr, zdtdz, zdrhodr, zdrhodz, zlr, zlz, zdrhodpsi

  integer :: inions, j, isp
  real(fp) :: zdpwr, zpwr_max, zrtmp, zztmp, zpwrtot, zpwr0tot, zfrnm

  character*200 zfile
  integer :: ilunt, iopstat, ll, lup

  !--------------------------------------------------------------------------
  ! for each antenna and nphi, fetch wavefield data and prepare
  call infolog(' orbrz_set_rf2: start orbrz_set_rf_info')

#ifdef __MPI
  call mpi_bcast(mcrf_nbuf, 3, MPI_INTEGER, 0, nbi_comm, mpierr)
#endif
  if(mcrf_nbuf(1).le.0.or.mcrf_nbuf(2).le.0) then
     mcrf_on = .false. 
     call warnlog(' orbrz_set_rf2: mcrf is forced to be disabled')
     call warnlog('                mcrf_nbuf(1)=',mcrf_nbuf(1))
     call warnlog('                mcrf_nbuf(2)=',mcrf_nbuf(2))
     if(myidd.eq.0) then
       write(nonlin,*)' !orbrzv_set_rf2 : mcrf is forced to be disabled'
     endif
     return
  endif
  if(mcrf_nbuf(3).eq.0) then
     mcrf_on = .false. 
     call infolog(' orbrz_set_rf2: mcrf is disabled')
     if(myidd.eq.0) write(nonlin,*) ' %orbrzv_set_rf2 : mcrf is disabled'
     return
  else
     mcrf_on = .true.
  endif

#ifdef __MPI
  if(myidd.ne.0) then
     if(allocated(mcrf_ibuf)) deallocate(mcrf_ibuf)
     if(allocated(mcrf_rbuf)) deallocate(mcrf_rbuf)
     allocate(mcrf_ibuf(mcrf_nbuf(1)))
     mcrf_ibuf=0
     allocate(mcrf_rbuf(mcrf_nbuf(2)))
     mcrf_rbuf=0.0_fp
  endif
  
  call mpi_bcast(mcrf_ibuf, mcrf_nbuf(1), MPI_INTEGER, &
       & 0, nbi_comm, mpierr)
  call mpi_bcast(mcrf_rbuf, mcrf_nbuf(2), MPI_REAL8, &
       & 0, nbi_comm, mpierr)
#endif
  
  ! free previous allocations -----------------------------------------------
  if(allocated(orbrzv_rfw)) call dealloc_orbrzv_rfw(ierr)

  if(allocated(orbrzv_rf_nnphi)) deallocate(orbrzv_rf_nnphi)

  if(allocated(orbrzv_lmin))  deallocate(orbrzv_lmin)
  if(allocated(orbrzv_lmax))  deallocate(orbrzv_lmax)
  if(allocated(orbrzv_trmax)) deallocate(orbrzv_trmax)
  if(allocated(orbrzv_nrf))   deallocate(orbrzv_nrf)

  ! rfdbg_p* -------------------------------------------------
  if(allocated(rfdbg_pdiag))   deallocate(rfdbg_pdiag)
  ! rfdbg_p* -------------------------------------------------

!mg 11/16 cleanup  if(allocated(orbrzv_my_rfdist))     deallocate(orbrzv_my_rfdist)
!mg 11/16 cleanup  if(allocated(orbrzv_my_rfdist_tmp)) deallocate(orbrzv_my_rfdist_tmp)

  if(allocated(orbrzv_fbm_dvl)) deallocate(orbrzv_fbm_dvl)
  if(allocated(orbrzv_fbm_dvp)) deallocate(orbrzv_fbm_dvp)
  if(allocated(orbrzv_fbm)) deallocate(orbrzv_fbm)
  ! -------------------------------------------------------------------------

  ! set number of antenna
  !call rfxqlo_dim_get(zdims, ierr)
  !orbrzv_rf_nant = zdims%n_antenna
  !orbrzv_rf_max_nphi = zdims%max_nphi    
  orbrzv_rf_nant     = mcrf_ibuf(1)
  orbrzv_rf_max_nphi = mcrf_ibuf(2)
  inant    = mcrf_ibuf(1)
  imaxnphi = mcrf_ibuf(2)

  allocate(orbrzv_rf_nnphi(1:orbrzv_rf_nant))
  call alloc_orbrzv_rfw(0,0,0)
  orbrzv_rf_nnphi=0
  ! fetch wavefield array dimensions
  !call rfxqlo_wavefield_dims(inx, inth, ierr)
  !call rfxqlo_ion_get(zions, ierr)
  inx  = mcrf_ibuf(3)
  inth = mcrf_ibuf(4)
  inth = inth - 1

  ! number of ion species
  inions = mcrf_ibuf(5)

  ! setup for index function
  idx_nx  = inx
  idx_nth = inth+1
  idx_nant  = inant
  idx_nnphi = imaxnphi

  iwfdim = inx*idx_nth

  orbrzv_rf_nx  = inx    
  orbrzv_rf_nth = inth            
  zdx  = 1.0_fp/real(orbrzv_rf_nx,fp)
  zdth = twopi/real(orbrzv_rf_nth,fp)
  orbrzv_rf_dx  = zdx
  orbrzv_rf_dth = zdth

  !orbrzv_rf_nsp = zions%n_ions
  orbrzv_rf_nsp = mcrf_ibuf(5)

  !print *, "inx, inth ", inx, inth
  allocate(zx(1:inx+1), zth(1:inth+1))
  allocate(zkp1(inx+1,inth+1), zep1(inx+1,inth+1))
  allocate(zcrfpol1(inx+1,inth+1))
  zx=0.0_fp
  zth=0.0_fp
  zkp1=0.0_fp
  zep1=0.0_fp
  zcrfpol1=0.0_fp
  bc_x = (/ 0, 0 /)
  bc_th = (/ -1, -1 /)

  do ii = 1, inx+1
    zx(ii) = zdx*real(ii-1,fp)
  enddo

  do jj = 1, inth+1
    zth(jj) = zdth*real(jj-1,fp)
  enddo

  orbrzv_toric_prftot = 0.0_fp
  orbrzv_toric_prf = 0.0_fp

  ioffset = 0
  orbrzv_kp2_max = 0.0_fp

  ! now loop over antenna and nphi to initialize and setup wavefield
  do iant = 1, orbrzv_rf_nant

    orbrzv_rf_nnphi(iant) = mcrf_ibuf(5+iant)

    do inphi = 1, orbrzv_rf_nnphi(iant)

      orbrzv_rfw(iant,inphi)%freq  = mcrf_rbuf(ioffset+1)
      orbrzv_rfw(iant,inphi)%power = mcrf_rbuf(ioffset+2)
      orbrzv_rfw(iant,inphi)%nphi  = mcrf_rbuf(ioffset+3)
      orbrzv_rfw(iant,inphi)%kl    = mcrf_rbuf(ioffset+3)
      zksym                        = mcrf_rbuf(ioffset+4)
      if(zksym > 1.0_fp) then
        orbrzv_rfw(iant,inphi)%ksym  = 2
      else
        orbrzv_rfw(iant,inphi)%ksym  = 1
      endif
!mg 04/07/16 symmetry added
!mg1      orbrzv_rfw(iant,inphi)%ksym  = 2
      ioffset = ioffset + 4

      ! kperp,  1/m not kperp^2
      do jj = 1, inth+1
        do ii = 1, inx
          zkp1(ii+1,jj) = mcrf_rbuf(ioffset+mcrf_idxf(ii,jj))
        enddo
        zkp1(1,jj) = mcrf_rbuf(ioffset+mcrf_idxf(1,jj))
      enddo
      ioffset = ioffset + iwfdim

      ! e+
      do jj = 1, inth+1
        do ii = 1, inx
          zep1(ii+1,jj) = mcrf_rbuf(ioffset+mcrf_idxf(ii,jj))
        enddo
        zep1(1,jj) = mcrf_rbuf(ioffset+mcrf_idxf(1,jj))
      enddo
      ioffset = ioffset + iwfdim

      ! real part of e+/e-
      do jj = 1, inth+1
        do ii = 1, inx
          zcrfpol1(ii+1,jj) = dcmplx(mcrf_rbuf(ioffset+mcrf_idxf(ii,jj)))
        enddo
        zcrfpol1(1,jj) = dcmplx(mcrf_rbuf(ioffset+mcrf_idxf(1,jj)))
      enddo
      ioffset = ioffset + iwfdim

      ! imaginary part of e+/e-
      do jj = 1, inth+1
        do ii = 1, inx
          zctmp = dcmplx(0.0_fp, mcrf_rbuf(ioffset+mcrf_idxf(ii,jj)))
          zcrfpol1(ii+1,jj) = zcrfpol1(ii+1,jj) + zctmp
        enddo
        zctmp = dcmplx(0.0_fp, mcrf_rbuf(ioffset+mcrf_idxf(1,jj)))
        zcrfpol1(1,jj) = zcrfpol1(1,jj) + zctmp
      enddo
      ioffset = ioffset + iwfdim

      do isp = 1, nsbeam 
        do j = lcentr, ledge
        !do j = lcentr+1, ledge ! neglect centreal power peak !
          zdpwr = mcrf_rbuf(ioffset+mcrf_idxp(j,isp))
          orbrzv_toric_prf(j, isp) = orbrzv_toric_prf(j, isp) + zdpwr
          orbrzv_toric_prftot(isp) = orbrzv_toric_prftot(isp) + zdpwr
        enddo
      enddo
      ioffset = ioffset + mj*nsbeam

      ! mks to cgs
      ! m^-1 ==> cm^-1
      zkp1(1:inx+1,1:inth+1) = 1.0e-2_fp*zkp1(1:inx+1,1:inth+1)
      ! v/m  ==> statvolt/cm
      zep1(1:inx+1,1:inth+1) = 1.0e-4_fp*zep1(1:inx+1,1:inth+1)/3.0_fp

      ! find maximum kp2 for later use
      do ii = 1, inx+1
        do jj = 1, inth+1
          orbrzv_kp2_max = max(zkp1(ii,jj), orbrzv_kp2_max)
        enddo
      enddo

      ! initialize spline objects
      call ezlinear_init(orbrzv_rfw(iant,inphi)%kp_spo,inx+1,inth+1,ierr)
      call ezspline_error(ierr)

      !call ezspline_init(orbrzv_rfw(iant,inphi)%ep_spo, &
      !                 & inx+1,inth+1,bc_x,bc_th,ierr)
      call ezlinear_init(orbrzv_rfw(iant,inphi)%ep_spo,inx+1,inth+1,ierr)
      call ezspline_error(ierr)

      !call ezspline_init(orbrzv_rfw(iant,inphi)%rp_spo, &
      !                 & inx+1,inth+1,bc_x,bc_th,ierr)
      call ezlinear_init(orbrzv_rfw(iant,inphi)%rp_spo,inx+1,inth+1,ierr)
      call ezspline_error(ierr)

      !call ezspline_init(orbrzv_rfw(iant,inphi)%ip_spo, &
      !                 & inx+1,inth+1,bc_x,bc_th,ierr)
      call ezlinear_init(orbrzv_rfw(iant,inphi)%ip_spo,inx+1,inth+1,ierr)
      call ezspline_error(ierr)

      ! setup spline objects
      ! x1 = sqrt(psi/psix)
      orbrzv_rfw(iant,inphi)%kp_spo%x1(1:inx+1)  = zx(1:inx+1)
      orbrzv_rfw(iant,inphi)%ep_spo%x1(1:inx+1)  = zx(1:inx+1)
      orbrzv_rfw(iant,inphi)%rp_spo%x1(1:inx+1)  = zx(1:inx+1)
      orbrzv_rfw(iant,inphi)%ip_spo%x1(1:inx+1)  = zx(1:inx+1)
      ! x2 = theta
      orbrzv_rfw(iant,inphi)%kp_spo%x2(1:inth+1) = zth(1:inth+1)
      orbrzv_rfw(iant,inphi)%ep_spo%x2(1:inth+1) = zth(1:inth+1)
      orbrzv_rfw(iant,inphi)%rp_spo%x2(1:inth+1) = zth(1:inth+1)
      orbrzv_rfw(iant,inphi)%ip_spo%x2(1:inth+1) = zth(1:inth+1)

      call ezspline_setup(orbrzv_rfw(iant,inphi)%kp_spo, zkp1*zkp1, ierr)
      call ezspline_error(ierr)
      call ezspline_setup(orbrzv_rfw(iant,inphi)%ep_spo, zep1, ierr)
      call ezspline_error(ierr)

      zkp1(1:inx+1,1:inth+1) = real(zcrfpol1(1:inx+1,1:inth+1),fp)
      zep1(1:inx+1,1:inth+1) = dimag(zcrfpol1(1:inx+1,1:inth+1))
      call ezspline_setup(orbrzv_rfw(iant,inphi)%rp_spo, zkp1, ierr)
      call ezspline_error(ierr)
      call ezspline_setup(orbrzv_rfw(iant,inphi)%ip_spo, zep1, ierr)
      call ezspline_error(ierr)

    enddo ! inphi
  enddo ! iant

  ! find the species subject to maximum rf-heating
  ! sum total rf power in this time step and previous time step
  zpwrtot = 0.0_fp
  zpwr0tot = 0.0_fp
  zpwr_max = 0.0_fp
  do isp = 1, nsbeam
     call infolog(' orbrz_set_rf2: toric prf(isp),isp=',isp)
     call infolog('                toric prf(isp)=',orbrzv_toric_prftot(isp))
     zpwrtot = zpwrtot + orbrzv_toric_prftot(isp)
     zpwr0tot = zpwr0tot + orbrzv_toric_prftot0(isp)
  enddo

  ! adjust frnm by comparing zpwrtot & zpwr0tot
  !if(min(zpwrtot,zpwr0tot) >= 1.0e4_fp) then
  !  zfrnm = sqrt(zpwrtot/zpwr0tot)
  !  zfrnm = 0.5_fp*(1.0_fp + zfrnm)
  !  print *, ' %orbrzv_set_rf : frnm updated by rf-input power change :', &
  !        & orbrzv_toric_frnm, zfrnm*orbrzv_toric_frnm
  !  orbrzv_toric_frnm = zfrnm*orbrzv_toric_frnm
  !endif

  !if(zpwr_max < 1.0e4) then
  !  print *, ' %orbrzv_set_rf : too low power, skip this time mcrf'
  !  orbrzv_rf_option = 0
  !endif


  !done with mcrf_... burffer, free memory
  if(allocated(mcrf_ibuf)) deallocate(mcrf_ibuf)
  if(allocated(mcrf_rbuf)) deallocate(mcrf_rbuf)

  deallocate(zx, zth)
  deallocate(zkp1, zep1)
  deallocate(zcrfpol1)

  !--------------------------------------------------------------------------
  ! allocate min & max of allowable harmonic number for each species
  allocate(orbrzv_lmin(1:nsbeam), orbrzv_lmax(1:nsbeam))
  allocate(orbrzv_trmax(1:nsbeam))
  orbrzv_lmin=0
  orbrzv_lmax=0
  orbrzv_trmax=0.0_fp
  !!allocate(orbrzv_prftot(1:nsbeam)); 
  !!allocate(orbrzv_prf(1:nsbeam,1:mj)); 
  allocate(orbrzv_nrf(1:nsbeam,1:3,1:mj)); orbrzv_nrf = 0.0_fp

  orbrzv_prftot = 0.0_fp
  orbrzv_prf = 0.0_fp

!mg 11/16 cleanup  ! my_rfdist_diag
!mg 11/16 cleanup  orbrzv_my_nr = 10
!mg 11/16 cleanup  orbrzv_my_ne = 10
!mg 11/16 cleanup  orbrzv_my_nx = 12
!mg 11/16 cleanup  allocate(orbrzv_my_rfdist(1:orbrzv_my_nr,&
!mg 11/16 cleanup                            1:orbrzv_my_ne,&
!mg 11/16 cleanup                            1:orbrzv_my_nx))
!mg 11/16 cleanup  allocate(orbrzv_my_rfdist_tmp(1:orbrzv_my_nr,&
!mg 11/16 cleanup                                1:orbrzv_my_ne,&
!mg 11/16 cleanup                                1:orbrzv_my_nx))
!mg 11/16 cleanup  orbrzv_my_rfdist = 0.0_fp
!mg 11/16 cleanup  orbrzv_my_rfdist_tmp = epslon


  ! rfdbg_p* ----------------------------------------------------------
  allocate(rfdbg_pdiag(1:9,1:mj))
  rfdbg_pdiag  = 0.0_fp
  ! rfdbg_p* ----------------------------------------------------------

  ! distribution function
  allocate(orbrzv_fbm_dvl(1:nsbeam),orbrzv_fbm_dvp(1:nsbeam))
  allocate(orbrzv_fbm(1:nsbeam,1:mj,&
                    & -orbrzv_fbm_nvl:orbrzv_fbm_nvl,&
                    & 0:orbrzv_fbm_nvp)); orbrzv_fbm = 0.0_fp

  ! for each species, find min & max harmonic numbers expected to exist

  ! find maximum freq, kl (nphi)
  zkl = orbrzv_rfw(1,1)%kl
  zfreq = orbrzv_rfw(1,1)%freq
  zfreq_min = orbrzv_rfw(1,1)%freq
  do iant = 1, orbrzv_rf_nant
    do inphi = 1, orbrzv_rf_nnphi(iant)
      zkl = dmax1(zkl, orbrzv_rfw(iant, inphi)%kl)
      zfreq = dmax1(zfreq, orbrzv_rfw(iant, inphi)%freq)
      zfreq_min = dmin1(zfreq_min, orbrzv_rfw(iant, inphi)%freq)
    enddo
  enddo

  ! set information on gradient scale length of erf
  ! lrf_spo stores the inverse of the gradient scale of erf (cm^-1)
  ! pick up maximum lrf_inv among the spectrum
  allocate(zlrfi(1:orbrzv_nr, 1:orbrzv_nz))
  zlrfi = 0.0_fp

  ! free previous setup
  call ezspline_free(orbrzv_lrfi_spo, ierr)

  call ezlinear_init(orbrzv_lrfi_spo, orbrzv_nr, orbrzv_nz, ierr)
  do ii = 1, orbrzv_nr
    orbrzv_lrfi_spo%x1(ii) = orbrzv_rmin + orbrzv_dr*real(ii-1,fp)
  enddo
  do jj = 1, orbrzv_nz
    orbrzv_lrfi_spo%x2(jj) = orbrzv_zmin + orbrzv_dz*real(jj-1,fp)
  enddo

  do iant = 1, orbrzv_rf_nant
    do inphi = 1, orbrzv_rf_nnphi(iant)

      if(orbrzv_rfw(iant,inphi)%freq.le.0.0_fp) cycle !mg, no data for this antenna

      do ii = 1, orbrzv_nr
        zr = orbrzv_rmin + (real(ii-1,fp))*orbrzv_dr
        do jj = 1, orbrzv_nz
          zz = orbrzv_zmin + (real(jj-1,fp))*orbrzv_dz

          call ezspline_interp2(orbrzv_psi_spo, zr, zz, zpsi, ierr)
          call ezspline_gradient2(orbrzv_psi_spo, zr, zz, zdpsi, ierr)
          zpsi = abs(zpsi)
          zdrhodpsi = 0.5_fp/(sqrt(zpsi*orbrzv_psix)+epslon)

          zrho = sqrt(zpsi/orbrzv_psix)
          
          zrtmp = zr - orbrzv_rc
          zztmp = zz - orbrzv_zc

          if(zrtmp == 0.0_fp .and. zztmp == 0.0_fp) then
            ztheta = 0.0_fp
          else
            ztheta = atan2(zztmp, zrtmp) 
          endif

          if(ztheta < 0.0_fp) then
            ztheta = ztheta + twopi
          else if(ztheta >= twopi) then 
            ztheta = ztheta - twopi
          endif

          if(zrho < 1.0_fp) then
            call ezspline_interp2(orbrzv_rfw(iant,inphi)%ep_spo, &
                                  zrho, ztheta, zerf, ierr)
            call ezspline_gradient2(orbrzv_rfw(iant,inphi)%ep_spo, &
                                  zrho, ztheta, zderf, ierr)

            zerf = abs(zerf) + epslon
            zrr = sqrt((zr-orbrzv_rc)**2+(zz-orbrzv_zc)**2)
            zdtdr = -sin(ztheta)/(zrr+epslon)
            zdtdz =  cos(ztheta)/(zrr+epslon)
            zdrhodr = zdrhodpsi*zdpsi(1)
            zdrhodz = zdrhodpsi*zdpsi(2)

            zlr = abs(zdrhodr*zderf(1) + zdtdr*zderf(2))/zerf
            zlz = abs(zdrhodz*zderf(1) + zdtdz*zderf(2))/zerf
            zlrf_inv = max(zlr, zlz)

            zlrfi(ii,jj) = max(zlrfi(ii,jj), zlrf_inv)
          endif

        enddo
      enddo

    enddo
  enddo

  call ezspline_setup(orbrzv_lrfi_spo, zlrfi, ierr)
  call ezspline_error(ierr)
  deallocate(zlrfi)


  ! estimate max and min harmonic number for each species
  zfreq = zfreq*twopi ! convert to angular freq
  zfreq_min = zfreq_min*twopi ! convert to angular freq
  zxe2cgs = xe*1.0e7_fp ! energy conversion factor ev --> erg
  do ii = 1, nsbeam
    zms = zmp*abeams(ii)
    zqs = zel*xzbeams(ii)

    ! maximum velocity of the species
    zvmax = dsqrt(2.0_fp*zxe2cgs*18.0*ebdmaxa(ii)/zms) !added factor 18, see ebatom if rf is on,  mg 04/07/16 
    zkvmax = zkl*zvmax/orbrzv_prmin ! max n*vl/r
!mg    zkvmax = dabs(zkvmax)
    zkvmin = zkl*vmins(ii)/orbrzv_prmax ! min n*vl/r, mg 04/07/16 

    ! fbm-distribution function dimension
    orbrzv_fbm_dvl(ii) = zvmax*5.0_fp/real(orbrzv_fbm_nvl,fp)
    orbrzv_fbm_dvp(ii) = zvmax*5.0_fp/real(orbrzv_fbm_nvp,fp)

    ! rf related quantities
    ! max and min of gyro-frequency within the plasma boundary
    zomg_max = zqs*bmaxm(ledge)*1.0e4_fp/(zms*zc) ! cgs
    zomg_min = zqs*bminm(ledge)*1.0e4_fp/(zms*zc) ! cgs

    ! maximum and minimum harmonic number present within plasma boundary
    orbrzv_tmp(1)=(zfreq - zkvmax)/zomg_max
    orbrzv_tmp(2)=(zfreq - zkvmax)/zomg_min
    orbrzv_tmp(3)=(zfreq + zkvmax)/zomg_max
    orbrzv_tmp(4)=(zfreq + zkvmax)/zomg_min
    orbrzv_tmp(5)=(zfreq_min - zkvmax)/zomg_max
    orbrzv_tmp(6)=(zfreq_min - zkvmax)/zomg_min
    orbrzv_tmp(7)=(zfreq_min + zkvmax)/zomg_max
    orbrzv_tmp(8)=(zfreq_min + zkvmax)/zomg_min
    orbrzv_tmp(9)=(zfreq - zkvmin)/zomg_max
    orbrzv_tmp(10)=(zfreq - zkvmin)/zomg_min
    orbrzv_tmp(11)=(zfreq + zkvmin)/zomg_max
    orbrzv_tmp(12)=(zfreq + zkvmin)/zomg_min
    orbrzv_tmp(13)=(zfreq_min - zkvmin)/zomg_max
    orbrzv_tmp(14)=(zfreq_min - zkvmin)/zomg_min
    orbrzv_tmp(15)=(zfreq_min + zkvmin)/zomg_max
    orbrzv_tmp(16)=(zfreq_min + zkvmin)/zomg_min
!mg 04/07/16
    orbrzv_lmin(ii) =dint(minval(orbrzv_tmp))
    orbrzv_lmin(ii) =max(-lup_max,orbrzv_lmin(ii))-1
    orbrzv_lmax(ii) = dint(maxval(orbrzv_tmp))
    orbrzv_lmax(ii) = min(lup_max,orbrzv_lmax(ii))+1
!mg 04/07/16   orbrzv_lmin(ii) = max(orbrzv_lmin(ii), 1) !by the definition- resonance number >=1
!mg 04/07/16   orbrzv_lmax(ii) = max(orbrzv_lmax(ii), 1)

    !print *, "zfreq, zkvmax", zfreq, zkvmax
    !print *, "zomg_max, zomg_min", zomg_min, zomg_max

    ! maximum resonance time allowed
    ! kjm-constrf
    !zbdbq_m = dabs(1.0e7_fp*orbrzv_bdb_min*zqs/(zms*zc))
    zbdbq_m = dabs(1.0e6_fp*orbrzv_bdb_min*zqs/(zms*zc))
    orbrzv_trmax(ii) = dsqrt(twopi/zbdbq_m)
  enddo
  ll=minval(orbrzv_lmin)
  lup=maxval(orbrzv_lmax)
  
  call alloc_orbrzv_rfw(1,ll,lup)
  if(allocated(orbrzv_zbesl)) deallocate(orbrzv_zbesl)
  orbrzv_llmax=max(abs(ll),abs(lup))+2
  allocate(orbrzv_zbesl(1:orbrzv_llmax))
  orbrzv_zbesl=0.0_fp

  ! rf-field strength renormalizer
  !orbrzv_toric_frnm = 3.3_fp
  !orbrzv_toric_frnm = 4.7_fp
  ! initial time setting from input
  !if(orbrzv_toric_frnm < 0.0_fp) then
  !  orbrzv_toric_frnm = orbrzv_toric_frnm0
  !endif

  ! controler to set onoff rf-realted goosing
  nlmcrfgcon = .false.

  !orbrzv_out = 314159
  !open(unit=orbrzv_out, file="rfinfo.txt")
!  call infolog(' orbrz_set_rf2: toric_frnm =',orbrzv_toric_frnm)
  do ii = 1, nsbeam
     call infolog(' orbrz_set_rf2: trmax(ii)=',orbrzv_trmax(ii))
     call infolog('                species ii=',ii)
     call infolog('                species ebdmaxa',ebdmaxa(ii))
     call infolog('                species abeama',abeams(ii))
     call infolog('                species xzbeama',xzbeams(ii))
     call infolog('                l-range : 1, <= ',orbrzv_lmin(ii))
     call infolog('                ,',orbrzv_lmax(ii))
     call infolog('                <=',orbrzv_llmax)
!mg    orbrzv_lmin(ii) = max(orbrzv_lmin(ii), 1)
!mg    orbrzv_lmax(ii) = min(orbrzv_lmax(ii), orbrzv_llmax)
  enddo
  do ii = 1, ng
     call infolog(' orbrz_set_rf2: thermal ii=',ii)
     call infolog('                thermal aplasm=',aplasm(ii))
     call infolog('                thermal backz=',backz(ii))
  enddo

  ! print out rfw interpolation for debugging ----------------------------- !
!#ifdef __DEBUG
#ifndef __MPI
  orbrzv_out = 314159
  open(unit=orbrzv_out, file="rfw.txt")

  do ii = 1, nsbeam
    write (unit=orbrzv_out, fmt=*) "# ion", ii, ":", &
          & ebdmaxa(ii),abeams(ii),xzbeams(ii)
  enddo

  write (unit=orbrzv_out, fmt=*) "# sim box r : ", orbrzv_rmin, orbrzv_rmax
  write (unit=orbrzv_out, fmt=*) "# sim box z : ", orbrzv_zmin, orbrzv_zmax

  do ii = 1, orbrzv_nr
    zr = orbrzv_rmin + (real(ii-1,fp))*orbrzv_dr
    do jj = 1, orbrzv_nz
      zz = orbrzv_zmin + (real(jj-1,fp))*orbrzv_dz

      call ezspline_interp2(orbrzv_psi_spo, zr, zz, zpsi, ierr)
      zpsi = abs(zpsi)
      zrho = sqrt(zpsi/orbrzv_psix)

      zrtmp = zr - orbrzv_rc
      zztmp = zz - orbrzv_zc

      if(zrtmp == 0.0_fp .and. zztmp == 0.0_fp) then
        ztheta = 0.0_fp
      else
        ztheta = atan2(zztmp, zrtmp) 
      endif

      if(ztheta < 0.0_fp) then
        ztheta = ztheta + twopi
      else if(ztheta >= twopi) then 
        ztheta = ztheta - twopi
      endif

      if(zrho < 1.0_fp) then
        call ezspline_interp2(orbrzv_rfw(1,1)%kp_spo, &
                              zrho, ztheta, zomg_max, ierr)
        call ezspline_interp2(orbrzv_rfw(1,1)%ep_spo, &
                              zrho, ztheta, zomg_min, ierr)
        call ezspline_interp2(orbrzv_lrfi_spo, zr, zz, zlrf_inv, ierr)
      else
        zomg_max = 0.0_fp
        zomg_min = 0.0_fp
      endif

      write (unit=orbrzv_out, fmt=*) zr, zz, zomg_max, zomg_min, &
                                     zrho, ztheta, zlrf_inv
    enddo
    write (unit=orbrzv_out, fmt=*) " "
  enddo
  close(unit=orbrzv_out)
  call infolog(' orbrz_set_rf_info2:orbrz_rf_out completed')
#endif
!#endif
  !read *, zkl
  ! ----------------------------------------------------------------------- !

  !call rfxqlo_free(iwarn)

  contains

    integer function mcrf_idxf(ii, jj)
      integer, intent(in) :: ii, jj
      mcrf_idxf = jj + idx_nth*(ii-1)
    end function mcrf_idxf

    integer function mcrf_idxp(j, isp)
      integer, intent(in) :: j, isp
      mcrf_idxp = isp + nsbeam*(j-1)
    end function mcrf_idxp

    subroutine alloc_orbrzv_rfw(istat,ll,lup)
      implicit none
      integer, intent(in) :: istat,ll,lup
      integer :: iant, inphi
      if(istat.eq.0) then
         allocate(orbrzv_rfw(1:orbrzv_rf_nant, 1:orbrzv_rf_max_nphi))
      elseif(istat.eq.1) then
         do iant = 1, orbrzv_rf_nant
            do inphi = 1, orbrzv_rf_max_nphi
               allocate(orbrzv_rfw(iant,inphi)%lam0(ll:lup,2),orbrzv_rfw(iant,inphi)%lam(ll:lup,2))
               allocate(orbrzv_rfw(iant,inphi)%dsgn(ll:lup,2),orbrzv_rfw(iant,inphi)%time(ll:lup,2))
               allocate(orbrzv_rfw(iant,inphi)%dt(ll:lup,2),orbrzv_rfw(iant,inphi)%nstep(ll:lup,2))
               allocate(orbrzv_rfw(iant,inphi)%time0(ll:lup,2),orbrzv_rfw(iant,inphi)%dt0(ll:lup,2))
               allocate(orbrzv_rfw(iant,inphi)%ztf(ll:lup,2),orbrzv_rfw(iant,inphi)%zzm(ll:lup,2,4))
               allocate(orbrzv_rfw(iant,inphi)%kp2(ll:lup,2),orbrzv_rfw(iant,inphi)%ep(ll:lup,2))
               allocate(orbrzv_rfw(iant,inphi)%cpol(ll:lup,2))
               orbrzv_rfw(iant,inphi)%kp2(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%ep(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%cpol(ll:lup,1:2)=(0.0_fp,0.0_fp)
               orbrzv_rfw(iant,inphi)%ztf(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%zzm(ll:lup,1:2,1:4)=0.0_fp
               orbrzv_rfw(iant,inphi)%time0(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%dt0(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%nstep(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%dt(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%dsgn(ll:lup,1:2)=0
               orbrzv_rfw(iant,inphi)%time(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%lam0(ll:lup,1:2)=0.0_fp
               orbrzv_rfw(iant,inphi)%lam(ll:lup,1:2)=0.0_fp
            enddo
         enddo
      endif
    end subroutine alloc_orbrzv_rfw

    subroutine dealloc_orbrzv_rfw(ierr)
      implicit none
      integer, intent(out) :: ierr
      integer :: iant, inphi

      do iant = 1, orbrzv_rf_nant
         do inphi = 1, orbrzv_rf_nnphi(iant)
            call ezspline_free(orbrzv_rfw(iant,inphi)%kp_spo, ierr)
            call ezspline_free(orbrzv_rfw(iant,inphi)%ep_spo, ierr)
            call ezspline_free(orbrzv_rfw(iant,inphi)%rp_spo, ierr)
            call ezspline_free(orbrzv_rfw(iant,inphi)%ip_spo, ierr)
            deallocate(orbrzv_rfw(iant,inphi)%lam0,orbrzv_rfw(iant,inphi)%lam)
            deallocate(orbrzv_rfw(iant,inphi)%dsgn,orbrzv_rfw(iant,inphi)%time)
            deallocate(orbrzv_rfw(iant,inphi)%dt,orbrzv_rfw(iant,inphi)%nstep)
            deallocate(orbrzv_rfw(iant,inphi)%time0,orbrzv_rfw(iant,inphi)%dt0)
            deallocate(orbrzv_rfw(iant,inphi)%ztf,orbrzv_rfw(iant,inphi)%zzm)
            deallocate(orbrzv_rfw(iant,inphi)%kp2,orbrzv_rfw(iant,inphi)%ep)
            deallocate(orbrzv_rfw(iant,inphi)%cpol)
         enddo
      enddo

      deallocate(orbrzv_rfw)
    end subroutine dealloc_orbrzv_rfw

end subroutine orbrz_set_rf_info2
