module simulation_domain
use precision, only: lk
implicit none

public  frc_domain_setup, &
        Zsol_domain_width, psi0_sol, psi1_sol, psi0_core, psi1_core, &
        R0sol, R1sol, R0core, R1core, R_sep, &
        sfcore, sfsol, couplefrc, &
        Z_corner_left, Z_corner_right, Z_corner_width
private



integer  sfcore, sfsol, couplefrc
real(lk) Zsol_domain_width, &
         Z_corner_left, Z_corner_right, Z_corner_width, & 
         R0core, R1core, psi0_core, psi1_core, &
         R0sol, R1sol, psi0_sol, psi1_sol, R_sep
save


contains

!========================================================================
subroutine frc_domain_setup
    use global_parameters,only:gtcout,mype
    use read_in_parameters, only: numerfrc, Zsol_domain0, Zsol_domain1, Rdomain0, Rdomain1
    use spline_function, only: spline2d,sppsiR1d
    use equilibrium, only: psirz_Eq, ne_Eq, r_Eq, z_Eq,&
                           psi_sep, lsr, lsz, spdR, spdZ, zmax, rmax
    use trace_fieldlines, only: trace_sep
    implicit none

    integer :: i,ngtmp,mpara0tmp
    real(lk) :: psifrctmp,deltas
    real(lk),dimension(:),allocatable :: Rfasp,Zfasp,tfasp

    deltas=0.001_lk
    ngtmp=int(3.0*zmax/deltas)
    allocate(Rfasp(ngtmp),Zfasp(ngtmp),tfasp(ngtmp))
    Rfasp=0.0_lk
    Zfasp=0.0_lk
    tfasp=0.0_lk
    call sppsiR1d(1,lsr,lsz,spdR,spdZ,r_Eq,z_Eq,psirz_Eq,zmax,psi_sep,R_sep,psifrctmp)
! Notice that first trace field line towards right, and then towards left!    
! get the values of Z_corner_right, Z_corner_left, Z_corner_width
    call trace_sep(ngtmp,R_sep,0.0_lk,deltas,mpara0tmp,Rfasp,Zfasp,tfasp,&
                   Z_corner_right,Z_corner_left,Z_corner_width)
    deallocate(Rfasp,Zfasp,tfasp)

    if(mype==0)then
        write(gtcout,*)"R_sep=",R_sep
        write(gtcout,*)"Z_corner_left=",Z_corner_left
        write(gtcout,*)"Z_corner_right=",Z_corner_right
        write(gtcout,*)"Z_corner_width=",Z_corner_width
    endif

! safety guard for Rdomain0 and Rdomain1 values
    if(Rdomain0<0.0_lk .or. Rdomain1<0.0_lk .or. Rdomain0>rmax .or. Rdomain1>rmax)then
        if(mype==0)write(gtcout,*)"Rdomain0 and Rdomain1 are out of equilibrium range"
        stop
    endif

    if(numerfrc==1)then           ! SOL and CORE coupled simulation
        if(.not. (Rdomain0<R_sep .and. Rdomain1>R_sep))then
            if(mype==0)write(gtcout,*)"wrong Rdomain0 and Rdomain1 values for numerfrc=1"
            stop
        endif
        sfcore=1
        sfsol=1
        couplefrc=1

        R0core=Rdomain0
        R1core=R_sep
        R0sol=R_sep
        R1sol=Rdomain1
        
        psi0_core=spline2d(R0core,zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        psi1_core=psi_sep
        psi0_sol=psi_sep
        psi1_sol=spline2d(R1sol,zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

        if(mype==0)then
            write(gtcout,*)'Global coupled simulation'
            write(gtcout,*)"sfcore,sfsol,couplefrc=",sfcore,sfsol,couplefrc
            write(gtcout,*)"R0core,R1core=",R0core,R1core
            write(gtcout,*)"psi0_core,psi1_core=",psi0_core,psi1_core
            write(gtcout,*)"R0sol,R1sol=",R0sol,R1sol
            write(gtcout,*)"psi0_sol,psi1_sol=",psi0_sol,psi1_sol
        endif

    elseif(numerfrc==2)then    ! Independent core simulation
        if(.not. (Rdomain0<R_sep .and. Rdomain1<R_sep .and. Rdomain0<Rdomain1))then
            if(mype==0)write(gtcout,*)"wrong Rdomain0 and Rdomain1 values for numerfrc=2"
            stop
        endif
        sfcore=1
        sfsol=0
        couplefrc=0

        R0core=Rdomain0
        R1core=Rdomain1
        psi0_core=spline2d(R0core,zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        psi1_core=spline2d(R1core,zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

        if(mype==0)then
            write(gtcout,*)'Core simulation'
            write(gtcout,*)"sfcore,sfsol,couplefrc=",sfcore,sfsol,couplefrc
            write(gtcout,*)"R0core,R1core=",R0core,R1core
            write(gtcout,*)"psi0_core,psi1_core=",psi0_core,psi1_core
        endif

    elseif(numerfrc==3)then    ! Independent SOL simulation
        if(.not.(Rdomain0>R_sep .and. Rdomain1>R_sep .and. Rdomain0<Rdomain1))then
            if(mype==0)write(gtcout,*)"wrong Rdomain0 and Rdomain1 values for numerfrc=3"
            if(mype==0)write(gtcout,*)"R0",Rdomain0,"R1",Rdomain1,"R_sep",R_sep
            stop
        endif
        sfcore=0
        sfsol=1
        couplefrc=0

        R0sol=Rdomain0
        R1sol=Rdomain1
        psi0_sol=spline2d(R0sol,zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)
        psi1_sol=spline2d(R1sol,zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

        if(mype==0)then
            write(gtcout,*)'SOL simulation'
            write(gtcout,*)"sfcore,sfsol,couplefrc=",sfcore,sfsol,couplefrc
            write(gtcout,*)"R0sol,R1sol=",R0sol,R1sol
            write(gtcout,*)"psi0_sol,psi1_sol=",psi0_sol,psi1_sol
        endif

    elseif(numerfrc==4)then    ! Independent SOL simulation with separatrix
        sfcore=0
        sfsol=1
        couplefrc=0

        R0sol=R_sep
        R1sol=Rdomain1
        psi0_sol = psi_sep
        psi1_sol = spline2d(R1sol,zmax,0,lsr,lsz,spdR,spdZ,psirz_Eq,0,0)

        if(mype==0)then
            write(gtcout,*)'SOL simulation starts at separatrix'
            write(gtcout,*)"sfcore,sfsol,couplefrc=",sfcore,sfsol,couplefrc
            write(gtcout,*)"R0sol,R1sol=",R0sol,R1sol
            write(gtcout,*)"psi0_sol,psi1_sol=",psi0_sol,psi1_sol
        endif
    endif

    Zsol_domain_width=Zsol_domain1-Zsol_domain0
    if(mype==0)write(gtcout,*)"Zsol_domain_width=",Zsol_domain_width

end subroutine frc_domain_setup

end module simulation_domain
