!-----------------------------------------------------------------------
!  jbshld -- classical/neoclassical beam current shielding factor
!
!   taken out of fokker and put here dmc sept 1990
!
subroutine jbshld
  !
  use nbi_com
  !
  !============
  implicit none
  integer :: j
  !============
  real(fp) :: zrzon,zrboun,zrbounp,zrmajor,zcurr,zvstae,zdelta,zd1m
  real(fp) :: zft,zef,zd32,za11,zb11,zc11,zd11,zk11,za12,zb12,zc12
  real(fp) :: zd12,zk12
  real(fp) :: zvstaeh,zxt,zl31e,zzf1
  !============
  logical ilspiz
  !
  !-------------------------------------------------------------
  !  d. mc cune impliment neoclassical beam current option
  !    collisionality, trapping fraction formulae from "resis.for"
  !    a,b, k0 coefficients from
  !
  !  neoclassically shielded beam driven current: options:
  !
  !  nmcurb=0 no fast ion driven current -- beam current drive suppressed.
  !    (logic to suppress elsewhere in nubeam)
  !
  !  nmcurb=1 (neoclassical shielding vsn #1) references:
  ! [1] hirshman, phys. fluids vol 21 no. 8 aug 1978, p.1295.
  ! [2] hirshman, phys. fluids vol 21 no. 2 feb. 1978, p.224.
  ! [3] s.p.hirshman, d.j.sigmar, nuclear fusion, vol.21, no.9(1981), p.1079
  !
  !  nmcurb=2 (classical okhawa current shielding): use xjbfac = (1.0 - zb/zeff)
  !     (no neoclassical correction)
  !
  !  nmcurb=3 (neoclassical shielding vsn #2) reference:
  !     y.r. lin-liu, f.l. hinton, phys. plasmas 4 (1997) 4179
  !
  !  study by marina gorelenkova (mgorelen@pppl.gov) in response to
  !  physicists' queries -- mg april 27,2011: verified work of michael kraus
  !  (link to his write up, below).
  !
  !  in summary:
  !
  !     nmcurb=1 is good in banana and plateau regimes, low to middling
  !     collisionality, and limited aspect ratio;
  !
  !     nmcurb=3 is only good in banana regime (low collisionality) but
  !     there is no validity limit on aspect ratio.
  !
  !     current drive algorithms by hirshman and sigmar [1] are used for
  !     nmcurb=1
  !
  !     analytical function of nu_e^* and zeff (eq.40,[1]) is programmed
  !     in jbshld if nmcurb=1, where nu_e^* is the normalized electron
  !     collisionality (eq.37,[1]).
  !
  !     the hirshman algorithm is valid in the banana as well as plateau
  !     regimes, and for the range
  !
  !     10^(-2) < (zeff * nu_e^*) <10^2 * eps^(-3/2) ... (p.1300,[1])
  !
  !     eps=r/r  is the inverse aspect ratio, midplane definition
  !
  !     the formula nmcurb=1 formula is only good for low aspect ratio.
  !     hirshman's estimation is 0.001 =< eps =< 0.15  with an error +/- 6%
  !     (1st page of [1]).
  !
  !     more details on derivation of the transp algorithm for nmcurb=1 can
  !     be found in michael kraus (michael.kraus@ipp.mpg.de) paper
  !     "transp current drive algorithm" -- nc references verified and
  !     equations checked by marina 27 apr 2011
  !
  !     for michael kraus's write up:
  !     see http://w3.pppl.gov/transp -- papers & documents section
  !       "transp/nubeam current drive algorithm" -- pdf:
  !       http://w3.pppl.gov/transp/papers/kraus_en.pdf

  !     nmcurb  = 4, shielding calculation by honda m., kikuchi m., azumi m.
  !     "collisionality dependence on shielding factor of beam driven current"
  !     nuclear fusion, submitted by october 2011. martix inversion method are used
  !     to include collisionality dependence on a shielding factor.
  !     s. kaye (11/2011)
  !     -----------------------------
  !     added option nmcurb=4 from model fit in honda, 2011 article (submitted
  !     nuclear fusion), eqs. 16 - 18. the model fit takes into account the
  !     collisionality dependence of the shielding factor at arbitrary aspect ratio.
  !     this model fit is easier to implement than the matrix inversion method
  !     described in the paper.
  !

  !
  !----------------

  ilspiz=.false.
  if(nmcurb.eq.2) then
     ilspiz=.true.
  else if(nmcurb.eq.3) then
     ilspiz=.false.
     call nc_lin_liu        ! n.c lin-liu model phys.plasmas 4(11),p.4179, 1997
     return
  end if
  !
  do j=lcentr,ledge
     xjbfac(j)=(1.0_fp-xzbeami/xzeffc(j,2))
     !
     !     if resistivity is classical, use classical beam current
     !     (nlspiz=.true. ==> classical, spitzer resistivity s.r. resis)
     !     *** mod dmc oct 1989 ***
     !     can force classical/n.c. beam driven current independently
     !     see namelist control nmcurb, and defn. of ilspiz above ***
     !
     if(ilspiz) then  !nmcurb=2
         xjbfac(j)=max(1.e-5_fp,xjbfac(j))
         xjbfacs(j,lsbeam)=xjbfac(j) ! multi-species
         cycle
     endif
     !
     zrzon=rmnrmp(j,2)
     zrboun=rmnrmp(j,1)
     zrbounp=rmnrmp(j+1,1)
     zrmajor=rmjrmp(j,2)
     !
     !     neoclassical jb
     zcurr=2.5e4_fp*bzxr/zrmajor**2 &
          *zrzon*(zrboun/qgeo(j)+zrbounp/qgeo(j+1))
     !     collisionality nu*e
     zvstae=3.46e-9_fp*rhoel(j,2)*cloge(j)*bzxr* &
          sqrt(zrmajor*zrzon)/(zcurr*te(j,2)*te(j,2))
     zdelta=zrzon/zrmajor
     zvstaeh=zvstae*xzeffc(j,2)
     zd1m=1._fp-zdelta
     !     trapping fraction
     zft=1._fp-zd1m*zd1m/(sqrt(1._fp-zdelta*zdelta)* &
          (1._fp+1.46_fp*sqrt(zdelta)))

     if (nmcurb.eq.4) then

        zxt=zft/(1._fp+(1._fp-0.1_fp*zft)*sqrt(zvstaeh)+0.5_fp* &
             (1._fp-zft)*zvstaeh/xzeffc(j,2))
        zzf1=xzeffc(j,2)+1._fp
        zl31e=(1._fp+1.4_fp/zzf1)*zxt - (1.9_fp/zzf1)*zxt*zxt + &
             (0.3_fp/zzf1)*zxt*zxt*zxt + (0.2_fp/zzf1)* &
             zxt*zxt*zxt*zxt
        xjbfac(j)=1._fp-(xzbeami/xzeffc(j,2))*(1._fp-zl31e)

     else !nmcurb=1, or nmcurb >4

        !     k11 coefficient
        zef=xzeffc(j,2)
        zd32=zdelta*sqrt(zdelta)
        za11=3.44_fp*zef+(.52_fp-.42_fp*zef)/(1._fp+1.35_fp*zef)
        zb11=.56_fp+.96_fp*zef
        zc11=.25_fp*zef+(.14_fp+.55_fp*zef)/(1._fp+5.0_fp*zef)
        zd11=.51_fp*zef+(.7_fp+.78_fp*zef)/(1._fp+zef)
        zk11=(.53_fp+zef)/ &
             ((1._fp+sqrt(za11*zvstae)+zb11*zvstae)* &
             (1._fp+sqrt(zc11*zvstae*zd32)+zd11*zvstae*zd32))
        !     k12 coefficient
        za12=.31_fp*zef+(.1_fp+.084_fp*zef)/(1._fp+1.3_fp*zef)
        zb12=.26_fp+.35_fp*zef
        zc12=.081_fp*zef+(.072_fp+.15_fp*zef)/(1._fp+3._fp*zef)
        zd12=.29_fp*zef+(.42_fp+.62_fp*zef)/(1._fp+1.42_fp*zef)
        zk12=(.71_fp+zef)/ &
             ((1._fp+sqrt(za12*zvstae)+zb12*zvstae)* &
             (1._fp+sqrt(zc12*zvstae*zd32)+zd12*zvstae*zd32))
        !     neoclassical beam current correction
        xjbfac(j)=xjbfac(j)+xzbeami/zef*zft* &
             (1.5_fp*zef*(zk12-2.5_fp*zk11)+ &
             (1.414214_fp+3.25_fp*zef)*zk11) &
             /(zef*(1.414214_fp+3.25_fp*zef)-2.25_fp*zef*zef)
     end if
     !
     xjbfac(j)=max(1.e-5_fp,xjbfac(j))
     xjbfacs(j,lsbeam)=xjbfac(j) ! multi-species
     !
  end do
  !
  !--------------------------------------------------------
  !
contains
  subroutine nc_lin_liu
    use iso_c_binding, only: fp => c_double
    ! neoclassical beam current shielding factor, nmcurb=3.
    ! lin-liu, hinton analytic expressions for the bootstap coefficient valid for
    ! arbitrary aspect ratio tokamaksin the banana regime and in the limit
    ! that electron thermal velocity greatly exceeds the fast ion velocity
    ! for electrical current driven by neutral beams.
    ! references:
    !------------
    ! 1. n.c lin-liu, hinton phys.plasmas 4(11),p.4179, 1997
    ! 2. n.c lin_liu,miller phys. plasmas 2(5), p.1666, 1995
    !
    !  the lin-liu, hinton paper identifies the shielding factor in the form
    !  [1-(zion/zeff)*(1-g)]
    !
    !  where g=l_31 and approximated by
    !
    !   g = x*[(0.754+2.21*zeff+zeff**2) + x*(0.348+1.243*zeff+zeff**2)] / d
    !
    !   d = 1.414*zeff+zeff**2 + x*(0.754+2.657*zeff+2*zeff**2) +
    !       x**2*(0.348+1.243*zeff+zeff**2)
    !
    !   x = ft/(1-ft), ft is the effective trapped ptcl.fraction
    !
    !
    !  (1-g)/zeff={1.414+zeff-ft*(2.381+zeff)+ft**2*0.967)}/
    !             {1.414*zeff+zeff**2+ft*(0.754-0.171*zeff)-ft**2*0.406}
    !
    !
    real(fp) :: zft,zftplus,zftminus ! the effective trapped ptcl.fraction zft,
    !  upper limit - zftplus, lower limit - zftminus
    ! zft = 3/4*zftplu+1/4*zftminus [2]

    integer :: jzb     ! =1:  compute at inner bdy of zone j-zone; avg btw 2 neighboring zones
    ! =2:  compute at center of zone j-zone
    real(fp) :: zef2 !zef^2
    real(fp) :: zft2 !zft^2
    real(fp) :: one_g_zef  ! simlified coefficient (1-g)/zeff from eqs(21,22) [1]
    real(fp) :: ag,cg !local
    jzb = 2
    do j=lcentr,ledge !loop over zone centers
       zef=xzeffc(j,2)
       zef2=zef*zef
       call nbi_ncift(j,jzb,zft,zftminus,zftplus) !find the effective trapped ptcl.fraction zft
       zft2=zft*zft
       cg = 1.414_fp+zef
       ag = 0.754_fp-0.171_fp*zef
       one_g_zef =(cg-zft*(2.381+zef)+zft2*0.967_fp)/ &
            (cg*zef+ag*zft-0.406_fp*zft2) !simplified (1-g)/zef  see eqs. (21,22) [1]

       xjbfac(j)=1.0_fp-xzbeami*one_g_zef
       xjbfac(j)=max(1.e-5_fp,xjbfac(j))
       xjbfacs(j,lsbeam)=xjbfac(j) ! multi-species
    end do
  end subroutine nc_lin_liu


end subroutine jbshld
