module grids_fields
use precision
implicit none

public frc_mesh_perturb_field_initial, gradient, periodicity, toroidal_filter, smooth, &
! common fields       
       phi,phi0,phi0_real,phi0_imag, phi_real, phi_imag, gradphi,gradphi0, & 
       phi_r_zeta_core, phi_r_zeta_sol, sol_parallel_filter, phi_rms, phi0_rms
!private sol_parallel_filter


! common fields
real(lk),dimension(:),allocatable ::  phi_real, phi_imag,phi0_real, phi0_imag, &
                                      phi_r_zeta_core, phi_r_zeta_sol
real(lk),dimension(:,:),allocatable :: phi, phi0
real(lk),dimension(:,:,:),allocatable :: gradphi, gradphi0
real(lk) :: phi_rms, phi0_rms
!$acc declare create(gradphi, gradphi0)

save

contains

!===============================================================
subroutine frc_mesh_perturb_field_initial
    use grids_setup, only: mgrid, mpsi_core, mpsi_sol, mpsi
    use global_parameters, only: mtoroidal
 
! allocate memory for common fields
    allocate(phi_real(mgrid),phi_imag(mgrid),phi(0:1,mgrid),gradphi(3,0:1,mgrid),&
             phi0(0:1,mgrid),phi0_real(mgrid),phi0_imag(mgrid),gradphi0(3,0:1,mgrid))
   
    phi=0.0_lk
    phi_real=0.0_lk
    phi_imag=0.0_lk
    phi0=0.0_lk
    phi0_real=0.0_lk
    phi0_imag=0.0_lk
    gradphi=0.0_lk
    gradphi0=0.0_lk
! diagnositic array

    allocate(phi_r_zeta_core(mpsi_core*mtoroidal), phi_r_zeta_sol(mpsi_sol*mtoroidal))
    phi_r_zeta_core=0.0_lk
    phi_r_zeta_sol=0.0_lk
end subroutine frc_mesh_perturb_field_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine gradient(scalar,vector)
    use global_parameters
    use read_in_parameters
    use grids_setup
    use simulation_domain, only: sfcore, sfsol
  
    integer :: i,ii,ij,j,k,icount,idest,isource,isendtag,irecvtag,ierror,jt,&
        istatus(MPI_STATUS_SIZE),ijsol,ij1,ijn
    real(lk) :: diffz,pleft(mgrid),pright(mgrid),&
        sendl(mgrid),recvr(mgrid),sendr(mgrid),recvl(mgrid)
    real(lk),dimension(mgrid),intent(in) :: scalar
    real(lk),dimension(3,0:1,mgrid),intent(out) :: vector
    real(lk) :: h0,h1,wp0,wp1,wt10,wt00,wt11,wt01,sdum0,sdum1,pdum0,pdum1,&
             diff_s1,diff_s2,diff_s,weight1,weight2,vectmp1,vectmp2,&
             x0x1,x0x2,x1x2,L0prime,L1prime,L2prime
    real(lk) :: dp0,dp1
    integer :: ijlow,ijhigh
  
! finite difference for e-field in equilibrium unit
  
    diffz=0.5_lk/deltaz
  
!$omp parallel do private(i,j,k)
    do i=1,mgrid
       do j=0,1
          do k=1,3
            vector(k,j,i)=0.0_lk
          enddo
       enddo
    enddo
  
! d_scalar/d_para
    if(sfcore==1)then
!$omp parallel do private(i,j,ij,diff_s1,diff_s2,vectmp1,vectmp2,weight1,weight2, ij1,ijn,diff_s)
        do i=1,mpsi_core
            do j=2,mpara_core(i)-1
                ij=igrid_core(i)+j
    
                diff_s1=1.0_lk/deltas_core1d(i)
                diff_s2=1.0_lk/deltas_core1d(i)

                vectmp1=diff_s1*(scalar(ij)-scalar(ij-1))
                vectmp2=diff_s2*(scalar(ij+1)-scalar(ij))
    
                weight1=0.5_lk
                weight2=0.5_lk
    
                vector(3,1,ij)=weight1*vectmp1+weight2*vectmp2
            enddo
            ij1=igrid_core(i)+1
            ijn=igrid_core(i)+mpara_core(i)
       
            diff_s=1.0_lk/(2.0_lk*deltas_core1d(i))

            vector(3,1,ij1)=diff_s*(scalar(ij1+1)-scalar(ijn-1))
            vector(3,1,ijn)=diff_s*(scalar(ij1+1)-scalar(ijn-1))
        enddo
    endif
  
    if(sfsol==1)then
