subroutine mcrf_toric_pack2(ierr)

  use iso_c_binding, only: c_double_complex
  use rfxqlo_interface
  use orbrzv
  use nbi_com, only : nsbeam,xzbeams,abeams,nlfprod,nonlin,nzones,mj,orbrzv_zzerr_con
  use nbi_mpi_com
  use logmod

  implicit none

  integer, intent(out) :: ierr  ! status code on exit (0=ok)
  integer :: ierr0 

  !-----------------------------
  !  fetch data from rfxqlo; pack into buffers for mpi broadcast
  !-----------------------------

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

  integer ::  j, isp,is_p, iant, inphi
  integer :: inx, inth, inions, inant, imaxnphi, inanp
  integer :: nsize(1:2), iwfdim, inbb, ioffset
  integer :: iznphi, ksym, imode
  integer :: idx_nx, idx_nth
  integer, allocatable, dimension(:) :: inum_nphi
  integer :: innap
  integer :: ilunt, iopstat

  real(fp) :: zfreq, zpower
  real(fp), allocatable, dimension(:,:) :: zkp, zep
  complex(c_double_complex), allocatable, dimension(:,:) :: zcrfpol
  real(fp), allocatable, dimension(:,:) :: zpwr
  integer :: ii,jj,iz,ia,iz0,ia0,imatch,iname,icount
  character*120 :: icname
  integer, dimension(:), allocatable :: imap
  logical :: ilname
  integer :: isave_freq

  ierr=0
  call nbi_mpi_com_init
  if(myidd.gt.0) return ! only myid=0 go there, to fill in mcrf_...buffer
  !-----------------------------------------------------
  ! get wavefield dimension information
  call rfxqlo_wavefield_dims(inx, inth, ierr)
  if(ierr.ne.0) return
  allocate(zkp(inx,inth), zep(inx,inth))
  zkp=0.0_fp
  zep=0.0_fp
  allocate(zcrfpol(inx,inth))
  zcrfpol=0.0_fp

  !-----------------------------------------------------
  ! get the number of ion species
  call rfxqlo_ion_get(zions, ierr)
  if(ierr.ne.0) return
  inions = zions%n_ions
  if(allocated(imap)) deallocate(imap)
  allocate(imap(inions+2))
  imap=0
  icount=0
  do ii=1,inions
     iz=zions%z_ions(ii) + 0.1d0
     ia=zions%a_ions(ii) + 0.1d0
     iname=len_trim(zions%name_ions(ii))
     icname=zions%name_ions(ii)
     if(iname.lt.4) cycle ! as all species have 'rfmi', 'beam', 'fusn' or none if thermal
     if(icname(iname-3:iname).eq.'rfmi'.or.icname(iname-3:iname).eq.'    ') cycle
     ! only beam and fusion needed
     ilname=(icname(iname-3:iname).eq.'fusn')
     imatch = 0
     do jj=1,nsbeam
        iz0=xzbeams(jj)+0.1d0
        ia0=abeams(jj)+0.1d0
        if((ia0.eq.ia).and.(iz0.eq.iz).and.(ilname.eqv.nlfprod(jj))) then !eqv -to satisfy gcc4
           imatch=jj
           icount=icount+1
           exit
        endif
     enddo
     if(imatch.gt.0) then
        imap(ii)=imatch
     endif
  enddo
  if(icount.gt.nsbeam)then
     call errorlog(' mcrf_toric_pack2: inconsistency in # of fast species found. nsbeam=',nsbeam)
     call errorlog(' mcrf_toric_pack2: found icount= ',icount)
     write(nonlin,*)' ? mcrf_toric_pack2:inconsistency in # of fast species found. nsbeam=',&
          &nsbeam,' found icount=',icount
     ierr=1
     return
  endif
  !-----------------------------------------------------
  ! get antenna and nphi spectrum number
  call rfxqlo_dim_get(zdims, ierr)
  if(ierr.ne.0) return
  inant    = zdims%n_antenna
  imaxnphi = zdims%max_nphi
  allocate(inum_nphi(1:inant))
  inum_nphi=0

  !-----------------------------------------------------
  ! get toroidal spectrum for each antenna
  innap = 0
  do iant = 1, inant
    call rfxqlo_select_antenna(iant, ierr)
    if(ierr.ne.0) return
    call rfxqlo_antenna_get(zant, ierr)
    if(ierr.ne.0) return
    inum_nphi(iant) = zant%num_nphi
    innap = innap + inum_nphi(iant)
  enddo

  !-----------------------------------------------------
  ! allocate array to hold profiles (heating, current by rf)
  if(allocated(mcrf_prof)) deallocate(mcrf_prof)
  allocate(mcrf_prof(1:inant, 1:imaxnphi))
  
  do iant = 1, inant
      do inphi=1,imaxnphi
         call rfxqlo_outpro_blank(mcrf_prof(iant,inphi))
      enddo
   enddo
  allocate(zpwr(1:nzones+1, inions+2))
  zpwr = 0.0_fp

  !-----------------------------------------------------
  ! setup for index function
  idx_nx  = inx
  idx_nth = inth

  !-----------------------------------------------------
  ! size of mcrf_ibuf
  nsize(1) = 5+inant

  !-----------------------------------------------------
  iwfdim = inx*inth
  ! real array size for each wave spectrum
  inbb = 4 + 4*iwfdim + mj*nsbeam

  !-----------------------------------------------------
  ! size of mcrf_rbuf
  nsize(2) = innap*inbb !nsize(2) = imaxnphi*inant*inbb

  !-----------------------------------------------------
  !  mcrf_nbuf is statically allocated, size=3

  ! size of mcrf_ibuf
  mcrf_nbuf(1) = nsize(1)
  ! size of mcrf_rbuf
  mcrf_nbuf(2) = nsize(2)

  mcrf_nbuf(3) = 1  ! in old code this meant rf data is present

  !-----------------------------------------------------
  ! free and allocate buffer to hold rf-information
  if(allocated(mcrf_ibuf)) deallocate(mcrf_ibuf)
  if(allocated(mcrf_rbuf)) deallocate(mcrf_rbuf)
  allocate(mcrf_ibuf(mcrf_nbuf(1)))
  allocate(mcrf_rbuf(mcrf_nbuf(2)))
  mcrf_rbuf = 0.0_fp
  mcrf_ibuf = 0
  !-----------------------------------------------------
  ! fill up integer buffer
  mcrf_ibuf(1) = inant
  mcrf_ibuf(2) = imaxnphi
  mcrf_ibuf(3) = inx
  mcrf_ibuf(4) = inth
  mcrf_ibuf(5) = inions
  do iant = 1, inant
    mcrf_ibuf(5+iant) = inum_nphi(iant)
  enddo

  !-----------------------------------------------------
  ! fill up real buffer
  ioffset = 0
  ! information of each nphi spectrum
  do iant = 1, inant
    call rfxqlo_select_antenna(iant, ierr)
    if(ierr.ne.0) return

    do inphi = 1, inum_nphi(iant)
      call rfxqlo_select_nphi(inphi, ierr)
      if(ierr.ne.0) return

      ! fetch basic information for this spectrum
      call rfxqlo_rfbasic(zfreq, zpower, iznphi, ksym, ierr)
      if(ierr.ne.0) return

      ! fetch wavefield information
      imode = 0
      call rfxqlo_wavefield_get1(imode, zkp, zep, zcrfpol, ierr)
      if(ierr.ne.0) return


      mcrf_rbuf(ioffset+1) = zfreq
      isave_freq = ioffset+1 !save an address for a frequency storage
      mcrf_rbuf(ioffset+2) = zpower/2._fp
      mcrf_rbuf(ioffset+3) = real(iznphi,fp)
      mcrf_rbuf(ioffset+4) = real(ksym,fp)
      ioffset = ioffset + 4

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

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

      do ii = 1, inx
        do jj = 1, inth
          mcrf_rbuf(ioffset+mcrf_idxf(ii, jj)) = real(zcrfpol(ii,jj),fp)
        enddo
      enddo
      ioffset = ioffset + iwfdim

      do ii = 1, inx
        do jj = 1, inth
          mcrf_rbuf(ioffset+mcrf_idxf(ii, jj)) = dimag(zcrfpol(ii, jj))
        enddo
      enddo
      ioffset = ioffset + iwfdim

      ! fetch profiles : power, current, etc.
      ierr0=0
      call rfxqlo_outpro_fetch(mcrf_prof(iant,inphi), ierr0)
      if(ierr0.ne.0) then
         mcrf_rbuf(isave_freq)=0.0_fp ! set 0 frequency if there is no profile saved
         ioffset = ioffset + mj*nsbeam
         cycle
      endif
      do isp = 1, inions+2
         is_p=imap(isp)
         if(is_p.eq.0) cycle
         if(mcrf_prof(iant,inphi)%id_pwrs(isp).eq.0) then
            call errorlog(' mcrf_toric_pack2: id_pwrs = 0, antenna=',iant)
            cycle
         endif
        call mcrf_fetch_pwr(mcrf_prof(iant,inphi)%id_pwrs(isp), &
                          & zpwr(1:nzones+1,isp), ierr)
        if(ierr.ne.0) then
           call errorlog(' mcrf_toric_pack2: error getting pwr prof from xplasma, antenna=',iant)
           goto 999
        endif

        do j = 1, nzones+1
          mcrf_rbuf(ioffset+mcrf_idxp(j,is_p)) = zpwr(j, isp)/2._fp
        enddo
        !do j = lcentr, ledge
        !  jj = j - lcentr + 1 
        !  mcrf_rbuf(ioffset+mcrf_idxp(j,isp))=zpwr(jj+1,isp)-zpwr(jj,isp)
        !enddo

      enddo
      ioffset = ioffset + mj*nsbeam

    enddo
  enddo

  !-----------------------------------------------------
  if(ioffset .ne. nsize(2)) then
    print *, ' ?mcrf_toric_pack : size not match - ', ioffset, nsize(2)
    ierr=999
  endif

