subroutine orbrz_init_spline_coeffs

  use ezspline, only : ezspline_free,ezspline_init,ezspline_error,&
       ezspline_setup,ezspline_interp2,ezspline_gradient2,&
       ezspline_derivative2
  use orbrzv
  use nbi_com,only : orbrzv_nr,orbrzv_nz,id_psi,id_bphi,id_epot,&
       nsbeam,lsbeam,ebdmaxa,nsnccw
  use xplasma_obj_instance
  use logmod

  implicit none
  integer bcs_r(1:2), bcs_z(1:2), ier, i, j, flist(1:3)
  real(fp) :: fval(1:1,1:3), r, z, dumr, psi, g
  real(fp), allocatable, dimension(:,:) :: psi_arry, g_arry, phi_arry
!mg 11/16 cleanup  real(fp), allocatable, dimension(:) :: zbuf1
  real(fp) :: psiv(1:6), gv(1:3), dpsi(2), dpdrr, dpdzz, dpdrz, dg(2), &
            & br, bphi, bz, b, b_inv, dbdr, dbdz, r_inv, phi, &
            & zdrhop, zrhop, zrhot
  integer :: inrhop
        
! mg 11/2016 i intensively checked and compared orbg.for(flux coordinate) 
! and orbrzg.f90(rz integrator). all looks good.
  ! kjm-rfdbg --------------------------------------------- 
  !orbrzv_nr = 160
  !orbrzv_nz = 160
  ! kjm-rfdbg --------------------------------------------- 

  ! (orbrzv_rmin, orbrzv_rmax) x (orbrzv_zmin, orbrzv_zmax)
  ! got from 'xstraln/xstraln_init.f90'
  ! call eq_bdlims(itype, orbrzv_rmin, orbrzv_rmax, &
  !                    & orbrzv_zmin, orbrzv_zmax, ier)
  call eq_rzlims(orbrzv_rmin, orbrzv_rmax, &
               & orbrzv_zmin, orbrzv_zmax, ier)

  ! conversion from mks to cgs
  orbrzv_rmin = 1.0e2_fp*orbrzv_rmin
  orbrzv_rmax = 1.0e2_fp*orbrzv_rmax
  orbrzv_zmin = 1.0e2_fp*orbrzv_zmin
  orbrzv_zmax = 1.0e2_fp*orbrzv_zmax

  orbrzv_dr = (orbrzv_rmax-orbrzv_rmin)/(real(orbrzv_nr-1,fp))
  orbrzv_dz = (orbrzv_zmax-orbrzv_zmin)/(real(orbrzv_nz-1,fp))


  ! allocate grid array and temporary function value arrays
  allocate(psi_arry(1:orbrzv_nr, 1:orbrzv_nz), &
         & g_arry(1:orbrzv_nr, 1:orbrzv_nz), &
         & phi_arry(1:orbrzv_nr, 1:orbrzv_nz))


  ! free previously allocated spline object
  call ezspline_free(orbrzv_psi_spo, ier)
  call ezspline_free(orbrzv_g_spo, ier)
  call ezspline_free(orbrzv_phi_spo, ier)
  call ezspline_free(orbrzv_bdb_spo, ier)
  call ezspline_free(orbrzv_b_spo, ier)
