module nbi_aep_mod
  use iso_c_binding, only: fp=>c_double
  use physconst_mod, only: aee
  use nbi_aep_com
  use nbi_aep_com, only : npt_mode=>ntaep,nk_aep=>nkaep
  use nbi_com, only: minb,lsbeam,xksid,xmuay,pmechay,xion,th,tdepay,rmjion,yion
  use nbi_com, only: tflux,tfluxay,pphiay,xiay,thay,vay,xksidy,rmjionay,xzionay
  use nbi_com, only: bmaxm,v2toev,nptcls,nbndex,nlfprod,abeams,xzbeams,tirabds,phiprg,mj
  use nbi_com, only: plflxg,ledge,lcentr,ngc,zmp
  use nbi_com, only: wghtay,pbaepms,nstep_aep
  use nbi_com, only: tstop_in=>tbm2,tstart_in=>tbm1
  use nbi_com, only : naeptot,taep_fac,aep_ascale,nltest_output
  use nbi_com, only : xminbm,lep1,nzones,nznbmr,nthzsm,nlsym2b,thbdy0
  use nbi_com, only : bimodb,bixi,bithst,bijac,xbmbnd,vion,wght,nsjdotb,nsnccw,nlbflr,nflr
  use nbi_com, only : big,btori,bi1ovb,nfast_map,twopi,aep_uscale,workpath,runid,delt,nsbeam
  use nbi_mpi_com
  use nbi_random_aep
  use logmod
  use ezspline_obj, only : ezspline_allocated
  use ezspline_type, only :  ezspline1,ezspline2 
  use ezspline, only : ezspline_interp,ezspline_error,ezspline_setup,&
       ezspline_free,ezspline_init,ezlinear_init
 
  implicit none
  private :: fp
  integer :: id_ptcl !mc particle ind (from orball)
  real(fp),allocatable,dimension(:,:,:) :: sp_f2 ! function values
  integer :: ier_a !error related to  spl_amode
  integer :: ier_psi !error related to  spl_rhophi
  integer :: ier_p !error related to spline2_o
  integer :: cnt_ax
  real(fp), dimension(1)  :: rfbm,zfbm,polfbm,bfbm,mufbm,pphfbm,efbm,tdepfbm,pfbm
  real(fp), dimension(1)  :: thfbm,newvalfbm    ! initialize flag newvalfbm, 
                                                         ! which at output indicates whether 
                                                         ! particle's variables have been 
                                                         ! modified by this script 
                                                         !(=1) or not (=0)
  real(fp) :: psi_lcfs,psiw
  real(fp) :: gout !g function
  real(fp) :: baxis,xiayp,pmechayp
  real(fp) :: b_min_lcfs,b_max_lcfs,gfun_w,gfun_ax, r_axis
  real(fp) :: fct_pz_xa ! normalization factor used in get_femu
  real(fp) :: fct_pz ! normalization factor used in get_femu, 
                     ! fct_pz=fct_pz_xa*xzbeams/sqrt(abeams)
  real(fp) :: pestep,pvar_e1,pvar_en
  real(fp), parameter :: tol_max=1.0e-8_fp ! [nominal: 1d-8]
  integer, parameter :: npsi=10001
  integer, parameter :: mnkaep=10 !must be consistent with nbspec.dat
  integer, dimension(:,:),allocatable  :: imap_nkaep_use !mapping array
  real(fp),dimension(mnkaep) ::  dtstep_mc 
  real(fp) ::  min_ene,max_ene
  integer,dimension(mnkaep) :: nstep_mc
  real(fp), dimension(:,:), allocatable :: astep_mc
  real(fp), dimension(:,:), allocatable :: tstep_mc
  integer, dimension(:,:), allocatable :: valid
  real(fp) :: tstart,zmp_zelonzc,tstop
  integer,dimension(mnkaep) :: nstep_mc_in,nstep_mc_out
  integer :: ntry_evlv_max    ! number of attempt to extract de,dp
  logical, save :: il_old=.false.
  real(fp) :: u_ax,u_lcfs
  real(fp) :: epot0 !potential energy before the kick, ev
  real(fp),dimension(:,:), allocatable :: aepd2 ! aepd(1:npde,1:npdp) interpolated in (e, pphi, mu)
  integer :: npde_max,npdp_max
  integer :: indmin(1),indmax(1)
  type(ezspline1) :: spl_rhophi ! 1d interpolation
  type(ezspline1), allocatable, dimension(:) :: spl_amode ! 1d interpolation
  type(ezspline2), allocatable, dimension(:) :: spline2_o ! 2d interpolation


contains

  subroutine aep_driver(bi_time)
    use iso_c_binding, only: fp => c_double
    real(fp) :: bi_time ! time step duirng orbiting
    integer :: i,ierr,k,k_use
    integer,save :: id_ptclp=0 ! id ptcl from prev time step
    real(fp),save :: tdep_ay=0.0_fp
    ierr=0
    call aep_init_ptcl(bi_time,ierr)
    if(ierr.ne.0) then
       return !bad particle
    endif
    if(id_ptclp.ne.id_ptcl) then !new ptcl
       !get first kick
       id_ptclp=id_ptcl
       tdep_ay=1.e3_fp*max(tstart_in,tdepay(id_ptcl,lsbeam))
       do k=1,nkaep_use
          k_use=imap_nkaep_use(k,lsbeam)
          if(k_use.eq.0) cycle
          if(nstep_mc(k_use).eq.0)cycle
          do i=1,nstep_mc(k_use)
             if(tdep_ay.lt.tstep_mc(i,k)) then
                nstep_mc_in(k_use)=i
                exit
             endif
          enddo
          do i=nstep_mc_in(k_use),nstep_mc(k_use)
             if(tdepfbm(1).lt.tstep_mc(i,k)) then
                nstep_mc_out(k_use)=i
                exit
             endif
          enddo
       enddo
    else
       do k=1,nkaep_use
          k_use=imap_nkaep_use(k,lsbeam)
          if(k_use.eq.0) cycle
          if(nstep_mc(k_use).eq.0)cycle
          nstep_mc_in(k_use)=nstep_mc_out(k_use)+1
          do i=nstep_mc_in(k_use)-1,nstep_mc(k_use)
             if(tdepfbm(1).lt.tstep_mc(i,k)) then
                nstep_mc_out(k_use)=i
                exit
             endif
          enddo
        enddo
    endif
    tdep_ay=tdepfbm(1)
   
    call evolve_fbm(nstep_mc_in,nstep_mc_out)
    il_old=.false.
    
  end subroutine aep_driver

  subroutine aep_init(ierr)
    use iso_c_binding, only: fp => c_double
    integer :: k,ibs
    real(fp) :: temp(1)
    integer :: ierr
    real(fp) :: max_mamp
    integer, dimension(:,:), allocatable :: noccur
    integer :: nstep_mc0
    real(fp) :: max_amp, min_estep, min_pzstep, max_destep, max_dpzstep, min_ntsteps
    real(fp) :: ppstep,pvar_de,pvar_dp,scale_face,scale_facpz,tdum1
    integer :: kt,j1,k_use
    real(fp), parameter :: min_ene_par=1.e12_fp    ! min energy
    real(fp) :: chi_b_min,chi_b_max
    integer :: isize, imatch,in_aep,isb
    logical :: il_match

!   call infolog('entering aep_init, mode # =',nk_aep)
    ierr=1
    np_e=0
    np_p=0
    np_m=0
    naeptot=0
    if(minb.eq.0) return
    if(nk_aep.le.0) return
    ierr=0
    if(allocated(imap_nkaep_use)) deallocate(imap_nkaep_use)
    allocate(imap_nkaep_use(nk_aep,nsbeam))
    imap_nkaep_use=0
    nkaep_use=0
    indmin=count(tstart_in.ge.aept)
    indmin=max(indmin(1),1)
    indmax=count(tstop_in.gt.aept)+1
    indmax=min(indmax(1),size(aept))
    if(tstart_in.ge.aept(indmax(1))) then
       ierr=2
       return
    endif
    if(tstop_in.le.aept(indmin(1))) then
       ierr=2
       return
    endif
    do k=1,nk_aep !# of mode
       if(aep_ascale(k).gt.0.0_fp)then
          if(maxval(aepa(indmin(1):indmax(1),k)).gt.0.0_fp) then
             nkaep_use=nkaep_use+1
             if(int(aepsp(k)).eq.0) then
                imap_nkaep_use(nkaep_use,1:nsbeam)=k
                cycle
             endif
             do isb=1,nsbeam
                call match_fi(k,isb,il_match)
                if(il_match) then
                   imap_nkaep_use(nkaep_use,isb)=k
                endif
             enddo
          endif
       endif
    enddo
    if (maxval(imap_nkaep_use).eq.0) then
       ierr=1
       return !no aep activities
    endif
    if(nltest_output) then
       if(allocated(pbaep_ms)) deallocate(pbaep_ms)
       if(allocated(pbaep_pol)) deallocate(pbaep_pol)
       if(allocated(pph_ms)) deallocate(pph_ms)
       if(allocated(aep_count)) deallocate(aep_count)
       np_e=maxval(aepe(:,1))
       np_p=maxval(aepp(:,1))
       np_m=maxval(aepmu(:,1))
       allocate(pbaep_ms(np_e,np_p,np_m,nkaep_use,1))
       allocate(pbaep_pol(n_norm_pol,nkaep_use,1))
       allocate(pph_ms(np_e,np_p,np_m,nkaep_use,1))
       allocate(aep_count(np_e,np_p,np_m,nkaep_use,1))
       pbaep_ms=0.0_fp
       pbaep_pol=0.0_fp
       pph_ms=0.0_fp
       aep_count=0.0_fp
    endif
    call create_spla(ierr)
    if(ierr.ne.0) then
       call errorlog(' aep_init failed in create_spla')
       call bad_exit
    endif
    call create_splp(ierr)
    if(ierr.ne.0) then
       call errorlog(' aep_init failed in create_splp')
       call bad_exit
    endif

    baxis = bmaxm(lcentr)       !  tesla
    ! get min/max b at lcfs, here chi angle location of min/max(b)
    call eq_glimb1(1.0_fp,0.0_fp,b_min_lcfs,chi_b_min,b_max_lcfs,chi_b_max,ierr)
    if(ierr.ne.0) then
       call errorlog(' aep_init failed: cannot get b_min, b_max at lcfs')
       return
    endif
    call get_g(1.0_fp,gfun_w)
    call get_g(0.0_fp,gfun_ax)


    ! define normalization factors, fct_pz, used in get_femu
    zmp_zelonzc=zmp/aee
!mg    fnorm=(xzbeams(lsbeam)*baxis)**2
! mario's (1d19*zmp/(zel/zc)) = 1.0e-5_fp*1.6726_fp/4.8032_fp*2.9979_fp
!mg    fct_pz=sqrt(fnorm/(2._fp*abeams(lsbeam)*zmp_zelonzc))
!mg     same as above
    fct_pz_xa=baxis*sqrt(0.5_fp/zmp_zelonzc)
 
    ! get max flux (at lcfs)
    temp(1)=plflxg(lcentr)
    psi_lcfs=plflxg(ledge)

    if (temp(1).gt.psi_lcfs) psi_lcfs = temp(1)
    if(psi_lcfs.eq.0.0_fp) then
       call errorlog(' aep_init psi_lcfs=0')
       call bad_exit
    endif
    call create_splpsi(ierr)
    if(ierr.ne.0) then
       call errorlog(' aep_init failed in create_splpsi')
       call bad_exit
    endif
    psiw = psi_lcfs/abs(baxis)

    tstart=max(aept(indmin(1)),tstart_in)
    tstop=min(aept(indmax(1)),tstop_in)
    nstep_mc(1:nk_aep)=0
    max_amp=0.0_fp
    min_estep=1.e12_fp
    min_pzstep=1.e12_fp
    max_destep=0._fp
    max_dpzstep=0._fp
    min_ntsteps=5.e2_fp
    do k=1,nkaep_use
       k_use=maxval(imap_nkaep_use(k,1:nsbeam))
       if(abs(taep_fac(k_use)).gt.0.0_fp) then 
          ! re-define time step [s] for inner montecarlo loop
          ! based on the ntaep_fac provided as input.
          ! in this case, no check on the minimum number
          ! of mc steps is performed.
          nstep_mc(k_use)=int((tstop-tstart)*abs(taep_fac(k_use))/&
               aeptstep(k_use)) !tstop,tstart [s]
          cycle
       endif
       !if(taep_fac.eq.0.0_fp)
       ! assign default value 
       ! compute max mode amplitude for each mode at this
       ! nubeam step. also compute max (de,dpz) grid steps and
       ! minimum (e,pz) grid steps.
       nstep_mc0=int((tstop-tstart)/aeptstep(k_use)) !tstop,tstart [s]
       do j1=indmin(1),indmax(1)
          if(aepa(j1,k_use).ge.max_amp) then
             max_amp=aepa(j1,k_use)
          endif
       enddo
       pestep=(aepe(k_use,3)-aepe(k_use,2))*1.e-3_fp/(aepe(k_use,1)-1) !step in energy, [kev]
       ppstep=(aepp(k_use,3)-aepp(k_use,2))/(aepp(k_use,1)-1)
       pvar_de=(aepde(k_use,3)-aepde(k_use,2))*1.e-3_fp/(aepde(k_use,1)-1)
       pvar_dp=(aepdp(k_use,3)-aepdp(k_use,2))/(aepdp(k_use,1)-1)
       if(pestep.le.min_estep) min_estep=pestep
       if(ppstep.le.min_pzstep) min_pzstep=ppstep
       if(pvar_de.ge.max_destep) max_destep=pvar_de
       if(pvar_dp.ge.max_dpzstep) max_dpzstep=pvar_dp
       ! define new multiplier based on max mode
       ! amplitude and (e,de) grid steps.
       scale_face=max_destep/min_estep
       scale_facpz=max_dpzstep/min_pzstep
       ! re-define time step [ms] for inner montecarlo loop
       ! based on grid steps, max mode amplitude.
       nstep_mc(k_use)=max_amp*min_ntsteps*max(scale_face,scale_facpz)
       nstep_mc(k_use)=max(nstep_mc0,nstep_mc(k_use)) 
    enddo
    if(maxval(nstep_mc).eq.0) then
       ierr=1 !aep not started yet
       return
    endif
    dtstep_mc=0.0_fp
    in_aep=0
    do k=1,nk_aep
       if(nstep_mc(k).le.0) cycle
       dtstep_mc(k)=(tstop-tstart)/nstep_mc(k) !tstop,tstart [s]
       in_aep=in_aep+1
    enddo
    if(in_aep.gt.0) then
       ntry_evlv_max=sum(nstep_mc(1:nk_aep))/in_aep*10
    else
       ierr=1
       return
    endif
!c --------------------------------------------------------
!c     compute average mode amplitude for each mode at
!c     each mc step
      
!c     initialize arrays 
    if(allocated(valid)) deallocate(valid)
    allocate(valid(maxval(nstep_mc),nkaep_use))
    valid=0
    if(allocated(astep_mc)) deallocate(astep_mc)
    allocate(astep_mc(maxval(nstep_mc),nkaep_use))
    astep_mc=0.0_fp
    if(allocated(noccur)) deallocate(noccur)
    allocate(noccur(maxval(nstep_mc),nkaep_use))
    noccur=0
    if(allocated(tstep_mc)) deallocate(tstep_mc)
    allocate(tstep_mc(maxval(nstep_mc),nkaep_use))
    do k=1,nkaep_use
       k_use=maxval(imap_nkaep_use(k,1:nsbeam))
       if(nstep_mc(k_use).eq.0) cycle
       do kt=1,nstep_mc(k_use)
          tstep_mc(kt,k)=1.e3*(tstart+(kt-1)*dtstep_mc(k_use)+&
               dtstep_mc(k_use)/2.0_fp) ! [ms]
       enddo
    enddo
!c     compute average amplitude
    ! if modes are specified on a time grid which has
    ! time step smaller than that of the inner mc loop,
    ! compute amplitude as the average of the input
    ! values at each mc sub-step
    do k=1,nkaep_use
       k_use=maxval(imap_nkaep_use(k,1:nsbeam))
       if((aept(indmin(1)+1)-aept(indmin(1))).lt.dtstep_mc(k_use)) then ! use method #1
          do j1=indmin(1),indmax(1)
             kt=int((aept(j1)-tstart)/dtstep_mc(k_use)) ! aept(j1) [s]
             if(kt.ge.1.and.kt.le.nstep_mc(k_use)) then
                astep_mc(kt,k)=astep_mc(kt,k)+aepa(j1,k_use)*aep_ascale(k_use)
                noccur(kt,k)=noccur(kt,k)+1
                if(astep_mc(kt,k).gt.0.0_fp)then
                   valid(kt,k)=1 ! added 03/24/2014, mp
                endif
             endif
          enddo
             
       ! if modes are specified on a time grid which has
       ! time step larger than that of the inner mc loop,
       ! compute amplitude by interpolating the input
       ! values at each mc sub-step
       
       else
          do j1=1,nstep_mc(k_use)
!             indmax=minloc(aept,aept.ge.tstep_mc(j1,k)/1.e3_fp) !right bndy
!             indmin=maxloc(aept,aept.le.tstep_mc(j1,k)/1.e3_fp) !left bndy
                                                             ! could be same
             if(maxval(aepa(indmin(1):indmax(1),k)).gt.0.0_fp) then
             ! if original amode gt 0 then get astep_mc from spline
                call ezspline_interp(spl_amode(k),tstep_mc(j1,k)/1.e3_fp,astep_mc(j1,k),ier_a)
                call ezspline_error(ier_a)	
                if (ier_a.ne.0) then 
                   call errorlog(' evolve_mc failed interpolate a')
                   call bad_exit
                endif
             else
                ! original amode eq 0
                astep_mc(j1,k)=0.0_fp
             endif
             astep_mc(j1,k)=astep_mc(j1,k)*aep_ascale(k_use)
             noccur(j1,k)=noccur(j1,k)+1
             if(astep_mc(j1,k).gt.0.0_fp)then
                valid(j1,k)=1 
             endif
          enddo
       endif
    enddo
    !    write(0,'(a,i3,a)') '	using ',nk_aep,' mode(s):'
    do k=1,nkaep_use
       k_use=maxval(imap_nkaep_use(k,1:nsbeam))
       if(nstep_mc(k_use).eq.0) cycle
       if(maxval(valid(1:nstep_mc(k_use),k)).eq.0) cycle
       max_mamp=0.
       do j1=1,nstep_mc(k_use)
          if(noccur(j1,k).gt.1)then
             astep_mc(j1,k)=astep_mc(j1,k)/noccur(j1,k)
          endif
          if(astep_mc(j1,k).ge.max_mamp) then 
             max_mamp=astep_mc(j1,k)
          endif
       enddo
    enddo
