module mapping_splines
use precision
implicit none

public frc_coordinate_spline, er_init, &
       lstcore, lspcore, lszsol, lspsol, &
       del_s_core, del_theta_core, del_psi_core, &
       del_s_sol, del_z_sol, del_psi_sol, &
       msp_psi_core, msp_psi_sol, &  ! psi value on each mapping spline surface
       aspcore_psi, aspsol_psi, &  ! arc length of each field line as a function of psi
       rspcore_psi_theta, zspcore_psi_theta, sspcore_psi_theta, &
       rspcore_psi_s, zspcore_psi_s, &
       sspsol_psi_z, rspsol_psi_s, zspsol_psi_s, &
       phip00_sp, lsp_er, dsp_er   ! Er
private


! grid number in parallel direction for mapping_splines: lstcore, lszsol
! grid number in perpendicular direction for mapping_splines: lspcore, lspsol, lsp_er
integer, parameter :: lstcore=300, lspcore=100, lszsol=300, lspsol=100, lsp_er=100

real(lk) del_s_core, del_theta_core, del_psi_core, &
         del_s_sol, del_z_sol, del_psi_sol
real(lk), dimension(:), allocatable :: msp_psi_core, msp_psi_sol

real(lk),dimension(:,:),allocatable :: aspcore_psi,aspsol_psi

real(lk),dimension(:,:,:),allocatable :: rspcore_psi_theta, zspcore_psi_theta, sspcore_psi_theta, &
                                         rspcore_psi_s, zspcore_psi_s, &
                                         sspsol_psi_z, rspsol_psi_s, zspsol_psi_s
!Er
real(lk) dsp_er
real(lk), dimension(:,:), allocatable :: phip00_sp

save
#ifndef GPU_UM
!$acc declare create(msp_psi_core, msp_psi_sol)
!$acc declare create(rspcore_psi_theta, zspcore_psi_theta, sspcore_psi_theta)
!$acc declare create(rspcore_psi_s,zspcore_psi_s)
!$acc declare create(rspsol_psi_s, zspsol_psi_s, sspsol_psi_z)
#endif

contains
!======================================================================
subroutine er_init
    use global_parameters, only: mype,gtcout
    use simulation_domain, only: psi0_sol, psi1_sol, psi0_core, psi1_core, sfcore, sfsol
    use spline_function, only: construct_spline1d,spline2d
    use equilibrium, only: psirz_Eq,zmax,lsr,lsz,spdR,spdZ
    use read_in_parameters, only: numerfrc,Rdomain0,Rdomain1,gamma0,diag_region, &
                                  mpsicore_input,mpsisol_input,iflux_core,iflux_sol
