module grids_setup
use precision
implicit none

public
private core_initial_mype0, sol_initial_mype0, sep_initial, &
        irreg2reg_interp_sol, reg2irreg_interp_sol, orthogonal_grid_sep, grids_dataout_mype0
!==================================================================
integer mpsi_sol, mpsi_core, mgrid_core, mgrid_sol, mpsi, mgrid, &
        mpara_sep_sol, & ! grid number on the separatrix in sol region
        mpara_share, & ! the shared grid on separatrix for core and sol regions
        iflux ! diag & uniform profile reference fluxsuface

integer, dimension(:), allocatable :: igrid_sol, mpara_sol, igrid_core, mpara_core, &   
                                      mpara_core_share_index, mpara_sol_share_index   ! mpara_core_share_index(mpara_share) for grid index in core region, mpara_sol_share_index(mpara_share) for grid index in sol region   
real(lk) deltar_sol, deltar_core   ! interval of different flux surfaces of simulation grids
real(lk), dimension(:), allocatable :: rg1d_core, zg1d_core, psimesh_core, & ! mpsi_core length
                                       rg1d_sol, zg1d_sol, psimesh_sol, &   ! mpsi_sol length
                                       Rcore, Zcore, Score, Thetacore, &  ! igrid_core length
                                       Rsol, Zsol, Ssol,&  ! igrid_sol length
                                       totals_core, deltas_core1d, & ! length of each fieldline of simulation flux surface/ parallel interval between each simulation grid
                                       totals_sol, deltas_sol

! mpara_share length
integer,dimension(:),allocatable :: interp_index_core_share,interp_index_sol_share
real(lk),dimension(:),allocatable :: interp_weight_core_share,interp_weight_sol_share
real(lk),dimension(:),allocatable :: Rsol_share,Zsol_share,Ssol_share,&
                                     Rcore_share,Zcore_share,Score_share,&
                                     dcore_orthogonal,dsol_orthogonal

! mpara_sep_sol length
integer,dimension(:),allocatable :: irreg2reg_index_sol, reg2irreg_index_sol
real(lk),dimension(:),allocatable :: irreg2reg_weight_sol, reg2irreg_weight_sol, &
                                     Rsol_sep_reg, Zsol_sep_reg, Ssol_sep_reg, &
                                     Rsol_sep_irreg, Zsol_sep_irreg, Ssol_sep_irreg
save
#ifndef GPU_UM
!$acc declare create(psimesh_core, psimesh_sol, Ssol, igrid_core, igrid_sol, mpara_sol, mpara_core)
#endif
contains
!==============================================================
subroutine frc_mesh_create
    use global_parameters,only:gtcout,mype
    use spline_function, only: spline2d
    use equilibrium, only: psi_sep
    use simulation_domain, only: sfcore, sfsol, couplefrc
    use read_in_parameters, only: mpsisol_input, mpsicore_input, mparasol_input, mparacore_input, diag_region, iflux_sol, iflux_core

    integer :: i,mtest,icount,ierror
 
    if(sfcore==1)then
        
        mpsi_core=mpsicore_input

        allocate(mpara_core(mpsi_core))
        do i=1,mpsi_core
            mpara_core(i)=mparacore_input
        enddo
        mgrid_core=0
        do i=1,mpsi_core
           mgrid_core = mgrid_core + mpara_core(i)
        enddo

        allocate(rg1d_core(mpsi_core),zg1d_core(mpsi_core),psimesh_core(mpsi_core))
        allocate(totals_core(mpsi_core),deltas_core1d(mpsi_core))
        allocate(Rcore(mgrid_core),Zcore(mgrid_core))
        allocate(igrid_core(mpsi_core))
        allocate(Score(mgrid_core),Thetacore(mgrid_core))
    
        if(mype==0)then
            call core_initial_mype0
        endif
        call MPI_Barrier(MPI_COMM_WORLD,ierror)
    