!c --------------------------------------------------------
!c     find minimum & maximum energy over which
!c     p(de,dp) is 'active' in this time window
    min_ene=min_ene_par
    max_ene=0._fp
    do k=1,nkaep_use
       k_use=maxval(imap_nkaep_use(k,1:nsbeam))
       if(nstep_mc(k_use).eq.0) cycle
       pestep=(aepe(k_use,3)-aepe(k_use,2))*1.e-3_fp/(aepe(k_use,1)-1) !step in energy, [kev]
       pvar_e1=aepe(k_use,2)*1.e-3_fp
       pvar_en=aepe(k_use,3)*1.e-3_fp
       do j1=1,nstep_mc(k_use)
          if(noccur(j1,k).gt.0.and.pvar_e1-pestep/2._fp.le.min_ene) then
             min_ene=pvar_e1-pestep/2._fp
          endif
          if(noccur(j1,k).gt.0.and.pvar_en+pestep/2._fp.ge.max_ene) then
             max_ene=pvar_en+pestep/2._fp
          endif
       enddo
    enddo
    if(allocated(noccur)) deallocate(noccur)
    if(allocated(aepd2)) deallocate(aepd2)
    npde_max=maxval(aepde(:,1))
    npdp_max=maxval(aepdp(:,1))
    allocate(aepd2(npde_max,npdp_max)) !do it once for all modes
    
    
  end subroutine aep_init
  
  subroutine aep_init_ptcl(bi_time,ierr)
    use iso_c_binding, only: fp => c_double
    use nbi_com, only : id_psi
    integer :: ivalid,k,i
    real(fp) :: temp(1),rhodum,bi_time
    real(fp) :: vdum,dumj(2,2),knorm,eprime
    real(fp) :: fnbint
    integer :: ierr,idummy,ngc_wrk
!    call infolog('entering aep_init_ptcl')
    ierr=0

    ivalid=0 ! counter for 'valid' particles in particle list (fbm)

    ! initialize flag newvalfbm, which at output indicates whether particle's
    ! variables have been modified by this script (=1) or not (=0)
    newvalfbm=0
    ! initialize deposition time for current particle
    tdepfbm=0._fp
    
    ! estimate # of 'valid' mc particles 
    cnt_ax=0

    k=id_ptcl
    if(k.eq.0.or.wght.le.0.0_fp.or.(tdepay(k,lsbeam).eq.0.0_fp).or.&
         &(vion.le.0.0_fp).or.(xion.ge.1.0_fp).or.&
         &((ngc.lt.lcentr).or.(ngc.gt.ledge))) then
       ierr=-1
       return !bad particle
    endif

    fct_pz=fct_pz_xa*xzbeams(lsbeam)/sqrt(abeams(lsbeam))

    cnt_ax=1

    u_ax=xzbeams(lsbeam)*fnbint(0.0_fp,phiprg,1,mj,1,idummy,ngc_wrk)*1.e-3_fp
    u_lcfs=xzbeams(lsbeam)*fnbint(1.0_fp,phiprg,1,mj,1,idummy,ngc_wrk)*1.e-3_fp
    u_ax=aep_uscale*u_ax
    u_lcfs=aep_uscale*u_lcfs
    do i=1,cnt_ax
       thfbm(i)=th
       rfbm(i)=rmjion*1.e-2_fp	! r [m]
       zfbm(i)=yion*1.e-2_fp	! z [m]
       bfbm(i)=bimodb(1) !tesla
       tdepfbm(i)=1.e3_fp*(tstart_in+bi_time) ! [ms]
       vdum=vion ! [cm/s]
       pfbm(i)=xksid	! pitch in b see errata in roscoe white book
       efbm(i)=v2toev*vdum*vdum*1.e-3_fp ! energy [kev]
!               v2toev*1.e-3_fp = .5e-3_fp*m2cf*mass2charge/1.e4_fp 
!               v2toev = .5 * m2cf*mass2charge/1.e4_fp              
!               v2toev*1.e4_fp = .5 * m2cf*mass2charge              
       call eq_rgetf(1,xion,id_psi,0,temp(1),ierr)
       if(ierr.ne.0)then
          return
       endif
       polfbm(i)=temp(1)/psi_lcfs ! pol. flux coord. normalized to psi @ lcfs
! compute magnetic moment, mu=e/b*(1-p^2) 
       mufbm(i)=efbm(i)/bfbm(i)*(1.0_fp-pfbm(i)*pfbm(i)) ! un-normalized definition
!mufbm(k)=baxis/bfbm(k)*(1.d0-pfbm(k)*pfbm(k)) ! mu'=mu bo/e
! note:  the value used for p(de,dp) mu-axis is the
!	normalized mu, which is:
!  mu_prime=baxis/bfbm(k)*(1.d0-pfbm(k)*pfbm(k))=mufbm(k)*bo/efbm(k)

!	       knorm=a*(1d19*zmp/(zel/zc))/(zcharge*baxis)**2
! (1d19*zmp/(zel/zc)) = 1.0e-5_fp*1.6726_fp/4.8032_fp*2.9979_fp
!              use orbit definition for p_zeta
       knorm=abeams(lsbeam)*zmp_zelonzc/(xzbeams(lsbeam)*baxis)**2
       eprime=knorm*efbm(i) ! energy in normalized units
       rhodum=pfbm(i)*sqrt(2._fp*eprime)*baxis/bfbm(i)
       gout=1.e-6_fp*big(1) ! now  b in tesla r in m
! p_zeta
       pphfbm(i)=-polfbm(i)+(rhodum*gout/baxis)/psiw	! toroidal angular momentum
	       ! note:  in the above definition of p_phi, the value psi_w of psi at the lcfs
	       !	is normalized to the value of b_0 to be consistent with orbit.
       epot0=aep_uscale*xzbeams(lsbeam)*fnbint(xion,phiprg,1,mj,1,idummy,ngc) ! potential energy [ev]
       efbm(i)=efbm(i)+epot0*1.e-3_fp ! energy [kev]
       if(efbm(i).le.0.0_fp) then
          call errorlog(" nbi_aep_mod: efbm(i) less than 0, check aep_uscale")
          call bad_exit(" nbi_aep_mod: efbm(i) less than 0")
       endif

    enddo
!    call infolog('exiting aep_init_ptcl')
    
  end subroutine aep_init_ptcl

  subroutine match_fi(k,isb,ilmatch)
    logical  :: ilmatch
    integer  :: k,isb
    real(fp) :: xz,xa
    
    ilmatch=.false.

    if(int(aepsp(k)).eq.1) then !"D"
       xz=1.0_fp
       xa=2.0_fp
    elseif(int(aepsp(k)).eq.2) then !"T"
       xz=1.0_fp
       xa=3.0_fp
    elseif(int(aepsp(k)).eq.3) then !"HE3"
       xz=2.0_fp
       xa=3.0_fp
    elseif(int(aepsp(k)).eq.4) then !"HE4"
       xz=2.0_fp
       xa=4.0_fp
    elseif(int(aepsp(k)).eq.5) then !"p"
       xz=1.0_fp
       xa=1.0_fp
    endif

    if(max(abs(xz-xzbeams(isb)),abs(xa-abeams(isb))/&
         max(1.0_fp,max(xa,abeams(isb)))) &
         .le.0.05_fp) then
       ilmatch=.true.
    endif

    return
    
  end subroutine match_fi
  
  

  subroutine evolve_fbm(nstep_mc_in,nstep_mc_out)
!ccc- this routine returns a pair of random kicks
!ccc- in (e,p) based on the p(de,dp|e,p,mu) probability distribution
!ccc- at given (e,p,mu) coordinates (or "bins")
    use iso_c_binding, only: fp => c_double
    integer :: imaxsamp    ! set max number of micro-step after which
                                          ! new de,dp are sampled
    real(fp), parameter :: lowecutoff=5.           ! [kev]  
    integer :: samp_err   !sampling error
    integer :: kf,k,kt,kp,k_use,i
    real(fp) ::  pdelt_e,pdelt_p
    real(fp) ::  b_orig
    integer :: npe,npp,npm
    real(fp) ::  pvar_p1,pvar_pn,pvar_mu1,pvar_mun
    real(fp) ::  pp_com
    real(fp) ::  dtsc
    integer :: ind(1),ier
    real(fp) :: rnd1(1)
    real(fp) :: dum1,dum2,dum3,dum4,dum5
    integer :: k_ptcl
    real(fp) :: dumj(2,2),zepp,zep
!mgr    integer, allocatable :: seeds(:) 
    integer :: idum
    integer,dimension(mnkaep) :: nstep_mc_in,nstep_mc_out
    real(fp) ::  e_orig,pph_orig,th_orig,pol_orig,mu_orig
    integer :: do_update
    integer,save :: isamp,is_ok
    ! define 2 nearest integer bins in (e,p,mu) phase space
    integer,dimension(2,mnkaep) :: ije,ijph,ijmu,ijpol
    integer,dimension(2,mnkaep),save :: ijeold,ijphold,ijmuold
    integer,dimension(mnkaep),save :: do_samp
    real(fp),dimension(3,mnkaep) :: dist !1- dist_e, 2 dist_ph,3 -dist_mu
    real(fp),dimension(2) :: dist_e,dist_ph,dist_mu, dist_pol!1- dist_e, 2 dist_ph,3 -dist_mu
    integer,save :: ntry_evlv 
    integer,save :: id_ptcl_p=0 ! id ptcl from prev time step
    integer :: ntry_evlv_p
    real(fp),dimension(mnkaep),save :: delt_e,delt_p,p_com,sgn,e_fbm
    real(fp) :: zbdir(4)
    real(fp) :: btgcp, btgci !b_tor before and after aep kick
    logical,save :: il_first
    real(fp) :: xion1,th1
    real(fp) :: ztflux,zmu,zpphi,zmo,ztfluxp,zmup,zpphip,zmop
    real(fp) :: psi_guess
    real(fp) :: thp,rmjionp,yionp
    integer:: ngcp,nflrp,idummy,ie,ip,imu,ipol
    real(fp) :: fnbint

    real(fp) :: f_damp ! "damp" kicks for particles near a turning point
    integer :: nstep_mc_in_min, nstep_mc_out_max
!    call infolog('entering evolve_fbm, mode # =',nkaep_use)
!mgr    call random_seed(size=i_seeds) !this is added to fix gcc
!mgr    allocate(seeds(i_seeds))
!c ------------------------------------------------------------
!c     start main mc loop over particles
    
    imaxsamp=1e4 
      ! define/init counter of total number of particles processed
    do kp=1,1       ! store original values
       e_orig=efbm(kp)
       e_fbm=e_orig
       pph_orig=pphfbm(kp)
       th_orig=thfbm(kp)
       pol_orig=polfbm(kp)
       mu_orig=mufbm(kp)
       b_orig=bfbm(kp)
       ! skip this particle if its energy is outside the
       ! range over which p(de,dp|e,p,mu) is defined
       if(efbm(kp).gt.0.0_fp.and. &
            &   efbm(kp).ge.min_ene.and. &
            &   efbm(kp).le.max_ene.and. &
            &   rfbm(kp).gt.0.0_fp.and. &
            &   polfbm(kp).le.1.0_fp) then
          ntry_evlv_p=0
          if(id_ptcl_p.ne.id_ptcl) then
             call infolog(' evolve_fbm: ntry_evlv= ',ntry_evlv)
             ntry_evlv=1 ! new ptcl
             il_first=.true.
          else
             if(ntry_evlv.eq.ntry_evlv_max) then
                return !sorry, this ptcl used all reserve
             endif
           endif

111       continue ! entry_evlv point for mc iterations   
    
          aepd2=0.0_fp
          ! start loop over time steps for each particle
          if((id_ptcl_p.ne.id_ptcl).or.ntry_evlv.eq.ntry_evlv_p) then
             do_samp=1
             isamp=0
             is_ok=0
             ijeold=0
             ijphold=0
             ijmuold=0
             delt_e=0.0_fp
             delt_p=0.0_fp
             p_com =0.0_fp
             sgn=0.0_fp
             id_ptcl_p=id_ptcl
          endif
          nstep_mc_in_min=minval(nstep_mc_in,mask=nstep_mc_in.gt.0)
          nstep_mc_out_max=maxval(nstep_mc_out)
          if(nstep_mc_out_max.lt.nstep_mc_in_min) then
             return
          endif
          do  kt=nstep_mc_in_min,nstep_mc_out_max
             do_update=0 ! will be set =1 if any parameter changes
             do kf=1,nkaep_use ! loop over possible p(de,dp) functions
                k_use=imap_nkaep_use(kf,lsbeam)
                if(k_use.eq.0) cycle
                if(lno_tstep(kt,nstep_mc_in(k_use),nstep_mc_out(k_use))) cycle
                npe=aepe(k_use,1)
                npp=aepp(k_use,1)
                npm=aepmu(k_use,1)
                pvar_e1=aepe(k_use,2)*1.e-3_fp
                pvar_en=aepe(k_use,3)*1.e-3_fp
                pvar_p1=aepp(k_use,2)
                pvar_pn=aepp(k_use,3)
                pvar_mu1=aepmu(k_use,2)
                pvar_mun=aepmu(k_use,3)
                
                !c         find particle's "bin" in (e,p,mu) - these are float type
                call get_bin(efbm(kp),npe,pvar_e1,pvar_en,ije(1:2,k_use),dist(1,k_use))
                call get_bin(pphfbm(kp),npp,pvar_p1,pvar_pn,ijph(1:2,k_use),dist(2,k_use))
                call get_bin(mufbm(kp)/efbm(kp)*baxis,npm,pvar_mu1,pvar_mun,&
                     ijmu(1:2,k_use),dist(3,k_use))
             enddo
                
             do kf=1,nkaep_use ! loop over possible p(de,dp) functions
                k_use=imap_nkaep_use(kf,lsbeam)
                if(k_use.eq.0) cycle

                if(lno_tstep(kt,nstep_mc_in(k_use),nstep_mc_out(k_use))) cycle
                if(taep_fac(k_use).lt.0) then
                   imaxsamp=1 !independent kicks
                endif
                ! skip this time step if mode amplitude is zero or
                ! if particle has not been 'deposited' yet
                if(valid(kt,kf).gt.0.and. &
                     &       efbm(kp).gt.0.) then    !check if it lost
                   npe=aepe(k_use,1)
                   npp=aepp(k_use,1)
                   npm=aepmu(k_use,1)
                   pvar_e1=aepe(k_use,2)*1.e-3_fp
                   pvar_en=aepe(k_use,3)*1.e-3_fp
                   pvar_p1=aepp(k_use,2)
                   pvar_pn=aepp(k_use,3)
                   pvar_mu1=aepmu(k_use,2)
                   pvar_mun=aepmu(k_use,3)
                   if ((ije(1,k_use).ge.1.and.ije(1,k_use).le.npe).and. &
                        & (ijph(1,k_use).ge.1.and.ijph(1,k_use).le.npp).and. &
                        & (ijmu(1,k_use).ge.1.and.ijmu(1,k_use).le.npm)) then ! ok, valid bin - proceed
                      do_update=1
!c	    get new de,dp kicks if particle moves to a new/different bin in (e,p,mu)	  
                      if (kt.eq.1.or.ije(1,k_use)-ijeold(1,k_use).ne.0.or. &
                           & ije(2,k_use)-ijeold(2,k_use).ne.0.or.   &
                           & ijph(1,k_use)-ijphold(1,k_use).ne.0.or. &
                           & ijph(2,k_use)-ijphold(2,k_use).ne.0.or. &
                           & ijmu(1,k_use)-ijmuold(1,k_use).ne.0.or. &
                           & ijmu(2,k_use)-ijmuold(2,k_use).ne.0.or. &
                           & do_samp(k_use).eq.1.or.isamp.ge.imaxsamp) then
                         
                         !update de, dp and "skew" parameter for current particle 	      
112                      continue
!mg                         call sample_dedp(ije(k_use),ijph(k_use),ijmu(k_use),kf,pdelt_e,pdelt_p,pp_com,samp_err)
                         call sample_dedp_cell(ije(1:2,k_use),ijph(1:2,k_use),ijmu(1:2,k_use), &
                              & dist(1:3,k_use), kf,pdelt_e,pdelt_p,pp_com,samp_err)
                         if (samp_err.gt.0.and.ntry_evlv.le.ntry_evlv_max) then
                            ntry_evlv=ntry_evlv+1
                            ntry_evlv_p=ntry_evlv
                            goto 112 ! error during sampling process, retry
                         endif
                         if (samp_err.gt.0.and.ntry_evlv.gt.ntry_evlv_max) then
                            e_fbm=e_orig
                            efbm(kp)=e_orig
                            pphfbm(kp)=pph_orig
                            thfbm(kp)=th_orig
                            polfbm(kp)=pol_orig
                            mufbm(kp)=mu_orig
                            bfbm(kp)=b_orig
                            newvalfbm(kp)=0
                            goto 764 ! error during sampling process, exit
                         endif
                         ! reset flag for initial mc sampling
                         do_samp(k_use)=0
                         
                         ! reset counter
                         if(isamp.ge.imaxsamp) then
                            isamp=0 
                            do_samp(k_use)=1
                         endif
                         
                         delt_e(k_use)=pdelt_e
                         delt_p(k_use)=pdelt_p
                         p_com(k_use)=pp_com
                      endif
                        
                      isamp=isamp+1
! extract random sign for the kicks -
! this provides the 'biased sampling' which
! mimics the 'correlated random walk' experienced 
! by (resonant) particles.
                      rnd1(1)=nbrandm_aep(idum)
                      sgn(k_use)=dsign(1.0_fp,rnd1(1)-p_com(k_use))
                      
                   else ! invalid bin - added 02/17/2014, mp
                      
                      delt_e(k_use)=0.0_fp
                      delt_p(k_use)=0.0_fp
                      p_com(k_use)=0.5_fp
                      sgn(k_use)=0.0_fp
                   endif !  update kicks, finite amplitude & inside domain

                   !   now (e,pz,mu) bin indexes are arrays
                   ijeold(1:2,k_use)=ije(1:2,k_use)
                   ijphold(1:2,k_use)=ijph(1:2,k_use)
                   ijmuold(1:2,k_use)=ijmu(1:2,k_use)
                  
                else ! out-of-domain (e,p,mu)
                   
                   delt_e(k_use)=0.0_fp
                   delt_p(k_use)=0.0_fp
                   p_com(k_use)=0.5_fp
                   sgn(k_use)=0.0_fp
                endif ! valid domain (e,p,mu)
                
             enddo ! loop over possible p(de,dp)
!c         update particle's coordinates	(e,pz)
             if (do_update.eq.1) then 
                
                ! update energy, pz
                do kf=1,nkaep_use
                   k_use=imap_nkaep_use(kf,lsbeam)
                   if(k_use.eq.0) cycle

                   if(lno_tstep(kt,nstep_mc_in(k_use),nstep_mc_out(k_use))) &
                        cycle
                   if(valid(kt,kf).gt.0) then
                      dtsc=dtstep_mc(k_use)*sqrt(abs(taep_fac(k_use)))/aeptstep(k_use) !kick is ~sqrt(abs(teap_fac))
                      e_fbm(k_use) = e_fbm(k_use) + sgn(k_use)*delt_e(k_use)*astep_mc(kt,kf)*dtsc  
                      efbm(kp) = efbm(kp) + sgn(k_use)*delt_e(k_use)*astep_mc(kt,kf)*dtsc  
                      pphfbm(kp)=pphfbm(kp)+sgn(k_use)*delt_p(k_use)*astep_mc(kt,kf)*dtsc
                   endif
                enddo
             endif
	  
          enddo  ! kt, loop over mc time steps
