Within a simple SO(8) algebraic model, the coexistence between isoscalar and isovector pairing modes can be successfully described using a mean-field method plus restoration of broken symmetries. In order to port this methodology to real nuclei, we need to employ realistic density functionals in the pairing channel. In this article, we present an analytical derivation of matrix elements of a separable pairing interaction in Cartesian coordinates and we correct errors of derivations available in the literature. After implementing this interaction in the code hfodd, we study evolution of pairing gaps in the chain of deformed Erbium isotopes, and we compare the results with a standard density-dependent contact pairing interaction.