The need for flexible and efficient multi-porosity, reservoir-simulation capabilities has never been greater. Much of the world’s oil reserves are contained in highly variable fractured reservoirs. Moreover, there is now unprecedented interest in the simulation of unconventional gas reservoirs, where up to four porosity types may be required. This paper discusses the design of an N-porosity, full-featured reservoir simulator capable of handling these diverse scenarios. Essential to the design is the programming data-structure paradigm of the SubGrid. SubGrids represent collections of simulation nodes within spatial regions of interest and encapsulate all data and variables required. Any number of associated SubGrids may coexist, representing N-porosities. Inter- and intra-SubGrid connections may be arbitrarily specified. Because the linear solver operates on the associated SubGrids as distinct entities, the design has tremendous flexibility. For example, regions with simple connectivity, such as in dual-porosity, single-permeability models, can have simple pre-solver elimination performed while other regions with uniform intra-porosity connectivity require a merged SubGrid solution. The scheme efficiently treats regions with different numbers of porosities, levels of interconnectivity, and levels of implicitness. Applications may be applied across multiple reservoirs simultaneously. Treatment of missing "fracture" zones has always been an important feature in the dual-porosity, single-permeability simulation. The traditional approach is to treat the matrix within missing fracture zones effectively as though it were an interconnected portion of the fracture porosity. This preserves connectivity within the fracture-free zones and the more computationally efficient single-permeability solver pre-elimination step. A second missing "fracture" scheme is introduced, which differs in the way that the connections between the fracture and matrix blocks are computed at the boundary of fracture-free zones. Flexible transfer-function application and porosity dependency during equilibrium initialization are discussed.