! This file is part of GTC version 4.4
! GTC version 4.4 is released under the 3-Clause BSD license:

! Copyright (c) 2002,2010,2016, GTC Team (team leader: Zhihong Lin, zhihongl@uci.edu)
! All rights reserved.

! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:

! 1. Redistributions of source code must retain the above copyright notice,
!    this list of conditions and the following disclaimer.

! 2. Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation
!    and/or other materials provided with the distribution.

! 3. Neither the name of the GTC Team nor the names of its contributors may be
!    used to endorse or promote products derived from this software without
!    specific prior written permission.
! ==============================================================================

!> \file
!> Load or restart all particle species
subroutine load
  use system_env
  use precision, only: lk, mpi_Rsize
  use global_parameters
  use particle_array
  use field_array
  use equilibrium
  use spline_function, only: dxdp, dxdt, dzdp, dzdt, spx, construct_spline1d, invert_spline0
  use antenna, only: particle_initialization
  implicit none

  integer i,j,ij
  real(lk) dmarki(0:1,mgrid),pmarkit(0:mpsi),w_initial_ion,w_initial_fast,w_initial_electron,w_initial_faste
  !$acc declare create(dmarki)

  interface
    subroutine gkLoad(species_name,w_initial)
      use precision
      implicit none
      character(*),intent(in) :: species_name
      real(lk) w_initial
    end subroutine
  end interface

!iload=0: ideal MHD with nonuniform pressure profile (eload=3) or uniform profile (eload=1)
!iload=1: Marker uniform temperature & density: n0,te0,ti0 (two-scale expansion)
!iload=2: Marker non-uniform temperature & uniform density: n0,te(r),ti(r), marker weight corrected
!iload=3: Physical PDF for marker with non-uniform temperature & density: n(r),te(r),ti(r)
!iload>99: non-perturbative (full-f) simulation

! initialize particle position, velocity and weight: thermal ions, fast ions, electrons and fast electrons
  w_initial_ion=1.0e-5
  w_initial_fast=1.0e-5
  w_initial_electron=1.0e-5
  w_initial_faste=1.0e-5

! If antenna on, option to initialize to zero
  if(antenna_on>0)then
    if(particle_initialization==0)then
        w_initial_ion=0.0
        w_initial_fast=0.0
        w_initial_electron=0.0
        w_initial_faste=0.0
    endif
  endif

  call gkLoad("thermal-ion",w_initial_ion)
  if(fload>0) call gkLoad("fast-ion",w_initial_fast)
  if(nhybrid>0)call gkLoad("thermal-electron",w_initial_electron)
  if(feload>0) call gkLoad("fast-electron",w_initial_faste)
  !> @todo
  !> It might be more straight forward to check for restart first, and skip all
  !> the initialization if it is a restart run. L.Shi 03/31/2017
  if(irun>0)call restart_load

! Tag particle for diagnostics
  if(track_particles==1 .and. irun==0)CALL TAG
  call MPI_BARRIER(MPI_COMM_WORLD,ierror)

! For non-perturbative simulation:
! step 1: irun=0, nonlinear=0 for Monte-Carlo run to obtain PDF as constants of motion
! step 2: irun=1 Monte-Carlo run to calculate time-averaged # of marker particles per cell
! step 3: nonlinear=1: turbulence run by using time-averaged # of marker particles per cell
! irun>1: continue turbulence run without changing markeri and pmarki

  if(irun==1 .and. nonlinear==1 .and. iload>99)then
! smoothing # of marker particles per cell
     dmarki(1,:)=markerit
     call smooth_cpu(dmarki)
     markeri=1.0/dmarki(1,:)

! # of marker particles per annulus
     pmarkit=0.0
!$omp parallel do private(i,j,ij)
     do i=0,mpsi
        do j=1,mtheta(i)
           ij=igrid(i)+j
           pmarkit(i)=pmarkit(i)+dmarki(1,ij)
        enddo
     enddo

! global sum for # of marker particles per cell
     call MPI_ALLREDUCE(pmarkit,pmarki,mpsi+1,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)

! reset field and particle weight for perturbative (delta-f) simulation iload >1
!     zion(5,:)=0.0
!     phi=0.0
!     zonali=0.0
!     call random_number(zion(5,1:mi))
!     zion(5,1:mi)=0.001*zion(5,1:mi)

  endif

  markerit=0.0 ! re-set time averaging

end subroutine load

subroutine restart_load
  use system_env
  use utility, only: check_allocation, mstat
  use global_parameters
  use particle_tracking,only:ptrackedi,ptrackede
  use particle_array,only:mimax,feload,mfemax,&
      zion,zelectron,zfast,zfaste,mi,me,mf,mfe,fload,&
      markeri,markerf,markere,markerfe,zonali,zonalci,&
      zonalf,zonalcf,zonale,zonalce,zonalfe,zonalcfe,&
      pmarki,pmarke,pmarkf,pmarkfe
  implicit none

  integer iorestart
  character(len=9)cdum1
  if(track_particles==1) then
    allocate(ptrackedi(nparami,max(mimax,1)), STAT=mstat)
    call check_allocation(mstat, "ptrackedi in load.F90")
  endif

#ifdef _FRC
  if(track_particles==1 .and. feload>0) then
    allocate(ptrackede(nparamfe,max(mfemax,1)))
    call check_allocation(mstat, "ptrackede in load.F90")
  endif
#endif
  if(mype==0)then
!!get restart directory info
     iorestart=345
     open(iorestart,file="FileExit.out",status="old")
     read(iorestart,"(A9,i1)")cdum1,FileExit
     read(iorestart,"(A9,i5)")cdum1,irest
     close(iorestart)
     irest=irest-1
  endif
  call mpi_bcast(irest,1,mpi_integer,0,mpi_comm_world,ierror)
  call MPI_BARRIER(MPI_COMM_WORLD,ierror)
  call restart_io("read")
  !$acc update device(zion(:,1:mi),markeri,zonali,zonalci,pmarki)
  if(fload>0)then
    !$acc update device(zfast(:,1:mf),markerf,zonalf,zonalcf,pmarkf)
  endif
  if(nhybrid>0)then
    !$acc update device(zelectron(:,1:me),markere,zonale,zonalce,pmarke)
  endif
  if(feload>0)then
    !$acc update device(zfaste(:,1:mfe),markerfe,zonalfe,zonalcfe,pmarkfe)
  endif
  return

end subroutine restart_load

subroutine gkLoad(species_name,w_initial)
  use system_env
  use utility, only: check_allocation, mstat
  use global_parameters,only: mype,gtcout,nparami,nparame,nparamf,nhybrid,&
                              ncyclee,nparamfe,ncyclefe,irun,mgrid,mpsi,deltabpara
  use particle_array
  use field_array,only: nmodes,nmode
  use equilibrium,only: nipp,tipp,nepp,tepp,nfpp,tfpp,nfepp,tfepp
  implicit none

  integer ngyroi0, ngyrof0
  real(lk) w_initial,mpartfactor
  character(*),intent(in) :: species_name

  interface
    subroutine gkLoadParticle(meshn,tppp,nppp,zpart,zpart0,wppart,&
        wtpart0,wtpart1,density,flow,data1d,jtpart0,jtpart1,wzpart,&
        marker,markert,pmark,pflux,rdtem,rdden,zonal,zonalc,mesht,qpart,apart,&
        pload,nparam,ngyro,mp,mpmax,w_initial,ncyclep,zpart1,&
        trap,trapfracn,trapfrac,hybrid,pressurepara,pressureperp,phit,dnt,&
        phisave,dnsave)

      use precision
      use spline_function, only: construct_spline1d, invert_spline0
      implicit none
      integer pload,nparam,ngyro,mp,mpmax
      integer,optional :: ncyclep,trap,hybrid
      integer,dimension(:,:) :: jtpart0,jtpart1
      real(lk) qpart,apart,w_initial
      real(lk),optional :: trapfracn
      real(lk),dimension(0:) :: meshn,mesht,pmark,pflux,rdtem,rdden,zonal,zonalc
      real(lk),dimension(:) :: wzpart,marker,markert
      real(lk),dimension(:),allocatable,optional :: trapfrac
      real(lk),dimension(:,:) :: tppp,nppp
      real(lk),dimension(:,:) :: zpart,zpart0,wppart,wtpart0,wtpart1
      real(lk),dimension(0:,:) :: density,flow,data1d
      real(lk),dimension(:,:),optional :: zpart1
      real(lk),dimension(0:,:),optional :: pressurepara,&
        pressureperp,phit,dnt
      real(lk),dimension(0:,:,:),optional :: phisave,dnsave
    end subroutine gkLoadParticle
  end interface

