!******************** start file fokker.for ; group fokker *************
!***********************************************************
!
!	this subroutine sets up fokker-planck coefficients
!
subroutine nbi_fokker(ierr)
  !
  !
  !  evaluate fokker-planck coefficients
  !
  !  dmc jan 1992 -- support heavy beam ions injected by heavy neutral
  !  beams (neon,krypton,argon,xenon beams experiment -- tftr 1992)
  !
  !  upgrade dmc sept 1990 -- supporting fusion product ion slowing down
  !  for each of the several fast ion species, compute fokker planck
  !  coefficients at several energy levels -- taking into account
  !  variation in log(lambda) during slowing down process.
  !
  !  also taking into account possibility that z of fast ion .ne. 1
  !
  !  also the old bbslo hack is deleted.
  !
  !	modified by m.c.zarnstorff	april 1991
  !	to call general routine to get coulomb logarithm
  !
  use nbi_com
  use physconst_mod, only : evptocm_sec,me_kg,mp_kg
  !
  !============
  implicit none
  integer :: ig,isb,iengy,iheavy,j,izmax,iz
  !============
  integer, parameter :: migmax=50
  real(fp),dimension (migmax) :: zallhy
  real(fp),dimension (migmax+2) :: ztemp,zdens,za,zz,zall
  real(fp) :: zalle,zallim,zeinj
  real(fp) :: zabeam,zzbeam,zb2,zgam,zenll,zsqrte,zsqrne,zbz,zcfall
  real(fp) :: zdfall,zef,zdf,zn,zzai,zvc,zbi2,zcffac,zdffac,zcfa
  real(fp) :: zdfa,zefa,zbi
  !============
  integer, intent(out) :: ierr   ! status code (0=ok)
  !
  logical :: ilstix
  !
  data ilstix/.false./
  !-------------------------------------
  !   jan 1995 -- anomalous diffusivity setup code moved to new subroutine,
  !    gdifb.
  !   nov 2001 -- gdifb moved out of beam code -- driver routine now.
  !
  !xx      ztime=tbm1
  !xx      call gdifb(ztime)
  !
  !-------------------------------------
  !
  !  factor for energy indexing of fokker planck coefficients:
  !     iengy=1:  full energy
  !     iengy=2:  fac_e * full energy
  !     iengy=3:  fac_e**2 * full energy
  !     iengy=4:  fac_e**3 * full energy
  !
  ierr=0
  !
  if(migmax.lt.mig) then
     write(nonlin,*) ' mig = ',mig,' migmax = ',migmax
     write(nonlin,*) ' ?fokker: local parameter migmax too small!'
     ierr=1
     return
  end if
  !
  !    save the charge and mass of each species for calculating coloumb logs
  !
  za(1) = me_kg/mp_kg
  zz(1) = -1.0_fp

  do ig=1,ng
     za(2+ig) = aplasm(ig)
     zz(2+ig) = backz(ig)

  end do

  fac_e=0.25_fp
  !
  do isb=1,nsbeam
     !
     zeinj=einjs(isb)
     zabeam=abeams(isb)
     zzbeam = xzbeams(isb)
     !
     zb2=xzbeams(isb)**2
     !
     zgam=2.39e11_fp/(zabeam*zabeam)
     zenll=zeinj/fac_e
     !
     do iengy=1,4
        !
        zenll=zenll*fac_e
        !
        if(nlheavy(isb).and.(iengy.eq.1)) then
           iheavy=1
        else
           iheavy=0
        end if
        !
        do j=lcentr,ledge
           za(2) = aimpj(j)
           zz(2) = xzimpj(j)
           !
           if(iengy.eq.1) then
              vmin(j,isb)=evptocm_sec* &
                   sqrt(1.5_fp*max(ti(j,2),te(j,2))/zabeam)
              vmin3(j,isb)=vmin(j,isb)**3
           end if
           !
           zsqrte=sqrt(te(j,2))
           zsqrne=sqrt(rhoel(j,2))
           !
           !     save temperatures and densities in mks-kev for calculating
           !     coloumb logarithms
           !
           ztemp(1) = te(j,2)/1000._fp
           zdens(1) = rhoel(j,2)*1.e6_fp

           ztemp(2) = ti(j,2)/1000._fp
           zdens(2) = rhi(j,2)*1.e6_fp

           do ig=1,ng
              ztemp(ig+2) = ti(j,2)/1000._fp
              zdens(ig+2) = rhob(j,ig,2)*1.e6_fp

           end do

           zbz=bzxr*fbz(j,2)*grizc(j) ! approx. b, tesla (bz close enuf)

           !  compute coulomb logs ...
           call r8_coulog(zenll/1000._fp, zabeam, zzbeam, &
                zdens, ztemp, za, zz, zall, &
                ng+2, zbz, ilstix )

           zalle=zall(1)
           zallim=zall(2)
           zallhy=zall(3)

           if(iengy.eq.1) cloge(j) = zalle

           zcfall=zb2*zalle
           cfa(j,isb,iengy)=zcfall*6.62e-21_fp*rhoel(j,2)*zabeam/ &
                (te(j,2)*sqrt(te(j,2)))*zgam
           !
           dfa(j,isb,iengy)=zb2*zabeam*zgam*rhimom(j,2,5)*zallim
           zdfall=zb2*zallim
           !
           efa(j,isb,iengy)=zb2*zgam*0.5_fp*rhimom(j,2,2)*zallim
           !
           zef=efa(j,isb,iengy)
           zdf=dfa(j,isb,iengy)
           !
           do ig=1,ng
              zn=rhob(j,ig,2)
              dfa(j,isb,iengy)=dfa(j,isb,iengy)+ &
                   zb2*zabeam*zgam*zn*backz(ig)*backz(ig)* &
                   zallhy(ig)/aplasm(ig)
              efa(j,isb,iengy)=efa(j,isb,iengy)+ &
                   zb2*zgam*0.5_fp*zn*backz(ig)*backz(ig)*zallhy(ig)
              zdfall=zdfall+zb2*zallhy(ig)

           end do
           !
           if(iengy.eq.1) then
              zzai=dfa(j,isb,iengy)*zabeam/ &
                   (zb2*rhoel(j,2)*zalle*2.39e11_fp)
              !
              zvc=5.33e6_fp*zsqrte*zzai**0.3333333_fp
              vcrit(j,isb)=zvc
              vcrit3(j,isb)=zvc**3
           end if
           !-------------------------------------------------------------
           !  dmc jan 1992 -- correction factors for use with heavy beam ions
           !   which vary in charge as their electrons are stripped away ...
           !
           !  cfa, dfa, and efa were computed for fully stripped ions
           !
           !  see subroutine hvsigma which evaluated log(lambda) tables --
           !  note that the log lambda arrays are reused!
           !
           if((iengy.eq.1).and.(iheavy.eq.1)) then
              izmax=xzbeams(isb)+0.5_fp
              do iz=1,izmax
                 zbi=iz
                 zbi2=zbi*zbi
                 !
                 zcffac=zbi2*allame(j,iz,isb)/zcfall
                 allame(j,iz,isb)=zcffac
                 !
                 zdffac=zbi2*allamx(j,iz,isb)
                 do ig=1,ng
                    zdffac=zdffac+zbi2*allamh(ig,j,iz,isb)
                 end do
                 zdffac=zdffac/zdfall
                 allamx(j,iz,isb)=zdffac
                 !
              end do
           end if
           !-------------------------------------------------------------
           !
           !  ... still in j loop ...
           !
           !  ** jun 5 80 d. mc cune calculate fraction of pbi going
           !         to impurity species
           !
           if(iengy.eq.1) then
              paxfrc(j,isb)=zef/efa(j,isb,iengy)
              pbxfrc(j,isb)=zdf/dfa(j,isb,iengy)
           end if
           !

        end do

     end do

  end do
  !
  !-----------------------------------------------------------------
  !  dmc mar 1995 -- renormalization option
  !
  do iengy=1,4
     do isb=1,nsbeam
        do j=lcentr,ledge
           if(nlfprod(isb)) then
              zcfa=max(epslon,xcfafus)
              zdfa=max(epslon,xdfafus)
              zefa=max(epslon,xefafus)
           else
              zcfa=max(epslon,xcfanbi)
              zdfa=max(epslon,xdfanbi)
              zefa=max(epslon,xefanbi)
           end if
           cfa(j,isb,iengy)=zcfa*cfa(j,isb,iengy)
           dfa(j,isb,iengy)=zdfa*dfa(j,isb,iengy)
           efa(j,isb,iengy)=zefa*efa(j,isb,iengy)
        end do
     end do
  end do
  !-----------------------------------------------------------------
  !  dmc sept 1990 -- driven current shielding factor calculation has
  !  been moved to new subroutine jbshld, where the zb dependence has
  !  also bee added ...
  !
  !  dmc jan 1995 -- anomalous diffusivity setup stuff moved to new
  !  routine, gdifb
  !
  !-------------------------------------------------------------
  !
  return
end subroutine nbi_fokker
!
!******************** end file fokker.for ; group fokker ***************
