Using dt(u) as the time derivative term dictates that the units of the equation be acceleration (length/time^2). Consistent with this demand, the equation will be dimensionally correct if F is interpreted as (force/mass) (which is equivalent to length/time^2), and P is interpreted as kinematic pressure (pressure/rho). In this latter term, it would be more transparent if we interpreted P as standard pressure and divided by rho, as you have indicated.
In the Boussinesq approximation, rho=rho0*[1+alpha(T-T0)], there is always some latitude as to when the variation of pressure can be ignored relative to the other terms. In this context, one can argue that moving the density in and out of the grad is justifiable.
However, FlexPDE does not have any built-in interpretation of units of measure. As in all FlexPDE scripts, it is up to the user to be sure that his interpretation of the units of the parameters of the equation result in dimensionally consistent equations.
FlexPDE example scripts are meant as demonstrations of the techniques of equation description, and not as definitive formulations of physical processes. Users are at all times free to formulate equations in a manner consistent with their understanding of their own area of application.