!mg 11/16 cleanup  call ezspline_free(orbrzv_rhot_spo, ier)

  ! set boundary condition for spline interlation
  bcs_r = (/ 0, 0 /)
  bcs_z = (/ 0, 0 /)

  ! initialize spline object
  ! for rz-integrator
  call ezspline_init(orbrzv_psi_spo, orbrzv_nr, orbrzv_nz, bcs_r, bcs_z, ier)
  call ezspline_error(ier)
  call ezspline_init(orbrzv_g_spo, orbrzv_nr, orbrzv_nz, bcs_r, bcs_z, ier)
  call ezspline_error(ier)
  call ezspline_init(orbrzv_phi_spo, orbrzv_nr, orbrzv_nz, bcs_r, bcs_z, ier)
  call ezspline_error(ier)
  !read *, dumr
  ! for rf calculation
  call ezspline_init(orbrzv_bdb_spo, orbrzv_nr, orbrzv_nz, bcs_r, bcs_z, ier)
  call ezspline_error(ier)
  call ezspline_init(orbrzv_b_spo, orbrzv_nr, orbrzv_nz, bcs_r, bcs_z, ier)
  call ezspline_error(ier)

  ! set rz-domain
  do i=1,orbrzv_nr
    orbrzv_psi_spo%x1(i) = orbrzv_rmin + orbrzv_dr*(real(i-1,fp))
    orbrzv_g_spo%x1(i)   = orbrzv_rmin + orbrzv_dr*(real(i-1,fp))
    orbrzv_phi_spo%x1(i) = orbrzv_rmin + orbrzv_dr*(real(i-1,fp))
  end do
  do j=1,orbrzv_nz
    orbrzv_psi_spo%x2(j) = orbrzv_zmin + orbrzv_dz*(real(j-1,fp))
    orbrzv_g_spo%x2(j)   = orbrzv_zmin + orbrzv_dz*(real(j-1,fp))
    orbrzv_phi_spo%x2(j) = orbrzv_zmin + orbrzv_dz*(real(j-1,fp))
  end do
  !read *, dumr
  
  ! get psi & g values on the 2d grid points
  flist = (/ id_psi, id_bphi, id_epot /)
  ! kjm-rfdbg -------------------------------
  !call eq_gfnum('g', id_g)
  !call eqm_rzfun(id_g, 1.0_fp, 2, 0.0_fp, ier)
  !flist = (/ id_psi, id_g /)
  ! kjm-rfdbg -------------------------------

  do i=1,orbrzv_nr
    r = orbrzv_rmin + (real(i-1,fp))*orbrzv_dr
    do j=1,orbrzv_nz
      z = orbrzv_zmin + (real(j-1,fp))*orbrzv_dz

      call eq_frz(1, r*1.0e-2_fp, z*1.0e-2_fp, 3, flist, 1, fval, ier)

      psi_arry(i,j) = fval(1, 1)
      g_arry(i,j)   = fval(1, 2)*r   ! g = r*btor , t*cm
      phi_arry(i,j) = fval(1, 3)

      ! conversion from mks to cgs
      psi_arry(i,j) = psi_arry(i,j)*1.0e8_fp
      g_arry(i,j)   = g_arry(i,j)*1.0e4_fp !gauss *cm
      phi_arry(i,j) = phi_arry(i,j)/3.0e2_fp !volts tp statvolts

    enddo
  enddo
  !print *,"orbrzv : psi(mid) ", psi_arry(orbrzv_nr/2,orbrzv_nz/2)
  !print *,"orbrzv : g(mid) ", g_arry(orbrzv_nr/2,orbrzv_nz/2)
  !read *, dumr
  ! setup spline coefficients 
  call ezspline_setup(orbrzv_psi_spo, psi_arry, ier)
  call ezspline_error(ier)
  call ezspline_setup(orbrzv_g_spo, g_arry, ier)
  call ezspline_error(ier)
  call ezspline_setup(orbrzv_phi_spo, phi_arry, ier)
  call ezspline_error(ier)



  ! some global information from xplasma ------------------------------------
  call eq_rgetf(1, 1.0_fp, id_psi, 0, orbrzv_psix, ier)
  call eq_sgaxis(orbrzv_bc, orbrzv_rc, orbrzv_zc)
  call xplasma_rzminmax_plasma(s, orbrzv_prmin, orbrzv_prmax, &
                             & orbrzv_pzmin, orbrzv_pzmax, ier)
  ! conversion from mks to cgs
  orbrzv_psix = orbrzv_psix*1.0e8_fp
  orbrzv_bc = orbrzv_bc*1.0e4_fp
  orbrzv_rc = orbrzv_rc*1.0e2_fp
  orbrzv_zc = orbrzv_zc*1.0e2_fp
  orbrzv_prmin = orbrzv_prmin*1.0e2_fp
  orbrzv_prmax = orbrzv_prmax*1.0e2_fp
  orbrzv_pzmin = orbrzv_pzmin*1.0e2_fp
  orbrzv_pzmax = orbrzv_pzmax*1.0e2_fp

  ! determine bphi sign
  call ezspline_interp2(orbrzv_g_spo, orbrzv_rc, orbrzv_zc, g, ier)
  orbrzv_bpsign = sign(1.0_fp, g) 
  ! orbrzv_bpsign and nsnccw must be same
  if(orbrzv_bpsign*nsnccw.lt.0.0_fp) then
     ! redefine the sign of 'g'
     call ezspline_setup(orbrzv_g_spo, -g_arry, ier)
     call ezspline_error(ier)
     call ezspline_interp2(orbrzv_g_spo, orbrzv_rc, orbrzv_zc, g, ier)
     orbrzv_bpsign = sign(1.0_fp, g) 
  endif

  ! geometric quantities derived from psi, g, phi ---------------------------
  do i= 1, orbrzv_nr 
    orbrzv_bdb_spo%x1(i) = orbrzv_rmin + orbrzv_dr*(real(i-1,fp))
    orbrzv_b_spo%x1(i) = orbrzv_rmin + orbrzv_dr*(real(i-1,fp))
  end do
  do j= 1, orbrzv_nz
    orbrzv_bdb_spo%x2(j) = orbrzv_zmin + orbrzv_dz*(real(j-1,fp))
    orbrzv_b_spo%x2(j) = orbrzv_zmin + orbrzv_dz*(real(j-1,fp))
  end do

  ! store minimum of b dot db
  orbrzv_bdb_min = orbrzv_bc/orbrzv_rc

  do i=1,orbrzv_nr
    r = orbrzv_rmin + (real(i-1,fp))*orbrzv_dr
    r_inv = 1.0_fp/r

    do j=1,orbrzv_nz
      z = orbrzv_zmin + (real(j-1,fp))*orbrzv_dz

      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)

      br =  dpsi(2)*r_inv*orbrzv_bpsign
      bz = -dpsi(1)*r_inv*orbrzv_bpsign
      bphi = g*r_inv

      b = dsqrt(br*br + bphi*bphi + bz*bz)
      b_inv = 1.0_fp/b

      dbdr = r_inv*r_inv*b_inv*(dpsi(1)*dpdrr+dpsi(2)*dpdrz+g*dg(1)) - b*r_inv
      dbdz = r_inv*r_inv*b_inv*(dpsi(2)*dpdzz+dpsi(1)*dpdrz+g*dg(2))

      ! b dot db
      psi_arry(i,j) = (br*dbdr + bz*dbdz)*b_inv
      orbrzv_bdb_min = dmin1(orbrzv_bdb_min, dabs(psi_arry(i,j)))

      ! mod(b)
      phi_arry(i,j) = b
    end do
  end do
  call ezspline_setup(orbrzv_bdb_spo, psi_arry, ier)
  call ezspline_error(ier)
  call ezspline_setup(orbrzv_b_spo, phi_arry, ier)
  call ezspline_error(ier)

  orbrzv_bdb_min = max(1.0e-6*orbrzv_bc/orbrzv_rc, orbrzv_bdb_min)
  !orbrzv_bdb_min = max(1.0e-6*orbrzv_bc/orbrzv_rc, orbrzv_bdb_min)

  ! kjm-constrf
  !orbrzv_bdb_min = 1.0e-1*orbrzv_bc/orbrzv_rc


  ! setup rhot(rhop) interpolation
  inrhop = 64
  bcs_r = (/ 0, 0 /)
