You can construct   examples a dime a dozen. $\newcommand{\bR}{\mathbb{R}}$

 Here  is a first  simple way. Fix $N$ smooth functions $f_1,f_2,\dotsc, f_N:M\to\bR$ and $N$ independent Gaussian random variables. Then

$$f(x)=\sum_{k=1}^N X_k f_k(x) $$

is a Gaussian random field and the sample functions are a.s. smooth.

You can   allow  infinitely many functions  in the above examples, but then  you need to make some assumptions on the functions and the variables $X$ to guarantee  the convergence of the resulting series. 

We know from Kolmogorov that the convergence is a $0-1$ phenomenon. The the three-series theorem  tells us what these conditions should be.

Here is an example of this kind when $M$ is compact. Fix a  Riemann metric, denote by $\Delta$ the resulting Laplacian. Its eigenvalues (multiplicities included) are

$$  0=\lambda_0<\lambda_1\leq \lambda_2\leq \cdots. $$

Fix an orthonormal eigenbasis of $L^2(M, dV_g)$ $(\psi_k)_{k\geq 0}$

$$\Delta \psi_k=\lambda_k,\;\; \Vert\psi_k\Vert_{L^2}=1. $$

Next choose  independent  Gaussian random variables $(X_k)_{k\geq 0}$. We denote by $v_k$ the variance of $X_k$.   If $v_k$ goes to $0$ sufficiently fast, then the random series 

$$f(x)=\sum_{k\geq 0} X_k \psi_k(x) $$

defines a  Gaussian random field on $M$.  The regularity of the sample functions  of this random field   depends on the decay rate of  $v_k$. The faster $v_k$ decays as $k\to\infty$, the more regular is the random function. For example, if 

$$\lim_{k\to \infty}k^\alpha v_k=0,\;\;\forall \alpha>0, $$

then the random function $f(x)$ is a.s. smooth.

Ultimately, the most general construction  of a Gaussian random function on a manifold  is via Gaussian measures on the space of   distributions (i.e. generalized functions) on $M$.  I refer to   [*this paper*][1]  for  more details and additional references. I would start with reference [6] in this paper.


  [1]: http://www3.nd.edu/~lnicolae/GB2.pdf