subroutine nbdrive_load(ss,ierr)
  use iso_c_binding, only: fp => c_double

  ! load nbdrive_common, using data from passed Plasma State object,
  ! as well as data from NUBEAM.  It is assumed that NUBEAM has completed
  ! a time step and has output.

  use plasma_state_mod
  use nbdrive_common

  implicit none

  type (plasma_state) :: ss
  integer, intent(out) :: ierr  ! return code, 0=OK

  !--------------------------------
  integer :: ii,jj,iz,ia,itype,indx,idum,inum,iertmp
  integer :: index_nbfi

  real(fp) :: zta,zdta

  integer :: iobj_s(ps_container_size)
  !--------------------------------

  ierr = 0

  !---------------------
  ! RUNID

  runid = ss%runid

  !---------------------
  ! number of thermal species (non-impurities)

  ngmax = ss%nspec_gas
  if(ngmax.eq.0) then
     write(lun_msgs,*) ' ?nbdrive_load: passed state has ss%nspec_gas = 0.'
     ierr=1
     return
  end if

  !---------------------
  ! find the number of fast ion species that are known to NUBEAM

  nsfast = 0

  itype = 1
  do ii=1,ss%nspec_beam
     iz = ss%qatom_snbi(ii)/ps_xe + 0.1_fp
     ia = ss%m_snbi(ii)/ps_mp + 0.1_fp

     indx = index_nbfi(iz,ia,itype)
     if(indx.gt.0) nsfast = nsfast + 1
  end do

  itype = 3
  do ii=1,ss%nspec_fusion
     iz = ss%qatom_sfus(ii)/ps_xe + 0.1_fp
     ia = ss%m_sfus(ii)/ps_mp + 0.1_fp

     indx = index_nbfi(iz,ia,itype)
     if(indx.gt.0) nsfast = nsfast + 1
  end do

  if(nsfast.eq.0) then
     write(lun_msgs,*) ' ?nbdrive_load: passed state and NUBEAM share no fast species.'
     ierr=1
     return
  end if

  nsbfp = nsfast  ! for now assume no RF minorities

  !---------------------
  ! allocate species-indexed arrays in nbdrive_common

  if(.not.allocated(indx_fi)) then
     call alloc_species(ierr)
     if(ierr.ne.0) return
  end if

  !---------------------
  ! set A and Z arrays and labels of thermal species

  do ii=1,ss%nspec_gas
     backz(ii) = ss%qatom_sgas(ii)/ps_xe
     aplasm(ii) = ss%m_sgas(ii)/ps_mp
     call zlabel(backz(ii),aplasm(ii),'thermal plasma ion',label_thions(ii))
  end do

  !---------------------
  ! set A and Z arrays and labels of thermal species
  ! also get the energy grid sizes

  jj = 0  ! index into combined fast species list
  itype = 1
  do ii=1,ss%nspec_beam
     iz = ss%qatom_snbi(ii)/ps_xe + 0.1_fp
     ia = ss%m_snbi(ii)/ps_mp + 0.1_fp

     indx = index_nbfi(iz,ia,itype)
     if(indx.gt.0) then
        jj = jj + 1
        indx_fi(jj)=indx
        xzbeams(jj) = ss%qatom_snbi(ii)/ps_xe
        abeams(jj) = ss%m_snbi(ii)/ps_mp
        call zlabel(xzbeams(jj),abeams(jj),'beam ion',label_fastions(jj))
        nlfprod(jj)=.FALSE.

        ! number of dist. fcn. energy bins, this fast specie:
        call nbi_fbm_nzons_e(indx,nznbmea(jj))
     end if
  end do

  itype = 3
  do ii=1,ss%nspec_fusion
     iz = ss%qatom_sfus(ii)/ps_xe + 0.1_fp
     ia = ss%m_sfus(ii)/ps_mp + 0.1_fp

     indx = index_nbfi(iz,ia,itype)
     if(indx.gt.0) then
        jj = jj + 1
        indx_fi(jj)=indx
        xzbeams(jj) = ss%qatom_sfus(ii)/ps_xe
        abeams(jj) = ss%m_sfus(ii)/ps_mp
        call zlabel(xzbeams(jj),abeams(jj),'fusion product ion',label_fastions(jj))
        nlfprod(jj)=.TRUE.

        ! number of dist. fcn. energy bins, this fast specie:
        call nbi_fbm_nzons_e(indx,nznbmea(jj))
     end if
  end do

  !---------------------
  ! max no. of energy zones for any fast specie:

  nznbme = maxval(nznbmea)

  !---------------------
  ! no. of vpll/v zones -- for now, same for all species

  idum = 0
  call nbi_fbm_nzons_vpv(idum,nznbma)

  !---------------------
  ! allocate energy and vpll/v grids in nbdrive_common

  if(.not.allocated(xxksid)) then
     call alloc_vgrids(ierr)
     if(ierr.ne.0) return
  end if

  !---------------------
  ! load pitch and energy grids

  call nbi_fbm_vpvzons(idum,nznbma,xxksid,iertmp)
  call nbi_fbm_vpvbdys(idum,nznbma+1,xxksidb,iertmp)

  do jj = 1,nsfast
     indx = indx_fi(jj)
     inum=nznbmea(jj)
     call nbi_fbm_ezons(indx,inum,efbm(1:inum,jj),iertmp)
     call nbi_fbm_ebdys(indx,inum+1,efbmb(1:inum+1,jj),iertmp)
  end do

  !---------------------
  ! MCgrid info (2d irregular spatial grid for NUBEAM fbm outputs)
  ! easiest way to get this info is to use old_xplasma calls; NUBEAM
  ! step execution loads xplasma data.

  call mcgrid_find(id_mcgrid)
  if(id_mcgrid.le.0) then

     ierr=2
     write(lun_msgs,*) ' ?nbdrive_load: mcgrid_find returned id_mcgrid = ', &
          id_mcgrid
     return
  end if

  call mcgrid_getsym(id_mcgrid,idum)
  nlsym = (idum.eq.1)  ! .TRUE. means: up-down symmetric geometry

  call mcgrid_getnumr(id_mcgrid,idum,nzones_fb)

  call mcgrid_getnumz(id_mcgrid,idum,nzones_mc)

  if(.not.allocated(nthzns)) then
     call alloc_nthzns(ierr)
     if(ierr.ne.0) return
  end if

  ! use -id_mcgrid to get only the zone counts inside the plasma
  ! +id_mcgrid would also fill in zone counts beyond the plasma, which
  ! might cause a write off the end of an array...

  call mcgrid_getnumzns(-id_mcgrid,nthzns)

  ! add up...

  nthzsum(0)=0
  do jj=1,nzones_fb
     nthzsum(jj) = nthzsum(jj-1) + nthzns(jj)
  end do

  !---------------------------
  ! gfbm_set calls...

  zta = 0.5_fp*(ss%t0 + ss%t1)
  zdta = 0.5_fp*(ss%t1 - ss%t0)

  call gfbm_set_timestep(zta,zdta)

  call gfbm_set_tokamak(ss%tokamak_id)

  call gfbm_set_shot(ss%shot_number)

  call gfbm_set_trgrid(ss%nrho,ss%rho,ss%nmom) ! plasma grid

  call gfbm_set_omega((ss%nrho-1),ss%omegat)   ! plasma toroidal angular veloc.

  !  get old_xplasma "s" pointer as integer :: array
  call itrans_ptr(ps_container_size,iobj_s)

  call gfbm_set_rzmoms(nlsym,ps_container_size,iobj_s) ! MHD equilibrium moments

CONTAINS
  subroutine zlabel(zz,aa,suffix,label)
    use iso_c_binding, only: fp => c_double

    use periodic_table_mod

    !  generate labels for ion species

    real(fp), intent(in) :: zz,aa   ! Z & A of species
    character(len=*), intent(in) :: suffix   ! label suffix
    character(len=*), intent(out) :: label   ! generated label.

    real aaa
    integer :: iz,ia

    character(len=10) :: prefix

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

    aaa = aa
    iz = zz + 0.1
    ia = aa + 0.1

    prefix='?'
    if(iz.eq.1) then
       if(ia.eq.1) then
          prefix='H'
       else if(ia.eq.2) then
          prefix='D'
       else if(ia.eq.3) then
          prefix='T'
       end if
    else if(iz.eq.2) then
       if(ia.eq.3) then
          prefix='He3'
       else if(ia.eq.4) then
          prefix='He4'
       end if
    else

       prefix = to_periodic_table(iz,aaa,-1,0)

    end if

    label = trim(prefix)//' '//trim(suffix)

  end subroutine zlabel

end subroutine nbdrive_load