!mg 11/16 cleanup  call ezspline_init(orbrzv_rhot_spo, inrhop, bcs_r, ier)
!mg 11/16 cleanup  call ezspline_error(ier)

  zdrhop = 1.0_fp/real(inrhop-1,fp)
!mg 11/16 cleanup  allocate(zbuf1(1:inrhop))
!mg 11/16 cleanup
!mg 11/16 cleanup  do i = 1, inrhop
!mg 11/16 cleanup    zrhot = real(i-1,fp)*zdrhop
!mg 11/16 cleanup
!mg 11/16 cleanup    call xplasma_ctrans(s,ier,rho_in=zrhot,theta_in=0.0_fp,r_out=r,z_out=z)
!mg 11/16 cleanup    r = r*1.0e2_fp
!mg 11/16 cleanup    z = z*1.0e2_fp
!mg 11/16 cleanup    call ezspline_interp2(orbrzv_psi_spo, r, z, psi, ier)
!mg 11/16 cleanup    zrhop = psi
!mg 11/16 cleanup
!mg 11/16 cleanup    orbrzv_rhot_spo%x1(i) = zrhop
!mg 11/16 cleanup    zbuf1(i) = zrhot
!mg 11/16 cleanup  enddo
!mg 11/16 cleanup  call ezspline_setup(orbrzv_rhot_spo, zbuf1, ier)
!mg 11/16 cleanup  call ezspline_error(ier)


  ! free local allocation
  deallocate(psi_arry, g_arry, phi_arry)
