!-----------------------------------------------------------------------
!  depfus -- deposit monte carlo fusion product fast ions
!
subroutine depfus(ierr)
  !
  use xstraln_calls
  use nbi_com
  use nbi_random_dep
  !
  !  this routine is similar in spirit to subroutine dep, and is called
  !  in the same way, but, it does not have to deal with neutral tracks.
  !  fusion product fast ion source functions (vs. 2d beam spatial grid
  !  and reaction index) are stored in common array "sfusion".
  !
  !  first version dmc 14 sept 1990.  fpp support deferred
  !
  !  uda modification for theta allocation pms 14 july 1994
  !
  !------------------------------------------------------------------
  !%      use sglib_mod
  !
  !============
  implicit none
  integer :: idum,iz_save,ireact_save,ireact,jz2,iz2,izr,iz,jz1,jz
  integer :: jz2p1,iz1,iz_nz,ireact_nz,idummy,igc,iwarn,inoflr
  integer :: iuncon,itrap,jiav
  !============
  real(fp) :: zwghta,zpt,zsum,zsump,ztolp1,nbi_d1mach,zefusn
  real(fp) :: zvfusn,zcth,zv1,zv1pp,zvth2,zv2,zv3,zrth,zdth,zomega
  real(fp) :: fnbint,zr,zy,zvrot,zvt,zvr,zvy,zvsq,zphi0,zphigc,ze0
  real(fp) :: zegc,zbt0,zbtgc,zrgc,zygc,zvpll,zxiav,xjadep,zvpov
  real(fp) :: zeiav,zphiav,zxhw
  real(fp) :: zxiflr,zthflr
  !============
  integer, intent(out) :: ierr  ! status code, 0=ok
  !
  !  debugging:
  !
  !  ztdep common replaced by nbi_com variables...
  !
  !  common block also used in depnrm -- testing fpp first orbit loss
  !    accounting
  !   xn_dep,xn_orb-- #/sec deposited, left orbitting
  !   xp_dep,xp_orb-- power deposited, left orbitting
  !   xph_dep,xph_orb-- ang. mo./sec deposited, left orbitting
  !
  !   differences to be accounted by:  prompt orbit losses
  !   jxb torque and associated work term
  !
  !  plotting
  !
  !%      real(fp) :: zzr(mj),zzy(mj),zzr2(2),zzy2(2),zrflr,zyflr,zth,zrlar
  !%      integer :: inpts,ith,iad
  !
  real(fp) :: zdummy(2,2)
  !------------------------------------------------------------------
  !
  !  local:
  !
  real(fp) :: zjac(2,2)
  !
  !  set .true. for reactions that produce a desired fusion product:
  logical ilreact(mireact)
  !
  !  initial energy (ev) of fusion product so produced:
  real(fp) :: zereact(mireact)
  !
  real(fp) :: zvol_randm,zvincr,zx1,zx2
  real(fp), dimension(:), allocatable :: zvols
  !
  !-------------------------------------------------------------
  !
  !
  !%      zzr=0; zzy=0; zzr2=0; zzy2=0
  !%      if(nbplot.gt.0) then
  !%        zzr(1)=mj-1
  !%        zzy(1)=mj-1
  !%        zzr2(1)=1.0_fp
  !%        zzy2(1)=1.0_fp
  !%      end if
  !
  !
  !---------------------------------------------------------------------
  !
  ierr=0
  !
  !  debug code
  !
  kdebug1=0
  kdebug2=0
  !
  if(kdebug1.eq.1) then
     if(ewghtdep(lsbeam).eq.0.0_fp) then
        xn_dep(lsbeam)=0.0_fp
        xn_orb(lsbeam)=0.0_fp
        xp_dep(lsbeam)=0.0_fp
        xp_orb(lsbeam)=0.0_fp
        xph_dep(lsbeam)=0.0_fp
        xph_orb(lsbeam)=0.0_fp
     end if
     zwghta=1.0_fp
  end if
  !
  !---------------------------------------------------------------------
  !  error check
  !
  if(.not.nlfprod(lsbeam)) then
     write(nonlin,*) ' ?? depfus called for neutral beam deposition'
     write(nonlin,*) '    should call subroutine dep'
     ierr=1; return
  end if
  !---------------------------------------------------------------------
  !
  !  deposition
  !
  if(nlfbmfpp(lsbeam)) then
     write(nonlin,*) ' ?? depfus not ready to deposit for fppmod'
     ierr=1; return
  else
     !
     !  monte carlo launch.
     !  1.  choose random start time with timestep
     !
     bitime(1)=nbrandm_dep(idum)*delt
     !
     !  nbsce=0 indicates source is a fusion reaction, not a neutral beam
     !  set nbien=1 as there is only one "energy fraction" for fusion
     !  products
     !
     nbsce=0
     nbien=1
     !
     !  2.  find source point and energy.
     !      search through source rate array by position and reaction index
     !
     call fuslis(lsbeam,ilreact,zereact)
     !
     !      stop search at this random fraction of the total source:
     !      (see subroutine fuspwr which summed the source)
     !
     zpt=nbrandm_dep(idum)*xninjs(lsbeam)
     !
     iz_save=0
     ireact_save=0
     !
     continue190 :     do ! re-enter in case of close "source pt not found"
        zsum=0.0_fp
        !
        do ireact=1,nreact
           if(ilreact(ireact)) then
              !
              jz2=lcm1
              iz2=nthzsm(lcm1)
              izr=lcm1
              !
              do iz=1,nfbzns
                 if(iz.gt.iz2) then
                    !  new zone row
                    jz1=min(lep1,(jz2+1))
                    jz2=min(lep1,(jz1+nxskpb-1))
                    jz2p1=min(lep1,(jz2+1))
                    !  new theta indexing
                    izr=izr+1
                    iz1=iz2+1
                    iz2=nthzsm(izr)
                 end if
                 
                 !
                 zsump=zsum
                 zsum=zsum+sfusion(iz,ireact)
                 if((zsump.le.zpt).and.(zpt.le.zsum)) exit continue190
                 if((iz.eq.iz_save).and.(ireact.eq.ireact_save)) then
                    exit continue190           ! recovery path
                 end if
                 if(sfusion(iz,ireact).gt.0.0_fp) then
                    iz_nz=iz            ! save in case we fall through...
                    ireact_nz=ireact
                 end if
              end do
           end if

        end do
     !
        ztolp1=1.0_fp+10.0_fp*nbi_d1mach(4)
        if(zpt.gt.ztolp1*zsum) then
           write(nonlin,*) ' ===> error in depfus:'
        else
           write(nonlin,*) ' ===> warning in depfus:'
        end if
        write(nonlin,*) ' random point zpt = ',zpt,' is past the end of'
        write(nonlin,*)' the source distribution (sfusion sum=',zsum,')'
        if(zpt.gt.ztolp1*zsum) then
           write(nonlin,*) ' ?depfus: error tolerance exceeded.'
           ierr=1
           return
        endif
        write(nonlin,*) ' (fixup applied).'
        iz_save=iz_nz
        ireact_save=ireact_nz
        cycle 
     enddo continue190
     !
     !
     !  source point located.
     !  fusion product energy and velocity in plasma frame:
     !
     zefusn=zereact(ireact)
     zvfusn=sqrt(zefusn/v2toev)
     !
     !  velocity components (plasma frame):
     !  bugfix ** dmc 16 july 1993
     !bad        zvth1=xpi*nbrandm_dep(idum)
     !bad        zv1=zvfusn*cos(zvth1)
     !bad        zv1pp=zvfusn*sin(zvth1)
     zcth=2.0_fp*nbrandm_dep(idum)-1.0_fp
     zv1=zvfusn*zcth
     zv1pp=zvfusn*sqrt(1.0_fp-zcth**2)
     !  end bugfix
     zvth2=twopi*nbrandm_dep(idum)
     zv2=zv1pp*sin(zvth2)
     zv3=zv1pp*cos(zvth2)
     !
     !  fusion source point within the beam 2d spatial zone:
     !  xi:
     !xxxx        xiflr=xi(jz1,1)+nbrandm_dep(idum)*(xi(jz2p1,1)-xi(jz1,1))
     !
     !     ** take volumes into account ** (bugfix dmc apr 2007)
     !
     allocate(zvols(jz1:jz2p1))
     zvols(jz1)=0.0_fp
     do jz=jz1,jz2
        zvols(jz+1)=zvols(jz)+dvolzc(jz)
     end do

     zvol_randm = zvols(jz2p1)*nbrandm_dep(idum)

     xiflr = xi(jz2p1,1)  ! just in case...
     do jz=jz1,jz2
        if(zvol_randm.le.zvols(jz+1)) then

           zvincr = (zvol_randm-zvols(jz))/(zvols(jz+1)-zvols(jz))
           zx1=xi(jz,1)
           zx2=xi(jz+1,1)
           xiflr = max(zx1, &
                sqrt(zx1*zx1 + zvincr*(zx2*zx2-zx1*zx1)) )

           exit
        end if
     end do

     deallocate(zvols)
     !
     !  theta: including non-symmetric case
     !
     xistate=xiflr      ! for fast ion radial current tracking
     !
     zrth=2.0_fp*nbrandm_dep(idum)-1.0_fp
     zdth=(thbdy1-thbdy0)/(iz2-iz1+1)
     thflr=thbdy0+(iz-iz1)*zdth+abs(zrth)*zdth
     if(nlsym2b) then
        if(zrth.lt.0.0_fp) thflr=-thflr
     end if
     !
     !  local angular velocity
     !
     if(jz1.eq.lep1) then
        zomega=0.0_fp
     else
        zomega = fnbint(xiflr,omegag,1,mj,1,idummy,igc)
     end if
     !
     !  (r,y) position at source point
     !
     !test        call xmomry(xiflr,thflr,0,zr,zy,zjac)
     call xs_rzget(0.0_fp,thflr,xiflr,100.0_fp,zr,zy,iwarn)
     !
     !  local plasma rotation velocity
     !
     zvrot=zr*zomega
     !
     !  lab frame velocity
     !
     zvt=zvrot+zv1
     zvr=zv2
     zvy=zv3
     zvsq=zvt*zvt+zvr*zvr+zvy*zvy
     vion=sqrt(zvsq)         ! converted to gc value in depsub
     !
     !  get initial orbit parameters
     !
     if(nlbflr) then
        inoflr=2
     else
        inoflr=3
     end if
     !
     !  depsub may move (xiflr,thflr) slightly...
     !
     call depsub(zr,zy,zvt,zvr,zvy,inoflr, &
          zphi0,zphigc,ze0,zegc,zbt0,zbtgc)
     !
     !
     !  initial weight of ion and energy weighted sum
     !
     wght=1.0_fp
     ewghtdep(lsbeam)=ewghtdep(lsbeam)+wght*ze0
     wghtdep(lsbeam)=wghtdep(lsbeam)+wght
     !
     if((lcentr.le.ngc).and.(ngc.le.ledge)) then
        mc_newdens(ngc,lsbeam) = mc_newdens(ngc,lsbeam) + wght
     end if
     !
     !  momentum deposition
     bphdepsgc(lsbeam)=bphdeps(lsbeam)+wght*zphigc
     bphdeps(lsbeam)=bphdeps(lsbeam)+wght*zphi0
     bphdepsgc(lsbeam)=bphdeps(lsbeam)+wght*zphigc
     !
     !  flr jxb effect
     !
     call nb0jxb(wght,ze0,zegc,zphi0,zphigc,zbt0,zbtgc, &
          nflr,ngc,xiflr,xion,1,.false.)
     !
     !  estimated average orbit location
     !
     !test        call xmomry(xion,th,0,zrgc,zygc,zdummy)
     call xs_rzget(0.0_fp,th,xion,100.0_fp,zrgc,zygc,iwarn)
     zvpll=xksid*vion
     !
     if(xion.gt.1.0_fp) then
        iuncon=1
        jiav=lep1 ! unconfined
     else
        zxiav=xjadep(xion,th,zrgc,zygc,vion,zvpll, &
             itrap,iuncon,zvpov,zeiav,zphiav,jiav,zxhw)
        !  trapping fraction sums
        call addtrap0(wght,itrap,zeiav)
     end if
     sdep_oa(jiav,lsbeam)=sdep_oa(jiav,lsbeam)+wght
     zidep0(lsbeam)=zidep0(lsbeam)+1 !counter of #deposited ptcls on given cpu
     !
     if(iuncon.eq.1) then
        xjav=real(lep1)             ! probable lost orbit
     else
        xjav=lcentr+zxiav*nzones     ! probable confined orbit
     end if
     !
  end if
  !
  einj0=v2toev*vion*vion            ! lab frame birth energy
  !
  !%      if(nbplot.gt.0) then
  !%        call xmomry(xiflr,thflr,0,zrflr,zyflr,zdummy)
  !%        zzr2(2)=zrflr
  !%        zzy2(2)=zyflr
  !%        call symbl(2)
  !%        call line(-1)
  !%        call r8_cplot(zzr2,zzy2)
  !
  !%        call line(0)
  !%        call symbl(0)
  !
  !%        zrlar=sqrt((zrgc-zrflr)**2+(zygc-zyflr)**2)
  !
  !%        inpts=mj-1
  !%        do 999 ith=1,inpts
  !%          zth=(ith-1)*twopi/(inpts-1)
  !%          iad=ith+1
  !%          zzr(iad)=zrgc+zrlar*cos(zth)
  !%          zzy(iad)=zygc+zrlar*sin(zth)
  !% 999    continue
  !%        call r8_cplot(zzr,zzy)
  !%        call tsend
  !%      end if
  return
end subroutine depfus
! 11sep2004 fgtok
! 11sep2004 fgtok
