This work presents a finite element implementation to treat the Hydromechanical Coupling (HM) in fractured rock masses under the framework of the so-called ‘equivalent continuum’ approach. The multilaminar concept, introduced by Zienkiewicz and Pande, is used to simulate the mechanical behaviour of both the intact rock and the families of fractures. In that concept, the non-linearities in the constitutive relations are dealt by means of fictitious viscoplasticity. In the present implementation, the mechanical behaviour of the fractures is modelled by means of Barton–Bandis model. The shear stress/shear displacement/dilatancy relationship is modelled as viscoplastic and the normal stress/normal displacement as non-linear viscoelastic. Flow along fractures is considered to occur as a sequence of permanent states. The permeability tensor of the equivalent continuum is determined from the hydraulic apertures, in accordance of Barton et al. From the numerical point of view, the basic aim of the work is the implementation of an efficient scheme to solve the above described problem. This is done by designing a self-adaptive time step control, transparent to the user, which determines the highest possible time step while assuming the conditions of precision, stability and convergence. The paper presents the numerical details of such scheme together with validation/comparative examples and the results obtained on the analysis of the fractured rock foundation of a hypothetical dam. © 1998 John Wiley & Sons, Ltd.