module laplacian
use precision, only: lk
implicit none

public lapperp_initial, laplacian_perp,&
       mlapperp, nlapperp, indexplapperp, lapmatperp, lapmatperp0
private

integer mlapperp
integer, dimension(:), allocatable :: nlapperp
integer, dimension(:,:), allocatable :: indexplapperp
real(lk), dimension(:,:), allocatable :: lapmatperp, lapmatperp0

save

contains
!====================================================================
subroutine lapperp_initial
    use global_parameters
    use grids_setup
    use simulation_domain, only: sfcore, sfsol, couplefrc
#ifdef _PETSc
#include "petscversion.h"
#if PETSC_VERSION_GT(3,8,0)
    use petsc_gtc, only: lapmat_pscreate_fperp, petsc_ps, petsc_ps0
#else
    use petsc_gtc, only: lapmat_pscreate_fperp, lapmat_pscreate_fperp0
#endif
    use petsc_array, only: luserafperp,luserafperp0, luserpfperp, lusersfperp
#endif

    integer :: ierror,mtest,i,j


    mlapperp=9
!!allocate array
    allocate(nlapperp(mgrid),indexplapperp(mlapperp,mgrid),&
        lapmatperp(mlapperp,mgrid),lapmatperp0(mlapperp,mgrid),stat=mtest)

    if (mtest /= 0) then
        if(mype==0)write(0,*)mype,'lapperp_initial: Cannot allocate lapmat'
        call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif

    lapmatperp=0.0_lk
    lapmatperp0=0.0_lk
    nlapperp=0
    indexplapperp=0

!====================================================================
    if(sfsol==1)call lapsol_initial
    if(sfcore==1)call lapcore_initial
    if(couplefrc==1)call lapcoup_initial
!  call lapdebug_initial

! lapmat coefficient index Checking
    if(mype==0)then
        do i=1,mgrid
            do j=1,nlapperp(i)
                if(indexplapperp(j,i)==0)then
                    write(gtcout,*)"warning: indexplapperp is not defined!!! at","mgrid=",i,"nlapperp=",j
                    stop
                endif
            enddo
        enddo
    endif

#ifdef _PETSc
    call lap2petsc_perp
#if PETSC_VERSION_GT(3,8,0)
    call lapmat_pscreate_fperp(petsc_ps, mlapperp,mgrid, luserafperp, luserpfperp, lusersfperp)
    call lapmat_pscreate_fperp(petsc_ps0, mlapperp, mgrid, luserafperp0, luserpfperp, lusersfperp)
#else
    call lapmat_pscreate_fperp(mlapperp,mgrid, luserafperp, luserpfperp, lusersfperp)
    call lapmat_pscreate_fperp0(mlapperp,mgrid, luserafperp0, luserpfperp, lusersfperp)
#endif
#endif

    if(mype==0)write(gtcout,*)"lapperp_initial is used"

end subroutine lapperp_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef _PETSc
subroutine lap2petsc_perp
!*Using nindex0,indexp0,ring0 generate matrices users, usera, and userp
    use precision
    use global_parameters
    use read_in_parameters, only: aion, qion, kinetice
    use simulation_domain, only: sfcore, sfsol, couplefrc
    use grids_setup
    use grids_eqdata
    use petsc_array,only: luserafperp,luserafperp0,lusersfperp,luserpfperp,luserbfperp,luserxfperp

    integer nnz,windexp,msort,nindex0max
    integer,dimension(mgrid) :: nindex0
    integer,dimension(mlapperp,mgrid) :: indexp0
    real(lk),dimension(mlapperp,mgrid) :: ring0, ring00

    integer :: ierror,mtest
    real(lk) :: wring,wring0,rhos2,coe_temp,dielectric_constant,diagonal, ne_ni_ratio
    integer :: i,ij,j,k

!*local copy of variables
    nindex0=nlapperp
    indexp0=indexplapperp
    ring0=lapmatperp
    ring00=lapmatperp0

    if(sfcore==1)then
        if(kinetice==0)then
!$omp parallel do private(i,j,ij,rhos2,ne_ni_ratio,coe_temp,dielectric_constant,k)
            do i=1,mpsi_core
                do j=1,mpara_core(i)
                    ij=igrid_core(i)+j
    
                    rhos2 = meshte(ij)*aion/(qion*qion*bmesh(ij)*bmesh(ij))
                    ne_ni_ratio = meshne(ij)/meshni(ij)
                    coe_temp = 1.0_lk + meshti(ij)/meshte(ij)/qion/qion*ne_ni_ratio
                    dielectric_constant = - rhos2*(qion*qion*meshni(ij)/meshte(ij))*coe_temp

                    do k=1,nindex0(ij)
                        ring0(k,ij)=dielectric_constant*ring0(k,ij)
                        ring00(k,ij)=dielectric_constant*ring00(k,ij)
                    enddo
                enddo
            enddo
        else
