subroutine global_pssolver
#include "petscversion.h"
    use precision
    use global_parameters
    use grids_fields
    use particle_array
    use petsc_array, only: luserbfperp, luserxfperp
    use read_in_parameters
    use simulation_domain, only: sfcore, sfsol, couplefrc
    use grids_setup
    use grids_eqdata
    use cell_volume, only: volume
#if PETSC_VERSION_GT(3,8,0)
    use petsc_gtc, only: lapmat_pssolver_fperp, petsc_ps, petsc_ps0
#else
    use petsc_gtc, only: lapmat_pssolver_fperp, lapmat_pssolver_fperp0
#endif
    use laplacian, only: laplacian_perp
    use ieee_arithmetic
    use system_env
    implicit none
    
    integer :: i, j, ij, ij1, ijn, ierror, icount
    real(lk) :: dentmp(mgrid), dielectric_constant, coe_temp, ne_ni_ratio, &
                dentmp0(mgrid),adum(0:1,mpsi),total_v,data1dtmp(mpsi,2),total_v1d(mpsi)
    
    real(lk) :: dentmpr(mgrid, mtoroidal/2+1), dentmpi(mgrid, mtoroidal/2+1), &
                phitmpr(mgrid, mtoroidal/2+1), phitmpi(mgrid, mtoroidal/2+1),&
                dentmpr0(mgrid, mtoroidal/2+1), dentmpi0(mgrid, mtoroidal/2+1), &
                phitmpr0(mgrid, mtoroidal/2+1), phitmpi0(mgrid, mtoroidal/2+1)

    real(lk) :: ffilter(mtoroidal/2+1), &
                xz(mtoroidal),&
                phifluxlz(mgrid,mtoroidal),&
                phifluxlz0(mgrid,mtoroidal),&
                eachzphifluxlz_core(mygrid_core,mtoroidal),&
                eachzphifluxlz0_core(mygrid_core,mtoroidal),&
                eachzphifluxlz_sol(mygrid_sol,mtoroidal),&
                eachzphifluxlz0_sol(mygrid_sol,mtoroidal),&
                eachtphifluxlz_core(mgrid_core),&
                eachtphifluxlz0_core(mgrid_core),&
                eachtphifluxlz_sol(mgrid_sol),&
                eachtphifluxlz0_sol(mgrid_sol)
    real(lk) :: dentmpsend_core(mygrid_core*2), dentmprecv_core(mygrid_core*2,mtoroidal),dentmpsend_sol(mygrid_sol*2), dentmprecv_sol(mygrid_sol*2,mtoroidal)
    complex(lk) :: yz(mtoroidal/2+1), yz2(mtoroidal/2+1)

    if(ifluxtube==1)then
       if(mype==0)then
         write(gtcout,*) "use local solver for fluxtube simulation! set global=0 in gtc.in!"
       endif
       stop
    endif
    phitmpr = 0.0_lk
    phitmpi = 0.0_lk
    dentmpr = 0.0_lk
    dentmpi = 0.0_lk
    phitmpr0 = 0.0_lk
    phitmpi0 = 0.0_lk
    dentmpr0 = 0.0_lk
    dentmpi0 = 0.0_lk

    icount = mgrid

    if(sfcore == 1)then
        if(ifullk==0)then
        if( kinetice==0 )then
!$omp parallel do private(i, j, ij, ne_ni_ratio, coe_temp) 
            do i = 1, mpsi_core
                do j = 1, mpara_core(i)
                    ij = igrid_core(i) + j
                    ne_ni_ratio = meshne(ij)/meshni(ij)
                    coe_temp = 1.0_lk + meshti(ij)/meshte(ij)/qion/qion*ne_ni_ratio
                    dentmp(ij) = qion*meshni(ij)*densityi(1,ij)/coe_temp
                    dentmp0(ij) = qion*meshni(ij)*densityi0(1,ij)
                enddo
            enddo
        else
!$omp parallel do private(i, j, ij, ne_ni_ratio, coe_temp) 
            do i = 1, mpsi_core
                do j = 1, mpara_core(i)
                    ij = igrid_core(i) + j
                    ne_ni_ratio = meshne(ij)/meshni(ij)
                    dentmp(ij) = qion*meshni(ij)*densityi(1,ij) + qelectron*meshne(ij)*densitye(1,ij)
                    dentmp0(ij) = qion*meshni(ij)*densityi0(1,ij) + qelectron*meshne(ij)*densitye0(1,ij)
                enddo
            enddo
        endif
        else