!mg 11/16 cleanup  deallocate(zbuf1)
  !read *, dumr

  ! print out geometric info for debugging ----------------------------- !
#ifndef __MPI
!#ifdef __DEBUG
  orbrzv_out = 314159
  open(unit=orbrzv_out, file="rzb.txt")

  do i = 1, nsbeam
    write (unit=orbrzv_out, fmt=*) "# emax ion ", i, " : ", ebdmaxa(i)
  end do

  write (unit=orbrzv_out, fmt=*) "# sim box r : ", orbrzv_rmin, orbrzv_rmax
  write (unit=orbrzv_out, fmt=*) "# sim box z : ", orbrzv_zmin, orbrzv_zmax
  write (unit=orbrzv_out, fmt=*) "# r min & max : ", orbrzv_prmin,orbrzv_prmax
  write (unit=orbrzv_out, fmt=*) "# z min & max : ", orbrzv_pzmin,orbrzv_pzmax
  write (unit=orbrzv_out, fmt=*) "# psix, rc, zc, bc : ", &
    & orbrzv_psix, orbrzv_rc, orbrzv_zc, orbrzv_bc
  write (unit=orbrzv_out, fmt=*) "# bdb_min : ", orbrzv_bdb_min

  do i=1,orbrzv_nr
    r = orbrzv_rmin + (real(i-1,fp))*orbrzv_dr
    do j=1,orbrzv_nz
      z = orbrzv_zmin + (real(j-1,fp))*orbrzv_dz
      call ezspline_interp2(orbrzv_psi_spo, r, z, psi, ier)
      call ezspline_interp2(orbrzv_g_spo, r, z, g, ier)
      call ezspline_interp2(orbrzv_phi_spo, r, z, phi, ier)
      call ezspline_interp2(orbrzv_b_spo, r, z, b, ier)
      call ezspline_interp2(orbrzv_bdb_spo, r, z, dbdr, ier)
      write (unit=orbrzv_out, fmt=*) r, z, psi, g/r, b, dbdr, phi
    enddo
    write (unit=orbrzv_out, fmt=*) " "
  enddo
  close(unit=orbrzv_out)
  call infolog(' orbrz_init_spline_coeffs:orbrzv_out completed')
!#endif
#endif
  ! -------------------------------------------------------------------- !

end subroutine orbrz_init_spline_coeffs
