
ddthUs =
      - fac::div(phis, Us)
      //- fac::grad(pb*h/(2.*friction->rho())))
      + gs*h;

gradp = fac::interpolate(pb*h/(2.*friction->rho()));


forAll(owner, edgei)
{
    label P = owner[edgei];
    label N = neighbour[edgei];

    if ((h[P] > hminflux.value() && gradp[edgei] < 0) ||
        (h[N] > hminflux.value() && gradp[edgei] > 0))
    {
        ddthUs[P] -= gradp[edgei]/S[N]*Le[edgei];
        ddthUs[N] += gradp[edgei]/S[P]*Le[edgei];
    }
}

ddthUs -= n*(n & ddthUs);

ddth =
      -fac::div(phis);

tau = -friction->tauSp()*Us;
