!> Particle mover with Newtonian/relativistic Guiding Center Approximation (GCA)
!> By Jannis Teunissen, Bart Ripperda, Oliver Porth, and Fabio Bacchini (2016-2020)
module mod_particle_gca
  use mod_particle_base

  private

  !> Variable index for gradient B, with relativistic correction 1/kappa
  !> where kappa = 1/sqrt(1 - E_perp^2/B^2)
  integer, protected, allocatable      :: grad_kappa_B(:)
  !> Variable index for (B . grad)B (curvature B drift)
  integer, protected, allocatable      :: b_dot_grad_b(:)

  ! ExB related drifts (vE = ExB/B^2)
  !> Variable index for curvature drift
  integer, protected, allocatable      :: vE_dot_grad_b(:)
  !> Variable index for polarization drift
  integer, protected, allocatable      :: b_dot_grad_vE(:)
  !> Variable index for polarization drift
  integer, protected, allocatable      :: vE_dot_grad_vE(:)

  public :: gca_init
  public :: gca_create_particles
  integer, parameter :: RK4=1, ARK4=2

  ! Variables
  public :: bp, ep, grad_kappa_B, b_dot_grad_b
  public :: vE_dot_grad_b, b_dot_grad_vE, vE_dot_grad_vE
  public :: vp, jp, rhop

