!******************** start file bparms.for ; group bparms *************
!---------------------------------------------------------------------
!  bparms
!  set up field description for beam code
!
!  rga oct00 - removed levgb = 1,2
!
!  rev  dmc feb 1987 -
!
!  pam stubberfield 12 july 1994 uda modifications
!
!   btiltc array definition modified
!   btiltc(j,iz):  zone iz on 2d beam grid
!    j=1:  mod(total field strength)
!    j=2:  mod(toroidal field strength) / mod(total field strength)
!    j=3:  mod(poloidal field strength) / mod(total field strength)
!        z=sqrt(dr/dth**2+dy/dth**2)
!    j=4:   (dr/dth)/z
!    j=5:   (dy/dth)/z
!
!    from these can be formed +/- the poloidal field vector
!    the correct sign is not knowable from available information
!    (the direction of the toroidal plasma current is not input to
!    transp)
!    btiltc(4,iz),btiltc(5,iz) also form the theta unit vector in (r,y)
!    tangent to the flux surface going around the flux surface in the
!    poloidal direction (ccw from outside midplane intercept)
!
!    btiltc(5,iz),-btiltc(4,iz) forms a unit vector in (r,y) in the
!    direction normal to the flux surface
!
!  rev dmc may 1990 --
!    for higher geometries, compute (de)enhancement factor
!      [ 1 / ( grad(xi)*rmajor ) ] / < 1 / ( grad(xi)*rmajor ) >
!    for the beam timestep controller, in each flux zone, including
!    beyond the plasma bdy.
!
!    this factor is empirically based; we have used q as a scaling
!    factor for the beam orbit timestep, this method introduces a
!    poloidal variation based on qlocal, where
!
!      qlocal = l*bt / (2*pi*r*bp)
!
!    where l is the distance around the flux surface of interest in
!    the poloidal plane.  we are normalizing this so we only care about
!    the part which varies vs. theta, which turns out to be
!
!           1 / ( grad(xi)*rmajor )
!
!------------------------------------------------------
!
subroutine bparms(icall,ierr)
  !
  !  icall=1 -- set up field description
  !  icall=2 -- evaluate field at midplane/flux surface intercepts
  !          -- find max field around each flux surface
  !          -- estimate max orbit width for full energy ptcls of
  !             each species near magnetic axis
  !
  use nbi_com
  use physconst_mod, only : aee_amp
  !
  implicit none
  integer, parameter :: inth=400
  real(fp) :: zra(inth),zza(inth),ztha(inth),zmodb(0:inth+1)
  real(fp) :: zjac(2,2),zbmaxmu,zbmaxml,zbsing,zthsing
  !
  integer :: icall,ierr,inumx,ixi,ic,iz,iz1,iz2,iersum
  integer :: i,j,jp1,jm1,ith,jx,izr,jxb1,jxb2,izz1,izz2
  real(fp) :: zr,zy,zd,zgradx,zbpmod,zth,zxi,ziotb,zdpsdx,zdxi
  real(fp) :: zh,zxb1,zxb2,zxc,zgc,zxibc,zthtot
  real(fp) :: zsum,znorm,zterm,zavg,zrmp,zymp,zbi
  real(fp) :: zrmaxm,zrminm,zxlar,zdxdr
  real(fp) :: zdum1,zdum2,zdum3,zdum4,zdum5
  !
  integer :: i2zp1
  !
  real(fp) :: zximin,zrlim,ze,zv,zvpll,zvprp,zb,zrlar,zrtes
  integer :: izrlim,isb
  !
  !  statement function:
  !
  real(fp) :: zf1d
  integer :: izth, izone
  real(fp) :: zdth
  !
  !--------------------------------------------------
  !  for orbit half width est. near axis
  !
  integer :: ireact,ibeam,ivpv,ivpvmax
  logical :: ilreact(mireact)  ! t for reactions that produce a given ion
  real(fp) :: zereact(mireact)   ! nominal fusion product birth energy
  real(fp) :: zengy,zepwr
  integer :: itrap,iuncon,ixiav
  real(fp) :: zxksav,zeiav,zphiav,zxhw,xjadep
  !
  !--------------------------------------------------
  !  for tau_bounce...
  !
  real(fp), parameter :: zcut = 1.0e-8_fp
  real(fp), parameter :: ztolb = 1.0e-4_fp
  !
  integer, parameter :: inum_tau = 7
  real(fp), dimension(inum_tau) :: zbarg = &
       (/ 0.24_fp,0.26_fp,0.49_fp,0.51_fp,0.74_fp,0.76_fp,0.99_fp /)
  real(fp), dimension(inum_tau,mj) :: zvtaub
  real(fp) :: zbmin,zbmax
  integer :: i1,i2,ierrx
  real(fp) :: zvtaub1,zvtaub2
  !
  real(fp) :: zomeg_exb(mj)
  real(fp) :: zdphdx,zq
  !
  !  in the moments arb. geometry code 1/mod(b) is evaluated as a
  !  function of flux coordinates (xi,theta)
  !
  !
  goto (10,20) icall
  !
  !-----------------------------------------------------------------
  !  set up field description
  !