!$omp parallel do private(i, j, ij, ne_ni_ratio, coe_temp) 
        do i = 1, mpsi_core
            do j = 1, mpara_core(i)
                ij = igrid_core(i) + j
                dentmp(ij) = qfki*meshni(ij)*densityfki(1,ij)
            enddo
        enddo
        endif
    endif

    if(sfsol == 1)then
        if(ifullk==0)then
        if(kinetice==0)then
!$omp parallel do private(i, j, ij, ne_ni_ratio, coe_temp) 
            do i = 1, mpsi_sol
                do j = 1, mpara_sol(i)
                    ij = mgrid_core + igrid_sol(i) + j
                    ne_ni_ratio = meshne(ij)/meshni(ij)
                    coe_temp = 1.0_lk + meshti(ij)/meshte(ij)/qion/qion*ne_ni_ratio
                    dentmp(ij) = qion*meshni(ij)*densityi(1,ij)/coe_temp
                    dentmp0(ij) = qion*meshni(ij)*densityi0(1,ij)
                enddo
            enddo
        else
!$omp parallel do private(i, j, ij, ne_ni_ratio, coe_temp) 
            do i = 1, mpsi_sol
                do j = 1, mpara_sol(i)
                    ij = mgrid_core + igrid_sol(i) + j
                    dentmp(ij) = qion*meshni(ij)*densityi(1,ij) + qelectron*meshne(ij)*densitye(1,ij)
                    dentmp0(ij) = qion*meshni(ij)*densityi0(1,ij) + qelectron*meshne(ij)*densitye0(1,ij)
                enddo
            enddo
        endif
        else
!$omp parallel do private(i, j, ij, ne_ni_ratio, coe_temp) 
        do i = 1, mpsi_sol
            do j = 1, mpara_sol(i)
                ij = mgrid_core + igrid_sol(i) + j
                dentmp(ij) = qfki*meshni(ij)*densityfki(1,ij)
            enddo
        enddo
        dentmp0=0.0_lk
        endif
    endif


if(sfcore==1)then
dentmpsend_core(1:mygrid_core)=dentmp(ij_min_core:ij_max_core)
dentmpsend_core(mygrid_core+1:mygrid_core*2)=dentmp0(ij_min_core:ij_max_core)
call MPI_ALLGATHER(dentmpsend_core,2*mygrid_core,mpi_Rsize,dentmprecv_core,2*mygrid_core,&
                       mpi_Rsize,toroidal_comm,ierror)
do ij=ij_min_core, ij_max_core
   phifluxlz(ij,:)=dentmprecv_core(ij-ij_min_core+1,:)
   phifluxlz0(ij,:)=dentmprecv_core(ij-ij_min_core+1+mygrid_core,:)
enddo
endif

if(sfsol==1)then
dentmpsend_sol(1:mygrid_sol)=dentmp(ij_min_sol:ij_max_sol)
dentmpsend_sol(mygrid_sol+1:mygrid_sol*2)=dentmp0(ij_min_sol:ij_max_sol)

call MPI_ALLGATHER(dentmpsend_sol,2*mygrid_sol,mpi_Rsize,dentmprecv_sol,2*mygrid_sol,&
                       mpi_Rsize,toroidal_comm,ierror)
do ij=ij_min_sol, ij_max_sol
   phifluxlz(ij,:)=dentmprecv_sol(ij-ij_min_sol+1,:)
   phifluxlz0(ij,:)=dentmprecv_sol(ij-ij_min_sol+1+mygrid_sol,:)
enddo
endif

!================== transform from real space to fourier space ==================

    if(sfcore == 1)then
        yz = 0.0_lk
        ffilter = 0.0_lk
        ffilter(n_modes+1) = 1.0_lk
!$omp parallel do private(ij, xz, yz)
        do ij = ij_min_tor_core, ij_max_tor_core
                xz = phifluxlz(ij,:)
                call fftr1d(1, mtoroidal, 1.0_lk, xz, yz, 2)
    ! linear run only keep a few modes
                yz = ffilter*yz
                dentmpr(ij, n_modes(1) + 1) = real(real(yz(n_modes(1) + 1)))
                dentmpi(ij, n_modes(1) + 1) = real(aimag(yz(n_modes(1) + 1)))
                xz = phifluxlz0(ij,:)
                call fftr1d(1, mtoroidal, 1.0_lk, xz, yz, 2)
    ! linear run only keep n=0 mode
                dentmpr0(ij, 1) = real(real(yz(1)))
                dentmpi0(ij, 1) = real(aimag(yz(1)))
        enddo
    endif
    
    if(sfsol == 1)then
        yz = 0.0_lk
        ffilter = 0.0_lk
        ffilter(n_modes + 1) = 1.0_lk
