! integrate guiding center orbit in (r,z) coordinate
subroutine orbrzg(iorb, istep, ierr)

  use orbrzv
  use nbi_com,only:biz,bia,bixi,bithst,bi1ovb,orbrzv_rf_option,orbrzv_zzerr_con,&
       &goose,ngoose,bivpar,bivprp,bitime,xbmbnd,dxizc,ledge,bidt,bim,biv,bike,bimu,&
       &biksid,rmjion,yion,bimodb,nsnccw,zel,zmp

  use xplasma_obj_instance
  use logmod

  implicit none
  integer, parameter :: ion = 1
  integer, intent(in) :: iorb, istep ! orbit index and time step from orbit.for
  integer, intent(out) :: ierr ! status code, 0 = ok
  integer :: icnt, ierr_sum
  real(fp) :: zximin, zximax, dumj(2,2), dpsi(1:2), zg, zvl0, zr0
  integer :: jx, iant, inphi
  real(fp) :: zx0, zh, ziotb

  ! initialize error flag
  ierr  = 0
  ierr_sum = 0

  ! charge, mass, time step
  orbrzv_qs = zel*biz(ion)   ! statcoulombs
  orbrzv_ms = zmp*bia(ion)   ! grams
  
  
  ! ----------------------------------------------------------------------
  ! at initial time step, initialize r,z and related geometric variables
  ! otherwise, use the values set by previous time step
  if(il_rz_updated) then !r,z moved due to anomalous diffusion, update gc parameters
     orbrzv_zz(1) = rmjion
     orbrzv_zz(2) = yion
     ! re-evaluate mod(b) at current (r,z)
     call ezspline_gradient2(orbrzv_psi_spo, &
          & orbrzv_zz(1), orbrzv_zz(2), dpsi, ierr)
     call ezspline_interp2(orbrzv_g_spo, &
          & orbrzv_zz(1), orbrzv_zz(2), zg, ierr)
     orbrzv_b = dsqrt(dpsi(1)*dpsi(1)+zg*zg+dpsi(2)*dpsi(2))/orbrzv_zz(1)
     
     orbrzv_zz(3) = bivpar(ion)*(-nsnccw) !mg,04/22/16 (doesn't need it here for rz integrator)
     
     orbrzv_zz(4) = (0.5_fp*orbrzv_ms*bivprp(ion)**2)/orbrzv_b

     ! call this to initialize (orbrzv_rho, orbrzv_th, orbrzv_pflag)
     call orbrz_check_bdy(iorb, istep)
     if(orbrzv_pflag < 0) then
        call errorlog(' orbrzg: lost particle at the start of orbrz ?')
        ierr = 0 !lost orbit
        return
     end if
     ! initialize time derivatives
     ! ---> geometric quantities are also initialized inside
     !      the derivative function
     call orbrz_derivs(1, orbrzv_zz, orbrzv_ms, orbrzv_qs, orbrzv_dzzdt)
     bi1ovb(ion) = 1.0_fp/orbrzv_b ! may be redundant ?
  endif

  if(istep.eq.1) then

    !call xmomry(bixi(ion), bithst(ion), 0, orbrzv_zz(1), orbrzv_zz(2), dumj)
     if(il_rz_updated) then
        orbrzv_zz(1) = rmjion
        orbrzv_zz(2) = yion
     else
        call xplasma_ctrans(s, ierr, rho_in = bixi(ion), theta_in = bithst(ion), &
             & r_out=orbrzv_zz(1), z_out=orbrzv_zz(2))
        orbrzv_zz(1) = orbrzv_zz(1)*1.0e2_fp
        orbrzv_zz(2) = orbrzv_zz(2)*1.0e2_fp
     endif
         
    ! re-evaluate mod(b) at current (r,z)
    call ezspline_gradient2(orbrzv_psi_spo, &
                             & orbrzv_zz(1), orbrzv_zz(2), dpsi, ierr)
    call ezspline_interp2(orbrzv_g_spo, &
                           & orbrzv_zz(1), orbrzv_zz(2), zg, ierr)
    orbrzv_b = dsqrt(dpsi(1)*dpsi(1)+zg*zg+dpsi(2)*dpsi(2))/orbrzv_zz(1)

    orbrzv_zz(3) = bivpar(ion)*(-nsnccw) !mg,04/22/16 (doesn't need it here for rz integrator)

    orbrzv_zz(4) = (0.5_fp*orbrzv_ms*bivprp(ion)**2)/orbrzv_b

    ! call this to initialize (orbrzv_rho, orbrzv_th, orbrzv_pflag)
    call orbrz_check_bdy(iorb, istep)
    if(orbrzv_pflag < 0) then
       call errorlog(' orbrzg: lost particle at the start of orbrz ?')
       ierr = 0 !lost orbit
       return
    end if

    ! initialize time derivatives
    ! ---> geometric quantities are also initialized inside
    !      the derivative function
    call orbrz_derivs(1, orbrzv_zz, orbrzv_ms, orbrzv_qs, orbrzv_dzzdt)
    bi1ovb(ion) = 1.0_fp/orbrzv_b ! may be redundant ?

    ! set initial trial time step
    !orbrzv_tdt = 10.0_fp*orbrzv_dt_min
    orbrzv_tdt = orbrzv_dt_min

    ! from zdtrob.for
!    bim(ion)    = orbrzv_ms
!    bimu(ion)   = orbrzv_zz(4)
!    jx  = int(bixi(ion)/dxi(lcentr,2)) + 1
!    zx0 = real(jx-1, r8)*dxi(lcentr,2)
!    zh  = bixi(ion) - zx0
!    ziotb = zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx),spliot(4,jx))
!    biq(ion) = 1.0_fp/ziotb
    !orbrzv_tdt = zdtrob(ierr)/5.0_fp