!$omp parallel do private(i,j,ij,ijsol,diff_s1,diff_s2,vectmp1,vectmp2,weight1,weight2)
        do i=1,mpsi_sol
            do j=2,mpara_sol(i)-1
                ij=mgrid_core+igrid_sol(i)+j
                ijsol=igrid_sol(i)+j
  
                diff_s1=1.0_lk/deltas_sol(ijsol-1)
                diff_s2=1.0_lk/deltas_sol(ijsol)
  
                vectmp1=diff_s1*(scalar(ij)-scalar(ij-1))
                vectmp2=diff_s2*(scalar(ij+1)-scalar(ij))
  
                weight1=deltas_sol(ijsol)/(deltas_sol(ijsol)+deltas_sol(ijsol-1))
                weight2=deltas_sol(ijsol-1)/(deltas_sol(ijsol)+deltas_sol(ijsol-1))
  
                vector(3,1,ij)=weight1*vectmp1+weight2*vectmp2
            enddo
        enddo
  
        if(sol_para_periodicity==0)then ! reflected field gradient boundary condition
            do i=1,mpsi_sol
                ij1=mgrid_core+igrid_sol(i)+1
                x0x1=-deltas_sol(ij1-mgrid_core) !x0-x1
                x0x2=-deltas_sol(ij1-mgrid_core)-deltas_sol(ij1-mgrid_core+1) !x0-x2
                x1x2=-deltas_sol(ij1-mgrid_core+1) !x1-x2
                L0prime=(x0x1+x0x2)/(x0x1*x0x2)
                L1prime=x0x2/(-x0x1*x1x2)
                L2prime=x0x1/((-x0x2)*(-x1x2))
                vector(3,1,ij1)=L0prime*scalar(ij1)+L1prime*scalar(ij1+1)+L2prime*scalar(ij1+2)
  
                ijn=mgrid_core+igrid_sol(i)+mpara_sol(i)
                x0x1=deltas_sol(ijn-mgrid_core-1) !x0-x1
                x0x2=deltas_sol(ijn-mgrid_core-1)+deltas_sol(ijn-mgrid_core-2) !x0-x2
                x1x2=deltas_sol(ijn-mgrid_core-2) !x1-x2
                L0prime=(x0x1+x0x2)/(x0x1*x0x2)
                L1prime=x0x2/(-x0x1*x1x2)
                L2prime=x0x1/((-x0x2)*(-x1x2))
                vector(3,1,ijn)=L0prime*scalar(ijn)+L1prime*scalar(ijn-1)+L2prime*scalar(ijn-2)
            enddo
        elseif(sol_para_periodicity==1)then ! periodic field gradient boundary condition
            do i=1, mpsi_sol
                ij1=mgrid_core+igrid_sol(i)+1
                ijn=mgrid_core+igrid_sol(i)+mpara_sol(i)
     
                diff_s=1.0_lk/(deltas_sol(ij1-mgrid_core)+deltas_sol(ijn-mgrid_core-1))
!                vector(3,1,ij1)=diff_s*(scalar(ij1+1)-scalar(ij1)+scalar(ijn)-scalar(ijn-1))
                vector(3,1,ij1)=diff_s*(scalar(ij1+1)-scalar(ijn-1))
                vector(3,1,ijn)=vector(3,1,ij1)
            enddo
        endif
    endif
  
! send scalar to right and receive from left
    sendr=scalar(:)
    recvl=0.0_lk
    icount=mgrid
    idest=right_pe
    isource=left_pe
    isendtag=myrank_toroidal
    irecvtag=isource
    call MPI_SENDRECV(sendr,icount,mpi_Rsize,idest,isendtag,&
        recvl,icount,mpi_Rsize,isource,irecvtag,toroidal_comm,istatus,ierror)
  