contains

  subroutine gca_init()
    use mod_global_parameters
    integer :: idir, nwx

    if (physics_type/='mhd') call mpistop(&
       "GCA particles need magnetic field!")
    if (ndir/=3) call mpistop("GCA particles need ndir=3!")

    nwx = 0
    allocate(bp(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      bp(idir) = nwx
    end do
    allocate(ep(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      ep(idir) = nwx
    end do
    allocate(grad_kappa_B(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      grad_kappa_B(idir) = nwx
    end do
    allocate(b_dot_grad_b(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      b_dot_grad_b(idir) = nwx
    end do
    allocate(vE_dot_grad_b(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      vE_dot_grad_b(idir) = nwx
    end do
    allocate(b_dot_grad_vE(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      b_dot_grad_vE(idir) = nwx
    end do
    allocate(vE_dot_grad_vE(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      vE_dot_grad_vE(idir) = nwx
    end do
    allocate(vp(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      vp(idir) = nwx
    end do
    allocate(jp(ndir))
    do idir = 1, ndir
      nwx = nwx + 1
      jp(idir) = nwx
    end do
    nwx = nwx + 1 ! density
    rhop = nwx

    ngridvars=nwx

    particles_fill_gridvars => gca_fill_gridvars

    if (associated(particles_define_additional_gridvars)) then
      call particles_define_additional_gridvars(ngridvars)
    end if

    select case(integrator_type_particles)
    case('RK4','Rk4','rk4')
      integrator = RK4
    case('ARK4','ARk4','Ark4','ark4')
      integrator = ARK4
    case default
      integrator = ARK4
    end select

    particles_integrate => gca_integrate_particles

  end subroutine gca_init

  subroutine gca_create_particles()
    ! initialise the particles
    use mod_global_parameters
    use mod_usr_methods, only: usr_create_particles, usr_update_payload,&
        usr_check_particle

    double precision :: b(ndir), u(ndir), magmom
    double precision :: bnorm, lfac, vnorm, vperp, vpar
    integer          :: igrid, ipe_particle
    integer          :: n, idir, nparticles_local
    double precision :: x(3, num_particles)
    double precision :: v(3, num_particles)
    double precision :: q(num_particles)
    double precision :: m(num_particles)
    double precision :: rrd(num_particles,ndir)
    double precision :: defpayload(ndefpayload)
    double precision :: usrpayload(nusrpayload)
    logical          :: follow(num_particles), check

    if (mype==0) then
      if (.not. associated(usr_create_particles)) then
        ! Randomly distributed
        do idir=1,ndir
          do n = 1, num_particles
            rrd(n,idir) = rng%unif_01()
          end do
        end do
        do n=1, num_particles
          x(1,n) = xprobmin1 + rrd(n+1,1) * (xprobmax1 - xprobmin1)
          x(2,n) = xprobmin2 + rrd(n+1,2) * (xprobmax2 - xprobmin2)
          x(3,n) = xprobmin3 + rrd(n+1,3) * (xprobmax3 - xprobmin3)
        end do
      else
        call usr_create_particles(num_particles, x, v, q, m, follow)
      end if
    end if

    call MPI_BCAST(x,3*num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(v,3*num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(q,num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(m,num_particles,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)
    call MPI_BCAST(follow,num_particles,MPI_LOGICAL,0,icomm,ierrmpi)

    nparticles_local = 0

    ! first find ipe and igrid responsible for particle
    do n = 1, num_particles
      call find_particle_ipe(x(:, n),igrid,ipe_particle)

      particle(n)%igrid = igrid
      particle(n)%ipe   = ipe_particle

      if(ipe_particle == mype) then
        check = .true.

        ! Check for user-defined modifications or rejection conditions
        if (associated(usr_check_particle)) call usr_check_particle(igrid, x(:,&
           n), v(:,n), q(n), m(n), follow(n), check)
        if (check) then
          call push_particle_into_particles_on_mype(n)
        else
          cycle
        end if

        nparticles_local = nparticles_local + 1

        call get_lfac_from_velocity(v(:, n), lfac)

        allocate(particle(n)%self)
        particle(n)%self%x      = x(:, n)
        particle(n)%self%q      = q(n)
        particle(n)%self%m      = m(n)
        particle(n)%self%follow = follow(n)
        particle(n)%self%index  = n
        particle(n)%self%time   = global_time
        particle(n)%self%dt     = 0.0d0

        call get_vec(bp, igrid, x(:, n), particle(n)%self%time, b)

        bnorm = norm2(b(:))
        vnorm = norm2(v(:, n))
        vpar  = sum(v(:, n) * b/bnorm)
        vperp = sqrt(vnorm**2 - vpar**2)

        ! The momentum vector u(1:3) is filled with the following components

        ! parallel momentum component (gamma v||)
        particle(n)%self%u(1) = lfac * vpar

        ! Mr: the conserved magnetic moment
        magmom = m(n) * (vperp * lfac)**2 / (2.0d0 * bnorm)
        particle(n)%self%u(2) = magmom

        ! Lorentz factor
        particle(n)%self%u(3) = lfac

        ! initialise payloads for GCA module
        allocate(particle(n)%payload(npayload))
        call gca_update_payload(igrid,x(:,n),particle(n)%self%u(:),q(n),m(n),&
           defpayload,ndefpayload,0.d0)
        particle(n)%payload(1:ndefpayload) = defpayload
        if (associated(usr_update_payload)) then
          call usr_update_payload(igrid,x(:,n),particle(n)%self%u(:),q(n),m(n),&
             usrpayload,nusrpayload,0.d0)
          particle(n)%payload(ndefpayload+1:npayload) = usrpayload
        end if
      end if
    end do

    call MPI_ALLREDUCE(nparticles_local,nparticles,1,MPI_INTEGER,MPI_SUM,icomm,&
       ierrmpi)

  end subroutine gca_create_particles

  subroutine gca_fill_gridvars
    use mod_global_parameters
    use mod_usr_methods, only: usr_particle_fields
    use mod_geometry

    integer                                   :: igrid, iigrid, idir, idim
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:ndir) :: beta
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:nw)   :: w,wold
    double precision                          :: current(ixGlo1:ixGhi1,&
       ixGlo2:ixGhi2,ixGlo3:ixGhi3,7-2*ndir:3)
    integer                                   :: idirmin
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
       1:ndir) :: vE, bhat
    double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
       ixGlo3:ixGhi3)        :: kappa, kappa_B, absB, tmp

    call fill_gridvars_default()

    do iigrid=1,igridstail; igrid=igrids(iigrid);
      w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         1:nw) = ps(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,1:nw)
      call phys_to_primitive(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,ixGlo1,&
         ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,w,ps(igrid)%x)

      ! grad(kappa B)
      absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) = &
         sqrt(sum(gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(:))**2,dim=ndim+1))
      vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         1) = gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(2)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(3)) - gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(3)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(2))
      vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         2) = gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(3)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(1)) - gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(1)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(3))
      vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         3) = gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(1)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(2)) - gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         ep(2)) * gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
         bp(1))
      do idir=1,ndir
        where (absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) .gt. 0.d0)
          vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             idir) = vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             idir) / absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3)**2
        elsewhere
          vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,idir) = 0.d0
        end where
      end do
      if (any(sum(vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,:)**2,&
         dim=ndim+1) .ge. c_norm**2) .and. relativistic) then
        call mpistop("GCA FILL GRIDVARS: vE>c! ABORTING...")
      end if
      if (any(vE .ne. vE)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN vE! ABORTING...")
      end if

      if (relativistic) then
        kappa(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
           ixGlo3:ixGhi3) = 1.d0/sqrt(1.0d0 - sum(vE(ixGlo1:ixGhi1,&
           ixGlo2:ixGhi2,ixGlo3:ixGhi3,:)**2,dim=ndim+1)/c_norm**2)
      else
        kappa(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) = 1.d0
      end if
      kappa_B(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) = absB(ixGlo1:ixGhi1,&
         ixGlo2:ixGhi2,ixGlo3:ixGhi3) / kappa(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
         ixGlo3:ixGhi3)

      if (any(kappa_B .ne. kappa_B)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN kappa_B! ABORTING...")
      end if

      tmp=0.d0
      do idim=1,ndim
        call gradient(kappa_B,ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
           ixGlo1+1,ixGlo2+1,ixGlo3+1,ixGhi1-1,ixGhi2-1,ixGhi3-1,idim,tmp)
        gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
           grad_kappa_B(idim)) = tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
           ixGlo3:ixGhi3)
      end do

      ! bhat
      do idir=1,ndir
        where (absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3) .gt. 0.d0)
          bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             idir) = gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,bp(idir)) / absB(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
        elsewhere
          bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,idir) = 0.d0
        end where
      end do
      if (any(bhat .ne. bhat)) then
        call mpistop("GCA FILL GRIDVARS: NaNs IN bhat! ABORTING...")
      end if

      do idir=1,ndir
        ! (b dot grad) b and the other directional derivatives
        do idim=1,ndim
          call gradient(bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,idir),&
             ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,ixGlo1+1,ixGlo2+1,&
             ixGlo3+1,ixGhi1-1,ixGhi2-1,ixGhi3-1,idim,tmp)
          gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             b_dot_grad_b(idir)) = gridvars(igrid)%w(ixGlo1:ixGhi1,&
             ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             b_dot_grad_b(idir)) + bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,idim) * tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
          gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             vE_dot_grad_b(idir)) = gridvars(igrid)%w(ixGlo1:ixGhi1,&
             ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             vE_dot_grad_b(idir)) + vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,idim) * tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
          call gradient(vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,idir),&
             ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,ixGlo1+1,ixGlo2+1,&
             ixGlo3+1,ixGhi1-1,ixGhi2-1,ixGhi3-1,idim,tmp)
          gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             b_dot_grad_vE(idir)) = gridvars(igrid)%w(ixGlo1:ixGhi1,&
             ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             b_dot_grad_vE(idir)) + bhat(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,idim) * tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
          gridvars(igrid)%w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             vE_dot_grad_vE(idir)) = gridvars(igrid)%w(ixGlo1:ixGhi1,&
             ixGlo2:ixGhi2,ixGlo3:ixGhi3,&
             vE_dot_grad_vE(idir)) + vE(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3,idim) * tmp(ixGlo1:ixGhi1,ixGlo2:ixGhi2,&
             ixGlo3:ixGhi3)
        end do
      end do

    end do

  end subroutine gca_fill_gridvars

  subroutine gca_integrate_particles(end_time)
    use mod_odeint
    use mod_global_parameters
    use mod_usr_methods, only: usr_create_particles, usr_update_payload
    double precision, intent(in)        :: end_time

    double precision                    :: lfac, absS
    double precision                    :: defpayload(ndefpayload)
    double precision                    :: usrpayload(nusrpayload)
    double precision                    :: dt_p, tloc, y(ndir+2),dydt(ndir+2),&
       ytmp(ndir+2), euler_cfl, int_factor
    double precision                    :: tk, k1(ndir+2),k2(ndir+2),&
       k3(ndir+2),k4(ndir+2)
    double precision, dimension(1:ndir) :: x, vE, e, b, bhat, x_new, vfluid,&
        current
    double precision, dimension(1:ndir) :: drift1, drift2
    double precision, dimension(1:ndir) :: drift3, drift4, drift5, drift6,&
        drift7
    double precision, dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
    double precision, dimension(1:ndir) :: bdotgradvE, vEdotgradvE
    double precision, dimension(1:ndir) :: gradBdrift, reldrift,&
        bdotgradbdrift
    double precision, dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
    double precision, dimension(1:ndir) :: vEdotgradvEdrift
    double precision                    :: kappa, Mr, upar, m, absb, gamma, q,&
        mompar, vpar, vEabs
    double precision                    :: gradBdrift_abs, reldrift_abs, epar
    double precision                    :: bdotgradbdrift_abs,&
        vEdotgradbdrift_abs
    double precision                    :: bdotgradvEdrift_abs,&
        vEdotgradvEdrift_abs
    double precision                    :: momentumpar1, momentumpar2,&
        momentumpar3, momentumpar4
    ! Precision of time-integration:
    double precision,parameter          :: eps=1.0d-6
    ! for odeint:
    double precision                    :: h1, hmin, h_old
    integer                             :: nok, nbad, ic11,ic12,ic13, ic21,&
       ic22,ic23, ierror, nvar
    integer                             :: ipart, iipart, seed, ic1,ic2,ic3,&
       igrid_particle, ipe_particle, ipe_working
    logical                             :: int_choice
    logical                             :: BC_applied

    nvar=ndir+2

    do iipart=1,nparticles_active_on_mype
      ipart                   = particles_active_on_mype(iipart)
      int_choice              = .false.
      dt_p                    = gca_get_particle_dt(particle(ipart), end_time)
      particle(ipart)%self%dt = dt_p

      igrid_working           = particle(ipart)%igrid
      ipart_working           = particle(ipart)%self%index
      tloc                    = particle(ipart)%self%time
      x(1:ndir)               = particle(ipart)%self%x(1:ndir)

      select case (integrator)
      case (RK4)
        ! Runge-Kutta order 4
        tk = tloc
        y(1:ndir) = x(1:ndir) ! position of guiding center
        y(ndir+1) = particle(ipart)%self%u(1) !parallel momentum component (gamma v||)
        y(ndir+2) = particle(ipart)%self%u(2) ! conserved magnetic moment Mr
        call derivs_gca_rk(tk,y,k1)
        tk = tloc + dt_p/2.d0
        y(1:ndir) = x + dt_p*k1(1:ndir)/2.d0
        y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k1(ndir+1)/2.d0
        call derivs_gca_rk(tk,y,k2)
        tk = tloc + dt_p/2.d0
        y(1:ndir) = x + dt_p*k2(1:ndir)/2.d0
        y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k2(ndir+1)/2.d0
        call derivs_gca_rk(tk,y,k3)
        tk = tloc + dt_p
        y(1:ndir) = x + dt_p*k3(1:ndir)
        y(ndir+1) = particle(ipart)%self%u(1) + dt_p*k3(ndir+1)
        call derivs_gca_rk(tk,y,k4)
        y(1:ndir) = x(1:ndir) ! position of guiding center
        y(ndir+1) = particle(ipart)%self%u(1) !parallel momentum component (gamma v||)
        y = y + dt_p/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
        particle(ipart)%self%x(1:ndir) = y(1:ndir)
        particle(ipart)%self%u(1)      = y(ndir+1)

      case (ARK4)
        ! Adaptive stepwidth RK4:
        ! initial solution vector:
        y(1:ndir) = x(1:ndir) ! position of guiding center
        y(ndir+1) = particle(ipart)%self%u(1) !parallel momentum component (gamma v||)
        y(ndir+2) = particle(ipart)%self%u(2) ! conserved magnetic moment Mr
      ! y(ndir+3) = particle(ipart)%self%u(3) ! Lorentz factor of particle
  
        ! we temporarily save the solution vector, to replace the one from the euler
        ! timestep after euler integration
        ytmp=y
  
        call derivs_gca(particle(ipart)%self%time,y,dydt)
  
        ! make an Euler step with the proposed timestep:
        ! factor to ensure we capture all particles near the internal ghost cells.
        ! Can be adjusted during a run, after an interpolation error.
        euler_cfl=2.5d0
  
        ! new solution vector:
        y(1:ndir+2) = y(1:ndir+2) + euler_cfl * dt_p * dydt(1:ndir+2)
        particle(ipart)%self%x(1:ndir) = y(1:ndir) !position of guiding center
        particle(ipart)%self%u(1)      = y(ndir+1) !parallel momentum component(gamma v||)
        particle(ipart)%self%u(2)      = y(ndir+2) ! conserved magnetic moment
  
        ! check if the particle is in the internal ghost cells
        int_factor =1.0d0
  
        if(.not. particle_in_igrid(ipart_working,igrid_working)) then
          ! if particle is not in the grid an euler timestep is taken instead of a RK4
          ! timestep. Then based on that we do an interpolation and check how much further
          ! the timestep for the RK4 has to be restricted.
          ! factor to make integration more accurate for particles near the internal
          ! ghost cells. This factor can be changed during integration after an
          ! interpolation error. But one should be careful with timesteps for i/o
  
          ! flat interpolation:
          ic1 = int((y(1)-rnode(rpxmin1_,igrid_working))/rnode(rpdx1_,&
             igrid_working)) + 1 + nghostcells
          ic2 = int((y(2)-rnode(rpxmin2_,igrid_working))/rnode(rpdx2_,&
             igrid_working)) + 1 + nghostcells
          ic3 = int((y(3)-rnode(rpxmin3_,igrid_working))/rnode(rpdx3_,&
             igrid_working)) + 1 + nghostcells
  
          ! linear interpolation:
          
          if (ps(igrid_working)%x(ic1,ic2,ic3,1) .lt. y(1)) then
             ic11 = ic1
          else
             ic11 = ic1 -1
          end if
          ic21 = ic11 + 1
          
          
          if (ps(igrid_working)%x(ic1,ic2,ic3,2) .lt. y(2)) then
             ic12 = ic2
          else
             ic12 = ic2 -1
          end if
          ic22 = ic12 + 1
          
          
          if (ps(igrid_working)%x(ic1,ic2,ic3,3) .lt. y(3)) then
             ic13 = ic3
          else
             ic13 = ic3 -1
          end if
          ic23 = ic13 + 1
          
  
          int_factor =0.5d0
  
          
          if (ic11 .le. ixGlo1-2 .or. ic21 .ge. ixGhi1+2) then
            int_factor = 0.05d0
          end if
          
          
          if (ic12 .le. ixGlo2-2 .or. ic22 .ge. ixGhi2+2) then
            int_factor = 0.05d0
          end if
          
          
          if (ic13 .le. ixGlo3-2 .or. ic23 .ge. ixGhi3+2) then
            int_factor = 0.05d0
          end if
          
  
          
          if (ic11 .eq. ixGlo1-1 .or. ic21 .eq. ixGhi1+1) then
            int_factor = 0.1d0
          end if
          
          
          if (ic12 .eq. ixGlo2-1 .or. ic22 .eq. ixGhi2+1) then
            int_factor = 0.1d0
          end if
          
          
          if (ic13 .eq. ixGlo3-1 .or. ic23 .eq. ixGhi3+1) then
            int_factor = 0.1d0
          end if
          
  
          dt_p=int_factor*dt_p
        end if
  
        ! replace the solution vector with the original as it was before the Euler timestep
        y(1:ndir+2) = ytmp(1:ndir+2)
  
        particle(ipart)%self%x(1:ndir) = ytmp(1:ndir) !position of guiding center
        particle(ipart)%self%u(1)      = ytmp(ndir+1) !parallel momentum component (gamma v||)
        particle(ipart)%self%u(2)      = ytmp(ndir+2) !conserved magnetic moment
  
        ! specify a minimum step hmin. If the timestep reaches this minimum, multiply by
        ! a factor 100 to make sure the RK integration doesn't crash
        h1 = dt_p/2.0d0; hmin=1.0d-9; h_old=dt_p/2.0d0
  
        if(h1 .lt. hmin)then
          h1=hmin
          dt_p=2.0d0*h1
        endif
  
        if (any(y .ne. y)) then
          call mpistop&
             ("NaNs DETECTED IN GCA_INTEGRATE BEFORE ODEINT CALL! ABORTING...")

        end if
  
        ! RK4 integration with adaptive stepwidth
        call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca_rk,&
           rkqs,ierror)
  
        if (ierror /= 0) then
           print *, "odeint returned error code", ierror
           print *, "1 means hmin too small, 2 means MAXSTP exceeded"
           print *, "Having a problem with particle", iipart
        end if
  
        if (any(y .ne. y)) then
          call mpistop&
             ("NaNs DETECTED IN GCA_INTEGRATE AFTER ODEINT CALL! ABORTING...")
        end if
  
        ! original RK integration without interpolation in ghost cells
  !       call odeint(y,nvar,tloc,tloc+dt_p,eps,h1,hmin,nok,nbad,derivs_gca,rkqs)
  
        ! final solution vector after rk integration
        particle(ipart)%self%x(1:ndir) = y(1:ndir)
        particle(ipart)%self%u(1)      = y(ndir+1)
        particle(ipart)%self%u(2)      = y(ndir+2)
        !particle(ipart)%self%u(3)      = y(ndir+3)
      end select

      ! now calculate other quantities, mean Lorentz factor, drifts, perpendicular velocity:
      call get_bfield(igrid_working, y(1:ndir), tloc+dt_p, b)
      call get_efield(igrid_working, y(1:ndir), tloc+dt_p, e)
