subroutine gkload(mp, mpmax, zpart, zpart0, wzpart, wppart, jtpart0, jtpart1,&
    wtpart0, wtpart1, jppart,rhogyropart, t_Eq, apart, qpart, marker, &
    pmark, meshn, mesht1d, nonuni, marknum, mpcore, mpsol, ngyro)
    use precision
    use grids_setup
    use grids_fields, only: phi
    use read_in_parameters, only: irun, frcfullf, uni_profile
    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 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
    real(lk), dimension(:) :: marker, meshn, mesht1d, nonuni, pmark
    integer, dimension(:) :: marknum
    integer mp, mpmax, mpcore, mpsol, ngyro
    real(lk) apart, qpart
! local parameters
    integer :: i,j,ij,ij0,ij1,ii,m,mtest,ierror,icount,mmin,mmax,nrsp,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,tp,&
             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


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
    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
!$acc update device(zpart)

    if(irun>0)return

    if(sfcore==1)then
! number of marker per grid, Jacobian=(gq+I)/B^2
        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
! number of marker per grid, Jacobian=(gq+I)/B^2
        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(mype==0)write(gtcout,*)"pdum",pdum
    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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! uni_profile==1: n_i = const., T_i = const., fullf = nonuni = 1.0.
    if(uni_profile==1)then
        pdum=sum(volume)
        if(sfcore==1)then
            do i=1,mpsi_core
                do j=1,mpara_core(i) - 1
                    ij=igrid_core(i)+j
                    marknum(ij)=max(1,int(volume(ij)*real(mp)/pdum)) ! load marker by volume
                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
                    marknum(ij)=max(1,int(volume(ij)*real(mp)/pdum)) ! load marker by volume
                enddo
            enddo
        endif
    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   ! Ghost grid of core region should be excluded
                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  ! Ghost grids in core region are also calculated
                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) ! # of marker per cell
            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(uni_profile==1)then
        nonuni=1.0_lk
    endif

    if(couplefrc==1)then
! Shared grids for SOL and core regions coupling
        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)  ! Periodicity boundary condition for ghost grid in core region
        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,1:mp))
    call random_number(zpart(2,1:mp))
    call random_number(zpart(3,1:mp))
    call random_number(zpart(4,1:mp))
    call random_number(zpart(5,1:mp))
    call random_number(zpart(6,1:mp))
    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*2.0_lk*(zpart(5,m)-0.5_lk)
                    zpart0(6,m)=nonuni(ij)
                    tdum=geotheta(zpart(1,m),zpart(3,m))
                    zpart(7,m)=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,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)))
!! 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,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*2.0_lk*(zpart(5,m)-0.5_lk)
                    zpart0(6,m)=nonuni(ij)
                    zpart(7,m)=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
                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
    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*mesht1d(ii+1)+(1-dp1)*mesht1d(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)
    
      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*mesht1d(ii+1+mpsi_core)+(1-dp1)*mesht1d(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)
      enddo
      if(mype==0)write(gtcout,*)"vth",sqrt(mesht1d(1+mpsi_core))/sqrt(apart)
    endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


    ! call locate_kernel(zpart,wzpart,wppart,wtpart0,wtpart1,jtpart0,&
    ! jtpart1,ngyro,mp,jppart,qpart,apart,rhogyropart, .false.)
!$acc update device(zpart)
    call locate_kernel(zpart,wzpart,wppart,wtpart0,wtpart1,jtpart0,&
    jtpart1,1,mp,jppart,qpart,apart,rhogyropart, .false.)
!$acc update self(wppart, wtpart0, wtpart1, jtpart0, jtpart1)

!$omp parallel do private(m,wp1,wp0,wt00,wt10,wt01,wt11,ij)
    do m=1,mp
        zpart0(6,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(6,m)=zpart0(6,m)+wt00*nonuni(ij)

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

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

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

    if(frcfullf==1)then
!$omp parallel do private(m)
        do m=1,mp
            zpart(5,m)=zpart0(6,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(zpart(1,m)) .or. isnan(zpart(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(zpart(1,m)) .or. isnan(zpart(3,m))))then
                    print *,"istep=",istep,"irk=",irk
                    print *,"mype",mype,"m=",m,"zpart(1,m)=",zpart(1,m),"zpart(3,m)=",zpart(3,m)
                    print *,"psidum=",psidum,"psidumtmp=",psidumtmp
                    print *,"tdum=",tdum,"rdum=",rdum,"zdum=",zdum
                    stop
                endif
            endif
        endif
    enddo

    if(mype==0)then
        write(gtcout,*)'mp=',mp
        write(gtcout,*)'mpcore=',mpcore
        write(gtcout,*)'mpsol=',mpsol
        write(gtcout,*)"mpmax=",mpmax
        call FLUSH(gtcout)
    endif

!$omp parallel do private(m)
    do m=1,mp
       zpart(5,m)=zpart(5,m)*zpart0(6,m)
       zpart(7,m)=spline2d(zpart(1,m),zpart(3,m)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    enddo
end subroutine gkload
