subroutine loadfki
  use particle_array
  use grids_eqdata
  use precision
  use spline_function
  use equilibrium
  use global_parameters
  use read_in_parameters
  implicit none

  integer mtest, ierror, ngyrofki

  interface
  subroutine loadfk_tracking(mp, mpmax, zpart, zpart0, wzpart, wppart, jtpart0, jtpart1, wtpart0, wtpart1, jppart,&
    qpart, apart)
    use precision
    use spline_function
    use equilibrium
    !use global_parameters, only:
    use read_in_parameters
 
    integer mtest, ierror, nparamfki
    real(lk) b, vth_norm, vthpart, psidum, particle_charge, paticle_mass, norm_energy, qpart, apart
    real(lk), dimension(:,:) :: zpart, zpart0, wppart, wtpart0, wtpart1
    integer, dimension(:,:) :: jtpart0, jtpart1, jppart
    real(lk), dimension(:) ::wzpart
  end subroutine
  
  subroutine loadfk(mp, mpmax, zpart, zpart0, wzpart, wppart, jtpart0, jtpart1,&
                  wtpart0, wtpart1, jppart,rhogyropart, t_Eq, apart, qpart, marker, &
                  pmark, meshn, mesht, nonuni, marknum, mpcore, mpsol)
  use precision
  use global_parameters, only: npartdom, mype, gtcout,numberpe,deltaz,zeta0,nmode,pi,&
                               pi2,istep,irk
  use mapping_splines, only: rspsol_psi_s, zspsol_psi_s, sspsol_psi_z
  use grids_setup
  use grids_eqdata
  use read_in_parameters
  use spline_function, only: spline2d, sppsiR1d
  use cylindrical_function, only: geotheta
  use equilibrium, only: psirz_Eq, B0Mag_Eq, r_Eq, z_Eq,&
                         psi_sep, lsr, lsz, spdR, spdZ, zmax, B0R_Eq, B0Z_Eq
  use simulation_domain, only: sfcore, sfsol, couplefrc, psi1_core, psi0_core
  use mapping_splines, only: msp_psi_core, msp_psi_sol,&
                             rspcore_psi_s, zspcore_psi_s, sspcore_psi_theta, &
                             rspcore_psi_theta, zspcore_psi_theta,&
                             rspsol_psi_s, zspsol_psi_s, &
                             lspcore, lstcore, del_psi_core, del_s_core, del_theta_core, &
                             lspsol, lszsol, del_psi_sol, del_s_sol
  use cell_volume, only: volume
  implicit none
  ! dummy parameters
  real(lk), dimension(:,:,:) :: t_Eq
  real(lk), dimension(:,:) ::zpart, zpart0, wppart, wtpart0, wtpart1,rhogyropart
  integer, dimension(:,:) :: jtpart0, jtpart1, jppart
  real(lk), dimension(:) :: wzpart, marker, meshn, mesht, nonuni, pmark
  integer, dimension(:) :: marknum
  integer mp, mpmax, mpcore, mpsol
  real(lk) apart, qpart

  ! local variables
  integer :: i, j, ij, ij0, ij1, ii, m, mtest, ierror, icount, mmin, mmax, nzsp, nz, np,&
             ipsi(1), iz0(1)
  real(lk) :: c0=2.515517, c1=0.802853, c2=0.010328, d1=1.432788, d2=0.189269, d3=0.001308
  real(lk) :: w_initial, rhop, b, umax, pdum, tdum, rmin1, rmax1, rmin2, rmax2, zzmin, zzmax,&
              rdum, zdum, dpr, dpz, del_r, del_z, rdum0, zdum0, vth, anonuni, rg, psi_cy0,&
              psi_cy1, deltapsi_cy, psidum, delta_V, sdum(mpsi), zdum2, dz2, psifrc, del_p, &
              jaco_frc, dr1, dr2, dr3, dr11, dr22, dr33, r0tmp, r1tmp, r2tmp, psidumup, &
              psidumlow, pmin, pmax, psidumtmp, dpsidrtmp, spdR_inv, spdZ_inv, &
              dpsi, dp1, del_s, ssmin, ssmax, ldum, wp1, wp0, wt00, wt10, wt01, wt11, &
              jj0, dz0, ndum, del_psidum, psi_inner_boundary, psi_outer_boundary, &
              vthp, vth_norm, proton_mass, proton_charge, norm_energy, tp, cmratio, &
              br, bz, bt, vecb(3), vec1(3), vec2(3), vecv(3), vecrho(3),vperp,vpara,omegac
  real(lk), dimension(:), allocatable:: alpha
  end subroutine loadfk
  end interface

  ngyrofki=1
  !(R,\zeta,Z,vR,weight,v\zeta,psi,vZ,energy,real \zeta,pos(1:3)[, tracking #1, tracking #2])
  nparamfki=13
  if(track_particles==1) nparamfki=nparamfki+2
  afki=aion
  qfki=qion
  mfkimax=1*mfki+500*ceiling(sqrt(real(mfki)))

  allocate(zfki(nparamfki,mfkimax),zfki0(nparamfki,mfkimax),wzfki(mfkimax),wpfki(ngyrofki,mfkimax),&
          jtfki0(ngyrofki,mfkimax),jtfki1(ngyrofki,mfkimax),wtfki0(ngyrofki,mfkimax),rhogyrofki(ngyrofki,mfkimax),rhogyrofki0(ngyrofki,mfkimax),&
          wtfki1(ngyrofki,mfkimax),jpfki(ngyrofki,mfkimax),pmarkfki(mpsi),markerfki(mgrid),&
          densityfki(0:1,mgrid),flowfki(0:1,mgrid),densityfki0(0:1,mgrid),&
          nonunifki(mgrid),marknumfki(mgrid), pressurefkipara(0:1,mgrid),pressurefkiperp(0:1,mgrid),&
          zonalfki(mpsi),zonalcfki(mpsi), pfluxfki(mpsi), rdtemfki(mpsi),&
          data1dfki(mpsi,mpdata1d),STAT=mtest)
  if(mtest /= 0) then
    write(0,*)mype,'*** Cannot allocate ion: mtest=',mtest
    call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
  endif
  pfluxfki=0.0_lk
  rdtemfki=0.0_lk
  data1dfki=0.0_lk

  call loadfk(mfki, mfkimax, zfki, zfki0, wzfki, wpfki, jtfki0, jtfki1,&
                  wtfki0, wtfki1, jpfki,rhogyrofki, ti_Eq, afki, qfki, markerfki, &
                  pmarkfki, meshni, meshti, nonunifki, marknumfki, mfkicore, mfkisol)
  if(track_particles==1) call loadfk_tracking(mfki, mfkimax, zfki, zfki0, wzfki, wpfki, jtfki0, jtfki1,wtfki0, wtfki1, jpfki,&
       qfki, afki)
      if(irun==0)then
#ifndef GPU_UM
      !$acc update device(zfki,zfki0,zfki1)
#endif
    endif

  if(mype==0) write(gtcout, *)'mfki=',mfki
end subroutine loadfki

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine loadfk_tracking(mp, mpmax, zpart, zpart0, wzpart, wppart, jtpart0, jtpart1, wtpart0, wtpart1, jppart,&
  qpart, apart)
  use precision
  use spline_function
  use equilibrium
  use global_parameters
  use read_in_parameters
  use simulation_domain
 
  real(lk) b, vth_norm, vthpart, psidum, particle_charge, particle_mass, norm_energy, qpart, apart, &
  vpara, vperp, vec1(3), vec2(3), vecb(3),br,bz,bt,alpha, vecv(3),rho(3), omegac
  real(lk), dimension(:,:) :: zpart, zpart0, wppart, wtpart0, wtpart1
  integer, dimension(:,:) :: jtpart0, jtpart1, jppart
  integer mp, mpmax,ierror
  real(lk), dimension(:) ::wzpart

  ! allocate particle memory
  fkitemp0=1.0_lk
  zpart(1,1)=2.0564_lk
  zpart(2,1)=zeta0+0.5_lk*deltaz
  zpart(3,1)=0.0_lk
  ! B field at guding center position
  b = spline2d(zpart(1,1),zpart(3,1)+zmax,0,lsr,lsz,spdR,spdZ,B0mag_Eq,0,0)
  br=spline2d(zpart(1,1),zpart(3,1)+zmax,0,lsr,lsz,spdR,spdZ,B0R_Eq,0,0)
  bz=spline2d(zpart(1,1),zpart(3,1)+zmax,0,lsr,lsz,spdR,spdZ,B0Z_Eq,0,0)
  bt=sqrt(b**2-br**2-bz**2)
  vth_norm = ulength*0.01_lk/utime !m/s
  vthpart = 9790.0_lk*sqrt(0.555_lk*etemp0*fkitemp0)/sqrt(apart) !m/s
  ! particle velocity and position
  vperp=vthpart/vth_norm
  vpara=vthpart/vth_norm
  ! alpha is the gyro phase angle
  call random_number(alpha)
  alpha = (alpha -0.5_lk) * pi 
  vecb=(/br/b, bt/b, bz/b/)
  vec1=(/bz/b*br/sqrt(br**2+bt**2), bz/b*bt/sqrt(br**2+bt**2), -sqrt(1-bz**2/b**2)/)
  call cross_product3(vecb, vec1, vec2)
  vecv=vperp*cos(alpha)*vec2+vperp*sin(alpha)*vec1+vpara*vecb
  ! rho = vecv x vecb / Omegac, x=X-rho
  omegac=qpart*b/apart
  rho=(vperp*cos(alpha)*vec1-vperp*sin(alpha)*vec2)/omegac
  zpart0(11,1)=zpart(1,1)-rho(1)
  zpart0(13,1)=zpart(3,1)-rho(3)
  zpart0(12,1)=zpart(2,1)
  zpart(10,1)=zpart0(12,1)
  zpart0(12,1)=modulo(zpart0(12,1),torbound)
  zpart(10,1)=modulo(zpart(10,1),pi2)
  zpart0(13,1)=modulo((zpart0(13,1)+zmax),(zmax*2.0_lk))-zmax
  psidum = spline2d(zpart0(11,1),zpart0(13,1)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
  zpart0(7,1)=1.0_lk
  zpart(5,1)=1.0e-3
  !initial particle position at -dt/2
  zpart(11,1)=zpart0(11,1)-vecv(1)*tstep*0.5
  zpart(12,1)=zpart0(12,1)-vecv(2)/zpart0(11,1)*tstep*0.5
  zpart(13,1)=zpart0(13,1)-vecv(3)*tstep*0.5
  zpart(13,1)=modulo((zpart(13,1)+zmax),(zmax*2.0_lk))-zmax
  zpart(4,1)=vecv(1)*cos(zpart0(12,1)-zpart(12,1))-vecv(2)*sin(zpart0(12,1)-zpart(12,1))
  zpart(6,1)=vecv(1)*sin(zpart0(12,1)-zpart(12,1))+vecv(2)*cos(zpart0(12,1)-zpart(12,1))
  zpart(6,1)=zpart(6,1)/zpart(11,1)
  zpart(8,1)=vecv(3)
  zpart(9,1)=0.5_lk*apart*(zpart(4,1)*zpart(4,1)+zpart(6,1)*zpart(6,1)*zpart(11,1)*zpart(11,1)+zpart(8,1)*zpart(8,1))
  particle_charge=qpart*1.602e-19 !C
  particle_mass=apart*1.673e-27 !kg
  norm_energy=particle_mass*vth_norm*vth_norm/particle_charge   ! eV

  if(mype==0)then
    write(gtcout,*) '----------info of tracking particle----------'
    write(gtcout,*) 'vperp', vperp,'vpara',vpara
    write(gtcout,*) 'vecv', vecv, 'v', sqrt(vecv(1)**2+vecv(2)**2+vecv(3)**2), 'energy', vecv(1)**2+vecv(2)**2+vecv(3)**2
    write(gtcout,*) 'B', b, 'rho', sqrt(rho(1)**2+rho(2)**2+rho(3)**2)
    write(gtcout,*) '-----------------------end-------------------'
  endif

end subroutine loadfk_tracking
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine cross_product3(A, B, C)
  use precision, only : lk
  implicit none
  real(lk) A(3), B(3), C(3)
  C(1)=A(2)*B(3)-A(3)*B(2)
  C(2)=A(3)*B(1)-A(1)*B(3)
  C(3)=A(1)*B(2)-A(2)*B(1)
end subroutine cross_product3

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine loadfk(mp, mpmax, zpart, zpart0, wzpart, wppart, jtpart0, jtpart1,&
                  wtpart0, wtpart1, jppart,rhogyropart, t_Eq, apart, qpart, marker, &
                  pmark, meshn, mesht, nonuni, marknum, mpcore, mpsol)
  use precision
  use global_parameters, only: npartdom, mype, gtcout,numberpe,deltaz,zeta0,nmode,pi,&
                               pi2,istep,irk,torbound
  use mapping_splines, only: rspsol_psi_s, zspsol_psi_s, sspsol_psi_z
  use grids_setup
  use grids_eqdata
  use read_in_parameters
  use spline_function, only: spline2d, sppsiR1d
  use cylindrical_function, only: geotheta
  use equilibrium, only: psirz_Eq, B0Mag_Eq, r_Eq, z_Eq,&
                         psi_sep, lsr, lsz, spdR, spdZ, zmax, B0R_Eq, B0Z_Eq
  use simulation_domain, only: sfcore, sfsol, couplefrc, psi1_core, psi0_core
  use mapping_splines, only: msp_psi_core, msp_psi_sol,& 
                             rspcore_psi_s, zspcore_psi_s, sspcore_psi_theta, &
                             rspcore_psi_theta, zspcore_psi_theta,&
                             rspsol_psi_s, zspsol_psi_s, &
                             lspcore, lstcore, del_psi_core, del_s_core, del_theta_core, &
                             lspsol, lszsol, del_psi_sol, del_s_sol
  use cell_volume, only: volume
  implicit none
  ! dummy parameters
  real(lk), dimension(:,:,:) :: t_Eq
  real(lk), dimension(:,:) ::zpart, zpart0, wppart, wtpart0, wtpart1, rhogyropart
  integer, dimension(:,:) :: jtpart0, jtpart1, jppart
  real(lk), dimension(:) :: wzpart, marker, meshn, mesht, nonuni, pmark
  integer, dimension(:) :: marknum
  integer mp, mpmax, mpcore, mpsol
  real(lk) apart, qpart

  ! local variables
  integer :: i, j, ij, ij0, ij1, ii, m, mtest, ierror, icount, mmin, mmax, nzsp, nz, np,&
             ipsi(1), iz0(1)
  real(lk) :: c0=2.515517, c1=0.802853, c2=0.010328, d1=1.432788, d2=0.189269, d3=0.001308
  real(lk) :: w_initial, rhop, b, umax, pdum, tdum, rmin1, rmax1, rmin2, rmax2, zzmin, zzmax,&
              rdum, zdum, dpr, dpz, del_r, del_z, rdum0, zdum0, vth, anonuni, rg, psi_cy0,&
              psi_cy1, deltapsi_cy, psidum, delta_V, sdum(mpsi), zdum2, dz2, psifrc, del_p, &
              jaco_frc, dr1, dr2, dr3, dr11, dr22, dr33, r0tmp, r1tmp, r2tmp, psidumup, &
              psidumlow, pmin, pmax, psidumtmp, dpsidrtmp, spdR_inv, spdZ_inv, &
              dpsi, dp1, del_s, ssmin, ssmax, ldum, wp1, wp0, wt00, wt10, wt01, wt11, &
              jj0, dz0, ndum, del_psidum, psi_inner_boundary, psi_outer_boundary, &
              vthp, vth_norm, proton_mass, proton_charge, norm_energy, tp, cmratio, &
              br, bz, bt, vecb(3), vec1(3), vec2(3), vecv(3), vecrho(3),vperp,vpara,omegac,tp_x,&
              tp_R
  real(lk), dimension(:), allocatable:: alpha
interface
    subroutine locate_kernel(zpart,wzpart,wppart,wtpart0,wtpart1,jtpart0,&
        jtpart1,ngyro,mp,jppart,qpart,apart,rhogyropart,fullk)
        use precision
        implicit none
    ! dummy parameters
        real(lk), dimension(:,:) :: zpart, wppart, wtpart0, wtpart1,rhogyropart
        integer, dimension(:,:) :: jtpart0, jtpart1, jppart
        real(lk), dimension(:) :: wzpart
        integer :: mp,ngyro
        logical :: fullk
        real(lk) :: qpart,apart
    end subroutine locate_kernel
end interface


  allocate(alpha(mpmax))

  zpart=0.0_lk
  zpart0=0.0_lk
  wzpart=0.0_lk
  wppart=0.0_lk
  jtpart0=0.0_lk
  jtpart1=0.0_lk
  wtpart0=0.0_lk
  wtpart1=0.0_lk
  jppart=0
  if(irun>0) return
  if(sfcore==1)then
! number of marker per grid, Jacobian
    do i=1, mpsi_core
      do j=1, mpara_core(i)-1
        ij=igrid_core(i)+j
        marker(ij)=meshn(ij)*volume(ij)
        nonuni(ij)=1.0_lk
      enddo
    enddo

    do j=1,mpara_core(1)-1
      ij=igrid_core(1)+j
      nonuni(ij)=0.5_lk*nonuni(ij)
    enddo
    do j=1,mpara_core(mpsi_core)-1
      ij=igrid_core(mpsi_core)+j
      nonuni(ij)=0.5_lk*nonuni(ij)
    enddo
  endif
  if(sfsol==1)then
    do i=1,mpsi_sol
      del_p=0.5_lk*(psimesh_sol(min(mpsi_sol,i+1))-psimesh_sol(i))+0.5_lk*(psimesh_sol(i)-psimesh_sol(max(1,i-1)))
      do j=1,mpara_sol(i)
        ij=mgrid_core+igrid_sol(i)+j
        ij0=igrid_sol(i)+j
        if(j==1)then
          del_s=0.5_lk*(Ssol(ij0+1)-Ssol(ij0))
        elseif(j==mpara_sol(i))then
          del_s=0.5_lk*(Ssol(ij0)-Ssol(ij0-1))
        else
          del_s=0.5_lk*(Ssol(ij0+1)-Ssol(ij0))+0.5_lk*(Ssol(ij0)-Ssol(ij0-1))
        endif
        marker(ij)=meshn(ij)*volume(ij)
        nonuni(ij)=1.0_lk
      enddo
    enddo

    do i=1,mpsi_sol
      nonuni(mgrid_core+igrid_sol(i)+1)=0.5_lk*nonuni(mgrid_core+igrid_sol(i)+1)
      nonuni(mgrid_core+igrid_sol(i)+mpara_sol(i))=0.5_lk*nonuni(mgrid_core+igrid_sol(i)+mpara_sol(i))
    enddo
    do j=1, mpara_sol(1)
      ij=mgrid_core+igrid_sol(1)+j
      nonuni(ij)=0.5_lk*nonuni(ij)
    enddo
    do j=1,mpara_sol(mpsi_sol)
      ij=mgrid_core+igrid_sol(mpsi_sol)+j
      nonuni(ij)=0.5_lk*nonuni(ij)
    enddo
  endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  pdum=sum(nonuni(1:mgrid))
  if(sfcore==1)then
!$omp parallel do private(i,j,ij)
    do i=1,mpsi_core
      do j=1,mpara_core(i)-1
        ij=igrid_core(i)+j
        marknum(ij)=max(1,int(nonuni(ij)*real(mp)/pdum)) !# of marker per cell
      enddo
    enddo
  endif
  if(sfsol==1)then
!$omp parallel do private(i,j,ij)
    do i=1,mpsi_sol
      do j=1,mpara_sol(i)
        ij=mgrid_core+igrid_sol(i)+j
        marknum(ij)=max(1, int(nonuni(ij)*real(mp)/pdum)) !# of marker per cell
      enddo
    enddo
  endif

  pmark=0.0_lk
  if(sfcore==1)then
    do i=1,mpsi_core
      do j=1,mpara_core(i)-1
        ij=igrid_core(i)+j
        pmark(i)=pmark(i)+real(marknum(ij))
      enddo
    enddo
  endif
  if(sfsol==1)then
    do i=1, mpsi_sol
      do j=1,mpara_sol(i)
        ij=mgrid_core+igrid_sol(i)+j
        pmark(mpsi_core+i)=pmark(mpsi_core+i)+real(marknum(ij))
      enddo
    enddo
  endif
  pmark=pmark*real(numberpe)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!marker number which is used for nonmalization, proportional to physcial density
  pdum=0.0_lk
  if(sfcore==1)then
    do i=1,mpsi_core
      do j=1,mpara_core(i)-1
        ij=igrid_core(i)+j
        pdum=pdum+marker(ij)
      enddo
    enddo
    if(mype==0) write(gtcout,*) "pdum_core=",pdum
  endif

  if(sfsol==1)then
    pdum=pdum+sum(marker(mgrid_core+1:mgrid))
    if(mype==0) write(gtcout,*) "pdum_total", pdum
  endif

  if(sfcore==1)then
!$omp parallel do private(i,j,ij)
  do i=1,mpsi_core
    do j=1,mpara_core(i)-1
       ij=igrid_core(i)+j
       marker(ij)=(marker(ij)*real(mp)/pdum) !# of marker per cell
    enddo
  enddo
  endif

  if(sfsol==1)then
!$omp parallel do private(i,j,ij)
  do i=1,mpsi_sol
    do j=1,mpara_sol(i)
      ij=mgrid_core+igrid_sol(i)+j
      marker(ij)=(marker(ij)*real(mp)/pdum)
    enddo
  enddo
  endif

  if(sfcore==1)then
!$omp parallel do private(i,j,ij)
    do i=1,mpsi_core
      do j=1,mpara_core(i) - 1
        ij=igrid_core(i)+j
        marker(ij)=real(npartdom)*marker(ij)
      enddo
    enddo
  endif

  if(sfsol==1)then
!$omp parallel do private(i,j,ij)
    do i=1,mpsi_sol
      do j=1,mpara_sol(i)
        ij=mgrid_core+igrid_sol(i)+j
        marker(ij)=real(npartdom)*marker(ij)
      enddo
    enddo
  endif

  if(sfcore==1)then
!$omp parallel do private(i,j,ij)
    do i=1,mpsi_core
      do j=1,mpara_core(i) - 1
        ij=igrid_core(i)+j
        nonuni(ij)=marker(ij)/(real(marknum(ij))*real(npartdom))
      enddo
    enddo
    do i=1,mpsi_core
      nonuni(igrid_core(i)+mpara_core(i))=nonuni(igrid_core(i)+1)
      marker(igrid_core(i)+mpara_core(i))=marker(igrid_core(i)+1)
    enddo
    ndum=0.0_lk
    do i=1,mpsi_core
      do j=1,mpara_core(i) - 1   ! Ghost grid of core region should be excluded
        ij=igrid_core(i)+j
        ndum=ndum+nonuni(ij)
      enddo
    enddo
    if(mype==0)write(gtcout,*)"average nonuni_core=",ndum/(real(mgrid_core-mpsi_core))    
  endif
  if(sfsol==1)then
!$omp parallel do private(i,j,ij)
    do i=1,mpsi_sol
      do j=1,mpara_sol(i)
         ij=mgrid_core+igrid_sol(i)+j
         nonuni(ij)=marker(ij)/(real(marknum(ij))*real(npartdom))
      enddo
    enddo
    if(mype==0) write(gtcout,*) "average nonuni_sol", sum(nonuni(mgrid_core+1:mgrid))/(real(mgrid_sol))
  endif
  
  if(couplefrc==1)then
    do i=1,mpara_share
      marker(igrid_core(mpsi_core)+mpara_core_share_index(i))=marker(igrid_core(mpsi_core)+mpara_core_share_index(i))+marker(mgrid_core+igrid_sol(1)+mpara_sol_share_index(i))
      marker(mgrid_core+igrid_sol(1)+mpara_sol_share_index(i))=marker(igrid_core(mpsi_core)+mpara_core_share_index(i))
    enddo
    marker(igrid_core(mpsi_core)+mpara_core(mpsi_core))=marker(igrid_core(mpsi_core)+1)
    if(mype==0) write(gtcout,*) "average nonuni", sum(nonuni)/(real(mgrid-mpsi_core))
  endif
!$omp parallel do private(i)
  do i=1,mgrid
    marker(i)=1.0_lk/marker(i) ! to avoid divide operation
  enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  mp=sum(marknum(1:mgrid))
  mpcore=sum(marknum(1:mgrid_core))
  mpsol=sum(marknum(mgrid_core+1:mgrid))
  w_initial=1.0e-3_lk
  umax=4.0
! random # uniformly distributed between 0 and 1
  call random_number(zpart(1:3,1:mp))
  call random_number(zpart(4,1:mp))
  call random_number(zpart(5,1:mp))
  call random_number(zpart(6,1:mp))
  call random_number(zpart(8,1:mp))
  call random_number(alpha)
  zpart(5,:)=zpart(5,:)-0.5_lk

  if(sfcore==1)then
    mmin=0
    mmax=0
    do i=1,mpsi_core
      do j=1,mpara_core(i)-1
        ij=igrid_core(i)+j
        del_s=deltas_core1d(i)/totals_core(i)
        ssmin=modulo(Score(ij)-0.5_lk*del_s, 1.0_lk)
        ssmax=Score(ij)+0.5_lk*del_s
        if(j/=1)del_s=ssmax-ssmin
        if(j==1)del_s=1.0_lk-ssmin+ssmax
        del_p=0.5_lk*(psimesh_core(min(mpsi_core,i+1))-psimesh_core(i))+0.5_lk*(psimesh_core(i)-psimesh_core(max(1,i-1)))
!! range of particle #
        mmin=mmax+1
        mmax=mmin-1+marknum(ij)
! load ion in radial direction: linear in psi within a grid cell
!!$omp parallel do private(m,zdum,ldum,rdum,psidum,tdum)
        do m=mmin,mmax
          zdum=zpart(3,m)
          ldum=modulo(ssmin+zdum*del_s,1.0_lk)
          rdum=zpart(1,m)
          psidum=psimesh_core(i)-0.5_lk*(psimesh_core(i)-psimesh_core(max(1,i-1)))+rdum*del_p
          zpart(1,m)=spline2d(psidum-msp_psi_core(1),ldum,0,lspcore,lstcore,del_psi_core,del_s_core,rspcore_psi_s,0,0)
          zpart(3,m)=spline2d(psidum-msp_psi_core(1),ldum,0,lspcore,lstcore,del_psi_core,del_s_core,zspcore_psi_s,0,0)
          if(zpart(1,m)<0) zpart(1,m)=-zpart(1,m)
          zpart(5,m)=w_initial*zpart(5,m)
          zpart0(7,m)=nonuni(ij)
          zpart0(1,m)=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
          tdum=geotheta(zpart(1,m),zpart(3,m))
          zpart(3,m)=spline2d(psidum-msp_psi_core(1),tdum,0,lspcore,lstcore,del_psi_core,del_theta_core,sspcore_psi_theta,0,0)
        enddo
      enddo
    enddo
  endif
  if(sfsol==1)then
    mmin=0
    mmax=mpcore
    do i=1,mpsi_sol
      do j=1,mpara_sol(i)
        ij=mgrid_core+igrid_sol(i)+j
        ij0=igrid_sol(i)+j

        ssmin=Ssol(ij0)-0.5_lk*(Ssol(ij0)-Ssol(max(ij0-1, igrid_sol(i)+1)))
        ssmax=Ssol(ij0)+0.5_lk*(Ssol(min(ij0+1, igrid_sol(i))+mpara_sol(i)))-Ssol(ij0)
       
        del_s=ssmax-ssmin
        del_p=0.5_lk*(psimesh_sol(min(mpsi_sol,i+1))-psimesh_sol(i))+0.5_lk*(psimesh_sol(i)-psimesh_sol(max(1,i-1)))
        mmin=mmax+1
        mmax=mmin-1+marknum(ij)
! load ion in radial direction: linear in psi within a grid cell
!!$omp parallel do private(m,rdum,ldum,zdum,psidum)
        do m=mmin,mmax
          zdum=zpart(3,m)
          ldum=ssmin+zdum*del_s
          rdum=zpart(1,m)
          psidum=psimesh_sol(i)-0.5_lk*(psimesh_sol(i)-psimesh_sol(max(1,i-1)))+rdum*del_p
          zpart(1,m)=spline2d(psidum-msp_psi_sol(1),ldum,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
          zpart(3,m)=spline2d(psidum-msp_psi_sol(1),ldum,0,lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)
          if(zpart(1,m)<0) zpart(1,m)=-zpart(1,m)

          zpart(5,m)=w_initial*zpart(5,m)
          zpart0(7,m)=nonuni(ij)
        enddo
      enddo
    enddo
  endif

! load ion uniformly in zeta
!$omp parallel do private(m,tdum)
    do m=1,mp
        tdum=zpart(2,m)*deltaz
        zpart(2,m)=zeta0+tdum
        zpart(10,m)=zpart(2,m)   ! zion(9,m) records actual zeta value for partial torus
    enddo

    if(mype==0)write(gtcout,*)"deltaz =",deltaz

!$omp parallel do private(m,zdum)
    do m=1,mp
        zdum=zpart(4,m)
        zpart(4,m)=zpart(4,m)-0.5
        zpart0(4,m)=sign(1.0_lk,zpart(4,m))
        zpart(4,m)=sqrt(max(1.0e-20_lk,log(1.0_lk/max(1.0e-20_lk,zpart(4,m)**2))))
        zpart(4,m)=zpart(4,m)-(c0+c1*zpart(4,m)+c2*zpart(4,m)**2)/&
          (1.0+d1*zpart(4,m)+d2*zpart(4,m)**2+d3*zpart(4,m)**3)
        if(zpart(4,m)>umax)zpart(4,m)=zdum
        zpart(4,m)=zpart0(4,m)*zpart(4,m)
        zpart(6,m)=max(1.0e-20_lk,min(umax*umax,-log(max(1.0e-20_lk,zpart(6,m)))))
    enddo

  if(sfcore==1)then
!$omp parallel do private(m,rdum,zdum,b,tp,psidum,ipsi,dpsi,ii,dp1,rhop,vth)
      do m=1,mpcore
         rdum=zpart(1,m)
         zdum=zpart(3,m)
         b = spline2d(rdum,zdum+zmax,0,lsr,lsz,spdR,spdZ,B0mag_Eq,0,0)
         tp = spline2d(rdum,zdum+zmax,0,lsr,lsz,spdR,spdZ,t_Eq,0,0)
         psidum = spline2d(rdum,zdum+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)


         if(.true.)then
           ipsi=minloc(abs(psimesh_core-psidum))
           dpsi=psidum-psimesh_core(ipsi(1))
           if(dpsi>0.0)then
             ii=min(ipsi(1),mpsi_core-1)
             psidum=min(psidum,psimesh_core(mpsi_core))
           else
             ii=max(ipsi(1)-1,1)
             psidum=max(psidum,psimesh_core(1))
           endif
           dp1=(psidum-psimesh_core(ii))/(psimesh_core(ii+1)-psimesh_core(ii)) !weight for outer flux surface
           tp=dp1*meshti1d(ii+1)+(1-dp1)*meshti1d(ii)
         endif

         rhop=sqrt(tp)*sqrt(apart)
         vth=sqrt(tp)/sqrt(apart)
         zpart(4,m)=vth*zpart(4,m)
         zpart(6,m)=vth*sqrt(apart*zpart(6,m)/b)

         zpart(8,m)=0.5*apart*zpart(4,m)*zpart(4,m)+zpart(6,m)*zpart(6,m)*b
      enddo
    endif
  if(sfsol==1)then
!$omp parallel do private(m,rdum,zdum,b,tp,psidum,ipsi,dpsi,ii,dp1,rhop,vth)
      do m=mpcore+1,mp
         rdum=zpart(1,m)
         zdum=zpart(3,m)
         b = spline2d(rdum,zdum+zmax,0,lsr,lsz,spdR,spdZ,B0mag_Eq,0,0)
         tp = spline2d(rdum,zdum+zmax,0,lsr,lsz,spdR,spdZ,t_Eq,0,0)
         psidum = spline2d(rdum,zdum+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

         if(.true.)then
         ipsi=minloc(abs(psimesh_sol-psidum))
         dpsi=psidum-psimesh_sol(ipsi(1))
         if(dpsi>0.0)then
         ii=min(ipsi(1),mpsi_sol-1)
         psidum=min(psidum,psimesh_sol(mpsi_sol))
         else
         ii=max(ipsi(1)-1,1)
         psidum=max(psidum,psimesh_sol(1))
         endif
         dp1=(psidum-psimesh_sol(ii))/(psimesh_sol(ii+1)-psimesh_sol(ii)) !weight for outer flux surface
         tp=dp1*meshti1d(ii+1+mpsi_core)+(1-dp1)*meshti1d(ii+mpsi_core)
         endif
         rhop=sqrt(tp)*sqrt(apart)
         vth=sqrt(tp)/sqrt(apart)
         zpart(4,m)=vth*zpart(4,m)
         zpart(6,m)=vth*sqrt(apart*zpart(6,m)/b)
         zpart(8,m)=0.5*apart*zpart(4,m)*zpart(4,m)+zpart(6,m)*zpart(6,m)*b
      enddo
      if(mype==0)write(gtcout,*)"vth",sqrt(meshti1d(1+mpsi_core))/sqrt(apart)
    endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call locate_kernel(zpart, wzpart, wppart, wtpart0, wtpart1, jtpart0, jtpart1, 1,&
      mp, jppart, qpart, apart,rhogyropart, .true.)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    if(.true.)then
!$omp parallel do private(m,wp1,wp0,wt00,wt10,wt01,wt11,ij)
    do m=1,mp
        zpart0(7,m)=0.0_lk

        wp1=wppart(1,m)       !outer flux surface
        wp0=1.0_lk-wp1              !inner flux surface

        wt10=wp0*wtpart0(1,m) !upper poloidal grid on inner flux surface
        wt00=wp0-wt10            !lower poloidal grid on inner flux surface

        wt11=wp1*wtpart1(1,m) !upper poloidal grid on outer flux surface
        wt01=wp1-wt11            !lower poloidal grid on outer flux surface  

        ij=jtpart0(1,m)       !lower poloidal grid on inner flux surface
        zpart0(7,m)=zpart0(7,m)+wt00*nonuni(ij)

        ij=ij+1                  !upper poloidal grid on inner flux surface
        zpart0(7,m)=zpart0(7,m)+wt10*nonuni(ij)

        ij=jtpart1(1,m)       !lower poloidal grid on outer flux surface
        zpart0(7,m)=zpart0(7,m)+wt01*nonuni(ij)

        ij=ij+1                  !upper poloidal grid on outer flux surface
        zpart0(7,m)=zpart0(7,m)+wt11*nonuni(ij)
    enddo
    endif

    alpha=(alpha-0.5_lk)*pi
    cmratio=qpart/apart
    do m=1,mp
        zpart(10,m)=zpart(9,m)
        zpart(9,m)=zpart(8,m)
        !b=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,B0mag_Eq,0,0)
        br=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,B0R_Eq,0,0)
        bz=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,B0Z_Eq,0,0)
        !b=max(b,sqrt(br**2+bz**2))
        !bt=sqrt(b**2-br**2-bz**2)
        bt=0.0_lk
        b=sqrt(bz**2+br**2)
        vecb=(/br/b,bt/b,bz/b/)
        vec1=(/bz/b*br/sqrt(br**2+bt**2), bz/b*bt/sqrt(br**2+bt**2), -sqrt(1-bz**2/b**2)/)
        call cross_product3(vecb, vec1, vec2)
        vperp=zpart(6,m)*sqrt(2*b/apart)
        vpara=zpart(4,m)
        vecv=vperp*cos(alpha(m))*vec2+vperp*sin(alpha(m))*vec1+vpara*vecb
        omegac=cmratio*b
        vecrho=(vperp*cos(alpha(m))*vec1-vperp*sin(alpha(m))*vec2)/omegac
        zpart0(11,m)=zpart(1,m)-vecrho(1)
        zpart0(13,m)=zpart(3,m)-vecrho(3)
        zpart0(12,m)=zpart(2,m)
        zpart0(13,m)=modulo((zpart0(13,m)+zmax),(zmax*2.0_lk))-zmax
        zpart(10,m)=zpart0(12,m)
        !initial particle position at -dt/2
        zpart(11,m)=zpart0(11,m)-vecv(1)*tstep*0.5
        zpart(12,m)=zpart0(12,m)-vecv(2)/zpart0(11,m)*tstep*0.5
        zpart(13,m)=zpart0(13,m)-vecv(3)*tstep*0.5
        zpart(13,m)=modulo((zpart(13,m)+zmax),(zmax*2.0_lk))-zmax
        zpart(4,m)=vecv(1)*cos(zpart0(12,m)-zpart(12,m))-vecv(2)*sin(zpart0(12,m)-zpart(12,m))
        zpart(6,m)=vecv(1)*sin(zpart0(12,m)-zpart(12,m))+vecv(2)*cos(zpart0(12,m)-zpart(12,m))
        zpart(6,m)=zpart(6,m)/zpart(11,m)
        zpart(8,m)=vecv(3)
    enddo

    if(frcfullf==1)then
!$omp parallel do private(m)
        do m=1,mp
            zpart(5,m)=zpart0(7,m)
        enddo
    endif


    del_psidum=0.1_lk*(psi1_core-psi0_core)
    psi_inner_boundary = psi0_core + del_psidum
    psi_outer_boundary = psi1_core - del_psidum
!$omp parallel do private(m,rdum,zdum,psidum,tdum,psidumtmp)
    do m=1,mp
        rdum=zpart(1,m)
        zdum=zpart(3,m)
        psidum = spline2d(rdum,zdum+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

        if(sfcore==1 .and. sfsol==0)then
!DONT CHANGE ANYTHING(repeat part should be in ifelse control)
            if(psidum<psi0_core)then
                psidum=max(psi_inner_boundary,min(psi_outer_boundary,psi_inner_boundary+(psi0_core-psidum)))
                tdum=geotheta(zpart(1,m),zpart(3,m))
                zpart(1,m)=spline2d(psidum-msp_psi_core(1),tdum,0,lspcore,lstcore,del_psi_core,del_theta_core,rspcore_psi_theta,0,0)
                zpart(3,m)=spline2d(psidum-msp_psi_core(1),tdum,0,lspcore,lstcore,del_psi_core,del_theta_core,zspcore_psi_theta,0,0)
                psidumtmp = spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

                if(psidumtmp < psi0_core)then
           !     if((isnan(zion(1,m)) .or. isnan(zion(3,m))))then
                    print *,"istep=",istep,"irk=",irk
                    print *,"mype",mype,"m=",m,"rdum=",zpart(1,m),"zdum=",zpart(3,m)
                    print *,"psidum=",psidum,"psidumtmp=",psidumtmp
                    print *,"tdum=",tdum,"rdum=",rdum,"zdum=",zdum
                    stop
                endif
            elseif(psidum>psi1_core)then
                psidum=min(psi_outer_boundary,max(psi_inner_boundary,psi_outer_boundary-(psidum-psi1_core)))
                tdum=geotheta(zpart(1,m),zpart(3,m))
                zpart(1,m)=spline2d(psidum-msp_psi_core(1),tdum,0,lspcore,lstcore,del_psi_core,del_theta_core,rspcore_psi_theta,0,0)
                zpart(3,m)=spline2d(psidum-msp_psi_core(1),tdum,0,lspcore,lstcore,del_psi_core,del_theta_core,zspcore_psi_theta,0,0)
                psidumtmp = spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

                if(psidumtmp > psi1_core)then
           !     if((isnan(zion(1,m)) .or. isnan(zion(3,m))))then
                    print *,"istep=",istep,"irk=",irk
                    print *,"mype",mype,"m=",m,"zion(1,m)=",zpart(1,m),"zion(3,m)=",zpart(3,m)
                    print *,"psidum=",psidum,"psidumtmp=",psidumtmp
                    print *,"tdum=",tdum,"rdum=",rdum,"zdum=",zdum
                    stop
                endif
            endif
        endif
    enddo
end subroutine loadfk