!      call get_vec(bp, igrid_working,y(1:ndir),tloc+dt_p,b)
!      call get_vec(vp, igrid_working,y(1:ndir),tloc+dt_p,vfluid)
!      call get_vec(jp, igrid_working,y(1:ndir),tloc+dt_p,current)
!      e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!      e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!      e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)

      absb         = sqrt(sum(b(:)**2))
      bhat(1:ndir) = b(1:ndir) / absb

      epar         = sum(e(:)*bhat(:))

      call cross(e,bhat,vE)

      vE(1:ndir) = vE(1:ndir) / absb
      vEabs = sqrt(sum(vE(:)**2))
      if (relativistic) then
        kappa = 1.d0/sqrt(1.0d0 - sum(vE(:)**2)/c_norm**2)
      else
        kappa = 1.d0
      end if
      Mr = y(ndir+2); upar = y(ndir+1); m=particle(ipart)%self%m
      q=particle(ipart)%self%q
      if (relativistic) then
        gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*Mr*absb/m/c_norm**2)*kappa
      else
        gamma = 1.d0
      end if

      particle(ipart)%self%u(3) = gamma

      ! Time update
      particle(ipart)%self%time = particle(ipart)%self%time + dt_p

      ! Update payload
      call gca_update_payload(igrid_working,particle(ipart)%self%x,&
         particle(ipart)%self%u,q,m,defpayload,ndefpayload,&
         particle(ipart)%self%time)
      particle(ipart)%payload(1:ndefpayload) = defpayload
      if (associated(usr_update_payload)) then
        call usr_update_payload(igrid_working,particle(ipart)%self%x,&
           particle(ipart)%self%u,q,m,usrpayload,nusrpayload,&
           particle(ipart)%self%time)
        particle(ipart)%payload(ndefpayload+1:npayload) = usrpayload
      end if

    end do

  end subroutine gca_integrate_particles

  subroutine derivs_gca_rk(t_s,y,dydt)
    use mod_global_parameters

    double precision                :: t_s, y(ndir+2)
    double precision                :: dydt(ndir+2)

    double precision,dimension(ndir):: vE, b, e, x, bhat, bdotgradb,&
        vEdotgradb, gradkappaB, vfluid, current
    double precision,dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1,&
        utmp2, utmp3
    double precision                :: upar, Mr, gamma, absb, q, m, epar,&
        kappa
    integer                         :: ic1,ic2,ic3

    ! Here the terms in the guiding centre equations of motion are interpolated for
    ! the RK integration. The interpolation is also done in the ghost cells such
    ! that the RK integration does not give an error

    q = particle(ipart_working)%self%q
    m = particle(ipart_working)%self%m

    x(1:ndir) = y(1:ndir)
    upar      = y(ndir+1) ! gamma v||
    Mr        = y(ndir+2)
    !gamma     = y(ndir+3)

    if (any(x .ne. x)) then
      write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN X OR Y!"
      write(*,*) "x",x
      write(*,*) "y",y(ndir+1:ndir+2)
      call mpistop("ABORTING...")
    end if

    call get_bfield(igrid_working, x, t_s, b)
    call get_efield(igrid_working, x, t_s, e)