!$omp parallel do private(ij, xz, yz)
        do ij = ij_min_tor_sol, ij_max_tor_sol
                xz = phifluxlz(ij,:)
                call fftr1d(1, mtoroidal, 1.0_lk, xz, yz, 2)
    ! linear run only keep a few modes
                yz = ffilter*yz
                dentmpr(ij, n_modes(1) + 1) = real(real(yz(n_modes(1) + 1)))
                dentmpi(ij, n_modes(1) + 1) = real(aimag(yz(n_modes(1) + 1)))
                xz = phifluxlz0(ij,:)
                call fftr1d(1, mtoroidal, 1.0_lk, xz, yz, 2)
    ! linear run only keep a few modes
                dentmpr0(ij, 1) = real(real(yz(1)))
                dentmpi0(ij, 1) = real(aimag(yz(1)))
        enddo
    endif
!=============================================================================
!   Important: interplate the fields on irregular grids to regular grids
!              for SOL and core region coupled simulation
!   Because the laplacian operator is discretized based on the regular grids
!=============================================================================
    if(couplefrc==1)then 
        call convert_irreg2reg(dentmpr(:, n_modes(1) + 1))
        call convert_irreg2reg(dentmpi(:, n_modes(1) + 1))
        call convert_irreg2reg(dentmpr0(:, n_modes(1) + 1))
        call convert_irreg2reg(dentmpi0(:, n_modes(1) + 1))
    endif
    CALL MPI_ALLREDUCE(dentmpr(:,n_modes(1)+1),dentmp,mgrid,MPI_RSIZE, MPI_SUM, mpi_comm_world,ierror)
    dentmpr(:,n_modes(1)+1)=dentmp
    CALL MPI_ALLREDUCE(dentmpi(:,n_modes(1)+1),dentmp,mgrid,MPI_RSIZE, MPI_SUM, mpi_comm_world,ierror)
    dentmpi(:,n_modes(1)+1)=dentmp
    CALL MPI_ALLREDUCE(dentmpr0(:,n_modes(1)+1),dentmp0,mgrid,MPI_RSIZE, MPI_SUM, mpi_comm_world,ierror)
    dentmpr0(:,n_modes(1)+1)=dentmp0
    CALL MPI_ALLREDUCE(dentmpi0(:,n_modes(1)+1),dentmp0,mgrid,MPI_RSIZE, MPI_SUM, mpi_comm_world,ierror)
    dentmpi0(:,n_modes(1)+1)=dentmp0

!====================== Zero boundary condition =======================
    if(sfcore == 1)then
!$omp parallel do private(j, ij)
        do j = 1, mpara_core(1)
            ij = igrid_core(1) + j
            dentmpr(ij, :) = 0.0_lk
            dentmpi(ij, :) = 0.0_lk
            dentmpr0(ij, :) = 0.0_lk
            dentmpi0(ij, :) = 0.0_lk
        enddo
      
        if(couplefrc == 0)then
!$omp parallel do private(j, ij)
            do j = 1,mpara_core(mpsi_core)
                ij = igrid_core(mpsi_core) + j
                dentmpr(ij, :) = 0.0_lk
                dentmpi(ij, :) = 0.0_lk
                dentmpr0(ij, :) = 0.0_lk
                dentmpi0(ij, :) = 0.0_lk
            enddo
        endif
    endif
    
    if(sfsol == 1)then
!$omp parallel do private(j, ij)
        do j = 1, mpara_sol(mpsi_sol)
            ij = mgrid_core + igrid_sol(mpsi_sol) + j
            dentmpr(ij, :) = 0.0_lk
            dentmpi(ij, :) = 0.0_lk
            dentmpr0(ij, :) = 0.0_lk
            dentmpi0(ij, :) = 0.0_lk
        enddo   
      
        if(couplefrc == 0)then 
!$omp parallel do private(j, ij)
            do j = 1, mpara_sol(1)
                ij = mgrid_core + igrid_sol(1) + j
                dentmpr(ij, :) = 0.0_lk
                dentmpi(ij, :) = 0.0_lk
                dentmpr0(ij, :) = 0.0_lk
                dentmpi0(ij, :) = 0.0_lk
            enddo
        endif
    endif

!====================  Real part solver  =========================
    luserbfperp = 0.0_lk
    luserxfperp = 0.0_lk
    if(ifullk==0)then
!$omp parallel do private(i)
    do  i = 1, mgrid
        luserbfperp(i-1) = dentmpr(i, n_modes(1)+1)
    enddo
#if PETSC_VERSION_GT(3,8,0)
    call lapmat_pssolver_fperp(petsc_ps, luserbfperp, luserxfperp)
