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, have been studied by numerical modelling in the case of a rhombohedral lead zirconate-titanate ferroelectric/ferroelastic thin film with an anisotropic misfit 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 in-plane 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. The flexoelectric effect exhibits a weak influence on the phase diagram of multi-domain states when compared with the phase diagram of single-domain states. Its effect, however, becomes significant in the case of skyrmion topological states, which spontaneously form near the film surface when compressive misfit strains are applied. 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 pronounced for topologically nontrivial structures.