subroutine nubeam_ctrl_psupdate(imem_xplasma,ss,pscon,ierr)

  !  set up update block based on plasma state time dependent data.
  !  this is called as part of the step execution sequence for the
  !  plasma state driven nubeam component.

  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: zero,one,half,twopi
  use nbi_types
  use nbi_pstypes
  use plasma_state_mod
  use xplasma_obj_instance

  implicit none

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

  integer, intent(in) :: imem_xplasma ! xplasma memory flag
               ! = 0: old_xplasma pointer "s" available for reuse
               !      (current object pointed to by "s" is replaced with
               !      copy from the plasma state).
               ! > 0: xplasma_obj_instance pointer "s" has been 
               !      prepositioned by caller (copy from plasma state not
               !      required).
               ! < 0: repoint "s" to old_xplasma s2 object, then copy from
               !      plasma state.  later: call this to restore "s" pointer:
               !      subroutine nubeam_ctrl_psupdate_cleanup(ierr)

  type (plasma_state) :: ss
  type (nbi_ps_update) :: pscon

  integer, intent(out) :: ierr
  integer :: ido,ido1,isize,r8size
  integer :: ibbtype


  !---------------------------------------------
  !  field extrapolator (in case needed...)

  external :: eqm_brz_adhoc

  !---------------------------------------------
  !  local:
 
  ! the xplasma ids returned for each profile fit are stored
  ! as elements of "nbi_ps_update" which can then be passed 
  ! conveniently as a block to the nubeam module...
 
  integer, parameter :: ispline   = 2  ! spline fit
  integer, parameter :: ihermite  = 1  ! hermite fit
  integer, parameter :: ilinear   = 0  ! piecewise linear

  integer, parameter :: ig00      = 1  ! f'=0 at axis, not-a-knot @edge
  integer, parameter :: iknot     = 0  ! not-a-knot bc, axis & edge

  integer, parameter :: izcen     = 2  ! zone centered quantity
  integer, parameter :: izbdy     = 1  ! bdy centered quanity
 
  integer :: inb,ing,infbe,iudsym,ii,ig,is,imp,id,iz,ia,ib,indx
  integer :: irho,inrho,inrhozc,inwk,inum,idrho,idrhozc
  integer :: inumr,inumz,itype,i,id_r,id_z,idlim
  integer :: iertmp

  character*32 zgxib,zgxizc,zgnam
  character*10 zid
  character*3 zid3

  real(fp), dimension(:), allocatable :: znsum,ztsum,zn1,zt1,zf1,zrgrid,zzgrid
  real(fp) :: zconv,znum
  real(fp) :: zrmin,zrmax,zzmin,zzmax,zrmid,zzmid,zdelr,zdelz
 
  integer :: ilun_dbg = 6

  !-----------------------------------------
  !  this item used to fetch grid size (nzones) if it is not already
  !  imbedded in the plasma state

  type(nbitype_grid) :: zgrid

  integer :: inzons,inzp1,inzp2
  real(fp), dimension(:), allocatable :: zrho_nbi,zrho_xzc

  real(fp), parameter :: ceps4 = 1.0e-4_fp
  integer, dimension(:), pointer :: idata,idata1
  real(fp), dimension(:), pointer :: r8data,r8data1
  integer :: iauth
  character*32 zauth

  !-----------------------------------------
  !  these items used if xplasma object needs to be extracted from plasma state

  type :: interp_object

     type (xplasma), pointer :: s_xpobj => null()

  end type interp_object

  type (interp_object) :: zobj
  type (xplasma), pointer :: sps

  !---------------------------------------------
  !  reset pscon structure elements to default marker values

  call nbi_ps_update_dflt(pscon)

  !---------------------------------------------
  !  get nubeam 1d gridsize

  inzp1 = 0
  inzp1 = ss%nrho_nbi
  if(inzp1.eq.0) inzp1 = ss%nrho_fus
  if(inzp1.eq.0) then
     call nbi_get_grid(zgrid)
     inzp1 = zgrid%nzones + 1
  endif

  inzons = inzp1 - 1
  inzp2 = inzp1 + 1

  !  rudimentary checks...

  if(inzp1.le.1) then
     write(ilun_dbg,*) ' ?nubeam_ctrl_psupdate: could not find radial grid.'
     ierr=1
     return
  endif
  
  if(ss%num_rzlim.eq.0) then
     write(ilun_dbg,*) ' ?nubeam_ctrl_psupdate: plasma state does not include a limiter or'
     write(ilun_dbg,*) '  vacuum vessel (r,z) contour (ss%num_rzlim = 0).'
     ierr=1
     return
  endif

  !  set up 1d grids: targets for profile interpolation to nubeam grid space
  !  and used for xplasma profile loading...

  allocate(zrho_nbi(inzp1),zrho_xzc(inzp2))

  ! regen these grids here, but, they will match what is in ss%rho_nbi
  ! which nubeam may have set up on an earlier call

  zrho_xzc(1)=zero
  zrho_xzc(inzp2)=one
  
  zrho_nbi(1)=zero
  do ii=2,inzp1
     zrho_nbi(ii) = (ii-1)*one/inzons
     zrho_xzc(ii) = half*(zrho_nbi(ii-1)+zrho_nbi(ii))
  enddo
  zrho_nbi(inzp1) = one  ! to be sure...

  ! zrho_nbi -- zone bdys; zrho_xzc -- zone ctrs + axis and boundary
  !---------------------------------------------
  !  scalar quantities

  !  time step
  if(ss%t1.le.ss%t0) then
     write(ilun_dbg,*) ' %nubeam_ctrl_psupdate: warning, positive time step required!'
     write(ilun_dbg,'(a,1pe11.4)') '  ss%t0 = ',ss%t0
     write(ilun_dbg,'(a,1pe11.4)') '  ss%t1 = ',ss%t1
     write(ilun_dbg,'(a,1pe11.4)') '  difference is: ',ss%t1-ss%t0
     !  allow to proceed-- nbstart will catch this...
     !!     ierr=1
     !!     return
  endif

  pscon%tbm1 = ss%t0
  pscon%tbm2 = ss%t1

  !  injection powers & voltages
  inb = ss%nbeam

  pscon%pinja(1:inb) = ss%power_nbi(1:inb)
  pscon%einja(1:inb) = ss%kvolt_nbi(1:inb)*1000.0_fp
  
  !  energy fractions (check against fraction tables if available)

  call ps_chk_nbfrac(0, ierr, state=ss)
  if(ierr.ne.0) return

  pscon%ffulla(1:inb) = ss%frac_full(1:inb)
  pscon%fhalfa(1:inb) = ss%frac_half(1:inb)

  !  handle trace elements in beams

  ii=inb
  do ib=1,inb
     if(ss%trace_flag(ib).ne.0) then
        ii=ii+1
        pscon%pinja(ii)=ss%power_nbi_trace(ib)
        pscon%einja(ii)=ss%kvolt_nbi(ib)*1000.0_fp
        if(ss%trace_flag(ib).eq.2) then
           pscon%ffulla(ii)=1.0_fp
           pscon%fhalfa(ii)=0.0_fp
        else
           pscon%ffulla(ii)=ss%frac_full(ib)
           pscon%fhalfa(ii)=ss%frac_half(ib)
        endif
     endif
  enddo

  !  sol model h0 neutral density

  pscon%dn0out = ss%dn0out * 1.0e-6_fp

  !  anom diffusivity energy dependence option

  pscon%nmdifb = ss%nmdifb
  infbe = ss%nefi_anom
  if(infbe.gt.0) then
     pscon%fdifbe(1:infbe) = ss%anom_evar(1:infbe)
  endif

  pscon%x_qlim0 = 0.0_fp  ! no x_qlim0 option when using ps
  ! instead the entire anomalous diffusivity profile is set for each
  ! class of fast ion species.

  ! set plhgt, if up-down symmetry is assumed in nubeam

  call nbstart_udsym(iudsym)

  if(iudsym.eq.1) then
     pscon%plhgt = ss%z_axis * 100.0_fp  ! -> cm
  else
     pscon%plhgt = 0.0_fp
  endif

  ! count number of non-impurity thermal ion species

  ing = 0
  do ii=1,ss%nspec_th
     if(ss%s_type(ii).ne.ps_impurity) then
        ing = ing + 1
     endif
  enddo

  pscon%ng = ing  ! in ps interface, inside nubeam, ng=ngmax always.

  !----------------------------------------------------------
  ! nthrf -- flag if all rf minority species are thermal
  !    (for now, nubeam has just one flag for all minority species; it is
  !    rare for there to be more than just one).

  pscon%nthrf = 1  ! assume thermal, initially
  ! (this value kept also if there are no rf minorities).

  do ii=1,ss%nspec_rfmin
     if(ss%isthermal(ii).ne.1) then
        pscon%nthrf = 2  ! non-thermal rf minority detected
     endif
  enddo

  if(imem_xplasma.le.0) then

     if(imem_xplasma.lt.0) then

        call xoi_init(ierr)  ! be sure "s" & "s2" pointers are initialized...

        call grab_f77_xplasma(s2, ierr)
        if(ierr.ne.0) then
           call nbi_warning( &
                ' ?nubeam_ctrl_psupdate: old_xplasma "s2" error (grab_f77_xplasma).')
           return
        endif
     endif

     !---------------------------------------------
     !  initialize xplasma instance; create id_mcgrid (2d irregular grid)
     !  for distribution function and related outputs.

     zobj = transfer(ss%iobj, zobj)
     sps => zobj%s_xpobj

     !  ...and copy this into old_xplasma instance "s"

     call xoi_init(ierr)  ! be sure "s" pointer is initialized...
     call xplasma_blkbxid(s,'aep',ido)
     call xplasma_blkbxid(sps,'aep',ido1)
     if(ido.gt.0.and.ido1.le.0) then
        call xplasma_blackbox_retrieve(s,ido,ierr,ibbtype, ia_ptr=idata, r8a_ptr=r8data)
        if(ierr.ne.0) then
           write(ilun_dbg,*) ' ?nubeam_ctrl_psupdate: unexpected error, blackbox aep:'
           call xplasma_error(s,ierr,ilun_dbg)
           return
        endif
        call xplasma_author_set(sps,'aep_transp',ierr)
        isize=idata(1)+1
        r8size=idata(9)
        call xplasma_create_blackbox(sps,'aep',ibbtype,ido1,ierr, &
             iasize=isize, r8asize=r8size,label='data for aep operator', &
             ia_ptr=idata1, r8a_ptr=r8data1)
        idata1=idata
        r8data1=r8data
        call xplasma_author_clear(sps,'aep_transp',ierr)
        
        if(ierr.ne.0) then
           write(ilun_dbg,*) &
                ' ?nubeam_ctrl_psupdate: unexpected error creating aep info object:'
           write(ilun_dbg,*) '  object name: aep'
           call xplasma_error(sps,ierr,ilun_dbg)
           return
        endif
        !  all done...
        nullify(idata1)  ! data has been given over to sps xplasma 
        nullify(r8data1) ! data has been given over to sps xplasma 

     endif
     call xplasma_copy(sps, s, ierr)
     if(ierr.ne.0) then
        call nbi_warning(' ?nubeam_ctrl_psupdate: unexpected error while copying plasma state xplasma data.')
        return
     endif
  endif

  pscon%id_mcgrid = 0   ! mcgrid to be produced in nbstart
  !  using xplasma data just copied.

  !---------------------------------------------
  !  ripple field(s) if any

  pscon%id_ripple = ss%id_tfrip_log
  pscon%id_ripple2 = ss%id_tfrip2_log

  !---------------------------------------------
  !  check [r,z] grids; if not provided from plasma state these need
  !  to be generated here...

  if(min(ss%nr,ss%nz).le.0) then

     !  (old xplasma calls can be used)

     inumr=max(81,3*(inzp1-1)+1)
     inumz=inumr

     !  (we also may need a limiter definition)

     call xplasma_find_item(s,'__limiter',idlim,ierr,nf_noerr=.true.)
     if(ierr.ne.0) then
        write(ilun_dbg,*) ' ?nubeam_ctrl_psupdate: unexpected xplasma error:'
        call xplasma_error(s,ierr,ilun_dbg)
        return
     endif

     if(idlim.eq.0) then

        call eqm_cbdy_grid(ss%num_rzlim, ss%rlim, ss%zlim, inumr, inumz, &
             id_r, id_z, ierr)
        if(ierr.ne.0) return

     else
        !  apparently limiter is present but not the [r,z] grids
        !  (code adapted from old_xplasma/eqm_cbdy_grid.f90)

        call eq_bdlims(itype,zrmin,zrmax,zzmin,zzmax,ierr)
        if(ierr.ne.0) return

        ! force z grid up-down symmetry
        zzmax = max(abs(zzmin),abs(zzmax))
        zzmin = -zzmax

        zrmid=half*(zrmin+zrmax)
        zzmid=half*(zzmin+zzmax)

        !  2a: expand slightly, but check...

        zdelr=(inumr-1)*half*(zrmax-zrmin)/(inumr-4)
        zdelz=(inumz-1)*half*(zzmax-zzmin)/(inumz-4)

        zrmin=max((zrmin/2),zrmid-zdelr)
        zrmax=zrmid+zdelr

        zzmin=zzmid-zdelz
        zzmax=zzmid+zdelz

        !  3.  generate the r,z grids

        allocate(zrgrid(inumr),zzgrid(inumz))

        do i=1,inumr
           zrgrid(i)=(zrmin*(inumr-i)+zrmax*(i-1))/(inumr-1)
        enddo
     
        do i=1,inumz
           zzgrid(i)=(zzmin*(inumz-i)+zzmax*(i-1))/(inumz-1)
        enddo

        call eqm_rzgrid(zrgrid,zzgrid,0,0,inumr,inumz,ceps4,id_r,id_z,ierr)
        if(ierr.ne.0) return

        deallocate(zrgrid,zzgrid)
     endif

     ! apply ad hoc extrapolation to define a field beyond plasma boundary

     call xplasma_brz(s,eqm_brz_adhoc,ierr)
     if(ierr.ne.0) then
        write(ilun_dbg,*) ' ?nubeam_ctrl_psupdate: unexpected xplasma error:'
        call xplasma_error(s,ierr,ilun_dbg)
        return
     endif

  endif
