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$.
Let $X$ and $Y$ be independent random variables each with normal distribution (the means and variances need not be the same). Let $W=XY$. Then
$$\text{Var}(W)=E((XY)^2)-(E(XY))^2.$$
We need to compute the two expectations on the right.
By independence, $E(XY)=E(X)E(Y)$. Also, $E(X^2Y^2)=E(X^2)E(Y^2)$.
Since $E(X^2)=\text{Var}(X) +(E(X))^2$, with a similar expression for $E(Y^2)$, once we know the mean and variance of $X$ and $Y$, we can use the above equations to find $\text{Var}(XY)$.
Added: Let the means be $\mu$ and $\nu$, and the variances be $\sigma^2$ and $\tau^2$. Then the variance of $XY$ is, by the above argument, equal to
$$(\sigma^2+\mu^2)(\tau^2+\nu^2)-\mu^2\nu^2.$$
This simplifies to $\sigma^2\tau^2+\sigma^2\nu^2+\tau^2\mu^2$.
Note that we did not use the normality of $X$ and $Y$.
Best Answer
Suppose $X = N(\mu_1,\sigma^2_1), Y = N(\mu_2, \sigma^2_2)$ independent . $E[Z] = E[X]E[Y] = \mu_1\mu_2$. $E[Z^2] = E[X^2]E[Y^2] = (\sigma_1^2+\mu_1^2)(\sigma_2^2+\mu_2^2)$
$Var(Z) = (\sigma_1^2+\mu_1^2)(\sigma_2^2+\mu_2^2) - \mu_1^2\mu_2^2$