We propose a new copula model for replicated multivariate spatial data. Unlike classical models that assume multivariate normality of the data, the proposed copula is based on the assumption that some factors exist that affect the joint spatial dependence of all measurements of each variable as well as the joint dependence among these variables. The model is parameterized in terms of a cross-covariance function that may be chosen from the many models proposed in the literature. In addition, there are additive factors in the model that allow tail dependence and reflection asymmetry of each variable measured at different locations, and of different variables to be modeled. The proposed approach can therefore be seen as an extension of the linear model of coregionalization widely used for modeling multivariate spatial data. The likelihood of the model can be obtained in a simple form and, therefore, the likelihood estimation is quite fast. The model is not restricted to the set of data locations, and using the estimated copula, spatial data can be interpolated at locations where values of variables are unknown. We apply the proposed model to temperature and pressure data, and we compare its performance with that of a popular model from multivariate geostatistics.