subroutine lappara_initial
     use precision
     use global_parameters
     use laplacian
     use grids_setup
     use simulation_domain, only: sfcore, sfsol
     use petsc, only: lapmat_pscreate_fpara
     use petsc_array, only: luserafpara, luserpfpara, lusersfpara
     implicit none

     integer :: i,j,ij,ijlow,ijhigh,ij1,ijn,ijsol,ierror,mtest
     real(lk) :: h1,h2,weight1,weight2
     real(lk) :: sh1,sh2,ws1,ws2,sh,ws


     mlappara=3
!!allocate array
     allocate(nlappara(mgrid_frc),indexplappara(mlappara,mgrid_frc),&
              lapmatpara(mlappara,mgrid_frc),stat=mtest)
     if (mtest /= 0) then
       if(mype==0)write(0,*)mype,'lapmat_initial: Cannot allocate lapmat'
       call MPI_ABORT(MPI_COMM_WORLD,1,ierror)
     endif

     lapmatpara=0.0_lk
     nlappara=0
     indexplappara=0


     if(sfcore==1)then
!$omp parallel do private(i,j,ij,ijlow,ijhigh,sh,ws)
        do i=1,mpsicore
           do j=1,mtheta_core(i)
         
              ij=igrid_core(i)+j

              if(j==1 .or. j==mtheta_core(i))then
                ijlow=igrid_core(i)+mtheta_core(i)-1
                ijhigh=igrid_core(i)+2
              else
                ijlow=ij-1
                ijhigh=ij+1
              endif

              sh=deltas_core(i)
              ws=0.5_lk

              nlappara(ij)=mlappara

              indexplappara(1,ij)=ijlow
              lapmatpara(1,ij)=(2.0_lk/(sh*sh))*ws

              indexplappara(2,ij)=ij
              lapmatpara(2,ij)=-2.0_lk/(sh*sh)

              indexplappara(3,ij)=ijhigh
              lapmatpara(3,ij)=(2.0_lk/(sh*sh))*ws

           enddo
        enddo
if(.false.)then
! The first grid
!$omp parallel do private(i,ij1)
        do i=1,mpsicore
           ij1=igrid_core(i)+1

           nlappara(ij1)=mlappara-2
           indexplappara(1,ij1)=ij1
           lapmatpara(1,ij1)=1.0_lk
        enddo

! The last grid
!$omp parallel do private(i,ijn)
        do i=1,mpsicore
           ijn=igrid_core(i)+mtheta_core(i)

           nlappara(ijn)=mlappara-2
           indexplappara(1,ijn)=ijn
           lapmatpara(1,ijn)=1.0_lk
        enddo
endif


     endif



     if(sfsol==1)then
!$omp parallel do private(i,j,ij,ijsol,sh1,sh2,ws1,ws2)
        do i=1,mpsisol
           do j=2,mtheta_sol(i)-1
              ij=mgrid_core+igrid_sol(i)+j
              ijsol=igrid_sol(i)+j

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

              nlappara(ij)=mlappara

              indexplappara(1,ij)=ij-1
              lapmatpara(1,ij)=(2.0_lk/(sh1*sh2))*ws1

              indexplappara(2,ij)=ij
              lapmatpara(2,ij)=-2.0_lk/(sh1*sh2)

              indexplappara(3,ij)=ij+1
              lapmatpara(3,ij)=(2.0_lk/(sh1*sh2))*ws2

              if(mype==0 .and. ijsol==120)write(gtcout,*)"deltas_sol=",deltas_sol(120)
              if(mype==0 .and. ijsol==120)write(gtcout,*)"Ssol(121)-Ssol(120)=",Ssol(121)-Ssol(120)

           enddo
        enddo

! The first grid
!$omp parallel do private(i,ij1)
        do i=1,mpsisol
           ij1=mgrid_core+igrid_sol(i)+1

           nlappara(ij1)=mlappara-2
           indexplappara(1,ij1)=ij1
           lapmatpara(1,ij1)=1.0_lk
        enddo



if(.true.)then
!$omp parallel do private(i,ij,ijsol,sh1,sh2,ws1,ws2)
        do i=1,mpsisol
              ij=mgrid_core+igrid_sol(i)+2
              ijsol=igrid_sol(i)+2

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

              nlappara(ij)=mlappara

              indexplappara(1,ij)=ij-1
              lapmatpara(1,ij)=-(2.0_lk/(sh1*sh2))*ws1*(sh1*sh2*(sh1+sh2)/(2.0_lk*sh1*sh2+sh2*sh2))

              indexplappara(2,ij)=ij
              lapmatpara(2,ij)=2.0_lk/(sh1*sh2)*(ws1*(sh1+sh2)*(sh1+sh2)/(2.0_lk*sh1*sh2+sh2*sh2)-1.0_lk)

              indexplappara(3,ij)=ij+1
              lapmatpara(3,ij)=(2.0_lk/(sh1*sh2))*(ws2-ws1*(sh1*sh1)/(2.0_lk*sh1*sh2+sh2*sh2))
        enddo
endif



if(.false.)then
!$omp parallel do private(i,ij,ijsol,sh1,sh2,ws1,ws2)
        do i=1,mpsisol
              ij=mgrid_core+igrid_sol(i)+mtheta_sol(i)-1
              ijsol=igrid_sol(i)+mtheta_sol(i)-1

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

              nlappara(ij)=mlappara

              indexplappara(1,ij)=ij-1
              lapmatpara(1,ij)=(2.0_lk/(sh1*sh2))*(ws1-ws2*(sh2*sh2)/(2.0_lk*sh1*sh2+sh1*sh1))

              indexplappara(2,ij)=ij
              lapmatpara(2,ij)=(2.0_lk/(sh1*sh2))*(ws2*(sh1+sh2)*(sh1+sh2)/(2.0_lk*sh1*sh2+sh1*sh1)-1.0_lk)

              indexplappara(3,ij)=ij+1
              lapmatpara(3,ij)=(2.0_lk/(sh1*sh2))*ws2*(sh1*sh2*(sh1+sh2)/(2.0_lk*sh1*sh2+sh1*sh1))
        enddo