! Broadcast the grid information of core region to each MPIs
        icount=mpsi_core
        call MPI_BCAST(rg1d_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(zg1d_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(psimesh_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(totals_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(deltas_core1d,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(igrid_core,icount,mpi_Integer,0,MPI_COMM_WORLD,ierror)
    
        icount=mgrid_core
        call MPI_BCAST(Rcore,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(Zcore,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(Score,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(Thetacore,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
#ifndef GPU_UM
        !$acc update device(psimesh_core,igrid_core,mpara_core)
#endif
    endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    if(sfsol==1)then
        mpsi_sol=mpsisol_input
    
        allocate(rg1d_sol(mpsi_sol),zg1d_sol(mpsi_sol),psimesh_sol(mpsi_sol),totals_sol(mpsi_sol))
        allocate(igrid_sol(mpsi_sol),mpara_sol(mpsi_sol))
        rg1d_sol=0.0_lk
        zg1d_sol=0.0_lk
        psimesh_sol=0.0_lk
        totals_sol=0.0_lk
        igrid_sol=0
        mpara_sol=0
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! give the grid number on separatrix in SOL region
        mpara_sol(1)=mparasol_input
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
        if(couplefrc==1)then
            call sep_initial
        endif
    
        mpara_sol(:)=mpara_sol(1)
    
        igrid_sol(1)=0
        do i=2,mpsi_sol
            igrid_sol(i)=igrid_sol(i-1)+mpara_sol(i-1)
        enddo
        mgrid_sol=0
        do i=1,mpsi_sol
           mgrid_sol=mgrid_sol+mpara_sol(i)
        enddo
    
        allocate(Rsol(mgrid_sol),Zsol(mgrid_sol),Ssol(mgrid_sol),deltas_sol(mgrid_sol))
        Rsol=0.0_lk
        Zsol=0.0_lk
        Ssol=0.0_lk
        deltas_sol=0.0_lk

        if(couplefrc==1)then
            allocate(Rsol_sep_reg(mpara_sep_sol),Zsol_sep_reg(mpara_sep_sol),Ssol_sep_reg(mpara_sep_sol))
            Rsol_sep_reg=0.0_lk
            Zsol_sep_reg=0.0_lk
            Ssol_sep_reg=0.0_lk
        endif
       
        if(mype==0)then
            call sol_initial_mype0
        endif
        call MPI_Barrier(MPI_COMM_WORLD,ierror)
! Broadcast the grid information of sol region to each MPIs
        icount=mgrid_sol
        call MPI_BCAST(Rsol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(Zsol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(Ssol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(deltas_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)

        icount=mpsi_sol
        call MPI_BCAST(rg1d_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(zg1d_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(psimesh_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(totals_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(igrid_sol,icount,mpi_Integer,0,MPI_COMM_WORLD,ierror)
    
        if(couplefrc==1)then
            icount=mpara_sep_sol
            call MPI_BCAST(Rsol_sep_reg,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
            call MPI_BCAST(Zsol_sep_reg,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
            call MPI_BCAST(Ssol_sep_reg,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        endif
    endif
    
    if(couplefrc==1)then
        call irreg2reg_interp_sol ! interpolate fields on regular grids of separatrix
        call reg2irreg_interp_sol ! interpolate fields on irregular grids of separatrix
        call orthogonal_grid_sep  ! perpendicularly trace field line for grids on sep
    endif
#ifndef GPU_UM
    !$acc update device(psimesh_sol,Ssol,mpara_sol,igrid_sol)
#endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    if(sfsol==0)then
        mgrid_sol=0
        mpsi_sol=0
        diag_region=1   ! safety guard
    elseif(sfcore==0)then
        mgrid_core=0
        mpsi_core=0
        diag_region=2   ! safety guard
    endif
    if(mype==0)write(gtcout,*)"grids_setup.F90, diag_region=",diag_region

! sum SOL and Core region grid number as mgrid and mpsi
    mgrid = mgrid_sol + mgrid_core
    mpsi = mpsi_sol + mpsi_core

    if(diag_region==1)then
        iflux = iflux_core
    elseif(diag_region==2)then
        iflux = mpsi_core + iflux_sol
    else
        if(mype==0)write(gtcout,*)"wrong diag_region,diag_region =",diag_region
        stop
    endif

    if(mype==0)then
        write(gtcout,101)mpsi,mpsi_core,mpsi_sol,mgrid,mgrid_core,mgrid_sol,iflux,iflux_core,iflux_sol
        call grids_dataout_mype0
    endif
101 format(' Simulation Grid Information',/,&
    1x, 'mpsi        =',i6/,&
    1x, 'mpsi_core    =',i6/,&
    1x, 'mpsi_sol     =',i6/,&
    1x, 'mgrid       =',i6/,&
    1x, 'mgrid_core  =',i6/,&
    1x, 'mgrid_sol   =',i6/,&
    1x, 'iflux       =',i6/,&
    1x, 'iflux_core  =',i6/,&
    1x, 'iflux_sol   =',i6)
end subroutine frc_mesh_create
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine core_initial_mype0
    use global_parameters, only:gtcout,pi
    use mapping_splines, only: rspcore_psi_theta, zspcore_psi_theta, rspcore_psi_s, zspcore_psi_s, msp_psi_core, aspcore_psi, &
                               lspcore, lstcore, del_psi_core, del_s_core, del_theta_core
    use spline_function, only: spline1d,spline2d,construct_spline1d,construct_spline2d
    use cylindrical_function, only: geotheta
    use equilibrium, only: psirz_Eq,&
                           lsr, lsz, spdR, spdZ, zmax
    use simulation_domain, only: R0core, R1core, psi1_core
    use read_in_parameters,only: min_r, para_diag_location_core

    integer :: i,j,ij
    real(lk) :: psidum

    write(gtcout,*)"CORE mesh"

    deltar_core=(R1core-R0core)/real(mpsi_core-1)

    zg1d_core(:)=0.0_lk
    do i=1,mpsi_core
       rg1d_core(i)=R0core+real(i-1)*deltar_core
       psimesh_core(i)=spline2d(rg1d_core(i),zg1d_core(i)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    enddo
    psimesh_core(mpsi_core)=psi1_core !Guarantee psi value on psi1_core

! Calculate the length of each closed field line, which is normalized by R0
    do i=1,mpsi_core
        totals_core(i)=spline1d(psimesh_core(i)-msp_psi_core(1),0,lspcore,del_psi_core,aspcore_psi,0)
    enddo
! Only allow equal-arc grid for core region, calculate the parallel grid size normalized by R0
    do i=1,mpsi_core
        deltas_core1d(i)=totals_core(i)/real(mpara_core(i)-1)
    enddo

    igrid_core(1)=0
    do i=2,mpsi_core
        igrid_core(i)=igrid_core(i-1)+mpara_core(i-1)
    enddo

! Score is the normalized field line length, normalized by length of local field line, [0,1]
    do i=1,mpsi_core
        do j=1,mpara_core(i)
            ij=igrid_core(i)+j
            Score(ij)=real(j-1)/real(mpara_core(i)-1)
        enddo
    enddo

    do i=1,mpsi_core
        do j=1,mpara_core(i)
            ij=igrid_core(i)+j
            Rcore(ij)=spline2d(psimesh_core(i)-msp_psi_core(1),Score(ij),0,lspcore,lstcore,del_psi_core,del_s_core,rspcore_psi_s,0,0)
            Zcore(ij)=spline2d(psimesh_core(i)-msp_psi_core(1),Score(ij),0,lspcore,lstcore,del_psi_core,del_s_core,zspcore_psi_s,0,0)
        enddo
    enddo

! quadratic spline with smooth of raw data has the risk that R < 0
! ensure that Rcore at the outest fluxsurface of core region is not negative   
    do j=1, mpara_core(mpsi_core)
        ij=igrid_core(mpsi_core) + j
            if(Rcore(ij) < min_r)then
                Rcore(ij)=0.0_lk  
            endif
    enddo

! geometric theta in core region   
    do i=1,mpsi_core
        do j=1,mpara_core(i)-1
            ij=igrid_core(i)+j
            Thetacore(ij)=geotheta(Rcore(ij),Zcore(ij))
        enddo
        Thetacore(igrid_core(i)+mpara_core(i))=2.0_lk*pi
    enddo

    write(gtcout,*) "Key information of core mesh at outer midplane:"
    write(gtcout,*)"parallel diag index/parallel grid number=",min(int(para_diag_location_core*real(mpara_core(1))) + 1, mpara_core(1)),"/",mpara_core(1)
    write(gtcout,110)
    do i=1,mpsi_core
        write(gtcout,120)i,rg1d_core(i),psimesh_core(i),totals_core(i),deltas_core1d(i),Rcore(igrid_core(i)+1),Zcore(igrid_core(i)+1),Thetacore(igrid_core(i)+1)
    enddo
110 format (1X,T4,"i",T8,"rg1d_core",T24,"psimesh_core",T40,"totals_core",T56,"deltas_core1d",T72,"Rcore(s=0)",T88,"Zcore(s=0)",T104,"Thetacore(s=0)")
120 format (1X,i4,7g16.8)
end subroutine core_initial_mype0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine sol_initial_mype0
    use global_parameters,only:gtcout,mype
    use mapping_splines, only: rspsol_psi_s, zspsol_psi_s, sspsol_psi_z, msp_psi_sol, aspsol_psi, &
                               lspsol, lszsol, del_psi_sol, del_z_sol, del_s_sol
    use spline_function, only: spline2d,spline1d
    use equilibrium, only: psirz_Eq, &
                           lsr, lsz, spdR, spdZ, zmax
    use simulation_domain, only: psi0_sol, R0sol, R1sol, &
                                 sfsol, sfcore, couplefrc
    use read_in_parameters,only: min_r, Zsol_domain0, para_diag_location_sol

    integer i,j,ij,ij1,j1,ij0,ii
    real(lk) psifrctmp,sdum

    write(gtcout,*)"SOL mesh"

    deltar_sol=(R1sol-R0sol)/real(mpsi_sol-1)
    if(mype==0)write(gtcout,*)"deltar_sol=",deltar_sol

    zg1d_sol=0.0_lk
    do i=1,mpsi_sol
        rg1d_sol(i)=R0sol+real(i-1)*deltar_sol
        psimesh_sol(i)=spline2d(rg1d_sol(i),zg1d_sol(i)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
    enddo
    psimesh_sol(1)=psi0_sol ! enforce the inner flux is on psi0_sol

    if(couplefrc==1)then
! set up separatrix grids in SOL region (irregular grid for grid-particle gather-scatter(need to use grid with same position on separatrix in core region), regular grid for finite difference) 
        do j=1,mpara_sep_sol
            ij=igrid_sol(1)+j
            Rsol(ij)=Rsol_sep_irreg(j)
            Zsol(ij)=Zsol_sep_irreg(j)
            Ssol(ij)=spline2d(psimesh_sol(1)-msp_psi_sol(1),Zsol(ij)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
        enddo

! notice that *_sep_reg should be same with *sol as following
        do j=1,mpara_sep_sol
            sdum=real(j-1)/real(mpara_sep_sol-1)
            Rsol_sep_reg(j)=spline2d(psimesh_sol(1)-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
            Zsol_sep_reg(j)=spline2d(psimesh_sol(1)-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)
            Ssol_sep_reg(j)=spline2d(psimesh_sol(1)-msp_psi_sol(1),Zsol_sep_reg(j)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
        enddo

    else

        do j=1,mpara_sol(1)
            ij=igrid_sol(1)+j
            sdum=real(j-1)/real(mpara_sol(1)-1)
            Rsol(ij)=spline2d(psimesh_sol(1)-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
            Zsol(ij)=spline2d(psimesh_sol(1)-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)
            Ssol(ij)=spline2d(psimesh_sol(1)-msp_psi_sol(1),Zsol(ij)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
        enddo

    endif


    do i=2,mpsi_sol
        do j=1,mpara_sol(i)
            ij=igrid_sol(i)+j
            sdum=real(j-1)/real(mpara_sol(i)-1)
            Rsol(ij)=spline2d(psimesh_sol(i)-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
            Zsol(ij)=spline2d(psimesh_sol(i)-msp_psi_sol(1),sdum,0,lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)
            Ssol(ij)=spline2d(psimesh_sol(i)-msp_psi_sol(1),Zsol(ij)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
        enddo
    enddo

! quadratic spline with smooth of raw data has the risk that R < 0
! ensure that Rsol at the innerest fluxsurface of SOL region is not negative   
    do j=1, mpara_sol(1)
        ij=igrid_sol(1) + j
        if(Rsol(ij) < min_r)then
            Rsol(ij)=0.0_lk  
        endif
    enddo
 
    do i=1,mpsi_sol
        totals_sol(i)=spline1d(psimesh_sol(i)-msp_psi_sol(1),0,lspsol,del_psi_sol,aspsol_psi,0)
    enddo

    do i=1,mpsi_sol
        do j=1,mpara_sol(i)
            ij=igrid_sol(i)+j
            ij1=min(ij+1,igrid_sol(i)+mpara_sol(i))
            deltas_sol(ij)=totals_sol(i)*(Ssol(ij1)-Ssol(ij))
        enddo
    enddo

    write(gtcout,*) "Key information of SOL mesh at outer midplane:"
    write(gtcout,*)"parallel diag index/parallel grid number=",&
                   min(int(para_diag_location_sol*real(mpara_sol(2))) + 1, mpara_sol(2)),"/",mpara_sol(2)
    write(gtcout,110)
    do i=1,mpsi_sol
        ij=igrid_sol(i) + int(0.5_lk*real(mpara_sol(i))) + 1 
        write(gtcout,120)i,rg1d_sol(i),psimesh_sol(i),totals_sol(i),&
                         deltas_sol(ij),Rsol(ij),Zsol(ij),Ssol(ij)
    enddo
110 format (1X,T4,"i",T8,"rg1d_sol",T24,"psimesh_sol",T40,"totals_sol",T56,"deltas_sol",T72,"Rsol",T88,"Zsol",T104,"Ssol")
120 format (1X,i4,7g16.8)
end subroutine sol_initial_mype0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine sep_initial
    use global_parameters, only: gtcout, mype
    use equilibrium, only: psi_sep
    use mapping_splines, only: sspsol_psi_z, msp_psi_sol, lspsol, lszsol, del_psi_sol, del_z_sol
    use spline_function, only: spline2d
    use read_in_parameters,only: min_r, Zsol_domain0, Zsol_domain1

    integer :: i,j,ij,ijsol,right_corner_index,left_corner_index, num_right_tmp, num_left_tmp
    real(lk) :: dz_right_tmp, dz_left_tmp, zleft_sol_tmp,zright_sol_tmp

    do j=1,mpara_core(mpsi_core)-1  ! Ghost grid in core region is excluded
        ij=igrid_core(mpsi_core)+j
        if(Rcore(ij) .le. min_r)then
! find the index of right corner grid in core region,right_corner_index is the beginning grid index which is not shared between sol and core regions
            right_corner_index=j 
            zright_sol_tmp=max(Zcore(ij-1),Zcore(ij))
            exit
        endif
    enddo
    do j=mpara_core(mpsi_core)-1,1,-1  ! Ghost grid in core region is excluded
        ij=igrid_core(mpsi_core)+j
        if(Rcore(ij) .le. min_r)then
! find the index of left corner grid in core region,left_corner_index is the end grid index which is not shared between sol and core regions
            left_corner_index=j  
            zleft_sol_tmp=min(Zcore(ij+1),Zcore(ij))
            exit
        endif
    enddo
 
    mpara_share=(mpara_core(mpsi_core)-1) - (left_corner_index-right_corner_index+1)
    allocate(mpara_core_share_index(mpara_share),mpara_sol_share_index(mpara_share))
    do i = 1,mpara_share/2
       mpara_core_share_index(i)=left_corner_index+i
    enddo
    do i = mpara_share/2+1, mpara_share
       mpara_core_share_index(i)=i-mpara_share/2
    enddo
    
    if(mpara_sol(1)-mpara_share .le. 0)then
        if(mype==0)write(gtcout,*)"increase mparasol_input to be further larger than mparacore_input"
        stop
    endif
   
    num_right_tmp=(mpara_sol(1)-mpara_share)/2
    num_left_tmp=(mpara_sol(1)-mpara_share)/2
    mpara_sep_sol=num_right_tmp+num_left_tmp+mpara_share

    allocate(Rsol_sep_irreg(mpara_sep_sol), Zsol_sep_irreg(mpara_sep_sol), Ssol_sep_irreg(mpara_sep_sol))
! Note that we need to devide the distance by actual grid number for interval calculation so that the grids are not in superposition at the corner
    dz_right_tmp=abs(Zsol_domain1-zright_sol_tmp)/real(num_right_tmp)
    dz_left_tmp=abs(zleft_sol_tmp-Zsol_domain0)/real(num_left_tmp)

    do i = 1, num_left_tmp
       Rsol_sep_irreg(i)=0.0_lk
       Zsol_sep_irreg(i)=Zsol_domain0 + dz_left_tmp*real(i-1)
    enddo
    do i = num_left_tmp+1, num_left_tmp+mpara_share
       Rsol_sep_irreg(i)=Rcore(igrid_core(mpsi_core)+mpara_core_share_index(i-num_left_tmp))
       Zsol_sep_irreg(i)=Zcore(igrid_core(mpsi_core)+mpara_core_share_index(i-num_left_tmp))
       mpara_sol_share_index(i-num_left_tmp)=i
    enddo
    do i = num_left_tmp+mpara_share+1, mpara_sep_sol
       Rsol_sep_irreg(i)=0.0_lk
       Zsol_sep_irreg(i)=min(Zsol_domain1, zright_sol_tmp+dz_right_tmp*real(i-num_left_tmp-mpara_share))
    enddo

    do i = 1, mpara_sep_sol
       Ssol_sep_irreg(i)=spline2d(psi_sep-msp_psi_sol(1),Zsol_sep_irreg(i)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
    enddo

    mpara_sol(1)=mpara_sep_sol

    if(mype==0)write(gtcout,*)"mpara_sol_share_index=",mpara_sol_share_index

end subroutine sep_initial
!=================================================================
! find the interpolation weight and index of neighbour irregular aligned grids for regular aligned grids on separatrix
! Get the value for regular grids from irregular grids
!=================================================================
subroutine irreg2reg_interp_sol
    use global_parameters

    integer :: j,jt(1)
    real(lk) :: ds

! allocate memory
    allocate(irreg2reg_index_sol(mpara_sep_sol))
    allocate(irreg2reg_weight_sol(mpara_sep_sol))

    irreg2reg_index_sol = 0
    irreg2reg_weight_sol = 0.0_lk

! find the lower neighbour irregular grid index for regular grid
    do j=1,mpara_sep_sol
        jt=minloc(abs(Ssol_sep_reg(j)-Ssol_sep_irreg))  ! jt is the index of Ssol_irregular closest to the target grid

        ds=Ssol_sep_reg(j)-Ssol_sep_irreg(jt(1))
        if(ds>0.0_lk)then
            irreg2reg_index_sol(j)=min(jt(1),mpara_sep_sol-1)
        else
            irreg2reg_index_sol(j)=max(jt(1)-1,1)
        endif
! higher neighbour irregular grid weight for regular grid
        irreg2reg_weight_sol(j)=(Ssol_sep_reg(j)-Ssol_sep_irreg(irreg2reg_index_sol(j)))/(Ssol_sep_irreg(irreg2reg_index_sol(j)+1)-Ssol_sep_irreg(irreg2reg_index_sol(j)))
    enddo

end subroutine irreg2reg_interp_sol
!=================================================================
! find the interpolation weight and index of neighbour regular aligned grids for irregular aligned grids on separatrix
! Get the value for irregular grids from regular grids
!=================================================================
subroutine reg2irreg_interp_sol
    use global_parameters

    integer :: j,jt(1)
    real(lk) :: ds

! allocate memory
    allocate(reg2irreg_index_sol(mpara_sep_sol))
    allocate(reg2irreg_weight_sol(mpara_sep_sol))

    reg2irreg_index_sol = 0
    reg2irreg_weight_sol = 0.0_lk

! find the lower neighbour irregular grid index for regular grid
    do j=1,mpara_sep_sol
        jt=minloc(abs(Ssol_sep_irreg(j)-Ssol_sep_reg))  ! jt is the index of Ssol_irregular closest to the target grid

        ds=Ssol_sep_irreg(j)-Ssol_sep_reg(jt(1))
        if(ds>0.0_lk)then
            reg2irreg_index_sol(j)=min(jt(1),mpara_sep_sol-1)
        else
            reg2irreg_index_sol(j)=max(jt(1)-1,1)
        endif
! higher neighbour irregular grid weight for regular grid
        reg2irreg_weight_sol(j)=(Ssol_sep_irreg(j)-Ssol_sep_reg(reg2irreg_index_sol(j)))/(Ssol_sep_reg(reg2irreg_index_sol(j)+1)-Ssol_sep_reg(reg2irreg_index_sol(j)))
    enddo

end subroutine reg2irreg_interp_sol
!==================================================================
! find the orthogonal neighbour grids in SOL and core regions for grids on separatrix
!==================================================================
subroutine orthogonal_grid_sep
    use global_parameters,only:gtcout,mype,pi
    use read_in_parameters, only: Zsol_domain0
    use spline_function, only: spline2d
    use cylindrical_function, only: geotheta
    use trace_fieldlines, only: trace_perp
    use mapping_splines, only: sspcore_psi_theta, sspsol_psi_z, msp_psi_core, msp_psi_sol, &
                               lspsol,lszsol,del_psi_sol,del_z_sol, &
                               lspcore,lstcore,del_psi_core,del_theta_core
    implicit none


    integer :: i,j,k,ij,trace_grid_num,jt(1),itrace,&
               maxgrid_num_sol_tmp, maxgrid_num_core_tmp
    real(lk) :: brtmp,bztmp,btmp,ds,bperp_r,bperp_z,pdum,tdum,&
                delta_pdum,deltascoretmp,deltassoltmp
    real(lk),dimension(:,:),allocatable :: rtrace_sol_tmp,ztrace_sol_tmp,rtrace_core_tmp,ztrace_core_tmp
    integer,dimension(:),allocatable :: step_num_sol_tmp,step_num_core_tmp
    real(lk),dimension(:),allocatable :: Ssoltmp,Scoretmp
 

    trace_grid_num=10000
! Towards SOL region
    deltassoltmp=(rg1d_sol(2)-rg1d_sol(1))/100.0_lk
    allocate(rtrace_sol_tmp(mpara_share,trace_grid_num),ztrace_sol_tmp(mpara_share,trace_grid_num))
    allocate(step_num_sol_tmp(mpara_share))
    rtrace_sol_tmp=0.0_lk
    ztrace_sol_tmp=0.0_lk
    step_num_sol_tmp=0  ! Must define initial value since the loop does not contain 1 and mpara_share

    do i=2,mpara_share-1
        call trace_perp(1,trace_grid_num,Rsol_sep_irreg(mpara_sol_share_index(i)),Zsol_sep_irreg(mpara_sol_share_index(i)),psimesh_sol(2),deltassoltmp,step_num_sol_tmp(i),rtrace_sol_tmp(i,:),ztrace_sol_tmp(i,:))
    enddo

!Towards core region
    deltascoretmp=(rg1d_core(mpsi_core)-rg1d_core(mpsi_core-1))/100.0_lk
    allocate(rtrace_core_tmp(mpara_share,trace_grid_num),ztrace_core_tmp(mpara_share,trace_grid_num))
    allocate(step_num_core_tmp(mpara_share))
    step_num_core_tmp=0  ! Must define initial value since the loop does not contain 1 and mpara_share
    rtrace_core_tmp=0.0_lk
    ztrace_core_tmp=0.0_lk

    do i=2,mpara_share-1
        call trace_perp(2,trace_grid_num,Rsol_sep_irreg(mpara_sol_share_index(i)),Zsol_sep_irreg(mpara_sol_share_index(i)),psimesh_core(mpsi_core-1),deltascoretmp,step_num_core_tmp(i),rtrace_core_tmp(i,:),ztrace_core_tmp(i,:))
    enddo


    allocate(Rsol_share(mpara_share),Zsol_share(mpara_share),Ssol_share(mpara_share),&
             Rcore_share(mpara_share),Zcore_share(mpara_share),Score_share(mpara_share),&
             dcore_orthogonal(mpara_share),dsol_orthogonal(mpara_share))
 
    Rsol_share=0.0_lk
    Zsol_share=0.0_lk
    Ssol_share=0.0_lk

    Rcore_share=0.0_lk
    Zcore_share=0.0_lk
    Score_share=0.0_lk

    dcore_orthogonal=0.0_lk
    dsol_orthogonal=0.0_lk
! find the orthogonal grid position in SOL region
    do i=2,mpara_share-1
        j=step_num_sol_tmp(i)
        Rsol_share(i)=rtrace_sol_tmp(i,j)
        Zsol_share(i)=ztrace_sol_tmp(i,j)
        Ssol_share(i)=spline2d(psimesh_sol(2)-msp_psi_sol(1),Zsol_share(i)-Zsol_domain0,0,lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
        dsol_orthogonal(i)=deltassoltmp*real(step_num_sol_tmp(i))
    enddo

! find the orthogonal grid position in core region
    do i=2,mpara_share-1
        j=step_num_core_tmp(i)
        if(mpara_core_share_index(i)==1)then
            Rcore_share(i)=rg1d_core(mpsi_core-1)   ! enforce the orthogonal grid is on outermidplane
            Zcore_share(i)=0.0_lk
            Score_share(i)=0.0_lk
        else
            Rcore_share(i)=rtrace_core_tmp(i,j)
            Zcore_share(i)=ztrace_core_tmp(i,j)
            tdum=geotheta(Rcore_share(i),Zcore_share(i))
            Score_share(i)=spline2d(psimesh_core(mpsi_core-1)-msp_psi_core(1),tdum,0,lspcore,lstcore,del_psi_core,del_theta_core,sspcore_psi_theta,0,0)
        endif
        dcore_orthogonal(i)=deltascoretmp*real(step_num_core_tmp(i))
    enddo
!========================================================================
!   calculate the interpolation grid index and weight
!========================================================================
    allocate(interp_weight_sol_share(mpara_share),interp_index_sol_share(mpara_share),interp_weight_core_share(mpara_share),interp_index_core_share(mpara_share))

    interp_weight_sol_share=0.0_lk
    interp_weight_core_share=0.0_lk
    interp_index_sol_share=0
    interp_index_core_share=0


    allocate(Ssoltmp(mpara_sol(2)),Scoretmp(mpara_core(mpsi_core-1)))
    do j=1,mpara_sol(2)
        ij=igrid_sol(2)+j
        Ssoltmp(j)=Ssol(ij)
    enddo

    do j=1,mpara_core(mpsi_core-1)
        ij=igrid_core(mpsi_core-1)+j
        Scoretmp(j)=Score(ij)
    enddo
! SOL region interpolation
    do j=2,mpara_share-1
        jt=minloc(abs(Ssol_share(j)-Ssoltmp))
        ds=Ssol_share(j)-Ssoltmp(jt(1))
        if(ds>0.0_lk)then
            interp_index_sol_share(j)=min(jt(1),mpara_sol(2)-1)
        else
            interp_index_sol_share(j)=max(jt(1)-1,1)
        endif
!weight for upper parallel grid
        interp_weight_sol_share(j)=(Ssol_share(j)-Ssoltmp(interp_index_sol_share(j)))/(Ssoltmp(interp_index_sol_share(j)+1)-Ssoltmp(interp_index_sol_share(j)))
    enddo

! core region interpolation
    do j=2,mpara_share-1
        jt=minloc(abs(Score_share(j)-Scoretmp))
        ds=Score_share(j)-Scoretmp(jt(1))
        if(ds>0.0_lk)then
            interp_index_core_share(j)=min(jt(1),mpara_core(mpsi_core-1)-1)
        else
            interp_index_core_share(j)=max(jt(1)-1,1)
        endif
!weight for upper parallel grid
        interp_weight_core_share(j)=(Score_share(j)-Scoretmp(interp_index_core_share(j)))/(Scoretmp(interp_index_core_share(j)+1)-Scoretmp(interp_index_core_share(j)))
    enddo

    maxgrid_num_sol_tmp=maxval(step_num_sol_tmp)
    maxgrid_num_core_tmp=maxval(step_num_core_tmp)
    if(mype==0)then
        write(gtcout,*)"maxgrid_num_core_tmp",maxgrid_num_core_tmp
        write(gtcout,*)"maxgrid_num_sol_tmp",maxgrid_num_sol_tmp
    endif

    if(mype==0)then
        itrace=1989
        open(itrace,file='trace_perp.out',status='replace')
        write(itrace,101)mpara_share,maxgrid_num_core_tmp,maxgrid_num_sol_tmp
        write(itrace,101)step_num_sol_tmp, step_num_core_tmp
        write(itrace,102)rtrace_sol_tmp(1:mpara_share,1:maxgrid_num_sol_tmp),ztrace_sol_tmp(1:mpara_share,1:maxgrid_num_sol_tmp)
        write(itrace,102)rtrace_core_tmp(1:mpara_share,1:maxgrid_num_core_tmp),ztrace_core_tmp(1:mpara_share,1:maxgrid_num_core_tmp)
        close(itrace)
    endif

    deallocate(Ssoltmp,Scoretmp)
    deallocate(rtrace_sol_tmp,ztrace_sol_tmp,rtrace_core_tmp,ztrace_core_tmp)
    deallocate(step_num_sol_tmp,step_num_core_tmp)

101 format(i6)
102 format(e13.6)
end subroutine orthogonal_grid_sep
!===================================================================
subroutine grids_dataout_mype0
    use read_in_parameters, only: numerfrc
    use mapping_splines, only: rspcore_psi_theta, zspcore_psi_theta, &
                               rspcore_psi_s, zspcore_psi_s, &
                               msp_psi_core, lspcore, lstcore, del_psi_core, &
                               del_s_core, del_theta_core
    use spline_function, only: spline2d
    use equilibrium, only: psirz_Eq, lsr, lsz, spdR, spdZ, zmax
    use simulation_domain, only: sfcore
    implicit none
  
    integer igrid, i, j, ij
    real(lk) psidum
    real(lk),dimension(:),allocatable :: r_psi_theta, z_psi_theta, r_psi_s, z_psi_s
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!             Spline check by using the simulation grids
!    calculate the \theta_s value for each simulation grids by using 
!                             1d spline function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if(sfcore==1)then
        allocate(r_psi_theta(mgrid_core),z_psi_theta(mgrid_core))
        
        do i=1,mpsi_core
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                psidum=spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
                r_psi_theta(ij)=spline2d(psidum-msp_psi_core(1),Thetacore(ij),0,lspcore,lstcore,del_psi_core,del_theta_core,rspcore_psi_theta,0,0)
                z_psi_theta(ij)=spline2d(psidum-msp_psi_core(1),Thetacore(ij),0,lspcore,lstcore,del_psi_core,del_theta_core,zspcore_psi_theta,0,0)
        
            enddo
        enddo
        
        allocate(r_psi_s(mgrid_core),z_psi_s(mgrid_core))
        do i=1,mpsi_core
            do j=1,mpara_core(i)
                ij=igrid_core(i)+j
                psidum=spline2d(Rcore(ij),Zcore(ij)+zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
                r_psi_s(ij)=spline2d(psidum-msp_psi_core(1),Score(ij),0,lspcore,lstcore,del_psi_core,del_s_core,rspcore_psi_s,0,0)
                z_psi_s(ij)=spline2d(psidum-msp_psi_core(1),Score(ij),0,lspcore,lstcore,del_psi_core,del_s_core,zspcore_psi_s,0,0)
            enddo
        enddo
    endif

    igrid=100
    open(igrid,file='grids.out',status='replace')
    if(numerfrc==1)then
        write(igrid,101)mpsi_core,mpara_core(1),mgrid_core, &
                        mpsi_sol,mpara_sol(1),mgrid_sol, &
                        mpara_sep_sol,mpara_share
        write(igrid,102)Rcore, Zcore, Score, Thetacore, &
                        Rsol, Zsol, Ssol, &
                        Rsol_sep_irreg, Zsol_sep_irreg, Ssol_sep_irreg, &
                        Rsol_sep_reg, Zsol_sep_reg, Ssol_sep_reg, &
                        Rsol_share, Zsol_share, Ssol_share, &
                        Rcore_share, Zcore_share, Score_share, &
                        dcore_orthogonal, dsol_orthogonal
        write(igrid,102)rg1d_core, zg1d_core, psimesh_core, &
                        rg1d_sol, zg1d_sol, psimesh_sol
        write(igrid,102)totals_core,deltas_core1d, &
                        totals_sol,deltas_sol
                  
        write(igrid,102)r_psi_theta,z_psi_theta,r_psi_s,z_psi_s
        deallocate(r_psi_theta, z_psi_theta, r_psi_s, z_psi_s)

    elseif(numerfrc==2)then
        write(igrid,101)mpsi_core,mpara_core(1),mgrid_core
        write(igrid,102)Rcore, Zcore, Score, Thetacore
        write(igrid,102)rg1d_core, zg1d_core, psimesh_core
        write(igrid,102)totals_core,deltas_core1d
                  
        write(igrid,102)r_psi_theta,z_psi_theta,r_psi_s,z_psi_s
        deallocate(r_psi_theta, z_psi_theta, r_psi_s, z_psi_s)

    elseif(numerfrc==3 .or. numerfrc==4)then
        write(igrid,101)mpsi_sol,mpara_sol(1),mgrid_sol
        write(igrid,102)Rsol, Zsol, Ssol
        write(igrid,102)rg1d_sol, zg1d_sol, psimesh_sol
        write(igrid,102)totals_sol,deltas_sol
           
    endif                 

    close(igrid)

101 format(i6)
102 format(e13.6)

end subroutine grids_dataout_mype0
!===================================================================
!   interplate the field from irregular grids on separatrix 
!   to the regular grids for SOL and core coupled simulation
!===================================================================
subroutine convert_irreg2reg(farray2d,farray1d)
    use global_parameters, only: mype, gtcout
    use simulation_domain,only: couplefrc
    real(lk), dimension(mgrid), intent(inout) :: farray2d
    real(lk), dimension(mpara_sep_sol), intent(inout), optional :: farray1d

    integer i
    real(lk) farraytmp(mpara_sep_sol)

    if(couplefrc==0)then
        if(mype==0)write(gtcout,*)"convert_irreg2reg is not applicable for couplefrc=0 case, it requires SOL and core region coupling"
        stop
    endif


    do i=1,mpara_sep_sol
        farraytmp(i)=(1.0_lk-irreg2reg_weight_sol(i))*farray2d(mgrid_core+igrid_sol(1)+irreg2reg_index_sol(i))&
                     +irreg2reg_weight_sol(i)*farray2d(mgrid_core+igrid_sol(1)+irreg2reg_index_sol(i)+1)
    enddo

    farray2d(mgrid_core+igrid_sol(1)+1 : mgrid_core+igrid_sol(1)+mpara_sep_sol) = farraytmp(1:mpara_sep_sol)

    if(present(farray1d))then
        farray1d=farraytmp
    endif
end subroutine convert_irreg2reg
!===================================================================
!   interplate the field from regular grids on separatrix 
!   to the irregular grids for SOL and core coupled simulation
!===================================================================
subroutine convert_reg2irreg(farray2d,farray1d)
    use global_parameters, only: mype, gtcout
    use simulation_domain,only: couplefrc
    real(lk), dimension(mgrid), intent(inout) :: farray2d
    real(lk), dimension(mpara_sep_sol), intent(inout), optional :: farray1d

    integer i
    real(lk) farraytmp(mpara_sep_sol)

    if(couplefrc==0)then
        if(mype==0)write(gtcout,*)"convert_irreg2reg is not applicable for couplefrc=0 case, it requires SOL and core region coupling"
        stop
    endif

    do i=1,mpara_sep_sol
        farraytmp(i)=(1.0_lk-reg2irreg_weight_sol(i))*farray2d(mgrid_core+igrid_sol(1)+reg2irreg_index_sol(i))&
                     +reg2irreg_weight_sol(i)*farray2d(mgrid_core+igrid_sol(1)+reg2irreg_index_sol(i)+1)
    enddo

    farray2d(mgrid_core+igrid_sol(1)+1 : mgrid_core+igrid_sol(1)+mpara_sep_sol) = farraytmp(1:mpara_sep_sol)

    if(present(farray1d))then
        farray1d=farraytmp
    endif
end subroutine convert_reg2irreg
!===================================================================
!   accumulate particle moments on the shared grids at separatrix
!===================================================================
subroutine particle_field_accumulation_sep(farray2d)
    use global_parameters, only: mype, gtcout
    use simulation_domain, only: couplefrc
 
    real(lk), dimension(0:1,mgrid), intent(inout) :: farray2d
    
    integer i, ijcore, ijsol

    if(couplefrc==0)then
        if(mype==0)write(gtcout,*)"particle_field_accumulation_sep is not applicable for couplefrc=0 case, it requires SOL and core region coupling"
        stop
    endif
! notice that mpara_core_share_index skip the end ghost grid in core
    do i=1,mpara_share
        ijcore = igrid_core(mpsi_core) + mpara_core_share_index(i)
        ijsol = mgrid_core + igrid_sol(1) + mpara_sol_share_index(i)

        farray2d(:,ijcore) = farray2d(:,ijcore) + farray2d(:,ijsol)
        farray2d(:,ijsol) = farray2d(:,ijcore)
    enddo
end subroutine particle_field_accumulation_sep
!===================================================================
!   core region parallel periodicity
!===================================================================
subroutine core_para_periodicity(farray2d)
    use global_parameters, only: mype, gtcout
    use simulation_domain, only: sfcore

    integer i, ij_begin, ij_end
    real(lk), dimension(0:1,mgrid), intent(inout) :: farray2d

    if(sfcore==0)then
        if(mype==0)write(gtcout,*)"core_para_periodicity is not applicable for sfcore=0 case"
        stop
    endif

    do i=1,mpsi_core
        ij_begin = igrid_core(i) + 1
        ij_end = igrid_core(i) + mpara_core(i)
        farray2d(:,ij_begin) = farray2d(:,ij_begin) + farray2d(:,ij_end)
        farray2d(:,ij_end) = farray2d(:,ij_begin)
    enddo

end subroutine core_para_periodicity

subroutine sol_para_period(farray2d)
    use simulation_domain, only: sfsol
    use global_parameters, only: mype, gtcout
    integer i, ij_begin, ij_end
    real(lk), dimension(0:1,mgrid), intent(inout) :: farray2d

    if(sfsol==0)then
        if(mype==0)write(gtcout,*)"sol_para_periodicity is not applicable for sfcore=0 case"
        stop
    endif

    do i=1,mpsi_sol
        ij_begin = mgrid_core + igrid_sol(i) + 1
        ij_end = mgrid_core + igrid_sol(i) + mpara_sol(i)
        farray2d(:,ij_begin) = farray2d(:,ij_begin) + farray2d(:,ij_end)
        farray2d(:,ij_end) = farray2d(:,ij_begin)
    enddo
end subroutine sol_para_period
end module grids_setup

subroutine radial_decay_1d(array)
use grids_setup
use read_in_parameters
use precision
implicit none

real(lk), dimension(mpsi), intent(inout) :: array

integer i

    do i=1,BC_nbound0
       array(i) = (i-1)*1.0_lk/BC_nbound0*array(i)
    enddo
    do i=mpsi-BC_nbound1+1, mpsi
       array(i) = (mpsi-i)*1.0_lk/BC_nbound1*array(i)
    enddo

end subroutine radial_decay_1d

subroutine radial_decay_2d(array)
use grids_setup
use read_in_parameters
use simulation_domain
use precision
implicit none

real(lk), dimension(mgrid), intent(inout) :: array

integer i,j,ij
    if(sfcore==1)then
      do i=1,BC_nbound0
         do j=1,mpara_core(i)
             ij = igrid_core(i) + j
             array(ij) = (i-1)*1.0_lk/BC_nbound0*array(ij)
         enddo
      enddo
      if(couplefrc==0)then
        do i=mpsi_core-BC_nbound1+1, mpsi_core
           do j=1,mpara_core(i)
              ij = igrid_core(i) + j
              array(ij) = (mpsi_core-i)*1.0_lk/BC_nbound1*array(ij)
           enddo
        enddo
      endif
    endif

    if(sfsol==1)then
    if(couplefrc==0)then
      do i=1,BC_nbound0
         do j=1,mpara_sol(i)
            ij = mgrid_core + igrid_sol(i) + j
            array(ij) = (i-1)*1.0_lk/BC_nbound0*array(ij)
         enddo
      enddo
    endif
    do i=mpsi-BC_nbound1+1, mpsi
       do j=1,mpara_sol(i)
          ij = mgrid_core + igrid_sol(i) + j
          array(ij) = (mpsi_sol-i)*1.0_lk/BC_nbound1*array(ij)
       enddo
    enddo
    endif

end subroutine radial_decay_2d