! 3D field has strongly changing Jacobian along zeta
! We enlarge the particle buffer by a factor of 3 to guarantee the particle
! shift can be successfully done.
#ifdef _TOROIDAL3D
  mpartfactor=3
#else
  mpartfactor=1
#endif
ngyroi0=ngyroi
if(deltabpara==1)ngyroi0=12
  if(species_name=='thermal-ion')then
    ! allocate particle memory
    mimax=mi+100*ceiling(sqrt(real(mi))) !particle array upper bound
    mimax=mimax*mpartfactor
    allocate(zion(nparami,mimax),zion0(nparami,mimax),wzion(mimax),&
      wpion(ngyroi0,mimax),jtion0(ngyroi0,mimax),jtion1(ngyroi0,mimax),&
      wtion0(ngyroi0,mimax),wtion1(ngyroi0,mimax),STAT=mstat)
    call check_allocation(mstat, "zion and others in load.F90")
    allocate(pmarki(0:mpsi),markeri(mgrid),densityi(0:1,mgrid),&
      flowi(0:1,mgrid),pfluxi(0:mpsi),rdtemi(0:mpsi),rddeni(0:mpsi),zonali(0:mpsi),&
      zonalci(0:mpsi),data1di(0:mpsi,mpdata1d),markerit(mgrid),STAT=mstat)
    call check_allocation(mstat, "pmarki and others in load.F90")
    if(deltabpara==1)then
      allocate(pressureiperp(0:1,mgrid),STAT=mstat)
      pressureiperp=0.0_lk
      call check_allocation(mstat, "pressureiperp in load.F90")
    endif

    call gkLoadParticle(meshni,tipp,nipp,zion,zion0,wpion,wtion0,wtion1,&
      densityi,flowi,data1di,jtion0,jtion1,wzion,markeri,markerit,pmarki,&
      pfluxi,rdtemi,rddeni,zonali,zonalci,meshti,qion,aion,iload,&
      nparami,ngyroi,mi,mimax,w_initial)

    if(iload==9)call gk2fk

    if(mype==0 .and. irun==0)then
      write(gtcout,*)'mi=',mi
      flush(gtcout)
    endif

    if(irun==0)then
      !$acc update device(zion,zion0,markeri,zonali,zonalci)
    endif

  elseif(species_name=='thermal-electron')then
    ! allocate particle memory
    memax=me+100*ceiling(sqrt(real(me))) !particle array upper bound
    memax=memax*mpartfactor
    allocate(zelectron(nparame,memax),zelectron0(nparame,memax),&
      wzelectron(memax),wpelectron(1,memax),jtelectron0(1,memax),&
      jtelectron1(1,memax),wtelectron0(1,memax),wtelectron1(1,memax),&
      STAT=mstat)
    call check_allocation(mstat, "zelectron and others in load.F90")
    allocate(pmarke(0:mpsi),markere(mgrid),densitye(0:1,mgrid),&
      flowe(0:1,mgrid),pfluxe(0:mpsi),rdteme(0:mpsi),rddene(0:mpsi),zonale(0:mpsi),&
      zonalce(0:mpsi),data1de(0:mpsi,mpdata1d),markeret(mgrid),STAT=mstat)
    call check_allocation(mstat, "pmarke and others in load.F90")
    allocate(pressureepara(0:1,mgrid),pressureeperp(0:1,mgrid),&
      phisave(0:1,mgrid,2*nhybrid),phit(0:1,mgrid),&
      dnesave(0:1,mgrid,2*nhybrid),dnet(0:1,mgrid),STAT=mstat)
    call check_allocation(mstat, "pressureepara and others in load.F90")

    if(ncyclee>0) then
      allocate(zelectron1(nparame,memax),STAT=mstat)
      call check_allocation(mstat, "zelectron1 in load.F90")
    endif

    call gkLoadParticle(meshne,tepp,nepp,zelectron,zelectron0,wpelectron,&
      wtelectron0,wtelectron1,densitye,flowe,data1de,jtelectron0,jtelectron1,&
      wzelectron,markere,markeret,pmarke,pfluxe,rdteme,rddene,zonale,zonalce,meshte,&
      qelectron,aelectron,eload,nparame,1,me,memax,&
      w_initial,ncyclep=ncyclee,zpart1=zelectron1,trap=etrap,&
      trapfracn=tfracn,trapfrac=tfrac,hybrid=nhybrid,&
      pressurepara=pressureepara,pressureperp=pressureeperp,phit=phit,&
      dnt=dnet,phisave=phisave,dnsave=dnesave)
    if(mype==0 .and. irun==0)then
       write(gtcout,*)'me=',me,',   trapped fraction',tfracn
       flush(gtcout)
    endif

    if(irun==0)then
       !$acc update device(zelectron,zelectron0,markere,zonale,zonalce,phit,dnet,phisave,dnesave)
    endif

  elseif(species_name=='fast-ion')then
    ngyrof0=ngyrof
    if(deltabpara==1)ngyrof0=12
    ! allocate particle memory
    mfmax=mf+100*ceiling(sqrt(real(mf))) !particle array upper bound
    mfmax=mfmax*mpartfactor
    allocate(zfast(nparamf,mfmax),zfast0(nparamf,mfmax),wzfast(mfmax),&
      wpfast(ngyrof0,mfmax),jtfast0(ngyrof0,mfmax),jtfast1(ngyrof0,mfmax),&
      wtfast0(ngyrof0,mfmax),wtfast1(ngyrof0,mfmax),STAT=mstat)
    call check_allocation(mstat, "zfast and orthers in load.F90")

    allocate(pmarkf(0:mpsi),markerf(mgrid),densityf(0:1,mgrid),&
      flowf(0:1,mgrid),pfluxf(0:mpsi),rdtemf(0:mpsi),rddenf(0:mpsi),zonalf(0:mpsi),&
      zonalcf(0:mpsi),data1df(0:mpsi,mpdata1d),markerft(mgrid),STAT=mstat)
    call check_allocation(mstat, "pmarkf and orthers in load.F90")
    if(deltabpara==1)then
      allocate(pressurefperp(0:1,mgrid),STAT=mstat)
      call check_allocation(mstat, "pressurefperp in load.F90")
    endif

    call gkLoadParticle(meshnf,tfpp,nfpp,zfast,zfast0,wpfast,wtfast0,&
      wtfast1,densityf,flowf,data1df,jtfast0,jtfast1,wzfast,markerf,markerft,&
      pmarkf,pfluxf,rdtemf,rddenf,zonalf,zonalcf,meshtf,qfast,afast,fload,&
      nparamf,ngyrof,mf,mfmax,w_initial)

    if(mype==0 .and. irun==0)then
      write(gtcout,*)'mf=',mf
      flush(gtcout)
    endif

    if(irun==0)then
      !$acc update device(zfast,zfast0,markerf,zonalf,zonalcf)
    endif

  elseif(species_name=='fast-electron')then
    ! allocate particle memory
    mfemax=mfe+100*ceiling(sqrt(real(mfe))) !particle array upper bound
    mfemax=mfemax*mpartfactor
    allocate(zfaste(nparamfe,mfemax),zfaste0(nparamfe,mfemax),wzfaste(mfemax),&
      wpfaste(ngyrofe,mfemax),jtfaste0(ngyrofe,mfemax),&
      jtfaste1(ngyrofe,mfemax),wtfaste0(ngyrofe,mfemax),&
      wtfaste1(ngyrofe,mfemax),STAT=mstat)
    call check_allocation(mstat, "zfaste and orthers in load.F90")
    allocate(pmarkfe(0:mpsi),markerfe(mgrid),densityfe(0:1,mgrid),&
      flowfe(0:1,mgrid),pfluxfe(0:mpsi),rdtemfe(0:mpsi),rddenfe(0:mpsi),zonalfe(0:mpsi),&
      zonalcfe(0:mpsi),data1dfe(0:mpsi,mpdata1d),markerfet(mgrid),STAT=mstat)
    call check_allocation(mstat, "pmarkfe and orthers in load.F90")

    if(ncyclefe>0) then
      allocate(zfaste1(nparamfe,mfemax),STAT=mstat)
      call check_allocation(mstat, "zfaste1 in load.F90")
    endif

    call gkLoadParticle(meshnfe,tfepp,nfepp,zfaste,zfaste0,wpfaste,wtfaste0,&
      wtfaste1,densityfe,flowfe,data1dfe,jtfaste0,jtfaste1,wzfaste,markerfe,&
      markerfet,pmarkfe,pfluxfe,rdtemfe,rddenfe,zonalfe,zonalcfe,meshtfe,qfaste,&
      afaste,feload,nparamfe,ngyrofe,mfe,mfemax,w_initial,&
      ncyclep=ncyclefe,zpart1=zfaste1,trap=fetrap,trapfracn=fetfracn,&
      trapfrac=fetfrac)

    if(mype==0 .and. irun==0)then
      write(gtcout,*)'mfe=',mfe,',   trapped fraction',fetfracn
      flush(gtcout)
    endif

    if(irun==0)then
      !$acc update device(zfaste,zfaste0,markerfe,zonalfe,zonalcfe)
      if(ncyclefe>0)then
        !$acc update device(zfaste1)
      endif
      continue
    endif

  else
    write(*,*)'load.F90: wrong choice'
    call mpi_barrier(mpi_comm_world,ierror)
  endif