!$omp parallel do private(i,j,ij,rhos2,ne_ni_ratio,coe_temp,dielectric_constant,k)
            do i=1,mpsi_core
                do j=1,mpara_core(i)
                    ij=igrid_core(i)+j

                    rhos2 = meshte(ij)*aion/(qion*qion*bmesh(ij)*bmesh(ij))
                    dielectric_constant = - rhos2*(qion*qion*meshni(ij)/meshte(ij))

                    do k=1,nindex0(ij)
                        ring0(k,ij)=dielectric_constant*ring0(k,ij)
                        ring00(k,ij)=dielectric_constant*ring00(k,ij)
                    enddo
                enddo
            enddo
        endif
    endif
    
    if(sfsol==1)then
        if(kinetice==0)then
!$omp parallel do private(i,j,ij,rhos2,ne_ni_ratio,coe_temp,dielectric_constant,k)
            do i=1,mpsi_sol
                do j=1,mpara_sol(i)
                    ij=mgrid_core+igrid_sol(i)+j
    
                    rhos2 = meshte(ij)*aion/(qion*qion*bmesh(ij)*bmesh(ij))
                    ne_ni_ratio = meshne(ij)/meshni(ij)
                    coe_temp = 1.0_lk + meshti(ij)/meshte(ij)/qion/qion*ne_ni_ratio
                    dielectric_constant = - rhos2*(qion*qion*meshni(ij)/meshte(ij))*coe_temp
    
                    do k=1,nindex0(ij)
                        ring0(k,ij)=dielectric_constant*ring0(k,ij)
                        ring00(k,ij)=dielectric_constant*ring00(k,ij)
                    enddo
                enddo
            enddo
        else
!$omp parallel do private(i,j,ij,rhos2,ne_ni_ratio,coe_temp,dielectric_constant,k)
            do i=1,mpsi_sol
                do j=1,mpara_sol(i)
                    ij=mgrid_core+igrid_sol(i)+j

                    rhos2 = meshte(ij)*aion/(qion*qion*bmesh(ij)*bmesh(ij))
                    dielectric_constant = - rhos2*(qion*qion*meshni(ij)/meshte(ij))

                    do k=1,nindex0(ij)
                        ring0(k,ij)=dielectric_constant*ring0(k,ij)
                        ring00(k,ij)=dielectric_constant*ring00(k,ij)
                    enddo
                enddo
            enddo
        endif
    endif

if(kinetice==0)then
!========================= add the diagonal term ==========================
! no diagonal zonal term for adiabatic e
    if(sfcore==1)then
!$omp parallel do private(i,j,ij,diagonal)
        do i=2,mpsi_core-1
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                diagonal = meshne(ij)/meshte(ij)
                ring0(5,ij) = ring0(5,ij) + diagonal
            enddo
        enddo
    
        do j=1,mpara_core(1)
            ij=igrid_core(1)+j
            diagonal = meshne(ij)/meshte(ij)
            ring0(1,ij) = ring0(1,ij) + diagonal
        enddo
    
        do j=1,mpara_core(mpsi_core)
            ij=igrid_core(mpsi_core)+j
            diagonal = meshne(ij)/meshte(ij)
            if(couplefrc==0)then
                ring0(1,ij) = ring0(1,ij) + diagonal
            elseif(couplefrc==1)then
                if((j.le.mpara_core_share_index(1)) .and. (j.ge.mpara_core_share_index(mpara_share)))then ! R=0 grids on axis, notice that the grid index increases along the clockwise direction starting at outer midplane in core region
                    ring0(1,ij) = ring0(1,ij) + diagonal
                else
                    ring0(3,ij) = ring0(3,ij) + diagonal
                endif
            else
                if(mype==0)write(gtcout,*)"wrong couplefrc value"
                stop
            endif
        enddo
    endif

! no diagonal zonal term for adiabatic e
    if(sfsol==1)then
