Given $\Omega:=[x_a,x_b]\times[y_a,y_b]\subset\mathbb{R}^2$, consider the problem $$ \max_{f\in \mathcal{C}^2(\Omega, \mathbb{R})} \iint_\Omega \left[f_y(x,y)\int_{x_a}^{x}f(z,y)\,dz\right]\, dx\, dy $$ subject to a constraint of the form $$ \iint_\Omega \left[\int_{x_a}^{x}f_y(z,y)\, dz -\int_{x_a}^{x_b} z\,f_y(z,y)\, dz\right]^2\, dx \, dy= c \in\mathbb{R} $$ where the subscript $y$ denotes the partial derivative w.r.t. the second variable, i.e., $y$.
I thought the first step would have been to write down the Euler-Lagrange equations, i.e., given the Langragian $\mathcal{L}$, $$ \frac{\partial}{\partial y}\left(\frac{\partial \mathcal{L}}{\partial f_y}\right) = \frac{\partial \mathcal{L}}{\partial f} $$
but the Lagrangian is pretty involved and I got stuck. Do you have any hint?
You can re-write: $$ F(x,y)=\int_{x_a}^xf(z,y)dz $$ It is easy to check that $F_x(x,y)=f(x,y)$. Your problem becomes:
$$ \int \int_{\Omega} F_{xy}(x,y)F(x,y)dxdy $$
Dealing with the constrain, it seems a lot more difficult. You can use differentiation under the integral sign to show that $$ \int_{x_a}^x f_y(z,y)dz=F_y(x,y) $$ As for the second term, you can use the integration by parts formula, which is (for $F'=f$): $$ \int_a^b xf(x)dx=xF(x)|_a^b-\int_a^b F(x)dx $$ Now we can get rid of that "annoying" $z$ in front of $f_y$: $$ \int_{x_a}^{x_b}zf_y(z,y)dz=\int_{x_a}^{x_b}zF_{zy}(z,y)dz=zF_y(z,y)|_{x_a}^{x_b}-\int_{x_a}^{x_b} F_y(z,y)dz $$ This seems all pretty good by now and your constraint becomes: $$ \int \int_{\Omega}(F_y(x,y)-zF_y(z,y)|_{x_a}^{x_b}+\int_{x_a}^{x_b} F_y(z,y)dz)^2dxdy=c $$
It might be possible to simplify the constrain even more by using Fubini, since all of the terms seem to involve a $\partial_y$, but that is a really lengthy calculation.
Now use that your function $f \in C^2$, which means that $F \in C^2$. It also makes exchanging the order of derivatives possible. You can now apply higher-order Euler-Lagrange equations to find the suitable PDE. See the Wikipedia page for a full formula (https://en.wikipedia.org/wiki/Euler–Lagrange_equation).
In your case, your Euler-Lagrange equation becomes (I used the different notations for a reason): $$ F_{xy}-\partial_{xy}F=0 $$ which is a Null-Lagrangian. That means it holds trivially for a all smooth functions. Now you just have to deal with the constraint. Note, however, you have only found a a candidate for a critical point so far.