endif

! The last grid
!$omp parallel do private(i,ijn)
        do i=1,mpsisol
           ijn=mgrid_core+igrid_sol(i)+mtheta_sol(i)
!
           nlappara(ijn)=mlappara-2
           indexplappara(1,ijn)=ijn
           lapmatpara(1,ijn)=1.0_lk
        enddo

     endif






if(.false.)then
! The first grid
!$omp parallel do private(i,ij1,ijsol,h1,h2,weight1,weight2)
    do i=1,mpsisol
       ij1=mgrid_core+igrid_sol(i)+1
       ijsol=igrid_sol(i)+1

       h1=1.0_lk/deltas_sol(ijsol)
       h2=1.0_lk/deltas_sol(ijsol+1)
! Notice h1 and h2 is the reciprocal of deltas_sol
       weight1=h1/(h1+h2)       ! weight of current(first) grid
       weight2=h2/(h1+h2)       ! weight of last(third) grid

       nlappara(ij1)=mlappara

       indexplappara(1,ij1)=ij1
       lapmatpara(1,ij1)=2.0_lk*h1*h2*weight1

       indexplappara(2,ij1)=ij1+1
       lapmatpara(2,ij1)=-2.0_lk*h1*h2

       indexplappara(3,ij1)=ij1+2
       lapmatpara(3,ij1)=2.0_lk*h1*h2*weight2

    enddo

! The last grid
!$omp parallel do private(i,ijn,ijsol,h1,h2,weight1,weight2)
    do i=1,mpsisol
       ijn=mgrid_core+igrid_sol(i)+mtheta_sol(i)
       ijsol=igrid_sol(i)+mtheta_sol(i)

       h1=1.0_lk/deltas_sol(ijsol-1)
       h2=1.0_lk/deltas_sol(ijsol-2)

! Notice h1 and h2 is the reciprocal of deltas_sol
       weight1=h1/(h1+h2)       ! weight of current(first) grid
       weight2=h2/(h1+h2)       ! weight of last(third) grid


       nlappara(ijn)=mlappara

       indexplappara(1,ijn)=ijn-2
       lapmatpara(1,ijn)=2.0_lk*h1*h2*weight2

       indexplappara(2,ijn)=ijn-1
       lapmatpara(2,ijn)=-2.0_lk*h1*h2

       indexplappara(3,ijn)=ijn
       lapmatpara(3,ijn)=2.0_lk*h1*h2*weight1

    enddo
endif

#ifdef _PETSc
     call lap2petsc_para
     call lapmat_pscreate_fpara(mlappara,mgrid_frc, luserafpara, luserpfpara, lusersfpara)
#endif

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

end subroutine lappara_initial
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine lap2petsc_para
!*Using nindex0,indexp0,ring0 generate matrices users, usera, and userp
     use precision
     use global_parameters
     use laplacian
     use grids_setup
     use petsc_array,only: luserafpara,lusersfpara,luserpfpara,luserbfpara,luserxfpara
     implicit none

     integer nnz,windexp,msort,nindex0max
     integer,dimension(mgrid_frc) :: nindex0
     integer,dimension(mlappara,mgrid_frc) :: indexp0
     real(lk),dimension(mlappara,mgrid_frc) :: ring0

     integer :: ierror,mtest
     real(lk) :: wring
     integer :: i,ij,j,k


!*local copy of variables
     nindex0=nlappara
     indexp0=indexplappara
     ring0=lapmatpara
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     ring0=-ring0

!*sort indexp and ring
     do i=1,mgrid_frc
       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)
             indexp0(j,i)=indexp0(j+1,i)
             ring0(j,i)=ring0(j+1,i)
             indexp0(j+1,i)=windexp
             ring0(j+1,i)=wring
           endif
         enddo
       enddo
     enddo

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

     allocate(luserafpara(nnz),luserpfpara(nnz),luserbfpara(0:mgrid_frc-1),lusersfpara(-1:mgrid_frc-1),&
              luserxfpara(0:mgrid_frc-1),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
     if(mype==0)write(gtcout,*)"nnz=",nnz

!*Make matrix lusera,luserp
     nnz=0
     lusersfpara(-1)=0
     do i=1,mgrid_frc
       lusersfpara(i-1)=lusersfpara(i-2)+nindex0(i)
       do j=1,nindex0(i)
         nnz=nnz+1
         luserafpara(nnz)=ring0(j,i)
         luserpfpara(nnz)=indexp0(j,i)-1
       enddo
     enddo

end subroutine lap2petsc_para
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine laplacian_para(scalarin,scalarout)
    use precision
    use global_parameters
    use grids_fields
    use grids_setup
    use laplacian
    implicit none

    integer :: i,j
    real(lk),dimension(mgrid_frc),intent(in) :: scalarin
    real(lk),dimension(0:1,mgrid_frc),intent(out) :: scalarout


    scalarout=0.0_lk

!$omp parallel do private(i)
    do i=1,mgrid_frc
       do j=1,nlappara(i)
         scalarout(1,i)=scalarout(1,i)+lapmatpara(j,i)*scalarin(indexplappara(j,i))
        enddo
    enddo

    call cy_periodicity(scalarout)

end subroutine laplacian_para
