The new formulation for conformation-tensor based viscoelastic fluid models, written in terms of the contravariant deformation tensor as introduced by Hutter et al. (2018), has been approached numerically. It has been implemented in a finite element method framework with the DEVSS-G/SUPG method for stabilisation. A stress-implicit as well as a stress-explicit formulation is presented that allows to integrate the non-linear flow equations, with or without a solvent contribution, forward in time with second-order accuracy. The time-dependent stability in shear flow is maintained using the contravariant deformation, which is tested by solving a perturbed planar Couette flow of an upper-convected Maxwell (UCM) fluid. The stability and accuracy of our implementation has furthermore been tested by solving the flow around a cylinder with confining walls. The new formulation turns out to be more stable as compared to the formulation in terms of the conformation tensor. It enables to simulate the flow of a Giesekus fluid with non-linear parameter alpha = 0.01 and viscosity ratio nu = eta(s/)eta(0) = 0.59 way beyond the Weissenberg number at which the High Weissenberg Number Problem manifests itself using the standard formulation. The stability appears to be similar to the log-conformation representation. The computed stress profiles up to a Weissenberg number of Wi = 0.6 for an Oldroyd-B fluid with viscosity ratio nu = 0.59, compare well with benchmark results from other studies, including other finite element, finite volume and spectral element methods.