Inequality constraints in calculus of variations $\def\d{\mathrm{d}}$It turns out that Yuri's answer to my earlier question, whilst correct (and I thank him for his effort), was not quite what I desired. I had not posed the question properly, so I have chosen to re-ask as I am still struggling to implement the strict inequality constraints. 
Let me begin by reposing the question. I wish to extremise 
$$Q = \int_0^h u \, \, \d y$$
subject to the constraints$$
B = B_0,\ M = M_0,$$
where
$$B = \int_0^h ug \, \d y, \ M = \int_0^h \left( u^2 + \int_0^y g(y^*) \, \d y^*\right) \, \d y.$$
An important new difference is that we have that $h$, the upper limit, can vary, and crucially we need
$$u(h) = 0, \ g(h) = 0, \\ \int_0^h u \, \d y> 0, \ \int_0^h g \, \d y = G(h) > 0.$$
Here is my attempt (mostly following @Yuri). We set $$ G(y) = \int_0^y g(y^*) \, \d y^*$$
such that $$B = \int_0^h uG' \, \d y, \ M = \int_0^h u^2 + G \, \d y$$ and hence the Lagrangian becomes
$$L(u,G,G',\lambda,\mu) = Q + \lambda(B-B_0) + \mu(M-M_0).$$ The variations of this (using Euler-Lagrange $\frac{\partial L}{\partial F} -  \frac{\d}{\d y}\frac{\partial L}{\partial F'} = 0$) become
$$1 + \lambda G' + 2\mu u = 0, \ \mu - \lambda u' = 0, \\ \int_0^h uG' = B_0, \ \int_0^h u^2 + G = M_0.$$
The second equation here tells us
$$u' = \frac{\lambda}{\mu} \Longrightarrow u = \frac{\lambda}{\mu}y + u_0.$$
Now since $h$ can vary we have two further transversality conditions 
$$\left.\frac{\partial L}{\partial G'}\right|_{y=h} = 0, \ \left.L - G' \frac{\partial L}{\partial G'}\right|_{y=h} = 0.$$
The first transversality condition tells us that 
$$\lambda u(h) = 0 \Longrightarrow u(h) = -\frac{\lambda}{\mu} \Longrightarrow u = -\frac{\mu}{\lambda}\left(h-y\right).$$
This is great as it satisifes $u(h) = 0$.
The second transversality condition tells us (using the first variation equation to give us $G$) that 
$$\mu G(h) = 0 \Longrightarrow h = \frac{\lambda}{\mu^2}.$$
(However I think this condition violates the requirement $G(h) > 0$)
We can now calculate everything we need,
$$B_0 = -\frac{1}{6\mu^3}, \ M_0 = \frac{\lambda}{2\mu^4} \\ \lambda = \frac{M_0}{3\left(6^{1/3}B_0^{4/3}\right)}, \ \mu = -\frac{1}{6^{1/3}B_0^{1/3}}, \ h = \frac{2^{1/3}M_0}{3^{2/3}B_0^{2/3}} \\ Q = \frac{M_0}{6^{1/3}B_0^{1/3}}.$$
If we now sub in some values $B_0 = 6, M_0 = 36.5$ I am especially happy with the answer for $Q \approx 11.05, h \approx 6.696$ as it verifies a separate result I have reached. The problem comes with the profile for $g$. If look at the profile for $u$ we have everything behaving as expected. $u(h) = 0$ and $\int_0^h u \, \, dy >0$.
 
However, the issue comes when looking at the profile for $g$. We can immediately see that $g(h) \neq 0$ (strangely $g(h/2) = 0$), and also we have $\int_0^h g \, \d y = 0$ which is not strictly bigger than 0.

So my question is, how can I build the machinery for the constraints on $g$ into the Lagrangian. I have read a little online about possibly using slack variables, but I am not sure how to implement them. Any pointers would be greatly appreciated.
Interestingly, if I were to transform the profile of $g$ to $\hat{g}$ such that $g(h/2) = \hat{g}(h) = 0$ and $g(0) = g^*, \hat{g}(0) = g^* / 2$ the solution for $\hat{g}$ would match the separate result I have reached identically.
EDIT
So I have a found a paper in which the author uses slack variables to impose a condition $\theta_z \geq 0$, my thinking is that this condition is equivalent to $G' >0$. Though in his derivation, his multiplier becomes a function of $z$. Something that I have never seen before.
EDIT 2
I am now not convinced the second transversality condition is correct. This condition $G(h) = \int_0^h g = 0$ is what is forcing the profile to be negative above $h/2$ (to keep the area 0). This actually contradicts the essential requirement that $\int_0^h g > 0$.
EDIT 3
I thought it would be helpful to add some further motivation. This method of calculus of variations as described above finds two solution profiles for $u$ and $g$ that are given by $$ u_1 = \frac{3 B_0 (h - y)}{M_0} \\ g_1 = \frac{3 \left(6 B_0^2 (h-y)-6^{1/3} B_0^{4/3} M_0\right)}{M_0^2} $$
which we can see satisfy the equations and maximise $Q$. 
I believe that the true solution I am seeking is given by profiles 
$$ u_2 = \frac{3 B_0 (h - y)}{M_0}, \ g_2 = \frac{9 B_0^2 (h-y)}{2 M_0^2}, \ h = \frac{2^{1/3}M_0}{3^{2/3}B_0^{2/3}}.$$
which can be shown to satisfy all the equations for $M,B$ and also give the same maximal $Q$. Furthermore, the profile for $g_2$ in particular satisfies the requirements that $g_2(h) = 0$ and $\int_0^h g_2 > 0$. 
It would appear to me that the $g_1, g_2$ profiles are equivalent solutions and I would expect that the calculus of variations to find both? I find it confusing that it doesn't.
I suppose one could frame the question in the way. What constraints needed to be added to the Lagrangian such that the calculus of variations gives the me the desired $g_2$ profile, rather than the unphysical $g_1$ profile?
 A: So, the task is to maximize the value of
$$Q=\int\limits_0^h udy,\qquad(1)$$
taking in account the integral relations of
$$\begin{cases}
B(uG') = \int\limits_0^h uG'dy = B_0\\
M(u^2+G) = \int\limits_0^h (u^2 + G)dy = M_0\\
G(y) = \int\limits_0^y g(y^*)dy^*\\
\int\limits_0^h g dy > 0, 
\end{cases}\qquad(2)$$
and boundary conditions of
$$\begin{cases}
u(h) = 0\\
g(h) = 0.
\end{cases}\qquad(3)$$
That task can be solved by using of Euler-Lagrange theorem of variations
$$\dfrac{\delta \mathcal L(y,f,f')}{\delta f} = \mathcal L'_f - \dfrac{d}{dy}\mathcal L'_{f'}.$$
Then, for the Lagrangian in the form of
$$\mathcal L(u, G, \lambda, \mu) = Q(u) + \lambda(B(uG')-B_0) + \mu(M(u^2+G)-M_0)$$ the requirements of zero variations with attention to $u$ and $G$ give the next complement equations:
$$\begin{cases}
1 + \lambda g + 2\mu u = 0\\
-\lambda u' + \mu = 0.\\
\end{cases}\qquad(4)$$
$(4.2)$ and $(3.1)$ leads to
$$u= - \frac\mu\lambda(h-y),\qquad(5)$$
then from $(4.1)$
$$g(y) = \dfrac{2\mu^2}{\lambda^2}(h-y) - \dfrac1\lambda.$$
At this point it turns out that condition $(3.2)$ leads to
$$\dfrac1\lambda = 0,\quad u=0,\quad g=0,$$
so it contradicts with $(2.4).$
At that time, the physical representation of the task means that $h$ isn't a constant and must be used as optimization variable in the form
$$\dfrac{\partial \mathcal L}{\partial h} = 0,$$
or
$$u(h) + \lambda u(h)g(h) + \mu (u^2(h)+G(h)) = 0,$$ 
$$G(h) = 0.$$ 
Then it's easy to get
$$G(y) = \dfrac{\mu^2}{\lambda^2}(2hy - y^2) - \dfrac y\lambda,$$
$$\dfrac{\mu^2}{\lambda^2}h^2 - \dfrac h\lambda = 0,$$
and that gives 
$$\begin{cases}
h = \dfrac\lambda{\mu^2}\\
g(y) = \dfrac{\mu^2}{\lambda^2}(h-2y)\\[4pt]
G(y) = \dfrac{\mu^2}{\lambda^2}y(h-y).
\end{cases}\qquad(6)$$
Then
$$B_0 = -\,\dfrac{\mu^3}{\lambda^3}\int\limits_0^h(h-y)(h-2y)\,dy = -\,\dfrac{\mu^3}{\lambda^3}\left(h^2y - \dfrac32hy^2 + \dfrac23y^3\right)\biggr|_0^h = -\,\dfrac{\mu^3h^3}{6\lambda^3},$$
$$M_0 = \,\dfrac{\mu^2}{\lambda^2}\int\limits_0^h\left((h-y)^2 + y(h-y)\right)\,dy = \,\dfrac{\mu^2}{\lambda^2}h\int\limits_0^h(h-y)\,dy  = \,\dfrac{\mu^2h^2}{2\lambda^2}.$$
This leads to
$$\dfrac\mu\lambda h = -\,\sqrt[3]{6B_0\,}$$
$$\,\sqrt[3]{6B_0\,} = \,\sqrt{2M_0\,},$$
$$Q = -\,\dfrac\mu\lambda\int\limits_0^h (h-y)\,dy = -\,\dfrac\mu\lambda\dfrac{ h^2}2 = \sqrt[3]{6B_0}\dfrac h2.$$
And the source data is not sufficient for an unambiguous answer
