We present a computationally efficient multiscale method for preparing equilibrated, isotropic long-chain model polymer melts. As an application, we generate Kremer-Grest melts of 1000 chains with 200 entanglements and 25 000-2000 beads/chain, which cover the experimentally relevant bending rigidities up to and beyond the limit of the isotropic-nematic transition. In the first step, we employ Monte Carlo simulations of a lattice model to equilibrate the large-scale chain structure above the tube scale while ensuring a spatially homogeneous density distribution. We then use theoretical insight from a constrained mode tube model to introduce the bead degrees of freedom together with random walk conformational statistics all the way down to the Kuhn scale of the chains. This is followed by a sequence of simulations with carefully parameterized force-capped bead-spring models, which slowly introduce the local bead packing while reproducing the larger-scale chain statistics of the target Kremer-Grest system at all levels of force-capping. Finally, we can switch to the full Kremer-Grest model without perturbing the structure. The resulting chain statistics is in excellent agreement with literature results on all length scales accessible in brute-force simulations of shorter chains.