Convolution of normal probability distribution and Gaussian filter

I'm interested in generating a random scalar field according to a probability density function with a correlation length. This might correspond to say a spatial distribution of strength of a material or a temporal distribution of a signal.

My plan for doing so is to generate normally-distributed random numbers $R$ at each sampling point $x_i$ and introduce the spatial correlation by convolving with a Gaussian filter. Specifically, I perform this in the frequency domain by multiplying the Fourier transforms element-wise and then performing the inverse transform to reconstruct the filtered variable. I'm using numpy's multi-dimensional Fourier transforms to do the calculation (numpy.fft.fftn, numpy.fft.ifftn).

$$R(x_i) \sim N(\mu=0,\sigma^2=1)$$

$$f(x_i) = \exp \left( \frac{-x_i^2}{\ell^2/2}\right)$$

$$\tilde{R} = \alpha \cdot IFFT_n (FFT_n(R) \cdots FFT_n(f))$$

Where I am not sure is about how to scale the filtered random field $\tilde{R}$ such that it is still normally distributed ($\sigma^2(\tilde{R})=1$). The mean is zero, but the variance is generally not one and nonlinearly varies with the correlation length $\ell$. What scaling factor $\alpha$ should I use, and how is this derived? I'm interested primarily in three-dimensional fields ($n=3$, $x_i \in R^3$).