10 continue
  !
  !  generalized code - define splines on fourier coefficients
  !
  call magcf(ierr)
  if(ierr.ne.0) return
  !
  !  dmc 13 june 1996 -- generate theta for a straight field line
  !   coordinate system, and the length along fieldline -- bicubic
  !   splines vs. xi(0 to 1) and theta(-pi to pi) periodic in the
  !   derivatives w.r.t. theta
  !
  !   no. of theta points:  nthspl
  !   stored at datspl(lsptheta...) the theta grid
  !
  !   stored at datspl(lspthets...):
  !     thets(xi,theta)  straight line (boozer) theta vs. transp theta
  !
  !   stored at datspl(lspphilb...):
  !     boozer nu(xi,theta) coordinate that allows evaluation of:
  !     philb(xi,theta)  length of field lines (start from -pi)
  !
  !   nb dmc -- not in current use; disabled unless nphipzc is set...
  !
  if(nphipzc.gt.0) then
     call gthets(ierr)
     if(ierr.ne.0) return
     !
     !  now on each surface find b field singularities and evaluate
     !  phi(p) vs. theta
     !
     call bphip(ierr)
     if(ierr.ne.0) return
     !
     !%        call chkspl     ! debug plots of numeric fits & splines
     !%        call eqdbg_plot ! debug plot of xplasma stored data
  end if
  !
  !  check that the poloidal flux spline fit is well behaved
  !
  ierr=0
  inumx=300
  do ixi=1,inumx
     zxi=ixi*(xbmbnd/real(inumx,fp))
     jx = int(zxi/dxi(lcentr,2)) + 1
     zh = zxi - real(jx-1,fp)*dxi(lcentr,2)
     ziotb  = zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx), &
          spliot(4,jx))
     zdpsdx = zxi*ziotb*tflux/xpi
     if(zdpsdx.le.0.0_fp) then
        write(nonlin,9925) zxi,zdpsdx

9925    format(' ?bparms:  at xi=',1pe11.4,' d(psi-pol)/d(xi)=', &
             1pe11.4)
        ierr=1
     end if

  end do
  if(ierr.eq.1) then

     write(nonlin,9926)
