  This paper considers a multivariate spatial random field, with each component
having univariate marginal distributions of the skew-Gaussian type. We assume
that the field is defined spatially on the unit sphere embedded in
$\mathbb{R}^3$, allowing for modeling data available over large portions of
planet Earth. This model admits explicit expressions for the marginal and cross
covariances. However, the $n$-dimensional distributions of the field are
difficult to evaluate, because it requires the sum of $2^n$ terms involving the
cumulative and probability density functions of a $n$-dimensional Gaussian
distribution. Since in this case inference based on the full likelihood is
computationally unfeasible, we propose a composite likelihood approach based on
pairs of spatial observations. This last being possible thanks to the fact that
we have a closed form expression for the bivariate distribution. We illustrate
the effectiveness of the method through simulation experiments and the analysis
of a real data set of minimum and maximum temperatures.
