We present a parallel domain decomposition method based on a hybrid direct-iterative solver for 3D modelling of visco-acoustic waves in the frequency domain. The modelling method was developed as a modelling engine for frequency-domain full waveform inversion. Frequency-domain seismic modelling reduces to the solution of a large and sparse system of linear equations, resulting from the discretization of the heterogeneous Helmholtz equation. Our approach to high-performance, scalable solution of large sparse linear systems in parallel scientific computing is to combine direct and iterative methods. Such a hybrid approach exploits the advantages of both direct and iterative methods. The iterative component allows us to use a small amount of memory and provides a natural way for parallelization. The direct part provides its favorable numerical properties. The domain decomposition is based on the algebraic Schur complement method, which allows for the iterative solution of a reduced system, the solution of which is the pressure wavefield at the interfaces between the subdomains. Once the interface unknowns have been computed, the wavefield at the interior of each subdomain is efficiently computed by local substitutions. The reduced Schur complement system is solved with the global minimum residual method (GMRES) and is preconditioned by an algebraic additive Schwarz preconditioner. A direct solver is used to factorize the local impedance matrices defined on each subdomain. Simulations are performed in the overthrust and the salt models for frequencies up to 12.5~Hz. The numerical experiments show that the number of iterations increases linearly with the number of subdomains for a given computational domain but that the elapsed time of the iterative resolution remains almost constant. The number of iterations also increases linearly with the frequencies, when the grid interval is adapted to the frequencies and the size of the subdomains is kept constant over frequency. Although the hybrid approach allows one to tackle larger problems than the direct-solver approach, further improvements are needed to mitigate the computational burden of the iterative component of the hybrid solver within the framework of multisource modelling. On the numerical side, the use of block iterative solvers and of incremental two-level deflating preconditioners, and on the parallel implementation side the use of two levels of parallelism in the domain decomposition method should allow us to mitigate this computational burden.