#else
    call lapmat_pssolver_fperp(luserbfperp,luserxfperp)
#endif
!$omp parallel do private(i)
    do i = 1, mgrid
        phitmpr(i, n_modes(1)+1)=luserxfperp(i-1)
    enddo
    !$omp parallel do private(i)
    do  i = 1, mgrid
        luserbfperp(i-1) = dentmpr0(i, 1)
    enddo
#if PETSC_VERSION_GT(3,8,0)
    call lapmat_pssolver_fperp(petsc_ps0, luserbfperp, luserxfperp)
#else
    call lapmat_pssolver_fperp0(luserbfperp,luserxfperp)
#endif
!$omp parallel do private(i)
    do i = 1, mgrid
        phitmpr0(i, 1)=luserxfperp(i-1)
    enddo
    else
!Fully-kinetic
!$omp parallel do private(i)
        do i = 1, mgrid
        phitmpr(i, n_modes(1)+1)=dentmpr(i, n_modes(1)+1)*meshte(i)/meshne(i)
        enddo
    endif
!=================== Imaginary part solver  ======================
    luserbfperp = 0.0_lk
    luserxfperp = 0.0_lk
    if(ifullk==0)then
!$omp parallel do private(i)
    do i = 1, mgrid
        luserbfperp(i-1)=dentmpi(i, n_modes(1)+1)
    enddo
#if PETSC_VERSION_GT(3,8,0) 
    call lapmat_pssolver_fperp(petsc_ps, luserbfperp, luserxfperp)
#else
    call lapmat_pssolver_fperp(luserbfperp,luserxfperp)
#endif
!$omp parallel do private(i)
    do i = 1, mgrid
        phitmpi(i, n_modes(1)+1) = luserxfperp(i-1)
    enddo
!$omp parallel do private(i)
    do i = 1, mgrid
        luserbfperp(i-1)=dentmpi0(i, 1)
    enddo
#if PETSC_VERSION_GT(3,8,0) 
    call lapmat_pssolver_fperp(petsc_ps0, luserbfperp, luserxfperp)
#else
    call lapmat_pssolver_fperp0(luserbfperp,luserxfperp)
#endif
!$omp parallel do private(i)
        do i = 1, mgrid
            phitmpi0(i, 1) = luserxfperp(i-1)
        enddo
    else
!Fully-kinetic
!$omp parallel do private(i)
        do i = 1, mgrid
            phitmpi(i, n_modes(1)+1)=dentmpi(i, n_modes(1)+1)*meshte(i)/meshne(i)
        enddo
    endif
    phi_real=0.0_lk
    phi_imag=0.0_lk
    phi0_real=0.0_lk
    phi0_imag=0.0_lk

    phifluxlz=0.0_lk
    phifluxlz0=0.0_lk
    phi_real=phitmpr(:, n_modes(1)+1)
    phi_imag=phitmpi(:, n_modes(1)+1)
    phi0_real=phitmpr0(:, 1)
    phi0_imag=phitmpi0(:, 1)


eachzphifluxlz_core=0.0_lk
eachzphifluxlz_sol=0.0_lk
eachtphifluxlz_core=0.0_lk
eachtphifluxlz_sol=0.0_lk
eachzphifluxlz0_core=0.0_lk
eachzphifluxlz0_sol=0.0_lk
eachtphifluxlz0_core=0.0_lk
eachtphifluxlz0_sol=0.0_lk

!================== transform back to real space ==================
    if(sfcore == 1)then
!$omp parallel private(xz,yz2)
        yz2 = 0.0_lk
        xz = 0.0_lk
!$omp do private(ij)
        do ij = ij_min_tor_core, ij_max_tor_core
                yz2(n_modes(1) + 1) = cmplx(phitmpr(ij, n_modes(1)+1), phitmpi(ij, n_modes(1)+1))
                yz2 = ffilter*yz2
                call fftr1d(-1, mtoroidal, 1.0_lk, xz, yz2, 2)
                eachzphifluxlz_core(ij-ij_min_core+1,:) = xz
        enddo
        yz2 = 0.0_lk
        xz = 0.0_lk
!$omp do private(ij)
        do ij = ij_min_tor_core, ij_max_tor_core
                yz2(1) = cmplx(phitmpr0(ij, 1), phitmpi0(ij, 1))
                call fftr1d(-1, mtoroidal, 1.0_lk, xz, yz2, 2)
                eachzphifluxlz0_core(ij-ij_min_core+1,:) = xz
        enddo
