subroutine orbrz_rk4(iorb, istep, zdt, ierr)

  use orbrzv
  use nbi_com,only:epslon,zc

  implicit none
  integer, intent(in) :: iorb, istep ! orbit index and time step from orbit.for
  real(fp), intent(in) :: zdt
  integer, intent(out) :: ierr  ! status code, 0 = ok

  real(fp) :: zms, zqs, zsp
  real(fp) :: zzk1(1:3), zzk2(1:3), zzk3(1:3), zztmp(1:3), zzscal(1:3)
  integer :: i

  ! initialize error flag and lost flag
  ierr  = 0

  zms = orbrzv_ms
  zqs = orbrzv_qs

  ! bphi sign
  zsp = orbrzv_bpsign

  ! backup initial status
  orbrzv_zz0(1:4) = orbrzv_zz(1:4)
  orbrzv_dzzdt0(1:3) = orbrzv_dzzdt(1:3)
  orbrzv_psiv0(1:6) = orbrzv_psiv(1:6)
  orbrzv_gv0(1:3) = orbrzv_gv(1:3)
  orbrzv_phiv0(1:3) = orbrzv_phiv(1:3)
  orbrzv_b0 = orbrzv_b

  ! conserved quantities at start
  orbrzv_teng0 = 0.5_fp*zms*orbrzv_zz0(3)**2 + orbrzv_zz0(4)*orbrzv_b0 &
             & + zqs*orbrzv_phiv0(3)
  orbrzv_tcmt0 = zms*orbrzv_zz0(3)*orbrzv_gv0(1)/orbrzv_b0 &
             & - zsp*zqs*orbrzv_psiv0(1)/zc

  ! set scaling facter for later error estimation
  !zzscal(1:3) = dabs(orbrzv_zz0(1:3)) + dabs(zdt*orbrzv_dzzdt(1:3))
  zzscal(1:3) = dabs(orbrzv_zz0(1:3) + zdt*orbrzv_dzzdt(1:3))
  zzscal(1:3) = zzscal(1:3) + epslon
  !zzscal(1:3) = dmax1(zzscal(1:3), orbrzv_zzscal_min(1:3))

  ! ------------ rk1 -----------------------------------------------------
  orbrzv_zz(1:3) = orbrzv_zz0(1:3) + 0.5_fp*zdt*orbrzv_dzzdt(1:3)
  zzk1(1:3) = orbrzv_dzzdt(1:3)
  ! check whether particle is inside simulation domain
  call orbrz_check_bdy(iorb, istep)

  ! ------------ rk2 -----------------------------------------------------
  if(orbrzv_pflag >= 0) then
    call orbrz_derivs(2, orbrzv_zz, zms, zqs, orbrzv_dzzdt)
    orbrzv_zz(1:3) = orbrzv_zz0(1:3) + 0.5_fp*zdt*orbrzv_dzzdt(1:3)
    zzk2(1:3) = orbrzv_dzzdt(1:3)
    ! check whether particle is inside simulation domain
    call orbrz_check_bdy(iorb, istep)
  end if

  ! ------------ rk3 -----------------------------------------------------
  if(orbrzv_pflag >= 0) then
    call orbrz_derivs(3, orbrzv_zz, zms, zqs, orbrzv_dzzdt)
    orbrzv_zz(1:3) = orbrzv_zz0(1:3) + zdt*orbrzv_dzzdt(1:3)
    zzk3(1:3) = orbrzv_dzzdt(1:3)
    ! check whether particle is inside simulation domain
    call orbrz_check_bdy(iorb, istep)
  end if

  ! ------------ rk4 -----------------------------------------------------
  if(orbrzv_pflag >= 0) then
    call orbrz_derivs(4, orbrzv_zz, zms, zqs, orbrzv_dzzdt)
    orbrzv_zz(1:3) = orbrzv_zz0(1:3) + zdt*(zzk1(1:3) + 2.0_fp*zzk2(1:3) &
                 & + 2.0_fp*zzk3(1:3) + orbrzv_dzzdt(1:3))/6.0_fp
    ! check whether particle is inside simulation domain
    call orbrz_check_bdy(iorb, istep)
  end if

  ! ------------ after rk ------------------------------------------------
  ! check ilost after completing rk4 step.
  ! if inside, prepare next time step and estimate integration err
  ! if lost, ...
  if(orbrzv_pflag >= 0) then
    call orbrz_derivs(1, orbrzv_zz, zms, zqs, orbrzv_dzzdt)

    zztmp(1:3) = orbrzv_zz0(1:3) + zdt*(zzk1(1:3) + 2.0_fp*zzk2(1:3) &
             & + 2.0_fp*zzk3(1:3) + orbrzv_dzzdt(1:3))/6.0_fp
    orbrzv_zzerr(1:3) = (orbrzv_zz(1:3) - zztmp(1:3))/6.0_fp
    orbrzv_zzerr(1:3) = dabs(orbrzv_zzerr(1:3))/zzscal(1:3)
  else ! lost particle
    ! guess recycling position and refine it
    !orbrzv_zz(1:4) = orbrzv_zz0(1:4)
    !orbrzv_zz(2) = -orbrzv_zz0(2)
    !call orbrz_findcrz(r, z, psi0, b0, 1.0e-5, crz_cnt, 10)

    !call orbrz_derivs(1, orbrzv_zz, ms, qs, orbrzv_dzzdt)
    ! set error equal to zero
    orbrzv_zzerr(1:3) = 0.0_fp
    orbrzv_teng = orbrzv_teng0
    orbrzv_tcmt = orbrzv_tcmt0

    ! leave the orbit integrator as bad orbit by retaining orbrzv_pflag = -1
  end if

  orbrzv_zzerr(4) = dabs(orbrzv_teng-orbrzv_teng0)
  orbrzv_zzerr(5) = dabs(orbrzv_tcmt-orbrzv_tcmt0)
  orbrzv_zzerr(4) = orbrzv_zzerr(4)/(dabs(orbrzv_teng0) + epslon)
  orbrzv_zzerr(5) = orbrzv_zzerr(5)/(dabs(orbrzv_tcmt0) + epslon)

  ! set maximum error to orbrzv_zzerr(0)
  orbrzv_zzerr(0) = 0.0_fp
  do i = 1, 5 ! kjm-rfdbg
    orbrzv_zzerr(0) = dmax1(orbrzv_zzerr(0), orbrzv_zzerr(i))
  end do

end subroutine orbrz_rk4
