The classical shallow water equations have a limited range of application because they cannot model vertical effects. The multi-layer approach allows to recover some vertical variations of the flow while keeping some simplification of the shallow water model. The computational costs obviously grow with the number of layers, which is often around 50 in ocean simulations. With a large computational domain and/or long time scales and a focus on the evolution on tracers rather than on the free surface, a strategy should be employed to reduce the computational cost. In this contribution we focus on the barotropic-baroclinic splitting which is employed in (numerical) ocean models. The fast barotropic gravity waves are treated in a vertically averaged manner. The slower baroclinic dynamic is fully multidimensional but has much larger time step. We reformulate this strategy in the nonlinear multilayer framework in terrain-following coordinates and write it as an exact splitting which behaves well regarding the total energy. No filters or correction are needed. The barotropic step gathers the evolution of the free surface and the mean velocity with an accurate and well-balanced one layer shallow water model. The baroclinic step includes the vertical exchange between the layers and an adjustment of the velocities around their mean vertical value. In low Froude simulation the computational cost is divided by the total number of layer without degrading the accuracy. Currently this work deals with the constant density case, but in ongoing work we are extending the barotropic-baroclinic splitting to the variable density case in order to model situations such as coastal upwelling.