! read snap, then import to pyton
module snap_parameters
    integer numerfrc,frcsnap, mpsi_core,mpara_core1,mgrid_core,mpsi_sol,&
    mpara_sol1, mgrid_sol,mtoroidal,mgrid

    !TODO: confirm the size of phi
    real(kind=8),dimension(:),ALLOCATABLE:: phi,phi_real,phi_imag,densityi,&
    flowi, densitye, flowe, pipara, piperp, pepara, peperp, phi_r_zeta_core,&
    phi_r_zeta_sol
    CHARACTER(len=999) fname
!f2py intent(in) fname

    contains
    subroutine filename(name)
        CHARACTER(*) name
        fname=trim(name)
    end subroutine filename
    subroutine deallocate
        implicit none
        if(numerfrc==1)then
            deallocate(phi, phi_real, phi_imag, densityi, &
                densitye, flowi, pipara, piperp, &
                flowe, pepara, peperp, phi_r_zeta_core,&
                phi_r_zeta_sol)
        elseif(numerfrc==3)then
            deallocate(phi, phi_real, phi_imag, densityi, &
                densitye, flowi, pipara, piperp, &
                flowe, pepara, peperp, &
                phi_r_zeta_sol)
        endif
    end subroutine deallocate
    subroutine snap_read
        frcsnap=321
        open(frcsnap,file=fname,status='old')
        if(numerfrc==1)then
           read(frcsnap,101) mpsi_core, mpara_core1, mgrid_core
           read(frcsnap,101) mpsi_sol, mpara_sol1, mgrid_sol
           read(frcsnap,101) mtoroidal
           mgrid=mgrid_core+mgrid_sol
           allocate(phi(mgrid), phi_real(mgrid), phi_imag(mgrid), densityi(mgrid), &
                densitye(mgrid), flowi(mgrid), pipara(mgrid), piperp(mgrid), &
                flowe(mgrid), pepara(mgrid), peperp(mgrid), phi_r_zeta_core(mpsi_core*mtoroidal),&
                phi_r_zeta_sol(mpsi_sol*mtoroidal))
           read(frcsnap,102) phi, phi_real, phi_imag, phi_r_zeta_core, phi_r_zeta_sol
           read(frcsnap,102) densityi, flowi, pipara, peperp
           read(frcsnap,102) densitye, flowe, pepara, peperp
        elseif(numerfrc==3)then
            read(frcsnap,101)mpsi_sol,mpara_sol1,mgrid_sol
            read(frcsnap,101)mtoroidal
            mgrid=mgrid_sol
            allocate(phi(mgrid), phi_real(mgrid), phi_imag(mgrid), densityi(mgrid), &
                densitye(mgrid), flowi(mgrid), pipara(mgrid), piperp(mgrid), &
                flowe(mgrid), pepara(mgrid), peperp(mgrid), &
                phi_r_zeta_sol(mpsi_sol*mtoroidal))
            read(frcsnap,102)phi,phi_real,phi_imag,phi_r_zeta_sol
            read(frcsnap,102)densityi,flowi,pipara,piperp
            read(frcsnap,102)densitye,flowe,pepara,peperp
        else
           write(*,*) 'Warning: in snap_read,numerfrc/=1 unsupported, exiting...'
           stop
        endif
        close(frcsnap)

    101 format(i6)
    102 format(e13.6)
    end subroutine snap_read
end module
