  subroutine nbi_interp_profiles(ierr)
 
    implicit none
    integer, intent(out) :: ierr
 
    !  interapolate profiles; xplasma ids have been loaded.
 
    call nbi_1d_grid   ! 1d "xi" grid  xi = "rho" = sqrt.norm.tor.flux
 
    call nbi_getprofiles(ierr)   ! interpolate profiles
 
    if(ierr.eq.0) then
       call nbi_profcomck   ! adjust nbi_com scalars based on profile data.
    endif
 
  end subroutine nbi_interp_profiles
 
  subroutine nbi_chk_inputs(ierr)

    !  simple basic checks on input data-- done just once per run
    !  because only time invariant things are checked

    logical, save :: ichecked = .false.

    ierr=0
    if(.not.ichecked) then
       call nbi_datchk(ierr)
       if(ierr.eq.0) ichecked=.true.
    endif

  end subroutine nbi_chk_inputs

  subroutine nbi_getprofiles(ierr)
    !
    !  extracted from nbi_profiles module (dmc 12/16/01) so all of
    !  nbi_com can be used
 
    use nbi_com
    use xplasma_obj_instance
 
    !
    ! local variables
    !
    implicit none
    integer :: inumth, iaccur, id_xyz=0
    integer :: id_g=0, id_tflx=0
    !
    real(fp)    :: zdum, zndv, zntns(6), znlim
    real(fp) :: zrho_min, zrat, zne
    !
    real(fp),  allocatable   :: zwork(:,:)
    real(fp), allocatable  :: zwork8(:,:)
    integer,allocatable  :: ilist(:)
    real(fp), allocatable    :: zdens_fusn(:,:),zdens_fusn_rf(:,:)  ! d, t, he3 thermal and rf densities
    real(fp), allocatable    :: ti_mini(:,:)  ! d, t, he3 rf tmeperatures
    real(fp)                 :: ztol
    integer              :: inum, ierr, j, iz, ig, impt, is, ii, id, inst
    integer              :: ibuf(1)
    integer              :: isnccwi,isnccwb
    real(fp)                 :: zbuf(1), zrhix, zth
    real(fp)                 :: z1, z2
    real(fp),   allocatable  :: zbvec(:, :)
    real(fp),   allocatable  :: zrho(:), zchi(:,:), zphi(:), zzdum(:)

    real(fp), allocatable :: zneprof(:,:),zne_fdiff(:)
    real(fp) :: znesum,znesum_orig,zne_mindiff,znfrac(nsbeam)

    integer, parameter   :: izrmode = 1            ! output b in (r,z,phi)
    integer, allocatable :: inregion(:)
    integer :: inrhix,inrhmin,imin
    real(fp) zmin,ztim1,ztim2
    integer :: idchk1,idchk2,ifound_d2d
 
    logical :: ilfdens_override  ! =t to use fast densities from
    ! input data rather than internal state.  this may be needed
    ! when input test data is provided to nubeam for which the fast
    ! specie densities are preset.  nbi_getprofiles decides based on
    ! a quasi-neutrality test.

    logical :: ilfdens_avail     ! =t if input data fast density override
    ! is available

    logical :: ilfdens_quasi     ! =t to estimate fast density from 
    ! quasi-neutrality equation (error recovery available if nubeam
    ! namelist switch quasi_check=.false.; normally, a quasineutrality
    ! error stops the run).

    real(fp) :: zquasi_errtol
    integer :: iwarn

    real(fp) :: zdq1,zdq2,zdq3,zdfmax
    !-----------------------------------

    ierr=0

    if(only_io) then
       write(nonlin,*) ' ?nbi_getprofiles: improper call, only_io = .true.'
       ierr=666
       return
    endif

    !-----------------------
    ! init cpu timings.
    ! as this is the first "major" call to nubeam facilities, clear cpu
    ! measures here.  dep/orb/fin calculations cannot run without setting
    ! up profiles first.

    cpbstart = zero
    cpbdep = zero
    cpborb = zero
    cpbfin = zero
    cpbmax = zero
    cpbmin = zero

    call cpu_time(ztim1)
    !-----------------------
    !  apply sawtooth correction first, if needed...

    if(sawflag.or.sawflag2) call sawnbi_prof

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

    call nbi_chk_inputs(ierr)
    if(ierr.ne.0) return
 
    !  (xplasma) flux surface label (rho) range
    !
    call eq_rholim(rho_axis,rho_bdy)
    !
    !-----------------------------------
    !  compute volumes and zonal volume averages
    !
    inumth=1                          ! integrate around whole zone
    iaccur=1                          ! accuracy control (1 recommended)
    !
    zrho_min=1.0e-7_fp
    call eq_flxint_init(0,nzp1,xi_r8(lcentr:lep1,1),zrho_min,ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint_init'
       return
    endif
    !...............
    call eq_flxint('dvol', iaccur, dvolzc_r8(lcentr:ledge), inumth, nzp1-1, ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint("dvol",...)'
       return
    endif
 
    dvolzc(lcentr:ledge) = dvolzc_r8(lcentr:ledge)
 
    allocate(zwork8(mj,5)); zwork8=zero
 
    !...............
    !  (dmc march 2002 -- use surface integrals -- for speed)
 
    call eq_flxint('<b^2>s', iaccur, zwork8(lcentr:lep1,1), inumth, nzp1, ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint("<b^2>",...)'
       return
    endif
    !...............
 
    call eq_flxint('<r^2>s', iaccur, zwork8(lcentr:lep1,2), inumth, nzp1, ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint("<r^2>",...)'
       return
    endif
    !...............
 
    call eq_flxint('<1/r>s', iaccur, zwork8(lcentr:lep1,3), inumth, nzp1, ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint("<1/r>",...)'
       return
    endif
    !...............
 
    call eq_flxint('<1/r^2>s', iaccur, zwork8(lcentr:lep1,4), inumth, nzp1, ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint("<1/r^2>",...)'
       return
    endif
    !...............
 
    call eq_flxint('<grad(rho)>s', iaccur, zwork8(lcentr:lep1,5), inumth, nzp1, ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint("<grad(rho)>",...)'
       return
    endif
 
    !...............
    ! shift to zone ctr
 
    do id=1,5
       do ii=lcentr,ledge
          zwork8(ii,id)=(zwork8(ii,id)+zwork8(ii+1,id))*half
       enddo
    enddo
 
    !...............
    ! copy
 
    ! note (dmc) -- code is all real(fp) now; these copies may no longer
    ! be needed...

    gb2zc(lcentr:ledge)=zwork8(lcentr:ledge,1)
    gr2zc(lcentr:ledge)=zwork8(lcentr:ledge,2)
    grizc(lcentr:ledge)=zwork8(lcentr:ledge,3)
    gr2izc(lcentr:ledge)=zwork8(lcentr:ledge,4)
    dravizc(lcentr:ledge)=zwork8(lcentr:ledge,5)
 
    dareazc(lcentr:ledge)=dvolzc(lcentr:ledge)*grizc(lcentr:ledge)/twopi
    dareazc_r8(lcentr:ledge)= &
         dvolzc_r8(lcentr:ledge)*zwork8(lcentr:ledge,3)/twopi
 
    deallocate(zwork8)
    allocate(zwork8(mj,2)); zwork8=zero

    !...............set nc 

    call eq_flxint('<b>s', iaccur, zwork8(lcentr:lep1,1), inumth, nzp1, ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint("<b>",...)'
       return
    endif
    !...............
    call eq_flxint('<f(h)nclass>s', iaccur, zwork8(lcentr:lep1,2), inumth, nzp1, ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  error in eq_flxint("<f(h)nclass>",...)'
       return
    endif
    !...............
    do id=1,2
       do ii=lcentr,ledge
          zwork8(ii,id)=(zwork8(ii,id)+zwork8(ii+1,id))*half
       enddo
    enddo
 
    !...............
    ! copy
 
    gb1zc(lcentr:ledge)=zwork8(lcentr:ledge,1)
    gncfhzc(lcentr:ledge)=zwork8(lcentr:ledge,2)
    deallocate(zwork8)

    !...............
    ! set edge value & fix units
 
    dvolzc(lep1)      = two*dvolzc(lep1-1)-dvolzc(lep1-2)
    dvolzc            = 1.e6_fp*dvolzc               ! m**3 -> cm**3
    dvolzc_r8(lep1)      = two*dvolzc_r8(lep1-1)-dvolzc_r8(lep1-2)
    dvolzc_r8            = 1.e6_fp*dvolzc_r8               ! m**3 -> cm**3
    !
    gr2zc(lep1)       = two*gr2zc(lep1-1)-gr2zc(lep1-2)
    gr2zc             = 1.e4_fp*gr2zc                ! m^2 -> cm^2; gr2zc(:,1) is not used by nubeam
    !
    grizc(lep1)       = two*grizc(lep1-1)-grizc(lep1-2)
    grizc             = 1.e-2_fp*grizc               ! 1/m -> 1/cm; grizc(:, 1) is not used by nubeam
    !
    dravizc(lcentr:ledge)= 1.0e-2_fp*dravizc(lcentr:ledge)/dxi(lcentr:ledge,2)
    dravizc(lep1)        = dravizc(ledge)
    !
    dareazc(lep1)     = two*dareazc(lep1-1)-dareazc(lep1-2)
    dareazc           = 1.e4_fp*dareazc              ! m**2 -> cm**2
    dareazc_r8(lep1)     = two*dareazc_r8(lep1-1)-dareazc_r8(lep1-2)
    dareazc_r8           = 1.e4_fp*dareazc_r8              ! m**2 -> cm**2
    !
    gr2izc(lep1)      = two*gr2izc(lep1-1)-gr2izc(lep1-2)
    gr2izc            = 1.e-4_fp*gr2izc              ! 1/m**2 -> 1/cm**2
 
    !-----------------------------------
    !  get toroidal flux from xplasma; define fbz
    !  get poloidal flux; define qgeo numerically: d[psi-pol]/d[psi-tor].
    !
    allocate(zwork(mj,2))
 
    call eq_gfnum('TFLUX', id_tflx)
    if(id_tflx.eq.0) then
       !  if tflux is undefined: define it.
       call eq_gfnum('G',id_g)
       call eq_rgetf(nzones, xi(lcentr:ledge,2), id_g, 0, &
            zwork(lcentr:ledge,2),   ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  error in eq_rgetf("g",...) 0'
          return
       endif
       tflux = 0         ! integrate (r*bt)*<1/r**2>*dvol/2pi (toroidal flux)
       zwork(lcentr,1)=0
       do j=lcentr,ledge
          tflux = tflux + zwork(j,2)*1.0e-2_fp*gr2izc(j)*dvolzc(j)/twopi
          zwork(j+1,1)=tflux
       enddo

       call nbo_load1_xizb  ! define bdy-oriented grid (id_rhozb)

       !  spline; d(tflux)/drho -> 0 on axis; default at edge.
       call eqm_rhofun(2,id_rhozb,'tflux',zwork(lcentr:lep1,1), &
            1, zero, 0, zero, id_tflx, ierr)
       if(ierr.ne.0) then
          write(nonlin,*) ' ?nbi_getprofiles: tflux construction failed.'
          return
       endif
    endif
 
    call eq_rgetf(nzp1, xi(lcentr:lep1,1), id_tflx, 0, &
         zwork(lcentr:lep1,1),    ierr)
    if(ierr.ne.0) then
       write(nonlin,*)  &
            ' ?nbi_getprofiles:  error in eq_rgetf("tflux",...)'
       return
    endif
    tflux            = zwork(lep1,1)
 
    call eq_gfnum('PSI',id_psi)
    call eq_rgetf(nzp1, xi(lcentr:lep1,1), id_psi, 0, &
         zwork(lcentr:lep1,2),    ierr)
    if(ierr.ne.0) then
       write(nonlin,*)  &
         ' ?nbi_getprofiles:  error in eq_rgetf("psi",...)'
       return
    endif
    plflxg(lcentr:lep1)=zwork(lcentr:lep1,2)
 
    qgeo(lcentr+1:ledge)= &
         (zwork(lcentr+2:ledge+1,1)-zwork(lcentr:ledge-1,1))/ &
         (twopi*(plflxg(lcentr+2:ledge+1)-plflxg(lcentr:ledge-1)))
    qgeo(lcentr)=qgeo(lcentr+1)-(qgeo(lcentr+2)-qgeo(lcentr+1))/3
 
    ! set q @ bdy, using a constrained extrapolation

    zdq1 = qgeo(ledge-1)-qgeo(ledge-2)
    zdq2 = qgeo(ledge)-qgeo(ledge-1)
    if(zdq2.le.zero) then 
       zdq3=zero
    else
       zdq3 = min(two*zdq2,max(zdq2,(two*zdq2-zdq1)))
    endif

    qgeo(lep1) = qgeo(ledge)+zdq3

    deallocate(zwork)
    !-----------------------------------
    !  get r*bt from xplasma, & fbz
    !
    call eq_gfnum('G', id_g)
    call eq_rgetf(nzp1, xi(lcentr:lep1,1), id_g, 0, fbz(lcentr:lep1,1), &
         ierr)
    if(ierr.ne.0) then
       write(nonlin,*)  &
         ' ?nbi_getprofiles:  error in eq_rgetf("g",...) 1'
       return
    endif
    bzxr              = fbz(lep1,1)
    call eq_rgetf(nzones, xi(lcentr:ledge,2), id_g, 0, fbz(lcentr:ledge,2),&
         ierr)
    if(ierr.ne.0) then
       write(nonlin,*)  &
         ' ?nbi_getprofiles:  error in eq_rgetf("g",...) 2'
       return
    endif
 
    ! normalize to edge value (bzxr, tesla-m at the moment)
 
    fbz(lcm1,2)=fbz(lcentr,2)
    fbz(lep1,2)=fbz(lep1,1)
    fbz = fbz/bzxr
    bzxr              = 1.0e2_fp*bzxr  ! convert to tesla*cm
    !-----------------------------------
    !  b & ip direction information
    !
    call eq_bchk_sign(isnccwb,isnccwi)
    nsnccw=isnccwb
    nsjdotb=isnccwb*isnccwi
    !
    !-----------------------------------
    !  get min(mod(b)), max(mod(b)), and 1/mod(b)
    !  and corresponding r locations (convert to cm) from xplasma
    !
 
    allocate(zrho(lep1),zchi(lep1,2),zphi(lep1),zzdum(lep1),inregion(lep1))
 
    call eq_gfnum('BMOD',id_xyz)
    zdum=zero
    call eq_frhochi(1,rho_axis,zdum,1,id_xyz,1,bminm(lcentr),ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  eq_frhochi(bmod) error.'
       return
    endif
    bmaxm(lcentr)=bminm(lcentr)  ! on axis
    zrho(lcentr)=zero; zchi=zero; zphi(lcentr)=zero
 
    !
    do is=lcp1,lep1
       zrho(is)=xi(is,1)
       zphi(is)=zero
       call eq_glimb1(zrho(is),zphi(is),bminm(is),zchi(is,1), &
            bmaxm(is),zchi(is,2),ierr)
       if(ierr.ne.0) then
          write(nonlin,*) ' ?nbi_getprofiles:  eq_glimb1 error.'
          return
       endif
    enddo
 
    call eq_rzget(nzones+1,zrho(lcentr:lep1),zchi(lcentr:lep1,1), &
         zphi(lcentr:lep1),rbminm(lcentr:lep1),zzdum,inregion,ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  eq_rzget (@bminm) error.'
       return
    endif
    rbminm=rbminm*100
 
    call eq_rzget(nzones+1,zrho(lcentr:lep1),zchi(lcentr:lep1,2), &
         zphi(lcentr:lep1),rbmaxm(lcentr:lep1),zzdum,inregion,ierr)
    if(ierr.ne.0) then
       write(nonlin,*) ' ?nbi_getprofiles:  eq_rzget (@bmaxm) error.'
       return
    endif
    rbmaxm=rbmaxm*100._fp
 
    thbminm(lcentr:lep1)=zchi(lcentr:lep1,1)
    thbmaxm(lcentr:lep1)=zchi(lcentr:lep1,2)
 
    deallocate(zrho,zchi,zphi,zzdum,inregion)
    !
    !--------------------------------------------------
    allocate(zwork(nzp1, 4), ilist(4))
    ilist(1) = id_te
    ilist(2) = id_ne
    ilist(3) = id_ti
    ilist(4) = id_tifusn
    inum=4
    if(ilist(4).eq.0) inum=3
    !
    call eq_rgetff(nzp1,xi(lcentr:lep1,1),inum,ilist,0,nzp1,zwork,ierr)
    if(ierr.ne.0) then
       write(nonlin,*)  &
         ' ?nbi_getprofiles:  unexpected:  eq_rgetff error (te/ti/ne) 1.'
       return
    endif
    te(lcentr:lep1, 1)      = 1000._fp*zwork(1:nzp1,1)  ! te
    rhoel(lcentr:lep1, 1)   = 1.e-6_fp*zwork(1:nzp1,2)  ! ne
    ti(lcentr:lep1, 1)      = 1000._fp*zwork(1:nzp1,3)  ! ti
    if(inum.eq.4) then
       if(maxval(ti(:,1)).le.zero) then
          ti(lcentr:ledge, 1)=1000._fp*zwork(1:nzp1-1,4)
       endif
    endif
    !
    call eq_rgetff(nzp1-1,xi(lcentr:ledge,2),inum,ilist,0,nzp1-1, &
         zwork(1:nzp1-1,1:inum),ierr)
    if(ierr.ne.0) then
       write(nonlin,*)  &
         ' ?nbi_getprofiles:  unexpected:  eq_rgetff error (te/ti/ne) 2.'
       return
    endif
    te(lcentr:ledge, 2)      = 1000._fp*zwork(1:nzp1-1,1)  ! te
    rhoel(lcentr:ledge, 2)   = 1.e-6_fp*zwork(1:nzp1-1,2)  ! ne
    ti(lcentr:ledge, 2)      = 1000._fp*zwork(1:nzp1-1,3)  ! ti
    if(inum.eq.4) then
       tmjsm(lcentr:ledge)       = 1000._fp*zwork(1:nzp1-1,4)  ! tmajority
       if(maxval(ti(:,2)).le.zero) then
          ti(lcentr:ledge, 2)=tmjsm(lcentr:ledge)
       endif
    else
       tmjsm=ti(1:mj,2)
    endif
    !
    te(lep1, 2)      = max(half*te(lep1, 1), two*te(lep1, 1)-te(ledge, 2))
    te(lcm1, 1)      = te(lcp1, 1)   ; te(lcm1, 2)	= te(lcentr, 2)
    rhoel(lep1, 2)   = max(half*rhoel(lep1, 1), two*rhoel(lep1, 1)-rhoel(ledge, 2))
    rhoel(lcm1, 1)   = rhoel(lcp1, 1); rhoel(lcm1, 2)	= rhoel(lcentr, 2)
    ti(lep1, 2)      = max(half*ti(lep1, 1), two*ti(lep1, 1)-ti(ledge, 2))
    ti(lcm1, 1)      = ti(lcp1, 1)   ; ti(lcm1, 2)	= ti(lcentr, 2)
    tmjsm(lcm1)           = tmjsm(lcentr)
    tmjsm(lep1)           = max(tmjsm(lep1-1), two*tmjsm(lep1-1)-tmjsm(ledge-1))
    if(id_aimp.gt.0) then
       ilist(1)=id_nimp(1)
       ilist(2)=id_zimp
       ilist(3)=id_aimp
       inum=3
       call eq_rgetff(nzp1,xi(lcentr:lep1,1),inum,ilist,0,nzp1,zwork,ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
            ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: nx/ax/zx.'
          return
       endif
       rhi(lcentr:lep1,1) = 1.e-6_fp*zwork(1:nzp1,1)
       call eq_rgetff(nzp1-1,xi(lcentr:ledge,2),inum,ilist,0,nzp1-1, &
            zwork(1:nzp1-1,1:inum),ierr)
       rhi(lcentr:ledge,2) = 1.e-6_fp*zwork(1:nzp1-1,1)
       rhi(lep1, 2)     = max(half*rhi(lep1, 1), two*rhi(lep1, 1)-rhi(ledge, 2))
       rhi(lcm1, 1)     = rhi(lcp1, 1)  ; rhi(lcm1, 2)	= rhi(lcentr, 2)
       xzimpj(lcentr:ledge)  = zwork(1:nzp1-1,2); xzimpj(lcm1) = xzimpj(lcentr)
       xzimpj(lep1)     = max(xzimpj(lep1-1), two*xzimpj(lep1-1)-xzimpj(ledge-1))
       xzimpj(lcm1)=xzimpj(lcentr)
       aimpj(lcentr:ledge)   = zwork(1:nzp1-1,3); aimpj(lcm1)  = aimpj(lcentr)
       aimpj(lep1)      = max(aimpj(lep1-1), two*aimpj(lep1-1)-aimpj(ledge-1))
       aimpj(lcm1)=aimpj(lcentr)
       if(nrhix.eq.1) then
          !
          ! tokamakium -- radially varying <a> and <z> impurity
          ! compute #-weighted profile averaged <<a>> and <<z>>, scalars.
          !
          aimpx(1)=zero
          xzimpx(1)=zero
          zdum=zero
          do ii=lcentr,ledge
             zndv=max(one,rhi(ii,2))*dvolzc(ii)
             zdum=zdum+zndv
             aimpx(1)=aimpx(1)+zndv*aimpj(ii)
             xzimpx(1)=xzimpx(1)+zndv*xzimpj(ii)
          enddo
          aimpx(1)=aimpx(1)/zdum
          xzimpx(1)=xzimpx(1)/zdum
       endif
    else
       if(nrhix.eq.0) then
          rhi=zero
          xzimpj=6.0_fp
          aimpj=12.0_fp
       endif
    endif
!
    deallocate(zwork, ilist)
    !--------------------------------------------------
    ! background plasma density at bdy/zone j
 
    ! here also: start density sum for thermonuclear fusion
    allocate(zdens_fusn(3,mj),zdens_fusn_rf(3,mj),ti_mini(3,mj))
    zdens_fusn=zero
    zdens_fusn_rf=zero
    ti_mini=zero
    rhob=zero
    if (ng>0) then
       allocate(zwork(nzp1, ng))
       call eq_rgetff(nzp1,xi(lcentr:lep1,1),ng,id_nh(1:ng),0, &
            nzp1,zwork,ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
            ' ?nbi_getprofiles:  unexpected:  eq_rgetff error (species ni).'
          return
       endif
       rhob(lcentr:lep1, 1:ng, 1)   = 1.e-6_fp*zwork(1:nzp1, 1:ng)
       call eq_rgetff(nzp1-1,xi(lcentr:ledge,2),ng,id_nh(1:ng),0, &
            nzp1-1,zwork(1:nzp1-1,1:ng),ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
            ' ?nbi_getprofiles:  unexpected:  eq_rgetff error (species ni) 2.'
          return
       endif
       rhob(lcentr:ledge, 1:ng, 2)  = 1.e-6_fp*zwork(1:nzp1-1, 1:ng)
       do ig = 1, ng
          rhob(lep1, ig, 2) = &
               max(half*rhob(lep1,ig,1), two*rhob(lep1, ig, 1)-rhob(ledge, ig, 2))
       end do
       rhob(lcm1, 1:ng, 1) = rhob(lcp1, 1:ng, 1)
       rhob(lcm1, 1:ng, 2) = rhob(lcentr, 1:ng, 2)
       deallocate(zwork)
       do ig=1,ng
          if(ig.eq.ndeut) zdens_fusn(1,lcentr:ledge)=rhob(lcentr:ledge,ig,2)
          if(ig.eq.ntrit) zdens_fusn(2,lcentr:ledge)=rhob(lcentr:ledge,ig,2)
          if(ig.eq.nhel3) zdens_fusn(3,lcentr:ledge)=rhob(lcentr:ledge,ig,2)
       enddo
      !
       den0mn=zero
       en0mn=zero
       omg0mn=zero
       inst=mis
       allocate(zwork(ledge,inst))
       do ig=1, ng
          do is=1, inst
             if(id_n0(is,ig).eq.0) then
                zwork(lcentr:ledge,is)=zero
             else
                call eq_rgetf(nzones,xi(lcentr:ledge,2), &
                     id_n0(is,ig),0,zwork(lcentr:ledge,is),ierr)
                if(ierr.ne.0) then
                   write(nonlin,*)  &
                        ' ?nbi_getprofiles:  unexpected:  eq_rgetf error: n0.'
                   return
                endif
             endif
             den0mn(ig, lcentr:ledge, is)   = 1.e-6_fp*zwork(lcentr:ledge, is)
          enddo
       enddo
       !
       do ig=1, ng
          do is=1, inst
             if(id_t0(is,ig).eq.0) then
                zwork(lcentr:ledge,is)=0
             else
                call eq_rgetf(nzones,xi(lcentr:ledge,2), &
                     id_t0(is,ig),0,zwork(lcentr:ledge,is),ierr)
                if(ierr.ne.0) then
                   write(nonlin,*)  &
                        ' ?nbi_getprofiles:  unexpected:  eq_rgetf error: t0.'
                   return
                endif
             endif
             en0mn(ig, lcentr:ledge, is)   = 1.e3_fp*zwork(lcentr:ledge, is)
          enddo
       enddo
       !
       do ig=1, ng
          do is=1, inst
             if(id_omeg0(is,ig).eq.0) then
                zwork(lcentr:ledge,is)=zero
             else
                call eq_rgetf(nzones,xi(lcentr:ledge,2), &
                     id_omeg0(is,ig),0,zwork(lcentr:ledge,is),ierr)
                if(ierr.ne.0) then
                   write(nonlin,*)  &
                       ' ?nbi_getprofiles: unexpected:  eq_rgetf error: omeg0.'
                   return
                endif
             endif
             omg0mn(ig, lcentr:ledge, is)   = one*zwork(lcentr:ledge, is)
          enddo
       enddo
       !
       deallocate(zwork)
    endif
    !--------------------------------------------------
    ! multiple impurity density at bdy/zone j used in code
    rhix=zero
    if ((nrhix.gt.1).or.((nrhix.eq.1).and.(id_aimp.eq.0))) then
       allocate(zwork(nzp1, nrhix))
       call eq_rgetff(nzp1,xi(lcentr:lep1,1),nrhix,id_nimp(1:nrhix), &
            0,nzp1,zwork,ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
            ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: nx.'
          return
       endif
       do ig = 1, nrhix
          rhix(lcentr:lep1, 1, ig)   = 1.e-6_fp*zwork(1:nzp1, ig)
       enddo
       call eq_rgetff(nzp1-1,xi(lcentr:ledge,2),nrhix,id_nimp(1:nrhix), &
            0,nzp1-1,zwork(1:nzp1-1,1:nrhix),ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
            ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: nx 2.'
          return
       endif
       do ig = 1, nrhix
          rhix(lcentr:ledge, 2, ig)   = 1.e-6_fp*zwork(1:nzp1-1, ig)
       enddo
      !
       do ig = 1, nrhix
          rhix(lep1, 2, ig) = max(half*rhix(lep1, 1, ig), &
               two*rhix(lep1, 1, ig)-rhix(ledge, 2, ig))
       enddo
       rhix(lcm1, 1, 1:nrhix) = rhix(lcp1, 1, 1:nrhix)
       rhix(lcm1, 2, 1:nrhix) = rhix(lcentr, 2, 1:nrhix)
       !
       !  <z> and <a>
       !
       do ii=lcentr,lep1
          rhi(ii,1)=sum(rhix(ii,2,1:nrhix))
          zrhix=sum(rhix(ii,2,1:nrhix))
          rhi(ii,2)=zrhix
          if(zrhix.gt.0) then
             aimpj(ii)=sum(rhix(ii,2,1:nrhix)*aimpx(1:nrhix))/zrhix
             xzimpj(ii)=sum(rhix(ii,2,1:nrhix)*xzimpx(1:nrhix))/zrhix
          else
             aimpj(ii)=aimpx(1)
             xzimpj(ii)=xzimpx(1)
          endif
       enddo
       rhi(lcm1,2)=rhi(lcentr,2)
       xzimpj(lcm1)=xzimpj(lcentr)
       aimpj(lcm1)=aimpj(lcentr)
       deallocate(zwork)
    endif
    !--------------------------------------------------
    rhmin=zero
    timin=zero
    if(nmini>0) then
       allocate(zwork(nzp1,nmini))
       call eq_rgetff(nzp1,xi(lcentr:lep1,1),nmini,id_nrf(1:nmini), &
            0,nzp1,zwork,ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: nmini.'
          return
       endif
       do ig=1,nmini
          rhmin(lcentr:lep1,1,ig)=1.e-6_fp*zwork(1:nzp1,ig)
       enddo
       call eq_rgetff(nzp1-1,xi(lcentr:ledge,2),nmini,id_nrf(1:nmini), &
            0,nzp1-1,zwork(1:nzp1-1,1:nmini),ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: nmini 2.'
          return
       endif
       do ig = 1, nmini
          rhmin(lcentr:ledge, 2, ig)   = 1.e-6_fp*zwork(1:nzp1-1, ig)
       enddo
       do ig = 1, nmini
          rhmin(lep1, 2, ig) = max(half*rhmin(lep1, 1, ig), &
               two*rhmin(lep1, 1, ig)-rhmin(ledge, 2, ig))
       enddo
       rhmin(lcm1, 1, 1:nmini) = rhmin(lcp1, 1, 1:nmini)
       rhmin(lcm1, 2, 1:nmini) = rhmin(lcentr, 2, 1:nmini)
! mg added 08/2019, get temperture of rf minority
       call eq_rgetff(nzp1,xi(lcentr:lep1,1),nmini,id_eperpmin(1:nmini), &
            0,nzp1,zwork,ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: eperpimini'
          return
       endif
       timin(lcentr:lep1,1,1:nmini)=zwork(1:nzp1,1:nmini) 
       call eq_rgetff(nzp1,xi(lcentr:lep1,1),nmini,id_epllmin(1:nmini), &
            0,nzp1,zwork,ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: epllmin.'
          return
       endif
       timin(lcentr:lep1,1,1:nmini)=(timin(lcentr:lep1,1,1:nmini)+zwork(1:nzp1,1:nmini))*1.e+3_fp*2._fp/3._fp
       call eq_rgetff(nzp1-1,xi(lcentr:ledge,2),nmini,id_eperpmin(1:nmini), &
            0,nzp1-1,zwork(1:nzp1-1,1:nmini),ierr)
       timin(lcentr:ledge,2,1:nmini)=zwork(1:nzp1-1,1:nmini) 
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: eperpmin 2.'
          return
       endif
       call eq_rgetff(nzp1-1,xi(lcentr:ledge,2),nmini,id_epllmin(1:nmini), &
            0,nzp1-1,zwork(1:nzp1-1,1:nmini),ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  unexpected:  eq_rgetff error: epllmin 2.'
          return
       endif
       timin(lcentr:ledge,2,1:nmini)= (timin(lcentr:ledge,2,1:nmini)+zwork(1:nzp1-1,1:nmini))*1.e+3_fp*2._fp/3._fp 
       do ig = 1, nmini
          timin(lep1, 2, ig) = max(half*timin(lep1, 1, ig), &
               two*timin(lep1, 1, ig)-timin(ledge, 2, ig))
       enddo
       timin(lcm1, 1, 1:nmini) = timin(lcp1, 1, 1:nmini)
       timin(lcm1, 2, 1:nmini) = timin(lcentr, 2, 1:nmini)
       
       deallocate(zwork)
      
       do ig=1,nmini
          if(ig.eq.mrfdeut) then
             zdens_fusn_rf(1,lcentr:ledge)=rhmin(lcentr:ledge,2,ig)
             ti_mini(1,lcentr:ledge)=timin(lcentr:ledge,2,ig)
          endif
          if(ig.eq.mrftrit) then
             zdens_fusn_rf(2,lcentr:ledge)=rhmin(lcentr:ledge,2,ig)
             ti_mini(2,lcentr:ledge)=timin(lcentr:ledge,2,ig)
          endif
          if(ig.eq.mrfhel3) then
             zdens_fusn_rf(3,lcentr:ledge)=rhmin(lcentr:ledge,2,ig)
             ti_mini(3,lcentr:ledge)=timin(lcentr:ledge,2,ig)
          endif
       enddo
    endif

    ! initially: assume fast density is taken from internal state

    ilfdens_override = .false.
    ilfdens_avail = .true.  ! this is tested
    ilfdens_quasi = .false.

    allocate(zneprof(mj,2),zne_fdiff(mj))
    do
       rhbs=zero
       if(abs(ntest_max).eq.1) then
          ebavps=zero
          ebapp2s=zero
          ebapl2s=zero
          do is=1,nsbeam
             call eq_rgetf(nzp1,xi(lcentr:lep1,1),id_eperpfast(is),0, &
                  ebapp2s(lcentr:lep1,is),ierr)
             call eq_rgetf(nzp1,xi(lcentr:lep1,1),id_epllfast(is),0, &
                  ebapl2s(lcentr:lep1,is),ierr)
             ebapp2s=ebapp2s*1.0e+3_fp !to ev
             ebapl2s=ebapl2s*1.0e+3_fp
             ebavps(lcentr:lep1,is)=ebapp2s(lcentr:lep1,is)+ebapl2s(lcentr:lep1,is)
          enddo
       endif
       do is=1,nsbeam
          if(nlfbmfpp(is).or.ilfdens_override) then
             call eq_rgetf(nzp1,xi(lcentr:lep1,1),id_nfast(is),0, &
                  rhbs(lcentr:lep1,1,is),ierr)
             if(ierr.ne.0) then
                write(nonlin,*)  &
                     ' ?nbi_getprofiles:  unexpected:  eq_rgetf error: nfast'
                return
             endif
             rhbs(lcentr:lep1,1,is)=1.e-6_fp*rhbs(lcentr:lep1,1,is)
             call eq_rgetf(nzones,xi(lcentr:ledge,2),id_nfast(is),0, &
                  rhbs(lcentr:ledge,2,is),ierr)
             if(ierr.ne.0) then
                write(nonlin,*)  &
                     ' ?nbi_getprofiles:  unexpected:  eq_rgetf error: nfast 2'
                return
             endif
             rhbs(lcentr:ledge,2,is)=1.e-6_fp*rhbs(lcentr:ledge,2,is)
          else
             ! check if input data could be used in place of bdenss(...)
             if(id_nfast(is).eq.0) ilfdens_avail=.false.

             rhbs(lcentr:ledge,2,is)=bdenss(lcentr:ledge,is)/ &
                  dvolzc(lcentr:ledge)
             rhbs(lcentr,1,is)=rhbs(lcentr,2,is)
             rhbs(lcp1:ledge,1,is)=half*(rhbs(lcentr:ledge-1,2,is)+ &
                  rhbs(lcp1:ledge,2,is))
             rhbs(lep1,1,is)=max(half*rhbs(ledge,2,is), &
                  two*rhbs(ledge,2,is)-rhbs(ledge,1,is))
          endif
       enddo
       do is = 1, nsbeam
          rhbs(lep1, 2, is) = max(half*rhbs(lep1, 1, is), &
               two*rhbs(lep1, 1, is)-rhbs(ledge, 2, is))
          rhbs(lcm1, 1, 1:nsbeam) = rhbs(lcp1, 1, 1:nsbeam)
          rhbs(lcm1, 2, 1:nsbeam) = rhbs(lcentr, 2, 1:nsbeam)
       enddo
       !--------------------------------------------------
       ! apply fixes to density profiles
       !
       ! thermal plasma -- max dynamic range 10**4 in density profiles
       ! fast species -- max dynamic range 10**8 in density profiles
       ! neutral atoms -- density.ge. 1 ptcl/cc
 
       znlim=1.0e-4_fp*max(zero,maxval(rhoel(lcm1:lep1,1:2)))
       rhoel(lcm1:lep1,1:2)=max(znlim,rhoel(lcm1:lep1,1:2))
 
       do ig=1,ng
          znlim=1.0e-4_fp*max(zero,maxval(rhob(lcm1:lep1,ig,1:2)))
          rhob(lcm1:lep1,ig,1:2)=max(znlim,rhob(lcm1:lep1,ig,1:2))
       enddo
 
       if(nrhix.gt.1) then
          do ig=1,nrhix
             znlim=1.0e-4_fp*max(zero,maxval(rhix(lcm1:lep1,1:2,ig)))
             rhix(lcm1:lep1,1:2,ig)=max(znlim,rhix(lcm1:lep1,1:2,ig))
          enddo
       else
          znlim=1.0e-4_fp*max(zero,maxval(rhi(lcm1:lep1,1:2)))
          rhi(lcm1:lep1,1:2)=max(znlim,rhi(lcm1:lep1,1:2))
       endif
 
       do ig=1,nmini
          znlim=1.0e-8_fp*max(zero,maxval(rhmin(lcm1:lep1,1:2,ig)))
          rhmin(lcm1:lep1,1:2,ig)=max(znlim,rhmin(lcm1:lep1,1:2,ig))
       enddo
 
       do ig=1,nsbeam
          znlim=1.0e-8_fp*max(zero,maxval(rhbs(lcm1:lep1,1:2,ig)))
          rhbs(lcm1:lep1,1:2,ig)=max(znlim,rhbs(lcm1:lep1,1:2,ig))
       enddo
 
       inst=mis
       do is=1,inst
          do ig=1,ng
             znlim=zero
             if(maxval(den0mn(ig,lcentr:ledge,is)).gt.zero) znlim=one
             den0mn(ig,lcentr:ledge,is)=max(znlim,den0mn(ig,lcentr:ledge,is))
          enddo
       enddo
 
       !--------------------------------------------------
       ! moments of the impurity density wrt
       ! z**1, z**2, a, rmass, z**2/a at bdy/zone j
       !
       do ii = lcm1, lep1
          rhimom(ii,2,1) = rhi(ii,2)*xzimpj(ii)
          if(nrhix.gt.1) then
             rhimom(ii,2,2) = sum(rhix(ii,2,1:nrhix)*xzimpx(1:nrhix)**2)
             rhimom(ii,2,5) = sum(rhix(ii,2,1:nrhix)*(xzimpx(1:nrhix)**2/ &
                  aimpx(1:nrhix)))
          else
             rhimom(ii,2,2) = rhi(ii,2)*xzimpj(ii)**2
             rhimom(ii,2,5) = rhi(ii,2)*xzimpj(ii)**2/aimpj(ii)
          endif
          rhimom(ii,2,3) = rhi(ii,2)*aimpj(ii)
          rhimom(ii,2,4) = rhi(ii,2)*aimpj(ii)*zmp
          if (ii > lcm1) then
             do iz = 1, 5
                rhimom(ii,1,iz) = half*(rhimom(ii-1,2,iz)+rhimom(ii,2,iz))
             end do
          endif
       enddo
       !--------------------------------------------------
       ! compute zeff; check quasi-neutrality against input ne profile (rhoel)

       zneprof = zero
       xzeffc = zero

       znesum = zero
       znesum_orig = zero

       iz = 2

       if(ilfdens_avail) then
          zquasi_errtol=0.001_fp
       else
          zquasi_errtol=0.05_fp
       endif

       iwarn = 0
       ierr = 0

       do ii = lcentr,lep1
          if(ii.le.ledge) then
             znesum_orig = znesum_orig + dvolzc(ii)*rhoel(ii,iz)
          endif

          zne = zero
          do ig = 1, ng
             zne = zne + rhob(ii, ig, iz) * backz(ig)
             xzeffc(ii, iz) = xzeffc(ii, iz) + rhob(ii, ig, iz) * backz(ig)**2
          enddo
          zne = zne + rhimom(ii,2,1)
          xzeffc(ii,iz) = xzeffc(ii,iz) + rhimom(ii,2,2)
          do ig = 1, nmini
             zne = zne + rhmin(ii, iz, ig) * xzmini(ig)
             xzeffc(ii, iz)= xzeffc(ii, iz) + rhmin(ii, iz, ig) * xzmini(ig)**2
          enddo

          zne_fdiff(ii) = rhoel(ii,iz) - zne

          do ig = 1, nsbeam
             zne = zne + rhbs(ii, iz, ig) * xzbeams(ig)
             xzeffc(ii, iz)= xzeffc(ii, iz) + rhbs(ii, iz, ig) * xzbeams(ig)**2
          enddo

          zneprof(ii,iz) = zne
          if(ii.le.ledge) then
             znesum = znesum + dvolzc(ii) * zne
          endif
          
          xzeffc(ii, iz) = xzeffc(ii, iz) / zne
       enddo

       write(nonlin,'(1x,a,1pe11.4)') &
            '   nbi_getprofiles ne*dvol sum (input):',znesum_orig
       write(nonlin,'(1x,a,1pe11.4)') &
            '   nbi_getprofiles ne*dvol sum (ions): ',znesum
       if(abs(one-znesum/znesum_orig).gt.zquasi_errtol) then
          if(ilfdens_avail) then
             iwarn = 1
             write(nonlin,*) ' %nbi_getprofiles: quasineutrality check => '
             write(nonlin,*) '  use externally provided fast ion densities.'
          else
             if(quasi_check) then
                write(nonlin,*) ' ?nbi_getprofiles: quasineutrality check error.'
                ierr = ierr + 1
             else
                iwarn = 2
                write(nonlin,*) ' %nbi_getprofiles: quasi_check=.false.:'
                write(nonlin,*) '  attempt to set fast fast density from quasi-neutrality deficit.'
                call qcfdens_set
                quasi_check=.true.
             endif
          endif
       endif

       if(iwarn.eq.1) then
          ilfdens_override=.true.
          ilfdens_avail=.false.
          write(nonlin,*) ' ==> retry using input data for fast densities.'
          cycle
       else if((iwarn.eq.2).and.(ierr.eq.0)) then
          ilfdens_override=.false.
          ilfdens_avail=.false.
          write(nonlin,*) ' ==> retry using patched data for fast densities.'
          cycle
       endif

       exit
    enddo
    if(ierr.ne.0) return

    xzeffc(lcm1,2)=xzeffc(lcentr,2)
    zneprof(lcm1,2)=zneprof(lcentr,2)
    do ii = lcentr,lep1
       xzeffc(ii,1)=half*(xzeffc(ii-1,2)+xzeffc(ii,2))
       zneprof(ii,1)=half*(zneprof(ii-1,2)+zneprof(ii,2))
    enddo
    !--------------------------------------------------
    !  ** rhoel replaced **
    !  use zi*ni sum

    rhoel = zneprof
    deallocate(zneprof,zne_fdiff)

    !--------------------------------------------------
    call eq_rgetf(nzp1-1,xi(lcentr:ledge,2), id_vpoh, 0, &
         vpoh(lcentr:ledge), ierr)
    if(ierr.ne.0) then
       write(nonlin,*)  &
            ' ?nbi_getprofiles:  eq_rgetf error:  vpoh (loop voltage).'
       return
    endif
    omegag=zero
    if(id_omega.gt.0) then
       call eq_rgetf(nzp1-1,xi(lcentr:ledge,2), id_omega, 0, &
            omegag(lcentr:ledge), ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  eq_rgetf error:  omegag (ang.velocity).'
          return
       endif
    endif
    curt=zero
    if(id_curt.gt.0) then
       call eq_rgetf(nzp1,xi(lcentr:lep1,1), id_curt, 0, &
            curt(lcentr:lep1), ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
            ' ?nbi_getprofiles:  eq_rgetf error: curt (area integral toroidal current).'
          return
       endif
    else
       ! data not provided; compute from equilibrium, b.dl loop integral
       call eq_flxint('itor',iaccur,curt(lcentr:lep1),inumth,nzp1,ierr)
       if(ierr.ne.0) then
          write(nonlin,*) &
               ' ?nbi_getprofiles: eq_flxint error, evaluating curt (itor(rho))'
          return
       endif
    endif

    ! dmc june 2009 -- constrain curt edge value
    zdq1 = curt(ledge-1)-curt(ledge-2)
    zdq2 = curt(ledge)-curt(ledge-1)
    zdq3 = two*zdq2-zdq1
    zdfmax = two/nzones
    curt(lep1)=curt(ledge)+zdq3
    if(curt(lep1).gt.(one+zdfmax)*curt(ledge)) then
       curt(lep1) = (one+zdfmax)*curt(ledge)
       write(nonlin,*) ' %note: constrained curt @ bdy to: ',curt(lep1),' a.'
    endif
    if(curt(lep1).lt.(one-zdfmax)*curt(ledge)) then
       curt(lep1) = (one+zdfmax)*curt(ledge)
       write(nonlin,*) ' %note: constrained curt @ bdy to: ',curt(lep1),' a.'
    endif

    phiprg=zero
    id_phi = id_epot   ! these are synonyms
    if(id_epot.gt.0) then
       call eq_rgetf(nzp1,xi(lcentr:lep1,1), id_epot, 0, &
            phiprg(lcentr:lep1), ierr)
       if(ierr.ne.0) then
          write(nonlin,*)  &
               ' ?nbi_getprofiles:  eq_rgetf error: phiprg (electrostatic potential).'
          return
       endif
    endif

    ! 1d anomalous transport

    idchk1 = max(id_velb,id_difb)
    idchk2 = max(id_velb_bi,id_difb_bi,id_velb_fi,id_difb_fi)
    if((idchk1.gt.0).and.(idchk2.gt.0)) then
       ierr=1
       write(nonlin,*) &
            ' ?nbi_getprofiles: if one or more of {id_velb,id_difb} is set,'
       write(nonlin,*) &
            '  none of {id_velb_fi,id_difb_fi,id_velb_bi,id_difb_bi} may be set.'
       return
    endif

    velb=zero
    if(id_velb.gt.0) then
       call vdget(id_velb,velb,1.e2_fp,0,"velb (anomalous velocity)")
       if(ierr.ne.0) return
    endif

    difb=zero
    if(id_difb.gt.0) then
       call vdget(id_difb,difb,1.e4_fp,1,"difb (anomalous diffusivity)")
       if(ierr.ne.0) return
    endif

    velb_fi=zero
    if(id_velb_fi.gt.0) then
       call vdget(id_velb_fi,velb_fi,1.e2_fp,0,"velb_fi (anomalous velocity)")
       if(ierr.ne.0) return
    endif

    difb_fi=zero
    if(id_difb_fi.gt.0) then
       call vdget(id_difb_fi,difb_fi,1.e4_fp,1,"difb_fi (anomalous diffusivity)")
       if(ierr.ne.0) return
    endif

    velb_bi=zero
    if(id_velb_bi.gt.0) then
       call vdget(id_velb_bi,velb_bi,1.e2_fp,0,"velb_bi (anomalous velocity)")
       if(ierr.ne.0) return
    endif

    difb_bi=zero
    if(id_difb_bi.gt.0) then
       call vdget(id_difb_bi,difb_bi,1.e4_fp,1,"difb_bi (anomalous diffusivity)")
       if(ierr.ne.0) return
    endif

    ! 2d anomalous transport: d(x,e) for orbit topology limits...

    if(ndifbe_2d.gt.1) then

       if(max(idchk1,idchk2).gt.0) then
          write(nonlin,*) &
               ' *** nbi_getprofiles: detected both 1d & 2d anomalous'
          write(nonlin,*) &
               '     transport data!  1d data ignored.'
       endif

       difbe_2d_fd0 = zero
       difbe_2d_fdb = zero
       difbe_2d_fdp = zero
       difbe_2d_fdq = zero
       difbe_2d_fdr = zero
       difbe_2d_fds = zero

       ifound_d2d = 0

       do is=1,nsbeam
          if(id_fd0(is).gt.0) then
             call d2load('fd0',is,id_fd0(is),difbe_2d_fd0(:,:,is))
             if(ierr.ne.0) exit
          endif
          if(id_fdb(is).gt.0) then
             call d2load('fdb',is,id_fdb(is),difbe_2d_fdb(:,:,is))
             if(ierr.ne.0) exit
          endif
          if(id_fdp(is).gt.0) then
             call d2load('fdp',is,id_fdp(is),difbe_2d_fdp(:,:,is))
             if(ierr.ne.0) exit
          endif
          if(id_fdq(is).gt.0) then
             call d2load('fdq',is,id_fdq(is),difbe_2d_fdq(:,:,is))
             if(ierr.ne.0) exit
          endif
          if(id_fdr(is).gt.0) then
             call d2load('fdr',is,id_fdr(is),difbe_2d_fdr(:,:,is))
             if(ierr.ne.0) exit
          endif
          if(id_fds(is).gt.0) then
             call d2load('fds',is,id_fds(is),difbe_2d_fds(:,:,is))
             if(ierr.ne.0) exit
          endif
       enddo

    endif

    if(ierr.ne.0) return

    !--------------------------------------------------
    ! apply fixes to temperature profiles
    !   charged ptcls: at least 10 ev
    !   neutral atoms: at least 2 ev
 
    te(lcentr:lep1,1:2)=max(10.0_fp,te(lcentr:lep1,1:2))
    ti(lcentr:lep1,1:2)=max(10.0_fp,ti(lcentr:lep1,1:2))
    tmjsm(lcentr:lep1)=max(10.0_fp,tmjsm(lcentr:lep1))
 
    inst=mis
    do is=1,inst
       do ig=1,ng
          en0mn(ig,lcentr:ledge,is)=max(2.0_fp,en0mn(ig,lcentr:ledge,is))
       enddo
    enddo
 
    !--------------------------------------------------
    ! compute thntns
    thntns=zero
    do ii=lcentr,ledge
       call nbi_thfusn(tmjsm(ii),zdens_fusn(1:3,ii),ti_mini(1:3,ii),zdens_fusn_rf(1:3,ii),zntns,&
            lev_ttnuk,nthrf)
       thntns(ii,1:6)=zntns(1:6)
    enddo
    !
    deallocate(zdens_fusn,zdens_fusn_rf,ti_mini)
    !

    call cpu_time(ztim2)
    cpbstart = cpbstart + (ztim2-ztim1)

  contains

    subroutine vdget(id,zprof,zconv,ipos,zlabel)

      use iso_c_binding, only: fp => c_double
      use physconst_mod, only: half, zero
      integer, intent(in) :: id      ! integer id
      real(fp), dimension(:) :: zprof  ! array to receive profile data
      real(fp), intent(in) :: zconv    ! convert, mks -> nubeam internal units
      integer :: ipos                ! =1 to force positive (diffusivities)
      character*(*), intent(in) :: zlabel   ! name & description of profile

      call eq_rgetf(nzp1,xi(lcentr:lep1,1), id, 0, &
            zprof(lcentr:lep1), ierr)

      if(ierr.ne.0) then
         write(nonlin,*)  &
              ' ?nbi_getprofiles:  eq_rgetf error: '//trim(zlabel)
         return
      endif
      zprof(lcentr:lep1)=zconv*zprof(lcentr:lep1)

      if(ipos.eq.0) return

      !  find minimum positive value  -- will enforce db > 0 everywhere
 
      imin=0
      zmin=epsinv
      do ii=lcentr,lep1
         if(zprof(ii).gt.zero) then
            if(zprof(ii).lt.zmin) then
               imin=ii
               zmin=zprof(ii)  ! minimum positive value
            endif
         endif
      enddo
      if(imin.eq.0) then
         zprof=zero  ! no positive values
      else
         zmin=min(0.01_fp,half*zmin)
         zprof(lcentr:lep1)=max(zmin,zprof(lcentr:lep1))
      endif

    end subroutine vdget

    subroutine d2load(zlbl,isb,id,arr2d)
      ! load 2d d(x,e) anom diffusivity data
      use iso_c_binding, only: fp => c_double
      use physconst_mod, only : z2ev=>zkev2ev, z2kev=>zev2kev
      character*(*), intent(in) :: zlbl  ! data label
      integer, intent(in) :: isb         ! species index
      integer, intent(in) :: id          ! xplasma id
 
      ! intent(inout):
      real(fp), dimension(:,:) :: arr2d    ! the data to load if id.gt.0

      !-------------------------
      real(fp) :: zetest(med2d),zetol,zevec(lcentr:lep1)
      integer :: idg1,idg2,irank,isize,icoord,icoord_e,idc_rho
      integer :: ii,ie,idg
      !-------------------------
      !  array arg check:

      ierr=0

      if(size(arr2d,1).ne.mj) then
         write(nonlin,*) ' ?nbi_getprofiles(d2load): 1st dim size .ne.mj:', &
              size(arr2d,1)
         ierr=1
         return
      endif

      if(size(arr2d,2).ne.med2d) then
         write(nonlin,*) ' ?nbi_getprofiles(d2load): 2nd dim size .ne.med2d:',&
              size(arr2d,2)
         ierr=1
         return
      endif

      !  id check: if zero, ok, no data for this array section...
      if(id.eq.0) return

      !-------------------------------
      !  non-zero data id provided...

      ifound_d2d = ifound_d2d + 1

      call xplasma_prof_info(s,id,ierr, rank=irank)
      if(ierr.ne.0) then
         call xplasma_error(s,ierr,nonlin)
         write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
              ' 2d diffusivity, specie is=',is
         return
      endif

      if(irank.ne.2) then
         write(nonlin,*) ' ?nbi_getprofiles: rank=2 expected, got rank=',irank
         write(nonlin,*) ' error with '//zlbl//' 2d diffusivity, specie is=',is
         ierr=1
         return
      endif

      call xplasma_prof_info(s,id,ierr, gridid1=idg1, gridid2=idg2)
      if(ierr.ne.0) then
         call xplasma_error(s,ierr,nonlin)
         write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
              ' 2d diffusivity, specie is=',is
         return
      endif

      idc_rho=0
      icoord_e=0

      do ii=1,irank
         if(ii.eq.1) idg=idg1
         if(ii.eq.2) idg=idg2

         call xplasma_grid_info(s,idg,ierr, size=isize, coord=icoord)
         if(ierr.ne.0) then
            call xplasma_error(s,ierr,nonlin)
            write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
                 ' 2d diffusivity, specie is=',is
            write(nonlin,*) '  grid lookup error.'
            return
         endif

         if(icoord.eq.xplasma_rho_coord) then
            idc_rho=icoord

         else
            icoord_e=icoord
            if(isize.ne.med2d) then
               write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
                    ' 2d diffusivity, specie is=',is
               write(nonlin,*) '  energy grid size ',med2d, &
                    ' expected, instead got: ',isize
               ierr=1
               return
            endif

            call xplasma_grid(s,idg,zetest,ierr)
            if(ierr.ne.0) then
               call xplasma_error(s,ierr,nonlin)
               write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
                    ' 2d diffusivity, specie is=',is
               write(nonlin,*) '  energy grid fetch error.'
               return
            endif

            zetest = zetest*z2ev  ! -> ev

            if(ifound_d2d.eq.1) then
               edifbe_2d(1:med2d) = zetest(1:med2d)

            else
               zetol = 1.0e-10_fp*maxval(abs(edifbe_2d))
               if(maxval(abs(zetest-edifbe_2d)).gt.zetol) then
                  write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
                    ' 2d diffusivity, specie is=',is
                  write(nonlin,*) '  energy grid value check error.'
                  write(nonlin,*) '  all energy grid must match, but: '
                  write(nonlin,*) '    #  reference  current_value'
                  do ie=1,med2d
                     write(nonlin,'(2x,i4,2(2x,1pe15.8))') ie, &
                          edifbe_2d(ie),zetest(ie)
                  enddo
                  ierr=1
                  return
               endif
            endif
         endif
      enddo

      if(idc_rho.eq.0) then
         write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
              ' 2d diffusivity, specie is=',is
         write(nonlin,*) '  rho coordinate not found.'
         ierr=1
         return
      endif

      if(icoord_e.eq.0) then
         write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
              ' 2d diffusivity, specie is=',is
         write(nonlin,*) '  energy coordinate not found.'
         ierr=1
         return
      endif

      ! ok: at last, interpolate the data...

      do ie=1,med2d
         zevec(lcentr:lep1) = edifbe_2d(ie)*z2kev  ! all the same value, kev...
         call xplasma_eval_prof(s,id, &
              idc_rho,  xizb(lcentr:lep1), &
              icoord_e, zevec(lcentr:lep1), &
              arr2d(lcentr:lep1,ie), &
              ierr)

         if(ierr.ne.0) then
            call xplasma_error(s,ierr,nonlin)
            write(nonlin,*) ' ?nbi_getprofiles: error with '//zlbl// &
                 ' 2d diffusivity, specie is=',is
            write(nonlin,*) '  data fetch error.'
            return
         endif

      enddo

      ! apply units conversion for diffusivity
      arr2d = arr2d*dif2d_conv

    end subroutine d2load

    subroutine qcfdens_set

      ! set fast density fractions from input power
      ! (quasi_check=f recovery for failed quasineutrality check)

      use iso_c_binding, only: fp => c_double
      use physconst_mod, only: zero
      integer :: ib,isb
      real(fp) :: zptot,zftot

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

      ierr = 0

      do j=lcentr,ledge
         if(zne_fdiff(j).lt.-1.0e-6_fp*zneprof(j,2)) then
            if(ierr.eq.0) then
               write(nonlin,*) ' ?negative quasineutrality gap cannot be filled with added fast ion density.'
               ierr=1
            endif
         else
            zne_fdiff(j) = max(zero,zne_fdiff(j))
         endif
      enddo
      if(ierr.ne.0) return

      znfrac = zero
      zptot = one  ! sum + 1.0

      ! beam power fractions

      do ib=1,nbeam
         zptot = zptot + pinja(ib)
         do isb=1,nsbeam
            if(max(abs(xzbeama(ib)-xzbeams(isb)), &
                 abs(abeama(ib)-abeams(isb))).lt.0.1_fp) then
               znfrac(isb)=znfrac(isb) + pinja(ib)
            endif
         enddo
      enddo

      ! nominal small fraction for fusion productions and/or beam species
      ! with no current injected power

      zftot = zero
      do isb=1,nsbeam
         znfrac(isb)=max(zptot*1.0e-4_fp,znfrac(isb))
         zftot=zftot+znfrac(isb)
      enddo
      znfrac = znfrac/zftot

      ! apply fractions to bridge quasineutrality gap

      do j=lcentr,ledge
         do isb=1,nsbeam
            bdenss(j,isb)=dvolzc(j)*znfrac(isb)*zne_fdiff(j)/xzbeams(isb)
         enddo
      enddo

    end subroutine qcfdens_set

  end subroutine nbi_getprofiles



  subroutine nbi_thfusn(zti,zdens,zti_rf,zdens_rf,zrdens,ilev_ttnuk,n_thrf)
 
    use iso_c_binding, only: fp => c_double
    use physconst_mod, only: zero,half
    implicit none
    integer,intent(in):: ilev_ttnuk,n_thrf
    real(fp), intent(in) :: zti           ! ion temperature, ev
    real(fp), intent(inout) :: zdens(3)      ! nd,nt,nhe3, cm**-3, respectively
    real(fp), intent(in) :: zti_rf(3)          !rf- ion temperature, ev
    real(fp), intent(in) :: zdens_rf(3)      ! rf - ion nd,nt,nhe3, cm**-3, respectively
    real(fp), intent(out) :: zrdens(6)    ! reactions/cm3/sec for:
 
    ! (dmc) 3 jan 2001 -- using preact calls get thermonuclear fusion
    ! .....(might want to move this to preact if precision options added)
 
 
                         ! dt,dhe3,dd(p),dd(n),tt, [t-he3 not available yet]
                         ! respectively
 
    !-------------------------
 
    real(fp) ztikev(1),zsigv
 
    !-------------------------
 
    zrdens=zero
 
    ztikev=0.001_fp*zti
    if(n_thrf.eq.1) then !treat minorities as thermal bulk
       zdens(1:3)= zdens(1:3)+zdens_rf(1:3)
    endif
    if(zdens(1)*zdens(2).gt.zero) then
       ! dt
       call fusn_thermo_dt(ztikev,zsigv,1)  ! returns zsigv in m**3/sec
       zrdens(1)=zdens(1)*zdens(2)*zsigv ! #/m3/sec
       if(n_thrf.ne.1) then
          if(zdens_rf(1)*zdens(2).gt.zero) then
             call fusn_btsigv_dt(zti_rf(1)*0.001_fp,ztikev,zsigv,1)
             zrdens(1)= zrdens(1)+zdens_rf(1)*zdens(2)*zsigv
          endif
          if(zdens_rf(1)*zdens_rf(2).gt.zero) then
             call fusn_btsigv_dt(zti_rf(1)*0.001_fp,zti_rf(2)*0.001_fp,zsigv,1)
             zrdens(1)= zrdens(1)+zdens_rf(1)*zdens_rf(2)*zsigv
          endif
       endif
       zrdens(1)= zrdens(1)*1.0e6_fp  ! #/cm3/sec
    endif
 
    if(zdens(1)*zdens(3).gt.zero) then
       ! dhe3 (gammas, not neutrons)
       call fusn_thermo_d3(ztikev,zsigv,1)  ! returns zsigv in m**3/sec
       zrdens(2)=zdens(1)*zdens(3)*zsigv ! #/m3/sec
       if(n_thrf.ne.1) then
          if(zdens_rf(1)*zdens(3).gt.zero) then
             call fusn_btsigv_d3(zti_rf(1)*0.001,ztikev,zsigv,1)
             zrdens(2)= zrdens(2)+zdens_rf(1)*zdens(3)*zsigv
          endif
          if(zdens_rf(1)*zdens_rf(3).gt.zero) then
             call fusn_btsigv_d3(zti_rf(1)*0.001,zti_rf(3)*0.001,zsigv,1)
             zrdens(2)= zrdens(2)+zdens_rf(1)*zdens_rf(3)*zsigv
          endif
       endif
       zrdens(2)= zrdens(2)*1.0e6_fp  ! #/cm3/sec
    endif

    if(zdens(1).gt.zero) then
       ! dd(p)
       call fusn_thermo_ddp(ztikev,zsigv,1) ! returns zsigv in m**3/sec
       zrdens(3)=zdens(1)*zdens(1)*zsigv! #/m3/sec
       zrdens(3)=half*zrdens(3)                 ! 0.5 due to self collision
       if(n_thrf.ne.1) then
          if(zdens_rf(1).gt.zero) then
             call fusn_btsigv_ddp(zti_rf(1)*0.001_fp,ztikev,zsigv,1)
             zrdens(3)= zrdens(3)+zdens_rf(1)*zdens(1)*zsigv
             call fusn_btsigv_ddp(zti_rf(1)*0.001_fp,zti_rf(1)*0.001_fp,zsigv,1)
             zrdens(3)= zrdens(3)+half*zdens_rf(1)*zdens_rf(1)*zsigv  ! 0.5 due to self collision
          endif
       endif
       zrdens(3)= zrdens(3)*1.0e6_fp ! #/cm3/sec
    endif
    if(zdens(1).gt.zero) then
       ! dd(n)
       call fusn_thermo_ddn(ztikev,zsigv,1) ! returns zsigv in m**3/sec
       zrdens(4)=zdens(1)*zdens(1)*zsigv ! #/m3/sec
       zrdens(4)=half*zrdens(4)                 ! 0.5 due to self collision
       if(n_thrf.ne.1) then
          if(zdens_rf(1).gt.zero) then
             call fusn_btsigv_ddn(zti_rf(1)*0.001_fp,ztikev,zsigv,1)
             zrdens(4)= zrdens(4)+zdens_rf(1)*zdens(1)*zsigv
             call fusn_btsigv_ddn(zti_rf(1)*0.001_fp,zti_rf(1)*0.001_fp,zsigv,1)
             zrdens(4)= zrdens(4)+half*zdens_rf(1)*zdens_rf(1)*zsigv ! 0.5 due to self collision
          endif
       endif
       zrdens(4)= zrdens(4)*1.0e6_fp  ! #/cm3/sec
    endif
 
    if(zdens(2).gt.zero) then
       ! tt
!  factor 0.5 from self-collision; factor 2.0 from 2 neutrons/ reaction
       if(ilev_ttnuk.eq.0) then
          call fusn_thermo_tt(ztikev,zsigv,1)  ! returns zsigv in m**3/sec
       elseif(ilev_ttnuk.eq.1) then
          call fusn_thermo_dpctt(ztikev,zsigv,1)  ! returns zsigv in m**3/sec
       endif
       zrdens(5)=zdens(2)*zdens(2)*zsigv ! #/m3/sec
       if(n_thrf.ne.1) then
          if(zdens_rf(2).gt.zero) then
             if(ilev_ttnuk.eq.0) then
                call fusn_btsigv_tt(ztikev*0.001_fp,zti_rf(2)*0.001_fp,zsigv,1)
                zrdens(5)=zrdens(5)+2.0_fp*zdens(2)*zdens_rf(2)*zsigv ! #/m3/sec
                call fusn_btsigv_tt(zti_rf(2)*0.001_fp,zti_rf(2)*0.001_fp,zsigv,1)
                zrdens(5)=zrdens(5)+zdens_rf(2)*zdens_rf(2)*zsigv ! #/m3/sec
             elseif(ilev_ttnuk.eq.1) then
                call fusn_btsigv_dpctt(ztikev*0.001_fp,zti_rf(2)*0.001_fp,zsigv,1)
                zrdens(5)=zrdens(5)+2.0_fp*zdens(2)*zdens_rf(2)*zsigv ! #/m3/sec
                call fusn_btsigv_dpctt(zti_rf(2)*0.001_fp,zti_rf(2)*0.001_fp,zsigv,1)
                zrdens(5)=zrdens(5)+zdens_rf(2)*zdens_rf(2)*zsigv ! #/m3/sec
             endif
          endif
       endif
       zrdens(5)=zrdens(5)*1.0e6_fp ! #/cm3/sec
       ! 2 neutrons / reaction * 0.5 (self collisions) = 1.0, so no correction
    endif
 
    zrdens(6)=zero  ! no t-he3 data in preact yet.
 
  end subroutine nbi_thfusn
 
  !-----------------------------------------------------------------------
 
  subroutine nbi_1d_grid_fetch(inzp1,zxizb,ierr)

    use nbi_com
    implicit none

    integer, intent(in) :: inzp1        ! #zones + 1 (expected)
    real(fp) :: zxizb(inzp1)              ! xi grid bdys (returned)
    integer, intent(out) :: ierr        ! 0=ok; 1=inzp1.ne.nzp1

    call nbi_1d_grid

    if(inzp1.ne.nzp1) then
       write(nonlin,*) ' ?nbi_1d_grid_fetch: inzp1.ne.nzp1: ',inzp1,nzp1
       ierr=1
    else
       ierr=0
       zxizb(1:inzp1) = xizb(lcentr:lep1)
    endif

  end subroutine nbi_1d_grid_fetch

  subroutine nbi_1d_grid
 
    use nbi_com
  
    !
    !----------------------------------------------!
    !set the standard 1-d grid                     !
    !..............................................!
 
    implicit none
    integer :: j, iz, im1, ith
 
    !-----------------------------------
    ! define fixed mesh array-- time invariant zone coordinate "xi"
    !
    lcentr=2
    ledge=lcentr+nzones-1
    !
    lcm1=lcentr-1
    !
    lcp1=lcentr+1
    lep1=ledge+1
    !
    mzones=lep1
    !
    nzp1=nzones+1
    !-----------------------------------
    !normalized flux label
    !
    do j=1,mj
       xizb_r8(j)=(j-lcentr)*one/nzones
       xizb(j)=xizb_r8(j)
       xi_r8(j,1)=xizb_r8(j)
       xi(j,1)=xizb(j)
       xizc_r8(j)=(j-lcentr+half)/nzones
       xizc(j)=xizc_r8(j)
       xi_r8(j,2)=xizc_r8(j)
       xi(j,2)=xizc(j)
    enddo
    !-----------------------------------
    ! differential elements of "xi"
    !
    do j=2,mzones
       im1=j-1
       dxi(im1,2)=xi(j,1)-xi(im1,1)
       dxi(j,1)=xi(j,2)-xi(im1,2)
    enddo
    dxi(mzones,2)   = dxi(ledge,2)
    dxizb(1:mzones) = dxi(1:mzones, 1)
    dxizc(1:mzones) = dxi(1:mzones, 2)
 
  end subroutine nbi_1d_grid
 
  !--------------------------------------------------------------
 
  subroutine nbi_profcomck
    use nbi_com
    implicit none
    integer j
 
    ! adjustments-- after all profiles and updates have been read in...
 
    xiliml=zero
    xibado=half
    if(xdepmod.ne.one) then
       ndepmod=1
    else
       ndepmod=0
    endif
 
    nlvphi=.false.
    if(nmdifb.le.3) then ! nmdifb from plasma-state file
       nmdifb=0
    endif
 
   ! if(ndifbe_2d.gt.1) then
   !    nmdifb=4  ! diffusivity vs. (x,e) data is present
   ! endif

    do j=lcentr,lep1
       if(omegag(j).ne.zero) nlvphi=.true.
       if(ndifbe_2d.le.1) then
          ! diffusivity & radial transport vs. (x) only...
          ! (namelist controlled energy variation possible...)
          if(difb(j).ne.zero) nmdifb=1
          if(velb(j).ne.zero) nmdifb=1
          if(difb_bi(j).ne.zero) nmdifb=2
          if(velb_bi(j).ne.zero) nmdifb=2
          if(difb_fi(j).ne.zero) nmdifb=2
          if(velb_fi(j).ne.zero) nmdifb=2
       endif
    enddo
 
  end subroutine nbi_profcomck
