subroutine nubeam_frantic(ss_plasma,ilun_msgs,iwall,ihalo,ireco, &
     wrkpath,ierr)

  ! use frantic subroutine to evaluate only halo neutral transport as indicated, using
  ! input data from two state objects:
  !   ss_plasma: mhd equilibrium, thermal plasma temperatures & densities, etc.
  !
  ! output data: same as frantic output, with same units as frantic data
  ! this routine based on nubeam_comp_sub/ps_nbi_frantic.f90
  ! 
  ! also, plasma_state will be loaded with output data for:
  !
  ! n0_halo(~nrho_nbi,nspec_gas) --  halo thermal neutral density, m^-3
  ! t0_halo(~nrho_nbi,nspec_gas) --  halo neutral temperature (2/3)<e0>, kev
  ! omeg0_halo(~nrho_nbi,nspec_gas)--  halo neutral toroidal angular velocity, rad/sec
  ! sbehalo(~nrho_nbi) --- cold electron sce due to
                         ! recapture of beam halo charge exchange thermal neutrals
                         ! = sum z(i)*sb0halo_recap(i), if electrons do not go into
                         ! an impurity charge state spectrum, #/sec 
  ! sb0halo_recap(~nrho_nbi,nspec_gas) --  cx neutral recapture
                                         ! recapture of thermal neutrals -- ion source, #/sec 
  ! pb0_halo(~nrho_nbi) --  power in halo neutral source (loss from thermal ions), w/zone
  ! psc_halo(~nrho_nbi) --  beam halo power: recapture - pb0_halo
                           !  (net source/sink to thermal ions)
                           !  (associated particle source/sink exists), w/zone
  ! pcx_halo(~nrho_nbi) --  beam halo driven cx power
                            !  (>0 means local source)
                            !  (no associated particle source/sink), w/zone
  ! tqb0_halo(~nrho_nbi) --  torque of halo neutral source (loss from thermal ions), nt*m/zone
  ! tqsc_halo(~nrho_nbi) --  beam halo torque: recapture - tqb0_halo
                            !  (net source/sink to thermal ions), nt*m/zone
                            !  (associated particle source/sink exists)
  ! tqcx_halo(~nrho_nbi) --  beam halo driven cx momentum torque
                            !  (positive value = momentum source)
                            !  (no associated particle source/sink), nt*m/zone

  ! sfsce(~nrho_fus,0:nspec_tha) -- net thermalization source of fusion ions
                                   ! sources due to fusion ions:  fusion ion thermalization;
                                   ! corrections due to charge exchange effects (usually small), #/sec
  ! sbsce(~nrho_nbi,0:nspec_tha) --  net source, thermal ions & electrons, #/sec


  use plasma_state_mod
  use ezspline_obj,only : ezspline1
  use ezspline,only : ezspline_init,ezspline_free,ezspline_setup
  use nbi_com
  use physconst_mod, only : z1000th=>zev2kev,zncv=>zm3tocm3,zvcv=>zcmtom
  implicit none

  type (plasma_state) :: ss_plasma ! state containing plasma & equilibrium
  ! state to contains gas nbi & fus component

  integer, intent(in) :: ilun_msgs

  logical, intent(in) :: iwall     ! (re)compute wall source neutral transport
  logical, intent(in) :: ihalo     ! (re)compute halo source neutral transport
  logical, intent(in) :: ireco     ! (re)compute recombination source transport

  character*(*), intent(in) :: wrkpath  ! directory for file output
  ! (only used if idebug=.true., below).

  integer, intent(out) :: ierr     ! status code returned, 0=ok

  !---------------------------------------
  integer :: ing,inx,intot,infast,inzons,iwarn
  integer :: iz,ia,iz0,ia0,ii,jj,ig,igp,igi,isc,ix,imatch,iok
  integer :: izf,iaf,izth,iath
  integer :: ik,ii_m
  real(fp) :: znorm0,znorm1,zcoeff
  integer, dimension(:), allocatable :: imapg,imapx,i0map
  real(fp) , dimension(:), allocatable :: zaa,zzz,zzfast
  real(fp) , dimension(:), allocatable :: zx,zxzcx,zrmj,zf,zvol,zdvol
  real(fp) , dimension(:), allocatable :: zrmjzc,zfzc,zsumzc
  real(fp) , dimension(:), allocatable :: zte,zresult,zresult2
  real(fp) , dimension(:), allocatable :: zti,zomeg
  real(fp) :: zflux,zengy
  type(ezspline1) :: zspline1,zspline2,zspline3

  logical :: change_grid

  logical, parameter :: idebug = .true.  ! set .true. for frantic test output
  character*180 :: debug_filename

  integer :: ind_l,ind_u,inx_ch
  !---------------------------------------

  ierr = 0

  call frantic_free
  ! first zero out quantities to be computed; 
  call nubeam_frantic_zero

  if(.not.(iwall.or.ihalo.or.ireco)) return

  if(iwall.or.ireco) then
     write(nonlin,*) ' ? nubeam_frantic: not support wall and/or recomb. source transport'
     return ! not support
  endif

 
  !-------------------
  ! #zones

  inzons = ss_plasma%nrho_nbi - 1  
  if(inzons.le.0) then
     ierr=1
     write(ilun_msgs,*) ' nubeam/frantic call attempt failed: #zones = ',inzons
     return
  endif

  !-------------------
  ! #neutral species

  ing = ss_plasma%nspec_gas
  if(ing.le.0) then
     ierr=1
     write(ilun_msgs,*) ' nubeam/frantic call attempt failed: #species = ',ing
     return
  endif

  !-------------------
  ! more sanity checks

  if(ss_plasma%nrho_eq.le.0) then
     ierr=1
     write(ilun_msgs,*) ' nubeam/frantic call attempt failed: null eq grid in ss_plasma.'
     return
  endif


  !-------------------
  ! #plasma species = #neutral species + #impurities
  !   in frantic, "impurities" refers to ion species for which no 
  !   neutral species are tracked.  note that the plasma state "abridged"
  !   species list is used: frantic does not deal with partially ionized
  !   impurity species and has only a simple model for impurity stopping
  !   of neutral species.  it is assumed that for each neutral species there
  !   is a corresponding (fully stripped light ion) plasma species present.
  !   frantic can treat isotopes of h, he, and li as "non-impurity" neutral
  !   and plasma species, but, the stopping cross sections for li are suspect.

  inx_ch = ss_plasma%nspec_th - ing ! number of impurity chages used in code
  if(nmimp.eq.0) then ! using single tokamakium impurity
     inx=1
  else
     inx = ss_plasma%nspec_imp0 ! number of impurity elements used in code
     if(inx_ch.ne.nrhix) then
        ierr=1
        write(ilun_msgs,*) ' nubeam/frantic call attempt failed: wrong number of impurity chages'
        return
     endif
  endif
  intot = ing + inx  ! may not equals to ss_plasma%nspec_tha
  allocate(imapg(intot),imapx(intot),zaa(intot),zzz(intot))
  allocate(i0map(ing))

  !-------------------
  ! #fast species

  infast = ss_plasma%nspec_beam + ss_plasma%nspec_fusion + &
       ss_plasma%nspec_rfmin
  allocate(zzfast(max(1,infast))); zzfast(1)=zero

  !-------------------
  ! set up species list for frantic

  imapg=0; imapx=0
  zaa(1:ing) = ss_plasma%m_sgas(1:ing)/ps_mp
  zzz(1:ing) = ss_plasma%qatom_sgas(1:ing)/ps_xe

  ix = 0
  do ii=1,intot
     if(ii.le.ing) then
        ii_m=ss_plasma%sgas_to_s(ii) ! map sgas index to "s" list
     else
