!-------------------------------------------------------------------
!  pre_sawnbi - evaluate effect of kadomtsev sawtooth on fast ion
!  population.  store results "on the side", optionally available for
!  later use.
!
!------------------------
!  dmc nov 2005 -- for parallelization, this single routine has been
!  broken up into serial & parallel parts:
!
!  (serial) pre_sawnbi_setup  -- initialization (call from nbstart)
!  (parallel) pre_sawnbi(ind) -- application of operator to ptcl ind
!                                (call from orball)
!  (serial) pre_sawnbi_finish -- completion (call from nbfinish)
!------------------------
!
subroutine pre_sawnbi_setup
  use physconst_mod, only: zero
  !
  !-------------------------------------------------------------------
  !  f90 modules...
  !
  use xstraln_calls
  use nbi_com
  use nbi_random_pll
  !
  !-------------------------------------------------------------------c
  implicit none
  !
  !-------------------------------------------------------------------
  !  set up the sawtooth calculation (serial code called from nbstart)
  !  (a) initialize sums; (b) determine which species participate;
  !  (c) call kdsaw_init to set up the mixing module.
  !
  !-------------------------------------------------------------------
  real(fp) :: ziotb(nzp1)                  ! temp iota(bar) profile
  integer :: iuse                     ! # of species in sawtooth
  integer :: isb,ierr
  real(fp) :: zxinv,zxmin
  !-------------------------------------------------------------------
  !
  !  clear sawtooth-associated sums and flags
  !
  wsum_kdspre = zero
  wsum_kdspost = zero
  !
  dxsm_kdsaw = zero
  dxwt_kdsaw = zero
  !
  kmix_kdsaw = 0
  xmix_kdsaw = zero
  !
  flag_cursaw = .false.
  !
  !  for now, zero the profile sums as well...
  !
  bdenssw = zero
  udensplw = zero
  udensppw = zero
  ebapl2sw = zero
  ebapp2sw = zero
  omegabsw = zero
  !
  !  determine what species if any are involved in sawtoothing...
  !
  iuse=0
  !
  do isb=1,nsfast
     if(.not.nlfbmfpp(isb)) then
        !  not an fpp specie.
        if(nlsawf.and.nlfprod(isb)) then
           !  fusion product...
           if(xswfrac_allfast*xswfrac_beam.gt.zero) then
              flag_cursaw(isb)=.true.
           end if
        end if
        if(nlsawb.and.(isb.le.nsbeam).and.(.not.nlfprod(isb))) then
           !  beam specie...
           if(xswfrac_allfast*xswfrac_fusn.gt.zero) then
              flag_cursaw(isb)=.true.
           end if
        end if
        if(flag_cursaw(isb)) then
           iuse=iuse+1
        end if
     end if
  end do
  !
  if(iuse.eq.0) return      ! no sawteeth
  !
  !-------------------------------------------------------------------
  !
  ziotb=zero
  ziotb(1:nzp1)=1.0_fp/qgeo(lcentr:lep1)
  !
  !  init the sawtooth (kdsaw)
  !
  call kdsaw_pre_init(fporcelli,nmix_kdsaw)
  call kdsaw_init(nonlin,nzones,xizb(lcentr:lep1),ziotb, &
       dvolzc(lcentr:ledge),kmix_kdsaw,ierr)
  !
  if((kmix_kdsaw.eq.0).or.(ierr.ne.0)) then
     if(ierr.ne.0) then
        write(nonlin,*) &
             ' %pre_sawnbi: kdsaw error, mixing suppressed.'
     end if
     kmix_kdsaw=0
     flag_cursaw=.false.
     return
  end if
  !
  kmix_kdsaw = kmix_kdsaw+lcentr-1 ! index to last mixed zone
  !
  !  get mixing radius...
  !
  call kdsaw_xmix(zxmin,zxinv,xmix_kdsaw)
  !
  return
