How to properly combine the diffuse and specular terms? As far as I understand, in a BRDF the Fresnel term is telling us the probability for a photon to be reflected or refracted when it hits a surface.

The reflected photons will contribute to the specular term, while the refracted ones will contribute to the diffuse term. Therefore when determining, in a physically based manner, the contribution of a light to the color of material, I feel tempted to just write:

// Assuming for example:
//   diffuse = dot(L, N);
//   specular = pow(dot(H, N), alpha) * (alpha + 2.0) / 8.0;
//   fresnel = f0 + (1.0 - f0) * pow(1.0 - dot(E, H), 5.0);
color = lightIntensity * Lerp(diffuse * albedo, specular, fresnel);


Yet, I don't think I've ever seen it written this way. I have seen the specular term being weighted according to the Fresnel term, but not the diffuse one. In his largely referenced article on PBR, Sébastien Lagarde even states that using $(1 - F)$ to weight the diffuse term is incorrect.

What am I missing?
I would much welcome an explanation that highlights in an obvious manner why this would be incorrect.
Using two Fresnel terms is correct in the sense that any given diffuse path will pass through the surface twice. If you're solving diffusion by tracing a path through the medium until it bounces out again then that you will get two (or more) Fresnel terms for that path as it interacts with the surface.

However, that's not what you're doing with a diffuse BRDF. A diffuse BRDF is intended to represent the average of all those possible diffusion paths. In the case of an Lambertian that average is modeled as uniform reflection and a single albedo value measuring the internal energy loss during diffusion, but more complicated models are possible. Crucially: a diffuse BRDF will already include the aggregate effect of some paths being reflected back into the medium to diffuse further and some passing out immediately. $1 - F_{out}$ is already "baked in" to the BRDF¹ and you do not need factor it in again.

What the Lambertian term does not include is the portion of energy that is lost by being reflected before light enters the medium of diffusion. This is view-dependent, but depends on the precise glossy lobe above it. There is no energy loss at a (non-metallic) surface interface so everything that isn't reflected will be refracted, meaning that what you actually want is to integrate the total energy loss at the surface over all outgoing directions, i.e. $1 - \int \texttt{glossy_bsdf}(\texttt{in}, \texttt{out}) \: \mathrm{d} \texttt{out}$.

It's possible to pre-compute approximations to that integral for specific BRDFs. The end result will in general depend on view direction, material roughness and IOR at least. As a first approximation you can assume that the glossy lobe is a perfectly specular reflector. That gives a weighting of $1 - \int \texttt{glossy} \: \mathrm{d}\texttt{out} = 1 - F_{in}$, which is exactly what you first suggested.

Additionally, note that the Lambertian BRDF is the albedo divided by $\pi$ and that the cosine term is a measure of how attenuated the incoming light is on the surface; it applies to both glossy and diffuse reflectance.

So, roughly:

// Assuming for example:
//   diffuse = albedo / PI;
//   specular = my_favorite_glossy_brdf(in_dir, out_dir, roughness);
//   fresnel = f0 + (1.0 - f0) * pow(1.0 - dot(E, H), 5.0);
//   total_surface_reflection = fresnel
color = lightIntensity * dot(L, N) * Lerp(diffuse, specular, total_surface_reflection);


¹ Really the integral of $F$ over all possible incident internal diffusion exit paths that result in your out direction, but I digress.