Here is a partial solution.
If $(X,Y)$ is jointly normal, with $\sigma^2_x \sigma_y^2>0$ (non-trivial case)
and $\mu_x=0$ then $XY$ cannot be normal.
If it were, then the third cumulant of $XY$
$$\kappa_3= 2\, \sigma_x\, \sigma_y(\sigma_y^2\, \sigma_x^2\, \rho^3\,
+3\, \sigma_y^2\, \sigma_x^2\, \rho
+3\, \sigma_x^2\, \rho\, \mu_y^2
+3\, \sigma_y\, \sigma_x\, \rho^2\, \mu_y\, \mu_x
+3\, \sigma_y\, \sigma_x\, \mu_x\, \mu_y
+3\, \sigma_y^2\rho\, \mu_x^2)$$
would vanish. Since $\mu_x=0$ this leads to
$$0=\rho(\rho^2 \sigma_y^2+3\sigma_y^2+3\mu_y^2).$$
Then $\rho=0$ and we are back in the independent case.
Update:
The more I think about it, the more I realize
that the result is not even plausible.
Generally speaking, the product $XY$ is much too likely to take values near zero
to be normal, or any other nice distribution.
Suppose that $(X,Y)$ has any "nice" bivariate distribution, with
a joint density function that is continuous and non-zero at the origin.
We argue that $XY$ cannot have a "nice" density at zero.
For $0<u<1$ define $$A_u=\{(x,y):u^{1/2}<x<u^{1/4},\, 0<y<u/x\}.$$
These sets shrink towards the origin as $u\to 0$.
We have
$$\mathbb{P}(0<XY\leq u)\geq \mathbb{P}((X,Y)\in A_u),$$
and so letting $\lambda$ denote two-dimensional Lebesgue measure
also
$${\mathbb{P}(0<XY\leq u)\over u}
\geq {\mathbb{P}((X,Y)\in A_u)\over\lambda(A_u)}\cdot{\lambda(A_u)\over u}
= {\mathbb{P}((X,Y)\in A_u)\over\lambda(A_u)}\cdot{\log(1/u)\over 4}.$$
As $u\to 0$, the first factor on the right converges to $f_{(X,Y)}(0,0)$ and the
second one blows up. Therefore, the left hand side also blows up, showing that $XY$
cannot have a continuous density with finite value at $u=0$.
Yes, if you know the joint probability distribution of $X,Y$ then the distributions of $U$ and $V$ are, respectively:
$$\begin{align}
\mathsf P(U=u) & = \sum_\overbrace{x\in\Omega_X, u-x\in\Omega_Y} \mathsf P(X=x, Y=u-x)
\\
\mathsf P(V=v) & = \sum_\overbrace{x\in\Omega_X, v/x\in\Omega_Y} \mathsf P(X=x, Y=v/x)
\end{align}$$
The joint distribution on $U,V$, would similarly be $$\mathsf P(U=u, V=v) = \sum_\overbrace{x\in\Omega_X, u-x\in\Omega_Y, v/x\in\Omega_Y} \mathsf P(X=x, Y=u-x, Y=v/x)$$
However, for any given pair of $U=u,V=v$ there is in fact only two corresponding pair of $X=\Box,Y=\Box$.
$\begin{align}
(U=X+Y) \wedge (V= XY) & \iff
\left(X= \dfrac{U\pm\sqrt{U^2-4V}}{2}\right)
\wedge \left(Y=\dfrac{U\mp\sqrt{U^2-4V}}{2}\right)
\\[4ex]
\therefore \mathsf P(U=u, V=v)
& = \mathsf P\left(X= \dfrac{u\pm\sqrt{u^2-4v}}{2}, Y=\dfrac{u\mp\sqrt{u^2-4v}}{2}\right)
\\[1ex]
& = {\mathsf P\left(X= \tfrac{u+\sqrt{u^2-4v}}{2}, Y=\tfrac{u-\sqrt{u^2-4v}}{2}\right)
+ \mathsf P\left(X= \tfrac{u-\sqrt{u^2-4v}}{2}, Y=\tfrac{u+\sqrt{u^2-4v}}{2}\right)}
\end{align}$
Best Answer
Elaborating on alex.jordan's answer. Consider the one-dimensional case, where $X$ and $Z$ are independent ${\rm N}(0,1)$ variables. Then, $$ p_{XZ,Z} (xz,z) = p_{XZ|Z} (xz|z)p_Z (z). $$ Now, conditionally on $Z=z$, $XZ$ is normally distributed with mean $0$ and variance ${\rm Var}(Xz)=z^2$. Hence, $$ p_{XZ|Z} (xz|z) = \frac{1}{{\sqrt {2\pi z^2 } }}\exp \bigg[ - \frac{{(xz)^2 }}{{2z^2 }}\bigg] = \frac{1}{{\sqrt {2\pi z^2 } }}e^{ - x^2 /2} , $$ and in turn, $$ p_{XZ,Z} (xz,z) = \frac{1}{{\sqrt {2\pi z^2 } }}e^{ - x^2 /2} \frac{1}{{\sqrt {2\pi } }}e^{ - z^2 /2} = \frac{1}{{2\pi |z|}}e^{ - (x^2 + z^2 )/2} , $$ which is not Gaussian because of the $|z|$ in the denominator.