There are a number of approaches to simulating a set of dependent random variables. In the context of spatial random fields, these include sequential indicator methods, turning bands, and the Karhunen-Loeve expansion. See Christakos (1992, Chapter 8) and Deutsch and Journel (1992, Chapter 5) for details.
In addition, there is the LU decomposition method, a particularly simple and computationally efficient for normal or Gaussian variates. For a given covariance matrix, the decomposition is computed once, and the simulation proceeds by repeatedly generating a vector of independent random variables and multiplying by the matrix.
One problem with this technique is that memory is required to hold the covariance matrix of all the analysis and conditioning variables in core.