end subroutine gkLoad

! Concrete particle loading subroutine
!> @brief Create and initialize particle arrays

!> @param[in] meshn
!> Equilibrium number density of the physical particles on simulation flux
!> surfaces. When pload==1, meshn is flattened before loading markers; when
!> pload==2, meshn is real profile, but marker density is flattened. \n
!> Dummy for: particle_array::meshne, particle_array::meshni, particle_array::meshnfe, particle_array::meshne
!> @see global_parameters::inorm

!> @param[in] tppp
!> equilibrium temperature spline coefficients. \n
!>  Dummy for: equilibrium::tepp, equilibrium::tipp, equilibrium::tfepp,
!> equilibrium::tfpp

!> @param[in] nppp
!> equilibrium density spline coefficients. \n
!>  Dummy for: equilibrium::nepp, equilibrium::nipp, equilibrium::nfepp,
!> equilibrium::nfpp

!> @param[out] zpart
!> marker information array. \n
!> Dummy for particle_array::zelectron,
!> particle_array::zion, particle_array::zfast, particle_array::zfaste. \n
!> The first index indicates the different information for each marker. \n
!> Their meanings are:
!> index | meaning
!> ----- | ----------
!> 1     | psi
!> 2     | theta
!> 3     | zeta
!> 4     | rho_para
!> 5     | weight
!> 6     | sqrt(mu)
!> 7     | f/g
!> 8,9   | particle tracking labels

!> @param[out] zpart0
!> temporary container for initial particle informations before each
!> Runge-Kutta push. See zpart for the format. \n
!> Dummy for particle_array::zelectron0,
!> particle_array::zion0, particle_array::zfast0, particle_array::zfaste0.

!> @param[out] wppart
!> psi weight array used for perturbed field evaluation at particle location.
!> weight for outer flux surface. \n
!> Dummy for particle_array::wpelectron, particle_array::wpion,
!> particle_array:: wpfast, particle_array::wpfaste. \n
!> See also: locateParticle.

!> @param[out] wtpart0
!> theta weight on upper theta grid point on inner flux surface. \n
!> Dummy for particle_array::wtelectron0, particle_array::wtion0,
!> particle_array::wtfast0, particle_array::wtfaste0.

!> @param[out] wtpart1
!> theta weight on upper theta grid point on outer flux surface. \n
!> Dummy for particle_array:: wtelectron1, particle_array::wtion1,
!> particle_array::wtfast1, particle_array::wtfaste1.

!> @param[out] wzpart
!> zeta weight on upper zeta plane. \n
!> Dummy for particle_array::wzelectron, particle_array::wzion,
!> particle_array::wzfast, particle_array::wzfaste

!>
subroutine gkLoadParticle(meshn,tppp,nppp,zpart,zpart0,wppart,wtpart0,&
    wtpart1,density,flow,data1d,jtpart0,jtpart1,wzpart,marker,markert,pmark,&
    pflux,rdtem,rdden,zonal,zonalc,mesht,qpart,apart,pload,nparam,ngyro,mp,mpmax,&
    w_initial,ncyclep,zpart1,trap,trapfracn,trapfrac,hybrid,&
    pressurepara,pressureperp,phit,dnt,phisave,dnsave)
  use system_env
  use utility, only: check_allocation, mstat
  use precision, only: lk, mpi_Rsize
  use global_parameters,only: mgrid,mype,irun,npartdom,numberpe,pi,rg0,gtcout,&
    mpsi,psi0,psi1,ipoincare,inorm,istep,eq_flux,fielddir,pi2,left_pe,right_pe,&
    myrank_toroidal,toroidal_comm,toroidal_domain_location,mtoroidal,numberpe,partd_comm
  use particle_array,only: mpdata1d,sd_v0,sd_vc,sd_l0,sd_widthInv
  use field_array,only: deltap,psimesh,qmesh,mtheta,igrid,bmesh,deltat,deltar,&
    deltaz,zeta0,bmesh,nmodes,nmode,solvermethod,zeta1,qtinv
  use equilibrium,only: lsp,spdpsi_inv,spdpsi,gpsi,cpsi,spdtheta,spdtheta_inv,&
    bsp,rgpsi,lst,ropp,xsp,erpp
  use spline_function, only: spb, spx, spq, sptrot, spgpsi, spcpsi, sprgpsi, invert_spline0, construct_spline1d
  implicit none

  !declaration of the dummy arguments
  !character(*),intent(in) :: w_init_drive
  integer pload,nparam,ngyro,mp,mpmax
  integer,optional :: ncyclep,trap,hybrid
  integer,dimension(:,:) :: jtpart0,jtpart1
  real(lk) qpart,apart
  real(lk),optional :: trapfracn
  real(lk),dimension(0:) :: meshn,mesht,pmark,pflux,rdtem,rdden,zonal,zonalc
  real(lk),dimension(:) :: wzpart,marker,markert
  real(lk),dimension(:),allocatable,optional :: trapfrac
  real(lk),dimension(:,:) :: tppp,nppp
  real(lk),dimension(:,:) :: zpart,zpart0,wppart,wtpart0,wtpart1
  real(lk),dimension(0:,:) :: density,flow,data1d
  real(lk),dimension(:,:),optional :: zpart1
  real(lk),dimension(0:,:),optional :: pressurepara,&
    pressureperp,phit,dnt
  real(lk),dimension(0:,:,:),optional :: phisave,dnsave

  ! Maxwellian cutoff parameter, unit in standard deviation
  real(lk),parameter:: umax=4.0
  ! Maxwellian Inverse CDF expression with rational fraction function
  ! parameters. Reference?
  real(lk),parameter :: c0=2.515517,c1=0.802853,c2=0.010328,d1=1.432788,d2=0.189269,d3=0.001308
  !declaration of the local variables
  integer i,isp,j,jst,ii,ij,m,ir,marknum(0:mpsi),&
    icount,jt,mmin,mmax,mptmp,perturbmarknum(0:mpsi),init
  real(lk) w_initial,energy,rg,delr,pdum,functiong,functioni,dp1,vth,&
    b,bm,q,upara,eperp,cmratio,zdum,tdum,delpsi,delp(mpsi),&
    dtx,dt2,pmin,pmax,volt(3,lst),tvol(3,lst),spdt_inv,dpx,dp2,fullf,&
    xdum,dden(0:mpsi),ddum(0:mpsi),majorr,upara0,&
    markload(0:mpsi),tpref,tp_inv,perturbpmark(0:mpsi),randvarv,&
    randvarf,randvart,sd_ftmp,fmax,zpartsign