!$omp parallel do private(i,j,ij,diagonal)
        do i = 2,mpsi_sol-1
            do j = 2,mpara_sol(i)-1
                ij = mgrid_core+igrid_sol(i)+j
                diagonal = meshne(ij)/meshte(ij)
                ring0(5,ij) = ring0(5,ij) + diagonal
            enddo
            ring0(4,mgrid_core+igrid_sol(i)+1) = ring0(4,mgrid_core+igrid_sol(i)+1) + diagonal
            ring0(6,mgrid_core+igrid_sol(i)+mpara_sol(i)) = ring0(6,mgrid_core+igrid_sol(i)+mpara_sol(i)) + diagonal
        enddo
    
        do j = 1,mpara_sol(1)
            ij = mgrid_core + igrid_sol(1) + j
            diagonal = meshne(ij)/meshte(ij)
            ring0(1,ij) = ring0(1,ij) + diagonal
        enddo
    
        do j = 1,mpara_sol(mpsi_sol)
            ij = mgrid_core + igrid_sol(mpsi_sol) + j
            diagonal = meshne(ij)/meshte(ij)
            ring0(1,ij)=ring0(1,ij) + diagonal
        enddo
    endif
endif

!*sort indexp and ring
    do i=1,mgrid
        do msort=nindex0(i)-1,1,-1   !do i=min,max,step
            do j=1,msort
                if(indexp0(j,i)>indexp0(j+1,i)) then
                    windexp=indexp0(j,i)
                    wring=ring0(j,i)
                    wring0=ring00(j,i)
                    indexp0(j,i)=indexp0(j+1,i)
                    ring0(j,i)=ring0(j+1,i)
                    ring00(j,i)=ring00(j+1,i)
                    indexp0(j+1,i)=windexp
                    ring0(j+1,i)=wring
                    ring00(j+1,i)=wring0
                endif
            enddo
        enddo
    enddo

!*Count nonzero and allocate
    nindex0max=0
    nnz=0
    do i=1,mgrid
        if(nindex0(i).gt.nindex0max)nindex0max=nindex0(i)
        do j=1,nindex0(i)
            nnz=nnz+1
        enddo
    enddo

    allocate(luserafperp(nnz),luserpfperp(nnz),luserbfperp(0:mgrid-1),lusersfperp(-1:mgrid-1),&
             luserxfperp(0:mgrid-1),luserafperp0(nnz),stat=mtest)

    if (mtest /= 0) then
        if(mype==0)write(0,*)mype,'matrix: Cannot allocate userX'
        call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
    endif

!*Make matrix lusera,luserp
    nnz=0
    lusersfperp(-1)=0
    do i=1,mgrid
        lusersfperp(i-1)=lusersfperp(i-2)+nindex0(i)
        do j=1,nindex0(i)
            nnz=nnz+1
            luserafperp(nnz)=ring0(j,i)
            luserafperp0(nnz)=ring00(j,i)
            luserpfperp(nnz)=indexp0(j,i)-1
        enddo
    enddo

end subroutine lap2petsc_perp

#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapsol_initial
    use precision
    use read_in_parameters, only: n_modes
    use global_parameters
    use grids_setup
    use simulation_domain, only: couplefrc
    use grids_eqdata
    integer :: i,j,ij,ijlow,ijhigh,ijsol,ierror,mtest,n,k,kn,xtmp(6),ytmp(6),ntmp(6),tmp,idex(1)

    real(lk) :: h1,h2,wp1,wp2,wp,ws1,ws2,ws,ph1,ph2,sh1,sh2,&
              aa,bb,cc,dd,ee,ff,xx,yy,laptmp(6)
    real(lk) :: jdpgpp(mgrid),jdsgps(mgrid),jdpgps(mgrid),jdsgss(mgrid)

! equilibrium jacobian and metric tensor derivatives are set to be zero. Todo: calculate those equilibrium derivatives using spline function
    do i=1,mpsi_sol
        do j=1,mpara_sol(i)
            ij=mgrid_core+igrid_sol(i)+j
            jdpgpp(ij)=0.0_lk
            jdsgps(ij)=0.0_lk
            jdpgps(ij)=0.0_lk
            jdsgss(ij)=0.0_lk
        enddo
    enddo

    do i=2,mpsi_sol-1
        do j=2,mpara_sol(i)-1
            ijsol=igrid_sol(i)+j
            ij=mgrid_core+igrid_sol(i)+j
            ijlow=mgrid_core+igrid_sol(i-1)+j
            ijhigh=mgrid_core+igrid_sol(i+1)+j
    
            ph1=psimesh_sol(i)-psimesh_sol(i-1)
            ph2=psimesh_sol(i+1)-psimesh_sol(i)
            wp1=ph2/(ph1+ph2)       ! weight of previous grid
            wp2=ph1/(ph1+ph2)       ! weight of later grid