!c     -------------------------------------------------------------------	
          ! make sure everything is ok and everybody is happy       
          ! if something is not ok with new orbit, try new mc run
          
          ! first, check final energy, make sure it is >0
          if(ntry_evlv.le.ntry_evlv_max.and.efbm(kp).lt.0.0_fp) then   
             e_fbm=e_orig
             efbm(kp)=e_orig
             pphfbm(kp)=pph_orig
             thfbm(kp)=th_orig
             polfbm(kp)=pol_orig
             mufbm(kp)=mu_orig
             bfbm(kp)=b_orig
             newvalfbm(kp)=0
             ntry_evlv=ntry_evlv+1
             ntry_evlv_p=ntry_evlv
             goto 111 ! try again
          endif
          
          if (abs(efbm(kp)-e_orig).gt.tol_max.or.abs(pphfbm(kp)-pph_orig).gt.tol_max) then 
             ! estimate new psi based on (e,pz) kicks
             dum5=pfbm(kp)
             dum3=-pphfbm(kp)+pph_orig
             dum4=max(1.e-4_fp,pol_orig+dum3) ! limit psi_guess to [1d-4,0.99]
             dum4=min(0.99_fp,dum4)
             psi_guess=dum4
             ! new, sep. 2015 - mp
             ! weight psi_guess on original pitch,
             ! avoid wrong guess due to ~1/p^2 term in dum2
             f_damp=1.0_fp-exp(-(pfbm(kp)*pfbm(kp))/1.0e-3_fp)
             psi_guess=f_damp*dum4+(1.0_fp-f_damp)*pol_orig
             psi_guess=max(1.e-4_fp,psi_guess) ! limit psi_guess to [1d-4,0.99]
             psi_guess=min(0.99_fp,psi_guess)
             ! end modified code, sep. 2015

             ! check whether new orbit is expected to
             ! be "unphysical" [not preserving sign of the pitch]
             ! based on *estimated* new psi_pol. if so, skip
             ! the time-consuming 'orb_update' step and move
             ! on to a new mc attempt.
	! -----------------------------------------
	! note, sep. 15 - mp:
	! skip this check, as it may lead to
	! discarding 'good' orbits when the 
	! estimated psi_guess is far from the actual psi.
	! this may especially occurr for large kicks, since
	! psi_guess is based on a *linear* expansion and
	! neglects changes in b(psi) and g(psi).
!             if(pfbm(kp)*dsign(1.0_fp,pphfbm(kp)+psi_guess).lt.0.0_fp) then
!                is_ok=0 ! new pitch likely to be wrong, failure expected
!             else  
             call orb_update(kp,psi_guess,is_ok)
                ! update flag, values have been modified
             newvalfbm(kp)=1 
 
!             endif
          else ! nothing has changed
             is_ok=1	
             newvalfbm(kp)=0
          endif

             
          ! if something is not ok with new orbit, try new mc run
          if(is_ok.eq.0.and.ntry_evlv.le.ntry_evlv_max) then 
             e_fbm=e_orig
             efbm(kp)=e_orig
             pphfbm(kp)=pph_orig
             thfbm(kp)=th_orig
             polfbm(kp)=pol_orig
             mufbm(kp)=mu_orig
             bfbm(kp)=b_orig
             newvalfbm(kp)=0
             ntry_evlv=ntry_evlv+1
             ntry_evlv_p=ntry_evlv
             goto 111 ! try again
             
          endif
          
          if (is_ok.eq.0.and.ntry_evlv.gt.ntry_evlv_max) then ! no hope... 
             e_fbm=e_orig
             efbm(kp)=e_orig
             pphfbm(kp)=pph_orig
             thfbm(kp)=th_orig
             polfbm(kp)=pol_orig
             mufbm(kp)=mu_orig
             bfbm(kp)=b_orig
             newvalfbm(kp)=0
             call errorlog(' evolve_mc: faild to evolve ptcl.,efbm= ', e_orig)
             call errorlog(' evolve_mc: faild to evolve ptcl.,polfbm= ', pol_orig)
             
          endif
       endif ! energy within allowed range
       
764    continue 
       
       if (abs(efbm(kp)-e_orig).gt.tol_max.or.abs(pphfbm(kp)-pph_orig).gt.tol_max) then 
          k_ptcl=id_ptcl
          if(il_first) then
             naeptot(lsbeam)=naeptot(lsbeam)+1
             il_first=.false.
          endif
          !c     velocity, |v|      
          !      vmod=sqrt(2.e+3_fp/(m2cf*mass2charge)*edum)
          !mg use relation:  v2toev*1.e4_fp = .5 * m2cf*mass2charge
          !      vmod=sqrt(1.e-1_fp/v2toev*edum)
          xiayp=xion
          thp=th
          rmjionp=rmjion
          yionp=yion
          call ptclsave(xiayp,th,vion,xksid,rmjion,yion,ztfluxp,zmup,zpphip,zmop)
          zepp=vion**2*v2toev ! ev
          pmechayp=pmechay(k_ptcl,lsbeam)
          xksid=pfbm(kp)!back to nubeam convention
          ngcp=ngc
          
          btgcp=btori
          rmjion=rfbm(kp)*1.0e2_fp !update after aep kick
          yion=zfbm(kp)*1.0e2_fp   !update after aep kick
          call nblocb(rmjion,yion,xion,th,dum1,zbdir,0,0,ier)!get xion,th
          zep=efbm(kp)*1.e+3_fp-&
               aep_uscale*xzbeams(lsbeam)*fnbint(xion,phiprg,1,mj,1,idummy,ngc) ! ev
          if(zep.le.0.0_fp) then
             call errorlog(" nbi_aep_mod: ptcl. energy less than 0, check aep_usclae")
             call bad_exit(" nbi_aep_mod: ptcl. energy less than 0")
          endif
          dum1=sqrt(zep/v2toev) ! cm/s
          vion=dum1 !updated
          !
          dum2=(1.0_fp-xksid*xksid)*zep/v2toev ! (perp. v)^2 [cm^2/s^2]
          dum3=0.5_fp*zmp*2.0_fp/(1.0e+4_fp*bfbm(kp))
          xmuay(k_ptcl,lsbeam)=dum2*dum3
          
          dum2=zmp*2.0_fp*xksid*dum1*1.e+2_fp*rfbm(kp)
          pmechay(k_ptcl,lsbeam)=dum2

          bixi(1) = xion
          bithst(1) = th
!c
!c    recompute the magnetics at the new position
!c
          bimodb(1)=0.0_fp !  zeroed to evaluate b, 1/b from bifourier spline bintrp,  jintrp
                           !  also evaluates derivatives for first step in orbit eqn.
          call nbigc(1)  !  evaluate auxilliary parameters associated 
                         !  with the *guiding center* of a fast ion orbit
          btgci=btori    !  btor after aep kick
          call ptclsave(xion,th,vion,xksid,rmjion,yion,ztflux,zmu,zpphi,zmo) !updated
          rmjionay(k_ptcl,lsbeam)=rmjion
          xzionay(k_ptcl,lsbeam)=yion
          xiay(k_ptcl,lsbeam)=xion
          thay(k_ptcl,lsbeam)=th
          vay(k_ptcl,lsbeam)=vion
          xksidy(k_ptcl,lsbeam)=xksid
          pmechay(k_ptcl,lsbeam)=zmo
          pphiay(k_ptcl,lsbeam)=zpphi
          xmuay(k_ptcl,lsbeam)=zmu
          tfluxay(k_ptcl,lsbeam)=ztflux
          if(bijac(1).le.0.0_fp) then
             ier=1
             call errorlog(' evolve_fbm: bijac <=0')
          endif
!
!  if doing flr, evaluate random gyrophase angle
!
!          if(nlbflr) then
!             call gphase
!          endif

!  find the gyro position in (r,y), (xi,theta (approx)), and indices
! (if not doing flr, copy the g.c. coords.)

          nflrp=nflr
          call nbiflr(dum2,0) !dummy dum2          
          !    radial current sum
          !
          call radcur(xiayp,xion,wght,tirabds)
          !  count  difference as jxb torque for gc, nj should be  = ngc here
          call nb0jxb(wght,zepp,zep,pmechay(k_ptcl,lsbeam),pmechayp, &
               btgcp,btgci,ngcp,ngc,xiayp,xion,0,.false.) !iflag=0 dosen't matter. taken at gc
          
          do i=1,nkaep_use
             k_use=imap_nkaep_use(i,lsbeam)
             if(k_use.eq.0) cycle

             nstep_aep(k_use,lsbeam)=nstep_mc(k_use)
             if(maxval(valid(nstep_mc_in(k_use):nstep_mc_out(k_use),i)).gt.0) then
                kf=k_use

                if(abs(e_fbm(kf)-e_orig).gt.tol_max) then
                   zep=e_fbm(kf)*1.e+3_fp ! ev
                   pbaepms(nflr,kf,lsbeam)=pbaepms(nflr,kf,lsbeam)+(zep-zepp-epot0)*wght
                   if(nltest_output) then
                      dist_e=0.0_fp
                      dist_ph=0.0_fp
                      dist_mu=0.0_fp
                      dist_pol=0.0_fp
                      npe=aepe(kf,1)
                      npp=aepp(kf,1)
                      npm=aepmu(kf,1)
                      pvar_e1=aepe(kf,2)*1.e-3_fp
                      pvar_en=aepe(kf,3)*1.e-3_fp
                      pvar_p1=aepp(kf,2)
                      pvar_pn=aepp(kf,3)
                      pvar_mu1=aepmu(kf,2)
                      pvar_mun=aepmu(kf,3)
                      call get_bin(efbm(kp),int(aepe(kf,1)),pvar_e1,pvar_en,ije(1:2,kf),dist_e(1))
                      if(ije(1,kf).eq.-1) then
                         ije(1:2,kf)=1
                         dist_e(1)=1.0_fp
                      endif
                      dist_e(2)=1.0_fp-dist_e(1)
                      call get_bin(pphfbm(kp),npp,pvar_p1,pvar_pn,ijph(1:2,kf),dist_ph(1))
                      if(ijph(1,kf).eq.-1) then
                         ijph(1:2,kf)=1
                         dist_ph(1)=1.0_fp
                      endif
                      dist_ph(2)=1.0_fp-dist_ph(1)
                      call get_bin(mufbm(kp)/efbm(kp)*baxis,npm,pvar_mu1,pvar_mun,&
                           ijmu(1:2,kf),dist_mu(1))
                      if(ijmu(1,kf).eq.-1) then
                         ijmu(1:2,kf)=1
                         dist_mu(1)=1.0_fp
                      endif
                      dist_mu(2)=1.0_fp-dist_mu(1)
                      call get_bin(polfbm(kp),n_norm_pol,0.0_fp,1.0_fp,&
                           ijpol(1:2,kf),dist_pol(1))
                      if(ijpol(1,kf).eq.-1) then
                         ijpol(1:2,kf)=1
                         dist_pol(1)=1.0_fp
                      endif
                      dist_pol(2)=1.0_fp-dist_pol(1)
                     do ie=1,2
                         do ip=1,2
                            do imu=1,2
                               pbaep_ms(ije(ie,kf),ijph(ip,kf),ijmu(imu,kf),i,1)=pbaep_ms(ije(ie,kf),ijph(ip,kf),ijmu(imu,kf),i,1)& !power
                                    +(zep-zepp-epot0)*wght*dist_e(ie)*dist_ph(ip)*dist_mu(imu)
                               pph_ms(ije(ie,kf),ijph(ip,kf),ijmu(imu,kf),i,1)=pph_ms(ije(ie,kf),ijph(ip,kf),ijmu(imu,kf),i,1)+& !canonical mom
                                    (pph_orig-pphfbm(kp))*wght*dist_e(ie)*dist_ph(ip)*dist_mu(imu)
                               aep_count(ije(ie,kf),ijph(ip,kf),ijmu(imu,kf),i,1)=&
                                    aep_count(ije(ie,kf),ijph(ip,kf),ijmu(imu,kf),i,1)+& ! #of particle
                                    wght*dist_e(ie)*dist_ph(ip)*dist_mu(imu)
                            enddo
                         enddo
                      enddo
                      do ipol=1,2
                         pbaep_pol(ijpol(ipol,kf),i,1)= pbaep_pol(ijpol(ipol,kf),i,1)+& !power
                              (zep-zepp-epot0)*wght*dist_pol(ipol)
                      enddo
                   endif
                endif
             endif
          enddo
       endif
    enddo  ! kp, loop over particles
    
    !mgr    if(allocated(seeds)) deallocate(seeds)
    !    call infolog('extiting evolve_fbm')
    
  end subroutine evolve_fbm

  !  cccccccccccccccccccccccccccccccccccccc
  subroutine sample_dedp(ije,ijph,ijmu,ipde,deout,dpout,pcom_out,samp_err)
!ccc- this routine returns a pair of random kicks
!ccc- in (e,p) based on the p(de,dp|e,p,mu) probability distribution
!ccc- at given (e,p,mu) coordinates (or "bins")
!c============
!c idecl:  explicitize implicit integer declarations:
    use iso_c_binding, only: fp => c_double
    integer samp_err
    integer ije,ijph,ijmu,ipde	! input integer bins in (e,p,mu) phase space
    integer ntry_sampl
    integer ntry_sampl_max ! number of attempt to extract de,dp
    integer npe,npp
!c============
!c idecl:  explicitize implicit real declarations:
    real(fp) :: deout,dpout ! output kicks in e,pphi
    real(fp) :: pcom_out ! skewness factor
    real(fp) :: ftrgt,fval1,fval2 ! interpolated values of p(de,dp)
    real(fp) :: max_ice
    real(fp) :: de_strt,dp_strt,de_range,dp_range
    real(fp) :: rnd1(1),rnd2(2),ranx	! random number(s)
    integer jemn,jemx,jpmn,jpmx
!c============
    integer :: idum,i_use,jde,jde_max,jdp,jdp_max
    integer, dimension(mnkaep),save :: ije_old,ijph_old,ijmu_old,ipde_old
!c============

    if(.not.il_old) then ! zero first time
     ije_old=0
     ijph_old=0
     ijmu_old=0
     ipde_old=0 
     il_old=.true.
    endif
    ntry_sampl=0
    samp_err=1 ! error flag
    
    i_use=imap_nkaep_use(ipde,lsbeam)
    if(i_use.eq.0) return
    npe=aepde(i_use,1)
    npp=aepdp(i_use,1)
    ntry_sampl_max=50*npe*npp
    ! if only one bin in aepd with non zero value then return
    max_ice=0.0_fp
    do jde=1,npe
       do jdp=1,npp
	  ! modified 05/09/2014
          if (aepd(jde,jdp,ije,ijph,ijmu,i_use).gt.max_ice) then
             max_ice=aepd(jde,jdp,ije,ijph,ijmu,i_use)
             jde_max=jde
             jdp_max=jdp
          endif
       enddo
    enddo ! jde  
!c ----------------------------------------------------------          
!c     check whether p(de,dp) does not contain valid points.
!c     if it doesn't, just return de=0,dpz=0 and p_com=0.5.
!c     particles for which p(de,dp) has no finite values
!c     will not receive kicks.
    if (max_ice.eq.0.) then
       deout=0.0_fp
       dpout=0.0_fp
       pcom_out=0.5_fp
       samp_err=0
       return     
    endif
!c ----------------------------------------------------------          
!c     if p(de,dp) contains only a single
!c     value at de=dp=0, compute kicks' value directly.
!c     this requires the indexes jde_max,jdp_max=jdp
!c     computed above.
    if (max_ice.ge.sum(aepd(1:npe,1:npp,ije,ijph,ijmu,i_use))) then
       deout=spline2_o(ipde)%x1(jde_max)
       dpout=spline2_o(ipde)%x2(jdp_max)
       pcom_out=0.0_fp ! complementary probability will always be zero here
       samp_err=0
      return     
    endif
! ----------------------------------------------------------            
!     look for de,dpz range over which p(de,dp) has finite
!     [i.e., non-zero] values to speed up mc sampling below.
      
! ----------------------------------      
!     new code - june 2015, mp
    jemn=-1
    jemx=npe+1
    jpmn=-1
    jpmx=npp+1      
    do jde=1,npe
       if(maxval(aepd(jde,1:npp,ije,ijph,ijmu,i_use)).gt.0.and.&
            &     jemn.lt.0) jemn=jde
       if(maxval(aepd(npe-jde+1,1:npp,ije,ijph,ijmu,i_use)).gt.0.and.&
            &     jemx.gt.npe) jemx=npe-jde+1
    enddo
    do jdp=1,npp
       if(maxval(aepd(1:npe,jdp,ije,ijph,ijmu,i_use)).gt.0.and.&
            &     jpmn.lt.0) jpmn=jdp
       if(maxval(aepd(1:npe,npp-jdp+1,ije,ijph,ijmu,i_use)).gt.0.and.&
            &     jpmx.gt.npp) jpmx=npp-jdp+1
    enddo
    de_range=(aepde(i_use,3)-aepde(i_use,2))*1.e-3_fp/(npe-1)
    dp_range=(aepdp(i_use,3)-aepdp(i_use,2))/(npp-1)
    de_strt=aepde(i_use,2)*1.e-3_fp+(jemn-1)*de_range
    dp_strt=aepdp(i_use,2)+(jpmn-1)*dp_range
    de_range=(jemx-jemn)*de_range
    dp_range=(jpmx-jpmn)*dp_range
!     end new code - june 2015, mp
! ----------------------------------
            

    ! set interpolating function
    if(ije_old(ipde).ne.ije.or.ijph_old(ipde).ne.ijph.or.ijmu_old(ipde).ne.ijmu.or.ipde_old(ipde).ne.ipde) then
       sp_f2(1:npe,1:npp,ipde)=aepd(1:npe,1:npp,ije,ijph,ijmu,i_use)
       ! setting up spline/linear interpolation
       if(maxval(sp_f2(1:npe,1:npp,ipde)).gt.0) then
          ier_p=0
          call ezspline_setup(spline2_o(ipde),sp_f2(1:npe,1:npp,ipde),ier_p)
          call ezspline_error(ier_p)
          if(ier_p.ne.0)then
             call errorlog(' sample_dedp failed in ezspline_setup')
             deout=0.0_fp
             dpout=0.0_fp
             pcom_out=0.50_fp
             samp_err=1
             return  
          endif
          ije_old(ipde)=ije
          ijph_old(ipde)=ijph
          ijmu_old(ipde)=ijmu
          ipde_old(ipde)=ipde
       else
          deout=0.0_fp
          dpout=0.0_fp
          pcom_out=0.50_fp
          samp_err=1
          return  
       endif
    endif
    
    !c     start mc random sampling with rejection technique
       
       ! extract "target" value for p(de,dp)
    do while (ntry_sampl.le.ntry_sampl_max.and.samp_err.eq.1)
       ftrgt=-1.0_fp
       ntry_sampl=ntry_sampl+1
       
       ! extract random kicks de,dp
!mgr                 call random_number(rnd1(1))
!mgr                 call errorlog(' 1 rnd1=',rnd1(1))
!mgr                 call random_number(rnd2)
!mgr                 call errorlog(' 1 rnd2(1)=',rnd2(1))
!mgr                 call errorlog(' 1 rnd2(2)=',rnd2(2))
       do while(ftrgt.le.0.0_fp)
          rnd1(1)=nbrandm_aep(idum)
          ftrgt=(1.0_fp-rnd1(1))*max_ice
       enddo
       rnd2(1)=nbrandm_aep(idum)
       rnd2(2)=nbrandm_aep(idum)
! ----------------------------------      
!     new code - june 2015, mp
        deout=de_strt+rnd2(1)*de_range
        dpout=dp_strt+rnd2(2)*dp_range