!$omp end parallel

        CALL MPI_ALLREDUCE(MPI_IN_PLACE,eachzphifluxlz_core,mygrid_core*mtoroidal,MPI_RSIZE, MPI_SUM, toroidal_comm,ierror)
        CALL MPI_ALLREDUCE(MPI_IN_PLACE,eachzphifluxlz0_core,mygrid_core*mtoroidal,MPI_RSIZE, MPI_SUM, toroidal_comm,ierror)
        eachtphifluxlz_core(ij_min_core:ij_max_core)=eachzphifluxlz_core(:,myrank_toroidal+1)
        eachtphifluxlz0_core(ij_min_core:ij_max_core)=eachzphifluxlz0_core(:,myrank_toroidal+1)
        CALL MPI_ALLREDUCE(MPI_IN_PLACE,eachtphifluxlz_core,mgrid_core,MPI_RSIZE, MPI_SUM, partd_comm,ierror)
        CALL MPI_ALLREDUCE(MPI_IN_PLACE,eachtphifluxlz0_core,mgrid_core,MPI_RSIZE, MPI_SUM, partd_comm,ierror)

!$omp parallel do private(i,j,ij) 
        do i=1,mpsi_core
            do j=1,mpara_core(i)
               ij=igrid_core(i)+j
               phi(1,ij)=eachtphifluxlz_core(ij)
               phi0(1,ij)=eachtphifluxlz0_core(ij)
            enddo
        enddo
    endif

phifluxlz=0.0_lk
phifluxlz0=0.0_lk

    if(sfsol==1)then
!$omp parallel private(xz,yz2)
        yz2=0.0_lk
        xz=0.0_lk
!$omp do private(ij)
        do ij = ij_min_tor_sol, ij_max_tor_sol
                yz2(n_modes(1) + 1) = cmplx(phitmpr(ij, n_modes(1)+1), phitmpi(ij, n_modes(1)+1))
                yz2=ffilter*yz2
                call fftr1d(-1,mtoroidal,1.0_lk,xz,yz2,2)
                eachzphifluxlz_sol(ij-ij_min_sol+1,:)=xz
        enddo
        yz2=0.0_lk
        xz=0.0_lk
!$omp do private(ij)
        do ij = ij_min_tor_sol, ij_max_tor_sol
                yz2(1) = cmplx(phitmpr0(ij, 1), phitmpi0(ij, 1))
                call fftr1d(-1,mtoroidal,1.0_lk,xz,yz2,2)
                eachzphifluxlz0_sol(ij-ij_min_sol+1,:)=xz
        enddo
!$omp end parallel
        CALL MPI_ALLREDUCE(MPI_IN_PLACE,eachzphifluxlz_sol,mygrid_sol*mtoroidal,MPI_RSIZE, MPI_SUM, toroidal_comm,ierror)
        CALL MPI_ALLREDUCE(MPI_IN_PLACE,eachzphifluxlz0_sol,mygrid_sol*mtoroidal,MPI_RSIZE, MPI_SUM, toroidal_comm,ierror)
        eachtphifluxlz_sol(ij_min_sol:ij_max_sol)=eachzphifluxlz_sol(:,myrank_toroidal+1)
        eachtphifluxlz0_sol(ij_min_sol:ij_max_sol)=eachzphifluxlz0_sol(:,myrank_toroidal+1)
        CALL MPI_ALLREDUCE(MPI_IN_PLACE,eachtphifluxlz_sol,mgrid_sol,MPI_RSIZE, MPI_SUM, partd_comm,ierror)
        CALL MPI_ALLREDUCE(MPI_IN_PLACE,eachtphifluxlz0_sol,mgrid_sol,MPI_RSIZE, MPI_SUM, partd_comm,ierror)

!$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
                phi(1,ij)=eachtphifluxlz_sol(ij)
                phi0(1,ij)=eachtphifluxlz0_sol(ij)
            enddo
        enddo
    endif

!========================================================================
!   Copy the fields on the grids on separatrix of core region to the grids 
!   on separatrix of SOL region
!========================================================================
    if(couplefrc==1)then
! convert the fields on regular mesh back to the one on irregular mesh on separatrix
        call convert_reg2irreg(phi(1,:))
! evaluate the fields on shared grids in SOL region using the one on shared grids in core region 
        do i = 1, mpara_share
            phi(1, mgrid_core+igrid_sol(1)+mpara_sol_share_index(i)) = phi(1, igrid_core(mpsi_core)+mpara_core_share_index(i))
            phi0(1, mgrid_core+igrid_sol(1)+mpara_sol_share_index(i)) = phi0(1, igrid_core(mpsi_core)+mpara_core_share_index(i))
        enddo
    endif   
!========================================================================
    if(sfcore==1)then
        if(ifullk==0)then
