! calculate the derivatives of (r,z) variables using 
! spline interpolation
subroutine orbrz_derivs(rk_step, zz, ms, qs, dzzdt)

  use ezspline, only : ezspline_interp2,ezspline_gradient2,&
       ezspline_derivative2
  use orbrzv
  use physconst_mod, only: vc,zmtocm
  implicit none
  integer, intent(in) :: rk_step
  real(fp), intent(in) :: zz(1:4), ms, qs
  real(fp), intent(out), dimension(1:3) :: dzzdt
  real(fp), parameter :: zc = vc*zmtocm
  real(fp) :: zsp
  real(fp) :: r, z, vl, mu, r_inv, r2_inv
  real(fp) :: psi, dpsi(1:2), g, dg(1:2), phi, dphi(1:2)
  real(fp) :: dpdr, dpdz, dpdrr, dpdzz, dpdrz, ddp, dpdp, dbdr, dbdz
  real(fp) :: b, b_inv, rb_inv, rb2_inv, r2b_inv, r2b2_inv, bs, zcf
  real(fp) :: bv(1:3), gv(1:3), df(1:3), bsv(1:3), g12, g13, g23
  integer :: ier

  ! bphi sign
  zsp = orbrzv_bpsign

  r  = zz(1)
  z  = zz(2)
  vl = zz(3)
  mu = zz(4)
  r_inv = 1.0_fp/r
  r2_inv = r_inv*r_inv
  ! interpolate geometric quantities
  call ezspline_interp2(orbrzv_psi_spo, r, z, psi, ier)
  call ezspline_gradient2(orbrzv_psi_spo, r, z, dpsi, ier)
  call ezspline_derivative2(orbrzv_psi_spo, 2, 0, r, z, dpdrr, ier)
  call ezspline_derivative2(orbrzv_psi_spo, 0, 2, r, z, dpdzz, ier)
  call ezspline_derivative2(orbrzv_psi_spo, 1, 1, r, z, dpdrz, ier)
  call ezspline_interp2(orbrzv_g_spo, r, z, g, ier)
  call ezspline_gradient2(orbrzv_g_spo, r, z, dg, ier)
  call ezspline_interp2(orbrzv_phi_spo, r, z, phi, ier)
  call ezspline_gradient2(orbrzv_phi_spo, r, z, dphi, ier)

  ! store them for later usage
  orbrzv_psiv(1) = psi
  orbrzv_psiv(2) = dpsi(1)
  orbrzv_psiv(3) = dpsi(2)
  orbrzv_psiv(4) = dpdrr
  orbrzv_psiv(5) = dpdzz
  orbrzv_psiv(6) = dpdrz
  orbrzv_gv(1) = g
  orbrzv_gv(2) = dg(1)
  orbrzv_gv(3) = dg(2)
  orbrzv_phiv(1) = dphi(1)
  orbrzv_phiv(2) = dphi(2)
  orbrzv_phiv(3) = phi

  ! b
  bv(1) =  dpsi(2)*r_inv*zsp
  bv(2) = -dpsi(1)*r_inv*zsp
  bv(3) =  g*r_inv*r_inv

  ! mod b
  b = sqrt(dpsi(1)**2 + dpsi(2)**2 + g**2)*r_inv
  b_inv = 1.0_fp/b

  rb_inv   = r_inv*b_inv
  rb2_inv  = rb_inv*b_inv
  r2b_inv  = rb_inv*r_inv
  r2b2_inv = r2b_inv*b_inv

  dpdr = dpsi(1)
  dpdz = dpsi(2)
  ddp  = dpdrr + dpdzz
  dpdp = dpsi(1)**2 + dpsi(2)**2

  ! dbdr, dbdz
  dbdr = r2b_inv*(dpdr*dpdrr+dpdz*dpdrz+g*dg(1)) - b*r_inv
  dbdz = r2b_inv*(dpdz*dpdzz+dpdr*dpdrz+g*dg(2))

  ! store b & dbdr for later usage
  orbrzv_b = b
  orbrzv_dbdr = dbdr

  ! (c*m*vl/q) curl b
  zcf = zc*ms*vl/qs
  gv(1) = -dg(2) + g*b_inv*dbdz
  gv(2) =  dg(1) - g*b_inv*dbdr
  gv(3) =  (r_inv*(ddp-dpdr*r_inv) - rb_inv*(dpdr*dbdr+dpdz*dbdz))*zsp
  gv(1:3) = zcf*rb_inv*gv(1:3)

  ! mod bs : bs dot b
  bs = g*(ddp - r_inv*dpdr) - (dpdr*dg(1) + dpdz*dg(2))
  bs = zcf*r2b2_inv*zsp*bs
  bs = b + bs

  ! bs
  bsv(1:3) = (bv(1:3) + gv(1:3))/bs

  zcf = zc/(b*bs*qs)
  g12 = g*r_inv
  g13 = dpdr*r2_inv*zsp
  g23 = dpdz*r2_inv*zsp

  df(1) = qs*dphi(1) + mu*dbdr
  df(2) = qs*dphi(2) + mu*dbdz
  df(3) = 0.0_fp

  dzzdt(1) = vl*bsv(1) + zcf*( g12*df(2) + g13*df(3))
  dzzdt(2) = vl*bsv(2) + zcf*(-g12*df(1) + g23*df(3))
  dzzdt(3) = -(bsv(1)*df(1) + bsv(2)*df(2) + bsv(3)*df(3))/ms

  ! conserved energy and momentum at present position
  orbrzv_teng = 0.5_fp*ms*orbrzv_zz(3)*orbrzv_zz(3) + orbrzv_zz(4)*orbrzv_b &
            & + qs*orbrzv_phiv(3)
  orbrzv_tcmt = ms*orbrzv_zz(3)*orbrzv_gv(1)/orbrzv_b &
            & - zsp*qs*orbrzv_psiv(1)/zc

end subroutine orbrz_derivs
