Note that, with $y\neq 1/2$ fixed, $f(x,y) = 2y \neq 1$ for irrational $x$ and $f(x,y) = 1$ for rational $x$. Since the rationals are dense, the map $x \mapsto f(x,y)$ is everywhere discontinuous. Hence, $\int_0^1 f(x,y) \, dx$ fails to exist for almost every $y$ and the second iterated (Riemann) integral fails to exist.
Nonexistence of one iterated integral is not sufficient to prove that $\int_{[0,1]^2} f(x,y) \,d(x,y)$ fails to exist. A counterexample can be furnished upon request.
However, we can prove that $f$ is not Riemann integrable over $[0,1]^2$ directly by showing the infimum of upper Darboux sums (upper integral) does not equal the supremum of lower Darboux sums (lower integral).
Take any partition $P$ of $[0,1]^2$ with $m\cdot n$ subrectangles $R_{jk} =[x_{j-1},x_j]\times [y_{k-1},y_k]$ formed from points $0 = x_0 < x_1 < \ldots < x_m = 1$ and $0 = y_0 < y_1 < \ldots < y_n= 1.$
Let $q$ be the index such that $y_{q-1} < 1/2 \leqslant y_q$. Since each subrectangle contains points $(x,y)$ where $x$ is rational and where $x$ is irrational, we have, for all $1 \leqslant j \leqslant m$,
$$M_{jk} =\sup_{x \in R_{jk}}f(x,y) = \begin{cases}2y_k,& k \geqslant q \\1, & k < q \end{cases}\quad ,\,\,m_{jk} =\inf_{x \in R_{jk}}f(x,y) = \begin{cases}1,& k \geqslant q \\2y_k, & k < q \end{cases} $$
Noting that $M_{jk}$ is independent of $j$ and $\sum_{j=1}^m (x_j - x_{j-1}) = 1$, the upper sum is
$$\begin{align}U(P,f) &= \sum_{j=1}^m\sum_{k=1}^n M_{jk}(x_j- x_{j-1})(y_k - y_{k-1}) \\ &=\sum_{k=1}^{q-1} (1)(y_k - y_{k-1}) + \sum_{k=q}^{n}2y_k(y_k - y_{k-1})\\ &\geqslant y_{q-1} + \int_{y_{q-1}}^1 2y \, dy \\ &= y_{q-1} + 1 - y_{q-1}^2 \end{align} $$
We have $y_{q-1} \to 1/2$ as the partition norm $\|P \| \to 0$ and the upper Darboux integral must satisfy
$$\overline{\int}_{[0,1]^2} f(x,y) \, dy = \lim_{\|P\|\to 0}U(P,f) \geqslant 5/4$$
Similarly, we can show for lower sums that
$$L(P,f) \leqslant \int_{0}^{y_{q-1}} 2y \, dy + 1 - y_{q} = y_{q-1}^2 + 1 - y_{q-1},$$
and
$$ \underline{\int}_{[0,1]^2} f(x,y) \, d(x,y) = \lim_{\|P\|\to 0}L(P,f) \leqslant 3/4$$
Therefore, $f$ is not Riemann integrable over $[0,1]^2$ since the upper and lower Darboux integrals are not equal.
Best Answer
We are seeing the step discontinuity.
Let $x_{0}=\frac{1}{n}=\frac{1}{3}$. As we approach from the left, $$\color{blue}{2.9,\ 2.99,\ 2.999},\dots$$ The function value is $$\frac{1}{0!} + \frac{1}{1!} + \frac{1}{2!}$$ As we approach from the right, $$\color{red}{ 3.1,\ 3.01,\ 3.001}$$ The function value is $$\frac{1}{0!} + \frac{1}{1!} + \frac{1}{2!} + \boxed{\frac{1}{3!}}$$
The difference is $\frac{1}{n!} = \frac{1}{3!}$