The numerical implementation of a recently developed thermomechanical constitutive model for fine-grained soils based on hyperelasticity-hyperplasticity theory (Golchin et al. 2020), is presented. A new unconventional implicit stress return mapping algorithm, compatible with elasticity derived from Gibbs (complementary) energy potential, in strain invariant space, is designed and the consistent tangent operator for use in boundary value problems (such as in the finite element method) is derived. It is shown that the rate of convergence of the stress integration algorithm is quadratic. The numerical results are in good agreement with available data from thermomechanical element tests found in literature.