For the simulation of turbulent flows we require very fine grids to resolve all the relevant scales. This makes the involved computations expensive. Closure models alleviate this by modeling the effect of the small scales on the large scales, allowing for coarser grids. Recently, machine learning algorithms have shown promise to outperform the physics-based closure models. This hybrid physics/machine learning approach is typically formulated as d ̄u/dt = f ( ̄u) + NN( ̄u), where ̄u is a coarse grid representation of the solution, f a discretization of the known physics, and NN a neural network. However, this formulation does not guarantee stability. This can cause the system to explode. In this work we resolve this by introducing energy conservation into the neural network. We first note that the kinetic energy of the system E can be decomposed into two positive energies: E(u) = 1/2 ̄u^T ̄u + ∆( ̄u, u). However, the problem is that ∆ also depends on the true solution u which is not known. We resolve this by introducing a set of auxiliary variables s which represent ∆ on the coarse grid. The approximated total kinetic energy then becomes A = 1 2 a^T a, where a^T = [ ̄u s]. Finally, we present a novel neural network architecture which solely dissipates this energy. This guarantees stability and improves general performance of the closure model.