We describe recent progress in the modelling of realistic equilibrium configurations of rotating superfluid neutron stars, in a fully general relativistic framework. We compute stationary and axisymmetric configurations of neutron stars composed of two interpenetrating and interacting fluids, namely superfluid neutrons and charged particles (protons and electrons), rotating with different rotation rates around a common axis. Two different realistic equations of state are considered. As a first application, we propose a simple bulk model for pulsar glitches, seen as angular momentum transfers between the two fluids through mutual friction force. From a series of equilibrium states, we compute the evolution in time of the properties of a neutron star during the rise period of a glitch. This enables us to infer characteristic features relative to glitches, such as spin-up timescales, that could be compared with future accurate observations in order to put some constraints on the interior of neutron stars.