!    call get_vec(bp, igrid_working,x,t_s,b)
!    call get_vec(vp, igrid_working,x,t_s,vfluid)
!    call get_vec(jp, igrid_working,x,t_s,current)
!    e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!    e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!    e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
    call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
    call get_vec(vE_dot_grad_b, igrid_working,x,t_s,vEdotgradb)
    call get_vec(grad_kappa_B, igrid_working,x,t_s,gradkappaB)
    call get_vec(b_dot_grad_vE, igrid_working,x,t_s,bdotgradvE)
    call get_vec(vE_dot_grad_vE, igrid_working,x,t_s,vEdotgradvE)

    if (any(b .ne. b) .or. any(e .ne. e) .or. any(bdotgradb .ne. bdotgradb) &
       .or. any(vEdotgradb .ne. vEdotgradb) .or. any(gradkappaB .ne. &
       gradkappaB) .or. any(bdotgradvE .ne. bdotgradvE) .or. any(vEdotgradvE &
       .ne. vEdotgradvE)) then
      write(*,*) "ERROR IN DERIVS_GCA_RK: NaNs IN FIELD QUANTITIES!"
      write(*,*) "b",b
      write(*,*) "e",e
      write(*,*) "bdotgradb",bdotgradb
      write(*,*) "vEdotgradb",vEdotgradb
      write(*,*) "gradkappaB",gradkappaB
      write(*,*) "bdotgradvE",bdotgradvE
      write(*,*) "vEdotgradvE",vEdotgradvE
      call mpistop("ABORTING...")
    end if

    absb  = sqrt(sum(b(:)**2))
    if (absb .gt. 0.d0) then
      bhat(1:ndir) = b(1:ndir) / absb
    else
      bhat = 0.d0
    end if
    epar = sum(e(:)*bhat(:))

    call cross(e,bhat,vE)
    if (absb .gt. 0.d0) then
      vE(1:ndir) = vE(1:ndir) / absb
    else
      vE(1:ndir) = 0.d0
    end if

    if (relativistic) then
      kappa = 1.d0/sqrt(1.0d0 - sum(vE(:)**2)/c_norm**2)
      gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*Mr*absb/m/c_norm**2)*kappa
    else
      kappa = 1.d0
      gamma = 1.d0
    end if

    if (absb .gt. 0.d0) then
      utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
    else
      utmp1 = 0.d0
    end if
    utmp2(1:ndir) = Mr/(gamma*q)*gradkappaB(1:ndir) + m/q* &
       (upar**2/gamma*bdotgradb(1:ndir) + upar*vEdotgradb(1:ndir) + &
       upar*bdotgradvE(1:ndir) + gamma*vEdotgradvE(1:ndir))
    if (relativistic) then
      utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*vE(1:ndir)
    end if

    call cross(utmp1,utmp2,utmp3)
    u(1:ndir) = vE(1:ndir) + utmp3(1:ndir)

    ! done assembling the terms, now write rhs:
    dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
    dydt(ndir+1) = q/m*epar - Mr/(m*gamma) * sum(bhat(:)*gradkappaB(:)) + &
       sum(vE(:)*(upar*bdotgradb(:)+gamma*vEdotgradb(:)))
    dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved

  end subroutine derivs_gca_rk

  subroutine derivs_gca(t_s,y,dydt)
    use mod_global_parameters

    double precision                :: t_s, y(ndir+2)
    double precision                :: dydt(ndir+2)

    double precision,dimension(ndir):: vE, b, e, x, bhat, bdotgradb,&
        vEdotgradb, gradkappaB, vfluid, current
    double precision,dimension(ndir):: bdotgradvE, vEdotgradvE, u, utmp1,&
        utmp2, utmp3
    double precision                :: upar, Mr, gamma, absb, q, m, epar,&
        kappa

    ! Here the normal interpolation is done for the terms in the GCA equations of motion

    q = particle(ipart_working)%self%q
    m = particle(ipart_working)%self%m

    x(1:ndir) = y(1:ndir)
    upar      = y(ndir+1) ! gamma v||
    Mr        = y(ndir+2)
    !gamma     = y(ndir+3)

    if (any(x .ne. x)) then
      call mpistop("ERROR IN DERIVS_GCA: NaNs IN X! ABORTING...")
    end if

    call get_bfield(igrid_working, x, t_s, b)
    call get_efield(igrid_working, x, t_s, e)