#ifdef _TOROIDAL3D
  integer k_sub, interp,izeta,mm,mpi0,mpi1,mpiadd
  real(lk),dimension(:,:),allocatable :: jacobian, delvolume,delvolume_2
  real(lk),dimension(:),allocatable :: vol1, vol2, marker3d !, vdum
  real(lk),dimension(0:mpsi) :: pmark3d, pmarktmp, marknum0,deltamp
  real(lk) zetadum, volume_3d, vvdum,vdum, whyx,whyz,ran
  integer idest,isource,isendtag,irecvtag,istatus(MPI_STATUS_SIZE)
  real(lk) sendl(mgrid),recvr(mgrid)
  integer addmp(0:mpsi,0:npartdom-1)
#endif

  logical subcycle,trapped,fkhybrid
  subcycle=.false.
  if(present(ncyclep))then
    if(ncyclep>0) subcycle=.true.
  endif

  trapped=.false.
  if(present(trap))then
    if(trap>0) trapped=.true.
  endif

  fkhybrid=.false.
  if(present(hybrid))then
    if(hybrid>0) fkhybrid=.true.
  endif

  if(trapped)then
    allocate(trapfrac(0:mpsi),STAT=mstat)
    call check_allocation(mstat, "trapfrac in load.F90")
    dden=0.0_lk
    ddum=0.0_lk
    icount=mpsi+1
    trapfracn=1.0_lk
    trapfrac=1.0_lk
  endif

  pflux=0.0
  rdtem=0.0
  rdden=0.0_lk
  zonal=0.0
  zonalc=0.0
  pmark=0.0
  marker=0.0
  markert=0.0
  density=0.0
  flow=0.0
  data1d=0.0

  zpart=0.0
  zpart0=0.0
  if(subcycle) zpart1=0.0
  wzpart=0.0
  wppart=0.0
  jtpart0=0
  jtpart1=0
  wtpart0=0.0
  wtpart1=0.0

  if(fkhybrid)then
    pressurepara=0.0
    pressureperp=0.0
    phisave=0.0
    dnsave=0.0
  endif

  if(irun>0 .or. pload==0)return


  delp=1.0/deltap
#ifdef _TOROIDAL3D
  interp=10
  mm=int(interp/2+1)
  allocate(jacobian(mgrid,interp+2),delvolume(mgrid,interp+2),delvolume_2(mgrid,mm),STAT=mstat)
  allocate(vol1(mgrid),vol2(mgrid),marker3d(mgrid),STAT=mstat)
  !allocate(parnum(mgrid,interp),STAT=mstat)
  !allocate(vdum(mtoroidal),STAT=mstat)
  
  jacobian=0.0_lk
  delvolume=0.0_lk
  delvolume_2=0.0_lk
  volume_3d=0.0_lk !! whole machine volume            
#endif

  !========================================================
  !> @note
  !> Calculation of the proper normalization of the marker weights
  !> In the end, marker(ij) contains the ratio between 
  !>      total particle number: Np and total marker number: Nm, 
  !>      divided by local particle number
  !> i.e. marker(ij) = Np/Nm / (n(i)*J(ij)*dpsi*dtheta) \n
  !> In charge, the gathered particle weight multiplied with marker(ij) and
  !> meshn(i) gives the correct density, since that is the physical particle number
  !> divided by the volume of the cell.

  ! First, marker(ij) contains the local physical particle number for each cell
  ! Jacobian=(gq+I)/B^2


#ifndef _TOROIDAL3D
!2D loading part

!$omp parallel do private(i,delpsi,pdum,isp,dpx,functiong,functioni,q,j,ij,b)
  do i=0,mpsi
    delpsi=0.5*(psimesh(min(mpsi,i+1))-psimesh(max(0,i-1)))
    pdum=psimesh(i)
    isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))
    dpx=pdum-spdpsi*real(isp-1)
    functiong=gpsi(1,isp)+dpx*gpsi(2,isp)+dpx*dpx*gpsi(3,isp)
    functioni=cpsi(1,isp)+dpx*cpsi(2,isp)+dpx*dpx*cpsi(3,isp)
    q=qmesh(i)
    do j=1,mtheta(i)
      ij=igrid(i)+j
      b=bmesh(ij)
      marker(ij)=meshn(i)*delpsi*deltat(i)*(functiong*q+functioni)/(b*b)
    enddo
  enddo

  ! # of physical particles of per annulus
  ! cannot use OpenMP in this loop
  pmark=0.0
  do i=0,mpsi
    do j=1,mtheta(i)
      ij=igrid(i)+j
      pmark(i)=pmark(i)+marker(ij)
    enddo
  enddo

  ! multiply by the ratio of total marker # and total physical particle
  ! number
  pdum=sum(pmark)
  tdum=real(mp)*real(npartdom)
  !$omp parallel do private(i,j,ij)
  do i=0,mpsi
    do j=1,mtheta(i)
      ij=igrid(i)+j
      marker(ij)=marker(ij)*tdum/pdum ! # of marker per cell
    enddo
    ! igrid(i) and igrid(i)+mtheta(i) are duplicated grid
    marker(igrid(i))=marker(igrid(i)+mtheta(i))
  enddo

#ifdef _FRC
  if(solvermethod==0)then
    do j=1,mtheta(1)
       ij=igrid(1)+j
       marker(ij)=marker(ij)+marker(igrid(0)+j)+marker(igrid(mpsi)+j)
    enddo
  endif
#endif

  ! Finally take the inverse to get the desired formula
  !$omp parallel do private(i,j,ij)
  do i=0,mpsi
    do j=1,mtheta(i)
      ij=igrid(i)+j
      marker(ij)=1.0/marker(ij) !to avoid divide operation
    enddo
    marker(igrid(i))=marker(igrid(i)+mtheta(i))
  enddo

  marknum=int(real(mp)*pmark/sum(pmark)) !# of marker per annulus per MPI
  pmark=real(marknum*numberpe)  !# of marker per annulus

  if(pload==2)then !for loading uniform marker in space 
    ! Nm/Np * integral dtheta (J dpsi)
    markload=pmark/meshn
    !# of marker per annulus per MPI, proportional to integral dtheta (J dpsi)
    marknum=int(real(mp)*markload/sum(markload))
    ! pdum is the volume averaged particle density
    pdum=sum(marknum*meshn)/sum(marknum)
    ! pmark is now n(i)* number of marker per annulus 
    pmark=pmark*pdum
    ! In this case, marker = Volume/Nm / (n(i)*J*dpsi*dtheta) 
    ! particle weight will be normalized accordingly later in zpart7
    ! initialization
    marker=marker/pdum
  endif
