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<=Teff(K)<=4000, 0.50<=C/O<0.99, 0<=logg<=5, [Fe/H]=0., -0.5dex, 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 (1980MNRAS.190..441P) that ZrO bands in warm S stars (Teff > 3200K) 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 Teff, 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 Teff=3200K, a change of C/O from 0.5 to 0.99 has a strong impact on V-K (2mag). Conversely, a range of 2 mag in V-K corresponds to a 200K shift along the (Teff, V-K) relationship (for a fixed C/O value). Hence, the use of a (Teff, 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.