!$omp parallel do private(i, j, ij, ne_ni_ratio, coe_temp) 
        do i = 1, mpsi_core
            do j = 1, mpara_core(i)
                ij = igrid_core(i) + j
                
                ne_ni_ratio=meshne(ij)/meshni(ij)
                if(kinetice==0)then
                coe_temp = 1.0_lk + meshti(ij)/meshte(ij)/qion/qion*ne_ni_ratio
                phi(1,ij) = phi(1,ij) + meshti(ij)*densityi(1,ij)/coe_temp/qion
                phi0(1,ij) = phi0(1,ij) + meshti(ij)*densityi0(1,ij)/qion
                else
                phi(1,ij) = phi(1,ij) + (qion*meshni(ij)*densityi(1,ij)+qelectron*meshne(ij)*densitye(1,ij))&
                            *meshti(ij)/(qion*qion)/meshni(ij)
                phi0(1,ij) = phi0(1,ij) + (qion*meshni(ij)*densityi0(1,ij)+qelectron*meshne(ij)*densitye0(1,ij))&
                            *meshti(ij)/(qion*qion)/meshni(ij)
                endif
            enddo
         enddo
         endif
         phi(1,igrid_core(1) + 1 : igrid_core(1) + mpara_core(1))=0.0_lk
         phi0(1,igrid_core(1) + 1 : igrid_core(1) + mpara_core(1))=0.0_lk
         if(couplefrc==0)then
             phi(1,igrid_core(mpsi_core) + 1 : igrid_core(mpsi_core) + mpara_core(mpsi_core))=0.0_lk
             phi0(1,igrid_core(mpsi_core) + 1 : igrid_core(mpsi_core) + mpara_core(mpsi_core))=0.0_lk
         endif
    endif
    
    if(sfsol==1)then
        if(ifullk==0)then
!$omp parallel do private(i, j, ij, ne_ni_ratio, coe_temp) 
        do i=1,mpsi_sol
            do j=1,mpara_sol(i)
                ij=mgrid_core+igrid_sol(i)+j
                if(kinetice==0)then
                ne_ni_ratio=meshne(ij)/meshni(ij)
                coe_temp = 1.0_lk + meshti(ij)/meshte(ij)/qion/qion*ne_ni_ratio
                phi(1,ij) = phi(1,ij) + meshti(ij)*densityi(1,ij)/coe_temp/qion
                phi0(1,ij) = phi0(1,ij) + meshti(ij)*densityi0(1,ij)/qion
                else
                phi(1,ij) = phi(1,ij) + (qion*meshni(ij)*densityi(1,ij)+qelectron*meshne(ij)*densitye(1,ij))&
                            *meshti(ij)/(qion*qion)/meshni(ij)
                phi0(1,ij) = phi0(1,ij) + (qion*meshni(ij)*densityi0(1,ij)+qelectron*meshne(ij)*densitye0(1,ij))&
                            *meshti(ij)/(qion*qion)/meshni(ij)
                endif
            enddo
        enddo
        endif
        phi(1,mgrid_core + igrid_sol(mpsi_sol) + 1 : mgrid_core + igrid_sol(mpsi_sol) + mpara_sol(mpsi_sol))=0.0_lk
        phi0(1,mgrid_core + igrid_sol(mpsi_sol) + 1 : mgrid_core + igrid_sol(mpsi_sol) + mpara_sol(mpsi_sol))=0.0_lk
        if(couplefrc==0)then
            phi(1,mgrid_core + igrid_sol(1) + 1 : mgrid_core + igrid_sol(1) + mpara_sol(1))=0.0_lk
            phi0(1,mgrid_core + igrid_sol(1) + 1 : mgrid_core + igrid_sol(1) + mpara_sol(1))=0.0_lk
        endif
    endif


!====================== Zero boundary condition for phi =======================
    if(sfcore == 1)then
!$omp parallel do private(j, ij)
        do j = 1, mpara_core(1)
            ij = igrid_core(1) + j
            phi_real(ij) = 0.0_lk
            phi_imag(ij) = 0.0_lk
        enddo

        if(couplefrc == 0)then
!$omp parallel do private(j, ij)
            do j = 1,mpara_core(mpsi_core)
                ij = igrid_core(mpsi_core) + j
                phi_real(ij) = 0.0_lk
                phi_imag(ij) = 0.0_lk
            enddo
        endif
    endif

    if(sfsol == 1)then
!$omp parallel do private(j, ij)
        do j = 1, mpara_sol(mpsi_sol)
            ij = mgrid_core + igrid_sol(mpsi_sol) + j
            phi_real(ij) = 0.0_lk
            phi_imag(ij) = 0.0_lk
        enddo

        if(couplefrc == 0)then