!    call get_vec(bp, igrid_working,x,t_s,b)
!    call get_vec(vp, igrid_working,x,t_s,vfluid)
!    call get_vec(jp, igrid_working,x,t_s,current)
!    e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!    e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!    e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)
    call get_vec(b_dot_grad_b, igrid_working,x,t_s,bdotgradb)
    call get_vec(vE_dot_grad_b, igrid_working,x,t_s,vEdotgradb)
    call get_vec(grad_kappa_B, igrid_working,x,t_s,gradkappaB)
    call get_vec(b_dot_grad_vE, igrid_working,x,t_s,bdotgradvE)
    call get_vec(vE_dot_grad_vE, igrid_working,x,t_s,vEdotgradvE)

    absb         = sqrt(sum(b(:)**2))
    if (absb .gt. 0.d0) then
      bhat(1:ndir) = b(1:ndir) / absb
    else
      bhat = 0.d0
    end if

    epar         = sum(e(:)*bhat(:))
    call cross(e,bhat,vE)
    if (absb .gt. 0.d0) vE(1:ndir)   = vE(1:ndir) / absb

    if (relativistic) then
      kappa = 1.d0/sqrt(1.0d0 - sum(vE(:)**2)/c_norm**2)
      gamma = sqrt(1.0d0+upar**2/c_norm**2+2.0d0*Mr*absb/m/c_norm**2)*kappa
    else
      kappa = 1.d0
      gamma = 1.d0
    end if
    if (absb .gt. 0.d0) then
      utmp1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
    else
      utmp1 = 0.d0
    end if
    utmp2(1:ndir) = Mr/(gamma*q)*gradkappaB(1:ndir) + m/q* &
       (upar**2/gamma*bdotgradb(1:ndir) + upar*vEdotgradb(1:ndir) + &
       upar*bdotgradvE(1:ndir) + gamma*vEdotgradvE(1:ndir))
    if (relativistic) then
      utmp2(1:ndir) = utmp2(1:ndir) + upar*epar/(gamma)*vE(1:ndir)
    end if

    call cross(utmp1,utmp2,utmp3)
    u(1:ndir) = vE(1:ndir) + utmp3(1:ndir)

    ! done assembling the terms, now write rhs:
    dydt(1:ndir) = ( u(1:ndir) + upar/gamma * bhat(1:ndir) )
    dydt(ndir+1) = q/m*epar - Mr/(m*gamma) * sum(bhat(:)*gradkappaB(:)) + &
       sum(vE(:)*(upar*bdotgradb(:)+gamma*vEdotgradb(:)))
    dydt(ndir+2) = 0.0d0 ! magnetic moment is conserved

  end subroutine derivs_gca

  !> Update payload subroutine
  subroutine gca_update_payload(igrid,xpart,upart,qpart,mpart,mypayload,&
     mynpayload,particle_time)
    use mod_global_parameters
    integer, intent(in)           :: igrid,mynpayload
    double precision, intent(in)  :: xpart(1:ndir),upart(1:ndir),qpart,mpart,&
       particle_time
    double precision, intent(out) :: mypayload(mynpayload)
    double precision, dimension(1:ndir) :: vE, e, b, bhat, vfluid, current
    double precision, dimension(1:ndir) :: drift1, drift2
    double precision, dimension(1:ndir) :: drift3, drift4, drift5, drift6,&
        drift7
    double precision, dimension(1:ndir) :: bdotgradb, vEdotgradb, gradkappaB
    double precision, dimension(1:ndir) :: bdotgradvE, vEdotgradvE
    double precision, dimension(1:ndir) :: gradBdrift, reldrift,&
        bdotgradbdrift
    double precision, dimension(1:ndir) :: vEdotgradbdrift, bdotgradvEdrift
    double precision, dimension(1:ndir) :: vEdotgradvEdrift
    double precision                    :: kappa, upar, absb, gamma, vpar,&
        vEabs
    double precision                    :: gradBdrift_abs, reldrift_abs, epar
    double precision                    :: bdotgradbdrift_abs,&
        vEdotgradbdrift_abs
    double precision                    :: bdotgradvEdrift_abs,&
        vEdotgradvEdrift_abs
    double precision                    :: momentumpar1, momentumpar2,&
        momentumpar3, momentumpar4

    call get_bfield(igrid, xpart(1:ndir), particle_time, b)
    call get_efield(igrid, xpart(1:ndir), particle_time, e)
