An integral inequality (one variable)
This one is quite fun. Allow me to share some reasoning (not an answer yet!), which I believe might be helpful.
First, let's rewrite the ODE inequality in equivalent ways:
$$r\left(g''+\frac{g'}{x}\right)\geq g'^2$$
$$g''- \frac{g'^2}{r}+\frac{g'}{x} \geq 0$$
$$g''e^{-g/r}- \frac{g'^2}{r} e^{-g/r}+\frac{g'}{x} e^{-g/r} \geq 0 $$
Introducing a new function:
$$f(x)=e^{-g(x)/r}$$
We obtain an inequality:
$$x f''+f' \leq 0$$
Or (provided $f' > 0$ for $x \in (0,1)$):
$$( \ln f' )' \leq - \frac{1}{x}$$
Not sure how this might be useful yet, but we'll see.
Now we play with (*). Our inequality becomes:
$$\left( \int_0^1 f^r xdx \right) \left( \int_0^1 f^{-r} xdx \right) \leq \frac{1}{4(1-r)}$$
Changing the dummy variable in the second integral we have equivalent expression:
$$ \int_0^1 \int_0^1 f(x)^r f(y)^{-r} xy~dxdy \leq \frac{1}{4(1-r)}$$
Introducing:
$$u(x,y)=\ln f(x) - \ln f(y)$$
We have:
$$ \int_0^1 \int_0^1 e^{u(x,y)r} xy~dxdy \leq \frac{1}{4(1-r)}$$
Expand both sides as Taylor series in $r$ (since $ r \in (0,1)$ the geometric series converges):
$$\int_0^1 \int_0^1 \left(1+ur+\frac{u^2r^2}{2}+\dots+\frac{u^nr^n}{n!}+\dots \right) xy~dxdy \leq \frac{1}{4} \left(1+r+r^2+\dots +r^n+\dots \right)$$
While not necessary, it is sufficient if each term on the left is $\leq$ than each term on the right for respective powers of $r$. This might be used as proof, in case it turns out to follow from the ODE inequality. (This part is now obsolete, since this condition is indeed too strong, see Update 2 for further details).
First terms are equal:
$$\int_0^1 \int_0^1 x y dx dy=\frac{1}{4}$$
For the other terms we are supposed to somehow use the ODE inequality to show:
$$I_n=\int_0^1 \int_0^1 \left(\ln f(x)-\ln f(y) \right)^n x y ~dx dy \leq \frac{n!}{4} \tag{1}$$
Funny thing is, $r$ kind of disappeared from the problem statement (provided we can prove the inequality by term in the general case).
Update
Let's consider (1) when $n$ is odd. It is quite clear that:
$$I_{2k+1}=\int_0^1 \int_0^1 \left(\ln f(x)-\ln f(y) \right)^{2k+1} x y ~dx dy =0, \qquad k \in \mathbb{N}$$
Proof. Let's introduce a function $\ln f(x)=h(x)$:
$$I_{2k+1}=\int_0^1 \int_0^1 \left(h(x)-h(y) \right) \left(h(x)-h(y) \right)^{2k} x y ~dx dy =$$
$$=\int_0^1 \int_0^1 h(x) \left(h(x)-h(y) \right)^{2k} x y ~dx dy-\int_0^1 \int_0^1 h(y) \left(h(x)-h(y) \right)^{2k} x y ~dx dy$$
But both integrals are the same, since we can always change the dummy variable. Thus:
$$I_{2k+1}=I-I=0$$
Which means we only need to prove the even case of (1):
$$I_{2k}=\int_0^1 \int_0^1 \left(\ln f(x)-\ln f(y) \right)^{2k} x y ~dx dy \leq \frac{(2k)!}{4} \tag{2}$$
Which can be rewritten as:
$$\sum_{l=0}^{2k} \frac{(-1)^l}{(2k-l)!l!} \int_0^1 h(x)^{2k-l} x~dx \int_0^1 h(y)^l y ~ dy \leq \frac{1}{4}$$
Update 2
From a comment by Binjiu below, it is apparent that requrements $(1)$ and $(2)$ are too strong, and they are not true for the example from the OP (as could be checked directly for $k=1$).
I am still leaving that part of the post intact, as it shows which way we can't go.
Now the problem could be reformulated in another way (where we include the odd terms from the right hand side):
Prove that if for $f:\left(0,1\right)\rightarrow\mathbb{R}$ the inequality $x f''+f' \leq 0$ holds for all $x \in (0,1)$ then there exists $r \in (0,1)$ such that:
$$\int_0^1 \int_0^1 \left(\ln f(x)-\ln f(y) \right)^{2k} x y ~dx dy \leq \frac{(2k)!}{4} \left(1+r \right) \tag{3}$$
I believe this can be written without $r$ as just:
$$\int_0^1 \int_0^1 \left(\ln f(x)-\ln f(y) \right)^{2k} x y ~dx dy < \frac{(2k)!}{2}$$
Let $f(x)=g(e^{-x})$, then $$r\left(g''(x)+\frac{g'(x)}{x}\right)\ge g'^2(x)\iff rf''\ge f'^2$$ and \begin{align} I&=\left( \int_0^1 e^{-g(x)} xdx \right) \left( \int_0^1 e^{g(x)} xdx \right)\\ &=\left( \int_0^{\infty} e^{-f(t)} e^{-2t}dt \right) \left( \int_0^{\infty} e^{f(t)} e^{-2t}dt \right)\tag{$x=e^{-t}$} \end{align} Let $h(x)=e^{-f(x/2)/r}$, then $$rf''\ge f'^2\iff h''\le0$$ and \begin{align} I&=\frac14\left( \int_0^{\infty} h(s)^r e^{-s}ds \right) \left( \int_0^{\infty} h(s)^{-r} e^{-s}ds \right)\tag{$s=2t$}\\ &=\frac14\int_0^{\infty}\int_0^{\infty} \left(\frac{h(t)}{h(s)}\right)^r e^{-t}e^{-s}dtds \end{align} Now the problem becomes
Let $h:\Bbb R_+\to\Bbb R_+$ be a twice-differentiable and concave function, let $r\in(0,1)$, prove that$$J=\int_0^{\infty}\int_0^{\infty} \left(\frac{h(t)}{h(s)}\right)^r e^{-t}e^{-s}dtds\le\frac1{1-r}$$
Let $s>0$, since $h$ is concave, we have $$\forall t>0,\ h(t)\le h'(s)(t-s)+h(s)$$ As $t\to0$, since $h$ is positive, we have $$h'(s)\le\frac{h(s)}{s}$$ Note that $h$ is positive and concave so is increasing, we have \begin{align} \forall t>0,\ \frac{h(t)}{h(s)}&\le\begin{cases} 1 & \text{if $t\le s$} \\ t/s & \text{if $t>s$} \end{cases}\\ &\le\frac{t}{s}+1 \end{align} Therefore, \begin{align} J&\le\int_0^{\infty}\int_0^{\infty} (s+t)^r s^{-r}e^{-t}e^{-s}dtds\\ &=\int_0^{\infty}\Gamma(r+1,s) s^{-r}ds\tag{1}\\ &=\frac{1}{1-r}\tag{2} \end{align} $(1):$ Incomplete Gamma function, $$\int_0^{\infty} (s+t)^r e^{-(s+t)}dt=\int_s^{\infty} u^r e^{-u}du\tag{$u=s+t$}$$ $(2):$ Integration by parts, \begin{align} \int_0^{\infty}\Gamma(r+1,s) s^{-r}ds&=\left[\frac{\Gamma(r+1,s) s^{1-r}}{1-r}\right]_0^{\infty}+\int_0^{\infty}s^re^{-s}\frac{s^{1-r}}{1-r}ds\\ &=0+\frac{\Gamma(2)}{1-r} \end{align}
Update:
As mentioned by @Binjiu, we can give a sharper bound by using directly: $$\forall s,t>0,\ \frac{h(t)}{h(s)}\le \begin{cases} 1 & \text{if $t\le s$} \\ t/s & \text{if $t>s$} \end{cases}$$ This yields \begin{align} J&\le\int_0^{\infty}\int_0^{s} e^{-t}e^{-s}dtds+\int_0^{\infty}\int_s^{\infty} t^r s^{-r}e^{-t}e^{-s}dtds\\ &=\frac{1}{2}+\int_0^{\infty}\gamma(1-r,t)t^{r}e^{-t}dt\tag{3}\\ &=1+r\sum_{n=0}^{\infty}\frac{(-1)^n}{n+1-r}\\ &=1+\frac{r}{2}\left(\psi\left(1-\frac r2\right)-\psi\left(\frac {1-r}2\right)\right)\tag{4} \end{align} $(3):$ Lower incomplete Gamma function, $$\gamma (1-r,t)=\int_0^ts^{-r}e^{-s}ds=-r\gamma (-r,t)-t^{-r}e^{-t}$$ $$\gamma (-r,t)=t^{-r}\sum_{n=0}^{\infty}\frac{(-1)^{n}t^{n}}{n!(n-r)}$$ $(4):$ Digamma function, $$\psi\left(\frac {z+1}2\right)-\psi\left(\frac {z}2\right)=2\sum_{n=0}^{\infty}\frac{(-1)^n}{n+z}$$ Below is the comparison of these two bounds. The maximum difference is $\ln2$ as $r\to1$.