!******************** start file nbiflr.for ; group nbiflr ******************
!--------------------------------------------------------------------
!  nbiflr
!
!  evaluate the flr position (position out on the gyro) of the
!  orbitting beam ion
!
!  rga feb01 - added charge and reversed sign on ephflr calculation
!              added gyro based larmor displacement (map_larmor)
!
!  mod dmc march 1987 - for momentum conservation
!   ion velocities on gyro are not identical to the gyro-avg velocities
!   tracked in the orbit eqn
!
!   compute flr velocities based on conservation of
!    m*vpll*r and of vpll**2+vperp**2+electrostatic potential
!
!  mod pam stubberfield 15 july 1994 for uda coding
!   nlsym2b=f no negation of btiltc(4,..) and nflr2 2d beam zone index change
!
subroutine nbiflr(zangl,itrig)
  use map_larmor_mod   ! module for gyro based flr offset
  !
  !  zangl - gyrophase angle
  !  itrig - 1 to evaluate cos and sin of zangl
  !          0 if already done; if 0, zangl is not used
  !
  ! ****
  !  common variables cospho, sinpho store the cos and sin of the
  !    gyrophase angle
  ! ****
  !
  !-----------------------
  !
  use xstraln_calls
  use nbi_com
  !
  !  common input data pertaining to beam ion ***
  !    rmjion,yion   (r,y) of guiding ctr
  !    rlar          larmour radius
  !
  !    vion,xksid    velocity, vpll/v of guiding ctr
  !    ngc, xingc    guiding center zone and interpolation factor
  !
  !    map_larmor also uses
  !    xbjaco, vperp, vpar, cospho, sinpho
  !
  !  common output data pertaining to beam ion ***
  !    rmjflr,yflr   (r,y) of ion out on gyro
  !    xiflr,thflr   (xi,theta) of ion out on gyro
  !    nflr,nflr2    transp zone, beam zone indices of ion on gyro
  !
  !    ephflr        electrostatic potential out on gyro
  !
  !    vflr,vppflr,vplflr,xksflr
  !                  v, vperp, vpll, vpll/v out on gyro
  !
  !  the inverse fourier moments map (rmjflr,yflr) --> (xiflr,thflr)
  !    is done with an approximate bilinear interpolation
  !    care is taken near the xiflr=0 axial singularity
  !    the interpolant array is set up in subroutine parms
  !
  !  ptcls have been limlos checked!
  !    however, this does not gaurantee that the entire ion gyro lies
  !    inside the region covered by the beam dist. fcn. spatial grid.
  !    if the ptcl is found to lie outside, nflr2=0 is returned! and
  !    the caller should check for this...
  !
  ! 1/pi ...
  !============
  implicit none
  integer :: itrig,iflag,ioops,iregion,iwarn,ibr,inz0,inz
  !============
  real(fp) :: zangl,zpinv,z_vpar,z_vperp,zdrx,zdrth,ztop,zphiarg,zv
  real(fp) :: alpha_v ! poloidal aiming angle
  real(fp) :: zvec(3) !vr,vphi,vz velocity of the particle (work array here)
  !============
  data zpinv/0.3183098_fp/
  !-------------------------
  !
  if(itrig.eq.1) then
     cospho=cos(zangl)
     sinpho=sin(zangl)
  end if
  !
  !  evaluate the (r,y) flr position
  !
  iflag=0
  !
  if(nlbflr) then
     if (nlbgflr .and. rlar>max(epslon,gflr_min)) then
        ! --- use a gyro based displacement ---
        call map_larmor(ioops, iregion, rmjflr, yflr, &
             z_vpar, z_vperp,alpha_v, .false.)              ! accept all outcomes
        iflag=1
     else
        ! --- displacement from guiding center ---
        zdrx=rlar*cospho
        zdrth=rlar*sinpho*nsnccw*btiltc(2,ngc2)
        if(.not.nlsym2b) then
           ztop=1.0_fp
        else
           if(yion.ge.0.0_fp) then
              ztop=1.0_fp
           else
              ztop=-1.0_fp
           end if
        end if
        rmjflr=rmjion+zdrx*btiltc(5,ngc2)+ztop*zdrth*btiltc(4,ngc2)
        yflr=yion-ztop*zdrx*btiltc(4,ngc2)+zdrth*btiltc(5,ngc2)
     end if
     !
     !  evaluate the (xi,theta) flr position using the fast bilinear
     !  inverse map
     !
     !test        call nbimap(rmjflr,yflr,xiflr,thflr,ier)
     zphiarg=0.0_fp
     call xs_fxget(0.01_fp,zphiarg,rmjflr,yflr,xiflr,thflr,iwarn)
     !
     !-----------------------------------
     !  evaluate transp and beam dist fcn zone indices
     !
     !
     if(xiflr.gt.xminbm) then
        nflr2=0
        nflrx=nznblo+lcentr-1
        xinflrx=1.0_fp
        nflr=lep1
        xinflr=0.0_fp
     else
        nflrx=int(lcentr+xiflr*nzones)
        xinflrx=(xiflr-xiblo(nflrx))/(xiblo(nflrx+1)-xiblo(nflrx))
        xinflrx=max(0.0_fp,min(1.0_fp,xinflrx))
        !
        nflr=min(lep1,nflrx)
        if(nflr.le.ledge) then
           xinflr=xinflrx
        else
           xinflr=0.0_fp
        end if
        !
        ibr=int(lcentr+nznbmr*(xiflr/xminbm))
        inz0=nthzsm(ibr-1)
        inz=nthzsm(ibr)
        if(nlsym2b) then
           !     midplane symmetry assumed; theta ranging from -pi to pi
           nflr2=min(inz,int(inz0+1+(abs(thflr)*zpinv)*(inz-inz0)))
        else
           !     no midplane symmetry
           nflr2=min(inz,int(inz0+1+((thflr-thbdy0)*0.5_fp*zpinv) &
                *(inz-inz0)))
        end if
     end if
     !
     !-----------------------------------
     !  evaluate flr velocity vector
     !
     ephflr= &
          xzbeami*(phiprg(nflr)+xinflr*(phiprg(nflr+1)-phiprg(nflr)))
     !
     if(ephflr.ne.ephion) then
        vflr=sqrt(max(1.0e12_fp,(vion*vion+(ephion-ephflr)/v2toev)))
     else
        vflr=vion
     end if
     !
     !
     !
     if(iflag.eq.1) then
        !
        ! use the parallel velocity from map_larmor, this is the
        ! projection of the flr velocity vector on the guiding center
        ! magnetic field direction
        !
        ! when map_larmor fails and returns the classical larmor
        ! displacement, the energy will be wrong (which is why it
        ! failed) for the flr particle so scale the parallel velocity
        !  vflr/zv =approx= 1.04  has been seen
        !
        zv = sqrt(z_vpar*z_vpar+z_vperp*z_vperp)
        vplflr=z_vpar * vflr/zv
     else
        vplflr=rmjion*vpll/rmjflr   ! original method
     end if
     vplflr=max(-0.999999_fp*vflr,min(0.999999_fp*vflr,vplflr))
     !
     vppflr=sqrt(vflr*vflr-vplflr*vplflr)
     !
     xksflr=vplflr/vflr
     !
     !-----------------------------------
     !  no flr option ...
     !
  else
     rmjflr=rmjion
     yflr=yion
     xiflr=xion
     thflr=th
     nflr=ngc
     nflr2=ngc2
     !
     !  flr velocity = gc velocity
     ephflr=ephion
     vppflr=vperp
     vplflr=vpll
     vflr=vion
     xksflr=xksid
     !
  end if

  if (iflag.eq.0) then
     !     get alpha_v
     if(.not.nlfbmflr) then
        alpha_v=0.0_fp
     else
        alpha_v=0.0_fp
        if(nflr.ge.lcentr.and.nflr.le.ledge) then
           call nbi_aimingangle &
                (rmjflr,yflr,zangl,vflr,vppflr,zvec,alpha_v)
        end if
     end if
  end if
  !     find ion lgyro in there
  if(.not.nlfbmflr) then
     lgyro=1
  else
     call nbi_fbm_grid_indx(alpha_v,lgyro)
  end if
  alpha_gyro=alpha_v
  !
  !  flag if outside the plasma
  !
  if((nflr.eq.lep1).or.(nflr2.gt.nfbznsi)) then
     nlpout=.true.
  else
     nlpout=.false.
  end if
  !
  return
  !
end subroutine nbiflr
!******************** end file nbiflr.for ; group nbiflr ******************