end subroutine pre_sawnbi_setup
!-------------------------------------------------------------------
!
subroutine pre_sawnbi(ind)
  use iso_c_binding, only: fp => c_double
  !
  !-------------------------------------------------------------------
  !  f90 modules...
  !
  use physconst_mod, only: zero,rad2deg
  use xstraln_calls
  use nbi_com
  use nbi_random_pll
  !
  !
  !-------------------------------------------------------------------
  !
  !  mod dmc -- multispecies beam support
  !   fusion product support.  note bdens sums do not include fusion products.
  !   note also:  fpp can be used for beams instead of mc.  if so, it is up
  !   to the sawfppic subroutine to do the sawtooth mixing, and the
  !   algorithm is not quite the same (the fpp dist fcn is simply flattened
  !   in the entire mixing region).  sawfppic can also contribute to the
  !   various common sums dnnbsw etc.
  !
  !  mod dmc nov 2004 -- partial mixing fractions now supported,
  !   with different net mixing fractions allowed btw beam and fusion
  !   product fast ions
  !
  !  mod dmc nov 2005 -- breakup of routine for parallelization; removal
  !   of explicit & implied loops over the entire mc particle list.
  !
  !-------------------------------------------------------------------
  !  dmc feb 1986
  !
  !  the monte carlo beam population is rearranged in accordance with
  !  the kadomtsev mixing factors (cf sawmix subroutine) or porcelli
  !  mixing.
  !
  !  each m.c. ion in the mixing region is moved from its original
  !  pre-sawtooth source zone to one of the possible post-sawtooth
  !  target zones (a random choice is made to choose the target zone)
  !  the theta coordinate of the moved ion is randomized.
  !
  !  the perp velocity of the moved ion is modified for mu conservation
  !  (vperp**2/b a constant; the magnitude b of the magnetic field is
  !  a function of flux zone and theta);
  !  the parallel velocity is conserved
  !
  !  thus both the fast ion speed vion = sqrt(vperp**2+vpll**2)
  !  and the fast ion pitch vpll/vion are modified
  !
  !  modified copies of the monte carlo fast ion list arrays are created:
  !  (cf nbspec.dat):
  !      njavay -->   njavay_saw
  !      rmjionay --> rmjionay_saw
  !      xzionay -->  xzionay_saw
  !      tfluxay -->  tfluxay_saw
  !      xmuay -->    xmuay_saw
  !      thay -->     thay_saw
  !      pphiay -->   pphiay_saw
  !      pmechay -->  pmechay_saw
  !      vay -->      vay_saw
  !      xiay -->     xiay_saw
  !      xistay -->   xistay_saw
  !      xksidy -->   xksidy_saw
  !
  !   (some mc list arrays are unaffected and so do not have separate
  !   sawtooth copies).
  !
  !   and
  !
  !  in addition the following summary profiles are evaluated
  !
  !   bdenssw(:,isb) -- postsawtooth #/zone, each fast specie
  !   udensplw(:,isb)-- postsawtooth upll, jles/zone, each fast specie
  !   udensppw(:,isb)-- postsawtooth uperp, jles/zone, each fast specie
  !   ebapl2sw(:,isb)-- postsawtooth <epll>, ev, each specie
  !   ebapp2sw(:,isb)-- postsawtooth <epll>, ev, each specie
  !   omegabsw(:,isb)-- postsawtooth angular velocity, 1/sec, each specie
  !
  !---------------------------------------------------------------
  !
  implicit none
  !
  integer, intent(in) :: ind ! ptcl index, current ptcl
  !
  integer :: iwarn,idum
  integer :: itrap,iuncon,ixiav
  integer :: ibr,inz0,inz,igc2,jold,jnew
  !
  logical :: imove_sawtooth ! set t if sawtooth moves this ion
  !
  real(fp) :: zxhw
  real(fp) :: zw,zxi,zth,zv,zr,zksid,zvpll,zvperp
  real(fp) :: zvp,zxisold,zxisnew,zrmaj,zymaj
  real(fp) :: zfac,zbiold,zbinew,zmu,zrand,zfrac,zfracf
  real(fp) :: zxksav,zeiav,zphiav
  real(fp) :: xjadep,fmodbx
  real(fp),parameter :: zpinv=rad2deg/180.0_fp
  !
  !---------------------------------------------------------------
  !
  imove_sawtooth = .false.
  if(flag_cursaw(lsbeam)) then
     !
     if(nlfprod(lsbeam)) then
        zfracf=min(1.0_fp,xswfrac_allfast*xswfrac_fusn)
     else
        zfracf=min(1.0_fp,xswfrac_allfast*xswfrac_beam)
     end if
     zfracf=max(0.1_fp,zfracf) ! at least 10% presumed to participate
     !
     !  this factor converts m.c. ion weight * v**2 to joules
     !
     zfac=xe*v2toev
     !
     !  get current rel. flux location of ion
     zxi=xiay(ind,lsbeam)
     jold=lcentr+zxi*nzones
     !
     !  is it inside the sawtooth mixing region?
     if(jold.gt.kmix_kdsaw) goto 100
     !
     !  it is.  get parameters
     !
     !  add in weight (no. of physical ions/ monte carlo ion)
     zw=wghtay(ind,lsbeam)
     wsum_kdspost(lsbeam)=wsum_kdspost(lsbeam)+zw
     !
     !  poloidal coordinate
     zth=thay(ind,lsbeam)
     !
     !  major radius
     zr=rmjionay(ind,lsbeam)
     !  speed
     zv=vay(ind,lsbeam)
     !  pitch angle (vpll/v)
     zksid=xksidy(ind,lsbeam)
     !
     !  compute vperp,vpll of ion
     zvpll=zksid*zv
     zvperp=sqrt(max(zero,(zv**2-zvpll**2)))
     !
     !  although ptcl is within zone that is at least partially mixed, is
     !  it actually inside the mixing region?
     !
     jnew=jold              ! used if ptcl does not move
     if(zxi.gt.xmix_kdsaw) then
        goto 90
     end if
     !
     !  it is.
     !  does it participate, in case of partial mixing?
     !
     if(zfracf.lt.1.0_fp) then
        zrand=nbrandm_pll(idum)
        if(zrand.gt.zfracf) goto 90
     end if
     !
     !  ok: the sawtooth moves this ion...
     !
     !  evaluate 1/b at pre sawtooth location
     !  ** using b field info stored in the **previous** beam code timestep
     !  fmodbx routine returns 1/b ***
     !
     zbiold=fmodbx(zxi,zth)
     !
     !  conserved quantity thru sawtooth
     zmu=zvperp*zvperp*zbiold
     !
     !  compute new "xi" location for fast ion
     !  random choice
     !
     zrand=nbrandm_pll(idum)
     !
     call kdsaw_find(zrand,jold-lcentr+1,jnew,zfrac)
     jnew=jnew+lcentr-1
     !
     !  new "xi" location
     zxi=xi(jnew,1)+dxi(jnew,2)*zfrac
     xiay_saw(ind,lsbeam)=zxi
     !
     !  randomize "theta"
     zth=-xpi+twopi*nbrandm_pll(idum)
     thay_saw(ind,lsbeam)=zth
     !
     !  compute new velocity coordinates - mu and vpll conserved
     zbinew=fmodbx(zxi,zth)
     zvperp=sqrt(zmu/zbinew)
     !
     zvp=zv
     zv=sqrt(zvpll**2+zvperp**2)
     vay_saw(ind,lsbeam)=zv
     !
     zksid=zvpll/zv
     xksidy_saw(ind,lsbeam)=zksid
     !
     !  radial current calculation -- also get orbit half width
     !
     zxisold=xistay(ind,lsbeam)
     call xs_rzget(zero,zth,zxi,100.0_fp,zrmaj,zymaj,iwarn)
     zxisnew=xjadep(zxi,zth,zrmaj,zymaj,zv,zvpll, &
          itrap,iuncon,zxksav,zeiav,zphiav,ixiav,zxhw)
     call radcur(zxisold,zxisnew,zw,tirsaws)
     xistay_saw(ind,lsbeam)=zxisnew
     njavay_saw(ind,lsbeam)=lcentr+zxisnew*nzones
     !
     !  use orbit half width for smoothing
     !
     dxsm_kdsaw(jnew,lsbeam)=dxsm_kdsaw(jnew,lsbeam)+zw*zxhw
     dxwt_kdsaw(jnew,lsbeam)=dxwt_kdsaw(jnew,lsbeam)+zw
     !
     !  re-init conserved quantities
     !
     imove_sawtooth = .true.
     call ptclsave(xiay_saw(ind,lsbeam),thay_saw(ind,lsbeam), &
          vay_saw(ind,lsbeam),xksidy_saw(ind,lsbeam), &
          rmjionay_saw(ind,lsbeam),xzionay_saw(ind,lsbeam), &
          tfluxay_saw(ind,lsbeam), &
          xmuay_saw(ind,lsbeam),pphiay_saw(ind,lsbeam), &
          pmechay_saw(ind,lsbeam))
     !
     zr=rmjionay_saw(ind,lsbeam)
     !
     !  add up the new part of the density profiles
     !
