  Being able to account for the missing mixing in stellar radiative zones is a
key step toward a better understanding of stellar evolution. Zahn (1974) argued
that thermally diffusive shear-induced turbulence might be responsible for some
of this mixing. In Part I and Part II of this series of papers we showed that
Zahn's (1974, 1992) mixing model applies when the properties of the turbulence
are local. But we also discovered limitations of the model when this locality
condition fails, in particular near the edge of a turbulent region. In this
paper, we propose a second-order closure model for the transport of momentum
and chemical species by shear-induced turbulence in strongly stratified,
thermally diffusive environments (the so-called low P\'eclet number limit),
which builds upon the work of Garaud \& Ogilvie (2005). Comparison against
direct numerical simulations (DNSs) shows that the model is able to predict the
vertical profiles of the mean flow and of the stress tensor (including the
momentum transport) in diffusive shear flows, often with a reasonably good
precision, and at least within a factor of order unity in the worst case
scenario. The model is sufficiently simple to be implemented in stellar
evolution codes, and all the model constants have been calibrated against DNSs.
While significant limitations to its use remain (e.g. it can only be used in
the low P\'eclet number, slowly rotating limit), we argue that it is more
reliable than most of the astrophysical prescriptions that are used in stellar
evolution models today.