!$omp parallel do private(j, ij)
            do j = 1, mpara_sol(1)
                ij = mgrid_core + igrid_sol(1) + j
                phi_real(ij) = 0.0_lk
                phi_imag(ij) = 0.0_lk
            enddo
        endif
    endif
!========================================================================


    call SMOOTH(phi)
    call SMOOTH(phi0)
    call toroidal_filter(phi)

!smooth over parallel field line
adum=0.0_lk
if(sfcore==1)then
  do i=1,mpsi_core
     do j=1,mpara_core(i)
        ij=igrid_core(i)+j
        adum(:,i)=adum(:,i)+phi0(:,ij)/mpara_core(i)
     enddo
     do j=1,mpara_core(i)
        ij=igrid_core(i)+j
        phi0(:,ij)=adum(:,i)
     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
        adum(:,mpsi_core+i)=adum(:,mpsi_core+i)+phi0(:,ij)/mpara_sol(i)
     enddo
     do j=1,mpara_sol(i)
        ij=mgrid_core+igrid_sol(i)+j
        phi0(:,ij)=adum(:,mpsi_core+i)
     enddo
  enddo
endif

phi_rms=0.0_lk
phi0_rms=0.0_lk
total_v = 0.0_lk
data1dtmp = 0.0_lk
total_v1d = 0.0_lk
if(sfcore==1)then
!$omp parallel do private(ij) reduction(+:data1dtmp,total_v1d,phi_rms,phi0_rms,total_v)
  do i=1,mpsi_core
     do j=1,mpara_core(i)
        ij=igrid_core(i)+j
        data1dtmp(i,1) = data1dtmp(i,1)+phi(1,ij)**2*volume(ij)
        data1dtmp(i,2) = data1dtmp(i,2)+phi0(1,ij)*volume(ij)
        total_v1d(i) = total_v1d(i) + volume(ij)
     enddo
     phi_rms = phi_rms + data1dtmp(i,1)
     total_v = total_v + total_v1d(i)
  enddo
endif
if(sfsol==1)then
!$omp parallel do private(ij) reduction(+:data1dtmp,total_v1d,phi_rms,phi0_rms,total_v)
  do i=1,mpsi_sol
     do j=1,mpara_sol(i)
        ij=mgrid_core+igrid_sol(i)+j
        data1dtmp(mpsi_core+i,1) = data1dtmp(mpsi_core+i,1)+phi(1,ij)**2*volume(ij)
        data1dtmp(mpsi_core+i,2) = data1dtmp(mpsi_core+i,2)+phi0(1,ij)*volume(ij)
        total_v1d(mpsi_core+i) = total_v1d(mpsi_core+i) + volume(ij)
     enddo
     phi_rms = phi_rms + data1dtmp(mpsi_core+i,1)
     total_v = total_v + total_v1d(mpsi_core+i)
  enddo
endif
 !$omp parallel do private(ij) reduction(+:phi_rms,phi0_rms,total_v)
 do ij=1,mgrid
    phi0_rms = phi0_rms + phi0(1,ij)**2*volume(ij)
 enddo