! for impurity
        ii_m=ss_plasma%simp0_to_s(ii-ing) ! map simp0 index to "s" list
     endif
     iz = ss_plasma%qatom_s(ii_m)/ps_xe + 0.1_fp
     ia = ss_plasma%m_s(ii_m)/ps_mp + 0.1_fp

     imatch = 0
     do jj=1,ing
        iz0=zzz(jj)+0.1_fp
        ia0=zaa(jj)+0.1_fp
        if((ia0.eq.ia).and.(iz0.eq.iz)) then
           imatch=jj
           exit
        endif
     enddo

     if(imatch.gt.0) then
        imapg(ii)=imatch

     else
        ix=ix+1
        if(nmimp.eq.0) then !using single tokamakium impurity
           zaa(ing+ix)=aimp
           zzz(ing+ix)=xzimp
        else
           zaa(ing+ix)=ss_plasma%m_s(ii_m)/ps_mp
           zzz(ing+ix)=ss_plasma%qatom_s(ii_m)/ps_xe
        endif
        imapx(ii)=ing+ix !# impurities

     endif
  enddo

  ix = 0
  do ii = 1,ss_plasma%nspec_beam
     ix = ix + 1
     zzfast(ix) = ss_plasma%qatom_snbi(ii)/ps_xe
  enddo

  do ii = 1,ss_plasma%nspec_fusion
     ix = ix + 1
     zzfast(ix) = ss_plasma%qatom_sfus(ii)/ps_xe
  enddo

  do ii = 1,ss_plasma%nspec_rfmin
     ix = ix + 1
     zzfast(ix) = ss_plasma%qatom_rfmin(ii)/ps_xe
  enddo
  zaa=zaa

  !------------------------------------
  ! frantic calls...

  call frantic_lun_set(ilun_msgs)

  call frantic_species_za(ing,inx,zzz,zaa,infast,zzfast)

  !------------------------------------
  ! radial grid; #zones

  call alloc_prof(ss_plasma%rho_nbi)
  if(ierr.ne.0) return

  !------------------------------------
  ! plasma parameters

  call init_prof
  if(ierr.ne.0) return
  call isplfree(zspline1)
  call isplfree(zspline2)
  call isplfree(zspline3)

  !-----------------------------
  ! compute wall sources -- will be renormalized...

  znorm0 = 1.0e+14_fp  ! temporary normalization, nominally ptcls/m^2 influx
  ! an additional normalization will be applied after the frantic computation
  ! to produce state profiles in the desired units, i.e. per unit atom influx.

  ig = 0

  !! dmc: this call not needed
  !!   (only important for wall, not halo sources)
  !!
  !!      call ps_gsc0_species_map(ierr,state=ss_plasma)
  !!      if(ierr.ne.0) return

  call clear_nubeam_frantic
  call isplinit1(zspline1)
  call frantic_set_atomic(lev_nbidep)
    
  do ig=1,ing
     if(maxval(abs(ss_plasma%sb0halo(:,ig))).le.zero) cycle
     
     call frantic_neutral_specie(zzz(ig),zaa(ig),i0map)
     call load_bsink
     
     call fmap1('halo_sce',ss_plasma%sb0halo(:,ig)/zdvol)
     call frantic_volume_source(zero,zspline1)
     
     if(idebug) then
        debug_filename = trim(wrkpath)//'nubeam_frantic_src_halo_'// &
             trim(ss_plasma%sgas_name(ig))//'.cdf'
        call frantic_input_save(debug_filename,iwarn)
        if(iwarn.ne.0) write(ilun_msgs,*) ' %nubeam_frantic: '// &
             'debug file write failed: '//trim(debug_filename)
     endif
     !----------------------------
     ! perform frantic computation
     call frantic_neuts(ierr)
     if(ierr.ne.0) then
        write(ilun_msgs,*) ' ?? nubeam_frantic(halo): frantic status: ',ierr
        exit
     endif
     
     !----------------------------
     ! get results; each source (ig) leads to contributions to other species
     ! due to charge exchange; these are added together here...
     
     do igi=1,ing
        ik=i0map(igi)
        if(ik.gt.0) then
           call frantic_get_n0(ik,zresult)
           nbi_den0mn(igi,lcentr:ledge,ig) = zresult(1:inzons) !m**-3
           ss_plasma%n0_halo(:,igi) = ss_plasma%n0_halo(:,igi) + zresult ! #/m^3
           
           call frantic_get_t0(ik,zresult2)
           nbi_en0mn(igi,lcentr:ledge,ig) = zresult2(1:inzons)  ! kev
           ss_plasma%t0_halo(:,igi) = ss_plasma%t0_halo(:,igi) + &
                zresult*zresult2  ! density weighted sum *kev
           
           call frantic_get_vphi0(ik,zresult2)
           nbi_omg0mn(igi,lcentr:ledge,ig) = zresult2(1:inzons)  ! m/sec
           ss_plasma%omeg0_halo(:,igi) = ss_plasma%omeg0_halo(:,igi) + &
                zresult*zresult2/zrmjzc  ! density weighted sum *(rad/s)

           call frantic_get_si(ik,zresult)
           nbi_snimn(igi,lcentr:ledge,ig) = zresult(1:inzons)    !m**-3/sec
           ss_plasma%sb0halo_recap(:,igi) = ss_plasma%sb0halo_recap(:,igi) + &
                zresult*zdvol     ! #/zone/s
           
           ss_plasma%sbehalo = ss_plasma%sbehalo + zresult*zdvol*zzz(igi)
           
           call frantic_get_outflux(ik,zflux,zengy)
           nbi_outflx(igi,ig)=zflux
           nbi_eescav(igi,ig)=zengy
           nbi_cexflx(igi,ig)=nbi_outflx(igi,ig)
           
        endif
     enddo
     
     call frantic_get_vsor_en(zresult)
     nbi_pb0mn(lcentr:ledge,ig)=nbi_pb0mn(lcentr:ledge,ig)+zresult(1:inzons)
     nbi_pscmn(lcentr:ledge,ig)=nbi_pscmn(lcentr:ledge,ig)-zresult(1:inzons)

     ss_plasma%pb0_halo = ss_plasma%pb0_halo + zresult*zdvol  ! w/zone
     ss_plasma%psc_halo = ss_plasma%psc_halo - zresult*zdvol


     
     call frantic_get_vsor_mo(zresult)
     nbi_tqb0mn(lcentr:ledge,ig)=nbi_tqb0mn(lcentr:ledge,ig)+zresult(1:inzons)
     nbi_tqscmn(lcentr:ledge,ig)=nbi_tqscmn(lcentr:ledge,ig)-zresult(1:inzons)
     
     ss_plasma%tqb0_halo = ss_plasma%tqb0_halo + zresult*zdvol  ! nt*m/zone
     ss_plasma%tqsc_halo = ss_plasma%tqsc_halo - zresult*zdvol
     
     call frantic_get_se(zresult)
     nbi_snemn(1,lcentr:ledge,ig)=zresult(1:inzons)
     
     call frantic_get_seb(zresult)
     nbi_snemn(2,lcentr:ledge,ig)=zresult(1:inzons)
     
     call frantic_get_scx(zresult)
     nbi_scxmn(lcentr:ledge,ig)=zresult(1:inzons)
     
     call frantic_get_pcx(zresult)
     nbi_pcxmn(lcentr:ledge,ig)=zresult(1:inzons)
     ! pcxsf -- frictional term due to velocity differences btw zones
     
     call frantic_get_pcxsf(zresult2)! (reduces cx loss power)
     nbi_pfcxmn(lcentr:ledge,ig)=zresult2(1:inzons)
     
     ss_plasma%pcx_halo = ss_plasma%pcx_halo + (zresult-zresult2)*zdvol  ! w/zone
     
     call frantic_get_tqcx(zresult)
     nbi_tqcxmn(lcentr:ledge,ig)=zresult(1:inzons)
     
     ss_plasma%tqcx_halo = ss_plasma%tqcx_halo + zresult*zdvol    ! => nt*m/zone 
     
     call frantic_get_tqii(zresult)
     nbi_tqnimn(lcentr:ledge,ig)=zresult(1:inzons)
     
     ss_plasma%tqsc_halo = ss_plasma%tqsc_halo + zresult*zdvol ! => nt*m/zone
     
     call frantic_get_pii(zresult)
     nbi_pnimn(lcentr:ledge,ig)=zresult(1:inzons)
     
     ! piisf -- frictional term due to velocity differences btw zones
     call frantic_get_piisf(zresult2) ! (increases ionization source power)
     nbi_pfnimn(lcentr:ledge,ig)=zresult2(1:inzons)
     
     ss_plasma%psc_halo = ss_plasma%psc_halo + &
          (zresult+zresult2)*zdvol    ! w/zone
     

     call frantic_get_sflx01(zresult)
     nbi_sflx0m(1,lcentr:ledge,ig)=zresult(1:inzons)
     
     call frantic_get_sflx02(zresult)
     nbi_sflx0m(2,lcentr:ledge,ig)=zresult(1:inzons)
     
     call frantic_get_pflx01(zresult)
     nbi_pflx0m(1,lcentr:ledge,ig)=zresult(1:inzons)
     
     call frantic_get_pflx02(zresult)
     nbi_pflx0m(2,lcentr:ledge,ig)=zresult(1:inzons)
     
     call frantic_get_tqflx0(zresult)
     nbi_tqflx0(lcentr:ledge,ig)=zresult(1:inzons)
     
     
  enddo

     ! divide out summed n0 to get avg t0 & omeg0 from all halo sources

  do ig=1,ing
     do jj=1,inzons
        if(ss_plasma%n0_halo(jj,ig).gt.zero) then
           ss_plasma%t0_halo(jj,ig) = &
                ss_plasma%t0_halo(jj,ig)/ss_plasma%n0_halo(jj,ig)
           ss_plasma%omeg0_halo(jj,ig) = &
                ss_plasma%omeg0_halo(jj,ig)/ss_plasma%n0_halo(jj,ig)
        endif
     enddo
  enddo
  
     ! finally: complete sbsce & sfsce elements
     ! all halo sources attributed to beam

     if(ss_plasma%nrho_fus.gt.0) then
        ss_plasma%sfsce = zero
        do ii = 1,ss_plasma%nspec_fusion
           izf = ss_plasma%qatom_sfus(ii)/ps_xe + 0.2_fp
           iaf = ss_plasma%m_sfus(ii)/ps_mp + 0.2_fp
           do jj=1,ss_plasma%nspec_tha
              izth = ss_plasma%qatom_sa(ii)/ps_xe + 0.2_fp
              iath = ss_plasma%m_sa(ii)/ps_mp + 0.2_fp
              if((izth.eq.izf).and.(iath.eq.iaf)) then
                 ss_plasma%sfsce(:,jj) = ss_plasma%sftherm(:,ii)
                 exit
              endif
           enddo
        enddo
     endif

     ! since halo data is on the nbi grid, we know the nbi gridded arrays
     ! are available if we reach this code.

     ss_plasma%sbsce(:,0) = ss_plasma%sbedep + ss_plasma%sbehalo

     do ii=1,ss_plasma%nspec_tha
        ss_plasma%sbsce(:,ii) = zero
     enddo

     do ii=1,ss_plasma%nspec_beam
        izf = ss_plasma%qatom_snbi(ii)/ps_xe + 0.2_fp
        iaf = ss_plasma%m_snbi(ii)/ps_mp + 0.2_fp
        do jj=1,ss_plasma%nspec_tha
           izth = ss_plasma%qatom_sa(ii)/ps_xe + 0.2_fp
           iath = ss_plasma%m_sa(ii)/ps_mp + 0.2_fp
           if((izth.eq.izf).and.(iath.eq.iaf)) then
              ss_plasma%sbsce(:,jj) = ss_plasma%sbsce(:,jj) + ss_plasma%sbtherm(:,ii)
              exit
           endif
        enddo
     enddo

     do ii=1,ss_plasma%nspec_gas
        jj = ss_plasma%sgas_to_s(ii)
        jj = ss_plasma%s_to_sa(jj)
        ss_plasma%sbsce(:,jj) = ss_plasma%sbsce(:,jj) + &
             ss_plasma%sb0halo_recap(:,ii) - ss_plasma%sb0halo(:,ii)
     enddo


  call isplfree(zspline1)
  !-----------------------------
  ! cleanup

  deallocate(i0map)
  deallocate(imapg,imapx)
  deallocate(zaa,zzz,zzfast)

  call free_loc_prof

  call frantic_free

