Large-scale reservoir modeling and simulation of gas reservoir flows in fractured porous media is currently an important topic of interest in petroleum engineering. In this paper, the dual-porosity dual-permeability (DPDP) model coupled with the Peng-Robinson equation of state (PR-EoS) is used for the mathematical model of the gas reservoir flow in fractured porous media. We develop and study a parallel and highly scalable reservoir simulator based on an adaptive fully implicit scheme and an inexact Newton type method to solve this dual-continuum mathematical model. In the approach, an explicit-first-step, single-diagonal-coefficient, diagonally implicit Runge–Kutta (ESDIRK) method with adaptive time stepping is proposed for the fully implicit discretization, which is second-order and L-stable. And then we focus on the family of Newton–Krylov methods for the solution of a large sparse nonlinear system of equations arising at each time step. To accelerate the convergence and improve the scalability of the solver, a class of multilevel monolithic additive Schwarz methods is employed for preconditioning. Numerical results on a set of ideal as well as realistic flow problems are used to demonstrate the efficiency and the robustness of the proposed methods. Experiments on a supercomputer with several thousand processors are also carried out to show that the proposed reservoir simulator is highly scalable.