subroutine ps_frantic(ss_plasma,ss_sces,ilun_msgs,iwall,ihalo,ireco,ierr)

  ! use FRANTIC to evaluate neutral transport as indicated, using
  ! input data from two state objects:
  !
  !   ss_plasma: MHD equilibrium, thermal plasma temperatures & densities, etc.
  !
  !   ss_sces: NUBEAM output: halo neutral source, nb*<sigma*v> thermal neutral
  !            stopping rate coefficients (**if available **)

  ! output data: (a) profiles in state GAS component are updated in ss_sces;
  !              (b) (if ihalo .TRUE.) complete the partice source profiles
  !                  for the NBI component, in ss_sces.

  ! DMC -- Feb 2011: pad FRANTIC input species list with a manufactured
  ! impurity, if necessary

  use plasma_state_mod
  use ps_frantic_privmod  ! idebug & dbg_path defined here...
  use ezspline_obj
  use ezspline
  implicit NONE

  INTEGER, PARAMETER :: R8=SELECTED_REAL_KIND(12,100)
  type (plasma_state) :: ss_plasma ! state containing plasma & equilibrium
  type (plasma_state) :: ss_sces   ! state to contains GAS NBI & FUS component

  ! updates.  NOTE: ss_plasma & ss_sces could be the same object.  If they are
  ! not, it is assumed that the dimensioning of species lists are consistent.

  integer, intent(in) :: ilun_msgs

  integer, intent(in) :: iwall  ! .ne.0: compute wall source neutral transport
  integer, intent(in) :: ihalo  ! .ne.0: compute halo source neutral transport
  integer, intent(in) :: ireco  ! .ne.0: compute recombination source transport

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

  !---------------------------------------
  integer :: ing,inx,inxa,intot,intota,infast,inzons,iwarn
  integer :: iz,ia,iz0,ia0,ii,jj,ig,igp,igi,isc,ix,imatch,iok
  integer :: izf,iaf,izth,iath
  integer :: ik
  real*8 :: znorm0,znorm1,zcoeff
  integer, dimension(:), allocatable :: imapg,imapx,i0map
  real*8 , dimension(:), allocatable :: zaa,zzz,zzfast
  real*8 , dimension(:), allocatable :: zx,zxzcx,zrmj,zf,zvol,zdvol,zwk1,zwk2
  real*8 , dimension(:), allocatable :: zrmjzc,zfzc,zsumzc
  real*8 , dimension(:,:), allocatable :: zta,zna,zvpar
  real*8 , dimension(:), allocatable :: zte,zne,zvpare,zresult,zresult2
  real*8 , dimension(:), allocatable :: zr2,zmoi,zti,zomeg,zntot,zscx,zpcx,ztcx
  real*8 , dimension(:), allocatable :: zsigv

  logical :: add_impurity

  type(ezspline1) :: zspline1,zspline2,zspline3

  logical :: change_grid

  character*300 :: debug_filename

  real*8, parameter :: ZERO = 0.0d0
  real*8, parameter :: ONE  = 1.0d0
  integer :: ii_m
  !---------------------------------------

  ierr = 0

  call frantic_free

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

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

  inzons = ss_sces%nrho_gas - 1   ! grid may be in ss_sces only...
  if(inzons.le.0) then
     ierr=1
     write(ilun_msgs,*) ' 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,*) ' FRANTIC call attempt failed: #species = ',ing
     return
  endif

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

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

  if(minval(ss_plasma%gbr2).le.ZERO) then
     ierr=1
     write(ilun_msgs,*) ' FRANTIC call attempt failed: <B*R**2> .le. ZERO detected: equilibrium not ready?'
     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.

  intot = ss_plasma%nspec_tha ! wiil include impurities be charges
  inx  = intot - ing  
  intota = intot
  inxa = inx

  if(inx.eq.0) then
     inxa=1
     intota=intota+1
     add_impurity=.TRUE.
  else
     add_impurity=.FALSE.
  endif

  allocate(imapg(intota),imapx(intota),zaa(intota),zzz(intota))
  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
