(1) 
This first equation is basicly the same we used in the SEM method, applied to the vorticity field instead of the velocity field. What we need now to calculate the velocity field from the generated vorticity.
(2) 
Using our divergence free hypotesis () we are allowed to delete the first term in the right part in equation (2), and we get:
(3) 
This is the equation we have to solve to get our velocity field. To solve it we can use the Green's function and the BiotSavart law that relate circulation and velocity. All together we get as the final result:
(4) 
Using (1) this formula, after some passages, can be simplified into:
(5) 
(6) 
The main idea is to develop the old SEM method, forcing one component to satisfy the divergence free condition. To pursue such a scope, we decided to force the 3rd component to be Divergece Free, even if such a constraint implies loss of part of the Reynolds Stress information.
The first problem I dealt with is the stability one. In fact, since we are basically solving to calculate the third component, and since this is a differential equation, we have to pay a lot of care about it.
The first algorithm used is the following one (note the pedix N,S,E,W,T,B have the meaning of the neighbours point):


As we can see from the following graph, this new algorithm did not solve the stability problem. In fact the new error propagation now is backward.
I tried then to run a simulation considering all the 27 sourroungind points (basicalli all the points forming a cube around the point P), but with mainly the same results.
Another problem that i saw is that the instantaneous velocity calculated has a shape totally different from the other two random generated components.