#else
!!!!!!!! 3D LOADING Part 1: 
!!!!!!!! added by HONGYU WANG, MAY 14TH, 2018 
!!!!!!!! 1. calculate jacobian volume delvolume (for calculating paticle number) and delvolume_2 (for getting pmark to do normalization)
!!!!!!!! 2. calculate marker and pmark  
!!!!!!!! 3. calculate total volume size of stellarator
  
!$omp parallel do private(i,delpsi,pdum,functiong,functioni,q,j,ij,b,k_sub,tdum,zetadum)
   do k_sub=1, interp+2
      zetadum=zeta1-(zeta1-zeta0)/(interp+1)*(k_sub-1)
        do i=0,mpsi
          delpsi=0.5*(psimesh(min(mpsi,i+1))-psimesh(max(0,i-1)))
          pdum=psimesh(i)
          functiong = spgpsi(pdum)
          functioni = spcpsi(pdum)
          q=spq(pdum)
            do j=1,mtheta(i)
              ! zeta1 is used here, instead of zetadum or zeta0
              ! the 1:interp/2+1 part are used for normalization later
              if(fielddir==1 .or. fielddir==3)then
                tdum=modulo(deltat(i)*real(j)+(zeta1-pi2)*qtinv(i),pi2)
              else
                tdum=modulo(deltat(i)*real(j)+zeta1*qtinv(i),pi2)
              endif
              b=spb(pdum,tdum,zetadum)
              ij=igrid(i)+j
              jacobian(ij,k_sub)=(functiong*q+functioni)/(b*b)
              delvolume(ij,k_sub)=meshn(i)*jacobian(ij,k_sub)*delpsi*deltat(i)
           enddo
         enddo
   enddo


!$omp parallel do private(i,delpsi,pdum,functiong,functioni,q,j,ij,b,k_sub,tdum,zetadum)
  do k_sub=1,interp/2+1
    zetadum=zeta1-(zeta1-zeta0)/(interp+1)*(k_sub+interp/2+1)
      do i=0,mpsi
        delpsi=0.5*(psimesh(min(mpsi,i+1))-psimesh(max(0,i-1)))                                            
        pdum=psimesh(i)
        functiong = spgpsi(pdum)
        functioni = spcpsi(pdum)
        q=spq(pdum)
          do j=1,mtheta(i)
            ! zeta0 ia used here, and this part will be given to left mpi 
            if(fielddir==1 .or. fielddir==3)then
              tdum=modulo(deltat(i)*real(j)+(zeta0-pi2)*qtinv(i),pi2)
            else
              tdum=modulo(deltat(i)*real(j)+zeta0*qtinv(i),pi2)
            endif
            b=spb(pdum,tdum,zetadum)
            ij=igrid(i)+j
            jacobian(ij,k_sub)=(functiong*q+functioni)/(b*b)
            delvolume_2(ij,k_sub)=meshn(i)*jacobian(ij,k_sub)*delpsi*deltat(i)
          enddo
       enddo
    enddo

!!!! vol1 and vol2 are used for different mpi for normalization.
!!!! vol2 will be given to left mpi, and current mpi will receive vol2 from right mpi.
  vol1=0.0_lk
  vol2=0.0_lk
  !$omp parallel do private(i,j,ij,k_sub)
  do i=0,mpsi
     do j=1,mtheta(i)
        ij=igrid(i)+j
        do k_sub=1,interp/2+1
          if(k_sub==1)then
          vol1(ij)=vol1(ij)+delvolume(ij,k_sub)/2
          else
           vol1(ij)=vol1(ij)+delvolume(ij,k_sub)
          endif
        enddo
        do k_sub=interp/2+2,interp+2
          if(k_sub==interp+2)then
            vol2(ij)=vol2(ij)+delvolume_2(ij,k_sub-interp/2-1)/2
          else
            vol2(ij)=vol2(ij)+delvolume_2(ij,k_sub-interp/2-1)
          endif
        enddo
     enddo
  enddo

  sendl=vol2
  recvr=0.0_lk
  icount=mgrid
  idest=left_pe
  isource=right_pe
  isendtag=myrank_toroidal
  irecvtag=isource
  call MPI_SENDRECV(sendl,icount,mpi_Rsize,idest,isendtag,&
         recvr,icount,mpi_Rsize,isource,irecvtag,toroidal_comm,istatus,ierror)
  vol2=recvr ! right of zeta1

 !$omp parallel do private(i,j,ij)
  do i=0,mpsi
     do j=1,mtheta(i)
        ij=igrid(i)+j
        marker(ij)=vol1(ij)+vol2(ij) ! an update of 2d situation, considering 3d rotation, using grid-center 
     enddo
    marker(igrid(i))=marker(igrid(i)+mtheta(i))
  enddo
                                                      
  !$omp parallel do private(i,j,ij)
  do i=0,mpsi
    do j=1,mtheta(i)
      ij=igrid(i)+j
      pmark(i)=pmark(i)+marker(ij)
    enddo
  enddo

!!!! marker3d is used to calculate pmark3d, and then calculate particle number
!from zeta0 to zeta1.
  marker3d=0.0_lk
  !$omp parallel do private(i,j,ij,k_sub)
  do i=0,mpsi
     do j=1,mtheta(i)
        ij=igrid(i)+j
        do k_sub=1,interp+2
          if(k_sub==1 .or. k_sub==interp+2)then
          marker3d(ij)=marker3d(ij)+delvolume(ij,k_sub)/2
          else
          marker3d(ij)=marker3d(ij)+delvolume(ij,k_sub)
          endif
        enddo
     enddo
    marker3d(igrid(i))=marker3d(igrid(i)+mtheta(i))
  enddo
                                                      
  pmark3d=0.0_lk
!$omp parallel do private(i,j,ij)
  do i=0,mpsi
    do j=1,mtheta(i)
      ij=igrid(i)+j
      pmark3d(i)=pmark3d(i)+marker3d(ij)
    enddo
  enddo

