read_tov_sph_id Subroutine

  SUBROUTINE read_tov_sph_id(filename1, filename2)

    !***********************************************************
    !
    !# Read the two SPH TOV ID files produced with setup_TOV.x,
    !  and place them symmetrically on the \(x\) axis so that
    !  their distance is equal to the periastron given as input
    !
    !  FT 13.12.2022
    !
    !***********************************************************

    USE sph_variables,  ONLY: npart,  &  ! particle number
                              n1, n2, &  ! particle numbers for each star
                              pos_u,  &  ! particle positions
                              vel_u,  &  ! particle velocities in
                                         ! coordinate frame
                              nlrf,   &  ! baryon number density in
                                         ! local rest frame
                              nu,     &  ! canonical baryon number per
                                         ! particle
                              Theta,  &  ! Generalized Lorentz factor
                              h,      &  ! Smoothing length
                              Pr,     &  ! Pressure
                              u,      &  ! Internal energy in local rest
                                         ! frame (no kinetic energy)
                              Ye,     &  ! Electron fraction
                              allocate_sph_memory, deallocate_sph_memory
    USE input_output,   ONLY: set_units, read_sphincs_dump
    USE utility,        ONLY: scan_1d_array_for_nans

    IMPLICIT NONE

    CHARACTER(LEN=*), INTENT(INOUT):: filename1, filename2

    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: pos_u1, vel_u1, &
                                                    pos_u2, vel_u2
    DOUBLE PRECISION, DIMENSION(:),   ALLOCATABLE:: u1, nu1, h1, nlrf1, &
                                                    Pr1, Ye1, Theta1, &
                                                    u2, nu2, h2, nlrf2, &
                                                    Pr2, Ye2, Theta2

    CALL set_units('NSM')

    !
    !-- Read just the particle number, to be able to allocate needed memory
    !
    OPEN(10, file= filename1, form='UNFORMATTED')
    READ(10) npart
    CLOSE(10)

    CALL allocate_sph_memory()

    CALL read_sphincs_dump(filename1)

    !PRINT *, "npart=", npart

    ALLOCATE(pos_u1(3,npart))
    ALLOCATE(vel_u1(3,npart))
    ALLOCATE(u1    (npart))
    ALLOCATE(nu1   (npart))
    ALLOCATE(h1    (npart))
    ALLOCATE(nlrf1 (npart))
    ALLOCATE(Pr1   (npart))
    ALLOCATE(Ye1   (npart))
    ALLOCATE(Theta1(npart))

    n1    = npart
    pos_u1= pos_u(:,1:npart)
    vel_u1= vel_u(:,1:npart)
    u1    = u    (1:npart)
    nu1   = nu   (1:npart)
    h1    = h    (1:npart)
    nlrf1 = nlrf (1:npart)
    Pr1   = Pr   (1:npart)
    Ye1   = Ye   (1:npart)
    Theta1= Theta(1:npart)

    !PRINT *, "SIZE(nlrf)=", SIZE(nlrf1)

    CALL deallocate_sph_memory()

    !
    !-- Read just the particle number, to be able to allocate needed memory
    !
    OPEN(10, file= filename2, form='UNFORMATTED')
    READ(10) npart
    CLOSE(10)

    CALL allocate_sph_memory()

    CALL read_sphincs_dump(filename2)

    ALLOCATE(pos_u2(3,npart))
    ALLOCATE(vel_u2(3,npart))
    ALLOCATE(u2    (npart))
    ALLOCATE(nu2   (npart))
    ALLOCATE(h2    (npart))
    ALLOCATE(nlrf2 (npart))
    ALLOCATE(Pr2   (npart))
    ALLOCATE(Ye2   (npart))
    ALLOCATE(Theta2(npart))

    n2    = npart
    pos_u2= pos_u(:,1:npart)
    vel_u2= vel_u(:,1:npart)
    u2    = u    (1:npart)
    nu2   = nu   (1:npart)
    h2    = h    (1:npart)
    nlrf2 = nlrf (1:npart)
    Pr2   = Pr   (1:npart)
    Ye2   = Ye   (1:npart)
    Theta2= Theta(1:npart)

    CALL deallocate_sph_memory()

    npart= n1 + n2

    CALL allocate_sph_memory()

    pos_u(1,1:n1)  = pos_u1(1,:)
    pos_u(2:3,1:n1)= pos_u1(2:3,:)
    vel_u(:,1:n1)  = vel_u1
    u    (1:n1)    = u1
    nu   (1:n1)    = nu1
    h    (1:n1)    = h1
    nlrf (1:n1)    = nlrf1
    Pr   (1:n1)    = Pr1
    Ye   (1:n1)    = Ye1
    Theta(1:n1)    = Theta1

    pos_u(1,n1 + 1: npart)  = pos_u2(1,:)
    pos_u(2:3,n1 + 1: npart)= pos_u2(2:3,:)
    vel_u(:,n1 + 1: npart)  = vel_u2
    u    (n1 + 1: npart)    = u2
    nu   (n1 + 1: npart)    = nu2
    h    (n1 + 1: npart)    = h2
    nlrf (n1 + 1: npart)    = nlrf2
    Pr   (n1 + 1: npart)    = Pr2
    Ye   (n1 + 1: npart)    = Ye2
    Theta(n1 + 1: npart)    = Theta2

    DEALLOCATE(pos_u1)
    DEALLOCATE(vel_u1)
    DEALLOCATE(u1)
    DEALLOCATE(nu1)
    DEALLOCATE(h1)
    DEALLOCATE(nlrf1)
    DEALLOCATE(Pr1)
    DEALLOCATE(Ye1)
    DEALLOCATE(Theta1)
    DEALLOCATE(pos_u2)
    DEALLOCATE(vel_u2)
    DEALLOCATE(u2)
    DEALLOCATE(nu2)
    DEALLOCATE(h2)
    DEALLOCATE(nlrf2)
    DEALLOCATE(Pr2)
    DEALLOCATE(Ye2)
    DEALLOCATE(Theta2)

    !
    !-- Ensure that the ID does not contain NaNs or infinities
    !
    PRINT *, "** Ensuring that the SPH ID does not have any NaNs or", &
             "infinities..."

    CALL scan_1d_array_for_nans( npart, pos_u(1,:), "pos_u(:,1)" )
    CALL scan_1d_array_for_nans( npart, pos_u(2,:), "pos_u(:,2)" )
    CALL scan_1d_array_for_nans( npart, pos_u(3,:), "pos_u(:,3)" )

    CALL scan_1d_array_for_nans( npart, nlrf, "nlrf" )
    CALL scan_1d_array_for_nans( npart, nu, "nu" )
    CALL scan_1d_array_for_nans( npart, u, "u" )
    CALL scan_1d_array_for_nans( npart, h, "h" )
    CALL scan_1d_array_for_nans( npart, Pr, "Pr" )
    CALL scan_1d_array_for_nans( npart, Ye, "Ye" )
    CALL scan_1d_array_for_nans( npart, Theta, "Theta" )

    CALL scan_1d_array_for_nans( npart, vel_u(1,:), "vel_u(:,1)" )
    CALL scan_1d_array_for_nans( npart, vel_u(2,:), "vel_u(:,2)" )
    CALL scan_1d_array_for_nans( npart, vel_u(3,:), "vel_u(:,3)" )

    PRINT *, " * the SPH ID does not have NaNs or infinities."
    PRINT *

  END SUBROUTINE read_tov_sph_id