9926 format(/ &
          ' ?? bparms -- the poloidal flux spline fit has singularities'/ &
          '    the beam orbit code cannot cope with such field topology.')
     return
  end if
  !
  !  mag axis
  !
  if(nfast_map.eq.1) then
     call nbmomry(0.0_fp,0.0_fp,0,raxis,yaxis,zjac) ! use bilinear maps
  else
     call nbmomry_xs(0.0_fp,0.0_fp,0,raxis,yaxis,zjac) !use xplasma spline maps
  end if
  !
  !  field coords at each beam zone ctr -- * use magcf output *
  !
  zdxi=dxi(lcentr,2)
  do izr=1,nznbmr
     !
     jxb1=1+(izr-1)*nxskpb
     jxb2=jxb1+nxskpb
     !
     zxb1=(izr-1)*nxskpb*zdxi
     zxb2=zxb1+nxskpb*zdxi
     zxc=0.5_fp*(zxb1+zxb2)
     !
     zxibc=0.5_fp*(spliot(1,jxb1) + spliot(1,jxb2))
     zgc=0.5_fp*(splg(1,jxb1)+splg(1,jxb2))/(1.e4_fp*bzxr)
     !  poloidal flux derivative
     zdpsdx=zxibc*zxc*tflux/xpi
     !  zone indices
     izz2=izr+lcm1
     izz1=izz2-1
     iz1=nthzsm(izz1)+1
     iz2=nthzsm(izz2)
     !
     zthtot = thbdy1-thbdy0
     !
     do iz=iz1,iz2
        zth=thbdy0+ &
             (zthtot*(0.5_fp+(iz-iz1)))/(iz2-iz1+1)
        if(nfast_map.eq.1) then
           call nbmomry(zxc,zth,1,zr,zy,xbjaco) ! use bilinear maps
        else
           call nbmomry_xs(zxc,zth,1,zr,zy,xbjaco) !use xplasma spline maps
        end if
        !  1/determinant
        zd=1.0_fp/(xbjaco(1,1)*xbjaco(2,2)-xbjaco(1,2)*xbjaco(2,1))
        !  grad(xi) norm; grad(xi) vector = 1/det. (dy/dth, -dr/dth)
        zgradx=sqrt(xbjaco(2,2)**2+xbjaco(1,2)**2)*abs(zd)
        !  mod(poloidal field)
        zbpmod= abs(1.e4_fp*zdpsdx*zgradx/zr)
        btiltc(2,iz)=abs(bzxr*zgc/zr) ! toroidal component
        btiltc(3,iz)=zbpmod           ! poloidal component
        !  r component of poloidal field (proportional to dr/dth)
        btiltc(4,iz) = xbjaco(1,2) * zd / zgradx
        !  y component of poloidal field (proportional to dy/dth)
        btiltc(5,iz) = xbjaco(2,2) * zd / zgradx
        !  normalize; mod(b) total:
        btiltc(1,iz)=sqrt(btiltc(2,iz)**2+btiltc(3,iz)**2)
        do ic=2,3
           btiltc(ic,iz)=btiltc(ic,iz)/btiltc(1,iz)

        end do

     end do

  end do
  !
  !
  !----------------------------------
  !  dmc may 1990
  !
500 continue

  !      comment out resetr call and replace with f90 assing -jed 08/01
  !      call resetr(gxtabl,idim,1.0)

  gxtabl=1.0_fp
  !
  do j=lcentr,lcentr+nznblo-1
     !
     jp1=j+1
     zxi=0.5_fp*(xiblo(j)+xiblo(jp1))
     !
     !  flux zone average:
     !
     zsum=0.0_fp
     znorm=0.0_fp
     do ith=1,nthzons
        !
        zth=thzons(ith)
        !
        if(nfast_map.eq.1) then
           call nbmomry(zxi,zth,1,zr,zy,xbjaco) ! use bilinear maps
        else
           call nbmomry_xs(zxi,zth,1,zr,zy,xbjaco) !use xplasma spline maps
        end if
        !
        !  determinant
        zd=(xbjaco(1,1)*xbjaco(2,2)-xbjaco(1,2)*xbjaco(2,1))
        !  grad(xi) norm; grad(xi) vector = 1/det. (dy/dth, -dr/dth)
        zgradx=sqrt(xbjaco(2,2)**2+xbjaco(1,2)**2)/abs(zd)
        !
        zterm=1.0_fp/(zr*zgradx)
        gxtabl(ith,j)=zterm
        !
        zsum=zsum+zterm*zr*zd
        znorm=znorm+zr*zd
        !

     end do               ! theta
     zavg=zsum/znorm
     gxtabl(1:nthzons,j)=gxtabl(1:nthzons,j)/zavg

  end do                  ! xi
  !
  !----------------------------------
  !  dmc 10 apr 2010: compute omega_exb
  !  (toroidal rotation through bpol that causes a lorentz force that balances
  !  the electrostatic force from grad(phiprg)).
  !
  do j=lcentr,ledge
     jp1 = j+1
     zdphdx = (phiprg(jp1)-phiprg(j))/(xizb(jp1)-xizb(j))
     zq = (qgeo(j)+qgeo(jp1))*0.5_fp
     zomeg_exb(j)=-zq*xpi*zdphdx/(xizc(j)*tflux)
     if(j.gt.lcentr) then
        jm1=j-1
        omega_exb(j) = 0.5_fp*(zomeg_exb(jm1)+zomeg_exb(j))
     end if
  end do
  !
  omega_exb(1:lcentr) = omega_exb(lcp1) - &
       (omega_exb(lcp1+1)-omega_exb(lcp1))/3
  omega_exb(lep1:mj) = omega_exb(ledge) + &
       (omega_exb(ledge)-omega_exb(ledge-1))/3
  !----------------------------------
  !
  return
  !----------------------------------------------------------------
  !  icall=2...
  !
  !  code for finding bmidp, bminm and bmaxm and theta locations
  !  has been moved to nbi_getprofiles, with associated xplasma calls.
  !
