This work aims to calculate interlaminar stress distribution through the thickness of multilayered composite shell structures by employing a novel nonlinear layer-wise shell finite element formulation. Adapting the Mindlin– Reissner theory in each layer, the shear-deformable layer-wise shell element presents the interlaminar shear stress distributions by increasing the number of layers. The interlaminar normal stress distribution is then determined using the finite difference solution of the general form of equilibrium equation in the non-orthogonal curvilinear grid along the Gaussian points. Two boundary conditions at the bottom and the top surfaces are satisfied by adopting the linear Lagrange interpolation function. The developed formulation is assessed through some illustrative problems solved using a proprietary finite element computer program. The results compare very well with those available in the literature and those obtained by simulations with the commercial finite element software Ansys.