check_if_in_profiles_domain Subroutine

  subroutine check_if_in_profiles_domain(Y,flag)
    !! @note Subrotuine that checks if particles in the simulation are
    !! within the spatial domain where interpolants and plasma profiles
    !! are known. @endnote
    !!External plasma profiles and interpolants can have different spatial
    !! domains where they are defined. Therefore, it is necessary to check
    !! if a given particle has left these spatial domains to stop following
    !! it, otherwise this will cause an error in the simulation.
    REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN)      :: Y
    !! Particles' position in cylindrical coordinates,
    !! Y(1,:) = \(R\), Y(2,:) = \(\phi\), and Y(3,:) = \(Z\).
    INTEGER(is), DIMENSION(:), ALLOCATABLE, INTENT(INOUT)  :: flag
    !! Flag that determines whether particles are followed
    !! in the simulation (flag=1), or not (flag=0).
    INTEGER                                                :: IR
    !! @param IR Variable used to localize the grid cell in
    !! the \((R,\phi,Z)\) or \((R,Z)\) grid containing the fields data that
    !! corresponds to the radial position of the particles.
    INTEGER                                                :: IPHI
    !! @param IPHI Variable used to localize the grid cell in
    !! the \((R,\phi,Z)\) or \((R,Z)\) grid containing the fields data that
    !! corresponds to the azimuthal position of the particles.
    INTEGER                                                :: IZ
    !! @param IZ Variable used to localize the grid cell in the
    !! \((R,\phi,Z)\) or \((R,Z)\) grid containing the fields data that
    !! corresponds to the vertical position of the particles.
    INTEGER(ip)                                            :: pp
    !! @param pp Particle iterator.
    INTEGER(ip)                                            :: ss
    !! @param ss Species iterator.
    
    if (size(Y,1).eq.1) then
       ss = size(Y,1)
    else
       if (Y(2,1).eq.0) then
          ss=1_idef
       else
          ss = size(Y,1)
       end if
    endif

    if (ALLOCATED(profiles_domain%FLAG3D)) then
       !$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,IR,IPHI,IZ) &
       !$OMP& SHARED(Y,flag,profiles_domain,profiles_3d)
       do pp=1_idef,ss
          IR = INT(FLOOR((Y(pp,1)  - profiles_domain%Ro + &
               0.5_rp*profiles_domain%DR)/profiles_domain%DR) + 1.0_rp,idef)
          IPHI = INT(FLOOR((Y(pp,2)  + 0.5_rp*profiles_domain%DPHI)/ &
               profiles_domain%DPHI) + 1.0_rp,idef)
          IZ = INT(FLOOR((Y(pp,3)  + ABS(profiles_domain%Zo) + &
               0.5_rp*profiles_domain%DZ)/profiles_domain%DZ) + 1.0_rp,idef)

          if ((profiles_domain%FLAG3D(IR,IPHI,IZ).NE.1_is).OR. &
               ((IR.GT.profiles_3d%NR).OR.(IZ.GT.profiles_3d%NZ))) then
             flag(pp) = 0_is
          end if
       end do
       !$OMP END PARALLEL DO
    else
       !$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,IR,IZ) &
       !$OMP& SHARED(Y,flag,profiles_domain,profiles_2d)
       do pp=1_idef,ss
          IR = INT(FLOOR((Y(pp,1)  - profiles_domain%Ro + &
               0.5_rp*profiles_domain%DR)/profiles_domain%DR) + 1.0_rp,idef)
          IZ = INT(FLOOR((Y(pp,3)  + ABS(profiles_domain%Zo) + &
               0.5_rp*profiles_domain%DZ)/profiles_domain%DZ) + 1.0_rp,idef)

          if ((profiles_domain%FLAG2D(IR,IZ).NE.1_is).OR. &
               ((IR.GT.profiles_2d%NR).OR.(IZ.GT.profiles_2d%NZ))) then
             flag(pp) = 0_is
          end if
       end do
       !$OMP END PARALLEL DO
    end if
  end subroutine check_if_in_profiles_domain