!    call get_vec(bp, igrid,xpart(1:ndir),particle_time,b)
!    call get_vec(vp, igrid,xpart(1:ndir),particle_time,vfluid)
!    call get_vec(jp, igrid,xpart(1:ndir),particle_time,current)
!    e(1) = -vfluid(2)*b(3)+vfluid(3)*b(2) + particles_eta*current(1)
!    e(2) = vfluid(1)*b(3)-vfluid(3)*b(1) + particles_eta*current(2)
!    e(3) = -vfluid(1)*b(2)+vfluid(2)*b(1) + particles_eta*current(3)

    absb         = sqrt(sum(b(:)**2))
    bhat(1:ndir) = b(1:ndir) / absb
    epar         = sum(e(:)*bhat(:))
    call cross(e,bhat,vE)
    vE(1:ndir)   = vE(1:ndir) / absb
    vEabs = sqrt(sum(vE(:)**2))
    if (relativistic) then
      kappa = 1.d0/sqrt(1.0d0 - sum(vE(:)**2)/c_norm**2)
    else
      kappa = 1.d0
    end if
    vpar = upart(1)/upart(3)
    upar = upart(1)

    call get_vec(b_dot_grad_b, igrid,xpart(1:ndir),particle_time,bdotgradb)
    call get_vec(vE_dot_grad_b, igrid,xpart(1:ndir),particle_time,vEdotgradb)
    call get_vec(grad_kappa_B, igrid,xpart(1:ndir),particle_time,gradkappaB)
    call get_vec(b_dot_grad_vE, igrid,xpart(1:ndir),particle_time,bdotgradvE)
    call get_vec(vE_dot_grad_vE, igrid,xpart(1:ndir),particle_time,&
       vEdotgradvE)

    drift1(1:ndir) = bhat(1:ndir)/(absb/kappa**2)
    drift2(1:ndir) = upart(2)/(upart(3)*qpart)*gradkappaB(1:ndir)

    call cross(drift1,drift2,gradBdrift)
    gradBdrift_abs = sqrt(sum(gradBdrift(:)**2))

    drift3(1:ndir) = upar*epar/upart(3)*vE(1:ndir)
    call cross(drift1,drift3,reldrift)
    reldrift_abs = sqrt(sum(reldrift(:)**2))

    drift4(1:ndir) = mpart/qpart* ( upar**2/upart(3)*bdotgradb(1:ndir))
    call cross(drift1,drift4,bdotgradbdrift)
    bdotgradbdrift_abs = sqrt(sum(bdotgradbdrift(:)**2))

    drift5(1:ndir) = mpart/qpart* ( upar*vEdotgradb(1:ndir))
    call cross(drift1,drift5,vEdotgradbdrift)
    vEdotgradbdrift_abs = sqrt(sum(vEdotgradbdrift(:)**2))

    drift6(1:ndir) = mpart/qpart* ( upar*bdotgradvE(1:ndir))
    call cross(drift1,drift6,bdotgradvEdrift)
    bdotgradvEdrift_abs = sqrt(sum(bdotgradvEdrift(:)**2))

    drift7(1:ndir) = mpart/qpart* (upart(3)*vEdotgradvE(1:ndir))
    call cross(drift1,drift7,vEdotgradvEdrift)
    vEdotgradvEdrift_abs = sqrt(sum(vEdotgradvEdrift(:)**2))

    momentumpar1 = qpart/mpart*epar
    momentumpar2 = -(upart(2)/mpart/upart(3))*sum(bhat(:)*gradkappaB(:))
    momentumpar3 = upar*sum(vE(:)*bdotgradb(:))
    momentumpar4 = upart(3)*sum(vE(:)*vEdotgradb(:))

    ! Payload update
    if (mynpayload > 0) then
      ! current gyroradius
      mypayload(1) = sqrt(2.0d0*mpart*upart(2)*absb)/abs(qpart*absb)
    end if
    if (mynpayload > 1) then
      ! pitch angle
      mypayload(2) = atan2(sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2)),&
         vpar)
    end if
    if (mynpayload > 2) then
      ! particle v_perp
      mypayload(3) = sqrt((2.0d0*upart(2)*absb)/(mpart*upart(3)**2))
    end if
    if (mynpayload > 3) then
      ! particle parallel momentum term 1
      mypayload(4) = momentumpar1
    end if
    if (mynpayload > 4) then
      ! particle parallel momentum term 2
      mypayload(5) = momentumpar2
    end if
    if (mynpayload > 5) then
      ! particle parallel momentum term 3
      mypayload(6) = momentumpar3
    end if
    if (mynpayload > 6) then
      ! particle parallel momentum term 4
      mypayload(7) = momentumpar4
    end if
    if (mynpayload > 7) then
      ! particle ExB drift
      mypayload(8) = vEabs
    end if
    if (mynpayload > 8) then
      ! relativistic drift
      mypayload(9) = gradBdrift_abs
    end if
    if (mynpayload > 9) then
      ! gradB drift
      mypayload(10) = reldrift_abs
    end if
    if (mynpayload > 10) then
      ! bdotgradb drift
      mypayload(11) = bdotgradbdrift_abs
    end if
    if (mynpayload > 11) then
      ! vEdotgradb drift
      mypayload(12) = vEdotgradbdrift_abs
    end if
    if (mynpayload > 12) then
      ! bdotgradvE drift
      mypayload(13) = bdotgradvEdrift_abs
    end if
    if (mynpayload > 13) then
      ! vEdotgradvE drift
      mypayload(14) = vEdotgradvEdrift_abs
    end if

  end subroutine gca_update_payload

  function gca_get_particle_dt(partp, end_time) result(dt_p)
    use mod_odeint
    use mod_global_parameters
    type(particle_ptr), intent(inout) :: partp
    double precision, intent(in)      :: end_time
    double precision                  :: dt_p

    double precision            :: tout, dt_particles_mype, dt_cfl0, dt_cfl1,&
        dt_a
    double precision            :: dxmin, vp, a, gammap
    double precision            :: v(ndir), y(ndir+2),ytmp(ndir+2),&
        dydt(ndir+2), v0(ndir), v1(ndir), dydt1(ndir+2)
    double precision            :: ap0, ap1, dt_cfl_ap0, dt_cfl_ap1, dt_cfl_ap
    double precision            :: dt_euler, dt_tmp
    ! make these particle cfl conditions more restrictive if you are interpolating out of the grid
    double precision            :: cfl, uparcfl
    double precision, parameter :: uparmin=1.0d-6*const_c
    integer                     :: ipart, iipart, nout, ic1,ic2,ic3,&
        igrid_particle, ipe_particle, ipe
    logical                     :: BC_applied

    if (const_dt_particles > 0) then
      dt_p = const_dt_particles
      return
    end if

    cfl = particles_cfl
    uparcfl = particles_cfl

    igrid_working = partp%igrid
    ipart_working = partp%self%index
    dt_tmp = (end_time - partp%self%time)
    if(dt_tmp .le. 0.0d0) dt_tmp = smalldouble
    ! make sure we step only one cell at a time, first check CFL at current location
    ! then we make an Euler step to the new location and check the new CFL
    ! we simply take the minimum of the two timesteps.
    ! added safety factor cfl:
    dxmin  = min(rnode(rpdx1_,igrid_working),rnode(rpdx2_,igrid_working),&
       rnode(rpdx3_,igrid_working),bigdouble)*cfl
    ! initial solution vector:
    y(1:ndir) = partp%self%x(1:ndir) ! position of guiding center
    y(ndir+1) = partp%self%u(1) ! parallel momentum component (gamma v||)
    y(ndir+2) = partp%self%u(2) ! conserved magnetic moment
    ytmp=y
    !y(ndir+3) = partp%self%u(3) ! Lorentz factor of guiding centre

    call derivs_gca(partp%self%time,y,dydt)
    v0(1:ndir) = dydt(1:ndir)
    ap0        = dydt(ndir+1)

    ! guiding center velocity:
    v(1:ndir) = abs(dydt(1:ndir))
    vp = sqrt(sum(v(:)**2))

    dt_cfl0    = dxmin / max(vp, smalldouble)
    dt_cfl_ap0 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap0),&
        smalldouble))
    !dt_cfl_ap0 = min(dt_cfl_ap0, uparcfl * sqrt(abs(unit_length*dxmin/(ap0+smalldouble))) )

    ! make an Euler step with the proposed timestep:
    ! new solution vector:
    dt_euler = min(dt_tmp,dt_cfl0,dt_cfl_ap0)
    y(1:ndir+2) = y(1:ndir+2) + dt_euler * dydt(1:ndir+2)

    partp%self%x(1:ndir) = y(1:ndir) ! position of guiding center
    partp%self%u(1)      = y(ndir+1) ! parallel momentum component (gamma v||)
    partp%self%u(2)      = y(ndir+2) ! conserved magnetic moment

    ! first check if the particle is outside the physical domain or in the ghost cells
    if(.not. particle_in_igrid(ipart_working,igrid_working)) then
      y(1:ndir+2) = ytmp(1:ndir+2)
    end if

    call derivs_gca_rk(partp%self%time+dt_euler,y,dydt)
    !call derivs_gca(partp%self%time+dt_euler,y,dydt)

    v1(1:ndir) = dydt(1:ndir)
    ap1        = dydt(ndir+1)

    ! guiding center velocity:
    v(1:ndir) = abs(dydt(1:ndir))
    vp = sqrt(sum(v(:)**2))

    dt_cfl1    = dxmin / max(vp, smalldouble)
    dt_cfl_ap1 = uparcfl * abs(max(abs(y(ndir+1)),uparmin) / max(abs(ap1),&
        smalldouble))
    !dt_cfl_ap1 = min(dt_cfl_ap1, uparcfl * sqrt(abs(unit_length*dxmin/(ap1+smalldouble))) )

    dt_tmp = min(dt_euler, dt_cfl1, dt_cfl_ap1)

    partp%self%x(1:ndir) = ytmp(1:ndir) ! position of guiding center
    partp%self%u(1)      = ytmp(ndir+1) !parallel momentum component (gamma v||)
    partp%self%u(2)      = ytmp(ndir+2) ! conserved magnetic moment
    !dt_tmp = min(dt_cfl1, dt_cfl_ap1)

    ! time step due to parallel acceleration:
    ! The standard thing, dt=sqrt(dx/a) where we compute a from d(gamma v||)/dt and d(gamma)/dt
    ! dt_ap = sqrt(abs(dxmin*unit_length*y(ndir+3)/( dydt(ndir+1) - y(ndir+1)/y(ndir+3)*dydt(ndir+3) ) ) )
    ! vp = sqrt(sum(v(1:ndir)**))
    ! gammap = sqrt(1.0d0/(1.0d0-(vp/const_c)**2))
    ! ap = const_c**2/vp*gammap**(-3)*dydt(ndir+3)
    ! dt_ap = sqrt(dxmin*unit_length/ap)

    !dt_a = bigdouble
    !if (dt_euler .gt. smalldouble) then
    !   a = sqrt(sum((v1(1:ndir)-v0(1:ndir))**2))/dt_euler
    !   dt_a = min(sqrt(dxmin/a),bigdouble)
    !end if

    !dt_p = min(dt_tmp , dt_a)
    dt_p = dt_tmp

    ! Make sure we do not advance beyond end_time
    call limit_dt_endtime(end_time - partp%self%time, dt_p)

  end function gca_get_particle_dt

end module mod_particle_gca