90   continue
     ibr=int(lcentr+nznbmr*(zxi/xminbm))
     inz0=nthzsm(ibr-1)
     inz=nthzsm(ibr)
     if(nlsym2b) then
        !  midplane symmetry assumed; theta ranging from -pi to pi;
        !  grid ranges from 0 to pi
        igc2=min(inz,int(inz0+1+(abs(zth)*zpinv)*(inz-inz0)))
     else
        !  asymmetric theta grid
        igc2=min(inz, &
             int(inz0+1+((zth-thbdy0)*0.5_fp*zpinv)*(inz-inz0)))
     end if
     bdenssw(jnew,lsbeam)=bdenssw(jnew,lsbeam)+zw
     udensppw(jnew,lsbeam)=udensppw(jnew,lsbeam)+ &
          zfac*zw*zvperp*zvperp
     udensplw(jnew,lsbeam)=udensplw(jnew,lsbeam)+ &
          zfac*zw*zvpll*zvpll
     omegabsw(jnew,lsbeam)=omegabsw(jnew,lsbeam)+ &
          zw*btiltc(2,igc2)*zvpll/zr
     !
100  continue
  end if
  !
  if(.not.imove_sawtooth) then
     ! unaffected ions -- copy the data...
     rmjionay_saw(ind,lsbeam)=rmjionay(ind,lsbeam)
     njavay_saw(ind,lsbeam)=njavay(ind,lsbeam)
     rmjionay_saw(ind,lsbeam)=rmjionay(ind,lsbeam)
     xzionay_saw(ind,lsbeam)=xzionay(ind,lsbeam)
     tfluxay_saw(ind,lsbeam)=tfluxay(ind,lsbeam)
     xmuay_saw(ind,lsbeam)=xmuay(ind,lsbeam)
     thay_saw(ind,lsbeam)=thay(ind,lsbeam)
     pphiay_saw(ind,lsbeam)=pphiay(ind,lsbeam)
     pmechay_saw(ind,lsbeam)=pmechay(ind,lsbeam)
     vay_saw(ind,lsbeam)=vay(ind,lsbeam)
     xiay_saw(ind,lsbeam)=xiay(ind,lsbeam)
     xistay_saw(ind,lsbeam)=xistay(ind,lsbeam)
     xksidy_saw(ind,lsbeam)=xksidy(ind,lsbeam)
  end if
  !
  return
