{
    // Solve the enthalpy equation in total enthalpy formulation (see K)

    T.storePrevIter();
    
    K = 0.5*(magSqr(U));

    fvScalarMatrix hEqn
    (
        fvm::div(phi, h)
      + fvm::SuSp(-fvc::div(phi), h)
      + fvc::div(phi, K)
      - fvm::laplacian(turbulence->alphaEff(), h)
     ==
        // Viscous heating: note sign (devRhoReff has a minus in it)
      - (turbulence->devRhoReff() && fvc::grad(U))
    );

    hEqn.relax();
    hEqn.solve();

    // Bounding of enthalpy taken out
    thermo.correct();
    psis = thermo.psi()/thermo.Cp()*thermo.Cv();
}