Using the Landau-Ginzburg-Devonshire theory, an influence of the misfit strain and surface screening charges, as well as the role of the flexoelectric effect, has been studied by numerical modelling in the case of a rhombohedral lead zirconate-titanate ferroelectric/ferroelastic thin film with an anisotropic mismatch produced by a substrate. It was established that the magnitude and sign of the misfit strain influence the domain structure and predominant directions of the polarization vector, providing misfit-dependent phases with different favourable polarization components. Whilst strong enough compressive misfit strains favour a phase with an orthorhombic-like polarization directions, strong tensile misfits only yield inplane polarization components. The strength of surface screening is seen to condition the existence of closure domain structures and, by increasing, supports the single-domain state depending on the value of the misfit strain. A crucial role of the flexoelectric effect is revealed as it saves the polar state of the film at high tensile misfit strains. Without the flexoelectric effect, the growing tensile misfit strain eventually suppresses the polarization, and hence any domain structure. Cooperative influence of the misfit strain, surface screening charges and temperature can set a thin rhombohedral ferroelectric film into a number of different polar and structural states, whereby the role of the flexoelectric effect is significant for establishing its phase diagram.