!    use grids_setup, only: iflux,mpsi
    integer i, icount, ierror, iflux, mpsi
    real(lk) psitmp,R_diag,psi_diag
   
    allocate(phip00_sp(3,lsp_er))
    phip00_sp=0.0_lk
    
    if (numerfrc==1)then
        mpsi = mpsicore_input+mpsisol_input
        if (diag_region==1)then
            iflux = iflux_core
        elseif (diag_region==2)then
            iflux = mpsicore_input+iflux_sol
        endif
    elseif (numerfrc==2)then
        mpsi = mpsicore_input
        iflux = iflux_core
    elseif (numerfrc==3)then
        mpsi = mpsisol_input
        iflux = iflux_sol
    endif
    R_diag = Rdomain0+(Rdomain1-Rdomain0)*1.0_lk*iflux/mpsi
    psi_diag=spline2d(R_diag,zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

    if(mype==0)then
      write(gtcout,*)'R_diag, psi_diag, gamma0 = '
      write(gtcout,*)R_diag,psi_diag,gamma0
      if(numerfrc==3)then
          dsp_er=(psi1_sol-psi0_sol)/real(lsp_er-1)
          do i=1,lsp_er
            psitmp=real(i-1)*dsp_er+psi0_sol
! initially we estimate R^2*B=6.0_lk, this will be rescaled in grids_eqdata.F90
            phip00_sp(1,i)=(psitmp-psi_diag)*gamma0/6.0_lk
          enddo
      elseif(numerfrc==2)then
          dsp_er=(psi1_core-psi0_core)/real(lsp_er-1)
          do i=1,lsp_er
            psitmp=real(i-1)*dsp_er+psi0_core
            phip00_sp(1,i)=(psitmp-psi_diag)*gamma0/6.0_lk
          enddo
      elseif(numerfrc==1)then
          dsp_er=(psi1_sol-psi0_core)/real(lsp_er-1)
          do i=1,lsp_er
            psitmp=real(i-1)*dsp_er+psi0_core
            phip00_sp(1,i)=(psitmp-psi_diag)*gamma0/6.0_lk
          enddo
      endif

      call construct_spline1d(lsp_er,dsp_er,phip00_sp,0)
    endif
  
    icount=1
    call MPI_BCAST(dsp_er,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
    icount=3*lsp_er
    call MPI_BCAST(phip00_sp,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
end subroutine er_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine frc_coordinate_spline
    use global_parameters,only:gtcout,mype
    use spline_function, only: spline2d
    use simulation_domain, only: sfcore, sfsol

    integer :: mtest,icount,ierror,istatus(MPI_STATUS_SIZE)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!       Core region spline creation
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if(sfcore==1)then

        allocate(msp_psi_core(lspcore))
        allocate(aspcore_psi(3,lspcore))
        allocate(rspcore_psi_theta(9,lspcore,lstcore),zspcore_psi_theta(9,lspcore,lstcore),sspcore_psi_theta(9,lspcore,lstcore))
        allocate(rspcore_psi_s(9,lspcore,lstcore),zspcore_psi_s(9,lspcore,lstcore))
        msp_psi_core=0.0_lk
        aspcore_psi=0.0_lk
        rspcore_psi_theta=0.0_lk
        zspcore_psi_theta=0.0_lk
        sspcore_psi_theta=0.0_lk
        rspcore_psi_s=0.0_lk
        zspcore_psi_s=0.0_lk

        if(mype==0)then
            call frc_coordinate_splinecore
        endif
        call MPI_Barrier(MPI_COMM_WORLD,ierror)

! Broadcast the spline function for core region to each MPIs
        icount=1
        call MPI_BCAST(del_theta_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        icount=1
        call MPI_BCAST(del_psi_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        icount=1
        call MPI_BCAST(del_s_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
     
        icount=lspcore
        call MPI_BCAST(msp_psi_core,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
     
        icount=3*lspcore
        call MPI_BCAST(aspcore_psi,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
     
        icount=9*lspcore*lstcore
        call MPI_BCAST(rspcore_psi_theta,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(zspcore_psi_theta,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(sspcore_psi_theta,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(rspcore_psi_s,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(zspcore_psi_s,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
#ifndef GPU_UM
        !$acc update device(rspcore_psi_theta, zspcore_psi_theta, sspcore_psi_theta)
        !$acc update device(rspcore_psi_s,zspcore_psi_s, msp_psi_core)
#endif
    
    endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!       SOL region spline creation
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if(sfsol==1)then
    
        allocate(msp_psi_sol(lspsol))
        msp_psi_sol=0.0_lk
    
        allocate(aspsol_psi(3,lspsol),sspsol_psi_z(9,lspsol,lszsol),rspsol_psi_s(9,lspsol,lszsol),zspsol_psi_s(9,lspsol,lszsol))
        aspsol_psi=0.0_lk
        sspsol_psi_z=0.0_lk
        rspsol_psi_s=0.0_lk
        zspsol_psi_s=0.0_lk
    
        if(mype==0)then
            call frc_coordinate_splinesol
        endif
        call MPI_Barrier(MPI_COMM_WORLD,ierror)
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        icount=1
        call MPI_BCAST(del_psi_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        icount=1
        call MPI_BCAST(del_z_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        icount=1
        call MPI_BCAST(del_s_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
     
        icount=lspsol
        call MPI_BCAST(msp_psi_sol,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
     
        icount=3*lspsol
        call MPI_BCAST(aspsol_psi,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
     
        icount=9*lspsol*lszsol
        call MPI_BCAST(rspsol_psi_s,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(zspsol_psi_s,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
        call MPI_BCAST(sspsol_psi_z,icount,mpi_Rsize,0,MPI_COMM_WORLD,ierror)
#ifndef GPU_UM
        !$acc update device(rspsol_psi_s, zspsol_psi_s, sspsol_psi_z,msp_psi_sol)
#endif
    
    endif

end subroutine frc_coordinate_spline
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine frc_coordinate_splinecore
    use global_parameters,only:gtcout,mype,pi
    use spline_function, only: spline1d,spline2d,sppsiR1d,construct_spline1d,construct_spline2d
    use equilibrium, only: r_Eq, z_Eq, psirz_Eq, &
                           psi_sep, lsr, lsz, spdR, spdZ, zmax
    use simulation_domain, only: psi0_core, psi1_core
    use read_in_parameters, only: Zsol_domain0, Zsol_domain1
    use trace_fieldlines, only: trace, trace_sep
  
    integer :: i,j,ij,mtest,ngtmp,j1,k,ik
    real(lk) psifrctmp,deltas,tdum,tdum1,tdum2,z1sep,z2sep,asep
    integer,dimension(:),allocatable :: mpara0tmp
    real(lk),dimension(:,:),allocatable :: Rfasp,Zfasp,tfasp
    real(lk),dimension(:),allocatable :: z0tmp,r0tmp,dscoretmp,n_dscoretmp
  
    msp_psi_core(1)=psi0_core
    msp_psi_core(lspcore)=psi1_core
  
    del_theta_core=2.0_lk*pi/real(lstcore-1)
    del_psi_core=(msp_psi_core(lspcore)-msp_psi_core(1))/real(lspcore-1)
  
    allocate(z0tmp(lspcore),r0tmp(lspcore),mpara0tmp(lspcore))
    z0tmp=0.0_lk
    r0tmp=0.0_lk
    mpara0tmp=0
  
    do i=1,lspcore
        msp_psi_core(i)=msp_psi_core(1)+del_psi_core*real(i-1)
        call sppsiR1d(1,lsr,lsz,spdR,spdZ,r_Eq,z_Eq,psirz_Eq,z0tmp(i)+zmax,msp_psi_core(i),r0tmp(i),psifrctmp)
    enddo
    msp_psi_core(lspcore)=psi1_core !Guarantee psi value on the separatrix
  
    deltas=0.001_lk
    ngtmp=int(3.0*zmax/deltas)
  
    allocate(Rfasp(lspcore,ngtmp),Zfasp(lspcore,ngtmp),tfasp(lspcore,ngtmp))
    Rfasp=0.0_lk
    Zfasp=0.0_lk
    tfasp=0.0_lk
  
    do k=1,lspcore-1
        call trace(1,1,ngtmp,r0tmp(k),z0tmp(k),msp_psi_core(k),deltas, Zsol_domain0, Zsol_domain1, mpara0tmp(k),Rfasp(k,:),Zfasp(k,:),tfasp(k,:))
    enddo
    if(psi1_core .lt. psi_sep)then
        call trace(1,1,ngtmp,r0tmp(lspcore),z0tmp(lspcore),msp_psi_core(lspcore),deltas, Zsol_domain0, Zsol_domain1, mpara0tmp(lspcore),Rfasp(lspcore,:),Zfasp(lspcore,:),tfasp(lspcore,:))
    else
        call trace_sep(ngtmp,r0tmp(lspcore),z0tmp(lspcore),deltas,mpara0tmp(lspcore),Rfasp(lspcore,:),Zfasp(lspcore,:),tfasp(lspcore,:),z1sep,z2sep,asep)
    endif
  
    do i=1,mpara0tmp(lspcore)
        if(Rfasp(lspcore,i)<0.0_lk)then
            if(mype==0)write(gtcout,*)"i"
        endif
    enddo
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          R(\psi,\theta), Z(\psi,\theta),S(\psi,\theta),s(\theta)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
    do i=1,lspcore
        do k=1,lstcore-1
            do j=1,mpara0tmp(i)
  
                tdum1=(tfasp(i,j)-real(k-1)*del_theta_core)
                tdum2=(tfasp(i,min(j+1,mpara0tmp(i)))-real(k-1)*del_theta_core)
                if(tdum2*tdum1 .le. 0 .and. abs(tdum2) .gt. abs(tdum1) )then
                    rspcore_psi_theta(1,i,k)=Rfasp(i,j)
                    zspcore_psi_theta(1,i,k)=Zfasp(i,j)
                    sspcore_psi_theta(1,i,k)=real(j-1)/real(mpara0tmp(i))
                    exit
                elseif(tdum2*tdum1 .le. 0 .and. abs(tdum2) .le. abs(tdum1))then
                    rspcore_psi_theta(1,i,k)=Rfasp(i,j+1)
                    zspcore_psi_theta(1,i,k)=Zfasp(i,j+1)
                    sspcore_psi_theta(1,i,k)=real(j)/real(mpara0tmp(i))
                    exit
                endif
  
            enddo
        enddo
        rspcore_psi_theta(1,i,lstcore)=Rfasp(i,1)
        zspcore_psi_theta(1,i,lstcore)=Zfasp(i,1)
        sspcore_psi_theta(1,i,lstcore)=1.0_lk
    enddo
  
    do i=1,lstcore
        if(rspcore_psi_theta(1,lspcore,i)<0.0_lk)then
            if(mype==0)write(gtcout,*)"1i",i
        endif
    enddo
  
    call construct_spline2d(lspcore,lstcore,del_psi_core,del_theta_core,rspcore_psi_theta,0,0)
    call construct_spline2d(lspcore,lstcore,del_psi_core,del_theta_core,zspcore_psi_theta,0,0)
    call construct_spline2d(lspcore,lstcore,del_psi_core,del_theta_core,sspcore_psi_theta,0,0)
  
    do i=1,lstcore
        if(rspcore_psi_theta(1,lspcore,i)<0.0_lk)then ! R.B.White's smooth algorithm for spline leads to the change of source value
           rspcore_psi_theta(1,lspcore,i)=0.0_lk
        endif
    enddo
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!       R(\psi,s),Z(\psi,s)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
    allocate(dscoretmp(lspcore),n_dscoretmp(lspcore))
  
    dscoretmp=0.0_lk
    n_dscoretmp=0.0_lk
  
    do i=1,lspcore
        aspcore_psi(1,i)=real(mpara0tmp(i))*deltas
    enddo
  
    do i=1,lspcore
      dscoretmp(i)=aspcore_psi(1,i)/real(lstcore-1)
      n_dscoretmp(i)=dscoretmp(i)/deltas
    enddo
   
    do i=1,lspcore
       do j=1,lstcore-1
        j1=nint(n_dscoretmp(i)*(j-1))+1
        rspcore_psi_s(1,i,j)=Rfasp(i,j1)
        zspcore_psi_s(1,i,j)=Zfasp(i,j1)
       enddo
        rspcore_psi_s(1,i,lstcore)=rspcore_psi_s(1,i,1)
        zspcore_psi_s(1,i,lstcore)=zspcore_psi_s(1,i,1)
    enddo
   
    del_s_core=1.0_lk/real(lstcore-1)
    call construct_spline1d(lspcore,del_psi_core,aspcore_psi,0)
    call construct_spline2d(lspcore,lstcore,del_psi_core,del_s_core,rspcore_psi_s,0,0)
    call construct_spline2d(lspcore,lstcore,del_psi_core,del_s_core,zspcore_psi_s,0,0)
  
    deallocate(Rfasp,Zfasp,tfasp,z0tmp,r0tmp,mpara0tmp,dscoretmp,n_dscoretmp)
  
end subroutine frc_coordinate_splinecore
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!         SOL region
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine frc_coordinate_splinesol
    use global_parameters,only:gtcout,mype,pi
    use simulation_domain, only: psi0_sol, psi1_sol
    use spline_function, only: spline1d,spline2d,sppsiR1d,construct_spline1d,construct_spline2d
    use equilibrium, only: r_Eq, z_Eq, psirz_Eq, &
                             psi_sep, lsr, lsz, spdR, spdZ, zmax
    use trace_fieldlines, only: trace, trace_sepsol
    use read_in_parameters, only: Zsol_domain0, Zsol_domain1
    
    integer :: i,j,ij,mtest,ngtmp,j1,k,ik
    real(lk) psifrctmp,deltas,tdum,tdum1,tdum2
    integer,dimension(:),allocatable :: mpara0tmp1,mpara0tmp2,mpara0tmp
    real(lk),dimension(:,:),allocatable :: Rfasp1,Zfasp1,tfasp1,&
                                           Rfasp2,Zfasp2,tfasp2,&
                                           Rfa,Zfa
    
    real(lk),dimension(:),allocatable :: z0tmp,r0tmp,num_deltassol
    
    msp_psi_sol(1)=psi0_sol
    msp_psi_sol(lspsol)=psi1_sol
    
    del_psi_sol=(msp_psi_sol(lspsol)-msp_psi_sol(1))/real(lspsol-1)
    del_s_sol=1.0_lk/real(lszsol-1)
    
    allocate(z0tmp(lspsol),r0tmp(lspsol))
    z0tmp=0.0_lk
    r0tmp=0.0_lk
    
    do i=1,lspsol
        msp_psi_sol(i)=msp_psi_sol(1)+del_psi_sol*real(i-1)
! notice iflag=1 for outmidplane rsol root in sol region(exclude core region on the separatrix)
        call sppsiR1d(1,lsr,lsz,spdR,spdZ,r_Eq,z_Eq,psirz_Eq,z0tmp(i)+zmax,msp_psi_sol(i),r0tmp(i),psifrctmp)
    enddo
    
    deltas=(Zsol_domain1-Zsol_domain0)/20000.0_lk
    ngtmp=int(3.0*abs(Zsol_domain1-Zsol_domain0)/deltas)
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                right part(positive Z)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    allocate(mpara0tmp1(lspsol))
    allocate(Rfasp1(lspsol,ngtmp),Zfasp1(lspsol,ngtmp),tfasp1(lspsol,ngtmp))
    Rfasp1=0.0_lk
    Zfasp1=0.0_lk
    tfasp1=0.0_lk
    mpara0tmp1=0
    
    do k=2,lspsol
        call trace(2,1,ngtmp,r0tmp(k),z0tmp(k),msp_psi_sol(k),deltas,Zsol_domain0,Zsol_domain1,mpara0tmp1(k),Rfasp1(k,:),Zfasp1(k,:),tfasp1(k,:))
    enddo
    if(psi0_sol .gt. psi_sep)then
        call trace(2,1,ngtmp,r0tmp(1),z0tmp(1),msp_psi_sol(1),deltas,Zsol_domain0,Zsol_domain1,mpara0tmp1(1),Rfasp1(1,:),Zfasp1(1,:),tfasp1(1,:))
    else
        call trace_sepsol(1,ngtmp,r0tmp(1),z0tmp(1),deltas,Zsol_domain0,Zsol_domain1,mpara0tmp1(1),Rfasp1(1,:),Zfasp1(1,:))
    endif
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!               Left part(negative Z)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    allocate(mpara0tmp2(lspsol))
    allocate(Rfasp2(lspsol,ngtmp),Zfasp2(lspsol,ngtmp),tfasp2(lspsol,ngtmp))
    
    Rfasp2=0.0_lk
    Zfasp2=0.0_lk
    tfasp2=0.0_lk
    mpara0tmp2=0
    
    do k=2,lspsol
        call trace(3,1,ngtmp,r0tmp(k),z0tmp(k),msp_psi_sol(k),deltas,Zsol_domain0,Zsol_domain1,mpara0tmp2(k),Rfasp2(k,:),Zfasp2(k,:),tfasp2(k,:))
    enddo
    if(psi0_sol .gt. psi_sep)then
        call trace(3,1,ngtmp,r0tmp(1),z0tmp(1),msp_psi_sol(1),deltas,Zsol_domain0,Zsol_domain1,mpara0tmp2(1),Rfasp2(1,:),Zfasp2(1,:),tfasp2(1,:))
    else
        call trace_sepsol(2,ngtmp,r0tmp(1),z0tmp(1),deltas,Zsol_domain0,Zsol_domain1,mpara0tmp2(1),Rfasp2(1,:),Zfasp2(1,:))
    endif
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    allocate(mpara0tmp(lspsol))
    mpara0tmp(:)=mpara0tmp1(:)+mpara0tmp2(:)-1
    ngtmp=maxval(mpara0tmp)
    
    allocate(Rfa(lspsol,ngtmp),Zfa(lspsol,ngtmp))
    Rfa=0.0_lk
    Zfa=0.0_lk
    
    do i=1,lspsol
        do j=1,mpara0tmp2(i)
            Rfa(i,j)=Rfasp2(i,mpara0tmp2(i)+1-j)
            Zfa(i,j)=Zfasp2(i,mpara0tmp2(i)+1-j)
        enddo
    enddo
    
    do i=1,lspsol
        do j=mpara0tmp2(i)+1,mpara0tmp(i)
            j1=j-mpara0tmp2(i)+1
            Rfa(i,j)=Rfasp1(i,j1)
            Zfa(i,j)=Zfasp1(i,j1)
       enddo
    enddo
     
    Zfa(:,1)=Zsol_domain0
    
    deallocate(Rfasp1,Zfasp1,tfasp1,mpara0tmp1,Rfasp2,Zfasp2,tfasp2,mpara0tmp2)
    
    del_z_sol=(Zsol_domain1-Zsol_domain0)/real(lszsol-1)
    
    do i=1,lspsol
        do k=1,lszsol-1
            do j=1,mpara0tmp(i)
    
                 tdum1=(Zfa(i,j)-real(k-1)*del_z_sol-Zsol_domain0)
                 tdum2=(Zfa(i,min(j+1,mpara0tmp(i)))-real(k-1)*del_z_sol-Zsol_domain0)
                 if(tdum2*tdum1 .le. 0 .and. abs(tdum2) .gt. abs(tdum1) )then
                   sspsol_psi_z(1,i,k)=real(j-1)/real(mpara0tmp(i))
                   exit
                 elseif(tdum2*tdum1 .le. 0 .and. abs(tdum2) .le. abs(tdum1))then
                   sspsol_psi_z(1,i,k)=real(j)/real(mpara0tmp(i))
                   exit
                 endif
    
            enddo
        enddo
        sspsol_psi_z(1,i,lszsol)=1.0_lk
    enddo
    
    call construct_spline2d(lspsol,lszsol,del_psi_sol,del_z_sol,sspsol_psi_z,0,0)
    
    do i=1,lspsol
        aspsol_psi(1,i)=real(mpara0tmp(i))*deltas
    enddo 
    
    
    allocate(num_deltassol(lspsol))
    num_deltassol=0.0_lk
    
    do i=1,lspsol
        num_deltassol(i)=real(mpara0tmp(i))/real(lszsol-1)
    enddo
     
    do i=1,lspsol
        do j=1,lszsol-1
            j1=nint(num_deltassol(i)*(j-1))+1
            rspsol_psi_s(1,i,j)=Rfa(i,min(j1,mpara0tmp(i)))
            zspsol_psi_s(1,i,j)=Zfa(i,min(j1,mpara0tmp(i)))
        enddo
        rspsol_psi_s(1,i,lszsol)=Rfa(i,mpara0tmp(i))
        zspsol_psi_s(1,i,lszsol)=Zsol_domain1
    enddo
    
    deallocate(num_deltassol)
    
    call construct_spline1d(lspsol,del_psi_sol,aspsol_psi,0)
    call construct_spline2d(lspsol,lszsol,del_psi_sol,del_s_sol,rspsol_psi_s,0,0)
    call construct_spline2d(lspsol,lszsol,del_psi_sol,del_s_sol,zspsol_psi_s,0,0)
   
    do i=1,lszsol
        if(rspsol_psi_s(1,1,i)<0.0_lk)then ! R.B.WHITE's smooth on spline leads the change of source value
           rspsol_psi_s(1,1,i)=0.0_lk
       endif
    enddo
    
end subroutine frc_coordinate_splinesol
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end module mapping_splines