! send scalar to left and receive from right
    sendl=scalar(:)
    recvr=0.0_lk
    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)
! unpack scalar_boundary and calculate E_zeta at boundaries toroidal
  
    pleft=recvl
    pright=recvr

    if(sfcore==1)then
!$omp parallel do private(i,j,ij)
        do i=1,mpsi_core
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                vector(2,1,ij)=(pright(ij)-pleft(ij))*diffz
            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
                vector(2,1,ij)=(pright(ij)-pleft(ij))*diffz
            enddo
        enddo
    endif

    if(sfcore==1)then
!$omp parallel do private(i,j,ij,wp0,wp1,dp0,dp1,ijlow,ijhigh)
       do i=2,mpsi_core-1
          do j=1,mpara_core(i)
             ij=igrid_core(i)+j
             ijlow=igrid_core(i-1)+j
             ijhigh=igrid_core(i+1)+j
             dp1=psimesh_core(i+1)-psimesh_core(i)
             dp0=psimesh_core(i)-psimesh_core(i-1)
             wp0=dp0/(dp0+dp1)
             wp1=1.0_lk-wp0
             vector(1,1,ij)=(scalar(ijhigh)-scalar(ij))/dp1*wp0+(scalar(ij)-scalar(ijlow))/dp0*wp1
          enddo
       enddo
    endif
    if(sfsol==1)then
!$omp parallel do private(i,j,ij,wp0,wp1,dp0,dp1,ijlow,ijhigh)
       do i=2,mpsi_sol-1
          do j=1,mpara_sol(i)
             ij=mgrid_core+igrid_sol(i)+j
             ijlow=mgrid_core+igrid_sol(i-1)+j
             ijhigh=mgrid_core+igrid_sol(i+1)+j
             dp1=psimesh_sol(i+1)-psimesh_sol(i)
             dp0=psimesh_sol(i)-psimesh_sol(i-1)
             wp0=dp0/(dp0+dp1)
             wp1=1.0_lk-wp0
             vector(1,1,ij)=(scalar(ijhigh)-scalar(ij))/dp1*wp0+(scalar(ij)-scalar(ijlow))/dp0*wp1
          enddo
       enddo

! boundary condition for gradient in psi direction
       do j=1,mpara_sol(1)
           ij1=mgrid_core+igrid_sol(1)+j
           x0x1=psimesh_sol(1)-psimesh_sol(2) !x0-x1
           x0x2=psimesh_sol(1)-psimesh_sol(3) !x0-x2
           x1x2=psimesh_sol(2)-psimesh_sol(3) !x1-x2
           L0prime=(x0x1+x0x2)/(x0x1*x0x2)
           L1prime=x0x2/(-x0x1*x1x2)
           L2prime=x0x1/((-x0x2)*(-x1x2))
           vector(1,1,ij1)=L0prime*scalar(ij1)+L1prime*scalar(ij1-igrid_sol(1)+igrid_sol(2))+L2prime*scalar(ij1-igrid_sol(1)+igrid_sol(3))
   
           ijn=mgrid_core+igrid_sol(mpsi_sol)+j
           x0x1=psimesh_sol(mpsi_sol-2)-psimesh_sol(mpsi_sol-1) !x0-x1
           x0x2=psimesh_sol(mpsi_sol-2)-psimesh_sol(mpsi_sol) !x0-x2
           x1x2=psimesh_sol(mpsi_sol-1)-psimesh_sol(mpsi_sol) !x1-x2
           L0prime=(x0x1+x0x2)/(x0x1*x0x2)
           L1prime=x0x2/(-x0x1*x1x2)
           L2prime=x0x1/((-x0x2)*(-x1x2))
           vector(1,1,ijn)=L0prime*scalar(ijn-igrid_sol(mpsi_sol)+igrid_sol(mpsi_sol-2)) +L1prime*scalar(ijn-igrid_sol(mpsi_sol)+igrid_sol(mpsi_sol-1)) +L2prime*scalar(ijn)
       enddo
    endif
  
    call periodicity(vector(1,:,:))
    call periodicity(vector(2,:,:))
    call periodicity(vector(3,:,:))
  
  