! find authorized author
  iauth=0
  zauth=' '
  if(ss%nrho_nbi.gt.0) then
     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
        call namels(iz,ia,zid)
        call xplasma_find_item(s,'eperp_'//trim(zid)//'_beam',id,ierr)
        if(id.gt.0) then
           call xplasma_get_item_info(s,id,ierr, author=zauth)
           exit
        endif
     enddo
  endif
  if(ss%nrho_fus.gt.0) then
     do ii=1,ss%nspec_fusion
        iz=ss%q_sfus(ii)/ps_xe + 0.1_fp
        ia=ss%m_sfus(ii)/ps_mp + 0.1_fp
        call namels(iz,ia,zid)
        call xplasma_find_item(s,'eperp_'//trim(zid)//'_fus',id,ierr)
        if(id.gt.0) then
           call xplasma_get_item_info(s,id,ierr, author=zauth)
           exit
        endif
     enddo
  endif

  if(zauth.ne.' ') then
     iauth = 1
     call xplasma_author_set(s,trim(zauth),ierr)
  endif
  if(iauth.eq.1) then
  !---------------------------------------------
  !  1d profile quantities

     inrho=inzp1
     inrhozc=inzp2
 
  !---------------------------
  !  make working grids

     zgxib = '_xib_nbdata'
     zgxizc= '_xizc_nbdata'


     call xplasma_gridid(s,zgxib,idrho)
     if(idrho.eq.0) then
        call xplasma_create_grid(s,trim(zgxib),xplasma_rho_coord,zrho_nbi, &
             idrho,ierr, label='rho bdy grid for zauth input profiles')
        if(ierr.ne.0) then
           call nbi_warning(' ?nubeam_ctrl_psupdate: unexpected xplasma grid creation error (1).')
           go to 999
        endif
     endif
     
     call xplasma_gridid(s,zgxizc,idrhozc)
     if(idrhozc.eq.0) then
        call xplasma_create_grid(s,trim(zgxizc),xplasma_rho_coord,zrho_xzc, &
          idrhozc,ierr, label='rho ctr grid for nubeam input profiles')
        if(ierr.ne.0) then
           call nbi_warning(' ?nubeam_ctrl_psupdate: unexpected xplasma grid creation error (2).')
           go to 999
        endif
     endif
     
  !---------------------------
  ! ps & xplasma units both mks so usually no conversion is needed...

     zconv = 1.0_fp
     ig=0
     do ii=1,ss%nspec_beam
        ig=ig+1
        iz=ss%qatom_snbi(ii)/ps_xe + 0.1_fp
        ia=ss%m_snbi(ii)/ps_mp + 0.1_fp
        call namels(iz,ia,zid)
        if(ss%nrho_nbi.gt.0) then
           call dload(ss%id_eperp_beami(ii),ihermite,ig00,izcen, &
                'eperp_'//trim(zid)//'_beam',pscon%id_eperpfast(ig))
           call dload(ss%id_epll_beami(ii),ihermite,ig00,izcen, &
                'epll_'//trim(zid)//'_beam',pscon%id_epllfast(ig))
        else
           call dload1(zn1,ihermite,ig00,izcen, &
                'eperp_'//trim(zid)//'_beam',pscon%id_eperpfast(ig))
           call dload1(zn1,ihermite,ig00,izcen, &
                'epll_'//trim(zid)//'_beam',pscon%id_epllfast(ig))
        endif
     enddo
     do ii=1,ss%nspec_fusion
        ig=ig+1
        iz=ss%q_sfus(ii)/ps_xe + 0.1_fp
        ia=ss%m_sfus(ii)/ps_mp + 0.1_fp
        call namels(iz,ia,zid)
        if(ss%nrho_fus.gt.0) then
           call dload(ss%id_eperp_fusi(ii),ihermite,ig00,izcen, &
                'eperp_'//trim(zid)//'_fus',pscon%id_eperpfast(ig))
           call dload(ss%id_epll_fusi(ii),ihermite,ig00,izcen, &
                'epll_'//trim(zid)//'_fus',pscon%id_epllfast(ig))
        else
           call dload1(zn1,ihermite,ig00,izcen, &
                'eperp_'//trim(zid)//'_fus',pscon%id_eperpfast(ig))
           call dload1(zn1,ihermite,ig00,izcen, &
             'epll_'//trim(zid)//'_fus',pscon%id_epllfast(ig))
        endif
     enddo
     if(iauth.eq.1) then
        call xplasma_author_clear(s,trim(zauth),ierr)
     endif

  endif
  !---------------------------------------------
  !  1d profile quantities

  inrho=inzp1
  inrhozc=inzp2
 
  !---------------------------
  !  make working grids

  zgxib = '_xib_nbdata'
  zgxizc= '_xizc_nbdata'

  call xplasma_author_set(s,'nubeam',iertmp)

  call xplasma_gridid(s,zgxib,idrho)
  if(idrho.eq.0) then
     call xplasma_create_grid(s,trim(zgxib),xplasma_rho_coord,zrho_nbi, &
          idrho,ierr, label='rho bdy grid for nubeam input profiles')
     if(ierr.ne.0) then
        call nbi_warning(' ?nubeam_ctrl_psupdate: unexpected xplasma grid creation error (1).')
        go to 999
     endif
  endif

  call xplasma_gridid(s,zgxizc,idrhozc)
  if(idrhozc.eq.0) then
     call xplasma_create_grid(s,trim(zgxizc),xplasma_rho_coord,zrho_xzc, &
          idrhozc,ierr, label='rho ctr grid for nubeam input profiles')
     if(ierr.ne.0) then
        call nbi_warning(' ?nubeam_ctrl_psupdate: unexpected xplasma grid creation error (2).')
        go to 999
     endif
  endif

  !---------------------------
  ! ps & xplasma units both mks so usually no conversion is needed...

  zconv = 1.0_fp

  !---------------------------
  if(iertmp.ne.0) then
     call xplasma_error(s,ierr,ilun_dbg)
     go to 999
  endif

  ! load te and ti
  call dload(ss%id_ts(0)  ,ihermite,ig00,izcen,'te_nbi',pscon%id_te)
  if(ierr.ne.0) go to 999

  call dload(ss%id_ti     ,ihermite,ig00,izcen,'ti_nbi',pscon%id_ti)
  if(ierr.ne.0) go to 999

  !---------------------------
  ! separate ti for fusion calculations -- from non-impurities only

  ierr=0
  allocate(znsum(inzons),ztsum(inzons)); znsum=0.0_fp; ztsum=0.0_fp
  allocate(zn1(inzons),zt1(inzons))
  do ii=1,ss%nspec_th
     if(ss%s_type(ii).ne.ps_impurity) then
        call ps_rho_rezone(zrho_nbi,ss%id_ns(ii),zn1,iertmp, state=ss)
        ierr=ierr+abs(iertmp)
        call ps_rho_rezone(zrho_nbi,ss%id_ts(ii),zt1,iertmp, state=ss)
        ierr=ierr+abs(iertmp)
        znsum=znsum+zn1
        ztsum=ztsum+zn1*zt1
     endif
  enddo

  if(minval(znsum).le.zero) then
     call ps_rho_rezone(zrho_nbi,ss%id_ti,zt1,iertmp, state=ss)
     ierr=ierr+abs(iertmp)
     ztsum=zt1          ! array assignment
  else
     ztsum=ztsum/znsum  ! array expression
  endif

  call dload1(ztsum,ihermite,ig00,izcen,'tmjsm_nbi',pscon%id_tifusn)
  deallocate(znsum,ztsum,zn1,zt1)
  if(ierr.ne.0) go to 999

  !---------------------------
  ! thermal plasma densities including impurities

  pscon%id_aimp=0  ! no avg z & a; full impurity list provided.
  pscon%id_zimp=0

  call dload(ss%id_ns(0),ihermite,ig00,izcen,'rhoel',pscon%id_ne)
  if(ierr.ne.0) go to 999

  imp=0
  ig=0
  do ii=1,ss%nspec_th
     if(ss%s_type(ii).ne.ps_impurity) then
        ig=ig+1
        iz=ss%q_s(ii)/ps_xe + 0.1_fp
        ia=ss%m_s(ii)/ps_mp + 0.1_fp
        call namels(iz,ia,zid)
        call dload(ss%id_ns(ii),ihermite,ig00,izcen, &
             'rhob_'//zid,pscon%id_nh(ig))
        if(ierr.ne.0) go to 999
     else
        imp=imp+1
        write(zid3,'(i3.3)') imp
        call dload(ss%id_ns(ii),ihermite,ig00,izcen, &
             'rhix'//zid3,pscon%id_nimp(imp))
        if(ierr.ne.0) go to 999
     endif
  enddo

  !---------------------------
  !  fast specie densities

  allocate(zn1(inzons)); zn1=zero

  do ii=1,ss%nspec_rfmin
     iz=ss%q_rfmin(ii)/ps_xe + 0.1_fp
     ia=ss%m_rfmin(ii)/ps_mp + 0.1_fp
     call namels(iz,ia,zid)
     if(ss%nrho_icrf.gt.0) then
        call dload(ss%id_nmini(ii),ihermite,ig00,izcen, &
             'rhmin_'//trim(zid),pscon%id_nrf(ii))
        call dload(ss%id_eperp_mini(ii),ihermite,ig00,izcen, &
             'eperp_'//trim(zid)//'_fppmin',pscon%id_eperpmin(ii))
        call dload(ss%id_epll_mini(ii),ihermite,ig00,izcen, &
             'epll_'//trim(zid)//'_fppmin',pscon%id_epllmin(ii))
     else
        call dload1(zn1,ihermite,ig00,izcen, &
             'rhmin_'//trim(zid),pscon%id_nrf(ii))
        call dload1(zn1,ihermite,ig00,izcen, &
             'eperp_'//trim(zid)//'_fppmin',pscon%id_eperpmin(ii))
        call dload1(zn1,ihermite,ig00,izcen, &
             'epll_'//trim(zid)//'_fppmin',pscon%id_epllmin(ii))
     endif
     if(ierr.ne.0) go to 999
  enddo

  ig=0
  do ii=1,ss%nspec_beam
     ig=ig+1
     iz=ss%qatom_snbi(ii)/ps_xe + 0.1_fp
     ia=ss%m_snbi(ii)/ps_mp + 0.1_fp
     call namels(iz,ia,zid)
     if(ss%nrho_nbi.gt.0) then
        call dload(ss%id_nbeami(ii),ihermite,ig00,izcen, &
             'bdens_'//trim(zid),pscon%id_nfast(ig))
     else
        call dload1(zn1,ihermite,ig00,izcen, &
             'bdens_'//trim(zid),pscon%id_nfast(ig))
     endif
     if(ierr.ne.0) go to 999

     ! dmc -- diffusivities vs. (x,e)

     if(allocated(ss%bidiff_dtrap)) then
        pscon%id_fd0(ig) = ss%id_bidiff_dtrap(ii)
        pscon%id_fdb(ig) = ss%id_bidiff_btrap(ii)
        pscon%id_fdp(ig) = ss%id_bidiff_bpass_co(ii)
        pscon%id_fdq(ig) = ss%id_bidiff_dpass_co(ii)
        pscon%id_fdr(ig) = ss%id_bidiff_bpass_ctr(ii)
        pscon%id_fds(ig) = ss%id_bidiff_dpass_ctr(ii)
     endif

  enddo

  do ii=1,ss%nspec_fusion
     ig=ig+1
     iz=ss%q_sfus(ii)/ps_xe + 0.1_fp
     ia=ss%m_sfus(ii)/ps_mp + 0.1_fp
     call namels(iz,ia,zid)
     if(ss%nrho_fus.gt.0) then
        call dload(ss%id_nfusi(ii),ihermite,ig00,izcen, &
             'fdens_'//trim(zid),pscon%id_nfast(ig))
     else
        call dload1(zn1,ihermite,ig00,izcen, &
             'fdens_'//trim(zid),pscon%id_nfast(ig))
     endif
     if(ierr.ne.0) go to 999

     ! dmc -- diffusivities vs. (x,e)

     if(allocated(ss%bidiff_dtrap)) then
        pscon%id_fd0(ig) = ss%id_fidiff_dtrap(ii)
        pscon%id_fdb(ig) = ss%id_fidiff_btrap(ii)
        pscon%id_fdp(ig) = ss%id_fidiff_bpass_co(ii)
        pscon%id_fdq(ig) = ss%id_fidiff_dpass_co(ii)
        pscon%id_fdr(ig) = ss%id_fidiff_bpass_ctr(ii)
        pscon%id_fds(ig) = ss%id_fidiff_dpass_ctr(ii)
     endif

  enddo

  deallocate(zn1)

  !---------------------------
  !  angular velocity
  
  call dload(ss%id_omegat,ihermite,iknot,izcen,'omegag',pscon%id_omega)
  if(ierr.ne.0) go to 999

  !---------------------------
  !  radial electrostatic potential

  zconv = 1000.0_fp  ! kvolts -> volts (nubeam expects volts) for this one
  !  zconv applied in dload

  call dload(ss%id_epot,ihermite,iknot,izbdy,'phiprg',pscon%id_epot)
  if(ierr.ne.0) go to 999
     
  zconv = 1.0_fp     ! restore

  !---------------------------
  !  toroidal loop voltage profile
 
  call dload(ss%id_v_loop,ihermite,iknot,izbdy,'vpoh',pscon%id_vpoh)
  if(ierr.ne.0) go to 999

  !---------------------------
  !  anomalous transport -- just use the pclin profiles from the plasma state

  pscon%id_difb=0
  pscon%id_velb=0

  pscon%id_difb_bi = ss%id_difb_nbi
  pscon%id_velb_bi = ss%id_velb_nbi

  pscon%id_difb_fi = ss%id_difb_fusi
  pscon%id_velb_fi = ss%id_velb_fusi

  !---------------------------
  !  neutral species:
  !  densities, temperatures = 3/2 <e0>, and avg toroidal angular velocities

  !  recycling and gasflow (edge) sources:

  pscon%id_n0 = 0
  pscon%id_t0 = 0
  pscon%id_omeg0 = 0

  allocate(zn1(inzons))

  do is=1,ss%ngsc0
     write(zid3,'(i3.3)') is      ! neutral source index
     znum = ss%sc0(is)            ! source magnitude (atoms/sec)

     do ig=1,ss%nspec_gas
        iz=ss%qatom_sgas(ig)/ps_xe + 0.1_fp
        ia=ss%m_sgas(ig)/ps_mp + 0.1_fp
        call namels(iz,ia,zid) ! neutral gas species

        zn1 = zero
        if(ss%nrho_gas.eq.0) then
           call dload1(zn1,ilinear,iknot,izcen, &
                'den0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_n0(is,ig))
           call dload1(zn1,ilinear,iknot,izcen, &
                't0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_t0(is,ig))
           call dload1(zn1,ilinear,iknot,izcen, &
                'omg0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_omeg0(is,ig))
        else

           call ps_rho_rezone(zrho_nbi,ss%id_n0norm(ig,is),zn1,ierr, state=ss)
           if(ierr.ne.0) go to 999
           zn1 = znum*zn1

           call dload1(zn1,ilinear,iknot,izcen, &
                'den0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_n0(is,ig))
           if(ierr.ne.0) go to 999

           call dload(ss%id_t0sc0(ig,is),ilinear,iknot,izcen, &
                't0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_t0(is,ig))
           if(ierr.ne.0) go to 999

           call dload(ss%id_omeg0sc0(ig,is),ilinear,iknot,izcen, &
                'omg0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_omeg0(is,ig))
           if(ierr.ne.0) go to 999
        endif

     enddo
  enddo

  !  beam halo

  indx=ss%ngsc0 + 1
  write(zid3,'(i3.3)') indx

  do ig=1,ss%nspec_gas
     iz=ss%qatom_sgas(ig)/ps_xe + 0.1_fp
     ia=ss%m_sgas(ig)/ps_mp + 0.1_fp
     call namels(iz,ia,zid) ! neutral gas species

     zn1 = zero
     if(ss%nrho_nbi.eq.0) then
        call dload1(zn1,ilinear,iknot,izcen, &
             'den0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_n0(is,ig))
        if(ierr.ne.0) go to 999
        call dload1(zn1,ilinear,iknot,izcen, &
             't0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_n0(is,ig))
        if(ierr.ne.0) go to 999
        call dload1(zn1,ilinear,iknot,izcen, &
             'omg0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_n0(is,ig))
        if(ierr.ne.0) go to 999
     else
        call dload(ss%id_n0_halo(ig),ilinear,iknot,izcen, &
             'den0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_n0(indx,ig))
        if(ierr.ne.0) go to 999

        call dload(ss%id_t0_halo(ig),ilinear,iknot,izcen, &
             't0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_t0(indx,ig))
        if(ierr.ne.0) go to 999

        call dload(ss%id_omeg0_halo(ig),ilinear,iknot,izcen, &
             'omg0mn_'//trim(zid)//'_'//trim(zid3),pscon%id_omeg0(indx,ig))
        if(ierr.ne.0) go to 999
     endif
  enddo

  deallocate(zn1)

999 continue

  if(allocated(zrho_nbi)) deallocate(zrho_nbi)
  if(allocated(zrho_xzc)) deallocate(zrho_xzc)
  call xplasma_author_clear(s,'nubeam',iertmp)

contains
  subroutine dload(id_ps,ifit,ibdy0,izc,name,id)

    ! load xplasma profile: idrho or idrhozc grid

    integer, intent(in) :: id_ps ! plasma state id in
    integer, intent(in) :: ifit  ! ispline, ihermite, or ilinear
    integer, intent(in) :: ibdy0 ! axial bc code
    integer, intent(in) :: izc   ! zone centered or bdy centered profile

    character*(*), intent(in) :: name
    integer, intent(out) :: id   ! xplasma id out

    if(izc.eq.izcen) then
       allocate(zf1(inzons))
       call ps_rho_rezone(zrho_nbi,id_ps,zf1,ierr, state=ss)
       if(ierr.eq.0) then
          call dload1(zconv*zf1,ifit,ibdy0,izc,name,id)
       else
          id=0
       endif
       deallocate(zf1)
    else
       allocate(zf1(inzp1))
       call ps_intrp_1d(zrho_nbi,id_ps,zf1,ierr, state=ss)
       if(ierr.eq.0) then
          call dload1(zconv*zf1,ifit,ibdy0,izc,name,id)
       else
          id=0
       endif
       deallocate(zf1)
    endif

  end subroutine dload

  subroutine dload1(zdata,ifit,ibdy0,izc,name,id)

    ! load xplasma profile: idrho or idrhozc grid

    real(fp), dimension(:) :: zdata  ! scaled input data
    integer, intent(in) :: ifit  ! ispline, ihermite, or ilinear
    integer, intent(in) :: ibdy0 ! axial bc code
    integer, intent(in) :: izc   ! zone centered or bdy centered profile

    character*(*), intent(in) :: name
    integer, intent(out) :: id   ! xplasma id out

    !------------------
    integer :: isize
    real(fp), dimension(:), allocatable :: zdatwk
    real(fp) :: zdmin,zdmax
    !------------------

    if(izc.eq.izcen) then
       isize = size(zdata)
       if(isize.ne.inzons) then
          write(ilun_dbg,*) ' ?nubeam_ctrl_psupdate(dload1) array size: ',isize
          write(ilun_dbg,*) '  it was expected to match inzons = ',inzons
          write(ilun_dbg,*) '  name: '//trim(name)
          ierr=1
          return
       endif

       allocate(zdatwk(inzp2))
       zdatwk(2:inzp1)=zdata

       zdmin = min(0.9_fp*zdatwk(2),1.1_fp*zdatwk(2))
       zdmax = max(0.9_fp*zdatwk(2),1.1_fp*zdatwk(2))

       zdatwk(1)=min(zdmax,max(zdmin, &
            zdatwk(2)-0.125_fp*(zdatwk(3)-zdatwk(2)) ))

       zdmin = min(0.1_fp*zdatwk(inzp1),2.0_fp*zdatwk(inzp1))
       zdmax = max(0.1_fp*zdatwk(inzp1),2.0_fp*zdatwk(inzp1))

       zdatwk(inzp2)=min(zdmax,max(zdmin, &
            zdatwk(inzp1) + 0.5_fp*(zdatwk(inzp1)-zdatwk(inzp1-1)) ))

       call dloadx(idrhozc,zdatwk,ifit,ibdy0,name,id)

       deallocate(zdatwk)

    else
       isize = size(zdata)
       if(isize.ne.inzp1) then
          write(ilun_dbg,*) ' ?nubeam_ctrl_psupdate(dload1) array size: ',isize
          write(ilun_dbg,*) '  it was expected to match inzons + 1 = ',inzp1
          write(ilun_dbg,*) '  name: '//trim(name)
          ierr=1
          return
       endif

       call dloadx(idrho,zdata,ifit,ibdy0,name,id)

    endif

  end subroutine dload1

  subroutine dloadx(idrho_use,zdata,ifit,ibdy0,name,id)

    ! load xplasma profile: any rho grid

    integer, intent(in) :: idrho_use
    real(fp), dimension(:) :: zdata
    integer, intent(in) :: ifit  ! ispline, ihermite, or ilinear
    integer, intent(in) :: ibdy0 ! axial bc code

    character*(*), intent(in) :: name
    integer, intent(out) :: id

    !-----------------
    integer :: ibc0
    !-----------------
    ! use zconv to scale the data

    ! bc (bdy cond) control has no effect for ifit=ilinear

    if((ibdy0.eq.ig00).and.(ifit.ne.ilinear)) then
       ibc0=1  ! f'=0 @ rho=0
    else
       ibc0=0  ! not a knot (or default for ilinear fit)
    endif

    call xplasma_create_1dprof(s,name,idrho_use,zdata, &
         id,ierr, &
         ispline=ifit, ibca=ibc0)

    if(ierr.ne.0) then
       write(ilun_dbg,*) '  error in nubeam_ctrl_psupdate(dloadx):'
       write(ilun_dbg,*) ' ?xplasma_create_1dprof failed on "',trim(name),'".'
       call xplasma_error(s,ierr,ilun_dbg)
    endif

  end subroutine dloadx

end subroutine nubeam_ctrl_psupdate

subroutine nubeam_ctrl_psupdate_cleanup(ierr)

  use xplasma_obj_instance
  implicit none

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

  integer, intent(out) :: ierr

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

  call release_f77_xplasma(ierr)
  if(ierr.ne.0) then
     call nbi_warning( &
          ' ?nubeam_ctrl_psupdate: old_xplasma error (release_f77_xplasma).')
  endif

end subroutine nubeam_ctrl_psupdate_cleanup
