ClearAll[x, t]
amp = 5787/1000;
tmax = 1000;
L = 20;
vdh[x_?NumericQ, t_?NumericQ] := 1/(1 + h[x, t]) - h[x, t] + Log[h[x, t]/(1 + h[x, t])];
[Lambda][t_Real] := -(1/L)*NIntegrate[vdh[x, t], {x, 0, L}];

myfunel = First[h /. NDSolve[{
 D[h[x, t], t] == D[h[x, t], {x, 2}] - vdh[x, t] - [Lambda][t],
 h[0, t] == h[L, t],
 Derivative[1, 0][h][0, t] == Derivative[1, 0][h][L, t],
 Derivative[2, 0][h][0, t] == Derivative[2, 0][h][L, t],
 h[x, 0] == 1 + 1/Sqrt[2*[Pi]]*Exp[-((x - 10)^2/2)],
 WhenEvent[h[L/2, t] >= amp, "StopIntegration"]},
 h, {x, 0, L}, {t, 0, tmax},
 Method -> {"MethodOfLines",
 "SpatialDiscretization" -> {"TensorProductGrid",
 "MinPoints" -> 101, "MaxPoints" -> 101,
 "DifferenceOrder" -> 4}}, AccuracyGoal -> 20,
 WorkingPrecision -> MachinePrecision, StepMonitor :> Print[t]]]

f[t_] := NIntegrate[h[x, t]^2 + 1/2*(D[h[x, t], x])^2, {x, 0, L}]
Plot[Evaluate[f[t]/.myfunel], {t, 0, tstop}, PlotRange -> All]

vdh = {x, t} [Function] 1/(1 + h[x, t]) - h[x, t] + Log[h[x, t]/(1 + h[x, t])];
λ = t [Function] -(1/L) int@t;
eqn = D[h[x, t], t] == D[h[x, t], {x, 2}] - vdh[x, t] - λ[t];

ic = h[x, 0] == 1 + 1/Sqrt[2*π]*Exp[-((x - 10)^2/2)];

amp = 5787/1000; tmax = 1000; L = 20;

points = 25; difforder = 4; domain = {0, L};

{nodes, weights} = Most[NIntegrate`GaussRuleData[points, MachinePrecision]];

midgrid = Rescale[nodes, {0, 1}, domain];

intrule = int@t -> -Subtract @@ domain weights.Map[Function[x, vdh[x, t]], midgrid] /.
   h[x_, t] :> h[x][t];

grid = Flatten[{0, midgrid, L}];

(*Definition of pdetoode isn't included in this post,
  please find it in the link above.*)
fdd[a__] := NDSolve`FiniteDifferenceDerivative[a, PeriodicInterpolation -> True];
ptoofunc = pdetoode[h[x, t], t, grid, difforder];

ode = ptoofunc[eqn] /. intrule;
odeic = ptoofunc@ic;
hmid = h@grid[[Length@grid/2 // Round]];

sollst = NDSolveValue[{ode, odeic,
    WhenEvent[hmid[t] >= amp // Evaluate, "StopIntegration"]}, h /@ grid, {t, 0, tmax}];
{{tmin, trealmax}} = sollst[[1]]["Domain"];
sol = ListInterpolation[
   Developer`ToPackedArray@#["ValuesOnGrid"] & /@ sollst, {grid,
    sollst[[1]]["Coordinates"][[1]]}];
Plot3D[sol[x, t], {x, 0, L}, {t, tmin, trealmax}, PlotRange -> All]