!!!!!!!! 3D LOADING Part 2: 
!!!!!!!! added by HONGYU WANG, MAY 14TH, 2018 
!!!!!!!! 1. calculate pmark for normalization
!!!!!!!! 2. calculate particle every partial_domaian: marknum
!!!!!!!! 3. calculate particle  correction every partial_domaian: addmp
   
  
  icount=1
  vdum=sum(pmark) !!!vdum is totalvolume of one toroidal plane in 3D.
  vvdum=0.0_lk
  call MPI_ALLREDUCE(vdum,vvdum,icount,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
  volume_3d=vvdum !!! volume_3d is the whole size of the stellarator
  
  
  pdum=volume_3d
  tdum=real(mp)*real(mtoroidal)*real(npartdom)
  !$omp parallel do private(i,j,ij)
  do i=0,mpsi
     do j=1,mtheta(i)
        ij=igrid(i)+j
        marker(ij)=marker(ij)*tdum/pdum
     enddo
    marker(igrid(i))=marker(igrid(i)+mtheta(i))
  enddo                                               
  
!$omp parallel do private(i,j,ij)
  do i=0,mpsi
     do j=1,mtheta(i)
        ij=igrid(i)+j
        marker(ij)=1.0/marker(ij)
     enddo
    marker(igrid(i))=marker(igrid(i)+mtheta(i))
  enddo    

  icount=mpsi+1
  pmark=pmark*real(mp)/volume_3d*mtoroidal*npartdom
  pmarktmp=0.0_lk
  call MPI_ALLREDUCE(pmark,pmarktmp,icount,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
  pmark=pmarktmp !! used for normalizetaion, from half to half.
  

  !!!! mp is total particle # of one partial of all mpi, from zeta0 to zeta1
  ! rough particle number of ONE MPI for a partdomain
  !# of marker per annulus per MPI. Using real number in mp to avoid integer overflow
  marknum=int(mtoroidal*real(mp)*pmark3d/volume_3d+0.5)
#endif    

  !initial density perturbation for fullk simulation
  init=0
  if(pload==9 .and. init==1)then
    do i=0,mpsi
      xdum=real(i)/real(mpsi)*PI
      xdum=1.0+1e-2*sin(xdum)
      perturbpmark(i)=pmark(i)*xdum
    enddo
    perturbmarknum=int(real(mp)*perturbpmark/sum(pmark))
    perturbpmark=real(perturbmarknum*numberpe)
    mp=sum(perturbmarknum)
  else
    perturbmarknum=marknum
  endif

  if(ipoincare==1)then
    call poincareLoad
    return
  endif

  mp=sum(perturbmarknum)
  delr=1.0/deltar
  spdt_inv=1.0/spdtheta

  ! random # uniformly distributed between 0 and 1
  call random_number(zpart(1,1:mp))  !psi
  call random_number(zpart(2,1:mp))  !theta
  call random_number(zpart(3,1:mp))  !zeta
  call random_number(zpart(4,1:mp))  !rho_para
  call random_number(zpart(5,1:mp))  !marker weight
  call random_number(zpart(6,1:mp))  !sqrt mu

  mmin=0
  mmax=0
  do ir=0,mpsi
    ! range of psi
    pmin=psimesh(max(0,ir-1))
    pmax=psimesh(min(mpsi,ir+1))

    ! range of particle #
    mmin=mmax+1
    mmax=mmin-1+perturbmarknum(ir)

    ! load particles in radial direction: linear in psi within a grid cell
    !$omp parallel do private(m,pdum)
    do m=mmin,mmax
      pdum=zpart(1,m)
      zpart(1,m)=(1.0-pdum)*pmin+pdum*pmax
    enddo

    ! Jacobian as function of theta on each flux surface & its inversion
    volt=0.0
    tvol=0.0
    ! cannot use OpenMP in this loop because volt(j+1) must be calculated
    ! after volt(j)
    do j=2,lst
      pdum=psimesh(ir)
      tdum=spdtheta*(real(j-1)-0.5)
#ifdef _TOROIDAL3D
      b = spb(pdum, tdum)
#else
      isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))
      dpx=pdum-spdpsi*real(isp-1)
      if(isp==1)dpx=sqrt(dpx)
      dp2=dpx*dpx

      jst=max(1,min(lst-1,ceiling(tdum*spdtheta_inv)))
      dtx=tdum-spdtheta*real(jst-1)
      dt2=dtx*dtx
      b=    bsp(1,isp,jst)+bsp(2,isp,jst)*dpx+bsp(3,isp,jst)*dp2+ &
           (bsp(4,isp,jst)+bsp(5,isp,jst)*dpx+bsp(6,isp,jst)*dp2)*dtx+ &
           (bsp(7,isp,jst)+bsp(8,isp,jst)*dpx+bsp(9,isp,jst)*dp2)*dt2
#endif
      ! volt is the CDF of a PDF proportional to the Jacobian
      volt(1,j)=volt(1,j-1)+1.0/(b*b)
    enddo
    ! Normalization of volt. Note the additional 2pi for inverse-method to
    ! generate a random variable in [0,2pi)
    volt=2.0*pi*volt/volt(1,lst)
    call construct_spline1d(lst,spdtheta,volt,0)
    ! The inverse-CDF method to generate a random variable according to the
    ! given PDF proportional to the jacobian
    call invert_spline0(0,lst,spdtheta,spdtheta,volt,tvol)

    ! load particles in theta: uniform density in theta
    ! evaluate tvol on a uniform random variable
    !$omp parallel do private(m,tdum,jt,dtx)
    do m=mmin,mmax
      tdum=2.0*pi*zpart(2,m)
      jt=max(1,min(lst-1,ceiling(tdum*spdt_inv)))
      dtx=tdum-spdtheta*real(jt-1)
      zpart(2,m)=tvol(1,jt)+dtx*tvol(2,jt)+dtx*dtx*tvol(3,jt)
    enddo
  enddo

  !> @todo
  !> Potential bug: In 3D magnetic equilibrium, Boozer coordinates Jacobian is
  !> in not uniform along zeta. Strictly speaking, it is not uniform in psi
  !> either. A more general way of generating the initial distribution is
  !> required. Maybe Von-Neumann Rejection method?

  ! load particles uniformly in zeta
  !$omp parallel do private(m,tdum)
  do m=1,mp
    tdum=zpart(3,m)*deltaz
    zpart(3,m)=zeta0+tdum
  ! initialize random weight
    zpart(5,m)=w_initial*2.0*(zpart(5,m)-0.5)
#ifdef _FRC
  ! uncomment to set particular n mode
    zpart(5,m)=zpart(5,m)*cos(zpart(3,m)*real(nmode))
#else
  ! set weight more heavily towards outer midplane
    zpart(5,m)=zpart(5,m)*(1.0+cos(zpart(2,m)))
#endif
  enddo

  if(pload==11)then ! Slowing down distribution
    fmax = sd_v0**2/(sd_v0**2+sd_vc**2)
  !$omp parallel do private(m,zdum,pdum,tdum,isp,dpx,dp2,jst,dtx,dt2,b)
    do m=1,mp
      randvarf = 1000.0_lk
      sd_ftmp  = 0.0001_lk
      do while(randvarf>sd_ftmp)
        call random_number(randvarv)
        call random_number(randvarf)
        call random_number(randvart)
        randvarf=fmax*randvarf
        randvarv=sd_v0*randvarv
        randvart=pi*randvart

        pdum=zpart(1,m)
        tdum=zpart(2,m)

        ! use spb function to automatically deal with both 2D and 3D case.
        ! in 3D, default zeta = zeta1 will be used.
        b = spb(pdum, tdum)

        sd_ftmp = exp(-((sin(randvart)**2/b-sd_l0)*sd_widthInv)**2)*randvarv**2/(sd_vc**3 +randvarv**3)*sin(randvart)
      enddo
        zpart(4,m)= randvarv*cos(randvart)
        zpart(6,m)= (randvarv*sin(randvart))**2
        zpart(4,m)=zpart(4,m)*apart/(qpart*b)
        zpart(6,m)=sqrt(0.5_lk*apart*zpart(6,m)/b)
        zpart0(4,m)=1.0_lk
        zpart0(6,m)=1.0_lk
    enddo
  else
  ! Maxwellian distribution in v_para, <v_para^2>=1.0, use zpartsign as temporary storage
  ! The generation of standard normal distributed random variable from uniform
  ! random variable in (0,1] is done by using a rational formula approximation
  ! of the inverse CDF of the normal distribution.

  ! Reference temperature are chosen to be at either axis or eq_flux
  ! consistent with the normalization scheme controlled by inorm
  if(inorm==0) then
    tpref = tppp(1,1)
  else
    tpref = mesht(eq_flux)
  endif
  tp_inv = 1.0_lk/tpref

#ifdef _TOROIDAL3D
  !$omp parallel do private(m,zdum,pdum,tdum,zetadum,isp,dpx,dp2,jst,dtx,dt2,zpartsign,b,majorr,upara0,rg,ii,dp1,vth)