! ----------------------------------      
!     end new code - june 2015, mp
!       deout=spline2_o(ipde)%x1(1)+rnd2(1)*(spline2_o(ipde)%x1(npe)-spline2_o(ipde)%x1(1))
!       dpout=spline2_o(ipde)%x2(1)+rnd2(2)*(spline2_o(ipde)%x2(npp)-spline2_o(ipde)%x2(1))
       ! compute p(de,dp) for random delt_e, delt_p point
       
       ! point interpolation    
       call ezspline_interp(spline2_o(ipde),deout,dpout,fval1,ier_p)
       fval1=max(0.0_fp,fval1) ! exclude error due to interpolation
       ! check for errors
       call ezspline_error(ier_p)
       if(ier_p.ne.0)then
          call errorlog(' sample_dedp failed in ezspline_interp')
          deout=0.0_fp
          dpout=0.0_fp
          pcom_out=0.50_fp
          samp_err=1
          return  
       endif
       
       if (fval1.ge.ftrgt) then
          samp_err=0 ! success
       endif
       
    enddo
    
    !c ------------------------------------------------------------------ 
    
    if (ntry_sampl.gt.ntry_sampl_max) then
       deout=0.0_fp
       dpout=0.0_fp
       pcom_out=0.50_fp
       samp_err=1
       return  
    endif
    
    ! point interpolation for complementary probability
    call ezspline_interp(spline2_o(ipde),-deout,-dpout,fval2,ier_p)
    fval2=max(0.0_fp,fval2) ! exclude error due to interpolation
    call ezspline_error(ier_p)
    if(ier_p.ne.0)then
       call errorlog(' sample_dedp failed in ezspline_interp for -deout')
       deout=0.0_fp
       dpout=0.0_fp
       pcom_out=0.50_fp
       samp_err=1
       return  
    endif
    
!c     ----------------------------------------------------------      

!c     make sure that both probabilities have finite value.
!c     if not, repeat sampling - this may be an indication of poor
!c     statistics in the p(de,dp) bin.     
!c      if(fval1.eq.0.d0.or.fval2.eq.0.d0) goto 997
      ! note: it was 'goto 998' until nov. 21, 2013 - mp
      
      
!cc    make sure that at least one probability value is finite.
!cc    if not, repeat sampling     
    if(fval1.eq.0.0_fp.and.fval2.eq.0.0_fp) then
       deout=0.0_fp
       dpout=0.0_fp
       pcom_out=0.50_fp
       samp_err=1
       return  
    endif
   
!c     compute skewness factor p_com
    pcom_out = fval2/(fval1+fval2)
    samp_err=0
    return
    
  end subroutine sample_dedp

  !  cccccccccccccccccccccccccccccccccccccc

  subroutine sample_dedp_cell(ije,ijph,ijmu,dist,ipde,deout,dpout,pcom_out,samp_err)
!ccc- this routine returns a pair of random kicks
!ccc- in (e,p) based on the p(de,dp|e,p,mu) probability distribution
!ccc- at given (e,p,mu) coordinates (or "bins")
!c============
!c idecl:  explicitize implicit integer declarations:
    use iso_c_binding, only: fp => c_double
    integer samp_err
    integer,intent(in),dimension(2) :: ije,ijph,ijmu	!  nearest integer bin in (e,p,mu) phase space
    real(fp),intent(in),dimension(3) :: dist !dist_e,dist_ph,dist_mu
    integer,intent(in) :: ipde
    integer ntry_sampl
    integer ntry_sampl_max ! number of attempt to extract de,dp
    integer npe,npp
!c============
!c idecl:  explicitize implicit real declarations:
    real(fp) :: deout,dpout ! output kicks in e,pphi
    real(fp) :: pcom_out ! skewness factor
    real(fp) :: ftrgt,fval1,fval2 ! interpolated values of p(de,dp)
    real(fp) :: max_ice
    real(fp) :: de_strt,dp_strt,de_range,dp_range
    real(fp) :: rnd1(1),rnd2(2),ranx	! random number(s)
    integer jemn,jemx,jpmn,jpmx
!c============
    integer :: idum,i_use,jde,jde_max,jdp,jdp_max
    integer, dimension(mnkaep),save :: ije_old,ijph_old,ijmu_old,ipde_old
    integer :: ict(8) !output request
    real(fp) :: hx,hy,hz
    logical :: il_e_ph_mu, il_e_ph, il_e_mu, il_e, il_ph_mu, il_ph, il_mu, il
!c============
    il_e_ph_mu=.false.
    il_e_ph=.false.
    il_e_mu=.false.
    il_e=.false.
    il_ph_mu=.false.
    il_ph=.false.
    il_mu=.false.
    il=.false.
    if(.not.il_old) then ! zero first time
     ije_old=0
     ijph_old=0
     ijmu_old=0
     ipde_old=0 
     il_old=.true.
    endif
    ntry_sampl=0
    samp_err=1 ! error flag
    
    i_use=imap_nkaep_use(ipde,lsbeam)
    if(i_use.eq.0) return
    npe=aepde(i_use,1)
    npp=aepdp(i_use,1)
    ntry_sampl_max=50*npe*npp
    ! if only one bin in aepd with non zero value then return
    max_ice=0.0_fp
    ict=0
    ict(1)=1 !get function only
    hx=(aepe(i_use,3)-aepe(i_use,2))*1.e-3_fp/aepe(i_use,1)
    hy=(aepp(i_use,3)-aepp(i_use,2))/aepp(i_use,1)
    hz=(aepmu(i_use,3)-aepmu(i_use,2))/aepmu(i_use,1)
    if(ije(1).eq.ije(2)) then
       if(ijph(1).eq.ijph(2)) then
          if(ijmu(1).eq.ijmu(2)) then
             il_e_ph_mu=.true.
          else
             il_e_ph=.true.
          endif
       else
          if(ijmu(1).eq.ijmu(2)) then
             il_e_mu=.true.
          else
             il_e=.true.
          endif
       endif
    else
       if(ijph(1).eq.ijph(2)) then
          if(ijmu(1).eq.ijmu(2)) then
             il_ph_mu=.true.
          else
             il_ph=.true.
          endif
       else
          if(ijmu(1).eq.ijmu(2)) then
             il_mu=.true.
          else
             il=.true.
          endif
       endif
    endif
    do jde=1,npe
       do jdp=1,npp
          if(il_e_ph_mu) then
             aepd2(jde,jdp)= aepd(jde,jdp,ije(1),ijph(1),ijmu(1),i_use)
          endif
          if(il_e_ph) then
             call pc1fcn(ict,1,1, &
               aepd2(jde,jdp),1,dist(3),hz,1.0_fp, &
               aepd(jde,jdp,ije(1),ijph(1),ijmu(1):ijmu(2),i_use),2)
          endif
          if(il_e_mu) then
             call pc1fcn(ict,1,1, &
               aepd2(jde,jdp),1,dist(2),hy,1.0_fp, &
               aepd(jde,jdp,ije(1),ijph(1):ijph(2),ijmu(1),i_use),2)
          endif
          if(il_e) then
             call pc2fcn(ict,1,1, &
                  aepd2(jde,jdp),1,1,dist(2),dist(3), &
                  hy,1.0_fp,hz,1.0_fp, &
                  aepd(jde,jdp,ije(1),ijph(1):ijph(2),ijmu(1):ijmu(2),i_use),2,2)
          endif
          if(il_ph_mu) then
             call pc1fcn(ict,1,1, &
                  aepd2(jde,jdp),1,dist(1),hx,1.0_fp, &
                  aepd(jde,jdp,ije(1):ije(2),ijph(1),ijmu(1),i_use),2)
          endif
          if(il_ph) then
             call pc2fcn(ict,1,1, &
                  aepd2(jde,jdp),1,1,dist(1),dist(3), &
                  hx,1.0_fp,hz,1.0_fp, &
                  aepd(jde,jdp,ije(1):ije(2),ijph(1),ijmu(1):ijmu(2),i_use),2,2)
          endif
          if(il_mu) then
             call pc2fcn(ict,1,1, &
                  aepd2(jde,jdp),1,1,dist(1),dist(2), &
                  hx,1.0_fp,hy,1.0_fp, &
                  aepd(jde,jdp,ije(1):ije(2),ijph(1):ijph(2),ijmu(1),i_use),2,2)
          endif
          if(il) then
             call pc3fcn(ict,1,1, &
                  aepd2(jde,jdp),1,1,1,dist(1),dist(2),dist(3), &
                  hx,1.0_fp,hy,1.0_fp,hz,1.0_fp, &
                  aepd(jde,jdp,ije(1):ije(2),ijph(1):ijph(2),ijmu(1):ijmu(2),i_use),2,2,2)
          endif
	  ! modified 05/09/2014
          if (aepd2(jde,jdp).gt.max_ice) then
             max_ice=aepd2(jde,jdp)
             jde_max=jde
             jdp_max=jdp
          endif
       enddo
    enddo ! jde  
!c ----------------------------------------------------------          
!c     check whether p(de,dp) does not contain valid points.
!c     if it doesn't, just return de=0,dpz=0 and p_com=0.5.
!c     particles for which p(de,dp) has no finite values
!c     will not receive kicks.
    if (max_ice.eq.0.) then
       deout=0.0_fp
       dpout=0.0_fp
       pcom_out=0.5_fp
       samp_err=0
       return     
    endif
!c ----------------------------------------------------------          
!c     if p(de,dp) contains only a single
!c     value at de=dp=0, compute kicks' value directly.
!c     this requires the indexes jde_max,jdp_max=jdp
!c     computed above.
    if (max_ice.ge.sum(aepd2(1:npe,1:npp))) then
       deout=spline2_o(ipde)%x1(jde_max)
       dpout=spline2_o(ipde)%x2(jdp_max)
       pcom_out=0.0_fp ! complementary probability will always be zero here
       samp_err=0
      return     
    endif
! ----------------------------------------------------------            
!     look for de,dpz range over which p(de,dp) has finite
!     [i.e., non-zero] values to speed up mc sampling below.
      
! ----------------------------------      
!     new code - june 2015, mp
    jemn=-1
    jemx=npe+1
    jpmn=-1
    jpmx=npp+1      
    do jde=1,npe
       if(maxval(aepd2(jde,1:npp)).gt.0.and.&
            &     jemn.lt.0) jemn=jde
       if(maxval(aepd2(npe-jde+1,1:npp)).gt.0.and.&
            &     jemx.gt.npe) jemx=npe-jde+1
    enddo
    do jdp=1,npp
       if(maxval(aepd2(1:npe,jdp)).gt.0.and.&
            &     jpmn.lt.0) jpmn=jdp
       if(maxval(aepd2(1:npe,npp-jdp+1)).gt.0.and.&
            &     jpmx.gt.npp) jpmx=npp-jdp+1
    enddo
    de_range=(aepde(i_use,3)-aepde(i_use,2))*1.e-3_fp/(npe-1)
    dp_range=(aepdp(i_use,3)-aepdp(i_use,2))/(npp-1)
    de_strt=aepde(i_use,2)*1.e-3_fp+(jemn-1)*de_range
    dp_strt=aepdp(i_use,2)+(jpmn-1)*dp_range
    de_range=(jemx-jemn)*de_range
    dp_range=(jpmx-jpmn)*dp_range
!     end new code - june 2015, mp
! ----------------------------------
            

    ! set interpolating function
    if(ije_old(ipde).ne.ije(1).or.ijph_old(ipde).ne.ijph(1).or.ijmu_old(ipde).ne.ijmu(1).or.ipde_old(ipde).ne.ipde) then
       sp_f2(1:npe,1:npp,ipde)=aepd2(1:npe,1:npp)
       ! setting up spline/linear interpolation
       if(maxval(sp_f2(1:npe,1:npp,ipde)).gt.0) then
          ier_p=0
          call ezspline_setup(spline2_o(ipde),sp_f2(1:npe,1:npp,ipde),ier_p)
          call ezspline_error(ier_p)
          if(ier_p.ne.0)then
             call errorlog(' sample_dedp failed in ezspline_setup')
             deout=0.0_fp
             dpout=0.0_fp
             pcom_out=0.50_fp
             samp_err=1
             return  
          endif
          ije_old(ipde)=ije(1)
          ijph_old(ipde)=ijph(1)
          ijmu_old(ipde)=ijmu(1)
          ipde_old(ipde)=ipde
       else
          deout=0.0_fp
          dpout=0.0_fp
          pcom_out=0.50_fp
          samp_err=1
          return  
       endif
    endif
    
    !c     start mc random sampling with rejection technique
       
       ! extract "target" value for p(de,dp)
    do while (ntry_sampl.le.ntry_sampl_max.and.samp_err.eq.1)
       ftrgt=-1.0_fp
       ntry_sampl=ntry_sampl+1
       
       ! extract random kicks de,dp
!mgr                 call random_number(rnd1(1))
!mgr                 call errorlog(' 1 rnd1=',rnd1(1))
!mgr                 call random_number(rnd2)
!mgr                 call errorlog(' 1 rnd2(1)=',rnd2(1))
!mgr                 call errorlog(' 1 rnd2(2)=',rnd2(2))
       do while(ftrgt.le.0.0_fp)
          rnd1(1)=nbrandm_aep(idum)
          ftrgt=(1.0_fp-rnd1(1))*max_ice
       enddo
       rnd2(1)=nbrandm_aep(idum)
       rnd2(2)=nbrandm_aep(idum)
! ----------------------------------      
!     new code - june 2015, mp
        deout=de_strt+rnd2(1)*de_range
        dpout=dp_strt+rnd2(2)*dp_range
! ----------------------------------      
!     end new code - june 2015, mp
!       deout=spline2_o(ipde)%x1(1)+rnd2(1)*(spline2_o(ipde)%x1(npe)-spline2_o(ipde)%x1(1))
!       dpout=spline2_o(ipde)%x2(1)+rnd2(2)*(spline2_o(ipde)%x2(npp)-spline2_o(ipde)%x2(1))
       ! compute p(de,dp) for random delt_e, delt_p point
       
       ! point interpolation    
       call ezspline_interp(spline2_o(ipde),deout,dpout,fval1,ier_p)
       fval1=max(0.0_fp,fval1) ! exclude error due to interpolation
       ! check for errors
       call ezspline_error(ier_p)
       if(ier_p.ne.0)then
          call errorlog(' sample_dedp failed in ezspline_interp')
          deout=0.0_fp
          dpout=0.0_fp
          pcom_out=0.50_fp
          samp_err=1
          return  
       endif
       
       if (fval1.ge.ftrgt) then
          samp_err=0 ! success
       endif
       
    enddo
    
    !c ------------------------------------------------------------------ 
    
    if (ntry_sampl.gt.ntry_sampl_max) then
       deout=0.0_fp
       dpout=0.0_fp
       pcom_out=0.50_fp
       samp_err=1
       return  
    endif
    
    ! point interpolation for complementary probability
    call ezspline_interp(spline2_o(ipde),-deout,-dpout,fval2,ier_p)
    fval2=max(0.0_fp,fval2) ! exclude error due to interpolation
    call ezspline_error(ier_p)
    if(ier_p.ne.0)then
       call errorlog(' sample_dedp failed in ezspline_interp for -deout')
       deout=0.0_fp
       dpout=0.0_fp
       pcom_out=0.50_fp
       samp_err=1
       return  
    endif
    
!c     ----------------------------------------------------------      

!c     make sure that both probabilities have finite value.
!c     if not, repeat sampling - this may be an indication of poor
!c     statistics in the p(de,dp) bin.     
!c      if(fval1.eq.0.d0.or.fval2.eq.0.d0) goto 997
      ! note: it was 'goto 998' until nov. 21, 2013 - mp
      
      
!cc    make sure that at least one probability value is finite.
!cc    if not, repeat sampling     
    if(fval1.eq.0.0_fp.and.fval2.eq.0.0_fp) then
       deout=0.0_fp
       dpout=0.0_fp
       pcom_out=0.50_fp
       samp_err=1
       return  
    endif
   
!c     compute skewness factor p_com
    pcom_out = fval2/(fval1+fval2)
    samp_err=0
    return
    
  end subroutine sample_dedp_cell
  
  
  subroutine orb_update(k,psi_guess,is_ok)
!ccc-called after a mc step to update particle's variables
!ccc-and make sure the mc step brought the particle
!ccc-into a valid region of phase space, were orbits exist.
!ccc-this routine updates the values of psi,r,z,b
!ccc-given the particle's index. if necessary, it
!ccc-assumes a random distribution in theta.
    use iso_c_binding, only: fp => c_double
    integer,intent(in) :: k
    integer,intent(out) :: is_ok
    integer :: ierr,j
    integer :: ntry_orb,ntry_orb_theta,ntry_orb_theta_max,ier
    integer :: valido(5),isol
    real(fp) :: pdumo(5),rdumo(5),zdumo(5),bdumo(5),ptchdumo(5)
    real(fp) :: thold,thmax,pol_in
    real(fp) :: dsol_t,dsol,pol_max,pol_min
    integer,parameter :: ntry_orb_max=5 ! increase search range at each iteration
    integer :: idum
    real(fp) :: psi_guess
    integer comp_sol
    real(fp) :: dthet,sgnthet,fct_thet
    
    is_ok=0 ! reset exit status flag, is_ok=1 if successfull
    comp_sol=0
    !c ------------------------------------------------------------------  
    !c     preliminary checks:
    
    ! if energy (or radius) are <=0, return - particle is not ok
    ierr=0
    if(efbm(k).le.0.0_fp.or.rfbm(k).le.0.0_fp) then
       is_ok=0
       return
    endif
    ! if particle is already lost (psi>=1), 
    ! then return
    if(polfbm(k).ge.1.0_fp) then
       is_ok=0
!mg       is_ok=1 ! nothing more can be done...
	! update some other parameters
       !	ptchdum=pfbm(k)
!mg       polfbm(k)=1.0_fp
!mg       xdum=1.0_fp
!mg       call get_b(xdum,thfbm(k),bfbm(k)) ! t
       return
    endif
    !c     find minimum & maximum energy over which
    !c     p(de,dp) is 'active' in this time window
    ! number of attempts to vary theta & find new orbit
    ! now depends on the value of |pitch|. this is
    ! done to improve accuracy for particles near a turning point,
    ! such as potatoes or trapped orbits near the tips.
    fct_thet=1.0e-3_fp/(pfbm(k)*pfbm(k)+1.0e-3_fp)
    ntry_orb_theta_max=2*int(fct_thet*3.6e+3_fp*abs(thfbm(k))/twopi)+1
    dthet=2.0_fp*thfbm(k)/ntry_orb_theta_max
    thold=thfbm(k)
    thmax=thfbm(k)
    ntry_orb_theta=1
    ntry_orb=1
    ier=1
    pdumo=0.0_fp
    rdumo=0.0_fp
    zdumo=0.0_fp
    bdumo=0.0_fp
    ptchdumo=0.0_fp
    valido=0

    
    do while(ntry_orb_theta.le.ntry_orb_theta_max.and.ier.gt.0)	  
       
       ntry_orb=1

      ! define search range, initialize counters/flags
       pol_min=0.0_fp
       pol_max=1.0_fp
       pol_in=psi_guess
      
       if(ntry_orb_theta.gt.1) then
	! define sign
          if(mod(ntry_orb_theta,2).eq.1) then
             sgnthet=1.0_fp
          else
             sgnthet=-1.0_fp
          endif
          
          ! vary theta
          thfbm(k)=thfbm(k)+ntry_orb_theta*sgnthet*dthet
          ! limit theta 
          if (thmax.ge.0) then
             thfbm(k)=max(thfbm(k),0.0_fp)
             thfbm(k)=min(thfbm(k),0.5_fp*twopi)
          else
             thfbm(k)=min(thfbm(k),0.0_fp)
             thfbm(k)=max(thfbm(k),-0.5_fp*twopi)
          endif
       endif