end subroutine gradient
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine periodicity(fgrid)
    use global_parameters
    use grids_setup
  
    integer :: i,ii,jt,icount,idest,isource,isendtag,irecvtag,ierror,istatus(MPI_STATUS_SIZE)
    real(lk) :: sendr(mgrid),recvl(mgrid)
    real(lk),dimension(0:1,mgrid),intent(inout) :: fgrid
  
! toroidal end point, send E to right and receive from left
!$omp parallel do private(i)
    do i=1,mgrid
       sendr(i)=fgrid(1,i)
       recvl(i)=0.0_lk
    enddo
    icount=mgrid
    idest=right_pe
    isource=left_pe
    isendtag=myrank_toroidal
    irecvtag=isource
    call MPI_SENDRECV(sendr,icount,mpi_Rsize,idest,isendtag,&
        recvl,icount,mpi_Rsize,isource,irecvtag,toroidal_comm,istatus,ierror)
!$omp parallel do private(i)
    do i=1,mgrid
        fgrid(0,i)=recvl(i)
    enddo
  
end subroutine periodicity
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine toroidal_filter(farray)
    use global_parameters
    use grids_setup
    use simulation_domain, only: sfcore, sfsol
    use read_in_parameters, only: n_modes, sol_pfilter
  
    integer i,j,k,ii,ij,ip,jt,kz,jpe,indp,indt,indp1,indt1,meachtheta,mteach,mzeach,jtp,&
         icount,ierror,idest,isource,isendtag,l,irecvtag,joffset,mteachpd,&
         istatus(MPI_STATUS_SIZE)
    real(lk) ::  wt,r,dt,wz,zdum,tdum,ffilter(mtoroidal/2+1),&
         xz(mtoroidal),&
         eachzphifluxlz_core(mygrid_core,mtoroidal),&
         eachtphifluxlz_core(mygrid_core,mtoroidal),&
         eachtphifluxlz_sol(mygrid_sol,mtoroidal),&
         eachzphifluxlz_sol(mygrid_sol,mtoroidal),&
         eachzeta_core(mygrid_core),eachzeta_sol(mygrid_sol)
    real(lk),dimension(0:1,mgrid),intent(inout) :: farray
    real(lk) :: farray_tmp(mgrid)
    complex(lk) yz(mtoroidal/2+1)
  
! periodic BC in poloidal and toroidal
!    CALL PERIODICITY(farray)

    eachzeta_core=0.0_lk
    eachzeta_sol=0.0_lk
    eachtphifluxlz_core=0.0_lk
    eachtphifluxlz_sol=0.0_lk
    icount=mgrid
    if(sfcore==1)then
!$omp parallel do private(ij) 
        do ij=ij_min_core,ij_max_core
                eachzeta_core(ij-ij_min_core+1)=farray(1,ij)
        enddo
    endif
    call MPI_ALLGATHER(eachzeta_core,mygrid_core,mpi_Rsize,eachzphifluxlz_core,mygrid_core,&
         mpi_Rsize,toroidal_comm,ierror)
    if(sfsol==1)then
!$omp parallel do private(ij)
        do ij=ij_min_sol,ij_max_sol
                eachzeta_sol(ij-ij_min_sol+1)=farray(1,ij)
        enddo
    endif
    call MPI_ALLGATHER(eachzeta_sol,mygrid_sol,mpi_Rsize,eachzphifluxlz_sol,mygrid_sol,&
            mpi_Rsize,toroidal_comm,ierror)
    farray_tmp=0.0_lk
    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=eachzphifluxlz_core(ij-ij_min_core+1,:)
                call fftr1d(1,mtoroidal,1.0_lk,xz,yz,2)
! linear run only keep a few modes
                yz=ffilter*yz