!MG bux fixed 09/27/2013
!        ii_m=ss_plasma%s_to_sa(ii) ! map simp0 index to "s" list
        ii_m=ss_plasma%simp0_to_s(ii-ing) ! map simp0 index to "s" list
     endif
!MG bux fixed 09/27/2013
!     iz = ss_plasma%qatom_sa(ii_m)/ps_xe + 0.1d0
     iz = ss_plasma%qatom_s(ii_m)/ps_xe + 0.1d0
     ia = ss_plasma%m_s(ii_m)/ps_mp + 0.1d0

     imatch = 0
     do jj=1,ing
        iz0=zzz(jj)+0.1d0
        ia0=zaa(jj)+0.1d0
        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
        zaa(ing+ix)=ss_plasma%m_s(ii_m)/ps_mp
!MG bux fixed 09/27/2013
!        zzz(ing+ix)=ss_plasma%qatom_sa(ii_m)/ps_xe
        zzz(ing+ix)=ss_plasma%qatom_s(ii_m)/ps_xe
        imapx(ii)=ing+ix
     endif
  enddo

  if(add_impurity) then 
     zaa(intota)= 12.0d0
     zzz(intota)=  6.0d0
     imapx(intota)= intota
  endif

  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

  !------------------------------------
  ! frantic calls...
  zaa=zaa

  call frantic_lun_set(ilun_msgs)

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

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

  call alloc_prof(ss_sces%rho_gas)
  if(ierr.ne.0) return

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

  call init_prof
  if(ierr.ne.0) return

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

  znorm0 = 1.0d14  ! 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

  call ps_gsc0_species_map(ierr,state=ss_sces)
  if(ierr.ne.0) return

  if((iwall.ne.0)) then
     do isc=1,ss_sces%ngsc0

        ! loop over edge neutral sources
        ! when the source species changes, update the fast ion nb*<sigma*v>
        ! stopping profile sets

        igp=ig
        ig=ss_sces%sc0_to_sgas(isc)

        if(ig.ne.igp) then
           call frantic_neutral_specie(zzz(ig),zaa(ig),i0map)
           call load_bsink
        endif  ! change in source neutral specie

        call frantic_edge_source(znorm0, ss_sces%e0_av(isc))
        call frantic_edge_vphi(ss_sces%vphi0_av(isc))
        
        if(idebug) then
           debug_filename = trim(dbg_path)//'frantic_src_'// &
                trim(ss_sces%gs_name(isc))//'.cdf'
           call frantic_input_save(debug_filename,iwarn)
           if(iwarn.ne.0) write(ilun_msgs,*) ' %ps_nbi_frantic: '// &
                'debug file write failed: '//trim(debug_filename)
        endif

        !----------------------------
        ! perform FRANTIC computation

        call frantic_neuts(ierr)
        if(ierr.ne.0) then
           write(ilun_msgs,*) ' ?? ps_nbi_frantic(wall): FRANTIC status: ',ierr
           exit
        endif

        ! output normalization: incident #atoms/sec @znorm0
        call frantic_get_inps(znorm1)
        
        ! fetch outputs

        do igi=1,ing
           ik=i0map(igi)
           if(ik.eq.0) then
              ss_sces%n0norm(:,igi,isc) = ZERO
              ss_sces%t0sc0(:,igi,isc) = ZERO
              ss_sces%omeg0sc0(:,igi,isc) = ZERO
              ss_sces%sprof0(:,igi,isc) = ZERO
           else
              call frantic_get_n0(ik,zresult)
              ss_sces%n0norm(:,igi,isc) = zresult/znorm1 ! #/m^3/influx

              call frantic_get_T0(ik,zresult)
              ss_sces%t0sc0(:,igi,isc) = zresult  ! KeV

              call frantic_get_si(ik,zresult)
              ss_sces%sprof0(:,igi,isc) = zresult*zdvol/znorm1 ! #/zone/s/influx

              call frantic_get_vphi0(ik,zresult)
              ss_sces%omeg0sc0(:,igi,isc) = zresult/zrmjzc  ! rad/s
           endif
        enddo

        call frantic_get_se(zresult)
        ss_sces%sprof0e(:,isc) = zresult*zdvol/znorm1  ! #/zone/s/influx

        ! CX data -- power and torque profiles

        call frantic_get_scx(zresult)
        zscx = zresult*zdvol/znorm1 ! CX reaction rate per zone per unit influx
        !  (used below for CX power and torque coefficients)

        !   pcxsf -- rotation associated frictional term

        call frantic_get_pcx(zresult)
        call frantic_get_pcxsf(zresult2) ! (reduces CX loss)
        zpcx = (zresult-zresult2)*zdvol/znorm1 ! W/zone/influx
        
        call frantic_get_tqcx(zresult)
        ztcx = zresult*zdvol/znorm1 ! Nt*m/zone/influx

        do jj=1,inzons
           if(zscx(jj).le.ZERO) then
              ss_sces%T0cx(jj,isc) = ZERO
              ss_sces%qqcx(jj,isc) = ZERO
              ss_sces%omeg0cx(jj,isc) = ZERO
              ss_sces%tqqcx(jj,isc) = ZERO
           else
              zcoeff = 1.5d3*ps_xe*zscx(jj)  ! 1e3*ps_xe = Jles/KeV; <E>=1.5*T
              ss_sces%qqcx(jj,isc) = zcoeff
              ss_sces%T0cx(jj,isc) = zti(jj)-zpcx(jj)/zcoeff  ! KeV
              zcoeff = zmoi(jj)*zscx(jj)
              ss_sces%tqqcx(jj,isc) = zcoeff
              ss_sces%omeg0cx(jj,isc) = zomeg(jj) - ztcx(jj)/zcoeff  ! rad/s
           endif
        enddo

        ! neutral ionization -- power and torque profiles
        !   piisf -- rotation associated frictional term

        call frantic_get_pii(zresult)
        call frantic_get_piisf(zresult2) ! (increases ionization source)
        ss_sces%qioniz(:,isc) = (zresult+zresult2)*zdvol/znorm1  ! W/zone/influx
        
        call frantic_get_tqii(zresult)
        ss_sces%tqioniz(:,isc) = zresult*zdvol/znorm1  ! Nt*m/zone/influx
        !----------------------------
     enddo  ! edge neutral source loop

  endif  ! ((iwall.ne.0))
  if(ierr.ne.0) return

  !-----------------------------
  ! compute volume sources: beam halo; recombination

  if((ireco.ne.0)) then
     ! first zero out PS quantities
     ss_sces%s0reco = ZERO
     ss_sces%s0reco_recap = ZERO
     ss_sces%s0reco_e = ZERO

     ss_sces%p0_reco = ZERO
     ss_sces%psc_reco = ZERO
     ss_sces%pcx_reco = ZERO

     ss_sces%tq0_reco = ZERO
     ss_sces%tqsc_reco = ZERO
     ss_sces%tqcx_reco = ZERO

     ss_sces%n0_reco = ZERO
     ss_sces%T0_reco = ZERO
     ss_sces%omeg0_reco = ZERO

     ! get recombination <sigma*v>(Te) -- isotope of H+ & e- -> H0
     ! at present no recombination computed for He...

     allocate(zsigv(inzons))
     call sigvte_reco(zte,inzons,zsigv)  ! zsigv in m^3/s

     do ig=1,ing
        if(abs(zzz(ig)-ONE).gt.0.1d0) cycle  ! non-H isotopes ignored...

        ss_sces%s0reco(:,ig) = zsigv*zne*zna(:,ig)*zdvol  ! #/zone/s
        if(maxval(abs(ss_sces%s0reco(:,ig))).le.ZERO) cycle

        ss_sces%s0reco_e = ss_sces%s0reco_e - ss_sces%s0reco(:,ig)

        call frantic_neutral_specie(zzz(ig),zaa(ig),i0map)
        call load_bsink

        call fmap1('reco_sce',ss_sces%s0reco(:,ig)/zdvol)
        call frantic_volume_source(ZERO,zspline1)

        if(idebug) then
           debug_filename = trim(dbg_path)//'frantic_src_reco_'// &
                trim(ss_sces%sgas_name(ig))//'.cdf'
           call frantic_input_save(debug_filename,iwarn)
           if(iwarn.ne.0) write(ilun_msgs,*) ' %ps_nbi_frantic: '// &
                'debug file write failed: '//trim(debug_filename)
        endif
        !----------------------------
        ! perform FRANTIC computation

        call frantic_neuts(ierr)
        if(ierr.ne.0) then
           write(ilun_msgs,*) ' ?? ps_nbi_frantic(reco): 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
              ! H isotope that CX's with source species (ig)
              call frantic_get_n0(ik,zresult)
              ss_sces%n0_reco(:,igi) = ss_sces%n0_reco(:,igi) + zresult ! #/m^3

              call frantic_get_T0(ik,zresult2)
              ss_sces%T0_reco(:,igi) = ss_sces%T0_reco(:,igi) + &
                   zresult*zresult2  ! density weighted sum *KeV

              call frantic_get_vphi0(ik,zresult2)
              ss_sces%omeg0_reco(:,igi) = ss_sces%omeg0_reco(:,igi) + &
                   zresult*zresult2/zrmjzc  ! density weighted sum *(rad/s)

              call frantic_get_si(ik,zresult)
              ss_sces%s0reco_recap(:,igi) = ss_sces%s0reco_recap(:,igi) + &
                   zresult*zdvol     ! #/zone/s

              ss_sces%s0reco_e = ss_sces%s0reco_e + zresult*zdvol
           endif
        enddo

        call frantic_get_vsor_en(zresult)
        ss_sces%p0_reco = ss_sces%p0_reco + zresult*zdvol  ! W/zone
        ss_sces%psc_reco = ss_sces%psc_reco - zresult*zdvol

        call frantic_get_vsor_mo(zresult)
        ss_sces%tq0_reco = ss_sces%tq0_reco + zresult*zdvol  ! Nt*m/zone
        ss_sces%tqsc_reco = ss_sces%tqsc_reco - zresult*zdvol

        call frantic_get_pii(zresult)
        ! piisf -- frictional term due to velocity differences btw zones
        call frantic_get_piisf(zresult2) ! (increases ionization source power)
        ss_sces%psc_reco = ss_sces%psc_reco + &
             (zresult+zresult2)*zdvol    ! W/zone

        call frantic_get_tqii(zresult)   ! => Nt*m/zone
        ss_sces%tqsc_reco = ss_sces%tqsc_reco + zresult*zdvol

        call frantic_get_pcx(zresult)
        ! pcxsf -- frictional term due to velocity differences btw zones
        call frantic_get_pcxsf(zresult2) ! (reduces CX loss power)
        ss_sces%pcx_reco = ss_sces%pcx_reco + (zresult-zresult2)*zdvol  ! W/zone

        call frantic_get_tqcx(zresult)   ! => Nt*m/zone
        ss_sces%tqcx_reco = ss_sces%tqcx_reco + zresult*zdvol

     enddo

     ! divide out summed n0 to get avg T0 & omeg0 from all reco sources

     do ig=1,ing
        do jj=1,inzons
           if(ss_sces%n0_reco(jj,ig).gt.ZERO) then
              ss_sces%T0_reco(jj,ig) = &
                   ss_sces%T0_reco(jj,ig)/ss_sces%n0_reco(jj,ig)
              ss_sces%omeg0_reco(jj,ig) = &
                   ss_sces%omeg0_reco(jj,ig)/ss_sces%n0_reco(jj,ig)
           endif
        enddo
     enddo

     deallocate(zsigv)
     if(ierr.ne.0) return

  endif

  if((ihalo.ne.0) .and. (ss_sces%nrho_nbi .gt. 0)) then

     change_grid=.FALSE.
     if(ss_sces%nrho_nbi .ne. ss_sces%nrho_gas) change_grid=.TRUE.
     if(.not.change_grid) then
        if(maxval(abs(ss_sces%rho_nbi - ss_sces%rho_gas)).gt.ZERO) then
           change_grid=.TRUE.
        endif
     endif

     change_grid = .TRUE.  ! *** TEMPORARY debug ***

     if(change_grid) then

        ! new grid... delete old grid-linked arrays here and in FRANTIC
        call free_loc_prof
        call frantic_free_profiles

        inzons = ss_sces%nrho_nbi - 1

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

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

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

        call init_prof
        if(ierr.ne.0) return
        call ezspline_free(zspline1, iok)
     endif

     ! first zero out PS quantities to be computed; 
     !  NUBEAM has provided ss_sces%sb0halo(<rho_nbi>,1:ing)

     ss_sces%sb0halo_recap = ZERO
     ss_sces%sbehalo = ZERO

     ss_sces%pb0_halo = ZERO
     ss_sces%psc_halo = ZERO
     ss_sces%pcx_halo = ZERO

     ss_sces%tqb0_halo = ZERO
     ss_sces%tqsc_halo = ZERO
     ss_sces%tqcx_halo = ZERO

     ss_sces%n0_halo = ZERO
     ss_sces%T0_halo = ZERO
     ss_sces%omeg0_halo = ZERO

     call isplinit1(zspline1)

     do ig=1,ing
        if(maxval(abs(ss_sces%sb0halo(:,ig))).le.ZERO) cycle

        call frantic_neutral_specie(zzz(ig),zaa(ig),i0map)
        call load_bsink

        call fmap1('halo_sce',ss_sces%sb0halo(:,ig)/zdvol)
        call frantic_volume_source(ZERO,zspline1)

        if(idebug) then
           debug_filename = trim(dbg_path)//'frantic_src_halo_'// &
                trim(ss_sces%sgas_name(ig))//'.cdf'
           call frantic_input_save(debug_filename,iwarn)
           if(iwarn.ne.0) write(ilun_msgs,*) ' %ps_nbi_frantic: '// &
                'debug file write failed: '//trim(debug_filename)
        endif
        !----------------------------
        ! perform FRANTIC computation

        call frantic_neuts(ierr)
        if(ierr.ne.0) then
           write(ilun_msgs,*) ' ?? ps_nbi_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)
              ss_sces%n0_halo(:,igi) = ss_sces%n0_halo(:,igi) + zresult ! #/m^3

              call frantic_get_T0(ik,zresult2)
              ss_sces%T0_halo(:,igi) = ss_sces%T0_halo(:,igi) + &
                   zresult*zresult2  ! density weighted sum *KeV

              call frantic_get_vphi0(ik,zresult2)
              ss_sces%omeg0_halo(:,igi) = ss_sces%omeg0_halo(:,igi) + &
                   zresult*zresult2/zrmjzc  ! density weighted sum *(rad/s)

              call frantic_get_si(ik,zresult)
              ss_sces%sb0halo_recap(:,igi) = ss_sces%sb0halo_recap(:,igi) + &
                   zresult*zdvol     ! #/zone/s

              ss_sces%sbehalo = ss_sces%sbehalo + zresult*zdvol*zzz(igi)
           endif
        enddo

        call frantic_get_vsor_en(zresult)
        ss_sces%pb0_halo = ss_sces%pb0_halo + zresult*zdvol  ! W/zone
        ss_sces%psc_halo = ss_sces%psc_halo - zresult*zdvol

        call frantic_get_vsor_mo(zresult)
        ss_sces%tqb0_halo = ss_sces%tqb0_halo + zresult*zdvol  ! Nt*m/zone
        ss_sces%tqsc_halo = ss_sces%tqsc_halo - zresult*zdvol

        call frantic_get_pii(zresult)
        ! piisf -- frictional term due to velocity differences btw zones
        call frantic_get_piisf(zresult2) ! (increases ionization source power)
        ss_sces%psc_halo = ss_sces%psc_halo + &
             (zresult+zresult2)*zdvol    ! W/zone

        call frantic_get_tqii(zresult)   ! => Nt*m/zone
        ss_sces%tqsc_halo = ss_sces%tqsc_halo + zresult*zdvol

        call frantic_get_pcx(zresult)
        ! pcxsf -- frictional term due to velocity differences btw zones
        call frantic_get_pcxsf(zresult2) ! (reduces CX loss power)
        ss_sces%pcx_halo = ss_sces%pcx_halo + (zresult-zresult2)*zdvol  ! W/zone

        call frantic_get_tqcx(zresult)   ! => Nt*m/zone
        ss_sces%tqcx_halo = ss_sces%tqcx_halo + zresult*zdvol

     enddo

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

     do ig=1,ing
        do jj=1,inzons
           if(ss_sces%n0_halo(jj,ig).gt.ZERO) then
              ss_sces%T0_halo(jj,ig) = &
                   ss_sces%T0_halo(jj,ig)/ss_sces%n0_halo(jj,ig)
              ss_sces%omeg0_halo(jj,ig) = &
                   ss_sces%omeg0_halo(jj,ig)/ss_sces%n0_halo(jj,ig)
           endif
        enddo
     enddo
     

     ! finally: complete sbsce & sfsce elements
     ! all halo sources attributed to beam

     if(ss_sces%nrho_fus.gt.0) then
        ss_sces%sfsce = ZERO
        do ii = 1,ss_sces%nspec_fusion
           izf = ss_sces%qatom_sfus(ii)/ps_xe + 0.2d0
           iaf = ss_sces%m_sfus(ii)/ps_mp + 0.2d0
           do jj=1,ss_sces%nspec_tha
              izth = ss_sces%qatom_sa(ii)/ps_xe + 0.2d0
              iath = ss_sces%m_sa(ii)/ps_mp + 0.2d0
              if((izth.eq.izf).and.(iath.eq.iaf)) then
                 ss_sces%sfsce(:,jj) = ss_sces%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_sces%sbsce(:,0) = ss_sces%sbedep + ss_sces%sbehalo

     do ii=1,ss_sces%nspec_tha
        ss_sces%sbsce(:,ii) = ZERO
     enddo

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

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

     if(ierr.ne.0) return
  endif

  !-----------------------------
  ! cleanup

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

  call free_loc_prof

  call frantic_free