431    continue ! entry point for new attempt with updated search range

 	  
       ! move to polfun_ezt function
       ! it has the same input/output parameters as the 
       ! previous polfun_ez, but contains some improvements
       ! to account for possible changes in the sign
       ! of the pitch for trapped, potatoes, stagnation orbits.
       call polfun_ezt(efbm(k),pphfbm(k),mufbm(k),thfbm(k),pol_min,pol_max,&
            pol_in,pfbm(k),pdumo,rdumo,zdumo,bdumo,ptchdumo,valido,comp_sol,ier)
       if (ier.eq.4) then !goto 987 ! unphysical orbit, return w/ error
          ! something wasn't ok - return with error
          is_ok=0
          return     
       endif
       if (ier.eq.0) goto 123 ! search is successfull, mova on
       
      ! if no solution is obtained, but a solution was
      ! bracketed at some point... try again within
      ! updated search range
       if (ier.eq.3.and.ntry_orb.lt.ntry_orb_max) then
          ier=1
          ntry_orb=ntry_orb+1
          pol_in=0.5_fp*(pol_min+pol_max) ! set psi_guess for next iteration
          goto 431
       endif
       ! comment following lines - not useful anymore
      ! avoid [most likely unsuccessful] iterations.
      ! this should mostly apply to trapped/potato/stagnation
      ! or cntr-passing that are put on (e,p,mu) regions with
      ! positive pitch - assumed to be unphysical.
!      if(pfbm(k).lt.0._fp.and.comp_sol.gt.0) then
!         ier=1
!         return
!      endif
      
      ! update counter
      
       ntry_orb_theta=ntry_orb_theta+1
    enddo ! vary theta
    
    
123 continue


    ! select solution that is closer to original psi
    if(ier.eq.0) then
       isol=0
       dsol=1.e+5_fp
       do j=1,5
          dsol_t=sqrt((rdumo(j)-rfbm(k))**2+(zdumo(j)-zfbm(k))**2+(ptchdumo(j)-pfbm(k))**2)
          if(valido(j).eq.1.and.dsol_t.lt.dsol) then
             isol=j
             dsol=dsol_t
          endif
       enddo
    endif
    
    ! if no valid solutions are found, return with
    ! error flag is_ok=0
!?mg    if((ier.ne.0.and.ier.ne.3).or.isol.eq.0) then
    if(ier.ne.0.or.isol.eq.0) then
       ! restore original theta
       thfbm(k)=thold
       is_ok=0
       return
    endif
    
    ! if a valid solution is found, update 
    ! particle's variables
    is_ok=1  
    polfbm(k)=pdumo(isol)
    rfbm(k)=rdumo(isol)
    zfbm(k)=zdumo(isol)
    bfbm(k)=bdumo(isol)
    pfbm(k)=ptchdumo(isol)
    return
  end subroutine orb_update
  
!cccccccccccccccccccccccccccccccccccccc
  subroutine bfield(nvec,psi_in,thet,rhoout,bout,ier)
!c============
!c idecl:  explicitize implicit integer declarations:
    use iso_c_binding, only: fp => c_double
    integer ier,nvec
!c============
!c idecl:  explicitize implicit real declarations:
    real(fp),intent(in) :: psi_in(nvec),thet
    real(fp),intent(out) :: bout(nvec),rhoout(nvec)
    real(fp) :: psi(nvec)
    real(fp) :: rho(nvec) 
    integer :: i
!c============
!ccc-this routine returns the value of the b field
!ccc-for a given psi/psi_w=[0...1] and theta pair.
!ccc-the first step is to get the radius rho that
!ccc-corresponds to the given psi, since this is
!ccc-the variable used in the plasma state profiles.
!ccc-then, the pair (rho,theta) is converted into (r,z).
!ccc-the latter pair is finally used to interpolate
!ccc-the 2d function b(r,z) at the desired location.
! limit psi to avoid bad behavior at lcfs
    ier=0
    bout=abs(baxis)
    rhoout=0.0_fp
    do i=1,nvec
       psi(i)=min(psi_in(i),1.0_fp)
       call ezspline_interp(spl_rhophi,psi(i),rho(i),ier_psi)
       if(ier_psi.ne.0) then
          call ezspline_error(ier_psi)
          call errorlog(' eq_rhopsi in bfield: failed to find rho for psi_pol=',psi(i))
          ier=1
       endif
    enddo
    !c       get magnetic field at required (r,z) location       
    do i=1,nvec
       call get_b(rho(i),thet,bout(i)) ! magnetic field at (rho,theta),  tesla
    enddo
    rhoout=rho
    return    
  end subroutine bfield
!cccccccccccccccccccccccccccccccccccccc

  subroutine polfun_ezt(edum,pz,mudum,thin,polf_min,polf_max,&
       polf_in,ptch_in,pdumo,rdumo,zdumo,bdumo,ptchdumo,valido,comp_sol,ier)
!c idecl:  explicitize implicit real declarations:
!cc find value of pol from (e,pz,mu,theta,pitch).
!cc note: 'polfun_ezt' function [loosely] inspired by
!cc 	 r. white's orbit code (see functions.f file
!cc	 in standard orbit distribution).
!c============
    use iso_c_binding, only: fp => c_double
    real(fp) :: pz,thin,edum,mudum
    real(fp) :: psi1,psi2,psiav,f1,f2,fav,dpsi
    real(fp) :: psq1,psq2,psqav
    real(fp) :: rdum,zdum,pdum,bdum,gdum,dum1
    real(fp) :: pdumo(5),rdumo(5),zdumo(5)
    real(fp) :: ptchdumo(5),bdumo(5)
    real(fp) :: parout(3,5)
    real(fp) :: par_out(5),rnd(1)
    real(fp) :: polf_min,polf_max,polf_in,ptch_in
    real(fp) :: ptchdum1,ptchdum2,ptchdumav
    real(fp) :: tolmin
    real(fp) :: pol2a,pol2b
    integer k,ier,ideb,kdeb,ierr
    integer cnd1,cnd2,comp_sol,id_sol
    integer cnt_lim,valido(5)
    integer retry
    integer :: idum
    real(fp),parameter :: max_tol=1.e-3_fp ! target tolerance for solution, 0<|f|<max_tol
    integer,parameter :: kiter_max=100
    integer,parameter :: n_dpsi=100
    real(fp) :: max_tol2,psiin(3),fout(3),tmp_min,tmp_max,rnd1(1)
    real(fp) :: pol(1),bmax_inv,bdum1(1),bdum2(1),rhowrk(1),max_tol_psi,dpsi_min,ptch
    integer otype ! flag for orbit type

!ccccccccccc-  simplified/faster (?) version of polfun_orb routine

!c --------------------------------------------------------------
!c     first, perform a basic test
    ierr=0
    ier=0
    if(edum.le.0.0_fp) then
       ier = 1	! some error occurred (?)
       return ! unhappy end  
    endif
    max_tol2=max_tol*max_tol
    !c --------------------------------------------------------------
    ier=1
    retry=0
    comp_sol=0 ! flag, set =1 if possible solutions with opposite sign are detected
    id_sol=0 ! flag, set =1 if a solution is bracketed [but not necessarily found]
    ! reset arrays
    do k=1,5 
       valido(k)=0
       pdumo(k)=0.0_fp
       rdumo(k)=0.0_fp
       zdumo(k)=0.0_fp
       ptchdumo(k)=0.0_fp
       bdumo(k)=0.0_fp
    enddo
    !c     set a minimum/maximum psi to avoid issues at
    !c     the magnetic axis and at lcfs
    polf_min=max(polf_min,plflxg(lcentr)/psi_lcfs) 
    polf_max=min(polf_max,1.0_fp) 
!c     compute max allowed b field (e.g. for trapped, stagnation orbits)  
    if(mudum.eq.0.0_fp) then
       return
    endif
    bmax_inv=edum/mudum
      ! get b field at the extrema
    pol(1)=polf_min
    call bfield(1,pol,thin,rhowrk,bdum1,ierr)
    if(ierr.ne.0)then
       ier = 1	! some error occurred (?)
       return
    endif
    pol(1)=polf_max
    call bfield(1,pol,thin,rhowrk,bdum2,ierr)
    if(ierr.ne.0)then
       ier = 1	! some error occurred (?)
       return
    endif
    if (bdum1(1).gt.bmax_inv.and.bdum2(1).gt.bmax_inv) then ! no solutions possible, return
       ier = 1	! some error occurred (?)
       return
    endif
    ! new choice of tolerance on psi - mp, sep. 2015
    max_tol_psi=1.0e-2_fp*polf_in ! target tolerance for solution, |psi2-psi1|<max_tol_psi
!c ------------------------------------------------------------------  
!c     start iterative search for f(e,p,mu)=0	    
    k=1
    dpsi=(polf_max-polf_min)/1.0e2_fp ! improve accuracy near axis, psi<0.1
    dpsi=min(dpsi,1.e-2_fp*sqrt(abs(ptch_in))) ! improve accuracy for trapped, stagnation near |pitch|=0
    dpsi=max(dpsi,1.e-4_fp)	! set lower limit for dpsi initial step
    dpsi_min=dpsi ! new, sep. 2015 - mp
    tolmin=1.e+5_fp
    f1=1.0_fp
    f2=1.0_fp
    fav=1.0_fp
    !psiav=(polf_min+polf_max)/2.0_fp ! mid point
    psiav=polf_in	    
    psi1=max(polf_min,psiav-dpsi)
    psi1=min(polf_max,psi1)
    psi2=min(polf_max,psiav+dpsi)
    psi2=max(polf_min,psi2)
    psqav=-1.0_fp
    cnt_lim=0
    
    do while (k.le.kiter_max)	!.and.tolmin.gt.max_tol.and.psqav.lt.0.0_fp)
       
       ! force early termination with error if
       ! search range is well beyond max input range
       if (dpsi.gt.1.5*(polf_max-polf_min))then
          ier = 1	! some error occurred (?)
          return ! unhappy end  
       endif
       
       ! reset flags
       cnd1=0
       cnd2=0
       
       ! get value of 'f' function at (psi,theta) pairs
       psiin(1)=min(1.0_fp,max(0.0_fp,psi1))
       psiin(2)=min(1.0_fp,max(0.0_fp,psi2))
       psiin(3)=min(1.0_fp,max(0.0_fp,psiav))
       call get_fepmu(edum,pz,mudum,psiin,thin, fout, parout,ierr)
       if(ierr.ne.0)then
          ier = 1	! some error occurred (?)
          return
       endif
       f1=fout(1)
       f2=fout(2)
       fav=fout(3)
       tolmin=abs(fav)
       pdum=psiav
       par_out(1:5)=parout(3,1:5)
       
       ! pitch^2 - must be consistent with definition in get_fepmu subroutine
       psq1 =parout(1,5)
       psq2 =parout(2,5)
       psqav=parout(3,5)
	! compute pitch at three searching locations.
        ! note: here below, the weird sqrt(max(0.,psi)) expression is required
	!       to prevent a 'floating exception' when calling sqrt(0.0) - mp
       ptchdum1=dsign(sqrt(max(0.,psq1)),pz+psi1) ! pitch corresponding to particle at psi=psi1
       ptchdum2=dsign(sqrt(max(0.,psq2)),pz+psi2) ! pitch corresponding to particle at psi=psi2
       ptchdumav=dsign(sqrt(max(0.,psqav)),pz+psiav) ! pitch corresponding to particle at psi=psiav

	! monitor possible presence of solutions with opposite sign of pitch.
	! this info is used in orb_update to skip additional, and likely unsuccessful
	! iterations at smaller theta for some classes of particles.
	if(psq1.gt.0.and.(ptch_in*ptchdum1).lt.0) comp_sol=1
	if(psq2.gt.0.and.(ptch_in*ptchdum2).lt.0) comp_sol=1
	if(psqav.gt.0.and.(ptch_in*ptchdumav).lt.0) comp_sol=1
	
       
       if (tolmin.le.max_tol.and.abs(psi2-psi1).le.max_tol_psi) then
          goto 96 ! solution at psiav is ok, exit from loop
       endif
       if (abs(f1).le.max_tol.and.abs(psiav-psi1).le.max_tol_psi) then
          pdum=psi1
          par_out(1:5)=parout(1,1:5)
	  goto 96 ! solution at psi1 is ok, exit from loop
       endif
       if (abs(f2).le.max_tol.and.abs(psi2-psiav).le.max_tol_psi) then
          pdum=psi2
          par_out(1:5)=parout(2,1:5)
          goto 96 ! solution at psi2 is ok, exit from loop
       endif
       ! else, continue search
	! ---------------------------------------------------
	! check whether a sign inversion is found
	! ---------------------------------------------------
	! new if...then conditions set june 2015 - mp
	! ---------------------------------------------------
	
	
       if (f1*fav.lt.0._fp) cnd1=1 ! root between psi1 and psiav
       if (f2*fav.lt.0._fp) cnd2=1 ! root between psiav and psi2
       
       ! check exit condition |psi2-psi1|<(max_tol)^2 when
       ! a solution is found in [psi1,psi2]
       if(((cnd1.eq.1.and.(psiav-psi1).lt.max_tol_psi).or.&
            &      (cnd2.eq.1.and.(psi2-psiav).lt.max_tol_psi)).and.&
            &       ptch_in*ptchdumav.gt.0) then
          pdum=psiav
          par_out(1:5)=parout(3,1:5)
          goto 96 ! solution at psiav is good enough, exit from loop	
       endif
       ! ---------------------------------------------------
       ! now decide how to update values for next iteration
       
       if (cnd2.eq.1) then ! root in [psiav,psi2]
          id_sol=1
          polf_min=psiav
          polf_max=psi2
          psi1=psiav
          pdum=psiav
          par_out(1:5)=parout(3,1:5)
          psiav=psi2-f2*(psi2-psiav)/(f2-fav)! approx newton's method
          if (psiav.lt.psi1.or.psiav.gt.psi2) then
             psiav=(psi1+psi2)/2.0_fp ! make sure new value is in between psi1,psi2
          endif
          dpsi=.5_fp*(polf_max-polf_min)
          psi2=min(polf_max,psiav+dpsi)
          psi2=max(polf_min,psi2)
          cnt_lim=0
          goto 246
       endif
       if (cnd1.eq.1) then ! root in [psi1,psiav]
          id_sol=1
          polf_min=psi1
          polf_max=psiav
          psi2=psiav	  
          pdum=psiav
          par_out(1:5)=parout(3,1:5)
          psiav=psi1-f1*(psi1-psiav)/(f1-fav)! approx newton's method
          if (psiav.lt.psi1.or.psiav.gt.psi2) then 
             psiav=(psi1+psi2)/2.0_fp ! make sure new value is in between psi1,psi2
          endif
          dpsi=.5_fp*(polf_max-polf_min)
          psi1=max(polf_min,psiav-dpsi)
          psi1=min(polf_max,psi1)
          cnt_lim=0
          goto 246
       endif
       ! if range is not changing for a while, return
       if (psi1.eq.polf_min.and.psi2.eq.polf_max.and.cnt_lim.gt.5) then
          ier = 2	! no solutions for f(e,p,mu)=0 at this theta
          return ! unhappy end   
       endif
       if (psi1.eq.polf_min.and.psi2.eq.polf_max.and.cnt_lim.le.5) then
          cnt_lim=cnt_lim+1
       endif
       if (psi1.eq.psi2.and.cnt_lim.gt.2) then
          ier = 2	! no solutions for f(e,p,mu)=0 at this theta
          return ! unhappy end   
       endif
       if (psi1.eq.psi2.and.cnt_lim.le.2) cnt_lim=cnt_lim+1
       ! else, try to expand range
       dpsi=1.5_fp*dpsi
       psi1=max(polf_min,psiav-dpsi)
       psi1=min(polf_max,psi1)
       psi2=min(polf_max,psiav+dpsi)
       psi2=max(polf_min,psi2)
       psiav=(psi2+psi1)/2.0_fp
       
246    continue	
       if (psi1.eq.polf_min.and.psi2.eq.polf_max.and.cnt_lim.gt.1) then ! try varying psiav
          rnd1(1)=nbrandm_aep(idum)
          psiav=psi1+rnd1(1)*(psi2-psi1)
       endif
      
       k=k+1 ! update counter
       
    enddo
    ! if no solutions are found, but a solution
    ! was bracketed at some point during the 
    ! iterations, exit with flag=3. the orb_update
    ! routine will then decide if it is worth to make
   ! another new attempt with the [polf_min,polf_max]
    if ((k.ge.kiter_max.or.tolmin.gt.max_tol).and.id_sol.eq.1) then
       ier = 3	! no solutions for f(e,p,mu)=0 at this theta, 
      		! but one solution was bracketed at some point.
		! there's still hope for this theta...
       return    ! partially unhappy end 
    endif

    
    ! if no solutions are found, retry or 
    ! return with error flag =1
    if (k.ge.kiter_max.or.tolmin.gt.max_tol) then
       ier = 2	! no solutions for f(e,p,mu)=0 at this theta:
       return    ! unhappy end 
    endif
    
96  continue
      
!c ------------------------------------------------------------------  

!c     now that psi is known, re-obtain (r,z) and pitch
    rdum=par_out(1)
    zdum=par_out(2)
    bdum=par_out(3)
    gdum=par_out(4)
    
!c       must also re-calculate pitch:
    dum1=par_out(5) ! get p^2 from get_fepmu() subroutine	
    ptch=dsign(sqrt(max(0.0_fp,dum1)),pz+pdum) ! picth with sign
    ! find orbit type
    ! note: orb_type function in debug.f file
    call orb_type(edum,pz,mudum,pdum,ptch,otype)
      
    ! --------------------------------------------------
    ! now decide what to do based on new orbit type
    ! and initial pitch.
    if (ptch_in*ptch.lt.0.and. & ! sign of pitch not preserved ...
         &    (otype.eq.1.or.otype.eq.3)) then ! ...and not lost, trapped, stagnation or potato...
	  ! if particle results as 'confined' but has the
	  ! wrong sign of pitch, the solution is assumed
	  ! to be unphysical for some orbit type (co-, cntr-passing).
	  ! for those, the routine will return with error flag =4.
	  ! for other types (trapped, potato, stagnation), 
	  ! there may actually be a change in pitch, as determined
	  ! from orbit simulations. 
	  
       ier=4 ! bad orbit, return
       return
    endif
    valido(1)=1 ! flag 
    ptchdumo(1)=ptch
      
!c     set remaining quantities
    pdumo(1)=pdum
    rdumo(1)=rdum
    zdumo(1)=zdum
    bdumo(1)=bdum 
    gout=gdum !g itself, not  a g/baxis
      
!c ----------------------------------------
100 continue
    ier = 0
!c      write(0,*) '	-> polfun_ezt ok'
    return ! happy end.

            
!c ----------------------------------------        
  end subroutine polfun_ezt

!cccccccccccccccccccccccccccccccccccccc
    subroutine get_fepmu(en,pz,mu,psi_in,thet,fout,parout,ierr)
      use iso_c_binding, only: fp => c_double
      real(fp),intent(in) :: en,pz,mu,psi_in(3),thet
      real(fp),intent(out) :: fout(3),parout(3,5)
      integer,intent(out) :: ierr
      real(fp) :: dum1,epot
      real(fp) :: bdum(3),g_out,rho_out(3)
      real(fp) :: f1,f2,p2out,r,z,dumj(2,2)
      real(fp) :: fnbint
      integer :: i,ierrx,idummy