! transform back to real space
                call fftr1d(-1,mtoroidal,1.0_lk,xz,yz,2)
                eachtphifluxlz_core(ij-ij_min_core+1,:)=xz
        enddo
        eachzphifluxlz_core=0.0_lk
        call mpi_allreduce(eachtphifluxlz_core,eachzphifluxlz_core,mygrid_core*mtoroidal,MPI_RSIZE, MPI_SUM,toroidal_comm,ierror)
!$omp parallel do private(ij)
        do ij=ij_min_core,ij_max_core
                farray_tmp(ij)=eachzphifluxlz_core(ij-ij_min_core+1,myrank_toroidal+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=eachzphifluxlz_sol(ij-ij_min_sol+1,:)
                call fftr1d(1,mtoroidal,1.0_lk,xz,yz,2)
! linear run only keep a few modes
                yz=ffilter*yz
! transform back to real space
                call fftr1d(-1,mtoroidal,1.0_lk,xz,yz,2)
                eachtphifluxlz_sol(ij-ij_min_sol+1,:)=xz
        enddo
        eachzphifluxlz_sol=0.0_lk
        call mpi_allreduce(eachtphifluxlz_sol,eachzphifluxlz_sol,mygrid_sol*mtoroidal,MPI_RSIZE, MPI_SUM,toroidal_comm,ierror)
!$omp parallel do private(ij)
        do ij=ij_min_sol,ij_max_sol
                farray_tmp(ij)=eachzphifluxlz_sol(ij-ij_min_sol+1,myrank_toroidal+1)
        enddo
    endif
    call mpi_ALLREDUCE(MPI_IN_PLACE,farray_tmp,mgrid,MPI_RSIZE, MPI_SUM, partd_comm,ierror)
    farray(1,:)=farray_tmp
    call periodicity(farray)
  
    if(sol_pfilter==1)call sol_parallel_filter(farray)
  
end subroutine toroidal_filter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine smooth(farray)
    use global_parameters
    use grids_setup
    use simulation_domain, only: sfcore, sfsol
    use read_in_parameters, only: BC_type
    integer :: i,j,ij,ijsol,ij1,ijn,icount,idest,isource,isendtag,irecvtag,ierror,jt,&
        istatus(MPI_STATUS_SIZE)
    real(lk) :: vectmp1,vectmp2,weight1,weight2,farray(0:1,mgrid),farraytmp(mgrid)
    real(lk) :: pleft(mgrid),pright(mgrid),&
        sendl(mgrid),recvr(mgrid),sendr(mgrid),recvl(mgrid)
  
!$omp parallel do private(i)
    do i=1,mgrid
       farraytmp(i)=farray(1,i)
    enddo
  
!Perpendicular smooth
! send scalar to right and receive from left
    sendr=farraytmp(:)
    recvl=0.0_lk
    icount=mgrid
    idest=right_pe
    isource=left_pe
    isendtag=myrank_toroidal
    irecvtag=isource
    call MPI_SENDRECV(sendr,icount,mpi_Rsize,idest,isendtag,&
        recvl,icount,mpi_Rsize,isource,irecvtag,toroidal_comm,istatus,ierror)
  
! send scalar to left and receive from right
    sendl=farraytmp(:)
    recvr=0.0_lk
    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)
! unpack scalar_boundary and calculate E_zeta at boundaries toroidal
  
    pleft=recvl
    pright=recvr
  
    if (sfcore==1) then