#else
!$omp parallel do private(m,zdum,pdum,tdum,isp,dpx,dp2,jst,dtx,dt2,zpartsign,b,majorr,upara0,rg,ii,dp1,vth)
#endif
   do m=1,mp
     ! initialize vth as thermal velocity at reference point
     ! vth may be modified to local thermal velocity in pload=2,3 cases
     vth=sqrt(tpref/apart)
     zdum=zpart(4,m)
     zpart(4,m)=zpart(4,m)-0.5
     zpartsign=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)
     ! The velocities outside the chosen range are reset to near the center.
     !> @todo
     !> Double check if this reset scheme is desired. All the out of range
     !> particles will stack near vpara=0.5
     if(zpart(4,m)>umax)zpart(4,m)=zdum
     zpart(4,m)=zpartsign*zpart(4,m)

    ! Maxwellian distribution in v_perp, <v_perp^2>=1.0
    ! v_perp^2 obey exponential distribution. The inverse-CDF is very easily
    ! obtained as -log(u)
     zpart(6,m)=max(1.0e-20_lk,min(umax*umax,-log(max(1.0e-20_lk,zpart(6,m)))))

    ! zpart(4,:) to rho_para, zpart(6,:) to sqrt(mu)
     pdum=zpart(1,m)
     tdum=zpart(2,m)
#ifdef _TOROIDAL3D
     zetadum=zpart(3,m)

     ! toroidal rotation
     upara0 = sptrot(pdum)
     ! magnetic field strength
     b = spb(pdum, tdum, zetadum)
     ! major radius
     majorr = spx(pdum, tdum, zetadum)
#else
     isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))
     dpx=pdum-spdpsi*real(isp-1)
     if(isp==1)dpx=sqrt(dpx)
     dp2=dpx*dpx

     jst=max(1,min(lst-1,ceiling(tdum*spdtheta_inv)))
     dtx=tdum-spdtheta*real(jst-1)
     dt2=dtx*dtx
     b=    bsp(1,isp,jst)+bsp(2,isp,jst)*dpx+bsp(3,isp,jst)*dp2+ &
          (bsp(4,isp,jst)+bsp(5,isp,jst)*dpx+bsp(6,isp,jst)*dp2)*dtx+ &
          (bsp(7,isp,jst)+bsp(8,isp,jst)*dpx+bsp(9,isp,jst)*dp2)*dt2
     majorr= xsp(1,isp,jst)+xsp(2,isp,jst)*dpx+xsp(3,isp,jst)*dp2+ &
          (xsp(4,isp,jst)+xsp(5,isp,jst)*dpx+xsp(6,isp,jst)*dp2)*dtx+ &
          (xsp(7,isp,jst)+xsp(8,isp,jst)*dpx+xsp(9,isp,jst)*dp2)*dt2                                      
     upara0=ropp(1,isp)  +ropp(2,isp)*dpx   +ropp(3,isp)*dp2
#endif

     upara0=majorr*upara0 !convert rotation frequency to rotation speed

     zpart(4,m)=apart*(zpart(4,m)*vth+upara0)/(qpart*b)   !zpart4 is rho_para
     zpart(6,m)=vth*sqrt(apart*zpart(6,m)/b)          !zpart6 is sqrt(mu)
     ! zpart7 is the f/g term, the ratio between total physical particle
     ! number and numerical marker number is included in marker(ij), which also has
     ! the volume size of each cell, so the physical "density" can be easily
     ! calculated
     ! For iload=1,3, markers are distributed proportionally to local physical
     ! particle number, so the ratio is constant, and zpart7=1
     zpart(7,m)=1.0_lk

     ! pload=2: marker distribution non-uniform in v, but uniform in space;
     !          marker weight zpart(7,:)=n_physical/n_marker
     ! pload=3: marker distribution non-uniform in both v & space;
     !          marker weight zpart(7,:)=1.0
     if(pload==2 .or. pload==3)then
       rg=  rgpsi(1,isp)  +rgpsi(2,isp)*dpx  +rgpsi(3,isp)*dp2
       !radial grid on inner flux surface
       ii=max(1,min(mpsi,ceiling((rg-rg0)*delr)))
       dp1=(psimesh(ii)-pdum)*delp(ii) !weight for inner flux surface

       vth=sqrt((dp1*mesht(ii-1)+(1.0-dp1)*mesht(ii))*tp_inv)
       zpart(4,m)=zpart(4,m)-apart*upara0/(qpart*b)
       zpart(4,m)=vth*zpart(4,m)+apart*upara0/(qpart*b)  !use global temperature
       zpart(6,m)=vth*zpart(6,m)                         !use global temperature

       ! The f/g term needs to be rescaled to be proportional to the local
       ! density, note that the volume/Nm factor is in marker(ij) expression
       if(pload==2)zpart(7,m)=(dp1*meshn(ii-1)+(1.0-dp1)*meshn(ii))
     endif
   enddo
  endif


! trap=1: load trapped particles only
! pick trapped particles using trapping condition
!
  if(trapped .and. trap==1)then
     cmratio=qpart/apart
     mptmp=0
     trapfracn=real(mp,lk)
     trapfrac=0.0

! do not use omp here
     do m=1,mp
        pdum=zpart(1,m)
        tdum=zpart(2,m)

#ifdef _TOROIDAL3D
        b = spb(pdum, tdum)
#else
        isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))
        dpx=pdum-spdpsi*real(isp-1)
        if(isp==1)dpx=sqrt(dpx)
        dp2=dpx*dpx

        jst=max(1,min(lst-1,ceiling(tdum*spdtheta_inv)))
        dtx=tdum-spdtheta*real(jst-1)
        dt2=dtx*dtx
        b=    bsp(1,isp,jst)+bsp(2,isp,jst)*dpx+bsp(3,isp,jst)*dp2+ &
             (bsp(4,isp,jst)+bsp(5,isp,jst)*dpx+bsp(6,isp,jst)*dp2)*dtx+ &
             (bsp(7,isp,jst)+bsp(8,isp,jst)*dpx+bsp(9,isp,jst)*dp2)*dt2
#endif

        upara=zpart(4,m)*b*cmratio
        energy=0.5*apart*upara*upara+zpart(6,m)*zpart(6,m)*b

#ifdef _TOROIDAL3D
        rg=  sprgpsi(pdum)
#else
        rg=  rgpsi(1,isp)  +rgpsi(2,isp)*dpx  +rgpsi(3,isp)*dp2
#endif
        ii=max(1,min(mpsi,int((rg-rg0)*delr+0.5)))
        bm=maxval(bmesh(igrid(ii):igrid(ii)+mtheta(ii))) !maximal B-field on flux-surface
        eperp=zpart(6,m)*zpart(6,m)*bm

! radial profile of particle and energy flux
        fullf=zpart(7,m)
        dp1=(psimesh(ii)-pdum)*delp(ii) !weight for inner flux surfac
        dden(ii-1)=dden(ii-1)+fullf*dp1
        dden(ii)  =dden(ii)+fullf*(1.0-dp1)

        if(eperp>energy)then
           mptmp=mptmp+1
           zpart(1,mptmp)=zpart(1,m)
           zpart(2,mptmp)=zpart(2,m)
           zpart(3,mptmp)=zpart(3,m)
           zpart(4,mptmp)=zpart(4,m)
           zpart(5,mptmp)=zpart(5,m)
           zpart(6,mptmp)=zpart(6,m)
           zpart(7,mptmp)=fullf !!zpart(7,m)

           trapfrac(ii-1)=trapfrac(ii-1)+fullf*dp1
           trapfrac(ii)  =trapfrac(ii)+fullf*(1.0-dp1)

        endif
     enddo
     mp=mptmp

     call MPI_ALLREDUCE(dden,ddum,icount,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)
     dden=ddum
     call MPI_ALLREDUCE(trapfrac,ddum,icount,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)
     trapfrac=ddum/dden
     call MPI_ALLREDUCE(trapfracn,tdum,1,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)
     zdum=real(mp,lk)
     call MPI_ALLREDUCE(zdum,trapfracn,1,mpi_Rsize,MPI_SUM,MPI_COMM_WORLD,ierror)
     trapfracn=trapfracn/tdum
   endif
   !$acc update device(pmark)
