S-type stars are late-type giants whose atmospheres are enriched in carbon and s-process elements because of either extrinsic pollution by a binary companion or intrinsic nucleosynthesis and dredge-up on the thermally-pulsing asymptotic giant branch. A grid of MARCS model atmospheres has been computed for S stars, covering the range 2700 ≤ T eff (K) ≤ 4000, 0.50 ≤ C/O ≤ 0.99, 0 ≤ log g ≤ 5, [Fe/H] = 0., −0.5 dex, and [s/Fe] = 0, 1, and 2 dex (where the latter quantity refers to the global overabundance of s-process elements). The MARCS models make use of a new ZrO line list. Synthetic spectra computed from these models are used to derive photometric indices in the Johnson and Geneva systems, as well as TiO and ZrO band strengths. A method is proposed to select the model best matching any given S star, a non-trivial operation since the grid contains more than 3500 models covering a five-dimensional parameter space. The method is based on the comparison between observed and synthetic photometric indices and spectral band strengths, and has been applied on a vast subsample of the Henize sample of S stars. Our results confirm the old claim by Piccirillo (1980, MNRAS, 190, 441) that ZrO bands in warm S stars (T eff > 3200 K) are not caused by the C/O ratio being close to unity, as traditionally believed, but rather by some Zr overabundance. The TiO and ZrO band strengths, combined with V − K and J − K photometric indices, are used to select T eff , C/O, [Fe/H] and [s/Fe]. The Geneva U − B 1 and B 2 − V 1 indices (or any equivalent) are good at selecting the gravity. The defining spectral features of dwarf S stars are outlined, but none is found among the Henize S stars. More generally, it is found that, at T eff = 3200 K, a change of C/O from 0.5 to 0.99 has a strong impact on V − K (2 mag). Conversely, a range of 2 mag in V − K corresponds to a 200 K shift along the (T eff , V − K) relationship (for a fixed C/O value). Hence, the use of a (T eff , V − K) calibration established for M stars will yield large errors for S stars, so that a specific calibration must be used, as provided in the present paper. Using the atmospheric parameters derived by our method for the sample of Henize S stars, we show that the extrinsic-intrinsic dichotomy among S stars reveals itself very clearly as a bimodal distribution in the effective temperatures. Moreover, the increase of s-process element abundances with increasing C/O ratios and decreasing temperatures is apparent among intrinsic stars, confirming theoretical expectations.