In this work, we develop a Reynolds stress model along the line of the approach presented by Commun. Nonlinear Sci. Numer. Simul. 9, 543 (2004)]], aiming to assess the role and contribution of the mean spin tensor in turbulence modeling. Here, the constitutive functional for the Reynolds stress depends on the mean spin tensor as well as the mean stretching tensor and its Jaumann derivative, the turbulent kinetic energy K , and the turbulent dissipation rate epsilon , which is at the complexity level of p=1,m=1 , and n=0 of a rate-type constitutive equation for the Reynolds stress proposed in the aforementioned paper. The explicit form for the Reynolds stress is obtained with recourse to the representation theorem and the theory of invariants developed in modern rational continuum mechanics, and, as an approximation, a nonlinear cubic K-epsilon model is worked out in which the model coefficients are analytically identified based on the experimental results of Tavoularis and Corrsin [J. Fluid Mech. 104, 311 (1981)]]. In addition, numerical results based on this model, in the forms of employing the Jaumann derivative and the Oldroyd derivative, respectively, for homogeneous turbulent shear flow and fully developed turbulent flow over a backward-facing step, are presented in comparison with those obtained based on a few previously proposed linear and nonlinear K-epsilon models, showing reasonably good agreement with the experimental results and the DNS data concerned and a better performance than the previously developed quadratic models.