The existing model for multivariate skew normal data does not cohere with the joint distribution of a random sample from a univariate skew normal distribution. This incoherence causes awkward interpretation for data analysis in practice, especially in the development of the sampling distribution theory. In this paper, we propose a refined model that is coherent with the joint distribution of the univariate skew :normal random sample, for multivariate skew normal data. The proposed model extends and strengthens the multivariate skew model described in Azzalini (1985, Scandinavian Journal of Statistics, 12, 171-178). We present a stochastic representation for the newly proposed model, and discuss a bivariate setting, which confirms that the newly proposed model is more plausible than the one given by Azzalini and Dalla Valle (1996, Biometrika, 83, 715-726).