!ccc-this routine returns the value of the
!ccc-function f(psi,theta|e,pz,mu) for a given (psi,theta) pair

!c ------------------------------------------------------------------
      ierr=0
      ierrx=0
! get b field at this location
      call bfield(3,psi_in,thet,rho_out,bdum,ierr)
      if(ierr.gt.0) then
          parout=0.0_fp
         return
      endif
! get value of 'g' function at this flux surface
      do i=1,3
         call get_g(rho_out(i),g_out) !g_out is not normalized
      ! compute f(psi,theta|e,pz,mu)
!mg         g_out=g_out/baxis ! normalized
!mg         f1=fct_pz*((pz+psi_in(i))*psiw)*bdum(i)/(g_out*baxis)
!mg same as above
         f1=fct_pz*((pz+psi_in(i))*psiw)*bdum(i)/g_out
         epot=xzbeams(lsbeam)*fnbint(rho_out(i),phiprg,1,mj,1,idummy,ngc)*1.e-3_fp ! energy [kev]
         epot=aep_uscale*epot
         dum1=(en-epot)-mu*bdum(i) 	! =>  dum1 = e p^2
!         p2out=dum1/en	! =>  p^2=(e-mub)/e=1-mub/e
         if(en.ne.epot) then
            p2out=1.0_fp-mu*bdum(i)/(en-epot)
         else
            p2out=-1.0_fp
         endif
         if(dum1.ge.0.0_fp) then
            f2=dsign(sqrt(dum1),pz+psi_in(i))  ! => sqrt(e) p, with correct sign
!c	f2=sqrt(dum1)
         else
!c	f2=dsign(sqrt(-dum1),pz+psi)
            f2=0.0_fp
         endif
      
!c      ! test - 04/11/2014, mp
         fout(i)=f1*f1-f2*f2
             
      
      ! extract additional output parameters from bfield.f

         if(nfast_map.eq.1) then
            call nbmomry(rho_out(i),thet,0,r,z,dumj)  ! use bilinear maps
         else
            call nbmomry_xs(rho_out(i),thet,0,r,z,dumj)  !use xplasma spline maps
         endif
         parout(i,1)=r*1.0e-2_fp !  m
         parout(i,2)=z*1.0e-2_fp !  m
         parout(i,3)=bdum(i)		! b
!mg         parout(i,4)=g_out*baxis 	! g
         parout(i,4)=g_out      	! g
        parout(i,5)=p2out		! pitch^2
      enddo
      return     
    end subroutine get_fepmu
!cccccccccccccccccccccccccccccccccccccc

    subroutine get_b(xi,th,bout)
! get b in tesla at a given (xi,th)
      use iso_c_binding, only: fp => c_double
      real(fp), intent(in):: xi !norm sqrt(tor)
      real(fp), intent(in):: th !theta
      real(fp), intent(out):: bout !magnetic field b in [t]

      bixi(1)=xi
      bithst(1)=th
      call bintrp
      bout=1.0_fp/bi1ovb(1) !  b in gauss 
      bout=1.e-4_fp*bout !  b in tesla
    end subroutine get_b

    subroutine get_g(xi,gout)
! get g_mhd at a given xi
      use iso_c_binding, only: fp => c_double
      real(fp), intent(in):: xi !norm sqrt(tor)
      real(fp), intent(out):: gout !magnetic field b in [t]
      bixi(1)=xi
      call gintrp(2)
      gout=big(1) !  b in gauss  r in cm
      gout=1.e-6_fp*gout ! now  b in tesla r in m
    end subroutine get_g

    subroutine create_splpsi(ierr)
      use iso_c_binding, only: fp => c_double
      integer :: ierr,i
      real(fp) :: dpsi
      integer :: bcs_sig(2)
      real(fp), allocatable,dimension(:) :: rho_vals,psi_vals !rho, psi wrk arraies

      ierr=0
      if(allocated(rho_vals)) deallocate(rho_vals,psi_vals)
      allocate(rho_vals(npsi),psi_vals(npsi))
      rho_vals=0.0_fp
      psi_vals(1)=0.0_fp
      dpsi=1.0_fp/(npsi-1)
      do i=2,npsi
         psi_vals(i)=psi_vals(1)+(i-1)*dpsi
      enddo
      psi_vals=psi_vals*psi_lcfs
      call eq_rhopsi(npsi,psi_vals,rho_vals,tol_max,ierr)
!      call eq_rhopsin(npsi,rho_vals,ierr)
      if(ierr.gt.0) then
         call errorlog(' create_splpsi: failed in eq_rhopsi')
         ierr=1
         return
      endif
      call ezspline_free(spl_rhophi,ier_psi)
      bcs_sig = (/0, 0/)     ! not-a-knot boundary conditions
!      call ezlinear_init(spl_rhophi,npsi,ier_psi)
      call ezspline_init(spl_rhophi,npsi,bcs_sig,ier_psi)      ! initialise spline       call ezspline_init(spl_rhophi,npsi,sp_ier) ! piecewise linear 
      call ezspline_error(ier_psi)
      if (ier_psi.ne.0) then 
         ierr=1   ! interrupt routine in case of error
         return
      endif
     spl_rhophi%x1=psi_vals(1:npsi)/psi_lcfs 
     call ezspline_setup(spl_rhophi,rho_vals(1:npsi),ier_psi) ! initialise spline 
     if (ier_psi.ne.0) then   ! interrupt routine in case of error
        ierr=1
        return
     endif
     if(allocated(rho_vals)) deallocate(rho_vals,psi_vals) !free memory
    end subroutine create_splpsi

    subroutine create_spla(ierr)
      use iso_c_binding, only: fp => c_double
      integer :: ierr,i,kf,nptmode
      integer :: bcs_sig(2)

      ierr=0
      if(allocated(spl_amode)) deallocate(spl_amode)
      allocate(spl_amode(1:nkaep_use))
      nptmode=indmax(1)-indmin(1)+1
      do i=1,nkaep_use
         kf=maxval(imap_nkaep_use(i,1:nsbeam))
         call ezspline_free(spl_amode(i),ier_a)
         bcs_sig = (/0, 0/)     ! not-a-knot boundary conditions
!      call ezlinear_init(spl_rhophi,npsi,ier_psi)
         call ezspline_init(spl_amode(i),nptmode,bcs_sig,ier_a)      ! initialise spline       call ezspline_init(spl_rhophi,npsi,sp_ier) ! piecewise linear 
         call ezspline_error(ier_a)
         if (ier_a.ne.0) then 
            ierr=1   ! interrupt routine in case of error
            return
         endif
         spl_amode(i)%x1=aept(indmin(1):indmax(1)) !time [s]
         call ezspline_setup(spl_amode(i),aepa(indmin(1):indmax(1),kf),ier_a) ! initialise spline 
         if (ier_a.ne.0) then   ! interrupt routine in case of error
            ierr=1
            return
         endif
      enddo
    end subroutine create_spla

    subroutine create_splp(ierr)
      use iso_c_binding, only: fp => c_double
      integer :: ierr
      integer i,kf,npe_max,npp_max,jde,jdp
      real(fp) :: pvar_de1,pvar_den,pvar_dp1,pvar_dpn
      integer :: npp_kf,npe_kf
      ier_p=0
      ierr=0
      ! allocate resources
      npe_kf=0
      npp_kf=0
      if(allocated(spline2_o)) deallocate(spline2_o)
      allocate(spline2_o(1:nkaep_use))
      npe_max=0
      npp_max=0
      do i=1,nkaep_use
         kf=maxval(imap_nkaep_use(i,1:nsbeam))
         npe_kf=aepde(kf,1)
         npp_kf=aepdp(kf,1)
         npe_max=max(npe_max,npe_kf)
         npp_max=max(npp_max,npp_kf)
         call ezlinear_init(spline2_o(i),npe_kf,npp_kf,ier_p)
         call ezspline_error(ier_p)
         if(ier_p.ne.0) then
            ierr=1
            return
         endif
         ! set axes
         pvar_de1=aepde(kf,2)*1.e-3_fp
         pvar_den=aepde(kf,3)*1.e-3_fp
         pvar_dp1=aepdp(kf,2)
         pvar_dpn=aepdp(kf,3)
         do jde=1,npe_kf
            spline2_o(i)%x1(jde) =pvar_de1*(npe_kf-jde)/(npe_kf-1)+pvar_den*(jde-1)/(npe_kf-1)
         enddo
         do jdp=1,npp_kf
            spline2_o(i)%x2(jdp) = pvar_dp1*(npp_kf-jdp)/(npp_kf-1)+pvar_dpn*(jdp-1)/(npp_kf-1)
         enddo
      enddo
      if(allocated(sp_f2)) deallocate(sp_f2)
      allocate(sp_f2(npe_max,npp_max,nkaep_use))
      sp_f2=0.0_fp
     
    end subroutine create_splp
    
    subroutine aep_free(ierr)
      use iso_c_binding, only: fp => c_double
      integer :: i,ierr
      ierr=0
      ier_p=0
      ier_psi=0
      if(allocated(spline2_o)) then
         do i=1,nkaep_use
            if(ezspline_allocated(spline2_o(i))) call ezspline_free(spline2_o(i),ier_p)
            call ezspline_error(ier_p)
            if(ier_p.ne.0) then
               ierr=1
               return
            endif
         enddo
         deallocate(spline2_o)
      endif
      
      if(allocated(spl_amode)) then
         do i=1,nkaep_use
            if(ezspline_allocated(spl_amode(i))) call ezspline_free(spl_amode(i),ier_p)
            call ezspline_error(ier_p)
            if(ier_p.ne.0) then
               ierr=1
               return
            endif
         enddo
         deallocate(spl_amode)
      endif

      if(ezspline_allocated(spl_rhophi)) call ezspline_free(spl_rhophi,ier_psi)
      call ezspline_error(ier_psi)
      if(ier_psi.ne.0) then
         ierr=1
         return
      endif

    end subroutine aep_free
!cccccccccccccccccccccccccccccccccccccc
    subroutine orb_type(edum,pdum,mu,pol,ptch_in,otype)
!c============
!cccc-  find orbit type for given (e,pz,mu,psi,p^2)
!c idecl:  explicitize implicit integer declarations:
      use iso_c_binding, only: fp => c_double
      integer otype
!c============
!c idecl:  explicitize implicit real declarations:
      real(fp) :: edum,pdum,gdum,ptch_in
      real(fp) :: dum,dum2,ptch
      real(fp) :: pzdum,pzdum2,elt,ert,eax,etp1,etp2,mu
      real(fp) :: evmin,evmax,ev0,mub
!mg?,elt2
      real(fp) :: pin1
      real(fp) :: pol,bfield_ptc,bfield_0,bfield_pi
      real(fp) :: bpar_out(2)
      real(fp) :: pol_wrk(1),rhowrk(1),bwrk(1)
      integer :: ierr,idummy,ngc_wrk
      real(fp) :: u_max
      real(fp) :: fnbint
  
!c============
!cccccc-  orbit classification.
!ccc   adopt conventions from r. white's orbit code,
!ccc   see rw's book for details.
!ccc   
!ccc   otp=1, co-passing confined, otp = 2, co-passing lost
!ccc   otp=3, ctr-passing confined, otp = 4, ctr-passing lost
!ccc   otp=5, trapped confined, otp = 6, trapped lost
!ccc   otp=7, stagnation, otp = 8, conf potato, otp = 9, trap potato
           
      ! init orbit type as 'unresolved'
      ierr=0
      otype=0
	   
      ! find pitch from input p^2, pz and psi
      ptch=ptch_in
!mg      ptch=dsign(sqrt(max(0.0_fp,psq)),pdum+pol)
	   
	   
           ! define factor for conversion from (pz b)^2/(2 g) -> e [kev]	   
!      fnorm=(zcharge*baxis)**2 
!      fct_pz=fnorm/(atnumb*1d19*zmp/zelonzc) already defined
           
	   	   
!c	   mu = mufbm(k) ! mu=e/b (1-p^2),  given as input
      mub = mu*baxis ! mu x baxis
!c          edum = efbm(k) ! [kev], given as input
!c	   pdum = pphfbm(k) !  given as input	
      pzdum = pdum*psiw	
      pzdum2 = (pdum+pol)*psiw
           
      pin1=max(-pdum,0.0_fp)
      pin1=min(pin1,1.0_fp)

      bwrk=0.0_fp
      rhowrk=0.0_fp
      pol_wrk(1)=pin1
      call bfield(1,pol_wrk,0.0_fp,rhowrk,bwrk,ierr)
      if(ierr.ne.0) then
         call errorlog(' orb_type: cannot get bfield_0 at pol=',pol_wrk(1))
         otype=0
         return
      endif
      bfield_0=bwrk(1)

      bwrk=0.0_fp
      rhowrk=0.0_fp
      pol_wrk(1)=pin1
      call bfield(1,pol_wrk,0.5_fp*twopi,rhowrk,bwrk,ierr)
      if(ierr.ne.0) then
         call errorlog(' orb_type: cannot get bfield_pi at pol=',pol_wrk(1))
         otype=0
         return
      endif
      bfield_pi=bwrk(1)

      u_max=xzbeams(lsbeam)*fnbint(rhowrk(1),phiprg,1,mj,1,idummy,ngc_wrk)*1.e-3_fp ! potential energy [kev]   
      u_max=aep_uscale*u_max
	   
      bwrk=0.0_fp
      rhowrk=0.0_fp
      pol_wrk(1)=pol
      call bfield(1,pol_wrk,0.0_fp,rhowrk,bwrk,ierr)
      if(ierr.ne.0) then
         call errorlog(' orb_type: cannot get bfield_ptc at pol=',pol_wrk(1))
         otype=0
         return
      endif
      bfield_ptc=bwrk(1)

      call get_g(rhowrk(1),gdum) !g_out is not normalized   
	   
       
      dum = (pzdum+psiw)**2/(2.0_fp*gfun_w**2)	 
      ert = fct_pz* dum*b_min_lcfs**2 + mu*b_min_lcfs + u_lcfs
      elt = fct_pz* dum*b_max_lcfs**2 + mu*b_max_lcfs + u_max
    
!did not use      dum = (pzdum2)**2/(2.0_fp*gdum**2)	
!did not use      elt2= fct_pz* dum*b_max_lcfs**2 + mu*b_max_lcfs + u_max !?mg u_pot
      
      dum = pzdum**2/(2.0_fp*gfun_ax**2)
      eax = fct_pz* dum*baxis**2 + mu*baxis + u_ax
      dum2= dum
      
      etp1 = mu*bfield_0 + u_max
      etp2 = mu*bfield_pi + u_max
      
      evmin = mu*b_min_lcfs + u_lcfs
      evmax = mu*b_max_lcfs + u_max
      ev0 = mu*baxis + u_ax
	   	   
!cccccccccccccccccccc-right
      if(pdum.ge.0.0_fp) then
         if(edum.ge.eax) otype = 1 ! co-passing confined
         if(edum.lt.eax) otype = 7 ! stagnation
         goto 999
      endif
!ccccccccccccc-  left
      if(pdum.le.-1.0_fp) then
         if(edum.lt.elt) otype = 4 ! cntr-passing lost
         if(edum.ge.elt) otype = 3 ! cntr-passing confined
         goto 999
      endif
!ccccccccccccccccccccc- middle
      if(pdum.gt.-1.0_fp.and.pdum.lt.0.0_fp) then
         if(edum.ge.elt) otype = 3  ! cntr-passing confined
         !c                 if(edum.gt.eax.and.edum.lt.etp2) otype = 8 ! conf. potato
         !c                 if(edum.gt.eax.and.edum.gt.ert)  otype = 2 ! co-passing lost
         if(edum.lt.eax.and.edum.lt.etp2) otype = 5 ! trapped confined
         if(edum.gt.ert.and.edum.lt.etp2) otype = 6 ! trapped lost
         if(edum.gt.etp2) then ! passing particle
            if(ptch.ge.0.0_fp) then
               if(edum.lt.elt.and.edum.lt.eax) otype = 2 ! co-passing lost
               if(edum.lt.ert.and.edum.lt.eax) otype = 1 ! co-passing confined
            endif
            if(ptch.lt.0.0_fp) then
               if(edum.lt.elt.and.edum.lt.eax) otype = 3 ! cntr-passing confined
            endif
            if(edum.lt.ert.and.edum.gt.eax) otype = 1 ! co-passing confined
         endif
         if(edum.lt.etp2) then
            if(edum.lt.ert.and.edum.lt.eax) otype = 5
         endif
         if(edum.gt.elt) otype = 3 ! cntr-passing confined
      endif
      !cccccccc-  stagnation
      if((edum.lt.etp1.and.edum.lt.eax).or.etp1.lt.u_max) otype = 7 ! stagnation
      !cccccccc- confined potato
      if(edum.gt.eax.and.edum.lt.etp2) otype = 8 ! conf. potato
      !cccccccc- lost potato
      if(edum.gt.ert) then
         if(edum.gt.eax.and.edum.lt.etp2) otype = 9 ! trapped potato
      endif
      !ccccccccccccccccccccccccccccccccccccccccc
      
999   continue   
      
      return
    end subroutine orb_type

  subroutine get_bin(val,ndim,arr1,arrn,indx,dist)
!     returns bin index from bin array arr
    use iso_c_binding, only: fp => c_double
    integer,intent(in) :: ndim !# of grid point
    integer,intent(out) :: indx(2) !returns index
    real(fp), intent(out) :: dist ! define distances from the corner of the bin one in (e,p,mu)  
    real(fp), intent(in) :: val,arr1,arrn !current value, region arr1<=val<=arrn
    integer :: k
    real(fp) darr1n,arrk
      
    indx=-1
    dist=0.0_fp
    darr1n=(arrn-arr1)/(ndim-1)
    do k=1,ndim
       arrk=arr1+darr1n*(k-1)
       if(arrk.le.val) then
          indx(1)=k
          dist=val-arrk
       else
          exit
       endif
    enddo
    if(indx(1).eq.-1) then
       return !error
    endif
    !normilized distance
    dist=dist/darr1n
    indx(2)=min(indx(1)+1,ndim)
    return
  end subroutine get_bin
