Over the last decade, increasing efforts have been directed towards the development of marine models using unstructured meshes. Such models allow for an accurate representation of the coastlines (e.g., islands, narrow straits) and the bathymetry. The mesh can easily be refined in regions of interest or coarsened in those regions where the dynamics is less demanding. Finally, unstructured meshes set up in spherical geometry allow to avoid singularities at the poles, rendering those techniques potentially useful for global ocean modelling. We present the latest version of {\it SLIM} (Second-generation Louvain-la-Neuve Ice-ocean Model, \textsf{http://www.climate.be/SLIM}), a three-dimensional, baroclinic, unstructured-mesh, finite-element ocean model. Herein, the emphasis is put on the mode splitting, which treats implicitly the two-dimensional mode while the three-dimensional part is treated explicitly. This method allows the representation of fast gravity waves with a reasonable time step, and generates a system much faster to solve than a fully implicit three-dimensional approach. It is well known that classical implementations of mode splitting resort to a weak coupling between two- and three-dimensional modes which can be the cause of inconsistencies. We introduce a new coupling method in which both modes are solved simultaneously. The new approach enables to ensure consistency: the mode splitting algorithm can be viewed as a specific time-stepping scheme, rather than a spatial approximation. Efficient implicit-explicit time discretizations ({\it IMEX}) schemes are used, which ensure stability and an optimal computational time. The model is validated by its application to a set of three-dimensional baroclinic simulations, such as the well-known overflow DOME test-case.