subroutine rfxqlo_iconc(inum,aconc,mainsp,zzeff,ierr)

  !  rfxqlo routine to satisfy legacy requirement of TORIC...
  !  get Zeff, "main" ion specie index, and relative concentrations
  !  n_i/n_e.

  use rfxqlo_private
  implicit NONE

  integer, intent(in) :: inum     ! dimension of aconc(...) must be large
  !                               ! enough to accomdate all ion species

  real(r8qlo), intent(out) :: aconc(inum)  ! n_i/n_e for each ion specie

  integer, intent(out) :: mainsp  ! index to specie with maximum n_i/n_e

  real(r8qlo), intent(out) :: zzeff   ! volume averaged Zeff
  !                                   ! sum (n_i*Z_i^2/n_e)

  integer, intent(out) :: ierr    ! completion code: 0=OK

  !-----------------------------------------

  real(r8qlo), dimension(:), allocatable :: zne
  real(r8qlo), dimension(:,:), allocatable :: zni
  real(r8qlo), dimension(:), allocatable :: zconc,zfracmax
  integer, dimension(:), allocatable :: imainsp
  logical, dimension(:), allocatable :: iscreen
  integer i,j,k,id,inumi,mainsp1
  real(r8qlo) :: znetot,znesub,zzi,zepsn,zfrac,zzmain
  integer :: iminsp(inum)

  !-----------------------------------------

  aconc=0
  mainsp=0
  mainsp1=0
  zzeff=0

  zepsn = rfx_options%epsn

  !-----------------------------------------

  inumi = rfx_ions%n_ions        ! no. of ion species

  allocate(zne(nx_qlo))
  allocate(zni(nx_qlo,inumi))
  allocate(zconc(nx_qlo),zfracmax(inumi))
  allocate(imainsp(nx_qlo))
  allocate(iscreen(inumi))

  do
     id=rfx_prof%id_ne
     call eq_rgetf(nx_qlo,x_qlo,id,0,zne,ierr)
     if(ierr.ne.0) exit

     do j=1,inumi
        id=rfx_ions%id_ni(j)
        call eq_rgetf(nx_qlo,x_qlo,id,0,zni(:,j),ierr)
        if(ierr.ne.0) exit
     enddo
     if(ierr.ne.0) exit

     zfracmax=0.0d0
     do j=1,inumi
        do i=1,nx_qlo
           zfrac = zni(i,j)/zne(i)
           zfracmax(j)=max(zfracmax(j),zfrac)
        enddo
     enddo

     iscreen=.FALSE.
     do j=1,inumi
        if(zfracmax(j).lt.zepsn) iscreen(j)=.TRUE.
     enddo

     znetot=0
     znesub=0
     iminsp=1
     do i=1,nx_qlo
        znetot=znetot+zne(i)*dv_qlo(i)
        do j=1,inumi
           zzi=rfx_ions%z_ions(j)
           if(iscreen(j)) then
              znesub = znesub + zzi*zni(i,j)*dv_qlo(i)
              zni(i,j)=0.0d0
           else
              aconc(j)=aconc(j)+zni(i,j)*dv_qlo(i)
              zzeff=zzeff+zni(i,j)*zzi*zzi*dv_qlo(i)
           endif
           if (zni(i,j)<zni(iminsp(j),j)) iminsp(j)=i
        enddo
     enddo

     mainsp=1
     do j=2,inumi
        if(aconc(j).gt.aconc(mainsp)) mainsp=j
     enddo

     do i=1,nx_qlo
        imainsp(i)=1
        zconc(i)=zni(i,imainsp(i))
        do j=2,inumi
           if(zni(i,j).gt.zni(i,imainsp(i))) then
              imainsp(i)=j
              zconc(i)=zni(i,j)
           endif
        enddo
     enddo

     mainsp1=1 
     do k=1,nx_qlo
        if(zconc(k).lt.zconc(mainsp1)) mainsp1=k
     enddo

     mainsp1=imainsp(mainsp1)
     if (mainsp1.ne.mainsp) then
        ! select the species with the largest minimum density
        if (zni(iminsp(mainsp1),mainsp1)>zni(iminsp(mainsp),mainsp)) then
           mainsp=mainsp1
        end if
     end if

     zzmain=rfx_ions%z_ions(mainsp)

     aconc(mainsp)=aconc(mainsp)+znesub/zzmain
     zzeff=zzeff+znesub*zzmain

     zzeff=zzeff/znetot
     aconc=aconc/znetot

     exit
  enddo
  deallocate(zne,zni,zconc,zfracmax,imainsp,iscreen)

end subroutine rfxqlo_iconc
