I proposed things to try (tighten the residual tolerance, switching to first order scheme) to help trouble shooting the issue in the other post.
As far as “an efficient way to quickly obtain the values of each term in the two-fluid model equations for every grid cell”, you will have to go in the code and spend quite a bit of time to decide which terms you want to look at. You can start with ‘iterate.f’ and 'solve_vel_star.f` and follow the iterative scheme. Then if the data is in an array with cell index IJK, you can save it into a vtk file by modifying vtk_out.f. However, I doubt this will help much, due to the instability occurring at random time and the fine mesh.
My suggestion is still to play with the numerics tolerances, try first order, or coarsen the mesh a bit.