The numerical simulation of geophysical flows has many applications and poses interesting questions of analysis and implementation. Here we are interested in the simulation of hydrostatic flows for which the vertical variation of the horizontal velocity can not be neglected. For such a flow, shallow water equations are not satisfactory and full 3d Euler or Navier-Stokes equations are unnecessary costly. Some years ago, we introduced a new multi-layer, or more precisely multi-velocities, model. In this talk, we describe a well balanced numerical scheme that is positive under a classical CFL condition. We will also construct the solution of the Riemann problem for the bilayer case.