A mesoscale numerical framework is presented for modeling progressive failure in composite laminates under high-cycle fatigue loading conditions. The recently proposed fatigue cohesive zone model (CZM) by Dávila, which takes into account initiation and propagation of cracks and requires limited fatigue input parameters, has been employed for simulating matrix cracking and delamination. In our contribution, the existing formulation of the fatigue CZM has been improved by using a fully implicit scheme for integrating the fatigue damage rate function. Furthermore, the traction update with fatigue damage is consistently linearized to obtain a consistent material tangent stiffness matrix. In addition, the fatigue CZM is combined with XFEM for modeling discrete matrix cracks that can initiate at arbitrary locations with a predefined crack spacing in each ply. The performance of the improved formulation of Dávila's fatigue CZM is demonstrated with double cantilever beam simulations. It is shown that the use of the implicit fatigue damage update and consistent tangent stiffness allows for larger cycle increments which significantly improves accuracy, efficiency and robustness of the numerical framework. The numerical framework is applied to the simulation of progressive failure of an open-hole [±45]-laminate. It is shown that the numerical model is capable of simulating the complete failure process of matrix cracking and delamination. Furthermore, the model is validated with experiments on notched laminates from literature. Moreover, the applicability of the model to simulate the mechanical response of ductile matrix material systems (thermoplastics) is explored.