20 continue
  !
  !-----------------------------
  !  compute "midplane" related quantities
  !
  !   dmc apr 2004 -- the "mnmx" arrays are really a compilation of info at
  !       the bmin and bmax locations on each flux surface-- there are
  !       two values (min,max) per non-singular surface, plus one value
  !       at the magnetic axis, hence the 2*nzones+1 grid.
  !
  i2zp1=nzp1+nzones
  !
  ympx=0.0_fp
  if(.not.nlsym2b) then
     ympx(lcentr)=yaxis
  end if
  rmjrmp(lcentr,1)=raxis
  rmnrmp(lcentr,1)=0.0_fp

  zh=0.0_fp
  rmnmx(nzp1)=raxis
  bmnmx(nzp1)=zf1d(zh,splg(1,1),splg(2,1),splg(3,1),splg(4,1))* &
       1.e-4_fp/raxis
  bpbmnmx(nzp1)=0.0_fp
  btbmnmx(nzp1)=1.0_fp
  psmnmx(nzp1)=0.0_fp
  !
  aratii(lcentr,1)=0.0_fp
  !
  do jm1=lcentr,lcentr+nznblo-1
     j=jm1+1
     do ith=1,inth-1
        zxi=xiblo(j)
        zth=(ith-1)*twopi/(inth-1)
        if(nfast_map.eq.1) then
           call nbmomry(zxi,zth,0,zra(ith),zza(ith),zjac) ! use bilinear maps
        else
           call nbmomry_xs(zxi,zth,0,zra(ith),zza(ith),zjac) !use xplasma spline maps
        end if
     end do
     zra(inth)=zra(1)
     zza(inth)=zza(1)
     !
     if(.not.nlsym2b) then
        call r8_plcentr(zra,zza,inth,zrmp,zymp)
        ympx(j)=zymp
     end if
     !
     zrmaxm=maxval(zra)
     zrminm=minval(zra)
     rmnrmp(j,1)=(zrmaxm-zrminm)/2
     rmjrmp(j,1)=(zrmaxm+zrminm)/2
     aratii(j,1)=rmnrmp(j,1)/rmjrmp(j,1)
     !
     rmnrmp(jm1,2)=(rmnrmp(jm1,1)+rmnrmp(j,1))/2
     rmjrmp(jm1,2)=(rmjrmp(jm1,1)+rmjrmp(j,1))/2
     aratii(jm1,2)=rmnrmp(jm1,2)/rmjrmp(jm1,2)
     !
     if(j.le.lep1) then
        i=j-lcentr+1
        !
        psmnmx(nzp1-(j-lcentr))=plflxg(j)
        psmnmx(nzp1+(j-lcentr))=plflxg(j)
        !
        bmnmx(nzp1-(j-lcentr))=bmaxm(j)
        bmnmx(nzp1+(j-lcentr))=bminm(j)
        !
        rmnmx(nzp1-(j-lcentr))=rbmaxm(j)
        rmnmx(nzp1+(j-lcentr))=rbminm(j)
        !
        zxc=xi(j,1)
        !
        zth=thbmaxm(j)
        if(nfast_map.eq.1) then
           call nbmomry(zxc,zth,1,zr,zy,xbjaco) ! use bilinear maps
        else
           call nbmomry_xs(zxc,zth,1,zr,zy,xbjaco) !use xplasma spline maps
        end if
        !  1/determinant
        zd=1.0_fp/(xbjaco(1,1)*xbjaco(2,2)-xbjaco(1,2)*xbjaco(2,1))
        !  grad(xi) norm; grad(xi) vector = 1/det. (dy/dth, -dr/dth)
        zgradx=sqrt(xbjaco(2,2)**2+xbjaco(1,2)**2)*abs(zd)
        !
        btbmnmx(nzp1-(j-lcentr)) = &
             zf1d(zh,splg(1,i),splg(2,i),splg(3,i),splg(4,i))* &
             1.e-4_fp/(rbmaxm(j)*bmaxm(j))
        btbmnmx(nzp1-(j-lcentr)) = &
             min(1.0_fp, btbmnmx(nzp1-(j-lcentr))) ! rga, allow some interp error
        bpbmnmx(nzp1-(j-lcentr)) = &
             sqrt(1.0_fp-btbmnmx(nzp1-(j-lcentr))**2)

        if(j.eq.lep1) then
           drmnmx(nzp1-(j-lcentr)) = 1.0_fp/ &
                (zgradx/dxi(ledge,2))
        else
           drmnmx(nzp1-(j-lcentr)) = 1.0_fp/ &
                (2*zgradx/(dxi(j,2)+dxi(jm1,2)))
        end if
        !
        zth=thbminm(j)
        if(nfast_map.eq.1) then
           call nbmomry(zxc,zth,1,zr,zy,xbjaco) ! use bilinear maps
        else
           call nbmomry_xs(zxc,zth,1,zr,zy,xbjaco) !use xplasma spline maps
        end if
        !  1/determinant
        zd=1.0_fp/(xbjaco(1,1)*xbjaco(2,2)-xbjaco(1,2)*xbjaco(2,1))
        !  grad(xi) norm; grad(xi) vector = 1/det. (dy/dth, -dr/dth)
        zgradx=sqrt(xbjaco(2,2)**2+xbjaco(1,2)**2)*abs(zd)
        !
        btbmnmx(nzp1+(j-lcentr)) = &
             zf1d(zh,splg(1,i),splg(2,i),splg(3,i),splg(4,i))* &
             1.e-4_fp/(rbminm(j)*bminm(j))
        btbmnmx(nzp1+(j-lcentr)) = &
             min(1.0_fp, btbmnmx(nzp1+(j-lcentr))) ! rga
        bpbmnmx(nzp1+(j-lcentr)) = &
             sqrt(1.0_fp-btbmnmx(nzp1+(j-lcentr))**2)
        if(j.eq.lep1) then
           drmnmx(nzp1+(j-lcentr)) = 1.0_fp/ &
                (zgradx/dxi(ledge,2))
        else
           drmnmx(nzp1+(j-lcentr)) = 1.0_fp/ &
                (2*zgradx/(dxi(j,2)+dxi(jm1,2)))
        end if
     end if
  end do
  drmnmx(nzp1)=(drmnmx(nzp1-1)+drmnmx(nzp1+1))/2
  !
  !---------------/start dmc 1-dec-1998 change
  !  dmc -- 1-dec-1998 -- check reasonableness of xiliml
  !  control settings  (moved here from nbstart 3-jan 2002 dmc).
  !
  zximin=1.0_fp
  izrlim=2*nzones+1
  zrlim=rmnmx(izrlim)
  !
  do isb=1,nsbeam
     lsbeam=isb
     if(.not.nlfbmfpp(lsbeam).and.(.not.nlheavy(lsbeam))) then
        call inispc
        ze=2.0_fp*einjs(lsbeam) ! test energy
        zv=sqrt(ze/v2toev)   ! test perp velocity

        do j=lcentr,ledge
           zxi=xi(j,1)
           if(zxi.lt.zximin) then
              izr=nzones+1+j-lcentr
              zb=bmnmx(izr)         ! local b on outer half midplane
              zr=rmnmx(izr)
              zrlar = zv*abeam/(aee_amp*xzbeami*zb)
              zrtes = zr+zrlar
              if(zrtes.gt.zrlim) then
                 zximin=zxi
              end if
           end if
        end do

     end if
  end do
  !
  if(zximin.lt.xiliml) then
     write(nonlin,4901) xiliml,zximin