!    orbrzv_tdt = zdtrob(ierr)
!    if(ierr.ne.0) then
!       call errorlog(' orbrzg error finding integration time step')
!       orbrzv_tdt =0.0_fp
!       return
!    endif
!    orbrzv_tdt =fdtnxy*orbrzv_tdt   !half the usual time step
!    orbrzv_tdt = min(orbrzv_dt_max,max(orbrzv_dt_min,orbrzv_tdt))

    ! non-goosed time counted by orbrzv_dt
    orbrzv_time0 = 0.0_fp

    ! initial rf-resonance info
    if(orbrzv_rf_option.eq.1 .and. mcrf_on) then
      do iant = 1, orbrzv_rf_nant
        do inphi = 1, orbrzv_rf_nnphi(iant)
          if(orbrzv_rfw(iant,inphi)%freq.le.0.0_fp) cycle !mg, no data for this antenna
          orbrzv_rfw(iant,inphi)%time  = bitime(1)
          orbrzv_rfw(iant,inphi)%time0 = orbrzv_time0
          orbrzv_rfw(iant,inphi)%dt    = 0.0_fp
          orbrzv_rfw(iant,inphi)%dt0   = 0.0_fp
          orbrzv_rfw(iant,inphi)%nstep = 0.0_fp
          orbrzv_rfw(iant,inphi)%ztf   = orbrzv_rf_dist
          orbrzv_rfw(iant,inphi)%dsgn=0
          orbrzv_rfw(iant,inphi)%lam0=0.0_fp
          orbrzv_rfw(iant,inphi)%lam=0.0_fp
        enddo
      enddo

      ! initialize fbm diag. per particle
      !orbrzv_fbm_cnt = 0
      !orbrzv_fbm_tmp = 0.0_fp

      ! accumulator for additional goose control for rf
      orbrzv_ztf_min  = 0.0_fp
      orbrzv_zrfd_min = 1.0_fp
      orbrzv_rfdv(-1:2) = 0.0_fp

      ! kjm-rfdbg-nogoose
      if(nlmcrfgcon) then
        goose = 1.0001_fp
        ngoose = 0
      endif
      orbrzv_nmidc = 0        ! mid-plane crossing counter
      orbrzv_bnw = 0.0_fp     ! orbit width
      orbrzv_bnw_lrf = 0.0_fp ! orbit width / rf-scale length
      orbrzv_trpas = 0.0_fp   ! flag to trap or passing

      ! backup error controler
      orbrzv_zzerr_con0 = orbrzv_zzerr_con
    endif

  endif !end of step.eq.1

  if(ierr .ne. 0) then
     call errorlog(' orbrzg: error in preparing orbrz')
     ierr_sum = ierr_sum + ierr
  endif

  ! if not the initial time, just reset velocity varialbes which might
  ! have been changed by cx, fpp...
  orbrzv_zz(3) = bivpar(ion)*(-nsnccw) !mg,04/22/16 (doesn't need it here for rz integrator)
  orbrzv_zz(4) = (0.5_fp*orbrzv_ms*bivprp(ion)**2)/orbrzv_b

  ! go one step
  call orbrz4(iorb, istep, 0, ierr)   ! rk4
  !call orbrz23(iorb, istep, 0, ierr) ! rk2 & rk3


  if(ierr .ne. 0) then
     call errorlog(' orbrzg: error after orbrz')
     ierr_sum = ierr_sum + ierr
  endif
 
  ! apply rf-heating !
  orbrzv_dbg_dt=0.0_fp ! use this var. to calculate the total resonance duration time
  if(orbrzv_rf_option.eq.1 .and. mcrf_on) then
    if(orbrzv_pflag .eq. 0) then
      ! check rf-resonance condition
      call orbrz_check_rfres(iorb, istep)
      ! apply quasi-linear heating
      call orbrz_qlrf(iorb, istep, ierr)
    endif
  endif

  !if(orbrzv_pflag .ne. 0) then
  !  print *, "bad orbit", iorb, istep, orbrzv_rho, orbrzv_zz(3), orbrzv_zz(4)
  !endif

  ! limiter of xi value : from orbxit.for
  zximin = 1.0e-5_fp
  zximax = xbmbnd + 0.5_fp*dxizc(ledge)

  ! store time step actually undergone in this step
  bidt(1) = orbrzv_dt+orbrzv_dbg_dt

  ! update non-goosed time elaspe
  orbrzv_time0 = orbrzv_time0 + orbrzv_dt

  ! set other legacy variables
  bim(ion)    = orbrzv_ms
  biv(ion)    = dsqrt(2.0_fp*(orbrzv_teng-orbrzv_qs*orbrzv_phiv(3))/orbrzv_ms)
  bike(ion)   = 0.5_fp*bim(ion)*biv(ion)**2
  bimu(ion)   = orbrzv_zz(4)
  bivpar(ion) = orbrzv_zz(3)*(-nsnccw) !mg,04/22/16 (back to orbpre.for definition)
  bivprp(ion) = dsqrt(2.0_fp*orbrzv_zz(4)*orbrzv_b/orbrzv_ms)
  biksid(ion) =  bivpar(ion)/biv(ion)
  rmjion      = orbrzv_zz(1)
  yion        = orbrzv_zz(2)
  bixi(ion)   = orbrzv_rho
  bixi(ion)   = dmax1(zximin,dmin1(bixi(ion),zximax))
  bithst(ion) = orbrzv_th

  bimodb(ion) = orbrzv_b*1.0e-4_fp ! in tesla
  bi1ovb(ion) = 1.0_fp/orbrzv_b

  if(ierr .ne. 0) then
     call errorlog(' orbrzg error after qlrf and finalizing')
    ierr_sum = ierr_sum + ierr
  endif

  ierr = ierr_sum

end subroutine orbrzg

! call old time step controller from orbrz*
! it is assumed that 
!
!   orbrzv_zz(1:4)
!   orbrzv_rho, orbrzv_th
!   orbrzv_b, orbrzv_ms, orbrzv_qs
!
! are initialized
subroutine orbrz_safe_call_zdtrob(zdtr)

  use orbrzv, only : orbrzv_ms,orbrzv_zz,orbrzv_b,orbrzv_rho,&
       orbrzv_th
  use nbi_com

  implicit none
  real(fp), intent(out) :: zdtr
  real(fp) :: zximin, zximax
  real(fp) :: zm0, zmu0, zb01, zv0, zrho0, zth0, zbq0
  integer :: ierr, jx
  real(fp) :: zx0, zh, ziotb
  real(fp) :: zf1d, zdtrob

  ! limiter of xi value : from orbxit.for
  zximin = 1.0e-5_fp
  zximax = xbmbnd + 0.5_fp*dxizc(ledge)

  ! prepare to call zdtrob to compare time step size
  zm0 = bim(1) 
  zmu0 = bimu(1) 
  zv0 = biv(1) 
  zb01 = bi1ovb(1)
  zrho0 = bixi(1) 
  zth0 = bithst(1)
  zbq0 = biq(1)

  bim(1)  = orbrzv_ms
  bimu(1) = orbrzv_zz(4)
  biv(1) = dsqrt(orbrzv_zz(3)**2 + 2.0_fp*orbrzv_zz(4)*orbrzv_b/orbrzv_ms)
  bi1ovb(1) = 1.0_fp/orbrzv_b
  bixi(1) = orbrzv_rho
  bithst(1) = orbrzv_th
  
  jx  = int(bixi(1)/dxi(lcentr,2)) + 1
  zx0 = real(jx-1,fp)*dxi(lcentr,2)
  zh  = bixi(1) - zx0
  ziotb = zf1d(zh,spliot(1,jx),spliot(2,jx),spliot(3,jx),spliot(4,jx))
  biq(1) = 1.0_fp/ziotb

  ! call zdtrob !
  zdtr = zdtrob(ierr)

  ! recover legacy variables
  bim(1) = zm0
  bimu(1) = zmu0
  biv(1) = zv0
  bi1ovb(1) = zb01
  bixi(1) = zrho0
  bithst(1) = zth0
  biq(1) = zbq0

end subroutine orbrz_safe_call_zdtrob