! For coupled SOL and core region case, the separatrix grid is converted to regular Ssol as other fluxsurfaces in SOL    
            sh1=Ssol(ijsol)-Ssol(ijsol-1)
            sh2=Ssol(ijsol+1)-Ssol(ijsol)

            ws1=sh2/(sh1+sh2)       ! weight of previous grid
            ws2=sh1/(sh1+sh2)       ! weight of later grid
    
! Define coefficients in notes for easy implementation
    
            aa=ph1/((ph1+ph2)*ph2)
            bb=-(ph1-ph2)/(ph1*ph2)
            cc=-ph2/((ph1+ph2)*ph1)
            
            dd=sh1/((sh1+sh2)*sh2)
            ee=-(sh1-sh2)/(sh1*sh2)
            ff=-sh2/((sh1+sh2)*sh1)
    
            xx=jmesh(ij)*(jdpgpp(ij)+jdsgps(ij))
            yy=jmesh(ij)*(jdpgps(ij)+jdsgss(ij))
    
    
            nlapperp(ij)=mlapperp
         
            indexplapperp(1,ij)=ijlow-1
            lapmatperp(1,ij)=2.0_lk*gups(ij)*cc*ff
     
            indexplapperp(2,ij)=ijlow
            lapmatperp(2,ij)=gupp(ij)*(2.0_lk/ph1/ph2)*wp1+2.0_lk*gups(ij)*cc*ee+cc*xx
            
            indexplapperp(3,ij)=ijlow+1
            lapmatperp(3,ij)=2.0_lk*gups(ij)*cc*dd
  
            indexplapperp(4,ij)=ij-1
            lapmatperp(4,ij)=guss(ij)*(2.0_lk/sh1/sh2)*ws1+2.0_lk*gups(ij)*bb*ff+ff*yy
  
            indexplapperp(5,ij)=ij
            lapmatperp(5,ij)=-gupp(ij)*(2.0_lk/ph1/ph2)-guss(ij)*(2.0_lk/sh1/sh2)+2.0_lk*gups(ij)*bb*ee+bb*xx+ee*yy-real(nmode)*real(nmode)*real(n_modes(1))*real(n_modes(1))*guzz(ij)
  
            indexplapperp(6,ij)=ij+1
            lapmatperp(6,ij)=guss(ij)*(2.0_lk/sh1/sh2)*ws2+2.0_lk*gups(ij)*bb*dd+dd*yy
  
  
            indexplapperp(7,ij)=ijhigh-1
            lapmatperp(7,ij)=2.0_lk*gups(ij)*aa*ff
          
            indexplapperp(8,ij)=ijhigh
            lapmatperp(8,ij)=gupp(ij)*(2.0_lk/ph1/ph2)*wp2+2.0_lk*gups(ij)*aa*ee+aa*xx
                    
            indexplapperp(9,ij)=ijhigh+1
            lapmatperp(9,ij)=2.0_lk*gups(ij)*aa*dd
            lapmatperp0(:,ij) = lapmatperp(:,ij)
            lapmatperp0(5,ij) = -gupp(ij)*(2.0_lk/ph1/ph2)-guss(ij)*(2.0_lk/sh1/sh2)+2.0_lk*gups(ij)*bb*ee+bb*xx+ee*yy
        enddo
    enddo

    do i=2,mpsi_sol-1
        ijsol=igrid_sol(i)+1
        ij=mgrid_core+igrid_sol(i)+1
        ijlow=mgrid_core+igrid_sol(i-1)+1
        ijhigh=mgrid_core+igrid_sol(i+1)+1

        ph1=psimesh_sol(i)-psimesh_sol(i-1)
        ph2=psimesh_sol(i+1)-psimesh_sol(i)
        wp1=ph2/(ph1+ph2)       ! weight of previous grid
        wp2=ph1/(ph1+ph2)       ! weight of later grid

        sh1=Ssol(ijsol+1)-Ssol(ijsol)
        sh2=Ssol(ijsol+2)-Ssol(ijsol+1)
        ws1=sh2/(sh1+sh2)       ! weight of previous grid
        ws2=sh1/(sh1+sh2)       ! weight of later grid

