conservation [@tt].\
For the discretization in space, we rely on the approach in [@acd]. We introduce the finite differences operators $$\label{eq:finite_diff}
(D_1^+ U )_{i,j} = \frac{ U_{i+1,j}-U_{i,j}   } {h} \quad \hbox{and }\quad  (D_2^+ U )_{i,j} = \frac{ U_{i,j+1}-U_{i,j}   } {h},$$ and define $$[D_h U]_{i,j} =\left((D_1^+ U )_{i,j} , (D_1^+ U )_{i-1,j}, (D_2^+ U )_{i,j}, (D_2^+ U )_{i,j-1}\right) ^T.\label{eq:discrete_gradient}$$ For the discretization of the Laplace operator, we consider a standard five points approximation $$(\Delta_h U)_{i,j} = -\frac{1}{h^2} (4 U_{i,j} - U_{i+1,j} - U_{i-1,j} - U_{i,j+1} - U_{i,j-1}).$$ In order to approximate the Hamiltonian $H$ in equation , we consider a numerical Hamiltonian $g: {\mathbb{T}}^2 \times {\mathbb{R}}^4\to {\mathbb{R}}$, $(x,q_1,q_2,q_3,q_4)\mapsto g\left(x,q_1,q_2,q_3,q_4\right)$ satisfying the following conditions:

-   $g$ is nonincreasing with respect to $q_1$ and $q_3$, and nondecreasing with respect to $q_2$ and $q_4$.

-   $g$ is consistent with the Hamiltonian $H$, i.e. $$g(x,q_1,q_1,q_2,q_2)=H(x,q), \quad \forall x\in {\mathbb{T}}^2,   q=(q_1,q_2)\in {\mathbb{R}}^2.$$

-   $g$ is of class $C^1$.

-   the function $(q_1,q_2,q_3,q_4) \mapsto g(x,q_1,q_2,q_3,q_4)$ is convex.

For the discretization of the coupling term $f[m]$ in the Hamilton-Jacobi-Bellman equation, denoted by $m_h$ the piecewise constant function with value $M_{i,j}$ in the square $|x-x_{i,j}|_\infty \le h/2$, we assume that $f[m_h]$ can be calculated in practice. Then, the approximation 