phi_rms = sqrt(phi_rms/total_v)
phi0_rms = sqrt(phi0_rms/total_v)
call MPI_ALLREDUCE(phi_rms,total_v,1,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
total_v = total_v / mtoroidal
phi_rms = total_v
call MPI_ALLREDUCE(phi0_rms,total_v,1,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
total_v = total_v / mtoroidal
phi0_rms = total_v

data1dtmp(:,1) = sqrt(data1dtmp(:,1)/total_v1d)
data1dtmp(:,2) = data1dtmp(:,2)/total_v1d
call MPI_ALLREDUCE(data1dtmp(:,1),total_v1d,mpsi,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
total_v1d = total_v1d / mtoroidal
data1dtmp(:,1) = total_v1d
call MPI_ALLREDUCE(data1dtmp(:,2),total_v1d,mpsi,mpi_Rsize,MPI_SUM,toroidal_comm,ierror)
total_v1d = total_v1d / mtoroidal
data1dtmp(:,2) = total_v1d

if(ifullk==0)then
data1di(:,4:5) = data1dtmp
else
data1dfki(:,4:5) = data1dtmp
endif

end subroutine global_pssolver

!=====================================================================
!   Not work, to be fixed
!=====================================================================
subroutine semi_spectrum_r2c(grid_poloidal_plane, grid_psi, grid_para_array, grid_index_array, farray_real, farray_spectrum_real, farray_spectrum_imag)
    use precision
    use read_in_parameters, only: n_modes, nfilter
    use global_parameters, only: toroidal_comm, mtoroidal, mype, gtcout
    implicit none

    integer, intent(in) :: grid_poloidal_plane, grid_psi
    integer, dimension(grid_psi), intent(in) :: grid_para_array, grid_index_array
    real(lk), dimension(grid_poloidal_plane), intent(in) :: farray_real

    real(lk), dimension(grid_poloidal_plane, mtoroidal/2+1), intent(out) :: farray_spectrum_real, farray_spectrum_imag

    integer icount, i, j, ij, ierror
    real(lk) :: farraytmp_1d(grid_poloidal_plane), farraytmp_2d(grid_poloidal_plane, mtoroidal), ffilter(mtoroidal/2+1)
    real(lk) :: xz(mtoroidal)
    complex(lk) :: yz(mtoroidal/2+1)

    farraytmp_1d = 0.0_lk
    farraytmp_2d = 0.0_lk

!$omp parallel do private(i)
    do i = 1, grid_poloidal_plane
        farraytmp_1d(i) = farray_real(i)
    enddo

    icount = grid_poloidal_plane
    call MPI_ALLGATHER(farraytmp_1d, icount, mpi_Rsize, farraytmp_2d, icount, &
                       mpi_Rsize, toroidal_comm, ierror)

    ffilter = 0.0_lk
    if(nfilter == 0)then
        ffilter = 1.0_lk

    elseif(nfilter == 1)then
        ffilter(n_modes+1) = 1.0_lk

    elseif(nfilter == 2)then
        ffilter(n_modes(1)+1) = 1.0_lk

    else
        if(mype==0)write(gtcout,*)"wrong nfilter value"
        stop
    endif


    yz = 0.0_lk
!$omp parallel do private(i, j, ij, xz, yz)
    do i = 1, grid_psi
        do j = 1, grid_para_array(i)
            ij = grid_index_array(i) + j
            xz = farraytmp_2d(ij,:)
            call fftr1d(1, mtoroidal, 1.0_lk, xz, yz, 2)
! linear run only keep a few modes
            yz = ffilter*yz
            farray_spectrum_real(ij, :) = real(real(yz(:)))
            farray_spectrum_imag(ij, :) = real(aimag(yz(:)))
        enddo
    enddo

end subroutine semi_spectrum_r2c
!=====================================================================
!   Not work, to be fixed
!=====================================================================
subroutine semi_spectrum_c2r(grid_poloidal_plane, grid_psi, grid_para_array, grid_index_array, farray_real, farray_spectrum_real, farray_spectrum_imag)
    use precision
    use read_in_parameters, only: n_modes
    use global_parameters, only: toroidal_comm, mtoroidal, myrank_toroidal
    implicit none

    integer, intent(in) :: grid_poloidal_plane, grid_psi
    integer, dimension(grid_psi), intent(in) :: grid_para_array, grid_index_array

    real(lk), dimension(grid_poloidal_plane, mtoroidal/2+1), intent(in) :: farray_spectrum_real, farray_spectrum_imag
    real(lk), dimension(grid_poloidal_plane), intent(out) :: farray_real

    integer  i, j, ij, k
    real(lk) :: farraytmp_2d(grid_poloidal_plane, mtoroidal), ffilter(mtoroidal/2+1)
    real(lk) :: xz(grid_poloidal_plane, mtoroidal)
    complex(lk) :: farraytmp_spectrum_2d(grid_poloidal_plane, mtoroidal/2+1)

!=========================================================================
    farraytmp_2d = 0.0_lk

!$omp parallel do private(i, j)
    do i = 1, grid_poloidal_plane
        do j = 1, mtoroidal/2 + 1
            farraytmp_spectrum_2d(i, j) = cmplx(farray_spectrum_real(i, j), farray_spectrum_imag(i, j))
        enddo
    enddo

    xz = 0.0_lk
!$omp parallel do private(i, j, ij)
    do i = 1, grid_psi
        do j = 1, grid_para_array(i)
            ij = grid_index_array(i) + j
            do k = 1, mtoroidal/2 + 1
                call fftr1d(-1, mtoroidal, 1.0_lk, xz(ij, :), farraytmp_spectrum_2d(ij, k), 2)
                farraytmp_2d(ij, :) = farraytmp_2d(ij, :) + xz(ij, :)
            enddo
        enddo
    enddo

!$omp parallel do private(i, j, ij) 
    do i = 1, grid_psi
        do j = 1, grid_para_array(i)
           ij = grid_index_array(i) + j
           farray_real(ij) = farraytmp_2d(ij, myrank_toroidal + 1)
        enddo
    enddo

end subroutine semi_spectrum_c2r
