The geostatistical analysis of multivariate data involves choosing and fitting theoretical models to the empirical matrix. This paper considers the specific case of the model of linear coregionalization, and describes an automated procedure for fitting models, that are adequate in the mathematical sense, using a least-squares like technique. It also describes how to decide whether the number of parameters of the cross-variogram matrix model should be reduced to improve stability of fit. The procedure is illustrated with an analysis of the spatial relations among the physical properties of an alluvial soil. The results show the main influence of the scale and the shape of the basic models on the goodness of fit. The choice of the number of basic models appears of secondary importance, though it greatly influences the resulting interpretation of the coregionalization analysis.