Multilayered plate structures, including composite laminates, are widely recognized as high-performance material solutions across numerous engineering applications, such as aerospace, naval, and automotive industries. The process of designing new allowables is, however, expensive and time consuming. For this reason achieving high-fidelity numerical simulations at a reduced computational cost is considered key for accelerating their design and market uptake. To tackle this challenge, the present study extends the multiscale 2D+ approach introduced in [1] to the nonlinear bending analysis of multilayered flat shells within an industrially relevant finite element framework. The proposed methodology employs dimensional reduction at both the macro- and meso-scales, substantially lowering computational expense while preserving the heterogeneous through-thickness response. The approach is based on a computational homogenization strategy in which flat shell structures with degenerated kinematics are used at the macro-scale and are coupled to a one-dimensional Representative Volume Element (RVE) through the thickness at the meso-scale. Representative numerical examples using a 2D+ flat shell formulation implemented in Abaqus demonstrate: (i) stress distribution predictions with accuracy comparable to full three-dimensional simulations, (ii) computational efficiency comparable to standard two-dimensional Equivalent Single Layer (ESL) theories, (iii) a fully non-intrusive implementation affecting only element level routines through user-defined elements (UEL), and (iv) the capability to naturally account for material nonlinearities.