! Define coefficients in notes for easy implementation

        aa=ph1/((ph1+ph2)*ph2)
        bb=-(ph1-ph2)/(ph1*ph2)
        cc=-ph2/((ph1+ph2)*ph1)

        dd=-sh1/((sh1+sh2)*sh2)
        ee=(sh1+sh2)/(sh1*sh2)
        ff=-(2.0_lk*sh1+sh2)/((sh1+sh2)*sh1)

        xx=jmesh(ij)*(jdpgpp(ij)+jdsgps(ij))
        yy=jmesh(ij)*(jdpgps(ij)+jdsgss(ij))


        nlapperp(ij)=mlapperp

        indexplapperp(1,ij)=ijlow
        lapmatperp(1,ij)=gupp(ij)*(2.0_lk/ph1/ph2)*wp1+2.0_lk*gups(ij)*cc*ff+cc*xx

        indexplapperp(2,ij)=ijlow+1
        lapmatperp(2,ij)=2.0_lk*gups(ij)*cc*ee

        indexplapperp(3,ij)=ijlow+2
        lapmatperp(3,ij)=2.0_lk*gups(ij)*cc*dd

        indexplapperp(4,ij)=ij
        lapmatperp(4,ij)=-gupp(ij)*(2.0_lk/ph1/ph2)+guss(ij)*(2.0_lk/sh1/sh2)*ws1+2.0_lk*gups(ij)*bb*ff+bb*xx+ff*yy-real(nmode)*real(nmode)*real(n_modes(1))*real(n_modes(1))*guzz(ij)

        indexplapperp(5,ij)=ij+1
        lapmatperp(5,ij)=-guss(ij)*(2.0_lk/sh1/sh2)+2.0_lk*gups(ij)*bb*ee+ee*yy

        indexplapperp(6,ij)=ij+2
        lapmatperp(6,ij)=guss(ij)*(2.0_lk/sh1/sh2)*ws2+2.0_lk*gups(ij)*bb*dd+dd*yy

        indexplapperp(7,ij)=ijhigh
        lapmatperp(7,ij)=gupp(ij)*(2.0_lk/ph1/ph2)*wp2+2.0_lk*gups(ij)*aa*ff+aa*xx

        indexplapperp(8,ij)=ijhigh+1
        lapmatperp(8,ij)=2.0_lk*gups(ij)*aa*ee

        indexplapperp(9,ij)=ijhigh+2
        lapmatperp(9,ij)=2.0_lk*gups(ij)*aa*dd
        lapmatperp0(:,ij) = lapmatperp(:,ij)
        lapmatperp0(4,ij) = -gupp(ij)*(2.0_lk/ph1/ph2)+guss(ij)*(2.0_lk/sh1/sh2)*ws1+2.0_lk*gups(ij)*bb*ff+bb*xx+ff*yy
    enddo

    do i=2,mpsi_sol-1
        ijsol=igrid_sol(i)+mpara_sol(i)
        ij=mgrid_core+ijsol
        ijlow=ij-mpara_sol(i) ! require mpara_sol are the same for different flux surfaces
        ijhigh=ij+mpara_sol(i)

        ph1=psimesh_sol(i)-psimesh_sol(i-1)
        ph2=psimesh_sol(i+1)-psimesh_sol(i)
        wp1=ph2/(ph1+ph2)       ! weight of previous grid
        wp2=ph1/(ph1+ph2)       ! weight of later grid

        sh1=Ssol(ijsol)-Ssol(ijsol-1)
        sh2=Ssol(ijsol-1)-Ssol(ijsol-2)
        ws1=sh2/(sh1+sh2)       ! weight of previous grid
        ws2=sh1/(sh1+sh2)       ! weight of later grid