contains

  subroutine clear_nubeam_frantic
    nbi_den0mn = zero
    nbi_en0mn = zero
    nbi_omg0mn = zero
    nbi_snimn = zero
    nbi_outflx = zero
    nbi_eescav = zero
    nbi_cexflx = zero
    nbi_snemn = zero
    nbi_scxmn = zero
    nbi_pcxmn = zero
    nbi_pnimn = zero
    nbi_tqcxmn = zero
    nbi_tqnimn = zero
    nbi_pfcxmn = zero
    nbi_pfnimn = zero
    nbi_sflx0m = zero
    nbi_pflx0m = zero
    nbi_tqflx0 = zero
     !  zero out ps quantities to be computed; 
     !  nubeam has provided ss_nbi%sb0halo(<rho_nbi>,1:ing)

     ss_plasma%sb0halo_recap = zero
     ss_plasma%sbehalo = zero

     ss_plasma%pb0_halo = zero
     ss_plasma%psc_halo = zero
     ss_plasma%pcx_halo = zero

     ss_plasma%tqb0_halo = zero
     ss_plasma%tqsc_halo = zero
     ss_plasma%tqcx_halo = zero

     ss_plasma%n0_halo = zero
     ss_plasma%t0_halo = zero
     ss_plasma%omeg0_halo = zero
    

  end subroutine clear_nubeam_frantic

  subroutine alloc_prof(zxb)

    real(fp), dimension(:) :: zxb

    call frantic_profile_alloc(inzons)

    allocate(zx(inzons+1),zf(inzons+1),zrmj(inzons+1),zvol(inzons+1),zxzcx(inzons+2))
    allocate(zrmjzc(inzons),zfzc(inzons),zsumzc(inzons),zresult(inzons),zresult2(inzons))
    allocate(zomeg(inzons))
    allocate(zte(inzons),zti(inzons),zdvol(inzons))

    zx = zxb

    zxzcx(1)=zx(1)
    zxzcx(inzons+2)=zx(inzons+1)
    zxzcx(2:inzons+1)=(zx(1:inzons)+zx(2:inzons+1))*0.5_fp

    call ps_intrp_1d(zx,ss_plasma%id_vol,zvol,ierr, state=ss_plasma)
    call ckintrp('vol')
    if(ierr.ne.0) return

    zdvol = zvol(2:inzons+1)-zvol(1:inzons)  ! zone volumes
    zrmj(1:inzons+1) =0.01_fp* rmjrmp(lcentr:lep1,1)
    zrmjzc(1:inzons) = 0.01_fp*rmjrmp(lcentr:lcentr+inzons-1,2)
    ! interpolation of (thermal species average) ti and omega
    zte(1:inzons)=z1000th*te(lcentr:ledge,2)
    zti(1:inzons)=z1000th*ti(lcentr:lcentr+inzons-1,2)
    zomeg(1:inzons)=zvcv*rmjrmp(lcentr:ledge,2)*0.5_fp*(omegag(lcentr:ledge)+omegag(lcentr+1:ledge+1))


    ! set zone volumes and major radii

    call frantic_radius(1,inzons+1,zx,zvol,zrmj)

  end subroutine alloc_prof

  subroutine init_prof


    call isplinit(zspline1)
    call isplinit(zspline2)
    call isplinit(zspline3)

    ix=0
    ind_l=1
    do ii=0,intot
       if(ii.eq.0) then
          call fmap('electrons',zncv*rhoel(lcentr:ledge,2), &
               &  zte(1:inzons),zomeg(1:inzons))

          call frantic_profile_set(0,zspline1,zspline2,zspline3)

       else if(imapg(ii).gt.0) then
          call fmap('thermal gas',zncv*rhob(lcentr:ledge,ii,2), &
               &  zti(1:inzons),zomeg(1:inzons))
          
          call frantic_profile_set(imapg(ii),zspline1,zspline2,zspline3)

       else if(imapx(ii).gt.0) then
          ix=ix+1
