We present in this work a new discretization of the multi-layer quasi-geostrophic (QG) model that relies on implicit dissipation rather than additional explicit dissipation. It is first based on the staggered discretization of the potential vorticity (PV) and the stream-function in order to solve the PV advection with a finite volume method. This ensures the exact material conservation of the PV. We compute PV fluxes with a WENO-5 interpolation whose implicit dissipation replaces the usual explicit (hyper-)viscous dissipation. We propose a new method for solving reversibly the elliptic equation which is nontrivial with this staggered discretization. The presented discretization does not require the tuning of any additional parameter, \textit{e.g.} additional hyper-viscosity. We test the proposed method on a challenging idealized wind-driven double-gyre configuration at eddy-resolving, eddy-permitting, and non-eddy-resolving resolutions. In the eddy-permitting and non-eddy-resolving resolutions, our method produces an eastward jet contrary to usual discretizations without parametrization. Moreover, our method produces statistics that have a stronger coherence across resolutions than usual discretizations. We release a very short, concise, and efficient PyTorch implementation of our method to facilitate future data assimilation or machine-learning developments upon this new discretization.