! Define coefficients in notes for easy implementation

        aa=ph1/((ph1+ph2)*ph2)
        bb=-(ph1-ph2)/(ph1*ph2)
        cc=-ph2/((ph1+ph2)*ph1)

        dd=sh1/((sh1+sh2)*sh2)
        ee=-(sh1+sh2)/(sh1*sh2)
        ff=(2.0_lk*sh1+sh2)/((sh1+sh2)*sh1)

        xx=jmesh(ij)*(jdpgpp(ij)+jdsgps(ij))
        yy=jmesh(ij)*(jdpgps(ij)+jdsgss(ij))

        nlapperp(ij)=mlapperp

        indexplapperp(1,ij)=ijlow-2
        lapmatperp(1,ij)=2.0_lk*gups(ij)*cc*dd

        indexplapperp(2,ij)=ijlow-1
        lapmatperp(2,ij)=2.0_lk*gups(ij)*cc*ee

        indexplapperp(3,ij)=ijlow
        lapmatperp(3,ij)=gupp(ij)*(2.0_lk/ph1/ph2)*wp1+2.0_lk*gups(ij)*cc*ff+cc*xx

        indexplapperp(4,ij)=ij-2
        lapmatperp(4,ij)=guss(ij)*(2.0_lk/sh1/sh2)*ws2+2.0_lk*gups(ij)*bb*dd+dd*yy

        indexplapperp(5,ij)=ij-1
        lapmatperp(5,ij)=-guss(ij)*(2.0_lk/sh1/sh2)+2.0_lk*gups(ij)*bb*ee+ee*yy

        indexplapperp(6,ij)=ij
        lapmatperp(6,ij)=-gupp(ij)*(2.0_lk/ph1/ph2)+guss(ij)*(2.0_lk/sh1/sh2)*ws1+2.0_lk*gups(ij)*bb*ff+bb*xx+ff*yy-real(nmode)*real(nmode)*real(n_modes(1))*real(n_modes(1))*guzz(ij)

        indexplapperp(7,ij)=ijhigh-2
        lapmatperp(7,ij)=2.0_lk*gups(ij)*aa*dd

        indexplapperp(8,ij)=ijhigh-1
        lapmatperp(8,ij)=2.0_lk*gups(ij)*aa*ee

        indexplapperp(9,ij)=ijhigh
        lapmatperp(9,ij)=gupp(ij)*(2.0_lk/ph1/ph2)*wp2+2.0_lk*gups(ij)*aa*ff+aa*xx
        lapmatperp0(:,ij) = lapmatperp(:,ij)
        lapmatperp0(6,ij) = -gupp(ij)*(2.0_lk/ph1/ph2)+guss(ij)*(2.0_lk/sh1/sh2)*ws1+2.0_lk*gups(ij)*bb*ff+bb*xx+ff*yy
    enddo

! Define the outer simulation boundary
    do j=1,mpara_sol(mpsi_sol)
        ij=mgrid_core+igrid_sol(mpsi_sol)+j
        nlapperp(ij)=1
        indexplapperp(1,ij)=ij
        lapmatperp(1,ij)=1.0_lk
        lapmatperp0(1,ij)=1.0_lk
    enddo

    do j=1,mpara_sol(1)
        ij=mgrid_core+igrid_sol(1)+j
        nlapperp(ij)=1
        indexplapperp(1,ij)=ij
        lapmatperp(1,ij)=1.0_lk
        lapmatperp0(1,ij)=1.0_lk
    enddo

    if(mype==0)write(gtcout,*)"lapsol_initial is used"

end subroutine lapsol_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapcore_initial
    use precision
    use read_in_parameters, only: n_modes
    use global_parameters
    use grids_setup
    use grids_eqdata
    integer :: i,j,ij,ijlow,ijhigh,ijlows,ijhighs,&
             ij0,ij2,ijlow0,ijlow2,ijhigh0,ijhigh2,ierror,mtest
    real(lk) :: h1,h2,wp1,wp2,wp,ws1,ws2,ws,ph1,ph2,sh1,sh2,&
              aa,bb,cc,dd,ee,ff,xx,yy
    real(lk) :: jdpgpp(mgrid),jdsgps(mgrid),jdpgps(mgrid),jdsgss(mgrid)

! equilibrium jacobian and metric tensor derivatives are set to be zero. Todo: calculate those equilibrium derivatives using spline function
    do i=1,mpsi_core
        do j=1,mpara_core(i)
            ij=igrid_core(i)+j          
            jdpgpp(ij)=0.0_lk         
            jdsgps(ij)=0.0_lk         
            jdpgps(ij)=0.0_lk          
            jdsgss(ij)=0.0_lk
        enddo
    enddo
! Require grid number to be the same in each flux surfaces!
    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

            if(j==1 .or. j==mpara_core(i))then
                ij2=igrid_core(i)+2
                ij0=igrid_core(i)+mpara_core(i)-1

                ijlow2=igrid_core(i-1)+2
                ijlow0=igrid_core(i-1)+mpara_core(i-1)-1

                ijhigh2=igrid_core(i+1)+2
                ijhigh0=igrid_core(i+1)+mpara_core(i+1)-1
            else
                ij2=ij+1
                ij0=ij-1
          
                ijlow2=ijlow+1
                ijlow0=ijlow-1

                ijhigh2=ijhigh+1
                ijhigh0=ijhigh-1
            endif


            ph1=psimesh_core(i)-psimesh_core(i-1)
            ph2=psimesh_core(i+1)-psimesh_core(i)
            wp1=ph2/(ph1+ph2)       ! weight of previous grid
            wp2=ph1/(ph1+ph2)       ! weight of later grid

            sh1=Score(ij)-Score(ij0)
            sh2=Score(ij2)-Score(ij)
            ws1=sh2/(sh1+sh2)       ! weight of previous grid
            ws2=sh1/(sh1+sh2)       ! weight of later grid