end subroutine gkLoadParticle

subroutine gk2fk
  use utility, only: check_allocation, mstat
  use global_parameters
  use particle_array
  use field_array,only: nmodes,nmode
  use equilibrium
  use spline_function, only: dxdp, dxdt, dzdp, dzdt, spx, spb, spgpsi,spcpsi
  implicit none


  integer m,isp,jst
  real(lk) pdum,tdum,dpx,dp2,dtx,dt2,gradt,gradz,gradp,b,ri,g,dxdptmp,dxdttmp,dzdptmp,dzdttmp,&
       pdot,tdot
  real(lk),dimension(:),allocatable:: alpha

  interface
    subroutine gkLoadParticle(meshn,tppp,nppp,zpart,zpart0,wppart,&
        wtpart0,wtpart1,density,flow,data1d,jtpart0,jtpart1,wzpart,&
        marker,markert,pmark,pflux,rdtem,rdden,zonal,zonalc,mesht,qpart,apart,&
        pload,nparam,ngyro,mp,mpmax,w_initial,ncyclep,zpart1,&
        trap,trapfracn,trapfrac,hybrid,pressurepara,pressureperp,phit,dnt,&
        phisave,dnsave)

      use precision, only: lk
      use spline_function, only: construct_spline1d, invert_spline0
      implicit none
      integer pload,nparam,ngyro,mp,mpmax
      integer,optional :: ncyclep,trap,hybrid
      integer,dimension(:,:) :: jtpart0,jtpart1
      real(lk) qpart,apart,w_initial
      real(lk),optional :: trapfracn
      real(lk),dimension(0:) :: meshn,mesht,pmark,pflux,rdtem,rdden,zonal,zonalc
      real(lk),dimension(:) :: wzpart,marker,markert
      real(lk),dimension(:),allocatable,optional :: trapfrac
      real(lk),dimension(:,:) :: tppp,nppp
      real(lk),dimension(:,:) :: zpart,zpart0,wppart,wtpart0,wtpart1
      real(lk),dimension(0:,:) :: density,flow,data1d
      real(lk),dimension(:,:),optional :: zpart1
      real(lk),dimension(0:,:),optional :: pressurepara,&
        pressureperp,phit,dnt
      real(lk),dimension(0:,:,:),optional :: phisave,dnsave
    end subroutine gkLoadParticle
  end interface

  allocate(alpha(mimax), STAT=mstat)
  call check_allocation(mstat, "alpha in load.F90")
  call random_number(alpha(mimax))
  alpha(:)=alpha(:)*2*pi

  !zion(1:3,m)->(psi,theta,zeta)
  !zion(4,m)->pdot
  !zion(5,m)->particle weight
  !zion(6,m)->(tdot)
  !zion(7,m)->
  !zion(8,m)->(zdot)
  !$omp parallel do private(m,pdum,isp,dpx,dp2,tdum,jst,dtx,dt2,&
  !$omp& gradt,gradz,gradp,b,ri,g,dxdptmp,dxdttmp,dzdptmp,dzdttmp)
  do m=1,mi
      pdum=zion(1,m)
      tdum=zion(2,m)
      dxdptmp=dxdp(pdum,tdum)
      dxdttmp=dxdt(pdum,tdum)
      dzdptmp=dzdp(pdum,tdum)
      dzdttmp=dzdt(pdum,tdum)
#ifdef _TOROIDAL3D
      b = spb(pdum, tdum)
      g = spgpsi(pdum)
      ri = spcpsi(pdum)
#else
      isp=max(1,min(lsp-1,ceiling(pdum*spdpsi_inv)))
      dpx=pdum-spdpsi*real(isp-1)
      if(isp==1)dpx=sqrt(dpx)
      dp2=dpx*dpx
      jst=max(1,min(lst-1,ceiling(tdum*spdtheta_inv)))
      dtx=tdum-spdtheta*real(jst-1)
      dt2=dtx*dtx
      b=bsp(1,isp,jst)+bsp(2,isp,jst)*dpx+bsp(3,isp,jst)*dp2+ &
        (bsp(4,isp,jst)+bsp(5,isp,jst)*dpx+bsp(6,isp,jst)*dp2)*dtx+ &
        (bsp(7,isp,jst)+bsp(8,isp,jst)*dpx+bsp(9,isp,jst)*dp2)*dt2
      g=gpsi(1,isp)   +gpsi(2,isp)*dpx   +gpsi(3,isp)*dp2
      ri=cpsi(1,isp)   +cpsi(2,isp)*dpx   +cpsi(3,isp)*dp2   
#endif
      gradp=sqrt(dxdttmp**2+dzdttmp**2)/abs(dxdptmp*dzdttmp-dxdttmp*dzdptmp)
      gradt=sqrt(dxdptmp**2+dzdptmp**2)/abs(dxdptmp*dzdttmp-dxdttmp*dzdptmp)
      gradz=1.0_lk/spx(pdum,tdum)
      pdot=sqrt(zion(6,m)**2*b*2.0_lk/aion)*sin(alpha(m))*gradp
      tdot=(zion(4,m)*qion/aion*ri*gradt*gradz+&
                sqrt(zion(6,m)**2*b*2.0_lk/aion)*g*gradz/b*cos(alpha(m)))*gradt
      zion(8,m)=(zion(4,m)*qion/aion*g*gradz-&
                sqrt(zion(6,m)**2*b*2/aion)*ri*gradt*gradz/b*cos(alpha(m)))*gradz
      zion(4,m)=pdot
      zion(6,m)=tdot
  enddo

end subroutine gk2fk


! subroutine poincareLoad
!> \brief load particles in the way of poincare mode
!> \details load particle anywhere you set in the simulation area

subroutine poincareLoad
  use precision
  use utility, only: check_allocation, mstat
  use global_parameters,only : ndiag,nhybrid,irestore,mype,mpsi,eq_flux
  use particle_array,only : mi,zion,meshti,aion,fload,feload
  use field_array,only : zeta0,psimesh
  implicit none

  integer m
  real(lk) vth,pmin,pmax,pdum
! reset default value when doing poincare
! ndiag=1,so that we wouldn't miss the particle when we aren't diagnosing them
   ndiag=1
! load ion only when doing poincare
   nhybrid=0
   fload=0
   feload=0
! the parameter could affect the weight until it overflow
   irestore=0
! load particle(ion) in first MPI(mype==0) only
  if(mype==0)then
! num of particle(ion)
! mi can be modified based on needs
     mi=100
  else
     mi=0
  end if
! load particle in radial direction: linear between 0 and 1
  if (mype==0)then
    do m=1,mi
      zion(1,m)=real(m,lk)/(1.0_lk+real(mi,lk))
    enddo
  end if
! range of psi
  pmin=psimesh(0)
  pmax=psimesh(mpsi)
! load particles in radial direction: linear in psi
  !$omp parallel do private(m,pdum)
  do m=1,mi
    pdum=zion(1,m)
    zion(1,m)=(1.0-pdum)*pmin+pdum*pmax
  enddo
! load particles in theta: theta=0
! load particles in zeta: zeta=zeta0
! load weight: weight=0.0 -- the value is of no use when doing poincare if we don't solve poisson eq.
  !$omp parallel do private(m)
  do m=1,mi
    zion(2,m)=0.0_lk
    zion(3,m)=zeta0
    zion(5,m)=0.0_lk
  enddo
! load rho_para: vth*aion,the sign is related to direction
! load sqrt(mu)=0
! load PDF=1

  !$omp parallel do private(m,vth)
   do m=1,mi
     vth=sqrt(meshti(eq_flux)/aion)
     zion(4,m)=vth*aion
     zion(6,m)=0.0_lk
     zion(7,m)=1.0_lk
   enddo

end subroutine poincareLoad