!ccccccccccccccccccccccccccccccc
    subroutine nbi_aep_clean
      integer :: ierr
      if(allocated(aepsp)) deallocate(aepsp)
      if(allocated(aept)) deallocate(aept)
      if(allocated(aepa)) deallocate(aepa)
      if(allocated(aeptstep)) deallocate(aeptstep)
      if(allocated(aepde)) deallocate(aepde)
      if(allocated(aepdp)) deallocate(aepdp)
      if(allocated(aepe)) deallocate(aepe)
      if(allocated(aepp)) deallocate(aepp)
      if(allocated(aepmu)) deallocate(aepmu)
      if(allocated(aepd)) deallocate(aepd)
      if(allocated(ilaep_init)) deallocate(ilaep_init)
      if(allocated(pbaep_ms)) then
         deallocate(pbaep_ms)
         deallocate(pph_ms)
         deallocate(aep_count)
      endif
      if(allocated(pbaep_pol)) deallocate(pbaep_pol)
      if(allocated(astep_mc)) deallocate(astep_mc)
      if(allocated(tstep_mc)) deallocate(tstep_mc)
      if(allocated(valid)) deallocate(valid)
      if(allocated(imap_nkaep_use)) deallocate(imap_nkaep_use)
      if(allocated(aepd2)) deallocate(aepd2)
      if(allocated(sp_f2)) deallocate(sp_f2)
      call aep_free(ierr)
    end subroutine nbi_aep_clean

    subroutine nbi_aep_write
      use iso_c_binding, only: fp => c_double
      use netcdf
      use execsystem
      implicit none
      integer :: i,kf
      character*2 :: cmode
      character*200 :: filename
      character*50  :: cvar,zlbl
      character*8 :: slabel
      integer :: ierr,icdf,ia,iz,id_tbm1,id_tbm2,ie,ip,im,id_tstart_aep,id_tstop_aep,id_n_norm_pol
      integer, dimension(nkaep_use) :: id_npe,id_npp,id_npmu,id_pvar_e1,id_pvar_en,npe,npp,npmu
      integer, dimension(nkaep_use) :: id_pvar_p1,id_pvar_pn,id_pvar_mu1,id_pvar_mun,id_pb,id_pph,id_cnt,id_ppol,id_nstep_aep
      integer, dimension(3) :: idims
      character*3 isotop(4,2)           ! (a,z)
      real(fp) :: wrk(1)
      data isotop /'h','d','t','?', '?','?','he3','he4'/

      iz=xzbeams(lsbeam) + 0.5_fp
      ia=abeams(lsbeam) + 0.5_fp
      if(nlfprod(lsbeam)) then
         slabel=trim(isotop(ia,iz))//'_fus'
      else
         slabel=trim(isotop(ia,iz))//'_beam'
      endif
      id_npe=0
      id_npp=0
      id_npmu=0

      filename = trim(adjustl(workpath))//&
           trim(adjustl(runid))//'_aep_'//trim(adjustl(slabel))//'.cdf'
      ierr=nf90_open(trim(adjustl(filename)),nf90_write,icdf)
      if(ierr.eq.0) then 
         !file exists
         !remove it
         ierr = nf90_close(icdf)
         ierr =jsystem('rm -rf '//trim(filename))
      endif
      call check(nf90_create(path=trim(adjustl(filename)),cmode=nf90_clobber,ncid=icdf),'create file')
      cvar='n_norm_pol'
      call check(nf90_def_dim(icdf,cvar,n_norm_pol,id_n_norm_pol),cvar)
      do i=1,nkaep_use
         do ie=1,np_e
            do ip=1,np_p
               do im=1,np_m
                  if(aep_count(ie,ip,im,i,1).gt.0) then
                     pph_ms(ie,ip,im,i,1)=pph_ms(ie,ip,im,i,1)/aep_count(ie,ip,im,i,1) !normalize canonical momentum
                  endif
               enddo
            enddo
         enddo
         pbaep_ms(1:np_e,1:np_p,1:np_m,i,1)=pbaep_ms(1:np_e,1:np_p,1:np_m,i,1)*aee/delt !watt
         pbaep_pol(1:n_norm_pol,i,1)=pbaep_pol(1:n_norm_pol,i,1)*aee/delt !watt
         kf=imap_nkaep_use(i,lsbeam)
         if(kf.eq.0) cycle
         cmode=' '
         write(cmode,'(i2)') kf
         npe(i)=int(aepe(kf,1))
         npp(i)=int(aepp(kf,1))
         npmu(i)=int(aepmu(kf,1))
         cvar='npe_'//trim(adjustl(cmode))
         call check(nf90_def_dim(icdf,cvar,npe(i),id_npe(i)),cvar)
!         zlbl='energy linear grid points for mode '//trim(adjustl(cmode))
!         call check(nf90_put_att(icdf,id_npe(i),cvar,trim(zlbl)),cvar)
         cvar='npp_'//trim(adjustl(cmode))
         call check(nf90_def_dim(icdf,cvar,npp(i),id_npp(i)),cvar)
!         zlbl='pphi linear grid points for mode '//trim(adjustl(cmode))
!         call check(nf90_put_att(icdf,id_npp(i),cvar,trim(zlbl)),cvar)
         cvar='npmu_'//trim(adjustl(cmode))
         call check(nf90_def_dim(icdf,cvar,npmu(i),id_npmu(i)),cvar)
!         zlbl='mu linear grid points for mode '//trim(adjustl(cmode))
!         call check(nf90_put_att(icdf,id_npmu(i),cvar,trim(zlbl)),cvar)
         cvar='pvar_e1_mode_'//trim(adjustl(cmode))
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_e1(i)),cvar)
         cvar='pvar_en_mode_'//trim(adjustl(cmode))
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_en(i)),cvar)
         zlbl='min energy for mode '//trim(adjustl(cmode))//' ,[ev]'
         call check(nf90_put_att(icdf,id_pvar_e1(i),cvar,trim(zlbl)),cvar)
         zlbl='max energy for mode '//trim(adjustl(cmode))//' ,[ev]'
         call check(nf90_put_att(icdf,id_pvar_en(i),cvar,trim(zlbl)),cvar)
         cvar='pvar_p1_mode_'//trim(adjustl(cmode))
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_p1(i)),cvar)
         cvar='pvar_pn_mode_'//trim(adjustl(cmode))
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_pn(i)),cvar)
         zlbl='min pphi for mode '//trim(adjustl(cmode))
         call check(nf90_put_att(icdf,id_pvar_p1(i),cvar,trim(zlbl)),cvar)
         zlbl='max pphi for mode '//trim(adjustl(cmode))
         call check(nf90_put_att(icdf,id_pvar_pn(i),cvar,trim(zlbl)),cvar)
         cvar='pvar_mu1_mode_'//trim(adjustl(cmode))
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_mu1(i)),cvar)
         cvar='pvar_mun_mode_'//trim(adjustl(cmode))
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_mun(i)),cvar)
         zlbl='min mu for mode '//trim(adjustl(cmode))
         call check(nf90_put_att(icdf,id_pvar_mu1(i),cvar,trim(zlbl)),cvar)
         zlbl='max mu for mode '//trim(adjustl(cmode))
         call check(nf90_put_att(icdf,id_pvar_mun(i),cvar,trim(zlbl)),cvar)
       
         cvar='pbaep_pol_mode_'//trim(adjustl(cmode))
         call check(nf90_def_var(icdf,cvar,nf90_double,(/ id_n_norm_pol /),id_ppol(i)),cvar)
         zlbl='power  exchange vs norm psi poloidal flux for mode '//trim(adjustl(cmode))//' ,watt'
         call check(nf90_put_att(icdf,id_ppol(i),cvar,trim(zlbl)),cvar)

         idims = (/id_npe(i),id_npp(i),id_npmu(i)/)
         cvar = 'pbaep_'//trim(slabel)//'_mode_'//trim(adjustl(cmode)) 
         call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_pb(i)),cvar)
         zlbl='power due to aep (mhd) activity for mode '//trim(adjustl(cmode))//' ,watt'
         call check(nf90_put_att(icdf,id_pb(i),cvar,trim(zlbl)),cvar)
         cvar = 'pph_'//trim(slabel)//'_mode_'//trim(adjustl(cmode)) 
         call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_pph(i)),cvar)
         zlbl='canonical momentum variation due to aep (mhd) activity for mode '//trim(adjustl(cmode))
         call check(nf90_put_att(icdf,id_pph(i),cvar,trim(zlbl)),cvar)
         cvar = 'aep_count_'//trim(slabel)//'_mode_'//trim(adjustl(cmode)) 
         call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_cnt(i)),cvar)
         zlbl=' number of particles contributing to aep (mhd) activit for mode '//trim(adjustl(cmode))
         call check(nf90_put_att(icdf,id_cnt(i),cvar,trim(zlbl)),cvar)
         cvar='nstep_mc_aep_mode_'//trim(adjustl(cmode))
         call check(nf90_def_var(icdf,cvar,nf90_int,idims,id_nstep_aep(i)),cvar)
         zlbl=' number of steps for aep activity for mode '&
              //trim(adjustl(cmode))
         call check(nf90_put_att(icdf,id_nstep_aep(i),cvar,trim(zlbl)),cvar)
      enddo
      cvar='tbm1'
      call check(nf90_def_var(icdf,cvar,nf90_double,id_tbm1),cvar)
      zlbl='beam start time step'//' ,sec'
      call check(nf90_put_att(icdf,id_tbm1,cvar,trim(zlbl)),cvar)
      cvar='tbm2'
      call check(nf90_def_var(icdf,cvar,nf90_double,id_tbm2),cvar)
      zlbl='beam finish time step'//' ,sec'
      call check(nf90_put_att(icdf,id_tbm2,cvar,trim(zlbl)),cvar)
      cvar='tstart_aep'
      call check(nf90_def_var(icdf,cvar,nf90_double,id_tstart_aep),cvar)
      zlbl='aep start time '//' ,sec'
      call check(nf90_put_att(icdf,id_tstart_aep,cvar,trim(zlbl)),cvar)
      cvar='tstop_aep'
      call check(nf90_def_var(icdf,cvar,nf90_double,id_tstop_aep),cvar)
      zlbl='aep stop time '//' ,sec'
      call check(nf90_put_att(icdf,id_tstop_aep,cvar,trim(zlbl)),cvar)
      ierr = nf90_enddef(icdf) ! out of def mode and go to data mode
       
      cvar='tbm1'
      wrk(1)=tstart_in
      call check(nf90_put_var(icdf,id_tbm1,wrk(1)),cvar)
      cvar='tbm2'
      wrk(1)=tstop_in
      call check(nf90_put_var(icdf,id_tbm2,wrk(1)),cvar)
      cvar='tstart_aep'
      wrk(1)=tstart
      call check(nf90_put_var(icdf,id_tstart_aep,wrk(1)),cvar)
      cvar='tstop_aep'
      wrk(1)=tstop
      call check(nf90_put_var(icdf,id_tstop_aep,wrk(1)),cvar)
      do i=1,nkaep_use
         kf=imap_nkaep_use(i,lsbeam)
         if(kf.eq.0) cycle
         cmode=' '
         write(cmode,'(i2)') kf
         cvar='nstep_mc_aep_mode_'//trim(adjustl(cmode))
         call check(nf90_put_var(icdf,id_nstep_aep(i),nstep_mc(kf)),cvar)
         cvar='pvar_e1_mode_'//trim(adjustl(cmode))
         call check(nf90_put_var(icdf,id_pvar_e1(i),aepe(kf,2)*1.e-3_fp ),cvar)
         cvar='pvar_en_mode_'//trim(adjustl(cmode))
         call check(nf90_put_var(icdf,id_pvar_en(i),aepe(kf,3)*1.e-3_fp ),cvar)
         cvar='pvar_p1_mode_'//trim(adjustl(cmode))
         call check(nf90_put_var(icdf,id_pvar_p1(i),aepp(kf,2) ),cvar)
         cvar='pvar_pn_mode_'//trim(adjustl(cmode))
         call check(nf90_put_var(icdf,id_pvar_pn(i),aepp(kf,3) ),cvar)
         cvar='pvar_mu1_mode_'//trim(adjustl(cmode))
         call check(nf90_put_var(icdf,id_pvar_mu1(i),aepmu(kf,2) ),cvar)
         cvar='pvar_mun_mode_'//trim(adjustl(cmode))
         call check(nf90_put_var(icdf,id_pvar_mun(i),aepmu(kf,3) ),cvar)
         cvar = 'pbaep_'//trim(slabel)//'_mode_'//trim(adjustl(cmode)) 
         call check(nf90_put_var(icdf,id_pb(i),pbaep_ms(1:npe(i),1:npp(i),1:npmu(i),i:i,1:1)),cvar)
         cvar = 'pph_'//trim(slabel)//'_mode_'//trim(adjustl(cmode)) 
         call check(nf90_put_var(icdf,id_pph(i),pph_ms(1:npe(i),1:npp(i),1:npmu(i),i:i,1:1)),cvar)
         cvar = 'aep_count_'//trim(slabel)//'_mode_'//trim(adjustl(cmode)) 
         call check(nf90_put_var(icdf,id_cnt(i),aep_count(1:npe(i),1:npp(i),1:npmu(i),i:i,1:1)),cvar)
         cvar='pbaep_pol_mode_'//trim(adjustl(cmode))
          call check(nf90_put_var(icdf,id_ppol(i),pbaep_pol(1:n_norm_pol,i:i,1:1)),cvar)
     enddo
      !all done
      ierr= nf90_close(icdf)
      
    end subroutine nbi_aep_write

    subroutine nbi_rbq_write
      use iso_c_binding, only: fp => c_double
      use netcdf
      use execsystem
      implicit none
      integer :: i,kf
      character*200 :: filename
      character*50  :: cvar,zlbl
      character*8 :: slabel
      integer :: ierr,icdf,ia,iz,id_tbm1,id_tbm2,ie,ip,im,id_tstart_aep,id_tstop_aep,id_n_norm_pol
      integer, dimension(nkaep_use) :: id_npe,id_npp,id_npmu,id_pvar_e1,id_pvar_en,npe,npp,npmu
      integer, dimension(nkaep_use) :: id_pvar_p1,id_pvar_pn,id_pvar_mu1,id_pvar_mun,id_pb,id_pph,id_cnt,id_ppol,id_nstep_aep
      integer, dimension(3) :: idims
      character*3 isotop(4,2)           ! (a,z)
      real(fp) :: wrk(1)
      data isotop /'h','d','t','?', '?','?','he3','he4'/
      kf=maxval(imap_nkaep_use(:,lsbeam))
      if(kf.ne.1) then
         !RBQ set kf=1 only
         call errorlog(" nbi_rbq_write not wrtiting data due to error: nk not eq 1")
         return
      endif
      iz=xzbeams(lsbeam) + 0.5_fp
      ia=abeams(lsbeam) + 0.5_fp
      if(nlfprod(lsbeam)) then
         slabel=trim(isotop(ia,iz))//'_fus'
      else
         slabel=trim(isotop(ia,iz))//'_beam'
      endif
      id_npe=0
      id_npp=0
      id_npmu=0

      filename = trim(adjustl(workpath))//&
           trim(adjustl(runid))//'_rbq_'//trim(adjustl(slabel))//'.cdf'
      ierr=nf90_open(trim(adjustl(filename)),nf90_write,icdf)
      if(ierr.eq.0) then 
         !file exists
         !remove it
         ierr = nf90_close(icdf)
         ierr =jsystem('rm -rf '//trim(filename))
      endif
      call check(nf90_create(path=trim(adjustl(filename)),cmode=nf90_clobber,ncid=icdf),'create file')
      do i=1,nkaep_use
         do ie=1,np_e
            do ip=1,np_p
               do im=1,np_m
                  if(pph_ms(ie,ip,im,i,1).gt.0) then
                     pbaep_ms(ie,ip,im,i,1)=pbaep_ms(ie,ip,im,i,1)/pph_ms(ie,ip,im,i,1)
                  endif
               enddo
            enddo
         enddo
         ! pph_ms(ie,ip,im,i,1) is a fast ion distribution function in (E, Pphi, lambda)
         ! lambda=(1-xksid^2)*B_axis*R_axis/B = mu/E *B_axis*R_axis/B
         ! Pphi -- noramlizes to Psi_lcfs
         npe(i)=int(aepe(kf,1))
         npp(i)=int(aepp(kf,1))
         npmu(i)=int(aepmu(kf,1))
         cvar='npe'
         call check(nf90_def_dim(icdf,cvar,npe(i),id_npe(i)),cvar)
!         zlbl='energy linear grid points'
!         call check(nf90_put_att(icdf,id_npe(i),cvar,trim(zlbl)),cvar)
         cvar='npp'
         call check(nf90_def_dim(icdf,cvar,npp(i),id_npp(i)),cvar)
!         zlbl='pphi linear grid points'
!         call check(nf90_put_att(icdf,id_npp(i),cvar,trim(zlbl)),cvar)
         cvar='npmu'
         call check(nf90_def_dim(icdf,cvar,npmu(i),id_npmu(i)),cvar)
!         zlbl='mu linear grid points'
!         call check(nf90_put_att(icdf,id_npmu(i),cvar,trim(zlbl)),cvar)
         cvar='pvar_e1'
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_e1(i)),cvar)
         cvar='pvar_en'
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_en(i)),cvar)
         zlbl='min energy'//' ,[ev]'
         call check(nf90_put_att(icdf,id_pvar_e1(i),cvar,trim(zlbl)),cvar)
         zlbl='max energy'//' ,[ev]'
         call check(nf90_put_att(icdf,id_pvar_en(i),cvar,trim(zlbl)),cvar)
         cvar='pvar_p1'
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_p1(i)),cvar)
         cvar='pvar_pn'
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_pn(i)),cvar)
         zlbl='min pphi'
         call check(nf90_put_att(icdf,id_pvar_p1(i),cvar,trim(zlbl)),cvar)
         zlbl='max pphi'
         call check(nf90_put_att(icdf,id_pvar_pn(i),cvar,trim(zlbl)),cvar)
         cvar='pvar_lambda1'
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_mu1(i)),cvar)
         cvar='pvar_lambdan'
         call check(nf90_def_var(icdf,cvar,nf90_double,id_pvar_mun(i)),cvar)
         zlbl='min lambda'
         call check(nf90_put_att(icdf,id_pvar_mu1(i),cvar,trim(zlbl)),cvar)
         zlbl='max lambda'
         call check(nf90_put_att(icdf,id_pvar_mun(i),cvar,trim(zlbl)),cvar)
       
         idims = (/id_npe(i),id_npp(i),id_npmu(i)/)
         cvar = 'polflux_av_'//trim(slabel)
         call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_pb(i)),cvar)
         zlbl='averaged polidal flux / psi_lcfs '
         call check(nf90_put_att(icdf,id_pb(i),cvar,trim(zlbl)),cvar)
         cvar = 'fbm_pph_'//trim(slabel)
         call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_pph(i)),cvar)
         zlbl=' fast ion distribution function in (E, Pphi, lambda)'
         call check(nf90_put_att(icdf,id_pph(i),cvar,trim(zlbl)),cvar)
         cvar = 'rbq_count_'//trim(slabel)
         call check(nf90_def_var(icdf,cvar,nf90_double,idims,id_cnt(i)),cvar)
         zlbl=' number of particles contributing to RBQ (mhd) activity'
         call check(nf90_put_att(icdf,id_cnt(i),cvar,trim(zlbl)),cvar)
      enddo
      cvar='tbm1'
      call check(nf90_def_var(icdf,cvar,nf90_double,id_tbm1),cvar)
      zlbl='beam start time step'//' ,sec'
      call check(nf90_put_att(icdf,id_tbm1,cvar,trim(zlbl)),cvar)
      cvar='tbm2'
      call check(nf90_def_var(icdf,cvar,nf90_double,id_tbm2),cvar)
      zlbl='beam finish time step'//' ,sec'
      call check(nf90_put_att(icdf,id_tbm2,cvar,trim(zlbl)),cvar)
     ierr = nf90_enddef(icdf) ! out of def mode and go to data mode

      cvar='tbm1'
      wrk(1)=tstart_in
      call check(nf90_put_var(icdf,id_tbm1,wrk(1)),cvar)
      cvar='tbm2'
      wrk(1)=tstop_in
      call check(nf90_put_var(icdf,id_tbm2,wrk(1)),cvar)
      do i=1,nkaep_use
         cvar='pvar_e1'
         call check(nf90_put_var(icdf,id_pvar_e1(i),aepe(kf,2)*1.e-3_fp ),cvar)
         cvar='pvar_en'
         call check(nf90_put_var(icdf,id_pvar_en(i),aepe(kf,3)*1.e-3_fp ),cvar)
         cvar='pvar_p1'
         call check(nf90_put_var(icdf,id_pvar_p1(i),aepp(kf,2) ),cvar)
         cvar='pvar_pn'
         call check(nf90_put_var(icdf,id_pvar_pn(i),aepp(kf,3) ),cvar)
         cvar='pvar_lambda1'
         call check(nf90_put_var(icdf,id_pvar_mu1(i),aepmu(kf,2) ),cvar)
         cvar='pvar_lambdan'
         call check(nf90_put_var(icdf,id_pvar_mun(i),aepmu(kf,3) ),cvar)
         cvar = 'fbm_pph_'//trim(slabel)
         call check(nf90_put_var(icdf,id_pph(i),pph_ms(1:npe(i),1:npp(i),1:npmu(i),i:i,1:1)),cvar)
         cvar = 'rbq_count_'//trim(slabel)
         call check(nf90_put_var(icdf,id_cnt(i),aep_count(1:npe(i),1:npp(i),1:npmu(i),i:i,1:1)),cvar)
         cvar='polflux_av_'//trim(slabel)
         call check(nf90_put_var(icdf,id_pb(i),pbaep_ms(1:npe(i),1:npp(i),1:npmu(i),i:i,1:1)),cvar)
      enddo
      !all done
      ierr= nf90_close(icdf)
      
    end subroutine nbi_rbq_write

    subroutine check(status,varnam)
      use netcdf
      implicit none
      integer, intent ( in) :: status 
      character*(*),intent(in) :: varnam ! dimension name input 
      
      if(status /= nf90_noerr) then 
         call errorlog(' nbi_aep/rbq_write: error:'& 
              &// trim(nf90_strerror(status))//' for var.= '//trim(varnam))
         call exit(status) 
      end if
    end subroutine check
  logical function lno_tstep(kt_time,nstep_in,nstep_out)
    integer, intent(in) :: kt_time,nstep_in,nstep_out
    lno_tstep=.false.
    if(nstep_in.eq.0) then
       lno_tstep=.true.
       return
    endif
    if(nstep_out.lt.nstep_in) then
       lno_tstep=.true.
       return
    endif
    if(kt_time.lt.nstep_in.or.kt_time.gt.nstep_out) then
       lno_tstep=.true.
    endif
    return
  end function lno_tstep

  subroutine rbq_init(ierr)
    use physconst_mod , only : mp_g, aee
    use nbi_com, only : nsbeam,plflxg,bmaxm,bixi,big
    use nbi_aep_com, only : nk_aep=>nkaep
    implicit none
    integer :: ierr,k,isb
    logical :: il_match
    real(fp) :: chi_b_min,chi_b_max !dummy
 
    ierr=1
    if(nk_aep.le.0) return
    ierr=0
    baxis=bmaxm(lcentr)
    if(allocated(imap_nkaep_use)) deallocate(imap_nkaep_use)
    allocate(imap_nkaep_use(nk_aep,nsbeam))
    imap_nkaep_use=0
    nkaep_use=0
    do k=1,nk_aep !# of mode
       if(aep_ascale(k).lt.0.0_fp)then
          nkaep_use=nkaep_use+1
          if(int(aepsp(k)).eq.0) then
             imap_nkaep_use(nkaep_use,1:nsbeam)=k
             cycle
          endif
          do isb=1,nsbeam
             call match_fi(k,isb,il_match)
             if(il_match) then
                imap_nkaep_use(nkaep_use,isb)=k
             endif
          enddo
       endif
    enddo
    if (maxval(imap_nkaep_use).eq.0) then
       ierr=1
       return !no RBQ activities
    endif
    psi_lcfs=plflxg(ledge)
    psiw = psi_lcfs/baxis
    bixi(1)=1.0_fp
    call gintrp(2)
    gfun_w=big(1)*1.e-6_fp
    bixi(1)=0.0_fp
    call gintrp(2)
    gfun_ax=big(1)*1.e-6_fp
    call eq_glimb1(1.0_fp,0.0_fp,b_min_lcfs,chi_b_min,b_max_lcfs,chi_b_max,ierr)
    u_ax=0.0_fp
    u_lcfs=0.0_fp
    zmp_zelonzc=mp_g/aee
    fct_pz_xa=baxis*sqrt(0.5_fp/zmp_zelonzc)
    if (plflxg(lcentr).gt.psi_lcfs) psi_lcfs = plflxg(lcentr)
    if(psi_lcfs.eq.0.0_fp) then
       call errorlog(' rbq_test failed:  psi_lcfs=0')
       call bad_exit
    endif
    call create_splpsi(ierr)
    if(ierr.ne.0) then
       call errorlog(' rbq_init failed in create_splpsi')
       call bad_exit
    endif
    if(nltest_output) then
       if(allocated(pbaep_ms)) deallocate(pbaep_ms)
       if(allocated(pph_ms)) deallocate(pph_ms)
       if(allocated(aep_count)) deallocate(aep_count)
       np_e=maxval(aepe(:,1))
       np_p=maxval(aepp(:,1))
       np_m=maxval(aepmu(:,1))
       allocate(pbaep_ms(np_e,np_p,np_m,nkaep_use,1))
       allocate(pph_ms(np_e,np_p,np_m,nkaep_use,1))
       allocate(aep_count(np_e,np_p,np_m,nkaep_use,1))
       pbaep_ms=0.0_fp !av poloidal normalized flux
       pph_ms=0.0_fp     ! FI distribution function in (E,Pphi,lambda)
       aep_count=0.0_fp ! # of particle responded to RBQ diffusion
    endif

   
  end subroutine rbq_init

  subroutine rbq_driver (zdt)
    use nbi_com, only: id_psi,vpll,vperp,bivpar
    use xplasma_definitions
    implicit none
    type (xplasma), pointer :: s
    integer  :: ierr
    real(fp) :: zdt
    integer, parameter :: ion=1,i=1
    real(fp) :: dumj(2,2) !dummy
    integer ::  is_ok, idum
    real(fp) :: diff_sign,zdifpphi,vpll2
    integer :: ie, ipphi, ilambda
    logical :: il_unknown
    integer :: nk
    real(fp) :: dum2,efbmp,pfbmp,bixip,bithstp,bigp,bimodbp,polfbmp,pphfbmp,rfbmp,zfbmp
    if(xion.gt.1.0_fp) return ! do not apply diffusion if ptcl is behind LCFS
    nk=1 ! this dimension for RBQ set to 1 on the input
    il_unknown=.false. ! set to .true. if Pphi for a ptcl not in the given range
    ierr=0
    zdifpphi=0.0_fp
    fct_pz=fct_pz_xa*xzbeams(lsbeam)/sqrt(abeams(lsbeam))
    efbm(i)=v2toev*vion*vion*1.e-3_fp ! energy [kev]
    pfbm(i)=xksid
    thfbm(i)=th    
    bixi(ion)=xion
    efbmp=efbm(i)
    pfbmp=pfbm(i)
    bixip=bixi(ion)
    bithstp=bithst(ion)
    call bintrp
    call gintrp(2)
    bigp=big(1)
    bimodbp=bimodb(1)
    if(nfast_map.eq.1) then
       call nbmomry(xion,th,0, &
            rmjion,yion,dumj)  ! use bilinear maps
    else
       call nbmomry_xs(xion,th,0, &
            rmjion,yion,dumj)  !use xplasma spline maps
    end if
    rfbm(i)=rmjion*1.e-2_fp	! r [m]
    zfbm(i)=yion*1.e-2_fp	! z [m]
    rfbmp=rfbm(i)
    zfbmp=zfbm(i)
    bfbm(i)=bimodb(1) !tesla
    mufbm(i)=efbm(i)/bfbm(i)*(1.0_fp-pfbm(i)*pfbm(i)) ! un-normalized definition
    call eq_rgetf(1, xion, id_psi, 0, polfbm(i), ierr)
   if(ierr.ne.0) then
       call errorlog(' rbq_driver failed in getting psi at xion=',xion)
       return
    endif
    polfbm(i)=polfbm(i)/psi_lcfs
    pphfbm(i)=-polfbm(i)+xksid*vion*big(1)*1.0e-6_fp*&
         sqrt(2.0_fp*v2toev*1.e-3_fp*abeams(lsbeam)*zmp_zelonzc)/ &
         psi_lcfs/bimodb(1)/xzbeams(lsbeam)
    polfbmp=polfbm(i)
    pphfbmp=pphfbm(i)
    call get_diff_rbq(zdifpphi,ie,ipphi,ilambda,ierr) !find diffusion coefficient
    zdifpphi=max(0.0_fp,zdifpphi)
    diff_sign=dsign(1.0_fp,1.0_fp -2.0_fp*nbrandm_aep(idum))
    zdifpphi=sqrt(2.0_fp*zdt*zdifpphi)
    pphfbm(i)=pphfbm(i) + diff_sign*zdifpphi
    call orb_update(i,polfbm(i),is_ok)
    if(is_ok.eq.1.and.polfbm(i).lt.1.0_fp.and.polfbm(i).gt.0.0_fp) then ! update ptcl characteristics
       call xplasma_ctrans(s, ierr, r_in =rfbm(i) , z_in = zfbm(i), &
            maptype=1, tol=1.0e-6_fp, i2pi=2, rho_out = xion, theta_out = th)     

       if(ierr.eq.0) then
          
          bixi(1) = xion
          bithst(1) = th
          call bintrp
          bfbm(i)=bimodb(1)
          rmjion=rfbm(i)*1.0e+2_fp
          yion=zfbm(i)*1.0e+2_fp
          xksid=pfbm(i)
          thfbm(i)=th
          vperp=sqrt(mufbm(i)/efbm(i)*bfbm(i))*vion
          vpll2=vion*vion*(1.0_fp- mufbm(i)/efbm(i)*bfbm(i))
          if(vpll2.lt.0.0_fp) then
             vpll=0.0_fp
             vperp=vion
          else
             vpll=sign(sqrt(vpll2),xksid)
             xksid=vpll/vion
          endif
          ! update pphi
          call gintrp(1)          
          call nbigc(-1)  !update ngc, ngc2          
          call nbiflr(0.0_fp,0) !nflr at theta=0.0    
          bivpar(1)=vpll
          if(aepp(nk,1).gt.1) then
             if(pphfbm(i).le.aepp(nk,2).or.pphfbm(i).gt.aepp(nk,3)) then
                il_unknown=.true.
             else
                il_unknown=.false.
                ipphi=1+ (pphfbm(i)-aepp(nk,2))*(aepp(nk,1)-1)/(aepp(nk,3)-aepp(nk,2))
             endif
          endif
          pphfbm(i)=-polfbm(i)+xksid*vion*big(1)*1.0e-6_fp*&
               sqrt(2.0_fp*v2toev*1.e-3_fp*abeams(lsbeam)*zmp_zelonzc)/ &
               psi_lcfs/bimodb(1)/xzbeams(lsbeam)
          if(nltest_output.and..not.il_unknown) then
             pph_ms(ie,ipphi,ilambda,nk,1)=pph_ms(ie,ipphi,ilambda,nk,1)+wght*zdt/delt
             aep_count(ie,ipphi,ilambda,nk,1)=aep_count(ie,ipphi,ilambda,nk,1)+zdt/delt
             pbaep_ms(ie,ipphi,ilambda,nk,1)=pbaep_ms(ie,ipphi,ilambda,nk,1)+polfbm(i)*wght*zdt/delt
          endif
       endif
    endif
  end subroutine rbq_driver
  
  subroutine get_diff_rbq(zdifpphi,ie_loc,ip_loc,imu_loc,ierr)
    use physconst_mod, only: aee,mp_g
    use nbi_aep_com, only : aepe,aepp,aepmu,aepd
    use nbi_com, only : raxis
   
    real(fp), intent(out) :: zdifpphi
    integer, intent(out)  :: ierr
    integer, intent(out)  :: ie_loc,ip_loc,imu_loc
    integer :: ik,nk,i
    real(fp) :: psi_lcfs,polflx(1)
    real(fp) :: z_e, z_pphi, z_mu,zlambda,zpphi
    real(fp) :: he_1,hpphi_1,hmu_1 ! inverse  step
    ierr=0
    i=1
    zdifpphi=0.0_fp
    nk=count(aep_ascale.lt.0.0_fp)
    if (nk.eq.0) then
       ierr=1
       return
    endif
    do ik=1,nk
       ie_loc=1
       if(aepe(ik,1).gt.1) then
          z_e=max(aepe(ik,2),min(aepe(ik,3),efbm(i)))
          he_1=(aepe(ik,1)-1)/(aepe(ik,3)-aepe(ik,2))
          ie_loc=1+ (z_e-aepe(ik,2))*he_1
       endif
       
       ip_loc=1
       if(aepp(ik,1).gt.1) then
          z_pphi=max(aepp(ik,2),min(aepp(ik,3),pphfbm(i)))
          hpphi_1=(aepp(ik,1)-1)/(aepp(ik,3)-aepp(ik,2))
          ip_loc=1+ (z_pphi-aepp(ik,2))*hpphi_1
       endif

       ! definition: lambda=(1-xksid^2)*B_axis*R_axis/B = mu/E *B_axis*R_axis/B 
       zlambda=(1.0_fp-pfbm(i)*pfbm(i))*baxis*raxis*1.0e-2_fp/bfbm(i)
       
       
       imu_loc=1
       if(aepmu(ik,1).gt.1) then
          z_mu=max(aepmu(ik,2),min(aepmu(ik,3),zlambda))
          hmu_1=(aepmu(ik,1)-1)/(aepmu(ik,3)-aepmu(ik,2))
          imu_loc=1+ (z_mu-aepmu(ik,2))*hmu_1
       endif
       
       if(ie_loc.eq.int(aepe(ik,1))) then
          if(ip_loc.eq.int(aepp(ik,1))) then
             if(imu_loc.eq.int(aepmu(ik,1))) then
                zdifpphi=aepd(1,1,ie_loc,ip_loc,imu_loc,ik)
             else
                call interp_linear1(aepmu(ik,2),hmu_1,z_mu,imu_loc,&
                     aepd(1:1,1:1,ie_loc:ie_loc,ip_loc:ip_loc,imu_loc:imu_loc+1,ik:ik),&
                     zdifpphi)
             endif
          else
             if(imu_loc.eq.int(aepmu(ik,1))) then
                call interp_linear1(aepp(ik,2),hpphi_1,z_pphi,ip_loc,&
                     aepd(1:1,1:1,ie_loc:ie_loc,ip_loc:ip_loc+1,imu_loc:imu_loc,ik:ik),&
                     zdifpphi)
             else
                call interp_linear2(aepp(ik,2),hpphi_1,z_pphi,ip_loc,&
                     aepmu(ik,2),hmu_1,z_mu,imu_loc,&
                     aepd(1:1,1:1,ie_loc:ie_loc,ip_loc:ip_loc+1,imu_loc:imu_loc+1,ik:ik),&
                     zdifpphi)
             endif
          endif !ip_loc
       else
          if(ip_loc.eq.int(aepp(ik,1))) then
             if(imu_loc.eq.int(aepmu(ik,1))) then
                call interp_linear1(aepe(ik,2),he_1,z_e,ie_loc,&
                     aepd(1:1,1:1,ie_loc:ie_loc+1,ip_loc:ip_loc,imu_loc:imu_loc,ik:ik),&
                     zdifpphi)
             else
                call interp_linear2(aepe(ik,2),he_1,z_e,ie_loc,&
                     aepmu(ik,2),hmu_1,z_mu,imu_loc,&
                     aepd(1:1,1:1,ie_loc:ie_loc+1,ip_loc:ip_loc,imu_loc:imu_loc+1,ik:ik),&
                     zdifpphi)
             endif
          else
             if(imu_loc.eq.int(aepmu(ik,1))) then
                call interp_linear2(aepe(ik,2),he_1,z_e,ie_loc,&
                     aepp(ik,2),hpphi_1,z_pphi,ip_loc,&
                     aepd(1:1,1:1,ie_loc:ie_loc+1,ip_loc:ip_loc+1,imu_loc:imu_loc,ik:ik),&
                     zdifpphi)
             else
                call interp_linear3(aepe(ik,2),he_1,z_e,ie_loc,&
                     aepp(ik,2),hpphi_1,z_pphi,ip_loc,&
                     aepmu(ik,2),hmu_1,z_mu,imu_loc,&
                     aepd(1:1,1:1,ie_loc:ie_loc+1,ip_loc:ip_loc+1,imu_loc:imu_loc+1,ik:ik),&
                     zdifpphi)
             endif
          endif
       endif
    enddo
    
  end subroutine get_diff_rbq

  subroutine interp_linear1(z1,hz_1,zz,iz,aep_d,zdiff)
    real(fp), intent(in) :: z1,hz_1,zz,aep_d(2)
    integer :: iz
    real(fp), intent(out) :: zdiff

    zdiff=aep_d(1)+((zz-z1)*hz_1-iz+1)*(aep_d(2)-aep_d(1))

  end subroutine interp_linear1

  subroutine interp_linear2(z1,hz_1,zz,iz,y1,hy_1,yy,iy,aep_d,zdiff)
    real(fp), intent(in) :: z1,hz_1,zz,y1,hy_1,yy,aep_d(2,2)
    integer, intent(in)  :: iz,iy
    real(fp), intent(out) :: zdiff
    real(fp) :: y_y1,y_y2,z_z1,z_z2

    y_y1=(yy-y1)*hy_1-iy+1.0_fp
    y_y2=y_y1-1.0_fp
    z_z1=(zz-z1)*hz_1-iz+1.0_fp
    z_z2=z_z1-1.0_fp

    zdiff=aep_d(1,1)*y_y2*z_z2-aep_d(2,1)*y_y1*z_z2-&
         aep_d(1,2)*y_y2*z_z1+aep_d(2,2)*y_y1*z_z1

  end subroutine interp_linear2

  subroutine interp_linear3(z1,hz_1,zz,iz,y1,hy_1,yy,iy,x1,hx_1,xx,ix,aep_d,zdiff)
    real(fp), intent(in) :: z1,hz_1,zz,y1,hy_1,yy,x1,hx_1,xx,aep_d(2,2,2)
    integer, intent(in)  :: iz,iy,ix
    real(fp), intent(out) :: zdiff
    real(fp) :: y_y1,y_y2,z_z1,z_z2,x_x1,x_x2

    x_x1=(xx-x1)*hx_1-ix+1.0_fp
    x_x2=x_x1-1.0_fp
    y_y1=(yy-y1)*hy_1-iy+1.0_fp
    y_y2=y_y1-1.0_fp
    z_z1=(zz-z1)*hz_1-iz+1.0_fp
    z_z2=z_z1-1.0_fp

    zdiff=-aep_d(1,1,1)*x_x2*y_y2*z_z2+aep_d(2,1,1)*x_x1*y_y2*z_z2+&
         aep_d(1,2,1)*x_x2*y_y1*z_z2+aep_d(1,1,2)*x_x2*y_y2*z_z1-&
         aep_d(2,1,2)*x_x1*y_y2*z_z1-aep_d(1,2,2)*x_x2*y_y1*z_z1-&
         aep_d(2,2,1)*x_x1*y_y1*z_z2+aep_d(2,2,2)*x_x1*y_y1*z_z1

  end subroutine interp_linear3

end module nbi_aep_mod
