A new method is presented to generate reduced order models (ROMs) in Fluid Dynamics problems. The method is based on the expansion of the flow variables on a Proper Orthogonal Decomposition (POD) basis, calculated from a limited number of snapshots, which are obtained via Computational Fluid Dynamics (CFD). Then, the POD-mode amplitudes are calculated as minimizers of a properly defined overall residual of the equations and boundary conditions. The residual can be calculated using only a limited number of points in the flow field, which can be scattered either all over the whole computational domain or over a smaller projection window. This means that the process is both computationally efficient (reconstructed flow fields require less than 1% of the time needed to compute a full CFD solution) and flexible (the projection window can avoid regions of large localized CFD errors). Also, various definitions of the residual are briefly discussed, along with the number and distribution of snapshots, the number of retained modes, and the effect of CFD errors, to conclude that the method is numerically robust. This is because the results are largely insensitive to the definition of the residual, to CFD errors, and to the CFD method itself, which may contain artificial stabilizing terms. Thus, the method is amenable for practical engineering applications.