Khawaja (Khawaja and Moatamedi 2018) wrote a clear explanation of this algorithm.
The continuous equations of looks like this:
\begin{aligned} & \frac{\partial u}{\partial t}=-\frac{u \partial u}{\partial x}-\frac{u \partial v}{\partial y}-\frac{\partial p}{\partial x}+\frac{\mu}{\rho}\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right) \\ & \frac{\partial v}{\partial t}=-\frac{v \partial v}{\partial y}-\frac{v \partial u}{\partial x}-\frac{\partial p}{\partial y}+\frac{\mu}{\rho}\left(\frac{\partial^2 v}{\partial y^2}+\frac{\partial^2 v}{\partial x^2}\right) \end{aligned}
The first step is to solve the discretized momentum equations at the current time step:
\begin{gathered} u^{t+1}(i, j)=u^t(i, j)+ \\ \Delta t\left(\text { function }\left(u^t, v^t, \rho, \mu, \Delta x, \Delta y\right)-\frac{(p(i, j)-p(i, j+1))}{\Delta x}\right) \\ v^{t+1}(i, j)=v^t(i, j)+ \\ \Delta t\left(\text { function }\left(u^t, v^t, \rho, \mu, \Delta x, \Delta y\right)-\frac{(p(i, j)-p(i+1, j))}{\Delta y}\right) \end{gathered}