!$omp parallel do private(i,j,ij)
        do i = 1,mpsi_core
            do j = 1,mpara_core(i)
                ij = igrid_core(i) + j
                farraytmp(ij) = (0.25_lk*pright(ij)+0.5_lk*farraytmp(ij)+0.25_lk*pleft(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
                farraytmp(ij)=(0.25_lk*pright(ij)+0.5_lk*farraytmp(ij)+0.25_lk*pleft(ij))
            enddo
        enddo
    endif
  
!Parallel smooth
    if (sfcore==1) then
!$omp parallel do private(i,j,ij,ij1,ijn,vectmp1,vectmp2,weight1,weight2)
        do i = 1, mpsi_core
            do j = 2, mpara_core(i) - 1
                ij = igrid_core(i) + j

                vectmp1 = 0.5_lk*(farraytmp(ij - 1) + farraytmp(ij))
                vectmp2 = 0.5_lk*(farraytmp(ij) + farraytmp(ij + 1))

                weight1 = 0.5_lk
                weight2 = 0.5_lk

                farray(1, ij) = weight1*vectmp1 + weight2*vectmp2
            enddo
            ij1 = igrid_core(i) + 1
            ijn = igrid_core(i) + mpara_core(i)
            vectmp1 = 0.5_lk*(farraytmp(ij1) + farraytmp(ij1 + 1))
            vectmp2 = 0.5_lk*(farraytmp(ijn - 1) + farraytmp(ijn))

            weight1 = 0.5_lk
            weight2 = 0.5_lk

            farray(1, ij1) = weight1*vectmp1 + weight2*vectmp2
            farray(1, ijn) = farray(1, ij1)
        enddo  
    endif

    if (sfsol==1) then
!$omp parallel do private(i,j,ij,ijsol,ij1,ijn,vectmp1,vectmp2,weight1,weight2)
        do i=1,mpsi_sol
            do j=2,mpara_sol(i)-1
                ij=mgrid_core+igrid_sol(i)+j
                ijsol=igrid_sol(i)+j
  
                vectmp1=0.5_lk*(farraytmp(ij)+farraytmp(ij-1))
                vectmp2=0.5_lk*(farraytmp(ij+1)+farraytmp(ij))
  
                weight1=deltas_sol(ijsol)/(deltas_sol(ijsol)+deltas_sol(ijsol-1))
                weight2=deltas_sol(ijsol-1)/(deltas_sol(ijsol)+deltas_sol(ijsol-1))
  
                farray(1,ij)=weight1*vectmp1+weight2*vectmp2
            enddo
            ij1=mgrid_core+igrid_sol(i)+1
            ijn=mgrid_core+igrid_sol(i)+mpara_sol(i)   
            vectmp1=0.5_lk*(farraytmp(ij1)+farraytmp(ij1+1))
            vectmp2=0.5_lk*(farraytmp(ijn-1)+farraytmp(ijn))
         
            weight1=deltas_sol(ijn-1-mgrid_core)/(deltas_sol(ijn-1-mgrid_core)+deltas_sol(ij1-mgrid_core))
            weight2=deltas_sol(ij1-mgrid_core)/(deltas_sol(ijn-1-mgrid_core)+deltas_sol(ij1-mgrid_core))
   
            farray(1,ij1)=weight1*vectmp1+weight2*vectmp2
            farray(1,ijn)=farray(1,ij1)
        enddo

        if(BC_type==1)call radial_decay_2d(farray(1,:))
    endif
    call periodicity(farray)
  
end subroutine smooth
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine sol_parallel_filter(farray)
    use global_parameters
    use grids_setup
    use read_in_parameters, only: m_modes
    integer i,j,ij
    real(lk) wt,r,dt,wz,zdum,tdum,ffilter(mpara_sol(1)/2+1),&
         xz(mpara_sol(1)),farray(0:1,mgrid),&
         eachr(mpsi_sol,mpara_sol(1))
    complex(lk) yz(mpara_sol(1)/2+1)
  
!$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
            eachr(i,j)=farray(1,ij)
        enddo
    enddo
  
    yz=0.0_lk
    ffilter=0.0_lk
    ffilter(m_modes+1)=1.0_lk
!$omp parallel do private(i,xz,yz)
    do i=1,mpsi_sol
        xz=eachr(i,:)
        call fftr1d(1,mpara_sol(1),1.0_lk,xz,yz,2)
! linear run only keep a few modes
        yz=ffilter*yz
! transform back to real space
        call fftr1d(-1,mpara_sol(1),1.0_lk,xz,yz,2)
        eachr(i,:)=xz
    enddo
  
!$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
            farray(1,ij)=eachr(i,j)
        enddo
    enddo
  
end subroutine sol_parallel_filter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module grids_fields
