My general question is how to construct an isotropic random vector field $\vec f: \mathbb{R}^3 \to \mathbb{R}^3$ with a given mean magnitude $\mathbb{E}[\|\vec f(\vec x)\|]=\mu$ and with vector magnitudes and direction correlated up to some length scales $l$ (beyond which the correlation goes to zero). I prefer a construction satisfying $\nabla \cdot \vec f=0$, but a construction without that condition is already interesting.
More precisely, we are given a vector $\vec{\mu}$ in $\mathbb{R^3}$ and a matrix-valued correlation function $C: \mathbb{R^3} \times \mathbb{R^3} \to M_3(\mathbb{R})$ which is isotropic, i.e $C(\vec v,\vec w)=C(\|\vec v-\vec w\|)$. One may define a Gaussian Process $f(\vec x) \sim GP(\vec{\mu},C)$ such that $\mathbb{E}[\vec f(\vec x)]=\vec\mu$ and $Cov(\vec f(v),\vec f(w))=C(\vec v,\vec w)$. I believe that constructions of such a random field $f$ are known. The equivalent problem for a scalar field $f$ seems well-studied: one can first draw the field in Fourier space by normalizing a white noise field with the appropriate power spectrum $P(k)$, and then transform back to real space. Here I am looking for a simple procedure for $\mathbb{R}^3$; I think this is known but I haven't found a clear description anywhere.
Question: How can one generate such a random field $f$ with a mean vector of zero ($\vec{\mu}=\vec 0$), and a mean magnitude of $\mu$ ($\mathbb{E}[\|\vec f(\vec x)\|]=\mu$)? And what if we also impose $\nabla \cdot \vec f=0$?
