A simple numerical scheme is presented for the construction of three-integral phase-space distribution functions for oblate galaxy models with a gravitational potential of St\"{a}ckel form, and an arbitrary axisymmetric luminous density distribution. The intrinsic velocity moments can be obtained simultaneously with little extra effort. The distribution of the inner and outer turning points of the short-axis tube orbits that are populated can be specified freely, and is chosen in advance. The entire distribution function is then derived from the density by an iterative scheme that starts from the explicitly known distribution function of the thin-orbit (maximum streaming) model, in which only the tubes with equal inner and outer turning points are populated. The versatility and limitations of this scheme are illustrated by the construction of a number of self-consistent three-integral flattened isochrone models of Kuzmin--Kutuzov type, and by investigation of special cases where the scheme is tractable analytically. This includes the behaviour of the distribution functions in the outer regions of the models. The scheme converges rapidly for models containing orbits with ratios of the outer to inner turning point as large as ten, and is particularly suited for the construction of tangentially anisotropic flattened models, self-consistent as well as non-consistent. The algorithm simplifies in the disk and spherical limit, and can be generalized to triaxial models.