! find low and upper indicies 
          if(ix.lt.inx) then
             ind_u=ind_l + ss_plasma%simp0_to_s(ix+1)-ss_plasma%simp0_to_s(ix)-1
          else
            ind_u=inx_ch
         endif
! rhix  is a multiple impurity density as a function of number of impurity charges used in code
! we need a multiple impurity density as a function of number of impurity element used in code
          call fmap('impurity',zncv*sum(rhix(lcentr:ledge,2,ind_l:ind_u),dim=2), &
               &  zti(1:inzons),zomeg(1:inzons))

          call frantic_profile_set(imapx(ii),zspline1,zspline2,zspline3)
          ind_l=ind_u + 1 

       endif
    enddo

    if(inx.eq.1) then
       !  impurity z and a profiles also
       
       call fmap('impurity z&a',&
            one*aimpj(lcentr:ledge),&
            one*xzimpj(lcentr:ledge),&
            one*xzimpj(lcentr:ledge))   ! zspline3 wasted.
       
       call frantic_impurity_zaprof(zspline2,zspline1,iok)
       if(iok.ne.0) call bad_exit  ! this would mean i got z & a order wrong.
    endif
    

  end subroutine init_prof

  subroutine free_loc_prof

    deallocate(zx,zf,zrmj,zvol,zxzcx,zrmjzc,zfzc,zsumzc,zresult,zresult2)
    deallocate(zdvol,zti,zomeg)
    deallocate(zte)

  end subroutine free_loc_prof

  subroutine ckintrp(str)
    ! error message
    character*(*), intent(in) :: str

    if(ierr.ne.0) then
       write(ilun_msgs,*) &
            ' ?nubeam_frantic: plasma state interpolation error: '//trim(str)
    endif

  end subroutine ckintrp

  subroutine isplinit(zspline)
    ! spline object init
    type(ezspline1) :: zspline
 
    call ezspline_init(zspline, inzons+2, (/1, 0/), iok)
    zspline%ishermite = 1
    zspline%bcval1min = 0        ! df/dx = 0 on axis
    zspline%x1 = zxzcx

  end subroutine isplinit

  subroutine isplinit1(zspline)
    ! spline object init
    type(ezspline1) :: zspline
 
    call ezspline_init(zspline, inzons+2, (/1, 1/), iok)
    zspline%ishermite = 1
    zspline%bcval1min = 0        ! df/dx = 0 on axis
    zspline%bcval1max = 0   ! df/dx = 0 at edge
    zspline%x1 = zxzcx

  end subroutine isplinit1


  subroutine isplfree(zspline)
    ! spline object free
    type(ezspline1) :: zspline
    call ezspline_free(zspline, iok)
    if(iok.ne.0) call errmsg_exit(' ?nubeam_frantic: ezspline_free error 1.')
    
  end subroutine isplfree

  subroutine fmap1(lbl,zf)
    ! set up 1 profile

    character*(*), intent(in) :: lbl
    real(fp), dimension(:), intent(in) :: zf

    real(fp) :: zbuf(inzons+2)

    call insert(zf,zbuf)
    call ezspline_setup(zspline1, zbuf, iok)
    if(iok.ne.0) write(ilun_msgs,*) &
         ' %nubeam_frantic: ezspline warning, fmap1, lbl=',trim(lbl)

  end subroutine fmap1

  subroutine fmap(lbl,zn,zt,zvpar)
    ! set up plasma specie splines

    character*(*), intent(in) :: lbl
    real(fp), dimension(:), intent(in) :: zn
    real(fp), dimension(:), intent(in) :: zt
    real(fp), dimension(:), intent(in) :: zvpar

    real(fp) :: zbuf(inzons+2)

    call insert(zn,zbuf)
    call ezspline_setup(zspline1, zbuf, iok)
    if(iok.ne.0) write(ilun_msgs,*) &
         ' %nubeam_frantic: ezspline warning on n, lbl=',trim(lbl)

    call insert(zt,zbuf)
    call ezspline_setup(zspline2, zbuf, iok)
    if(iok.ne.0) write(ilun_msgs,*) &
         ' %nubeam_frantic: ezspline warning on t, lbl=',trim(lbl)

    call insert(zvpar,zbuf)
    call ezspline_setup(zspline3, zbuf, iok)
    if(iok.ne.0) write(ilun_msgs,*) &
         ' %nubeam_frantic: ezspline warning on vpar, lbl=',trim(lbl)

  end subroutine fmap

  subroutine insert(zprof,zspl_data)
    ! insert profile with guarded extrapolation to axis & edge

    real(fp), dimension(:), intent(in) :: zprof
    real(fp), dimension(:), intent(out) :: zspl_data

    real(fp) :: zlim1,zlim2

    !----------------------------------
    zspl_data(2:inzons+1) = zprof(1:inzons)

    zlim1=min(0.8_fp*zprof(1),1.2_fp*zprof(1))
    zlim2=max(0.8_fp*zprof(1),1.2_fp*zprof(1))

    zspl_data(1) = max(zlim1, min(zlim2, (9*zprof(1)-zprof(2))/8 ))

    zlim1=min(0.8_fp*zprof(inzons),1.2_fp*zprof(inzons))
    zlim2=max(0.8_fp*zprof(inzons),1.2_fp*zprof(inzons))

    zspl_data(inzons+2) = max(zlim1, min(zlim2, &
         (3*zprof(inzons)-zprof(inzons-1))/2 ))

  end subroutine insert

  subroutine load_bsink

    ! load nb*<sigma*v> stopping profiles due to fast ions
    !   only for those neutral species that are present...

    do igi=1,ing
       ik=i0map(igi)
       if(ik.gt.0) then
          zsumzc = zero
          if(ss_plasma%nspec_beam.gt.0) then
             call ps_intrp_1d(zxzcx(2:inzons+1), &
                  ss_plasma%id_rate_sinb0i(igi),zfzc,ierr, state=ss_plasma)
             zsumzc=zsumzc+zfzc
          endif
          if(ss_plasma%nspec_fusion.gt.0) then
             call ps_intrp_1d(zxzcx(2:inzons+1), &
                  ss_plasma%id_rate_sinf0i(igi),zfzc,ierr, state=ss_plasma)
             zsumzc=zsumzc+zfzc
          endif
          call fmap1('sinb0i',zsumzc)
          call frantic_bsink_ii(ik,zspline1)
          ix=0
          do ii=1,ss_plasma%nspec_beam
             ix=ix+1
             call ps_intrp_1d(zxzcx(2:inzons+1), &
                  ss_plasma%id_rate_sinb0xs(igi,ii),zfzc,ierr, state=ss_plasma)
             call fmap1('sinb0xs',zfzc)
             call frantic_bsink_cx(ik,ix,zspline1)
          enddo
          do ii=1,ss_plasma%nspec_fusion
             ix=ix+1
             call ps_intrp_1d(zxzcx(2:inzons+1), &
                  ss_plasma%id_rate_sinf0xs(igi,ii),zfzc,ierr, state=ss_plasma)
             call fmap1('sinf0xs',zfzc)
             call frantic_bsink_cx(ik,ix,zspline1)
          enddo
       endif
    enddo  ! resultant neutral species loop

  end subroutine load_bsink
  
  subroutine nubeam_frantic_zero
    nbi_den0mn = zero
    nbi_en0mn = zero
    nbi_omg0mn = zero
    nbi_snimn = zero
    nbi_outflx = zero
    nbi_eescav = zero
    nbi_cexflx = zero
    nbi_snemn = zero
    nbi_scxmn = zero
    nbi_pcxmn = zero
    nbi_pnimn = zero
    nbi_tqcxmn = zero
    nbi_tqnimn = zero
    nbi_pfcxmn = zero
    nbi_pfnimn = zero
    nbi_sflx0m = zero
    nbi_pflx0m = zero
    nbi_tqflx0 = zero
    nbi_pb0mn = zero
    nbi_pscmn = zero
    nbi_tqb0mn = zero
    nbi_tqscmn = zero
  end subroutine nubeam_frantic_zero
end subroutine nubeam_frantic