CONTAINS
  subroutine alloc_prof(zxb)

    real*8, dimension(:) :: zxb

    call frantic_profile_alloc(inzons)

    allocate(zx(inzons+1),zrmj(inzons+1),zvol(inzons+1),zxzcx(inzons+2))
    allocate(zrmjzc(inzons),zsumzc(inzons),zresult(inzons),zresult2(inzons))
    allocate(zwk1(inzons+1),zwk2(inzons+1),zf(inzons+1),zfzc(inzons))
    allocate(zr2(inzons),zmoi(inzons),zscx(inzons),zpcx(inzons),ztcx(inzons))
    allocate(zti(inzons),zomeg(inzons),zntot(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.5d0

    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

    call ps_intrp_1d(zx,ss_plasma%id_R_midp_in,zwk1,ierr, state=ss_plasma)
    call ckintrp('R_midp_in')
    if(ierr.ne.0) return

    call ps_intrp_1d(zx,ss_plasma%id_R_midp_out,zwk2,ierr, state=ss_plasma)
    call ckintrp('R_midp_out')
    if(ierr.ne.0) return

    ! interpolation of <R**2> for moment-of-inertia estimate for CX torque coeff.
    call ps_intrp_1d(zxzcx(2:inzons+1),ss_plasma%id_gr2,zr2,ierr, &
         state=ss_plasma)
    call ckintrp('<R**2>')
    if(ierr.ne.0) return

    zrmj = 0.5d0*(zwk1+zwk2)
    zrmjzc = 0.5d0*(zrmj(1:inzons)+zrmj(2:inzons+1))

    ! interpolation of (thermal species average) Ti and omega
    call ps_ti_rezone(zx,ierr, state=ss_plasma, Ti=zTi, omegat=zomeg)
    call ckintrp('Ti & omega')
    if(ierr.ne.0) return

    ! set zone volumes and major radii

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

  end subroutine alloc_prof

  subroutine init_prof

    zmoi = ZERO
    zntot = ZERO

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

    allocate(zta(inzons,intota),zna(inzons,intota),zvpar(inzons,intota))
    allocate(zte(inzons),zne(inzons),zvpare(inzons))

    call ps_tha_rezone(zx,ierr, state=ss_plasma, &
         ns=zna(:,1:intot), ts=zta(:,1:intot), v_pars=zvpar(:,1:intot))

    if(add_impurity) then
       zna(:,intota) = 0.0d0
       zta(:,intota) = zta(:,intot)
       zvpar(:,intota) = zvpar(:,intot)
    endif

    call ckintrp('plasma_parameters')
    if(ierr.ne.0) return

    call ps_rho_rezone(zx,ss_plasma%id_ns(0),zne,ierr, state=ss_plasma)
    call ckintrp('ne')
    if(ierr.ne.0) return

    call ps_rho_rezone(zx,ss_plasma%id_Ts(0),zTe,ierr, state=ss_plasma)
    call ckintrp('Te')
    if(ierr.ne.0) return

    call ps_rho_rezone(zx,ss_plasma%id_v_pars(0),zvpare,ierr, state=ss_plasma)
    call ckintrp('vpar_e')
    if(ierr.ne.0) return

    ! to relate avg parallel velocity & angular velocity treated as flux fcn:
    !  <vpll> = <(B/B_phi)*v_phi> = (1/g)*<R*B*v_phi> = (omega/g)*<B*R**2>
    ! use state element R|units=T*m^2|pclin  gbr2(nrho_eq_geo)     ! <|B|*R^2>

    call ps_intrp_1d(zx,ss_plasma%id_gbr2,zwk1,ierr, state=ss_plasma)
    call ckintrp('gbr2')
    if(ierr.ne.0) return

    call ps_intrp_1d(zx,ss_plasma%id_g_eq,zwk2,ierr, state=ss_plasma)
    call ckintrp('g_eq')
    if(ierr.ne.0) return

    zf = zwk2/zwk1  ! (R*B_phi)/<R**2*B> ~ 1/R
    zfzc = 0.5d0*(zf(1:inzons)+zf(2:inzons+1))
    !  FRANTIC uses toroidal velocity in form R*omega ~ zrmjzc*zfzc*zvpar

    do ii=0,intota
       if(ii.eq.0) then
          call fmap('electrons',zne,zte,zvpare)
          call frantic_profile_set(0,zspline1,zspline2,zspline3)
       else if(imapg(ii).gt.0) then
          call fmap('thermal gas',zna(:,ii),zta(:,ii),zvpar(:,ii))
          call frantic_profile_set(imapg(ii),zspline1,zspline2,zspline3)
          zntot = zntot + zna(:,ii)
          zmoi  = zmoi  + zna(:,ii)*zr2*ss_plasma%m_sgas(ii)
       else if(imapx(ii).gt.0) then
          call fmap('impurity',zna(:,ii),zta(:,ii),zvpar(:,ii))
          call frantic_profile_set(imapx(ii),zspline1,zspline2,zspline3)
       endif
    enddo

    zmoi = zmoi/zntot  ! gas species density averaged moment of inertia per ion

  end subroutine init_prof

  subroutine free_loc_prof

    deallocate(zx,zf,zrmj,zvol,zwk1,zwk2,zxzcx,zrmjzc,zfzc,zsumzc,zresult)
    deallocate(zr2,zmoi,zscx,zpcx,ztcx,zntot,zdvol,zti,zomeg,zresult2)
    deallocate(zta,zna,zvpar)
    deallocate(zte,zne,zvpare)

    call ezspline_free(zspline1, iok)
    call ezspline_free(zspline2, iok)
    call ezspline_free(zspline3, iok)

  end subroutine free_loc_prof

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

    if(ierr.ne.0) then
       write(ilun_msgs,*) &
            ' ?ps_nbi_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 fmap1(lbl,zf)
    ! set up 1 profile

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

    real*8 :: zbuf(inzons+2)

    call insert(zf,zbuf)
    call ezspline_setup(zspline1, zbuf, iok)
    if(iok.ne.0) write(ilun_msgs,*) &
         ' %ps_nbi_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*8, dimension(:), intent(in) :: zn
    real*8, dimension(:), intent(in) :: zt
    real*8, dimension(:), intent(in) :: zvpar

    real*8 :: zbuf(inzons+2)

    call insert(zn,zbuf)
    call ezspline_setup(zspline1, zbuf, iok)
    if(iok.ne.0) write(ilun_msgs,*) &
         ' %ps_nbi_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,*) &
         ' %ps_nbi_frantic: EZspline warning on T, lbl=',trim(lbl)

    call insert(zrmjzc*zfzc*zvpar,zbuf)
    call ezspline_setup(zspline3, zbuf, iok)
    if(iok.ne.0) write(ilun_msgs,*) &
         ' %ps_nbi_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*8, dimension(:), intent(in) :: zprof
    real*8, dimension(:), intent(out) :: zspl_data

    real*8 :: zlim1,zlim2

    INTEGER, PARAMETER :: R8=SELECTED_REAL_KIND(12,100)

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

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

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

    zlim1=min(0.8_R8*zprof(inzons),1.2_R8*zprof(inzons))
    zlim2=max(0.8_R8*zprof(inzons),1.2_R8*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...

    integer :: idi

    do igi=1,ing
       ik=i0map(igi)
       if(ik.gt.0) then
          zsumzc = ZERO

          idi=0
          if(allocated(ss_sces%id_rate_sinb0i)) then
             idi=ss_sces%id_rate_sinb0i(igi)
          endif
          if(idi.gt.0) then
             call ps_intrp_1d(zxzcx(2:inzons+1),idi,zfzc,ierr, state=ss_sces)
             zsumzc=zsumzc+zfzc
          endif

          idi=0
          if(allocated(ss_sces%id_rate_sinf0i)) then
             idi=ss_sces%id_rate_sinf0i(igi)
          endif
          if(idi.gt.0) then
             call ps_intrp_1d(zxzcx(2:inzons+1),idi,zfzc,ierr, state=ss_sces)
             zsumzc=zsumzc+zfzc
          endif

          call fmap1('sinb0i',zsumzc)
          call frantic_bsink_ii(ik,zspline1)

          ix=0
          do ii=1,ss_sces%nspec_beam
             ix=ix+1
             idi=0
             if(allocated(ss_sces%id_rate_sinb0xs)) then
                idi=ss_sces%id_rate_sinb0xs(igi,ii)
             endif
             if(idi.ne.0) then
                call ps_intrp_1d(zxzcx(2:inzons+1), &
                     idi,zfzc,ierr, state=ss_sces)
             else
                zfzc=0.0d0
             endif
             call fmap1('sinb0xs',zfzc)
             call frantic_bsink_cx(ik,ix,zspline1)
          enddo

          do ii=1,ss_sces%nspec_fusion
             ix=ix+1
             idi=0
             if(allocated(ss_sces%id_rate_sinf0xs)) then
                idi=ss_sces%id_rate_sinf0xs(igi,ii)
             endif
             if(idi.ne.0) then
                call ps_intrp_1d(zxzcx(2:inzons+1), &
                     idi,zfzc,ierr, state=ss_sces)
             else
                zfzc=0.0d0
             endif
             call fmap1('sinf0xs',zfzc)
             call frantic_bsink_cx(ik,ix,zspline1)
          enddo
       endif
    enddo  ! resultant neutral species loop

  end subroutine load_bsink
end subroutine ps_frantic

subroutine ps_frantic_ccss(iobj_in,iobj_out,ilun_msgs,iwall,ihalo,ireco,ierr)

  ! c or c++ callable interface to ps_frantic(...)

  use plasma_state_mod
  implicit NONE

  integer :: iobj_in(*)   ! (ps_container_size) words used: input state object
  integer :: iobj_out(*)  ! (ps_container_size) words used: output state object

  integer, intent(in) :: ilun_msgs

  integer, intent(in) :: iwall  ! .ne.0: compute wall source neutral transport
  integer, intent(in) :: ihalo  ! .ne.0: compute halo source neutral transport
  integer, intent(in) :: ireco  ! .ne.0: compute recombination source transport

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

  !---------------------------------------
  type (plasma_state_container) :: ptr_in,ptr_out
  type (plasma_state), pointer :: ss_in,ss_sces
  !---------------------------------------

  call load_ps_container(iobj_in,ptr_in)
  ss_in => ptr_in%ps_ptr

  call load_ps_container(iobj_out,ptr_out)
  ss_sces => ptr_out%ps_ptr

  call ps_frantic(ss_in,ss_sces,ilun_msgs,iwall,ihalo,ireco,ierr)

end subroutine ps_frantic_ccss