! Define coefficients in notes for easy implementation

            aa=ph1/((ph1+ph2)*ph2)
            bb=-(ph1-ph2)/(ph1*ph2)
            cc=-ph2/((ph1+ph2)*ph1)
   
            dd=sh1/((sh1+sh2)*sh2)
            ee=-(sh1-sh2)/(sh1*sh2)
            ff=-sh2/((sh1+sh2)*sh1)
   
            xx=jmesh(ij)*(jdpgpp(ij)+jdsgps(ij))
            yy=jmesh(ij)*(jdpgps(ij)+jdsgss(ij))
   
            nlapperp(ij)=mlapperp
   
            indexplapperp(1,ij)=ijlow0
            lapmatperp(1,ij)=2.0_lk*gups(ij)*cc*ff
   
            indexplapperp(2,ij)=ijlow
            lapmatperp(2,ij)=gupp(ij)*(2.0_lk/ph1/ph2)*wp1+2.0_lk*gups(ij)*cc*ee+cc*xx
   
            indexplapperp(3,ij)=ijlow2
            lapmatperp(3,ij)=2.0_lk*gups(ij)*cc*dd
   
            indexplapperp(4,ij)=ij0
            lapmatperp(4,ij)=guss(ij)*(2.0_lk/sh1/sh2)*ws1+2.0_lk*gups(ij)*bb*ff+ff*yy
   
            indexplapperp(5,ij)=ij
            lapmatperp(5,ij)=-gupp(ij)*(2.0_lk/ph1/ph2)-guss(ij)*(2.0_lk/sh1/sh2)+2.0_lk*gups(ij)*bb*ee+bb*xx+ee*yy-real(nmode)*real(nmode)*real(n_modes(1))*real(n_modes(1))*guzz(ij)
   
            indexplapperp(6,ij)=ij2
            lapmatperp(6,ij)=guss(ij)*(2.0_lk/sh1/sh2)*ws2+2.0_lk*gups(ij)*bb*dd+dd*yy
   
            indexplapperp(7,ij)=ijhigh0
            lapmatperp(7,ij)=2.0_lk*gups(ij)*aa*ff
   
            indexplapperp(8,ij)=ijhigh
            lapmatperp(8,ij)=gupp(ij)*(2.0_lk/ph1/ph2)*wp2+2.0_lk*gups(ij)*aa*ee+aa*xx
   
            indexplapperp(9,ij)=ijhigh2
            lapmatperp(9,ij)=2.0_lk*gups(ij)*aa*dd
            lapmatperp0(:,ij)=lapmatperp(:,ij)
            lapmatperp0(5,ij)=-gupp(ij)*(2.0_lk/ph1/ph2)-guss(ij)*(2.0_lk/sh1/sh2)+2.0_lk*gups(ij)*bb*ee+bb*xx+ee*yy
        enddo
    enddo

!Define the inner boundary
    do j=1,mpara_core(1)
        ij=igrid_core(1)+j
        nlapperp(ij)=1
        indexplapperp(1,ij)=ij
        lapmatperp(1,ij)=1.0_lk
        lapmatperp0(1,ij)=1.0_lk
    enddo

!Define the outer boundary for couplefrc==0
    do j=1,mpara_core(mpsi_core)
        ij=igrid_core(mpsi_core)+j
        nlapperp(ij)=1
        indexplapperp(1,ij)=ij
        lapmatperp(1,ij)=1.0_lk
        lapmatperp0(1,ij)=1.0_lk
    enddo
end subroutine lapcore_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapcoup_initial
    use precision
    use read_in_parameters, only: n_modes
    use global_parameters
    use grids_setup
    use grids_eqdata

    integer :: i,j,ij,ij1,ijn
    real(lk) :: h1,h2,w1,w2,wsol,wcore

