subroutine gfbm_set_track_data

  ! from an XPLASMA object fetch NUBEAM track data if available
  ! also, get available fusion product deposition data
  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: zero,zrdg=>deg2rad,&
       z1=>one,z1000=>zkev2ev,z100=>zmtocm
  use get_fbm_module
  use deptrk_io_track
  use xplasma_definitions
  implicit NONE

  !-------------------------------
  integer :: ierr

  logical :: rfdeut = .FALSE.
  logical :: rftrit = .FALSE.
  logical :: rfhel3 = .FALSE.
  logical :: nlfst,nlfhe3
  integer :: itrit,iprot(2),ihel3,ihel4,itot,isb,iz,ia,ir,ireact,inum
  integer :: ibi,index_nbfi
  integer :: id_mcgrid,ist_th,ist_ph
  integer :: imax_ndep,imax_track,inztrack
  integer, dimension(:), allocatable :: inz_track !total # of zone per species
  integer :: inzmax,indep,ip,ip_tr
  integer :: imlost_data, imax_lost, imibs_use
  integer :: ii,nsfast_use,isfast_use  
  real(fp) :: zconv,zsftot

  character*30 :: cvar
  character*2 znum2
  character*32 zlbls(100),zlbl  ! fast ion labels
  character*150 fname_closed !NetCdf file name for track data
