pure subroutine find_intersec(nelem_in_1, nelem_in_2, arr_1, arr_2, intersec, nelem_out)
integer(int32), intent(in) :: nelem_in_1, nelem_in_2
integer(int32), intent(in) :: arr_1(1:), arr_2(1:)
integer(int32), intent(inout) :: intersec(1:)
integer(int32), intent(out) :: nelem_out
integer :: i, j
nelem_out = 0
i = 1
j = 1
do while (i <= nelem_in_1 .and. j <= nelem_in_2)
if (arr_1(i) < arr_2(j)) then
i = i + 1
else if (arr_1(i) > arr_2(j)) then
j = j + 1
else
nelem_out = nelem_out + 1
!intersec(nelem_out) = arr_1(i)
intersec(nelem_out) = i
i = i + 1
j = j + 1
end if
end do
end subroutine find_intersec