! Define the shared grid on separatrix which are not on axis (R=0)
    do i=2,mpara_share-1
        ij=igrid_core(mpsi_core)+mpara_core_share_index(i) 
        h1=psimesh_core(mpsi_core)-psimesh_core(mpsi_core-1)
        h2=psimesh_sol(2)-psimesh_sol(1)
        w1=h2/(h1+h2)
        w2=h1/(h1+h2)
        wsol=interp_weight_sol_share(i)
        wcore=interp_weight_core_share(i)
    
        nlapperp(ij)=5
  
        indexplapperp(1,ij)=igrid_core(mpsi_core-1)+interp_index_core_share(i)
        lapmatperp(1,ij)=(2.0_lk/(h1*h2))*w1*(1.0_lk-wcore)*gupp(ij)
  
        indexplapperp(2,ij)=igrid_core(mpsi_core-1)+interp_index_core_share(i)+1
        lapmatperp(2,ij)=(2.0_lk/(h1*h2))*w1*wcore*gupp(ij)
     
        indexplapperp(3,ij)=ij
        lapmatperp(3,ij)=-(2.0_lk/(h1*h2))*gupp(ij)-real(nmode)*real(nmode)*real(n_modes(1))*real(n_modes(1))*guzz(ij)
      
        indexplapperp(4,ij)=mgrid_core+igrid_sol(2)+interp_index_sol_share(i)
        lapmatperp(4,ij)=(2.0_lk/(h1*h2))*w2*(1.0_lk-wsol)*gupp(ij)
  
        indexplapperp(5,ij)=mgrid_core+igrid_sol(2)+interp_index_sol_share(i)+1
        lapmatperp(5,ij)=(2.0_lk/(h1*h2))*w2*wsol*gupp(ij)
        lapmatperp0(:,ij)=lapmatperp(:,ij)
        lapmatperp0(3,ij)=-(2.0_lk/(h1*h2))*gupp(ij)
    enddo

! Periodicity in core
    ij1=igrid_core(mpsi_core)+1
    ijn=igrid_core(mpsi_core)+mpara_core(mpsi_core)

    nlapperp(ijn)=nlapperp(ij1)
    indexplapperp(1:2,ijn)=indexplapperp(1:2,ij1)
    indexplapperp(4:5,ijn)=indexplapperp(4:5,ij1)
    indexplapperp(3,ijn)=ijn
    lapmatperp(:,ijn)=lapmatperp(:,ij1)
    lapmatperp0(:,ijn)=lapmatperp0(:,ij1)

    if(mype==0)write(gtcout,*)"lapcoup_initial is used"
end subroutine lapcoup_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine laplacian_perp(farray,scalarout)
    use global_parameters
    use simulation_domain, only: couplefrc
    use grids_setup
    use grids_fields, only: periodicity

    real(lk),dimension(mgrid),intent(in) :: farray
    real(lk),dimension(mgrid),intent(inout) :: scalarout

    integer :: i,j
    real(lk),dimension(mgrid) :: scalarin
    integer :: ierror
  
    scalarout=0.0_lk
    scalarin(:)=farray(:)
  
    if(couplefrc==1)then
! convert the fields on irregular mesh to the one on regular mesh on separatrix, because the laplacian operator is written based on regular mesh
        call convert_irreg2reg(scalarin)
    endif

!$omp parallel do private(i)
    do i=1,mgrid
        do j=1,nlapperp(i)
            scalarout(i)=scalarout(i)+lapmatperp(j,i)*scalarin(indexplapperp(j,i))
        enddo
    enddo
  
    if(couplefrc==1)then
! convert the fields on regular mesh back to the one on irregular mesh on separatrix
         call convert_reg2irreg(scalarout)
! evaluate the fields on shared grids in SOL region using the one on shared grids in core region 
         do i = 1, mpara_share
             scalarout(mgrid_core+igrid_sol(1)+mpara_sol_share_index(i)) = scalarout(igrid_core(mpsi_core)+mpara_core_share_index(i))
         enddo
    endif

end subroutine laplacian_perp
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lapdebug_initial
    use precision
    use global_parameters
    use grids_setup

    integer :: i,j,ij

    do i=1,mpsi_core
        do j=1,mpara_core(i)
            ij=igrid_core(i)+j

            nlapperp(ij)=1
            indexplapperp(1,ij)=ij
            lapmatperp(1,ij)=1.0_lk
        enddo
    enddo

    do i=1,mpsi_sol
        do j=1,mpara_sol(i)
            ij=mgrid_core+igrid_sol(i)+j

            nlapperp(ij)=1
            indexplapperp(1,ij)=ij
            lapmatperp(1,ij)=1.0_lk
        enddo
    enddo

end subroutine lapdebug_initial
end module laplacian