!-------------------------------
  ! initially, assume no data...
  ndep0_use = -1
  nldeptrk = .FALSE.

  call fuslis_mireact(mireact)
  if(mireact.le.0) then
     write(6,*) ' %gfbm_set_track_data: mireact=0, early exit.'
     return
  endif

  if(nznsi.le.0) then
     write(6,*) ' %gfbm_set_track_data: nznsi=0, early exit.'
     return
  endif

  call mcgrid_find(id_mcgrid)
  if(id_mcgrid.le.0) then
     write(6,*) ' %gfbm_set_track_data: id_mcgrid.le.0, early exit.'
     return
  endif

  !-------------------------------
  !  look at fusion products

  if(allocated(nlfusion_sce)) deallocate(nlfusion_sce)
  if(allocated(efusion_sce)) deallocate(efusion_sce)

  allocate(nlfusion_sce(mireact,nsfast)); nlfusion_sce = .FALSE.
  allocate(efusion_sce(mireact,nsfast)); efusion_sce = ZERO

  nlfst=.FALSE.
  nlfhe3=.FALSE.
  iprot = 0
  itrit = 0
  ihel3 = 0
  ihel4 = 0
  itot = 0
  nsfast_use = 0
  nldep_beam = .FALSE.
  ibi=1
  do isb=1,nsfast
     if(nlfprod(isb)) then
        iz = xzbeams(isb)+0.5
        ia = abeams(isb)+0.5

        if(iz.eq.1) then
           if(ia.eq.1) then
              iprot(ibi)=isb
              itot=itot+1
              if(ibi.gt.2) then
                 write(6,*) ' ?gfbm_set_track_data: only fusion p(3.02MeV) and p(14.7MeV) allowed'
                 exit
              endif
              ibi=ibi+1
           else if(ia.eq.3) then
              itrit=isb
              nlfst=.TRUE.
              itot=itot+1
           endif
        else if(iz.eq.2) then
           if(ia.eq.3) then
              ihel3=isb
              nlfhe3=.TRUE.
              itot=itot+1
           else if(ia.eq.4) then
              ihel4=isb
              itot=itot+1
           endif
        endif
     else
        nsfast_use = nsfast_use +1
        nldep_beam(isb) = .TRUE.
     endif
  enddo
  ierr = 0
  do isb=1,nsfast
     if(nlfprod(isb).and.(itot.gt.0)) then
        call fuslis_sub(mireact,nlfusion_sce(1,isb),efusion_sce(1,isb), &
             (isb.eq.itrit),(isb.eq.iprot(1)),(isb.eq.iprot(2)),(isb.eq.ihel3),(isb.eq.ihel4), &
             nlfst,nlfhe3,rfdeut,rftrit,rfhel3,ierr)
        if(ierr.ne.0) exit
     endif
  enddo

  if(ierr.ne.0) then
     write(6,*) ' ?gfbm_set_track_data: unexpected fuslis_sub error:'// &
          ' taking early exit.'
     return
  endif

  ! OK now get the sources

  if(allocated(fusion_sce)) deallocate(fusion_sce)
  allocate(fusion_sce(nznsi,mireact)); fusion_sce = ZERO

  if(nlsym) then
     ist_th=0
  else
     ist_th=-1
  endif
  ist_ph=0

  zconv = 1.0d0  ! sce in #/zone/sec

  ierr = 0
  if(itot.gt.0) then
     do ireact = 1,mireact
        write(znum2,'(i2.2)') ireact
        zlbl='sfusion_'//znum2

        call mcgrid_getobj(-id_mcgrid,zlbl, &
             ist_th,ist_ph,zconv,fusion_sce(1,ireact),nznsi, ierr)

        if(ierr.ne.0) exit
     enddo
  endif

  if(ierr.ne.0) then
     write(6,*) ' %gfbm_set_track_data: mcgrid_getobj(...,fusion_sce,...):'// &
          ' error, early exit.'
     return
  endif

  if(allocated(sftot)) deallocate(sftot)
  allocate(sftot(mireact))
  sftot(1:mireact) = ZERO
 
  if(itot.gt.0) then

     ! check if there are really fusion sources

     itot = 0
     do isb=1,nsfast
        zsftot = ZERO
        do ir=1,mireact
           if(nlfusion_sce(ir,isb)) then
              if(sftot(ir).eq.ZERO) then
                 do iz=1,nznsi
                    sftot(ir) = sftot(ir) + fusion_sce(iz,ir)
                 enddo
              endif
              zsftot=zsftot+sftot(ir)
           endif
        enddo
        nldep_fusion(isb) = (zsftot.gt.ZERO)
        if(nldep_fusion(isb)) itot=itot+1
     enddo
  endif

  if(itot.gt.0) then
     ! OK there are fusion product sources, at least...

     ndep0_use=0
     nldeptrk = .TRUE.
  endif

  !-------------------------------
  !  look for beam track data

  nstep_gtime = 1
  track_gtime1 = ta-dta
  track_gtime2 = ta+dta

  zlbls = ' '
  call nbfi_lbl(zlbls,100,inum)

  imax_track = 0
  inztrack = 0
  imax_ndep = 0

  if(nsfast_use.gt.0) then
     if(allocated(ntracks)) deallocate(ntracks)
     allocate(ntracks(nsfast_use)); ntracks = 0

     if(allocated(inz_track)) deallocate(inz_track)
     allocate(inz_track(nsfast_use)); inz_track = 0
  endif
  isfast_use=0
  do isb=1,nsfast
     ibi = 0
     if(.not.nlfprod(isb)) then
        iz = xzbeams(isb) + 0.5
        ia = abeams(isb) + 0.5
        ibi = index_nbfi(iz,ia,1)
     endif
     if(ibi.gt.0.and..not.nlfprod(isb).and.nldep_beam(isb)) then

        isfast_use=isfast_use+1
     endif
  enddo

  call deptrk_num(nsfast_use,imax_ndep,imax_track,inz_track(1:nsfast_use),imibs_use,iversion_deptrk)
  ntracks(1:isfast_use)=imax_ndep !same # of tracks is used for all species

  if(imax_ndep.eq.0) then
     write(6,*) ' %gfbm_set_track_data: no beam deposition track data found.'
     nldep_beam = .FALSE.
     return
  endif
  inztrack = sum(inz_track(1:isfast_use)) !total zone for all tracks and all beam ion species
  isfast_use=0
  do isb=1,nsfast
     ibi = 0
     if(ibi.gt.0.and..not.nlfprod(isb).and.nldep_beam(isb)) then

        isfast_use=isfast_use+1

        if(ntracks(isfast_use).gt.0) then
           nldep_beam(isb)=.TRUE.
        else
           nldep_beam(isb)=.FALSE.
        endif

     endif
  enddo

  if(iversion_deptrk.ne.3) then
     write(6,*) ' %gfbm_set_track_data: track data black box version=',iversion_deptrk
     write(6,*) '  (this routine is coded for version 3 only).'
     nldep_beam = .FALSE.
     return
  else
     nldeptrk = .TRUE.
  endif

  nzmax_track = imax_track
  mztrack = nzmax_track
  inzmax = nzmax_track

  ndep0_use = imax_ndep
  mntrack = ndep0_use
  indep = ndep0_use

  if(indep.gt.0.and.inztrack.gt.0) then
     ! allocate track data arrays (note: nfast = nsfast)

     if(allocated(nstat_track_xja)) deallocate(nstat_track_xja)
     allocate(nstat_track_xja(inztrack))
     nstat_track_xja = 0

     if(allocated(nztrack)) deallocate(nztrack)
     allocate(nztrack(mntrack,nsfast_use))
     nztrack = 0

     if(allocated(ntrack_nbsce)) deallocate(ntrack_nbsce)
     allocate(ntrack_nbsce(mntrack,nsfast_use))
     ntrack_nbsce = 0

     if(allocated(ntrack_nbien)) deallocate(ntrack_nbien)
     allocate(ntrack_nbien(mntrack,nsfast_use))
     ntrack_nbien = 0

     if(allocated(track_de_flr)) deallocate(track_de_flr)
     allocate(track_de_flr(inztrack))
     track_de_flr = ZERO

     if(allocated(track_dvpv_flr)) deallocate(track_dvpv_flr)
     allocate(track_dvpv_flr(inztrack))
     track_dvpv_flr = ZERO

     if(allocated(track_dR_flr)) deallocate(track_dR_flr)
     allocate(track_dR_flr(inztrack))
     track_dR_flr = ZERO

     if(allocated(track_dZ_flr)) deallocate(track_dZ_flr)
     allocate(track_dZ_flr(inztrack))
     track_dZ_flr = ZERO

     if(allocated(track_einj)) deallocate(track_einj)
     allocate(track_einj(mntrack,nsfast_use))
     track_einj = ZERO

     if(allocated(track_sinj)) deallocate(track_sinj)
     allocate(track_sinj(mntrack,nsfast_use))
     track_sinj = ZERO

     if(allocated(track_time)) deallocate(track_time)
     allocate(track_time(mntrack,nsfast_use))
     track_time = ZERO

     if(allocated(track_pdep)) deallocate(track_pdep)
     allocate(track_pdep(inztrack))
     track_pdep = ZERO

     if(allocated(track_xl)) deallocate(track_xl)
     allocate(track_xl(inztrack))
     track_xl = ZERO

     if(allocated(track_vpv_dep)) deallocate(track_vpv_dep)
     allocate(track_vpv_dep(inztrack))
     track_vpv_dep = ZERO

     if(allocated(track_rin)) deallocate(track_rin)
     allocate(track_rin(mntrack,nsfast_use))
     if(allocated(track_rout)) deallocate(track_rout)
     allocate(track_rout(mntrack,nsfast_use))
     track_rin = ZERO
     track_rout = ZERO

     if(allocated(track_phiin)) deallocate(track_phiin)
     allocate(track_phiin(mntrack,nsfast_use))
     if(allocated(track_phiout)) deallocate(track_phiout)
     allocate(track_phiout(mntrack,nsfast_use))
     track_phiin = ZERO
     track_phiout = ZERO

     if(allocated(track_zin)) deallocate(track_zin)
     allocate(track_zin(mntrack,nsfast_use))
     if(allocated(track_zout)) deallocate(track_zout)
     allocate(track_zout(mntrack,nsfast_use))
     track_zin = ZERO
     track_zout = ZERO
  endif

  ierr = 0
  mtrk0 = 0
  isfast_use = 0
  ip_tr = 1
  do isb=1,nsfast
     ibi = 0
     if(.not.nlfprod(isb)) then
        iz = xzbeams(isb) + 0.5
        ia = abeams(isb) + 0.5
        ibi = index_nbfi(iz,ia,1)
     endif
     if(.not.nlfprod(isb).and.nldep_beam(isb)) then
        isfast_use = isfast_use +1
        zlbl = zlbls(ibi)
        ip = 1
        ip_tr = ip_tr+ mtrk0

        !------------------------
        !  get INTEGER data


        cvar = 'xntrack'
        call deptrk_io(cvar,z1,nztrack(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/),ierr); call ckerr
        cvar = 'track_nbsce'
        call deptrk_io(cvar,z1,ntrack_nbsce(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/),ierr); call ckerr
        cvar = 'track_nbien'
        call deptrk_io(cvar,z1,ntrack_nbien(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/),ierr); call ckerr

        mtrk0 = inz_track(isfast_use) !#of tracks for a given species
        cvar = 'track_xjadep'
        call deptrk_io(cvar,z1,nstat_track_xja(ip_tr:ip_tr+mtrk0-1),(/ip_tr/), &
             ierr); call ckerr

        !------------------------
        !  get REAL(FP) 1d data

        !  cm, cm, rad
        cvar = 'track_rin'
        call deptrk_io(cvar,z1,track_rin(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        cvar = 'track_zin'
        call deptrk_io(cvar,z1,track_zin(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        cvar = 'track_phiin'
        call deptrk_io(cvar,zrdg,track_Phiin(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        
        !  cm, cm, rad
        cvar = 'track_rout'
        call deptrk_io(cvar,z1,track_rout(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        cvar = 'track_zout'
        call deptrk_io(cvar,z1,track_zout(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        cvar = 'track_phiout'
        call deptrk_io(cvar,zrdg,track_Phiout(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        
        !  s, eV, #ptcls
        cvar = 'track_time'
        call deptrk_io(cvar,z1,track_time(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        cvar = 'track_einj'
        call deptrk_io(cvar,z1,track_einj(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        cvar = 'track_sinj'
        call deptrk_io(cvar,z1,track_sinj(ip:ip+indep-1,isfast_use:isfast_use),(/ip,isfast_use/), &
             ierr); call ckerr
        
        !------------------------
        !  get REAL(FP) 2d data
        
        !  dimensionless (normalized distance), dimensionless (depo. probability)
        
        cvar = 'track_xl'
        call deptrk_io(cvar,z1,track_xL(ip_tr:ip_tr+mtrk0-1),(/ip_tr/), &
             ierr); call ckerr
        cvar = 'track_pdep'
        call deptrk_io(cvar,z1,track_pdep(ip_tr:ip_tr+mtrk0-1),(/ip_tr/), &
             ierr); call ckerr
        
        !  cm, cm, eV
        cvar = 'track_dr_flr'
        call deptrk_io(cvar,z1,track_dR_flr(ip_tr:ip_tr+mtrk0-1),(/ip_tr/), &
             ierr); call ckerr
        cvar = 'track_dz_flr'
        call deptrk_io(cvar,z1,track_dZ_flr(ip_tr:ip_tr+mtrk0-1),(/ip_tr/), &
             ierr); call ckerr
        cvar = 'track_de_flr'
        call deptrk_io(cvar,z1,track_dE_flr(ip_tr:ip_tr+mtrk0-1),(/ip_tr/), &
             ierr); call ckerr
        
        !  dimensionless (vpll/v @dep), dimensionless (delta(vpll/v) to GC)
        
        cvar = 'track_vpv_dep'
        call deptrk_io(cvar,z1,track_vpv_dep(ip_tr:ip_tr+mtrk0-1),(/ip_tr/), &
             ierr); call ckerr
        cvar = 'track_dvpv_flr'
        call deptrk_io(cvar,z1,track_dvpv_flr(ip_tr:ip_tr+mtrk0-1),(/ip_tr/), &
             ierr); call ckerr
        
     endif
  enddo
  call deptrk_io_close(fname_closed)
  if(ierr.ne.0) then
     write(6,*) ' ?gfbm_set_track_data: error reading track data.'
     nldep_beam = .FALSE.
  endif
  !remove <RUNID>_dpt.cdf file
  ierr=0
  call fdelete(fname_closed,ierr)
  if(ierr.ne.0) then
     write(6,*) ' ?gfbm_set_track_data: cannot delete file: '//trim(fname_closed)
     return
  endif
  
  if (allocated(inz_track)) deallocate(inz_track)
  !-------------------------------
  !  look for lost orbit data
  imax_lost=0
  ip = 1
  gtime1_lost = ta-dta
  gtime2_lost = ta+dta
  call lost_orbit_num(' ',imlost_data)
  imax_lost=max(imax_lost,imlost_data)
  mlost_data = imax_lost
  if(imax_lost.eq.0) then
     nllost = .false.
     return
  endif
  nllost = .TRUE.

  if(allocated(gtime_lost)) deallocate(gtime_lost)
  allocate(gtime_lost(imax_lost)); gtime_lost = zero
  if(allocated(ynbsce_lost)) deallocate(ynbsce_lost)
  allocate(ynbsce_lost(imax_lost)); ynbsce_lost = zero
  if(allocated(ynbien_lost)) deallocate(ynbien_lost)
  allocate(ynbien_lost(imax_lost)); ynbien_lost = zero
  if(allocated(wght_lost)) deallocate(wght_lost)
  allocate(wght_lost(imax_lost)); wght_lost = zero
  if(allocated(rmjion_lost)) deallocate(rmjion_lost)
  allocate(rmjion_lost(imax_lost)); rmjion_lost = zero
  if(allocated(yyion_lost)) deallocate(yyion_lost)
  allocate(yyion_lost(imax_lost)); yyion_lost = zero
  if(allocated(yze_lost)) deallocate(yze_lost)
  allocate(yze_lost(imax_lost)); yze_lost = zero
  if(allocated(xksid_lost)) deallocate(xksid_lost)
  allocate(xksid_lost(imax_lost)); xksid_lost = zero
  if(allocated(goose_lost)) deallocate(goose_lost)
  allocate(goose_lost(imax_lost)); goose_lost = zero
  if(allocated(species_lost)) deallocate(species_lost)
  allocate(species_lost(imax_lost)); species_lost = zero
  zlbl = ' '
  !  s
  call lost_orbit_r8get1(zlbl,'TIMELST',z1,imax_lost,gtime_lost(ip), &
       ierr)
  if(ierr.ne.0) then
     nllost = .false.
     write(6,*) ' ?gfbm_set_track_data: error reading lost orbit data.'
     return
  endif
  call lost_orbit_r8get1(zlbl,'NBSCELST',z1,imax_lost,ynbsce_lost(ip), &
       ierr) 
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  call lost_orbit_r8get1(zlbl,'NBIENLST',z1,imax_lost,ynbien_lost(ip), &
       ierr) 
  !   #ptcls
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  call lost_orbit_r8get1(zlbl,'WGHTLST',z1,imax_lost,wght_lost(ip), &
       ierr) 
  !  cm, cm
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  call lost_orbit_r8get1(zlbl,'RMJION',z100,imax_lost,rmjion_lost(ip), &
       ierr) 
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  call lost_orbit_r8get1(zlbl,'ZYION',z100,imax_lost,yyion_lost(ip), &
       ierr) 
  !  eV
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  call lost_orbit_r8get1(zlbl,'ZELST',z1000,imax_lost,yze_lost(ip), &
       ierr) 
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  call lost_orbit_r8get1(zlbl,'XKSIDLST',z1,imax_lost,xksid_lost(ip), &
       ierr) 
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  call lost_orbit_r8get1(zlbl,'GOOSELST',z1,imax_lost,goose_lost(ip), &
       ierr)  
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  call lost_orbit_r8get1(zlbl,'SPECIESLST',z1,imax_lost,species_lost(ip), &
       ierr)  
  if(ierr.ne.0) then
     nllost = .false.
     return
  endif
  if(nllost.or.nldeptrk) then
     ip=1
     call beam_info_num(' ',nbeam,nsjdotb,nsnccw)
     if(nbeam.le.0)return

     if(allocated(znlco)) deallocate(znlco)
     if(allocated(pinja)) deallocate(pinja)
     if(allocated(xbzeta)) deallocate(xbzeta)
     if(allocated(rtcena)) deallocate(rtcena)
     if(allocated(xlbapa)) deallocate(xlbapa)
     if(allocated(xlbtna)) deallocate(xlbtna)

     allocate(znlco(nbeam),pinja(nbeam),xbzeta(nbeam))
     allocate(rtcena(nbeam),xlbapa(nbeam),xlbtna(nbeam))

     zlbl = ' '
     call beam_info_r8get1(zlbl,'NLCO',z1,nbeam,znlco(ip),ierr)
     if(ierr.ne.0) then
             write(6,*) ' ?gfbm_set_track_data: error reading beam_info nlco.'
             return
     endif
     call beam_info_r8get1(zlbl,'PINJA',z1,nbeam,pinja(ip),ierr)
     if(ierr.ne.0) then
             write(6,*) ' ?gfbm_set_track_data: error reading beam_info pinja.'
             return
     endif
     call beam_info_r8get1(zlbl,'XBZETA',z1,nbeam,xbzeta(ip),ierr)
     if(ierr.ne.0) then
             write(6,*) ' ?gfbm_set_track_data: error reading beam_info xbzeta.'
             return
     endif
     call beam_info_r8get1(zlbl,'XLBAPA',z1,nbeam,xlbapa(ip),ierr)
     if(ierr.ne.0) then
             write(6,*) ' ?gfbm_set_track_data: error reading beam_info xlbapa.'
             return
     endif
     call beam_info_r8get1(zlbl,'XLBTNA',z1,nbeam,xlbtna(ip),ierr)
     if(ierr.ne.0) then
             write(6,*) ' ?gfbm_set_track_data: error reading beam_info xlbtna.'
             return
     endif
     call beam_info_r8get1(zlbl,'RTCENA',z1,nbeam,rtcena(ip),ierr)
     if(ierr.ne.0) then
             write(6,*) ' ?gfbm_set_track_data: error reading beam_info rtcena.'
             return
     endif
  endif
  
contains
  
  subroutine ckerr
    
    if(ierr.ne.0) then
       write(6,*) ' ?gfbm_set_track_data: an error was detected.'
       call bad_exit
    endif
    
  end subroutine ckerr
  
  
end subroutine gfbm_set_track_data
