We broaden the well-known Gneiting class of space-time covariance functions by introducing a very general parametric class of fully nonseparable direct and cross-covariance functions for multivariate random fields, where each component has a spatial covariance function from the Matern family with its own smoothness and scale parameters and, unlike most of currently available models, its own correlation function in time. We present sufficient conditions that result in valid models with varying degrees of complexity and we discuss the parameterization of those. Continuous-in-space and discrete-in-time simulation algorithms are also given, which are not limited by the number of target spatial coordinates and allow tens of thousands of time coordinates. The application of the proposed model is illustrated on a weather trivariate dataset over France. Our new model yields better fitting and better predictive scores in time compared to a more parsimonious model with a common temporal correlation function.