4901 format( &
          ' **nbstart:  due to large larmour radii, xiliml is reduced:'/ &
          '             from ',1pe11.4,' to ',1pe11.4)
     xiliml=zximin
  end if
  !
  !--------------------------------------------------
  !  calculate tau_bounce estimates for trapped orbits
  !  used to force early update of goose in case there are particles
  !  whose orbits do not reach the midplane and so cannot receive the
  !  normal midplane crossing goose updates.
  !
  do j=lcp1,lep1
     call eq_tau_bounce(xizb_r8(j),inum_tau,zbarg,ztolb,zcut, &
          zbmin,zbmax,zvtaub(1:inum_tau,j),ierr)
     if(ierr.ne.0) then
        write(nonlin,*) ' ?bparms: eq_tau_bounce: ierr=',ierr
        exit
     end if
  end do
  if(ierr.ne.0) return
  !
  zvtaub(1:inum_tau,lcentr)=zvtaub(1:inum_tau,lcp1)
  !
  !  choose lesser of paired values to avoid possible singularity
  !  convert v*tau_bounce to cm
  !
  do j=lcentr,ledge
     vtaub(4,j)=100.0_fp* &
          max(zvtaub(inum_tau,j),zvtaub(inum_tau,j+1))
     do i=1,3
        i1=2*(i-1)+1
        i2=2*i
        zvtaub1=min(zvtaub(i1,j),zvtaub(i2,j))
        zvtaub2=min(zvtaub(i1,j+1),zvtaub(i2,j+1))
        vtaub(i,j)=100.0_fp*max(zvtaub1,zvtaub2)
     end do
     vtaubmin(j)=minval(vtaub(1:4,j))
  end do
  vtaub(1:4,lep1)=vtaub(1:4,ledge)  ! for orbits beyond bdy...
  vtaubmin(lep1)=vtaubmin(ledge)
  !
  !--------------------------------------------------
  !  estimate max orbit widths near the magnetic axis
  !
  zxi=0.2_fp
  zth=0.0_fp
  if(nfast_map.eq.1) then
     call nbmomry(zxi,zth,0,zr,zy,zjac) ! use bilinear maps
  else
     call nbmomry_xs(zxi,zth,0,zr,zy,zjac) !use xplasma spline maps
  end if
  zdxdr=zxi/abs(zr-raxis)   ! delta(xi)/delta(r)
  !
  !  loop over species
  !
  do isb=1,nsbeam
     lsbeam=isb
     xhwdep0(lsbeam)=0.0_fp
     if(.not.nlfbmfpp(lsbeam).and.(.not.nlheavy(lsbeam))) then
        call inispc
        !
        !  set nominal energy
        !
        if(nlfprod(lsbeam)) then
           call fuslis(lsbeam,ilreact,zereact)
           zengy=50000.0_fp
           do ireact=1,mireact
              if(ilreact(ireact)) zengy=max(zengy,zereact(ireact))
           end do
        else
           zengy=0.0_fp
           zepwr=0.0_fp
           do ibeam=1,nbeam
              zengy=zengy+max(1.0_fp,pinja(ibeam))*einja(ibeam)
              zepwr=zepwr+max(1.0_fp,pinja(ibeam))
           end do
           zengy=max(50000.0_fp,zengy/max(1.0_fp,zepwr))
        end if
        !
        !  set v & loop over vpll/v; call xjadep; find max orbit half width estimate.
        !
        zb=bmnmx(nzp1)      ! b on axis, t
        zv=sqrt(zengy/v2toev)
        ivpvmax=50
        do ivpv=1,ivpvmax-1
           zvpll=ivpv*zv/ivpvmax
           zvprp=sqrt(zv*zv-zvpll*zvpll)
           zdum1=xjadep(zxi,zth,zr,zy,zv,zvpll, &
                itrap,iuncon,zxksav,zeiav,zphiav,ixiav,zxhw)
           zrlar=abeams(lsbeam)*zvprp/ &
                (aee_amp*xzbeams(lsbeam)*zb)
           zxlar=zrlar*zdxdr
           xhwdep0(lsbeam)=max((zxhw+zxlar),xhwdep0(lsbeam))
        end do
     end if
  end do
  !
  !--------------------------------------------------
  !
  return
end subroutine bparms
!******************** end file bparms.for ; group bparms ***************
