Questions regarding "integrable systems" 
Consider a smooth differential equation on the plane
  $$
x'=g(x,y),\quad y'=h(x,y).
$$
  Suppose there exists a function $D(x,y)$ such that
  $$
(Dg)_x+(Dh)_y=0. 
$$ 
  Then $D$ is an integrating factor and the system is integrable. 

A quick search for "integrable system" on Google returns results not satisfying. 
Could anyone explain what the last sentence in the argument above means?

[Added:]
The question is motivated by reading a paper about the Bendixson-Dulac Theorem. In particular,
$$
g(x,y)=ax+bx^2+cxy,\quad h(x,y)=dy+exy+fy^2
$$
and $D(x,y)=x^ry^s$ for some $r,s$.

[Added:] I asked this question in MO. I don't understand though, there is an answer there:  

$X=g\partial_x + h\partial_y$ is the vector field whose flow lines are wanted. 
  
  $\omega=hdx - gdy$ is a 1-form with kernel the span of $X$. 
  
  Also $D\omega$ has kernel the span of $X$ for any function $D$ which does not vanish anywhere. If $d(D\omega)=0$ (this is your condition) then $D\omega$ is a closed 1-form, thus exact on simply connected sets. So $D\omega = dF$ for a function $F$ which can easily be computed by line integrals.
  Thus the wanted flow lines are contained in the level sets of $F$. 
  
  Finally,
  the time dependence of the flow has to be computed extra. 

I would really appreciate it if anyone could explain what that answer means (in a more "elementary" way) here.
 A: The system $x'=gD$, $y'=hD$ has the same solution curves as the system $x'=g$, $y'=h$, since the velocity vector field $(gD,hD)$ is parallel (at every point) to the vector field $(g,h)$. (But the solutions to the two systems of course follow these solution curves with different speeds.)
Now, if your condition $(-hD)_y=(gD)_x$ is satisfied, then there is a function $F(x,y)$ such that 
$$
F_x=-hD
,\quad
F_y=gD
.
$$
(Recall that the condition for a system of the form $F_x=a$, $F_y=b$ to have a solution is the equality of mixed derivatives, $(F_x)_y=(F_y)_x$, i.e., $a_y=b_x$. This is always necessary, and it's sufficient if the domain is simply connected.)
This function $F$ will be a constant of motion of the system $x'=gD$, $y'=hD$.
In other words, it will have the property that it's constant on solution curves $(x(t),y(t))$,
as is easy to check using the chain rule:
$$
\frac{d}{dt}F(x(t),y(t))
= \frac{\partial F}{\partial x} \, \frac{dx}{dt} + \frac{\partial F}{\partial y} \, \frac{dy}{dt}
= F_x gD + F_y gD
= F_x F_y + F_y (-F_x) = 0
.
$$
If we know $F$, then we know its level curves, and thus we know the solution curves of the system $x'=gD$, $y'=hD$, and (as we noted initially) these are also the solution curves of the original system $x'=g$, $y'=h$. In this sense the system is "integrable"; it has a constant of motion (also called a "first integral"). By integrating the equations $F_x=-hD$, $F_y=gD$ (going from the derivatives of $F$ to $F$ itself), we solve our ODEs geometrically in the sense that we know what the solution curves look like.
(However, there are still some more computations to do if we want to find the exact time dependence, i.e., finding out exactly where on the curve the solution is at a given time $t$.)
