subroutine nubeam_load_state(ss,beam_present,fusi_present,ierr)

  ! load plasma state profiles output by nubeam
  ! use the traditional xplasma id output blocks

  ! it is assumed plasma state has been prepared for use by the nubeam
  ! component, and nubeam internal species lists have been set up to match
  ! the plasma state lists (i.e. this is part of the nubeam component
  ! interface).  nubeam uses a unified list for both beams and fusion
  ! products, ordered s.t. beam ions all come before fusion ions.

  ! note (sep 2008) frantic or other neutral gas model needed to complete
  ! the beam driven sources and compute ss%sb0halo_recap(...) & ss%sbsce(...)

  ! future version of nubeam will handle this internally
  !-------------------------------------------

  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: zero,one
  use nbi_types
  use plasma_state_mod

  implicit none

  !-------------------------------------------
  ! arguments:

  type (plasma_state) :: ss
  logical :: beam_present,fusi_present

  integer, intent(out) :: ierr

  !-------------------------------------------
  ! local:

  type (nbotype_density) :: zdens
  type (nbotype_average_energies) :: zengy
  type (nbotype_powers) :: zpower
  type (nbotype_currents) :: zcur
  type (nbotype_sources) :: zsce
  type (nbotype_neutral_sources) :: zsc0
  type (nbotype_neutral_sinks) :: zsn0
  type (nbotype_torques) :: ztq
  type (nbotype_hi_z_beams) :: ztzb

  integer :: ii,jj,ig,is,isb,inrho
  integer :: iaf,izf,iath,izth,izb
  real(fp), dimension(:), allocatable :: zx,zxzc,zwk,zsuma,zterma,zvol,zdvol
  real(fp), dimension(:), allocatable :: zsumx

  real(fp) :: zxrz(2),zans(1)  ! rezone-to-scalar
  integer, dimension(:), allocatable :: iqav

  integer, parameter :: idens = 1   ! integrated density profile
  integer, parameter :: iparam = 2  ! non-integrated parameter profile
  !-------------------------------------------

  ierr = 0

  inrho = 0

  if(beam_present) inrho = ss%nrho_nbi
  if(fusi_present) inrho = ss%nrho_fus  ! if both are present these will match

  if(inrho.eq.0) return

  zxrz(1)=zero
  zxrz(2)=one

  if(allocated(iqav)) deallocate(iqav)

  !------------------
  ! working arrays...

  allocate(zx(inrho),zwk(inrho))  ! bdy oriented grid & evaluation buffer
  allocate(zxzc(inrho-1),zsuma(inrho-1),zterma(inrho-1),zsumx(inrho-1))
  allocate(zvol(inrho),zdvol(inrho-1))

  zx = ss%rho_nbi
  zxzc = (zx(1:inrho-1)+zx(2:inrho))*0.5_fp

  call ps_intrp_1d_vec(zx,ss%id_vol,zvol,ierr,state=ss)
  if(ierr.ne.0) then
     call nbi_warning(' ??nubeam_load_state: fetch of volume elements failed.')
     return
  endif

  zdvol = zvol(2:inrho) - zvol(1:inrho-1)

  !-----------------
  ! fetch profile ids

  call nbo_get_density(zdens)
  call nbo_get_average_energies(zengy)
  call nbo_get_powers(zpower)
  call nbo_get_currents(zcur)
  call nbo_get_sources(zsce)
  call nbo_get_neutral_sources(zsc0)
  call nbo_get_neutral_sinks(zsn0)
  call nbo_get_torques(ztq)
  call nbo_get_hi_z_beams(ztzb)

  ! get the data...

  do
     !-----------------
     ! densities and average energies

     isb=0
     if(beam_present) then
        do ii=1,ss%nspec_beam
           isb=isb+1
           izb=ss%qatom_snbi(ii)/ps_xe + 0.1_fp

           call fetch1(zdens%io_bdenss(isb),idens,zterma,'bdens')
           if(ierr.ne.0) exit
           do jj=1,inrho-1
              if(zterma(jj).lt.zero) then
                 ss%nbeami(jj,ii) = zero
              else
                 ss%nbeami(jj,ii) = zterma(jj)/zdvol(jj)
              endif
           enddo

           call fetch1(zengy%io_ebapl2s(isb),iparam,ss%epll_beami(:,ii), &
                '<epll>beam')
           if(ierr.ne.0) exit

           call fetch1(zengy%io_ebapp2s(isb),iparam,ss%eperp_beami(:,ii), &
                '<eperp>beam')
           if(ierr.ne.0) exit

           if(izb.gt.2) then
              call fetch1(ztzb%io_xzbmav(isb),iparam,ss%qbeami(:,ii), &
                   '<q>impurity_beam')
              ss%qbeami(:,ii) = ss%qbeami(:,ii)*ps_xe  ! units conv. to c
              if(.not.allocated(iqav)) then
                 allocate(iqav(ss%nspec_beam)); iqav=0
              endif
              iqav(ii)=1
           else
              ss%qbeami(:,ii) = ps_xe*izb
              ss%q_snbi(ii) = ps_xe*izb
           endif

        enddo
        if(ierr.ne.0) exit
     endif

     if(fusi_present) then
        do ii=1,ss%nspec_fusion
           isb=isb+1

           call fetch1(zdens%io_bdenss(isb),idens,zterma,'fdens')
           if(ierr.ne.0) exit
           do jj=1,inrho-1
              if(zterma(jj).lt.zero) then
                 ss%nfusi(jj,ii) = zero
              else
                 ss%nfusi(jj,ii) = zterma(jj)/zdvol(jj)
              endif
           enddo

           call fetch1(zengy%io_ebapl2s(isb),iparam,ss%epll_fusi(:,ii), &
                '<epll>fusion')
           if(ierr.ne.0) exit

           call fetch1(zengy%io_ebapp2s(isb),iparam,ss%eperp_fusi(:,ii), &
                '<eperp>fusion')
           if(ierr.ne.0) exit

        enddo
        if(ierr.ne.0) exit
     endif

     !-----------------
     ! powers

     if(beam_present) then
        call fetch1(zpower%io_pbe,idens,ss%pbe,'pbe')
        if(ierr.ne.0) exit
        call fetch1(zpower%io_pbi,idens,ss%pbi,'pbi')
        if(ierr.ne.0) exit
     endif

     if(fusi_present) then
        call fetch1(zpower%io_pfe,idens,ss%pfuse,'pfe')
        if(ierr.ne.0) exit
        call fetch1(zpower%io_pfi,idens,ss%pfusi,'pfi')
        if(ierr.ne.0) exit
     endif

     isb=0

     if(beam_present) then
        zsuma=zero
        do ii=1,ss%nspec_beam
           isb=isb+1

           ! thermalization power
           call fetch1(zpower%io_pbths(isb),idens,zterma,'pbth')
           if(ierr.ne.0) exit
           zsuma = zsuma + zterma

           ! source friction
           call fetch1(zpower%io_pbthr2s(isb),idens,zterma,'pbthr2')
           if(ierr.ne.0) exit
           zsuma = zsuma + zterma
        enddo
        ss%pbth = zsuma
     endif

     if(fusi_present) then
        zsuma=zero
        do ii=1,ss%nspec_fusion
           isb=isb+1

           ! thermalization energy
           call fetch1(zpower%io_pbths(isb),idens,zterma,'pfth')
           if(ierr.ne.0) exit
           zsuma = zsuma + zterma

           ! source friction
           call fetch1(zpower%io_pbthr2s(isb),idens,zterma,'pfthr2')
           if(ierr.ne.0) exit
           zsuma = zsuma + zterma
        enddo
        ss%pfusth = zsuma
     endif

     !-----------------
     ! torques

     if(beam_present) then
        call fetch1(ztq%io_tqbe,idens,ss%tqbe,'tqbe')
        if(ierr.ne.0) exit
        call fetch1(ztq%io_tqbi,idens,ss%tqbi,'tqbi')
        if(ierr.ne.0) exit
        call fetch1(ztq%io_tqjxb,idens,ss%tqbjxb,'tqjxb')
        if(ierr.ne.0) exit
        call fetch1(ztq%io_tqbth,idens,ss%tqbth,'tqbth')
        if(ierr.ne.0) exit
     endif

     !-----------------
     ! particle sources

     ! note: current version of code needs frantic or other neutral gas model
     ! to complete ss%sbsce -- evaluation of beam halo recapture

     ! cx-related contributions are all treated as from the neutral beams

     isb=0

     ss%sfsce = zero
     ss%sbsce = zero
     ss%sb0halo_recap = zero

     zsuma=zero
     do ii=1,ss%nspec_beam
        isb=isb+1

        ! cold electron source
        call fetch1(zsce%io_sorbes(isb),idens,zterma,'sorbe')
        if(ierr.ne.0) exit
        zsuma = zsuma + zterma
     enddo

     do ii=1,ss%nspec_fusion
        isb=isb+1

        ! cold electron source -- cx recapture, attribute to beams
        call fetch1(zsce%io_sorbes(isb),idens,zterma,'sorbe')
        if(ierr.ne.0) exit
        zsuma = zsuma + zterma
     enddo

     ss%sbedep = zsuma

     ! halo neutral sources

     do ig=1,ss%nspec_gas

        isb = 0
        zsuma=zero

        do ii=1,ss%nspec_beam
           isb=isb+1

           ! thermal neutral source
           call fetch1(zsc0%io_sorbn0s(ig,isb),idens,zterma,'sorbn0')
           if(ierr.ne.0) exit
           zsuma = zsuma + zterma

        enddo

        do ii=1,ss%nspec_fusion
           isb=isb+1

           ! thermal neutral source -- cx, attribute to beams

           call fetch1(zsc0%io_sorbn0s(ig,isb),idens,zterma,'sorbn0')
           if(ierr.ne.0) exit
           zsuma = zsuma + zterma

        enddo

        ss%sb0halo(:,ig) = zsuma
     enddo

     ! thermal neutral sink rates

     do ig=1,ss%nspec_gas

        isb = 0
        zsuma = zero
        zsumx = zero

        do ii=1,ss%nspec_beam
           isb=isb+1

           ! impact ionization of thermal neutrals by beam fast ions

           call fetch1(zsn0%io_sinb0is(ig,isb),iparam,zterma,'sinb0i')
           if(ierr.ne.0) exit
           zsuma = zsuma + zterma

           call fetch1(zsn0%io_sinbn0s(ig,isb),iparam,zterma,'sinbn0')
           if(ierr.ne.0) exit
           zsumx = zsumx + zterma

           ss%rate_sinb0xs(:,ig,ii) = zterma
        enddo

        ss%rate_sinb0i(:,ig) = zsuma
        ss%rate_sinb0x(:,ig) = zsumx

        zsuma = zero
        zsumx = zero

        do ii=1,ss%nspec_fusion
           isb=isb+1

           ! impact ionization of thermal neutrals by fusion product fast ions

           call fetch1(zsn0%io_sinb0is(ig,isb),iparam,zterma,'sinb0i')
           if(ierr.ne.0) exit
           zsuma = zsuma + zterma

           call fetch1(zsn0%io_sinbn0s(ig,isb),iparam,zterma,'sinbn0')
           if(ierr.ne.0) exit
           zsumx = zsumx + zterma

           ss%rate_sinf0xs(:,ig,ii) = zterma
        enddo

        ss%rate_sinf0i(:,ig) = zsuma
        ss%rate_sinf0x(:,ig) = zsumx

     enddo

     isb = 0
     if(beam_present) then
        do ii=1,ss%nspec_beam
           isb=isb+1

           ! thermalization profile
           call fetch1(zsce%io_sorbths(isb),idens,ss%sbtherm(:,ii),'sbth')
           if(ierr.ne.0) exit
        enddo
     endif

     if(fusi_present) then
        do ii=1,ss%nspec_fusion
           isb=isb+1

           ! thermalization profile
           call fetch1(zsce%io_sorbths(isb),idens,ss%sftherm(:,ii),'sfth')
           if(ierr.ne.0) exit

           izf=ss%qatom_sfus(ii)/ps_xe + 0.1_fp
           iaf=ss%m_sfus(ii)/ps_mp + 0.1_fp

           do jj=1,ss%nspec_tha
              izth=ss%qatom_sa(ii)/ps_xe + 0.1_fp
              iath=ss%m_sa(ii)/ps_mp + 0.1_fp

              if((iath.eq.iaf).and.(izth.eq.izf)) then
                 ss%sfsce(:,jj) = ss%sftherm(:,ii)
              endif
           enddo

        enddo
     endif

     !-----------------
     ! driven current

     if(beam_present) then
        call fetch1(zcur%io_curjbeam,idens,ss%curbeam,'curbeam')
        if(ierr.ne.0) exit
     endif

     if(fusi_present) then
        call fetch1(zcur%io_curjfusi,idens,ss%curfusn,'curfusn')
        if(ierr.ne.0) exit
     endif

     exit
  enddo
  deallocate(zx,zxzc,zwk,zsuma,zsumx,zterma,zvol,zdvol)

  call ps_state_memory_update(ierr, state=ss, eqcheck=.false.)

  if(allocated(iqav)) then
     do ii=1,ss%nspec_beam
        if(iqav(ii).eq.1) then
           call ps_rho_rezone(zxrz, ss%id_qbeami(ii), zans, ierr, state=ss)
           ss%q_snbi(ii) = zans(1)
        endif
     enddo
     deallocate(iqav)
  endif

contains

  subroutine fetch1(id,itype,zarr,zitem)

    integer, intent(in) :: id     ! xplasma profile id
    integer, intent(in) :: itype  ! profile type (idens or iparam)

    real(fp), dimension(:) :: zarr  ! output data
    character*(*), intent(in) :: zitem   ! item label, for error message

    !---------------------------
    integer :: ix
    !---------------------------

    if(id.eq.0) then
       zarr = zero  ! xplasma id never set if data is all zero

    else if(itype.eq.idens) then
       call eq_rgetf(inrho,zx,id,0,zwk,ierr)
       do ix=1,inrho-1
          zarr(ix)=zwk(ix+1)-zwk(ix)
       enddo

    else
       call eq_rgetf(inrho-1,zxzc,id,0,zarr,ierr)
    endif

    if(ierr.ne.0) call nbi_warning( &
         ' ?? nubeam_load_state: error detected while fetching: '//trim(zitem))

  end subroutine fetch1

end subroutine nubeam_load_state
