    Info << "Reading field h" << endl;
    areaScalarField h
    (
        IOobject
        (
            "h",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        aMesh
    );


    Info << "Reading field hUs" << endl;
    areaVectorField hUs
    (
        IOobject
        (
            "hUs",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        aMesh
    );

    areaVectorField Us
    (
        IOobject
        (
            "Us",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        hUs/(h+dimensionedScalar("h_min", dimensionSet(0, 1, 0, 0, 0, 0, 0), 1e-6))
    );

    edgeScalarField phis
    (
        IOobject
        (
            "phis",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        fac::interpolate(hUs) & aMesh.Le()
    );

    edgeVectorField phi2s
    (
        IOobject
        (
            "phi2s",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        fac::interpolate(hUs*Us) & aMesh.Le()
    );

    edgeScalarField gradp
    (
        IOobject
        (
            "gradp",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        aMesh,
        dimensionedScalar("zero", dimensionSet(0, 3, -2, 0, 0, 0, 0), 0)
    );


    areaVectorField ddthUs
    (
        IOobject
        (
            "ddthUs",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        aMesh,
        hUs.dimensions()/dimTime
    );

    areaScalarField pb
    (
        IOobject
        (
            "pb",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        aMesh,
        dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0, 0, 0), 0)
    );

    areaVectorField pbn
    (
        IOobject
        (
            "pbn",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        aMesh,
        dimensionedVector("zero", dimensionSet(1, -1, -2, 0, 0, 0, 0), vector(0, 0, 0))
    );

    areaVectorField tau
    (
        IOobject
        (
            "tau",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        aMesh,
        hUs.dimensions()/dimTime
    );

    areaScalarField ddth
    (
        IOobject
        (
            "ddth",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        aMesh,
        h.dimensions()/dimTime
    );

    const areaVectorField& n = aMesh.faceAreaNormals();

    const areaScalarField gn = g & n;
    const areaVectorField gs = g - gn*n;

    const unallocLabelList& owner = aMesh.owner();
    const unallocLabelList& neighbour = aMesh.neighbour();
    const edgeVectorField& Le = aMesh.Le();
    const scalarField& S = aMesh.S();