999 continue
  if(allocated(mcrf_prof)) deallocate(mcrf_prof)
  if(allocated(inum_nphi)) deallocate(inum_nphi)
  if(allocated(zkp)) deallocate(zkp, zep)
  if(allocated(zcrfpol)) deallocate(zcrfpol)
  if(allocated(zpwr)) deallocate(zpwr)
  if(allocated(imap)) deallocate(imap)

  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

end subroutine mcrf_toric_pack2

!==========================================================================
! fetch power profile indicated by id_pwr

subroutine mcrf_fetch_pwr(id_pwr, zpwr, ierr)

  use nbi_com

  implicit none

  integer, intent(in) :: id_pwr ! xplasma id for pwr to fetch

  real(fp), intent(out) :: zpwr(1:nzones+1)
  integer, intent(out) :: ierr  ! status code on exit (0=ok)

  integer :: j, jj
  real(fp) :: zxis(1:nzones+1), zxsmoo, zpwr_tmp(1:nzones+1)

  zpwr = 0.0_fp
  zpwr_tmp = 0.0_fp

  ! note that mj = nzones+3
  zxis = xizb(lcentr:lep1)
  zxis(1) = 0.0_fp
  zxis(nzones+1) = 1.0_fp

  ! smoothing factor
  zxsmoo = 0.0_fp

  call eq_rget_binsm(nzones+1, zxis, id_pwr, zxsmoo, 2, &
                   & zpwr_tmp(1:nzones+1), ierr)
        
  do j = lcentr, ledge
    jj = j - lcentr + 1 
    zpwr(j) = zpwr_tmp(jj+1) - zpwr_tmp(jj)
  enddo

end subroutine mcrf_fetch_pwr