end subroutine pre_sawnbi
!-------------------------------------------------------------------
!
subroutine pre_sawnbi_finish
  !
  !-------------------------------------------------------------------
  !  f90 modules...
  !
  use xstraln_calls
  use nbi_com
  use nbi_random_pll
  !
  !-------------------------------------------------------------------c
  implicit none
  !
  !-------------------------------------------------------------------c
  real(fp) :: znorm
  integer :: iuse,j,isb,imix
  !-------------------------------------------------------------------
  !
  !  is there a species affected by sawteeth?
  !
  iuse=0
  do isb=1,nsfast
     if(flag_cursaw(isb)) iuse = iuse + 1
  end do
  !
  !  just copy profile data outside the mixing region; inside, sum the
  !  total ptcl count
  !
  if(iuse.eq.0) then
     imix=0
  else
     imix=kmix_kdsaw
  end if
  !
  do isb=1,nsfast
     do j=ledge,lcentr,-1
        if((.not.flag_cursaw(isb)).or.(j.gt.imix)) then
           omegabsw(j,isb)=omegabs(j,isb)
        else if(bdenssw(j,isb).gt.zero) then
           omegabsw(j,isb)=omegabsw(j,isb)/bdenssw(j,isb)
        else
           if(j.lt.ledge) then
              omegabsw(j,isb)=omegabsw(j+1,isb)
           else
              omegabsw(j,isb)=zero
           end if
        end if
     end do
  end do
  !
  do j=lcentr,ledge    ! outside mixing region, not everywhere !
     do isb=1,nsfast
        if(flag_cursaw(isb)) then
           if(j.gt.imix) then
              !  just copy outside mixing region
              bdenssw(j,isb)=bdenss(j,isb)
              ebapp2sw(j,isb)=ebapp2s(j,isb)
              ebapl2sw(j,isb)=ebapl2s(j,isb)
              udensppw(j,isb)=udenspp(j,isb)
              udensplw(j,isb)=udenspl(j,isb)
              omegabsw(j,isb)=omegabs(j,isb)
           else
              !  in sawtooth mixing region, ptcls will be moved
              !  and new density summed.  (angular velocity
              wsum_kdspre(isb)=wsum_kdspre(isb)+bdenss(j,isb)
           end if
        else
           !  just copy-- species not affected by sawtooth
           bdenssw(j,isb)=bdenss(j,isb)
           ebapp2sw(j,isb)=ebapp2s(j,isb)
           ebapl2sw(j,isb)=ebapl2s(j,isb)
           udensppw(j,isb)=udenspp(j,isb)
           udensplw(j,isb)=udenspl(j,isb)
           omegabsw(j,isb)=omegabs(j,isb)
        end if
     end do
  end do
  !
  do isb=1,nsfast
     if(flag_cursaw(isb)) then
        !
        !  compute orbit width-based smoothing factor
        !  fill in holes -- zones w 0 ptcls
        !
        do j=ledge,lcentr,-1
           if(dxwt_kdsaw(j,isb).gt.zero) then
              dxsm_kdsaw(j,isb)= &
                   dxsm_kdsaw(j,isb)/dxwt_kdsaw(j,isb)
           else
              if(j.eq.ledge) then
                 dxsm_kdsaw(j,isb)=zero
              else
                 dxsm_kdsaw(j,isb)=dxsm_kdsaw(j+1,isb)
              end if
           end if
        end do
        !
        !  smooth the smoothing coeffcient
        !
        call sawsmoo(dxsm_kdsaw(1:mj,isb),-1)

        if(wsum_kdspost(isb).gt.0) then
           znorm=wsum_kdspre(isb)/wsum_kdspost(isb)
           call sawsmoo(bdenssw(1:mj,isb),1)
           call sawsmoo(udensppw(1:mj,isb),1)
           call sawsmoo(udensplw(1:mj,isb),1)
           call sawsmoo(omegabsw(1:mj,isb),0)

           do j=lcentr,ledge
              !
              !  convert back the new average energies
              !
              if(bdenssw(j,isb).gt.zero) then
                 ebapp2sw(j,isb)=udensppw(j,isb)/(xe*bdenssw(j,isb))
                 ebapl2sw(j,isb)=udensplw(j,isb)/(xe*bdenssw(j,isb))
              else
                 ebapp2sw(j,isb)=zero
                 ebapl2sw(j,isb)=zero
              end if
           end do
        else
           !
           !  no ptcls found -- keep old profiles.
           !
           bdenssw(lcentr:imix,isb)=bdenss(lcentr:imix,isb)
           udensppw(lcentr:imix,isb)=udenspp(lcentr:imix,isb)
           udensplw(lcentr:imix,isb)=udenspl(lcentr:imix,isb)
           ebapp2sw(lcentr:imix,isb)=ebapp2s(lcentr:imix,isb)
           ebapl2sw(lcentr:imix,isb)=ebapl2s(lcentr:imix,isb)
           omegabsw(lcentr:imix,isb)=omegabs(lcentr:imix,isb)
        end if
        !
     end if                          !  mc species test
  end do
  !
contains
  !
  !  smoothing...
  !
  subroutine sawsmoo(zprof,inorm)
    use iso_c_binding, only: fp => c_double
    use physconst_mod, only: zbase=>zero,epsinv
    real(fp) :: zprof(mj)        ! profile to renorm & smooth...
    integer :: inorm        ! =1 for final renorm to conserve total #
    !  =0: no renorm; =-1: smooth the delta, passed as zprof

    real(fp) :: zeps(nzones),zdel(nzones),zwkin(nzones),zwksmoo(nzones)
    real(fp) :: zbc = 0.5_fp
    real(fp) :: zdum1,zsum,znorm2
    real(fp) :: zdel2,zdel1

    integer :: idum1,idum2,ibase

    !  renormalize mixing region

    zdel2 = min(0.1_fp,(2.0_fp/nzones))
    zdel1 = 1.0_fp/nzones

    if(inorm.eq.1) then
       zprof(lcentr:imix)=znorm*zprof(lcentr:imix)
       zsum=sum(zprof(lcentr:imix))
       if(zsum.le.zero) then
          write(nonlin,*) ' %pre_sawnbi(sawsmoo): profile = zero'
          zprof = zero
          return
       end if
    end if

    zeps=epsinv
    if(inorm.eq.1) then
       zdel=max(zdel1,dxsm_kdsaw(lcentr:ledge,isb))
    else
       zdel=zdel2
    end if
    if(inorm.eq.1) then
       zwkin=zprof(lcentr:ledge)/dvolzc_r8(lcentr:ledge)
       ibase=1
    else
       zwkin=zprof(lcentr:ledge)
       ibase=0
    end if
    zwksmoo=zero

    call r8filtr6(xizc(lcentr:ledge),zwkin,zwksmoo,nzones, &
         zeps,nzones,zeps,0,zdel,0,zeps(1),ibase,zbase, &
         zbc,zbc, zdum1,idum1,idum2)

    if(inorm.eq.1) then
       zprof(lcentr:ledge) = zwksmoo*dvolzc_r8(lcentr:ledge)
    else
       zprof(lcentr:ledge) = zwksmoo
    end if

    if(inorm.eq.1) then
       znorm2=zsum/sum(zprof(lcentr:imix))
       zprof(lcentr:imix) = znorm2*zprof(lcentr:imix)
    end if

  end subroutine sawsmoo

end subroutine pre_sawnbi_finish
