From 973be0cf7a9e0e13ec2f0acc6f3b2c176cc99ea6 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Tue, 7 Sep 2010 09:43:03 +0100
Subject: [PATCH] ENH: Merge recent surface film developments into main line

---
 src/surfaceFilmModels/Make/files              |   13 +
 src/surfaceFilmModels/Make/options            |   20 +-
 ...rectMappedFixedInternalValueFvPatchField.C |    6 +-
 .../directMappedNamedFixedValueFvPatchField.C |  334 +++
 .../directMappedNamedFixedValueFvPatchField.H |  164 ++
 ...directMappedNamedFixedValueFvPatchFields.C |   43 +
 ...directMappedNamedFixedValueFvPatchFields.H |   49 +
 ...ectMappedNamedFixedValueFvPatchFieldsFwd.H |   50 +
 ...ectMappedNamesFixedValueFvPatchFieldsFwd.H |   50 +
 .../htcConv/htcConvFvPatchScalarField.C       |  165 ++
 .../htcConv/htcConvFvPatchScalarField.H       |  159 ++
 ...alphatFilmWallFunctionFvPatchScalarField.C |  247 ++
 ...alphatFilmWallFunctionFvPatchScalarField.H |  172 ++
 .../mutFilmWallFunctionFvPatchScalarField.C   |  246 ++
 .../mutFilmWallFunctionFvPatchScalarField.H   |  172 ++
 .../cloudInjection/cloudInjection.C           |  109 +
 .../cloudInjection/cloudInjection.H           |  122 +
 .../injectionModel/injectionModel.C           |    3 +-
 .../injectionModel/injectionModel.H           |    2 +-
 .../injectionModel/injectionModelNew.C        |    2 +-
 .../injectionModel/noInjection/noInjection.C  |    4 +-
 .../removeInjection/removeInjection.C         |    4 +-
 .../removeInjection/removeInjection.H         |    2 +-
 .../constantHeatTransfer.C                    |  105 +
 .../constantHeatTransfer.H                    |  114 +
 .../heatTransferModel/heatTransferModel.C     |   69 +
 .../heatTransferModel/heatTransferModel.H     |  159 ++
 .../heatTransferModel/heatTransferModelI.H    |   44 +
 .../heatTransferModel/heatTransferModelNew.C  |   59 +
 .../mappedConvectiveHeatTransfer.C            |  121 +
 .../mappedConvectiveHeatTransfer.H            |  120 +
 .../noPhaseChange/noPhaseChange.C             |   10 +-
 .../noPhaseChange/noPhaseChange.H             |    2 +-
 .../phaseChangeModel/phaseChangeModel.C       |   15 +-
 .../phaseChangeModel/phaseChangeModel.H       |    8 +-
 .../phaseChangeModel/phaseChangeModelNew.C    |    2 +-
 .../standardPhaseChange/standardPhaseChange.C |  185 ++
 .../standardPhaseChange/standardPhaseChange.H |  121 +
 .../kinematicSingleLayer.C                    | 2364 +++++++++--------
 .../kinematicSingleLayer.H                    | 1040 ++++----
 .../kinematicSingleLayerI.H                   |   35 +
 .../kinematicSingleLayerTemplates.C           |  124 +-
 .../surfaceFilmModel/noFilm/noFilm.C          |   76 +-
 .../surfaceFilmModel/noFilm/noFilm.H          |   25 +-
 .../surfaceFilmModel/surfaceFilmModel.C       |   65 +-
 .../surfaceFilmModel/surfaceFilmModel.H       |   74 +-
 .../surfaceFilmModel/surfaceFilmModelI.H      |    7 +
 .../surfaceFilmModel/surfaceFilmModelNew.C    |   16 +-
 .../thermoSingleLayer/thermoSingleLayer.C     |  463 +++-
 .../thermoSingleLayer/thermoSingleLayer.H     |  167 +-
 .../thermoSingleLayer/thermoSingleLayerI.H    |   59 +-
 51 files changed, 5976 insertions(+), 1811 deletions(-)
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchField.C
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchField.H
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.C
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.H
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFieldsFwd.H
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamesFixedValueFvPatchFieldsFwd.H
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.C
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.H
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.H
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.C
 create mode 100644 src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.H
 create mode 100644 src/surfaceFilmModels/submodels/kinematic/injectionModel/cloudInjection/cloudInjection.C
 create mode 100644 src/surfaceFilmModels/submodels/kinematic/injectionModel/cloudInjection/cloudInjection.H
 create mode 100644 src/surfaceFilmModels/submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.C
 create mode 100644 src/surfaceFilmModels/submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.H
 create mode 100644 src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.C
 create mode 100644 src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.H
 create mode 100644 src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelI.H
 create mode 100644 src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelNew.C
 create mode 100644 src/surfaceFilmModels/submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.C
 create mode 100644 src/surfaceFilmModels/submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.H
 create mode 100644 src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
 create mode 100644 src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H

diff --git a/src/surfaceFilmModels/Make/files b/src/surfaceFilmModels/Make/files
index cbaf49b8100..8b84a94e146 100644
--- a/src/surfaceFilmModels/Make/files
+++ b/src/surfaceFilmModels/Make/files
@@ -10,17 +10,30 @@ surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
 submodels/kinematic/injectionModel/injectionModel/injectionModel.C
 submodels/kinematic/injectionModel/injectionModel/injectionModelNew.C
 submodels/kinematic/injectionModel/noInjection/noInjection.C
+submodels/kinematic/injectionModel/cloudInjection/cloudInjection.C
 submodels/kinematic/injectionModel/removeInjection/removeInjection.C
 
 submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.C
 submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModelNew.C
 submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C
+submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
+
+submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.C
+submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelNew.C
+submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.C
+submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.C
 
 
 /* Boundary conditions */
 derivedFvPatchFields/directMappedFixedInternalValue/directMappedFixedInternalValueFvPatchFields.C
+derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.C
 derivedFvPatchFields/directMappedFixedPushedInternalValue/directMappedFixedPushedInternalValueFvPatchFields.C
 derivedFvPatchFields/filmHeightInletVelocity/filmHeightInletVelocityFvPatchVectorField.C
+derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.C
+
+/* Wall functions for primary region */
+derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C
+derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.C
 
 
 LIB = $(FOAM_LIBBIN)/libsurfaceFilmModels
diff --git a/src/surfaceFilmModels/Make/options b/src/surfaceFilmModels/Make/options
index d27c95d033d..bd04197ec3c 100644
--- a/src/surfaceFilmModels/Make/options
+++ b/src/surfaceFilmModels/Make/options
@@ -1,7 +1,23 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude
 
 EXE_LIBS = \
+    -lSLGThermo \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lpdf \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/directMappedFixedInternalValue/directMappedFixedInternalValueFvPatchField.C b/src/surfaceFilmModels/derivedFvPatchFields/directMappedFixedInternalValue/directMappedFixedInternalValueFvPatchField.C
index 0f5c3db7c10..9c84101d853 100644
--- a/src/surfaceFilmModels/derivedFvPatchFields/directMappedFixedInternalValue/directMappedFixedInternalValueFvPatchField.C
+++ b/src/surfaceFilmModels/derivedFvPatchFields/directMappedFixedInternalValue/directMappedFixedInternalValueFvPatchField.C
@@ -111,10 +111,8 @@ void directMappedFixedInternalValueFvPatchField<Type>::updateCoeffs()
     directMappedFixedValueFvPatchField<Type>::updateCoeffs();
 
     // Get the coupling information from the directMappedPatchBase
-    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
-    (
-        this->patch().patch()
-    );
+    const directMappedPatchBase& mpp =
+        refCast<const directMappedPatchBase>(this->patch().patch());
     const polyMesh& nbrMesh = mpp.sampleMesh();
     const fvPatch& nbrPatch =
         refCast<const fvMesh>
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchField.C b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchField.C
new file mode 100644
index 00000000000..0be8a36bd1c
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchField.C
@@ -0,0 +1,334 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "directMappedNamedFixedValueFvPatchField.H"
+#include "directMappedPatchBase.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+directMappedNamedFixedValueFvPatchField<Type>::
+directMappedNamedFixedValueFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<Type>(p, iF),
+    fieldName_(iF.name()),
+    setAverage_(false),
+    average_(pTraits<Type>::zero)
+{}
+
+
+template<class Type>
+directMappedNamedFixedValueFvPatchField<Type>::
+directMappedNamedFixedValueFvPatchField
+(
+    const directMappedNamedFixedValueFvPatchField<Type>& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
+    fieldName_(ptf.fieldName_),
+    setAverage_(ptf.setAverage_),
+    average_(ptf.average_)
+{
+    if (!isA<directMappedPatchBase>(this->patch().patch()))
+    {
+        FatalErrorIn
+        (
+            "directMappedNamedFixedValueFvPatchField<Type>::"
+            "directMappedNamedFixedValueFvPatchField\n"
+            "(\n"
+            "    const directMappedNamedFixedValueFvPatchField<Type>&,\n"
+            "    const fvPatch&,\n"
+            "    const Field<Type>&,\n"
+            "    const fvPatchFieldMapper&\n"
+            ")\n"
+        )   << "\n    patch type '" << p.type()
+            << "' not type '" << directMappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+}
+
+
+template<class Type>
+directMappedNamedFixedValueFvPatchField<Type>::
+directMappedNamedFixedValueFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<Type>(p, iF, dict),
+    fieldName_(dict.lookupOrDefault<word>("fieldName", iF.name())),
+    setAverage_(readBool(dict.lookup("setAverage"))),
+    average_(pTraits<Type>(dict.lookup("average")))
+{
+    if (!isA<directMappedPatchBase>(this->patch().patch()))
+    {
+        FatalErrorIn
+        (
+            "directMappedNamedFixedValueFvPatchField<Type>::"
+            "directMappedNamedFixedValueFvPatchField"
+            "("
+            "    const fvPatch&, "
+            "    const DimensionedField<Type, volMesh>& iF, "
+            "    const dictionary&"
+            ")"
+        )   << "\n    patch type '" << p.type()
+            << "' not type '" << directMappedPatchBase::typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << fieldName_
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
+    //// Force calculation of schedule (uses parallel comms)
+    //const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
+    //(
+    //    this->patch().patch()
+    //);
+    //(void)mpp.map().schedule();
+}
+
+
+template<class Type>
+directMappedNamedFixedValueFvPatchField<Type>::
+directMappedNamedFixedValueFvPatchField
+(
+    const directMappedNamedFixedValueFvPatchField<Type>& ptf
+)
+:
+    fixedValueFvPatchField<Type>(ptf),
+    fieldName_(ptf.fieldName_),
+    setAverage_(ptf.setAverage_),
+    average_(ptf.average_)
+{}
+
+
+template<class Type>
+directMappedNamedFixedValueFvPatchField<Type>::
+directMappedNamedFixedValueFvPatchField
+(
+    const directMappedNamedFixedValueFvPatchField<Type>& ptf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<Type>(ptf, iF),
+    fieldName_(ptf.fieldName_),
+    setAverage_(ptf.setAverage_),
+    average_(ptf.average_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void directMappedNamedFixedValueFvPatchField<Type>::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    // Get the scheduling information from the directMappedPatchBase
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
+    (
+        directMappedNamedFixedValueFvPatchField<Type>::patch().patch()
+    );
+    const mapDistribute& distMap = mpp.map();
+
+    // Force recalculation of schedule
+    (void)distMap.schedule();
+
+    const fvMesh& nbrMesh = refCast<const fvMesh>(mpp.sampleMesh());
+
+    // Result of obtaining remote values
+    Field<Type> newValues;
+
+    switch (mpp.mode())
+    {
+        case directMappedPatchBase::NEARESTCELL:
+        {
+            if (mpp.sameRegion())
+            {
+                newValues = this->internalField();
+            }
+            else
+            {
+                newValues = nbrMesh.lookupObject<fieldType>
+                (
+                    fieldName_
+                ).internalField();
+            }
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                newValues
+            );
+
+            break;
+        }
+        case directMappedPatchBase::NEARESTPATCHFACE:
+        {
+            const label nbrPatchID =
+                nbrMesh.boundaryMesh().findPatchID(mpp.samplePatch());
+            if (nbrPatchID < 0)
+            {
+                FatalErrorIn
+                (
+                    "void directMappedNamedFixedValueFvPatchField<Type>::"
+                    "updateCoeffs()"
+                )<< "Unable to find sample patch " << mpp.samplePatch()
+                 << " in region " << mpp.sampleRegion()
+                 << " for patch " << this->patch().name() << nl
+                 << abort(FatalError);
+            }
+
+            const fieldType& nbrField =
+                nbrMesh.lookupObject<fieldType>(fieldName_);
+            newValues = nbrField.boundaryField()[nbrPatchID];
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                newValues
+            );
+
+            break;
+        }
+        case directMappedPatchBase::NEARESTFACE:
+        {
+            Field<Type> allValues(nbrMesh.nFaces(), pTraits<Type>::zero);
+
+            const fieldType& nbrField =
+                nbrMesh.lookupObject<fieldType>(fieldName_);
+            forAll(nbrField.boundaryField(), patchI)
+            {
+                const fvPatchField<Type>& pf =
+                    nbrField.boundaryField()[patchI];
+                label faceStart = pf.patch().patch().start();
+
+                forAll(pf, faceI)
+                {
+                    allValues[faceStart++] = pf[faceI];
+                }
+            }
+
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                allValues
+            );
+
+            newValues = this->patch().patchSlice(allValues);
+
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "directMappedNamedFixedValueFvPatchField<Type>::updateCoeffs()"
+            )<< "Unknown sampling mode: " << mpp.mode()
+             << nl << abort(FatalError);
+        }
+    }
+
+    if (setAverage_)
+    {
+        Type averagePsi =
+            gSum(this->patch().magSf()*newValues)
+           /gSum(this->patch().magSf());
+
+        if (mag(averagePsi)/mag(average_) > 0.5)
+        {
+            newValues *= mag(average_)/mag(averagePsi);
+        }
+        else
+        {
+            newValues += (average_ - averagePsi);
+        }
+    }
+
+    this->operator==(newValues);
+
+    if (debug)
+    {
+        Info<< "directMapped on field:" << fieldName_
+            << " patch:" << this->patch().name()
+            << "  avg:" << gAverage(*this)
+            << "  min:" << gMin(*this)
+            << "  max:" << gMax(*this)
+            << endl;
+    }
+
+    fixedValueFvPatchField<Type>::updateCoeffs();
+}
+
+
+template<class Type>
+void directMappedNamedFixedValueFvPatchField<Type>::write(Ostream& os) const
+{
+    fvPatchField<Type>::write(os);
+    os.writeKeyword("fieldName") << fieldName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl;
+    os.writeKeyword("average") << average_ << token::END_STATEMENT << nl;
+    this->writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchField.H b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchField.H
new file mode 100644
index 00000000000..711ed1ec98e
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchField.H
@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::directMappedNamedFixedValueFvPatchField
+
+Description
+    Variant of directMappedFixedValueFvPatch where the name of the field to
+    map is input.
+
+SourceFiles
+    directMappedNamedFixedValueFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedNamedFixedValueFvPatchField_H
+#define directMappedNamedFixedValueFvPatchField_H
+
+#include "fixedValueFvPatchField.H"
+#include "directMappedFvPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+          Class directMappedNamesFixedValueFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class directMappedNamedFixedValueFvPatchField
+:
+    public fixedValueFvPatchField<Type>
+{
+    // Private data
+
+        //- Name of field to sample - defaults to field associated with this
+        //  patch if not specified
+        word fieldName_;
+
+        //- If true adjust the mapped field to maintain average value average_
+        bool setAverage_;
+
+        //- Average value the mapped field is adjusted to maintain if
+        //  setAverage_ is set true
+        Type average_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("directMappedNamedFixedValue");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        directMappedNamedFixedValueFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        directMappedNamedFixedValueFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given directMappedNamedFixedValueFvPatchField
+        //  onto a new patch
+        directMappedNamedFixedValueFvPatchField
+        (
+            const directMappedNamedFixedValueFvPatchField<Type>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        directMappedNamedFixedValueFvPatchField
+        (
+            const directMappedNamedFixedValueFvPatchField<Type>&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type> > clone() const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new directMappedNamedFixedValueFvPatchField<Type>(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        directMappedNamedFixedValueFvPatchField
+        (
+            const directMappedNamedFixedValueFvPatchField<Type>&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type> > clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new directMappedNamedFixedValueFvPatchField<Type>(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "directMappedNamedFixedValueFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.C b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.C
new file mode 100644
index 00000000000..849cc558259
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.C
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "directMappedNamedFixedValueFvPatchFields.H"
+#include "volMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makePatchFields(directMappedNamedFixedValue);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.H b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.H
new file mode 100644
index 00000000000..2f815fb81aa
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFields.H
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedNamedFixedValueFvPatchFields_H
+#define directMappedNamedFixedValueFvPatchFields_H
+
+#include "directMappedNamedFixedValueFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(directMappedNamedFixedValue)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFieldsFwd.H b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFieldsFwd.H
new file mode 100644
index 00000000000..c7d6e5c11cf
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamedFixedValueFvPatchFieldsFwd.H
@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedNamedFixedValueFvPatchFieldsFwd_H
+#define directMappedNamedFixedValueFvPatchFieldsFwd_H
+
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type> class directMappedNamedFixedValueFvPatchField;
+
+makePatchTypeFieldTypedefs(directMappedNamedFixedValue)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamesFixedValueFvPatchFieldsFwd.H b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamesFixedValueFvPatchFieldsFwd.H
new file mode 100644
index 00000000000..c7d6e5c11cf
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/directMappedNamedFixedValue/directMappedNamesFixedValueFvPatchFieldsFwd.H
@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedNamedFixedValueFvPatchFieldsFwd_H
+#define directMappedNamedFixedValueFvPatchFieldsFwd_H
+
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type> class directMappedNamedFixedValueFvPatchField;
+
+makePatchTypeFieldTypedefs(directMappedNamedFixedValue)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.C b/src/surfaceFilmModels/derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.C
new file mode 100644
index 00000000000..c8ca479b6a9
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.C
@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "htcConvFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+htcConvFvPatchScalarField::htcConvFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    L_(1.0)
+{}
+
+
+htcConvFvPatchScalarField::htcConvFvPatchScalarField
+(
+    const htcConvFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    L_(ptf.L_)
+{}
+
+
+htcConvFvPatchScalarField::htcConvFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    L_(readScalar(dict.lookup("L")))
+{}
+
+
+htcConvFvPatchScalarField::htcConvFvPatchScalarField
+(
+    const htcConvFvPatchScalarField& htcpsf
+)
+:
+    fixedValueFvPatchScalarField(htcpsf),
+    L_(htcpsf.L_)
+{}
+
+
+htcConvFvPatchScalarField::htcConvFvPatchScalarField
+(
+    const htcConvFvPatchScalarField& htcpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(htcpsf, iF),
+    L_(htcpsf.L_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void htcConvFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField alphaEffw = rasModel.alphaEff()().boundaryField()[patchI];
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const vectorField& Uc = rasModel.U();
+    const vectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField& Tw = rasModel.thermo().T().boundaryField()[patchI];
+    const scalarField Cpw = rasModel.thermo().Cp(Tw, patchI);
+
+    const scalarField kappaw = Cpw*alphaEffw;
+    const scalarField Pr = muw*Cpw/kappaw;
+
+    scalarField& htc = *this;
+    forAll(htc, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar Re = rhow[faceI]*mag(Uc[faceCellI] - Uw[faceI])*L_/muw[faceI];
+
+        if (Re < 5.0E+05)
+        {
+            htc[faceI] = 0.664*sqrt(Re)*cbrt(Pr[faceI])*kappaw[faceI]/L_;
+        }
+        else
+        {
+            htc[faceI] = 0.037*pow(Re, 0.8)*cbrt(Pr[faceI])*kappaw[faceI]/L_;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void htcConvFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+    os.writeKeyword("L") << L_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    htcConvFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.H b/src/surfaceFilmModels/derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.H
new file mode 100644
index 00000000000..4b5585e578e
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/htcConv/htcConvFvPatchScalarField.H
@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::compressible::RASModels::htcConvFvPatchScalarField
+
+Description
+    Convective heat transfer boundary condition
+
+SourceFiles
+    htcConvFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleMutRoughWallFunctionFvPatchScalarField_H
+#define compressibleMutRoughWallFunctionFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class htcConvFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class htcConvFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+protected:
+
+    // Protected data
+
+        //- L Length scale [m]
+        const scalar L_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("htcConvection");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        htcConvFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        htcConvFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  htcConvFvPatchScalarField
+        //  onto a new patch
+        htcConvFvPatchScalarField
+        (
+            const htcConvFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        htcConvFvPatchScalarField
+        (
+            const htcConvFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new htcConvFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        htcConvFvPatchScalarField
+        (
+            const htcConvFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new htcConvFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C b/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C
new file mode 100644
index 00000000000..3a43ceaa7ae
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C
@@ -0,0 +1,247 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "alphatFilmWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "surfaceFilmModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "directMappedWallPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphatFilmWallFunctionFvPatchScalarField::
+alphatFilmWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    B_(5.5),
+    yPlusCrit_(11.05),
+    Cmu_(0.09),
+    kappa_(0.41),
+    Prt_(0.85)
+{}
+
+
+alphatFilmWallFunctionFvPatchScalarField::
+alphatFilmWallFunctionFvPatchScalarField
+(
+    const alphatFilmWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    B_(ptf.B_),
+    yPlusCrit_(ptf.yPlusCrit_),
+    Cmu_(ptf.Cmu_),
+    kappa_(ptf.kappa_),
+    Prt_(ptf.Prt_)
+{}
+
+
+alphatFilmWallFunctionFvPatchScalarField::
+alphatFilmWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    B_(dict.lookupOrDefault("B", 5.5)),
+    yPlusCrit_(dict.lookupOrDefault("yPlusCrit", 11.05)),
+    Cmu_(dict.lookupOrDefault("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault("kappa", 0.41)),
+    Prt_(dict.lookupOrDefault("Prt", 0.85))
+{}
+
+
+alphatFilmWallFunctionFvPatchScalarField::
+alphatFilmWallFunctionFvPatchScalarField
+(
+    const alphatFilmWallFunctionFvPatchScalarField& fwfpsf
+)
+:
+    fixedValueFvPatchScalarField(fwfpsf),
+    B_(fwfpsf.B_),
+    yPlusCrit_(fwfpsf.yPlusCrit_),
+    Cmu_(fwfpsf.Cmu_),
+    kappa_(fwfpsf.kappa_),
+    Prt_(fwfpsf.Prt_)
+{}
+
+
+alphatFilmWallFunctionFvPatchScalarField::
+alphatFilmWallFunctionFvPatchScalarField
+(
+    const alphatFilmWallFunctionFvPatchScalarField& fwfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(fwfpsf, iF),
+    B_(fwfpsf.B_),
+    yPlusCrit_(fwfpsf.yPlusCrit_),
+    Cmu_(fwfpsf.Cmu_),
+    kappa_(fwfpsf.kappa_),
+    Prt_(fwfpsf.Prt_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    bool ok =
+        db().objectRegistry::foundObject
+        <surfaceFilmModels::surfaceFilmModel>("surfaceFilmProperties");
+
+    if (!ok)
+    {
+        // do nothing on construction - film model doesn't exist yet
+        return;
+    }
+
+    const label patchI = patch().index();
+
+    // Retrieve phase change mass from surface film model
+    const surfaceFilmModels::surfaceFilmModel& filmModel =
+        db().objectRegistry::lookupObject
+            <surfaceFilmModels::surfaceFilmModel>("surfaceFilmProperties");
+
+    const directMappedWallPolyPatch& wpp =
+         refCast<const directMappedWallPolyPatch>(patch().patch());
+    const mapDistribute& distMap = wpp.map();
+    label filmPatchI = -1;
+    forAll(filmModel.primaryPatchIDs(), i)
+    {
+        if (filmModel.primaryPatchIDs()[i] == patchI)
+        {
+            filmPatchI = filmModel.filmBottomPatchIDs()[i];
+            break;
+        }
+    }
+
+    scalarField mDotFilm =
+        filmModel.massPhaseChangeForPrimary().boundaryField()[filmPatchI];
+    distMap.distribute(mDotFilm);
+
+    // Retrieve RAS turbulence model
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+
+    const scalarField& y = rasModel.y()[patchI];
+    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+    const scalarField& alphaw = rasModel.alpha().boundaryField()[patchI];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    // Populate alphat field values
+    scalarField& alphat = *this;
+    forAll(alphat, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar uTau = Cmu25*sqrt(k[faceCellI]);
+
+        scalar yPlus = y[faceI]*uTau/(muw[faceI]/rhow[faceI]);
+
+        scalar Pr = muw[faceI]/alphaw[faceI];
+
+        scalar factor = 0.0;
+        scalar mStar = mDotFilm[faceI]/(y[faceI]*uTau);
+        if (yPlus > yPlusCrit_)
+        {
+            scalar expTerm = exp(min(50.0, yPlusCrit_*mStar*Pr));
+            scalar yPlusRatio = yPlus/yPlusCrit_;
+            scalar powTerm = mStar*Prt_/kappa_;
+            factor =
+                mStar/(expTerm*(pow(yPlusRatio, powTerm)) - 1.0 + ROOTVSMALL);
+        }
+        else
+        {
+            scalar expTerm = exp(min(50.0, yPlus*mStar*Pr));
+            factor = mStar/(expTerm - 1.0 + ROOTVSMALL);
+        }
+
+        scalar dx = patch().deltaCoeffs()[faceI];
+
+        scalar alphaEff = dx*rhow[faceI]*uTau*factor;
+
+        alphat[faceI] = max(alphaEff - alphaw[faceI], 0.0);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphatFilmWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+    os.writeKeyword("B") << B_ << token::END_STATEMENT << nl;
+    os.writeKeyword("yPlusCrit") << yPlusCrit_ << token::END_STATEMENT << nl;
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphatFilmWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.H b/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.H
new file mode 100644
index 00000000000..c899e853da0
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::compressible::RASModels::alphatFilmWallFunctionFvPatchScalarField
+
+Description
+    Turbulent thermal diffusivity boundary conditions for use with surface
+    film models.
+
+SourceFiles
+    alphatFilmWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleMutRoughWallFunctionFvPatchScalarField_H
+#define compressibleMutRoughWallFunctionFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+          Class alphatFilmWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class alphatFilmWallFunctionFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+protected:
+
+    // Protected data
+
+        //- B Coefficient (default = 5.5)
+        scalar B_;
+
+        //- y+ value for laminar -> turbulent transition (default = 11.05)
+        scalar yPlusCrit_;
+
+        //- Turbulent Cmu coefficient (default = 0.09)
+        scalar Cmu_;
+
+        //- Von-Karman constant (default = 0.41)
+        scalar kappa_;
+
+        //- Turbulent Prandtl number (default = 0.85)
+        scalar Prt_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("alphatFilmWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        alphatFilmWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        alphatFilmWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  alphatFilmWallFunctionFvPatchScalarField
+        //  onto a new patch
+        alphatFilmWallFunctionFvPatchScalarField
+        (
+            const alphatFilmWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        alphatFilmWallFunctionFvPatchScalarField
+        (
+            const alphatFilmWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphatFilmWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        alphatFilmWallFunctionFvPatchScalarField
+        (
+            const alphatFilmWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphatFilmWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.C b/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.C
new file mode 100644
index 00000000000..f7a44c586ac
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.C
@@ -0,0 +1,246 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "mutFilmWallFunctionFvPatchScalarField.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "RASModel.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFilmModel.H"
+#include "directMappedWallPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<scalarField> mutFilmWallFunctionFvPatchScalarField::calcUTau
+(
+    const scalarField& magGradU
+) const
+{
+    tmp<scalarField> tuTau(new scalarField(patch().size(), 0.0));
+    scalarField& uTau = tuTau();
+
+    bool ok =
+        db().objectRegistry::foundObject
+        <surfaceFilmModels::surfaceFilmModel>("surfaceFilmProperties");
+
+    if (!ok)
+    {
+        // do nothing on construction - film model doesn't exist yet
+        return tuTau;
+    }
+
+    const label patchI = patch().index();
+
+    // Retrieve phase change mass from surface film model
+    const surfaceFilmModels::surfaceFilmModel& filmModel =
+        db().objectRegistry::lookupObject
+            <surfaceFilmModels::surfaceFilmModel>("surfaceFilmProperties");
+
+    const directMappedWallPolyPatch& wpp =
+         refCast<const directMappedWallPolyPatch>(patch().patch());
+    const mapDistribute& distMap = wpp.map();
+    label filmPatchI = -1;
+    forAll(filmModel.primaryPatchIDs(), i)
+    {
+        if (filmModel.primaryPatchIDs()[i] == patchI)
+        {
+            filmPatchI = filmModel.filmBottomPatchIDs()[i];
+            break;
+        }
+    }
+
+    scalarField mDotFilm =
+        filmModel.massPhaseChangeForPrimary().boundaryField()[filmPatchI];
+    distMap.distribute(mDotFilm);
+
+
+    // Retrieve RAS turbulence model
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const fvPatchScalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    forAll(uTau, faceI)
+    {
+        label faceCellI = patch().faceCells()[faceI];
+
+        scalar ut = Cmu25*sqrt(k[faceCellI]);
+
+        scalar yPlus = y[faceI]*ut/(muw[faceI]/rhow[faceI]);
+
+        scalar mStar = mDotFilm[faceI]/(y[faceI]*ut);
+
+        scalar factor = 0.0;
+        if (yPlus > yPlusCrit_)
+        {
+            scalar expTerm = exp(min(50.0, B_*mStar));
+            scalar powTerm = pow(yPlus, mStar/kappa_);
+            factor = mStar/(expTerm*powTerm - 1.0 + ROOTVSMALL);
+        }
+        else
+        {
+            scalar expTerm = exp(min(50.0, mStar));
+            factor = mStar/(expTerm*yPlus - 1.0 + ROOTVSMALL);
+        }
+
+        uTau[faceI] = sqrt(max(0, magGradU[faceI]*ut*factor));
+    }
+
+    return tuTau;
+}
+
+
+tmp<scalarField> mutFilmWallFunctionFvPatchScalarField::calcMut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField magGradU = mag(Uw.snGrad());
+    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+
+    return max
+    (
+        scalar(0),
+        rhow*sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - muw
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mutFilmWallFunctionFvPatchScalarField::mutFilmWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mutkWallFunctionFvPatchScalarField(p, iF),
+    B_(5.5),
+    yPlusCrit_(11.05)
+{}
+
+
+mutFilmWallFunctionFvPatchScalarField::mutFilmWallFunctionFvPatchScalarField
+(
+    const mutFilmWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mutkWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
+    B_(5.5),
+    yPlusCrit_(11.05)
+{}
+
+
+mutFilmWallFunctionFvPatchScalarField::mutFilmWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mutkWallFunctionFvPatchScalarField(p, iF, dict),
+    B_(dict.lookupOrDefault("B", 5.5)),
+    yPlusCrit_(dict.lookupOrDefault("yPlusCrit", 11.05))
+{}
+
+
+mutFilmWallFunctionFvPatchScalarField::mutFilmWallFunctionFvPatchScalarField
+(
+    const mutFilmWallFunctionFvPatchScalarField& wfpsf
+)
+:
+    mutkWallFunctionFvPatchScalarField(wfpsf),
+    B_(wfpsf.B_),
+    yPlusCrit_(wfpsf.yPlusCrit_)
+{}
+
+
+mutFilmWallFunctionFvPatchScalarField::mutFilmWallFunctionFvPatchScalarField
+(
+    const mutFilmWallFunctionFvPatchScalarField& wfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mutkWallFunctionFvPatchScalarField(wfpsf, iF),
+    B_(wfpsf.B_),
+    yPlusCrit_(wfpsf.yPlusCrit_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<scalarField> mutFilmWallFunctionFvPatchScalarField::yPlus() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+
+    return y*calcUTau(mag(Uw.snGrad()))/(muw/rhow);
+}
+
+
+void mutFilmWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+    writeLocalEntries(os);
+    os.writeKeyword("B") << B_ << token::END_STATEMENT << nl;
+    os.writeKeyword("yPlusCrit") << yPlusCrit_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, mutFilmWallFunctionFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.H b/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.H
new file mode 100644
index 00000000000..e6a40f5f78c
--- /dev/null
+++ b/src/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutFilmWallFunction/mutFilmWallFunctionFvPatchScalarField.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::compressible::RASModels::
+    mutFilmWallFunctionFvPatchScalarField
+
+Description
+    Wall function boundary condition for use with surface film models.
+
+SourceFiles
+    mutFilmWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleMutSpalartAllmarasWallFunctionFvPatchScalarField_H
+#define compressibleMutSpalartAllmarasWallFunctionFvPatchScalarField_H
+
+#include "mutkWallFunctionFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+           Class mutFilmWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class mutFilmWallFunctionFvPatchScalarField
+:
+    public mutkWallFunctionFvPatchScalarField
+{
+protected:
+
+    // Protected data
+
+        //- B Coefficient (default = 5.5)
+        scalar B_;
+
+        //- y+ value for laminar -> turbulent transition (default = 11.05)
+        scalar yPlusCrit_;
+
+
+    // Protected member functions
+
+        //- Calculate the turbulence viscosity
+        virtual tmp<scalarField> calcMut() const;
+
+        //- Calculate the friction velocity
+        virtual tmp<scalarField> calcUTau(const scalarField& magGradU) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("mutFilmWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        mutFilmWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        mutFilmWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  mutFilmWallFunctionFvPatchScalarField
+        //  onto a new patch
+        mutFilmWallFunctionFvPatchScalarField
+        (
+            const mutFilmWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        mutFilmWallFunctionFvPatchScalarField
+        (
+            const mutFilmWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutFilmWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        mutFilmWallFunctionFvPatchScalarField
+        (
+            const mutFilmWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new mutFilmWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Calculate and return the yPlus at the boundary
+            virtual tmp<scalarField> yPlus() const;
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/kinematic/injectionModel/cloudInjection/cloudInjection.C b/src/surfaceFilmModels/submodels/kinematic/injectionModel/cloudInjection/cloudInjection.C
new file mode 100644
index 00000000000..a1993a7e33e
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/kinematic/injectionModel/cloudInjection/cloudInjection.C
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "cloudInjection.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "Time.H"
+#include "mathematicalConstants.H"
+#include "Random.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace surfaceFilmModels
+    {
+        defineTypeNameAndDebug(cloudInjection, 0);
+        addToRunTimeSelectionTable(injectionModel, cloudInjection, dictionary);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::cloudInjection::cloudInjection
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    injectionModel(type(), owner, dict),
+    particlesPerParcel_(readScalar(coeffs_.lookup("particlesPerParcel"))),
+    rndGen_(label(0)),
+    parcelPDF_(pdfs::pdf::New(coeffs_.subDict("parcelPDF"), rndGen_)),
+    diameter_(owner.film().nCells(), 0.0)
+{
+    forAll(diameter_, faceI)
+    {
+        diameter_[faceI] = parcelPDF_->sample();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::cloudInjection::~cloudInjection()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::surfaceFilmModels::cloudInjection::correct
+(
+    scalarField& massToInject,
+    scalarField& diameterToInject
+)
+{
+    const scalar pi = constant::mathematical::pi;
+    const scalarField& rhoFilm = owner().rho();
+
+    // Collect the data to be transferred
+    forAll(massToInject, cellI)
+    {
+        scalar rho = rhoFilm[cellI];
+        scalar diam = diameter_[cellI];
+        scalar minMass = particlesPerParcel_*rho*pi/6*pow3(diam);
+
+        if (massToInject[cellI] > minMass)
+        {
+            // All mass can be injected - set particle diameter
+            diameterToInject[cellI] = diameter_[cellI];
+
+            // Retrieve new particle diameter sample
+            diameter_[cellI] = parcelPDF_->sample();
+        }
+        else
+        {
+            // Mass below minimum threshold - cannot be injected
+            massToInject[cellI] = 0.0;
+            diameterToInject[cellI] = -1.0;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/kinematic/injectionModel/cloudInjection/cloudInjection.H b/src/surfaceFilmModels/submodels/kinematic/injectionModel/cloudInjection/cloudInjection.H
new file mode 100644
index 00000000000..01c9f5eedb9
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/kinematic/injectionModel/cloudInjection/cloudInjection.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::cloudInjection
+
+Description
+    Cloud injection model
+
+SourceFiles
+    cloudInjection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cloudInjection_H
+#define cloudInjection_H
+
+#include "injectionModel.H"
+#include "pdf.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class cloudInjection Declaration
+\*---------------------------------------------------------------------------*/
+
+class cloudInjection
+:
+    public injectionModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        cloudInjection(const cloudInjection&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const cloudInjection&);
+
+
+protected:
+
+    // Protected data
+
+        //- Number of particles per parcel
+        scalar particlesPerParcel_;
+
+        //- Random number generator
+        Random rndGen_;
+
+        //- Parcel size PDF model
+        const autoPtr<pdfs::pdf> parcelPDF_;
+
+        //- Diameters of particles to inject into the cloud
+        scalarList diameter_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cloudInjection");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        cloudInjection(const surfaceFilmModel& owner, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~cloudInjection();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                scalarField& massToInject,
+                scalarField& diameterToInject
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModel.C b/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModel.C
index 96d59420ca4..fe3e0d9e771 100644
--- a/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModel.C
+++ b/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModel.C
@@ -41,8 +41,7 @@ namespace Foam
 
 Foam::surfaceFilmModels::injectionModel::injectionModel
 (
-    const surfaceFilmModel& owner,
-    const dictionary& dict
+    const surfaceFilmModel& owner
 )
 :
     owner_(owner),
diff --git a/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModel.H b/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModel.H
index bc140539bfe..a05d028ffde 100644
--- a/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModel.H
+++ b/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModel.H
@@ -97,7 +97,7 @@ public:
     // Constructors
 
         //- Construct null
-        injectionModel(const surfaceFilmModel& owner, const dictionary& dict);
+        injectionModel(const surfaceFilmModel& owner);
 
         //- Construct from type name, dictionary and surface film model
         injectionModel
diff --git a/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModelNew.C b/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModelNew.C
index 07bfb64cdf2..7bddebcc8ed 100644
--- a/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModelNew.C
+++ b/src/surfaceFilmModels/submodels/kinematic/injectionModel/injectionModel/injectionModelNew.C
@@ -34,7 +34,7 @@ Foam::surfaceFilmModels::injectionModel::New
     const dictionary& dict
 )
 {
-    const word modelType(dict.lookup("injectionModel"));
+    word modelType(dict.lookup("injectionModel"));
 
     Info<< "    Selecting injectionModel " << modelType << endl;
 
diff --git a/src/surfaceFilmModels/submodels/kinematic/injectionModel/noInjection/noInjection.C b/src/surfaceFilmModels/submodels/kinematic/injectionModel/noInjection/noInjection.C
index 39d2c347069..f696f1a1bcd 100644
--- a/src/surfaceFilmModels/submodels/kinematic/injectionModel/noInjection/noInjection.C
+++ b/src/surfaceFilmModels/submodels/kinematic/injectionModel/noInjection/noInjection.C
@@ -43,10 +43,10 @@ namespace Foam
 Foam::surfaceFilmModels::noInjection::noInjection
 (
     const surfaceFilmModel& owner,
-    const dictionary& dict
+    const dictionary&
 )
 :
-    injectionModel(owner, dict)
+    injectionModel(owner)
 {}
 
 
diff --git a/src/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.C b/src/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.C
index 58fbd9fe5af..97d9bb00aae 100644
--- a/src/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.C
+++ b/src/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.C
@@ -43,10 +43,10 @@ namespace Foam
 Foam::surfaceFilmModels::removeInjection::removeInjection
 (
     const surfaceFilmModel& owner,
-    const dictionary& dict
+    const dictionary&
 )
 :
-    injectionModel(owner, dict)
+    injectionModel(owner)
 {}
 
 
diff --git a/src/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.H b/src/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.H
index 4fc4a63fc4c..9c651560f82 100644
--- a/src/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.H
+++ b/src/surfaceFilmModels/submodels/kinematic/injectionModel/removeInjection/removeInjection.H
@@ -66,7 +66,7 @@ private:
 public:
 
     //- Runtime type information
-    TypeName("remove");
+    TypeName("removeInjection");
 
 
     // Constructors
diff --git a/src/surfaceFilmModels/submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.C b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.C
new file mode 100644
index 00000000000..682971f3a62
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.C
@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constantHeatTransfer.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace surfaceFilmModels
+    {
+        defineTypeNameAndDebug(constantHeatTransfer, 0);
+        addToRunTimeSelectionTable
+        (
+            heatTransferModel,
+            constantHeatTransfer,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::constantHeatTransfer::constantHeatTransfer
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    heatTransferModel(typeName, owner, dict),
+    c0_(readScalar(coeffs_.lookup("c0")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::constantHeatTransfer::~constantHeatTransfer()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::surfaceFilmModels::constantHeatTransfer::correct()
+{
+    // do nothing
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::surfaceFilmModels::constantHeatTransfer::h() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "htc",
+                owner_.time().timeName(),
+                owner_.film(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            owner_.film(),
+            dimensionedScalar
+            (
+                "c0",
+                dimEnergy/dimTime/sqr(dimLength)/dimTemperature,
+                c0_
+            ),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+}
+
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.H b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.H
new file mode 100644
index 00000000000..537e70fae20
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/constantHeatTransfer/constantHeatTransfer.H
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::constantHeatTransfer
+
+Description
+    Constant heat transfer model
+
+SourceFiles
+    constantHeatTransfer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantHeatTransfer_H
+#define constantHeatTransfer_H
+
+#include "heatTransferModel.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class constantHeatTransfer Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantHeatTransfer
+:
+    public heatTransferModel
+{
+private:
+
+    // Private data
+
+        //- Constant heat transfer coefficient [W/m2/K]
+        scalar c0_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        constantHeatTransfer(const constantHeatTransfer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const constantHeatTransfer&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constant");
+
+
+    // Constructors
+
+        //- Construct from surface film model and dictionary
+        constantHeatTransfer
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~constantHeatTransfer();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct();
+
+            //- Return the heat transfer coefficient [W/m2/K]
+            virtual tmp<volScalarField> h() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.C b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.C
new file mode 100644
index 00000000000..b7ee1852259
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.C
@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "heatTransferModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace surfaceFilmModels
+    {
+        defineTypeNameAndDebug(heatTransferModel, 0);
+        defineRunTimeSelectionTable(heatTransferModel, dictionary);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::heatTransferModel::heatTransferModel
+(
+    const surfaceFilmModel& owner
+)
+:
+    owner_(owner)
+{}
+
+
+Foam::surfaceFilmModels::heatTransferModel::heatTransferModel
+(
+    const word& type,
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    owner_(owner),
+    coeffs_(dict.subDict(type + "Coeffs"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::heatTransferModel::~heatTransferModel()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.H b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.H
new file mode 100644
index 00000000000..25caf051d91
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModel.H
@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::heatTransferModel
+
+Description
+    Base class for heat transfer models
+
+SourceFiles
+    heatTransferModelI.H
+    heatTransferModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef heatTransferModel_H
+#define heatTransferModel_H
+
+#include "surfaceFilmModel.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class heatTransferModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class heatTransferModel
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        heatTransferModel(const heatTransferModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const heatTransferModel&);
+
+
+protected:
+
+    // Protected data
+
+        //- Reference to the owner surface film model
+        const surfaceFilmModel& owner_;
+
+        //- Model coefficients dictionary
+        dictionary coeffs_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("heatTransferModel");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             heatTransferModel,
+             dictionary,
+             (
+                 const surfaceFilmModel& owner,
+                 const dictionary& dict
+             ),
+             (owner, dict)
+         );
+
+    // Constructors
+
+        //- Construct null
+        heatTransferModel(const surfaceFilmModel& owner);
+
+        //- Construct from type name, dictionary and surface film model
+        heatTransferModel
+        (
+            const word& type,
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected phase change model
+        static autoPtr<heatTransferModel> New
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~heatTransferModel();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the reference to the owner surface film model
+            inline const surfaceFilmModel& owner() const;
+
+            //- Return the model coefficients dictionary
+            inline const dictionary& coeffs() const;
+
+
+        // Evolution
+
+            //- Correct
+            virtual void correct() = 0;
+
+            //- Return the heat transfer coefficient [W/m2/K]
+            virtual tmp<volScalarField> h() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "heatTransferModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelI.H b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelI.H
new file mode 100644
index 00000000000..c9a99ffe298
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelI.H
@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "heatTransferModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline const Foam::surfaceFilmModels::surfaceFilmModel&
+Foam::surfaceFilmModels::heatTransferModel::owner() const
+{
+    return owner_;
+}
+
+
+inline const Foam::dictionary&
+Foam::surfaceFilmModels::heatTransferModel::coeffs() const
+{
+    return coeffs_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelNew.C b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelNew.C
new file mode 100644
index 00000000000..e16c40afaf0
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/heatTransferModel/heatTransferModelNew.C
@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "heatTransferModel.H"
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::surfaceFilmModels::heatTransferModel>
+Foam::surfaceFilmModels::heatTransferModel::New
+(
+    const surfaceFilmModel& model,
+    const dictionary& dict
+)
+{
+    word modelType(dict.lookup("heatTransferModel"));
+
+    Info<< "    Selecting heatTransferModel " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "heatTransferModel::New(const surfaceFilmModel&, const dictionary&)"
+        )   << "Unknown heatTransferModel type " << modelType << nl << nl
+            << "Valid heatTransferModel types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<heatTransferModel>(cstrIter()(model, dict));
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.C b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.C
new file mode 100644
index 00000000000..cfe6d577fe7
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.C
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "mappedConvectiveHeatTransfer.H"
+#include "zeroGradientFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "kinematicSingleLayer.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace surfaceFilmModels
+    {
+        defineTypeNameAndDebug(mappedConvectiveHeatTransfer, 0);
+        addToRunTimeSelectionTable
+        (
+            heatTransferModel,
+            mappedConvectiveHeatTransfer,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::mappedConvectiveHeatTransfer::
+mappedConvectiveHeatTransfer
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+//    heatTransferModel(typeName, owner, dict),
+    heatTransferModel(owner),
+    htcConvPrimary_
+    (
+        IOobject
+        (
+            "htcConv",
+            owner.time().timeName(),
+            owner.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        owner.mesh()
+    ),
+    htcConvFilm_
+    (
+        IOobject
+        (
+            htcConvPrimary_.name(), // must have same name as above for mapping
+            owner.time().timeName(),
+            owner.film(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        owner.film(),
+        dimensionedScalar("zero", dimMass/pow3(dimTime)/dimTemperature, 0.0),
+        refCast<const kinematicSingleLayer>(owner).pSp().boundaryField().types()
+    )
+{
+    // Update the primary-side convective heat transfer coefficient
+    htcConvPrimary_.correctBoundaryConditions();
+
+    // Pull the data from the primary region via direct mapped BCs
+    htcConvFilm_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::mappedConvectiveHeatTransfer::
+~mappedConvectiveHeatTransfer()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::surfaceFilmModels::mappedConvectiveHeatTransfer::correct()
+{
+    // Update the primary-side convective heat transfer coefficient
+    htcConvPrimary_.correctBoundaryConditions();
+
+    // Pull the data from the primary region via direct mapped BCs
+    htcConvFilm_.correctBoundaryConditions();
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::surfaceFilmModels::mappedConvectiveHeatTransfer::h() const
+{
+    return htcConvFilm_;
+}
+
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.H b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.H
new file mode 100644
index 00000000000..8f02cb61ba9
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveHeatTransfer.H
@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::mappedConvectiveHeatTransfer
+
+Description
+    Convective heat transfer model based on a re-working of a Nusselt number
+    correlation
+
+SourceFiles
+    mappedConvectiveHeatTransfer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef mappedConvectiveHeatTransfer_H
+#define mappedConvectiveHeatTransfer_H
+
+#include "heatTransferModel.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                Class mappedConvectiveHeatTransfer Declaration
+\*---------------------------------------------------------------------------*/
+
+class mappedConvectiveHeatTransfer
+:
+    public heatTransferModel
+{
+private:
+
+    // Private data
+
+        //- Heat transfer coefficient - primary region [W/m2/K]
+        volScalarField htcConvPrimary_;
+
+        //- Heat transfer coefficient - film region [W/m2/K]
+        //  Assumes that the primary regtion to film region boundaries are
+        //  described as directMappedPushed types
+        volScalarField htcConvFilm_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        mappedConvectiveHeatTransfer(const mappedConvectiveHeatTransfer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const mappedConvectiveHeatTransfer&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("mappedConvectiveHeatTransfer");
+
+
+    // Constructors
+
+        //- Construct from surface film model and dictionary
+        mappedConvectiveHeatTransfer
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~mappedConvectiveHeatTransfer();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct();
+
+            //- Return the heat transfer coefficient [W/m2/K]
+            virtual tmp<volScalarField> h() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C
index 551bfb2a443..67094e226dc 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.C
@@ -43,10 +43,10 @@ namespace Foam
 Foam::surfaceFilmModels::noPhaseChange::noPhaseChange
 (
     const surfaceFilmModel& owner,
-    const dictionary& dict
+    const dictionary&
 )
 :
-    phaseChangeModel(owner, dict)
+    phaseChangeModel(owner)
 {}
 
 
@@ -58,7 +58,11 @@ Foam::surfaceFilmModels::noPhaseChange::~noPhaseChange()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::surfaceFilmModels::noPhaseChange::correct()
+void Foam::surfaceFilmModels::noPhaseChange::correct
+(
+    const scalar,
+    scalarField&
+)
 {
     // do nothing
 }
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.H b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.H
index 487e3c2b3cf..22a2a17e475 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.H
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/noPhaseChange/noPhaseChange.H
@@ -84,7 +84,7 @@ public:
         // Evolution
 
             //- Correct
-            virtual void correct();
+            virtual void correct(const scalar dt, scalarField& dMass);
 };
 
 
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.C b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.C
index e2a0b93a754..7e7e310d611 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.C
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.C
@@ -41,8 +41,7 @@ namespace Foam
 
 Foam::surfaceFilmModels::phaseChangeModel::phaseChangeModel
 (
-    const surfaceFilmModel& owner,
-    const dictionary& dict
+    const surfaceFilmModel& owner
 )
 :
     owner_(owner),
@@ -59,17 +58,7 @@ Foam::surfaceFilmModels::phaseChangeModel::phaseChangeModel
 :
     owner_(owner),
     coeffs_(dict.subDict(type + "Coeffs"))
-{
-    WarningIn
-    (
-        "phaseChangeModel::phaseChangeModel"
-        "("
-            "const word&, "
-            "const surfaceFilmModel&, "
-            "const dictionary&"
-        ")"
-    )   << "Phase change models not implemented!" << endl;
-}
+{}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.H b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.H
index 87a08b10169..8421520d98d 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.H
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModel.H
@@ -25,6 +25,7 @@ Class
     Foam::phaseChangeModel
 
 Description
+    Phase change model for surface film modelling.
 
 SourceFiles
     phaseChangeModelI.H
@@ -37,6 +38,7 @@ SourceFiles
 
 #include "surfaceFilmModel.H"
 #include "runTimeSelectionTables.H"
+#include "scalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -96,7 +98,7 @@ public:
     // Constructors
 
         //- Construct null
-        phaseChangeModel(const surfaceFilmModel& owner, const dictionary& dict);
+        phaseChangeModel(const surfaceFilmModel& owner);
 
         //- Construct from type name, dictionary and surface film model
         phaseChangeModel
@@ -134,8 +136,8 @@ public:
 
         // Evolution
 
-            //- Correct // TODO
-            virtual void correct() = 0;
+            //- Correct
+            virtual void correct(const scalar dt, scalarField& dMass) = 0;
 };
 
 
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModelNew.C b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModelNew.C
index 1ad7bd0b717..7cc04c0027d 100644
--- a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModelNew.C
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/phaseChangeModel/phaseChangeModelNew.C
@@ -34,7 +34,7 @@ Foam::surfaceFilmModels::phaseChangeModel::New
     const dictionary& dict
 )
 {
-    const word modelType(dict.lookup("phaseChangeModel"));
+    word modelType(dict.lookup("phaseChangeModel"));
 
     Info<< "    Selecting phaseChangeModel " << modelType << endl;
 
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
new file mode 100644
index 00000000000..9011d02fcab
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "standardPhaseChange.H"
+#include "addToRunTimeSelectionTable.H"
+#include "thermoSingleLayer.H"
+#include "specie.H"
+#include "heatTransferModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace surfaceFilmModels
+    {
+        defineTypeNameAndDebug(standardPhaseChange, 0);
+        addToRunTimeSelectionTable
+        (
+            phaseChangeModel,
+            standardPhaseChange,
+            dictionary
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::surfaceFilmModels::standardPhaseChange::Sh
+(
+    const scalar Re,
+    const scalar Sc
+) const
+{
+    if (Re < 5.0E+05)
+    {
+        return 0.664*sqrt(Re)*cbrt(Sc);
+    }
+    else
+    {
+        return 0.037*pow(Re, 0.8)*cbrt(Sc);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::standardPhaseChange::standardPhaseChange
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    phaseChangeModel(typeName, owner, dict),
+    deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
+    L_(readScalar(coeffs_.lookup("L")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::standardPhaseChange::~standardPhaseChange()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::surfaceFilmModels::standardPhaseChange::correct
+(
+    const scalar dt,
+    scalarField& dMass
+)
+{
+    const thermoSingleLayer& film = refCast<const thermoSingleLayer>(owner_);
+
+    // set local thermo properties
+    const SLGThermo& thermo = film.thermo();
+    const label liqId = film.liquidId();
+    const liquid& liq = thermo.liquids().properties()[liqId];
+    const label vapId = thermo.carrierId(thermo.liquids().components()[liqId]);
+
+    // retrieve fields from film model
+    const scalarField& delta = film.delta();
+    const scalarField& YInf = film.YPrimary()[vapId];
+    const scalarField& pInf = film.pPrimary();
+    const scalarField& T = film.T();
+    const scalarField& Ts = film.Ts();
+    const scalarField& Tw = film.Tw();
+    const scalarField& TInf = film.TPrimary();
+    const scalarField& rho = film.rho();
+    const scalarField& mu = film.mu();
+    const scalarField& hs = film.hs();
+    const scalarField& magSf = film.magSf();
+    const scalarField hInf = film.htcs().h();
+    const scalarField hFilm = film.htcw().h();
+    const vectorField dU = film.UPrimary() - film.Us();
+
+    // Reynolds number
+    const scalarField Re = rho*mag(dU)*L_/mu;
+
+    // molecular weight of vapour
+    const scalar Wvap = thermo.carrier().W(vapId);
+
+    // molecular weight of liquid
+    const scalar Wliq = liq.W();
+
+    forAll(dMass, cellI)
+    {
+        if (delta[cellI] > deltaMin_)
+        {
+            // cell pressure
+            const scalar pc = pInf[cellI];
+
+            // saturation pressure
+            const scalar pSat = liq.pv(pc, Ts[cellI]);
+
+            // calculate mass transfer
+            if (pSat > pc)
+            {
+                // boiling
+                const scalar qDotInf = hInf[cellI]*(TInf[cellI] - T[cellI]);
+                const scalar qDotFilm = hFilm[cellI]*(T[cellI] - Tw[cellI]);
+                dMass[cellI] +=
+                    max
+                    (
+                        0.0,
+                        dt*magSf[cellI]/hs[cellI]*(qDotInf - qDotFilm)
+                    );
+            }
+            else
+            {
+                // vapour mass fraction at interface
+                const scalar Ys = Wliq*pSat/(Wliq*pSat + Wvap*(pc - pSat));
+
+                // bulk gas average density
+                const scalar rhoAve = pc/(specie::RR*Ts[cellI]);
+
+                // vapour diffusivity [m2/s]
+                const scalar Dab = liq.D(pc, Ts[cellI]);
+
+                // Schmidt number
+                const scalar Sc = mu[cellI]/(rho[cellI]*(Dab + ROOTVSMALL));
+
+                // Sherwood number
+                const scalar Sh = this->Sh(Re[cellI], Sc);
+
+                // mass transfer coefficient [m/s]
+                const scalar hm = Sh*Dab/L_;
+
+                // add mass contribution to source
+                dMass[cellI] =
+                    max
+                    (
+                        0.0,
+                        dt*magSf[cellI]*rhoAve*hm*(Ys - YInf[cellI])/(1.0 - Ys)
+                    );
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
new file mode 100644
index 00000000000..dd16488e03f
--- /dev/null
+++ b/src/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::standardPhaseChange
+
+Description
+
+
+SourceFiles
+    standardPhaseChange.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef standardPhaseChange_H
+#define standardPhaseChange_H
+
+#include "phaseChangeModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class standardPhaseChange Declaration
+\*---------------------------------------------------------------------------*/
+
+class standardPhaseChange
+:
+    public phaseChangeModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        standardPhaseChange(const standardPhaseChange&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const standardPhaseChange&);
+
+
+protected:
+
+    // Protected data
+
+        //- Minimum film height for model to be active
+        const scalar deltaMin_;
+
+        //- Length scalae / [m]
+        const scalar L_;
+
+
+    // Protected member functions
+
+        //- Return Sherwood number as a function of Reynolds and Schmidt numbers
+        scalar Sh(const scalar Re, const scalar Sc) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("standardPhaseChange");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        standardPhaseChange
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~standardPhaseChange();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct(const scalar dt, scalarField& dMass);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.C b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.C
index 36f28993533..cdfc81d774e 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.C
@@ -1,1062 +1,1302 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "kinematicSingleLayer.H"
-#include "fvm.H"
-#include "fvcDiv.H"
-#include "fvcLaplacian.H"
-#include "fvcSnGrad.H"
-#include "fvcReconstruct.H"
-#include "fvcVolumeIntegrate.H"
-#include "addToRunTimeSelectionTable.H"
-#include "directMappedWallPolyPatch.H"
-
-// Sub-models
-#include "injectionModel.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    namespace surfaceFilmModels
-    {
-        defineTypeNameAndDebug(kinematicSingleLayer, 0);
-        addToRunTimeSelectionTable
-        (
-            surfaceFilmModel,
-            kinematicSingleLayer,
-            mesh
-        );
-    }
-}
-
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-bool Foam::surfaceFilmModels::kinematicSingleLayer::read()
-{
-    if (surfaceFilmModel::read())
-    {
-        const dictionary& solution = filmRegion_.solutionDict().subDict("PISO");
-        solution.lookup("momentumPredictor") >> momentumPredictor_;
-        solution.lookup("nOuterCorr") >> nOuterCorr_;
-        solution.lookup("nCorr") >> nCorr_;
-        solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
-
-        coeffs_.lookup("Cf") >> Cf_;
-        coeffs_.lookup("deltaStable") >> deltaStable_;
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::initialise()
-{
-    if (debug)
-    {
-        Pout<< "kinematicSingleLayer::initialise()" << endl;
-    }
-
-    label nBoundaryFaces = 0;
-    DynamicList<label> primaryPatchIDs;
-    DynamicList<label> filmBottomPatchIDs;
-    const polyBoundaryMesh& bm = filmRegion_.boundaryMesh();
-    forAll(bm, patchI)
-    {
-        const polyPatch& pp = bm[patchI];
-        if (isA<directMappedWallPolyPatch>(pp))
-        {
-            if (debug)
-            {
-                Pout<< "found " << directMappedWallPolyPatch::typeName
-                    <<  " " << pp.name() << endl;
-            }
-
-            filmBottomPatchIDs.append(patchI);
-            const directMappedWallPolyPatch& dwpp =
-                refCast<const directMappedWallPolyPatch>(pp);
-
-            primaryPatchIDs.append
-            (
-                mesh_.boundaryMesh().findPatchID(dwpp.samplePatch())
-            );
-
-            const labelList& fCells = pp.faceCells();
-            nBoundaryFaces += fCells.size();
-
-            // Cache patch normals
-            UIndirectList<vector>(nHat_, fCells) = pp.faceNormals();
-
-            // Cache mesh face areas
-            UIndirectList<scalar>(magSf_, fCells) = mag(pp.faceAreas());
-        }
-    }
-    nHat_.correctBoundaryConditions();
-    magSf_.correctBoundaryConditions();
-
-    primaryPatchIDs_.transfer(primaryPatchIDs);
-    filmBottomPatchIDs_.transfer(filmBottomPatchIDs);
-
-    if (nBoundaryFaces == 0)
-    {
-        WarningIn("kinematicSingleLayer::initialise()")
-            << "Film model being applied without direct mapped boundary "
-            << "conditions" << endl;
-    }
-
-    if (nBoundaryFaces != filmRegion_.nCells())
-    {
-        FatalErrorIn("kinematicSingleLayer::initialise()")
-            << "Number of primary region coupled boundary faces not equal to "
-            << "the number of cells in the film region" << nl
-            << abort(FatalError);
-    }
-
-    scalarField topMagSf(magSf_.size(), 0.0);
-    filmTopPatchIDs_.setSize(filmBottomPatchIDs_.size(), -1);
-    forAll(filmBottomPatchIDs_, i)
-    {
-        const label patchI = filmBottomPatchIDs_[i];
-        const polyPatch& ppBottom = bm[patchI];
-        if (ppBottom.size() > 0)
-        {
-            label cellId = bm[patchI].faceCells()[0];
-            const cell& cFaces = filmRegion_.cells()[cellId];
-
-            label faceBottom = ppBottom.start();
-            label faceTop =
-                 cFaces.opposingFaceLabel(faceBottom, filmRegion_.faces());
-
-            label topPatchI = bm.whichPatch(faceTop);
-            filmTopPatchIDs_[i] = topPatchI;
-            const polyPatch& ppTop = bm[topPatchI];
-            UIndirectList<scalar>(topMagSf, ppTop.faceCells()) =
-                mag(ppTop.faceAreas());
-        }
-    }
-
-    Pstream::listCombineGather(filmTopPatchIDs_, maxEqOp<label>());
-    Pstream::listCombineScatter(filmTopPatchIDs_);
-
-    magSf_.field() = 0.5*(magSf_ + topMagSf);
-    magSf_.correctBoundaryConditions();
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::
-resetPrimaryRegionSourceTerms()
-{
-    rhoSpPrimary_ == dimensionedScalar("zero", rhoSp_.dimensions(), 0.0);
-    USpPrimary_ == dimensionedVector("zero", USp_.dimensions(), vector::zero);
-    pSpPrimary_ == dimensionedScalar("zero", pSp_.dimensions(), 0.0);
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::
-transferPrimaryRegionFields()
-{
-    // Update pressure and velocity from primary region via direct mapped
-    // (coupled) boundary conditions
-    UPrimary_.correctBoundaryConditions();
-    pPrimary_.correctBoundaryConditions();
-
-    // Retrieve the source fields from the primary region via direct mapped
-    // (coupled) boundary conditions
-    // - fields require transfer of values for both patch AND to push the
-    //   values into the first layer of internal cells
-    rhoSp_.correctBoundaryConditions();
-    USp_.correctBoundaryConditions();
-    pSp_.correctBoundaryConditions();
-
-    // Convert accummulated source terms into per unit area per unit time
-    // Note: boundary values will still have original (neat) values
-    const scalar deltaT = filmRegion_.time().deltaTValue();
-    rhoSp_.field() /= magSf_*deltaT;
-    USp_.field() /= magSf_*deltaT;
-    pSp_.field() /= magSf_*deltaT;
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::surfaceFilmModels::kinematicSingleLayer::pu()
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "pu",
-                filmRegion_.time().timeName(),
-                filmRegion_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            pPrimary_                  // pressure (mapped from primary region)
-          + pSp_                           // accumulated particle impingement
-          - fvc::laplacian(sigma_, delta_) // surface tension
-        )
-    );
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::surfaceFilmModels::kinematicSingleLayer::pp()
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "pp",
-                filmRegion_.time().timeName(),
-                filmRegion_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-           -rho_*gNormClipped() // hydrostatic effect only
-        )
-    );
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::correctDetachedFilm()
-{
-    massForPrimary_ == dimensionedScalar("zero", dimMass, 0.0);
-    diametersForPrimary_ == dimensionedScalar("zero", dimLength, -1.0);
-
-    const scalarField gNorm = this->gNorm();
-
-    forAll(gNorm, i)
-    {
-        if (gNorm[i] > SMALL)
-        {
-            scalar ddelta = max(0.0, delta_[i] - deltaStable_.value());
-            massForPrimary_[i] = ddelta*rho_[i]*magSf_[i];
-        }
-    }
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::updateSubmodels()
-{
-    correctDetachedFilm();
-
-    // Update injection model - mass returned is actual mass injected
-    injection_->correct(massForPrimary_, diametersForPrimary_);
-
-    // Update cumulative detached mass counter
-    detachedMass_ += sum(massForPrimary_.field());
-
-    // Push values to boundaries ready for transfer to the primary region
-    massForPrimary_.correctBoundaryConditions();
-    diametersForPrimary_.correctBoundaryConditions();
-
-    // Update source fields
-    const dimensionedScalar deltaT = filmRegion_.time().deltaT();
-    rhoSp_ -= massForPrimary_/magSf_/deltaT;
-    USp_ -= massForPrimary_*U_/magSf_/deltaT;
-}
-
-
-Foam::scalar
-Foam::surfaceFilmModels::kinematicSingleLayer::CourantNumber() const
-{
-    scalar CoNum = 0.0;
-    scalar meanCoNum = 0.0;
-
-    if (filmRegion_.nInternalFaces())
-    {
-        const scalar deltaT = time_.deltaTValue();
-
-        surfaceScalarField SfUfbyDelta =
-            filmRegion_.surfaceInterpolation::deltaCoeffs()*mag(phi_)
-           /fvc::interpolate
-            (
-                rho_*(delta_ + dimensionedScalar("SMALL", dimLength, SMALL))
-            );
-
-        CoNum = max(SfUfbyDelta/filmRegion_.magSf()).value()*deltaT;
-
-        meanCoNum = (sum(SfUfbyDelta)/sum(filmRegion_.magSf())).value()*deltaT;
-    }
-
-    Info<< "    Courant number mean: " << meanCoNum << " max: " << CoNum
-        << endl;
-
-    return CoNum;
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::continuityCheck()
-{
-    const volScalarField deltaRho0 = deltaRho_;
-
-    solveContinuity();
-
-    if (debug)
-    {
-        volScalarField mass = deltaRho_*magSf_;
-        dimensionedScalar totalMass =
-            fvc::domainIntegrate(mass)
-          + dimensionedScalar("SMALL", dimMass*dimVolume, ROOTVSMALL);
-
-        scalar sumLocalContErr =
-            (
-                fvc::domainIntegrate(mag(mass - magSf_*deltaRho0))/totalMass
-            ).value();
-
-        scalar globalContErr =
-            (
-                fvc::domainIntegrate(mass - magSf_*deltaRho0)/totalMass
-            ).value();
-
-        cumulativeContErr_ += globalContErr;
-
-        Info<< "Surface film: " << type() << nl
-            << "    time step continuity errors: sum local = "
-            << sumLocalContErr << ", global = " << globalContErr
-            << ", cumulative = " << cumulativeContErr_ << endl;
-    }
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::solveContinuity()
-{
-    if (debug)
-    {
-        Info<< "kinematicSingleLayer::solveContinuity()" << endl;
-    }
-
-    solve
-    (
-        fvm::ddt(deltaRho_)
-      + fvc::div(phi_)
-     ==
-        rhoSp_
-    );
-}
-
-
-Foam::tmp<Foam::fvVectorMatrix>
-Foam::surfaceFilmModels::kinematicSingleLayer::tau
-(
-    volVectorField& U
-) const
-{
-    DimensionedField<vector, volMesh> Uw
-    (
-        IOobject
-        (
-            "Uw",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        filmRegion_,
-        dimensionedVector("zero", dimVelocity, vector::zero)
-    );
-
-    forAll(filmBottomPatchIDs_, i)
-    {
-        label patchI = filmBottomPatchIDs_[i];
-        const polyPatch& pp = filmRegion_.boundaryMesh()[patchI];
-        UIndirectList<vector>(Uw, pp.faceCells()) =
-            UPrimary_.boundaryField()[patchI];
-    }
-    Uw -= nHat_*(Uw & nHat_);
-
-    volVectorField Us
-    (
-        IOobject
-        (
-            "Us",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        filmRegion_,
-        dimensionedVector("zero", dimVelocity, vector::zero)
-    );
-
-    // TODO: use wall functions
-    Us.dimensionedInternalField() = UPrimary_.dimensionedInternalField();
-    Us -= nHat_*(Us & nHat_);
-
-    volScalarField Cs("Cs", rho_*Cf_*mag(Us - U_));
-    volScalarField Cw
-    (
-        "Cw",
-        mu_/(0.5*(delta_ + dimensionedScalar("SMALL", dimLength, SMALL)))
-    );
-    Cw.min(1.0e+06);
-
-    return
-    (
-       - fvm::Sp(Cs, U) + Cs*Us
-       - fvm::Sp(Cw, U) + Cw*Uw
-    );
-}
-
-
-Foam::tmp<Foam::fvVectorMatrix>
-Foam::surfaceFilmModels::kinematicSingleLayer::solveMomentum
-(
-    const volScalarField& pu,
-    const volScalarField& pp
-)
-{
-    if (debug)
-    {
-        Info<< "kinematicSingleLayer::solveMomentum()" << endl;
-    }
-
-    // Momentum
-    tmp<fvVectorMatrix> tUEqn
-    (
-        fvm::ddt(deltaRho_, U_)
-      + fvm::div(phi_, U_)
-     ==
-        USp_
-      + tau(U_)
-    );
-
-    fvVectorMatrix& UEqn = tUEqn();
-
-    UEqn.relax();
-
-    if (momentumPredictor_)
-    {
-        solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-              - fvc::interpolate(delta_)
-              * (
-                    filmRegion_.magSf()
-                  * (
-                        fvc::snGrad(pu, "snGrad(p)")
-                      + fvc::snGrad(pp, "snGrad(p)")*fvc::interpolate(delta_)
-                      + fvc::snGrad(delta_)*fvc::interpolate(pp)
-                    )
-                  - (fvc::interpolate(rho_*gTan()) & filmRegion_.Sf())
-                )
-            )
-        );
-
-        // Remove any patch-normal components of velocity
-        U_ -= nHat_*(nHat_ & U_);
-        U_.correctBoundaryConditions();
-    }
-
-    return tUEqn;
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::solveThickness
-(
-    const volScalarField& pu,
-    const volScalarField& pp,
-    const fvVectorMatrix& UEqn
-)
-{
-    if (debug)
-    {
-        Info<< "kinematicSingleLayer::solveThickness()" << endl;
-    }
-
-    volScalarField rUA = 1.0/UEqn.A();
-    U_ = rUA*UEqn.H();
-
-    surfaceScalarField deltarUAf = fvc::interpolate(delta_*rUA);
-    surfaceScalarField rhof = fvc::interpolate(rho_);
-
-    surfaceScalarField phiAdd
-    (
-        "phiAdd",
-        filmRegion_.magSf()
-      * (
-            fvc::snGrad(pu, "snGrad(p)")
-          + fvc::snGrad(pp, "snGrad(p)")*fvc::interpolate(delta_)
-        )
-      - (fvc::interpolate(rho_*gTan()) & filmRegion_.Sf())
-    );
-    constrainFilmField(phiAdd, 0.0);
-
-    surfaceScalarField phid
-    (
-        "phid",
-        (fvc::interpolate(U_*rho_) & filmRegion_.Sf())
-      - deltarUAf*phiAdd*rhof
-    );
-    constrainFilmField(phid, 0.0);
-
-    surfaceScalarField ddrhorUAppf =
-        fvc::interpolate(delta_)*deltarUAf*rhof*fvc::interpolate(pp);
-//    constrainFilmField(ddrhorUAppf, 0.0);
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr_; nonOrth++)
-    {
-        // Film thickness equation
-        fvScalarMatrix deltaEqn
-        (
-            fvm::ddt(rho_, delta_)
-          + fvm::div(phid, delta_)
-          - fvm::laplacian(ddrhorUAppf, delta_)
-         ==
-            rhoSp_
-        );
-
-        deltaEqn.solve();
-
-        if (nonOrth == nNonOrthCorr_)
-        {
-            phiAdd +=
-                fvc::interpolate(pp)
-              * fvc::snGrad(delta_)
-              * filmRegion_.magSf();
-
-            phi_ == deltaEqn.flux();
-        }
-    }
-
-    // Bound film thickness by a minimum of zero
-    delta_.max(0.0);
-
-    // Update U field
-    U_ -= fvc::reconstruct(deltarUAf*phiAdd);
-
-    // Remove any patch-normal components of velocity
-    U_ -= nHat_*(nHat_ & U_);
-
-    U_.correctBoundaryConditions();
-
-    // Continuity check
-    continuityCheck();
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::surfaceFilmModels::kinematicSingleLayer::kinematicSingleLayer
-(
-    const word& modelType,
-    const fvMesh& mesh,
-    const dimensionedVector& g
-)
-:
-    surfaceFilmModel(modelType, mesh, g),
-    filmRegion_
-    (
-        IOobject
-        (
-            filmRegionName_,
-            time_.timeName(),
-            time_,
-            IOobject::MUST_READ_IF_MODIFIED
-        )
-    ),
-    nHat_
-    (
-        IOobject
-        (
-            "nHat",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        filmRegion_,
-        dimensionedVector("zero", dimless, vector::zero),
-        zeroGradientFvPatchVectorField::typeName
-    ),
-    magSf_
-    (
-        IOobject
-        (
-            "magSf",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        filmRegion_,
-        dimensionedScalar("zero", dimArea, 0.0),
-        zeroGradientFvPatchScalarField::typeName
-    ),
-    primaryPatchIDs_(0),
-    filmTopPatchIDs_(0),
-    filmBottomPatchIDs_(0),
-
-    momentumPredictor_(true),
-    nOuterCorr_(-1),
-    nCorr_(-1),
-    nNonOrthCorr_(-1),
-    cumulativeContErr_(0.0),
-
-    Cf_(0.0),
-    deltaStable_("deltaStable", dimLength, 0.0),
-
-    rho_
-    (
-        IOobject
-        (
-            "rhof",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        filmRegion_,
-        dimensionedScalar(coeffs_.lookup("rho"))
-    ),
-    mu_
-    (
-        IOobject
-        (
-            "muf",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        filmRegion_,
-        dimensionedScalar(coeffs_.lookup("mu"))
-    ),
-    sigma_
-    (
-        IOobject
-        (
-            "sigmaf",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        filmRegion_,
-        dimensionedScalar(coeffs_.lookup("sigma"))
-    ),
-
-    delta_
-    (
-        IOobject
-        (
-            "deltaf",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        filmRegion_
-    ),
-    U_
-    (
-        IOobject
-        (
-            "Uf",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        filmRegion_
-    ),
-    deltaRho_
-    (
-        IOobject
-        (
-            delta_.name() + "*" + rho_.name(),
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        delta_*rho_,
-        zeroGradientFvPatchScalarField::typeName
-    ),
-
-    phi_
-    (
-        IOobject
-        (
-            "phi",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        fvc::interpolate(delta_*rho_*U_) & filmRegion_.Sf()
-    ),
-
-    massForPrimary_
-    (
-        IOobject
-        (
-            "massForPrimary",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        filmRegion_,
-        dimensionedScalar("zero", dimMass, 0.0),
-        zeroGradientFvPatchScalarField::typeName
-    ),
-    diametersForPrimary_
-    (
-        IOobject
-        (
-            "diametersForPrimary",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        filmRegion_,
-        dimensionedScalar("zero", dimLength, -1.0),
-        zeroGradientFvPatchScalarField::typeName
-    ),
-
-    USp_
-    (
-        IOobject
-        (
-            "USpf",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        filmRegion_
-    ),
-    pSp_
-    (
-        IOobject
-        (
-            "pSpf",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        filmRegion_
-    ),
-    rhoSp_
-    (
-        IOobject
-        (
-            "rhoSpf",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        filmRegion_,
-        dimensionedScalar("zero", dimMass/dimTime/dimArea, 0.0),
-        pSp_.boundaryField().types()
-    ),
-
-    USpPrimary_
-    (
-        IOobject
-        (
-            USp_.name(), // must have same name as USp_ to enable mapping
-            time_.timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedVector("zero", USp_.dimensions(), vector::zero)
-    ),
-    pSpPrimary_
-    (
-        IOobject
-        (
-            pSp_.name(), // must have same name as pSp_ to enable mapping
-            time_.timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("zero", pSp_.dimensions(), 0.0)
-    ),
-    rhoSpPrimary_
-    (
-        IOobject
-        (
-            rhoSp_.name(), // must have same name as rhoSp_ to enable mapping
-            time_.timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("zero", rhoSp_.dimensions(), 0.0)
-    ),
-
-    UPrimary_
-    (
-        IOobject
-        (
-            "U", // must have same name as U on primary region to enable mapping
-            time_.timeName(),
-            filmRegion_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        filmRegion_
-    ),
-    pPrimary_
-    (
-        IOobject
-        (
-            "p", // must have same name as p on primary region to enable mapping
-            time_.timeName(),
-            filmRegion_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        filmRegion_
-    ),
-
-    injection_(injectionModel::New(*this, coeffs_)),
-
-    addedMass_(0.0),
-    detachedMass_(0.0)
-{
-    read();
-
-    initialise();
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::surfaceFilmModels::kinematicSingleLayer::~kinematicSingleLayer()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-inline const Foam::fvMesh&
-Foam::surfaceFilmModels::kinematicSingleLayer::film() const
-{
-    return filmRegion_;
-}
-
-
-const Foam::labelList&
-Foam::surfaceFilmModels::kinematicSingleLayer::filmBottomPatchIDs() const
-{
-    return filmBottomPatchIDs_;
-}
-
-
-const Foam::labelList&
-Foam::surfaceFilmModels::kinematicSingleLayer::filmTopPatchIDs() const
-{
-    return filmTopPatchIDs_;
-}
-
-
-const Foam::labelList&
-Foam::surfaceFilmModels::kinematicSingleLayer::primaryPatchIDs() const
-{
-    return primaryPatchIDs_;
-}
-
-
-bool Foam::surfaceFilmModels::kinematicSingleLayer::isFilmPatch
-(
-    const label patchI
-) const
-{
-    if (!active_)
-    {
-        return false;
-    }
-
-    forAll(primaryPatchIDs_, i)
-    {
-        if (primaryPatchIDs_[i] == patchI)
-        {
-            return true;
-        }
-    }
-
-    return false;
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::addSources
-(
-    const label patchI,
-    const label faceI,
-    const scalar massSource,
-    const vector& momentumSource,
-    const scalar pressureSource,
-    const scalar energySource
-)
-{
-    if (debug)
-    {
-        Info<< "\nSurface film: " << type() << ": adding to film source:" << nl
-            << "    mass     = " << massSource << nl
-            << "    momentum = " << momentumSource << nl
-            << "    pressure = " << pressureSource << endl;
-    }
-
-    rhoSpPrimary_.boundaryField()[patchI][faceI] += massSource;
-    USpPrimary_.boundaryField()[patchI][faceI] += momentumSource;
-    pSpPrimary_.boundaryField()[patchI][faceI] += pressureSource;
-
-    addedMass_ += massSource;
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::evolveFilm()
-{
-    transferPrimaryRegionFields();
-
-    updateSubmodels();
-
-    // Solve continuity for deltaRho_
-    solveContinuity();
-
-    for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
-    {
-        // Explicit pressure source contribution
-        tmp<volScalarField> pu = this->pu();
-
-        // Implicit pressure source coefficient
-        tmp<volScalarField> pp = this->pp();
-
-        // Solve for momentum for U_
-        tmp<fvVectorMatrix> UEqn = solveMomentum(pu(), pp());
-
-        // Film thickness correction loop
-        for (int corr=1; corr<=nCorr_; corr++)
-        {
-            // Solve thickness for delta_
-            solveThickness(pu(), pp(), UEqn());
-        }
-    }
-
-    // Update deltaRho_ with new delta_
-    deltaRho_ == delta_*rho_;
-
-    // Reset source terms for next time integration
-    resetPrimaryRegionSourceTerms();
-}
-
-
-const Foam::volVectorField&
-Foam::surfaceFilmModels::kinematicSingleLayer::U() const
-{
-    return U_;
-}
-
-
-const Foam::volScalarField&
-Foam::surfaceFilmModels::kinematicSingleLayer::rho() const
-{
-    return rho_;
-}
-
-
-const Foam::volScalarField&
-Foam::surfaceFilmModels::kinematicSingleLayer::T() const
-{
-    FatalErrorIn
-    (
-        "const volScalarField& kinematicSingleLayer::T() const"
-    )   << "T field not available for " << type() << abort(FatalError);
-
-    return volScalarField::null();
-}
-
-
-const Foam::volScalarField&
-Foam::surfaceFilmModels::kinematicSingleLayer::cp() const
-{
-    FatalErrorIn
-    (
-        "const volScalarField& kinematicSingleLayer::cp() const"
-    )   << "cp field not available for " << type() << abort(FatalError);
-
-    return volScalarField::null();
-}
-
-
-
-const Foam::volScalarField&
-Foam::surfaceFilmModels::kinematicSingleLayer::massForPrimary() const
-{
-    return massForPrimary_;
-}
-
-
-const Foam::volScalarField&
-Foam::surfaceFilmModels::kinematicSingleLayer::diametersForPrimary() const
-{
-    return diametersForPrimary_;
-}
-
-
-void Foam::surfaceFilmModels::kinematicSingleLayer::info() const
-{
-    Info<< "\nSurface film: " << type() << endl;
-
-    // Output Courant number for info only - does not change time step
-    CourantNumber();
-
-    Info<< indent << "added mass      = "
-        << returnReduce<scalar>(addedMass_, sumOp<scalar>()) << nl
-        << indent << "current mass    = "
-        << gSum((deltaRho_*magSf_)()) << nl
-        << indent << "detached mass   = "
-        << returnReduce<scalar>(detachedMass_, sumOp<scalar>()) << nl
-        << indent << "min/max(mag(U)) = " << min(mag(U_)).value() << ", "
-        << max(mag(U_)).value() << nl
-        << indent << "min/max(delta)  = " << min(delta_).value() << ", "
-        << max(delta_).value() << nl;
-}
-
-
-// ************************************************************************* //
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kinematicSingleLayer.H"
+#include "fvm.H"
+#include "fvcDiv.H"
+#include "fvcLaplacian.H"
+#include "fvcSnGrad.H"
+#include "fvcReconstruct.H"
+#include "fvcVolumeIntegrate.H"
+#include "addToRunTimeSelectionTable.H"
+#include "directMappedWallPolyPatch.H"
+
+// Sub-models
+#include "injectionModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace surfaceFilmModels
+    {
+        defineTypeNameAndDebug(kinematicSingleLayer, 0);
+        addToRunTimeSelectionTable
+        (
+            surfaceFilmModel,
+            kinematicSingleLayer,
+            mesh
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+bool Foam::surfaceFilmModels::kinematicSingleLayer::read()
+{
+    if (surfaceFilmModel::read())
+    {
+        const dictionary& solution = filmRegion_.solutionDict().subDict("PISO");
+        solution.lookup("momentumPredictor") >> momentumPredictor_;
+        solution.lookup("nOuterCorr") >> nOuterCorr_;
+        solution.lookup("nCorr") >> nCorr_;
+        solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
+
+        coeffs_.lookup("Cf") >> Cf_;
+        coeffs_.lookup("deltaStable") >> deltaStable_;
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::initialise()
+{
+    if (debug)
+    {
+        Pout<< "kinematicSingleLayer::initialise()" << endl;
+    }
+
+    label nBoundaryFaces = 0;
+    DynamicList<label> primaryPatchIDs;
+    DynamicList<label> filmBottomPatchIDs;
+    const polyBoundaryMesh& bm = filmRegion_.boundaryMesh();
+    forAll(bm, patchI)
+    {
+        const polyPatch& pp = bm[patchI];
+        if (isA<directMappedWallPolyPatch>(pp))
+        {
+            if (debug)
+            {
+                Pout<< "found " << directMappedWallPolyPatch::typeName
+                    <<  " " << pp.name() << endl;
+            }
+
+            filmBottomPatchIDs.append(patchI);
+            const directMappedWallPolyPatch& dwpp =
+                refCast<const directMappedWallPolyPatch>(pp);
+
+            primaryPatchIDs.append
+            (
+                mesh_.boundaryMesh().findPatchID(dwpp.samplePatch())
+            );
+
+            const labelList& fCells = pp.faceCells();
+            nBoundaryFaces += fCells.size();
+
+            // Cache patch normals
+            UIndirectList<vector>(nHat_, fCells) = pp.faceNormals();
+
+            // Cache mesh face areas
+            UIndirectList<scalar>(magSf_, fCells) = mag(pp.faceAreas());
+        }
+    }
+    nHat_.correctBoundaryConditions();
+    magSf_.correctBoundaryConditions();
+
+    primaryPatchIDs_.transfer(primaryPatchIDs);
+    filmBottomPatchIDs_.transfer(filmBottomPatchIDs);
+
+    if (nBoundaryFaces == 0)
+    {
+        WarningIn("kinematicSingleLayer::initialise()")
+            << "Film model being applied without direct mapped boundary "
+            << "conditions" << endl;
+    }
+
+    if (nBoundaryFaces != filmRegion_.nCells())
+    {
+        FatalErrorIn("kinematicSingleLayer::initialise()")
+            << "Number of primary region coupled boundary faces not equal to "
+            << "the number of cells in the film region" << nl
+            << abort(FatalError);
+    }
+
+    scalarField topMagSf(magSf_.size(), 0.0);
+    filmTopPatchIDs_.setSize(filmBottomPatchIDs_.size(), -1);
+    forAll(filmBottomPatchIDs_, i)
+    {
+        const label patchI = filmBottomPatchIDs_[i];
+        const polyPatch& ppBottom = bm[patchI];
+        if (ppBottom.size() > 0)
+        {
+            label cellId = bm[patchI].faceCells()[0];
+            const cell& cFaces = filmRegion_.cells()[cellId];
+
+            label faceBottom = ppBottom.start();
+            label faceTop =
+                 cFaces.opposingFaceLabel(faceBottom, filmRegion_.faces());
+
+            label topPatchI = bm.whichPatch(faceTop);
+            filmTopPatchIDs_[i] = topPatchI;
+            const polyPatch& ppTop = bm[topPatchI];
+            UIndirectList<scalar>(topMagSf, ppTop.faceCells()) =
+                mag(ppTop.faceAreas());
+        }
+    }
+
+    Pstream::listCombineGather(filmTopPatchIDs_, maxEqOp<label>());
+    Pstream::listCombineScatter(filmTopPatchIDs_);
+
+    magSf_.field() = 0.5*(magSf_ + topMagSf);
+    magSf_.correctBoundaryConditions();
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::correctThermoFields()
+{
+    if (thermoModel_ == tmConstant)
+    {
+        rho_ == dimensionedScalar(coeffs_.lookup("rho0"));
+        mu_ == dimensionedScalar(coeffs_.lookup("mu0"));
+        sigma_ == dimensionedScalar(coeffs_.lookup("sigma0"));
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "void Foam::surfaceFilmModels::kinematicSingleLayer::"
+            "correctThermo()"
+        )   << "Kinematic surface film must use "
+            << thermoModelTypeNames_[thermoModel_] << "thermodynamics" << endl;
+    }
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::
+resetPrimaryRegionSourceTerms()
+{
+    rhoSpPrimary_ == dimensionedScalar("zero", rhoSp_.dimensions(), 0.0);
+    USpPrimary_ == dimensionedVector("zero", USp_.dimensions(), vector::zero);
+    pSpPrimary_ == dimensionedScalar("zero", pSp_.dimensions(), 0.0);
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::
+transferPrimaryRegionFields()
+{
+    // Update pressure and velocity from primary region via direct mapped
+    // (coupled) boundary conditions
+    UPrimary_.correctBoundaryConditions();
+    pPrimary_.correctBoundaryConditions();
+
+    // Retrieve the source fields from the primary region via direct mapped
+    // (coupled) boundary conditions
+    // - fields require transfer of values for both patch AND to push the
+    //   values into the first layer of internal cells
+    rhoSp_.correctBoundaryConditions();
+    USp_.correctBoundaryConditions();
+    pSp_.correctBoundaryConditions();
+
+    // Convert accummulated source terms into per unit area per unit time
+    // Note: boundary values will still have original (neat) values
+    const scalar deltaT = time_.deltaTValue();
+    rhoSp_.field() /= magSf_*deltaT;
+    USp_.field() /= magSf_*deltaT;
+    pSp_.field() /= magSf_*deltaT;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::surfaceFilmModels::kinematicSingleLayer::pu()
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "pu",
+                time_.timeName(),
+                filmRegion_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            pPrimary_                  // pressure (mapped from primary region)
+          + pSp_                           // accumulated particle impingement
+          - fvc::laplacian(sigma_, delta_) // surface tension
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::surfaceFilmModels::kinematicSingleLayer::pp()
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "pp",
+                time_.timeName(),
+                filmRegion_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+           -rho_*gNormClipped() // hydrostatic effect only
+        )
+    );
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::correctDetachedFilm()
+{
+    massForPrimary_ == dimensionedScalar("zero", dimMass, 0.0);
+    diametersForPrimary_ == dimensionedScalar("zero", dimLength, -1.0);
+
+    const scalarField gNorm = this->gNorm();
+
+    forAll(gNorm, i)
+    {
+        if (gNorm[i] > SMALL)
+        {
+            scalar ddelta = max(0.0, delta_[i] - deltaStable_.value());
+            massForPrimary_[i] =
+                max
+                (
+                    0.0,
+                    ddelta*rho_[i]*magSf_[i] - massPhaseChangeForPrimary_[i]
+                );
+        }
+    }
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::updateSubmodels()
+{
+    correctDetachedFilm();
+
+    // Update injection model - mass returned is actual mass injected
+    injection_->correct(massForPrimary_, diametersForPrimary_);
+
+    // Update cumulative detached mass counter
+    detachedMass_ += sum(massForPrimary_.field());
+
+    // Push values to boundaries ready for transfer to the primary region
+    massForPrimary_.correctBoundaryConditions();
+    diametersForPrimary_.correctBoundaryConditions();
+
+    // Update source fields
+    const dimensionedScalar deltaT = time_.deltaT();
+    rhoSp_ -= massForPrimary_/magSf_/deltaT;
+    USp_ -= massForPrimary_*U_/magSf_/deltaT;
+}
+
+
+Foam::scalar
+Foam::surfaceFilmModels::kinematicSingleLayer::CourantNumber() const
+{
+    scalar CoNum = 0.0;
+    scalar meanCoNum = 0.0;
+
+    if (filmRegion_.nInternalFaces())
+    {
+        const scalar deltaT = time_.deltaTValue();
+
+        surfaceScalarField SfUfbyDelta =
+            filmRegion_.surfaceInterpolation::deltaCoeffs()*mag(phi_)
+           /fvc::interpolate
+            (
+                rho_*(delta_ + dimensionedScalar("SMALL", dimLength, SMALL))
+            );
+
+        CoNum = max(SfUfbyDelta/filmRegion_.magSf()).value()*deltaT;
+
+        meanCoNum = (sum(SfUfbyDelta)/sum(filmRegion_.magSf())).value()*deltaT;
+    }
+
+    Info<< "    Courant number mean: " << meanCoNum << " max: " << CoNum
+        << endl;
+
+    return CoNum;
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::continuityCheck()
+{
+    const volScalarField deltaRho0 = deltaRho_;
+
+    solveContinuity();
+
+    if (debug)
+    {
+        volScalarField mass = deltaRho_*magSf_;
+        dimensionedScalar totalMass =
+            fvc::domainIntegrate(mass)
+          + dimensionedScalar("SMALL", dimMass*dimVolume, ROOTVSMALL);
+
+        scalar sumLocalContErr =
+            (
+                fvc::domainIntegrate(mag(mass - magSf_*deltaRho0))/totalMass
+            ).value();
+
+        scalar globalContErr =
+            (
+                fvc::domainIntegrate(mass - magSf_*deltaRho0)/totalMass
+            ).value();
+
+        cumulativeContErr_ += globalContErr;
+
+        Info<< "Surface film: " << type() << nl
+            << "    time step continuity errors: sum local = "
+            << sumLocalContErr << ", global = " << globalContErr
+            << ", cumulative = " << cumulativeContErr_ << endl;
+    }
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::solveContinuity()
+{
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::solveContinuity()" << endl;
+    }
+
+    solve
+    (
+        fvm::ddt(deltaRho_)
+      + fvc::div(phi_)
+     ==
+        rhoSp_
+    );
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::updateSurfaceVelocities()
+{
+    // Push boundary film velocity values into internal field
+    for (label i=0; i<filmBottomPatchIDs_.size(); i++)
+    {
+        label patchI = filmBottomPatchIDs_[i];
+        const polyPatch& pp = filmRegion_.boundaryMesh()[patchI];
+        UIndirectList<vector>(Uw_, pp.faceCells()) =
+            U_.boundaryField()[patchI];
+    }
+    Uw_ -= nHat_*(Uw_ & nHat_);
+    Uw_.correctBoundaryConditions();
+
+    // TODO: apply quadratic profile to determine surface velocity
+    Us_ = U_;
+    Us_.correctBoundaryConditions();
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::surfaceFilmModels::kinematicSingleLayer::tau
+(
+    volVectorField& U
+) const
+{
+    // Calculate shear stress
+    volScalarField Cs("Cs", rho_*Cf_*mag(Us_ - U));
+    volScalarField Cw
+    (
+        "Cw",
+        mu_/(0.3333*(delta_ + dimensionedScalar("SMALL", dimLength, SMALL)))
+    );
+    Cw.min(1.0e+06);
+
+    return
+    (
+       - fvm::Sp(Cs, U) + Cs*Us_
+       - fvm::Sp(Cw, U) + Cw*Uw_
+    );
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix>
+Foam::surfaceFilmModels::kinematicSingleLayer::solveMomentum
+(
+    const volScalarField& pu,
+    const volScalarField& pp
+)
+{
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::solveMomentum()" << endl;
+    }
+
+    updateSurfaceVelocities();
+
+    // Momentum
+    tmp<fvVectorMatrix> tUEqn
+    (
+        fvm::ddt(deltaRho_, U_)
+      + fvm::div(phi_, U_)
+     ==
+        USp_
+      + tau(U_)
+      + fvc::grad(sigma_)
+    );
+
+    fvVectorMatrix& UEqn = tUEqn();
+
+    UEqn.relax();
+
+    if (momentumPredictor_)
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+              - fvc::interpolate(delta_)
+              * (
+                    filmRegion_.magSf()
+                  * (
+                        fvc::snGrad(pu, "snGrad(p)")
+                      + fvc::snGrad(pp, "snGrad(p)")*fvc::interpolate(delta_)
+                      + fvc::snGrad(delta_)*fvc::interpolate(pp)
+                    )
+                  - (fvc::interpolate(rho_*gTan()) & filmRegion_.Sf())
+                )
+            )
+        );
+
+        // Remove any patch-normal components of velocity
+        U_ -= nHat_*(nHat_ & U_);
+        U_.correctBoundaryConditions();
+    }
+
+    return tUEqn;
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::solveThickness
+(
+    const volScalarField& pu,
+    const volScalarField& pp,
+    const fvVectorMatrix& UEqn
+)
+{
+    if (debug)
+    {
+        Info<< "kinematicSingleLayer::solveThickness()" << endl;
+    }
+
+    volScalarField rUA = 1.0/UEqn.A();
+    U_ = rUA*UEqn.H();
+
+    surfaceScalarField deltarUAf = fvc::interpolate(delta_*rUA);
+    surfaceScalarField rhof = fvc::interpolate(rho_);
+
+    surfaceScalarField phiAdd
+    (
+        "phiAdd",
+        filmRegion_.magSf()
+      * (
+            fvc::snGrad(pu, "snGrad(p)")
+          + fvc::snGrad(pp, "snGrad(p)")*fvc::interpolate(delta_)
+        )
+      - (fvc::interpolate(rho_*gTan()) & filmRegion_.Sf())
+    );
+    constrainFilmField(phiAdd, 0.0);
+
+    surfaceScalarField phid
+    (
+        "phid",
+        (fvc::interpolate(U_*rho_) & filmRegion_.Sf())
+      - deltarUAf*phiAdd*rhof
+    );
+    constrainFilmField(phid, 0.0);
+
+    surfaceScalarField ddrhorUAppf =
+        fvc::interpolate(delta_)*deltarUAf*rhof*fvc::interpolate(pp);
+//    constrainFilmField(ddrhorUAppf, 0.0);
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr_; nonOrth++)
+    {
+        // Film thickness equation
+        fvScalarMatrix deltaEqn
+        (
+            fvm::ddt(rho_, delta_)
+          + fvm::div(phid, delta_)
+          - fvm::laplacian(ddrhorUAppf, delta_)
+         ==
+            rhoSp_
+        );
+
+        deltaEqn.solve();
+
+        if (nonOrth == nNonOrthCorr_)
+        {
+            phiAdd +=
+                fvc::interpolate(pp)
+              * fvc::snGrad(delta_)
+              * filmRegion_.magSf();
+
+            phi_ == deltaEqn.flux();
+        }
+    }
+
+    // Bound film thickness by a minimum of zero
+    delta_.max(0.0);
+
+    // Update U field
+    U_ -= fvc::reconstruct(deltarUAf*phiAdd);
+
+    // Remove any patch-normal components of velocity
+    U_ -= nHat_*(nHat_ & U_);
+
+    U_.correctBoundaryConditions();
+
+    // Continuity check
+    continuityCheck();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::kinematicSingleLayer::kinematicSingleLayer
+(
+    const word& modelType,
+    const fvMesh& mesh,
+    const dimensionedVector& g
+)
+:
+    surfaceFilmModel(modelType, mesh, g),
+    filmRegion_
+    (
+        IOobject
+        (
+            filmRegionName_,
+            time_.timeName(),
+            time_,
+            IOobject::MUST_READ
+        )
+    ),
+    nHat_
+    (
+        IOobject
+        (
+            "nHat",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        filmRegion_,
+        dimensionedVector("zero", dimless, vector::zero),
+        zeroGradientFvPatchVectorField::typeName
+    ),
+    magSf_
+    (
+        IOobject
+        (
+            "magSf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimArea, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    primaryPatchIDs_(0),
+    filmTopPatchIDs_(0),
+    filmBottomPatchIDs_(0),
+
+    momentumPredictor_
+    (
+        filmRegion_.solutionDict().subDict("PISO").lookup("momentumPredictor")
+    ),
+    nOuterCorr_
+    (
+        readLabel
+        (
+            filmRegion_.solutionDict().subDict("PISO").lookup("nOuterCorr")
+        )
+    ),
+    nCorr_
+    (
+        readLabel(filmRegion_.solutionDict().subDict("PISO").lookup("nCorr"))
+    ),
+    nNonOrthCorr_
+    (
+        readLabel
+        (
+            filmRegion_.solutionDict().subDict("PISO").lookup("nNonOrthCorr")
+        )
+    ),
+    cumulativeContErr_(0.0),
+
+    Cf_(readScalar(coeffs_.lookup("Cf"))),
+    deltaStable_(coeffs_.lookup("deltaStable")),
+
+    initialisedThermo_(false),
+    rho_
+    (
+        IOobject
+        (
+            "rhof",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimDensity, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    mu_
+    (
+        IOobject
+        (
+            "muf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimPressure*dimTime, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    sigma_
+    (
+        IOobject
+        (
+            "sigmaf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimMass/sqr(dimTime), 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+
+    delta_
+    (
+        IOobject
+        (
+            "deltaf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_
+    ),
+    U_
+    (
+        IOobject
+        (
+            "Uf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_
+    ),
+    Us_
+    (
+        IOobject
+        (
+            "Usf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        U_,
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    Uw_
+    (
+        IOobject
+        (
+            "Uwf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        U_,
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    deltaRho_
+    (
+        IOobject
+        (
+            delta_.name() + "*" + rho_.name(),
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", delta_.dimensions()*rho_.dimensions(), 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+
+    phi_
+    (
+        IOobject
+        (
+            "phi",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_,
+        dimLength*dimMass/dimTime
+    ),
+
+    massForPrimary_
+    (
+        IOobject
+        (
+            "massForPrimary",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimMass, 0.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    diametersForPrimary_
+    (
+        IOobject
+        (
+            "diametersForPrimary",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimLength, -1.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    massPhaseChangeForPrimary_
+    (
+        IOobject
+        (
+            "massPhaseChangeForPrimary",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimMass, 0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+
+    USp_
+    (
+        IOobject
+        (
+            "USpf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_
+    ),
+    pSp_
+    (
+        IOobject
+        (
+            "pSpf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_
+    ),
+    rhoSp_
+    (
+        IOobject
+        (
+            "rhoSpf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar("zero", dimMass/dimTime/dimArea, 0.0),
+        pSp_.boundaryField().types()
+    ),
+
+    USpPrimary_
+    (
+        IOobject
+        (
+            USp_.name(), // must have same name as USp_ to enable mapping
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector("zero", USp_.dimensions(), vector::zero)
+    ),
+    pSpPrimary_
+    (
+        IOobject
+        (
+            pSp_.name(), // must have same name as pSp_ to enable mapping
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero", pSp_.dimensions(), 0.0)
+    ),
+    rhoSpPrimary_
+    (
+        IOobject
+        (
+            rhoSp_.name(), // must have same name as rhoSp_ to enable mapping
+            time_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero", rhoSp_.dimensions(), 0.0)
+    ),
+
+    UPrimary_
+    (
+        IOobject
+        (
+            "U", // must have same name as U on primary region to enable mapping
+            time_.timeName(),
+            filmRegion_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_
+    ),
+    pPrimary_
+    (
+        IOobject
+        (
+            "p", // must have same name as p on primary region to enable mapping
+            time_.timeName(),
+            filmRegion_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_
+    ),
+
+    injection_(injectionModel::New(*this, coeffs_)),
+
+    addedMass_(0.0),
+    detachedMass_(0.0)
+{
+    initialise();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::surfaceFilmModels::kinematicSingleLayer::~kinematicSingleLayer()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+inline const Foam::fvMesh&
+Foam::surfaceFilmModels::kinematicSingleLayer::film() const
+{
+    return filmRegion_;
+}
+
+
+const Foam::labelList&
+Foam::surfaceFilmModels::kinematicSingleLayer::filmBottomPatchIDs() const
+{
+    return filmBottomPatchIDs_;
+}
+
+
+const Foam::labelList&
+Foam::surfaceFilmModels::kinematicSingleLayer::filmTopPatchIDs() const
+{
+    return filmTopPatchIDs_;
+}
+
+
+const Foam::labelList&
+Foam::surfaceFilmModels::kinematicSingleLayer::primaryPatchIDs() const
+{
+    return primaryPatchIDs_;
+}
+
+
+bool Foam::surfaceFilmModels::kinematicSingleLayer::isFilmPatch
+(
+    const label patchI
+) const
+{
+    if (!active_)
+    {
+        return false;
+    }
+
+    forAll(primaryPatchIDs_, i)
+    {
+        if (primaryPatchIDs_[i] == patchI)
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::addSources
+(
+    const label patchI,
+    const label faceI,
+    const scalar massSource,
+    const vector& momentumSource,
+    const scalar pressureSource,
+    const scalar energySource
+)
+{
+    if (debug)
+    {
+        Info<< "\nSurface film: " << type() << ": adding to film source:" << nl
+            << "    mass     = " << massSource << nl
+            << "    momentum = " << momentumSource << nl
+            << "    pressure = " << pressureSource << endl;
+    }
+
+    rhoSpPrimary_.boundaryField()[patchI][faceI] += massSource;
+    USpPrimary_.boundaryField()[patchI][faceI] += momentumSource;
+    pSpPrimary_.boundaryField()[patchI][faceI] += pressureSource;
+
+    addedMass_ += massSource;
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::preEvolveFilm()
+{
+    if (!initialisedThermo_)
+    {
+        correctThermoFields();
+
+        deltaRho_ == delta_*rho_;
+        phi_ = fvc::interpolate(deltaRho_*U_) & filmRegion_.Sf();
+        initialisedThermo_ = true;
+    }
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::evolveFilm()
+{
+    transferPrimaryRegionFields();
+
+    updateSubmodels();
+
+    // Solve continuity for deltaRho_
+    solveContinuity();
+
+    for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
+    {
+        // Explicit pressure source contribution
+        tmp<volScalarField> tpu = this->pu();
+
+        // Implicit pressure source coefficient
+        tmp<volScalarField> tpp = this->pp();
+
+        // Solve for momentum for U_
+        tmp<fvVectorMatrix> UEqn = solveMomentum(tpu(), tpp());
+
+        // Film thickness correction loop
+        for (int corr=1; corr<=nCorr_; corr++)
+        {
+            // Solve thickness for delta_
+            solveThickness(tpu(), tpp(), UEqn());
+        }
+    }
+
+    // Update deltaRho_ with new delta_
+    deltaRho_ == delta_*rho_;
+
+    // Update film wall and surface velocities
+    updateSurfaceVelocities();
+
+    // Reset source terms for next time integration
+    resetPrimaryRegionSourceTerms();
+}
+
+
+const Foam::volVectorField&
+Foam::surfaceFilmModels::kinematicSingleLayer::U() const
+{
+    return U_;
+}
+
+
+const Foam::volVectorField&
+Foam::surfaceFilmModels::kinematicSingleLayer::Us() const
+{
+    return Us_;
+}
+
+
+const Foam::volVectorField&
+Foam::surfaceFilmModels::kinematicSingleLayer::Uw() const
+{
+    return Uw_;
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::rho() const
+{
+    return rho_;
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::T() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& kinematicSingleLayer::T() const"
+    )   << "T field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::Ts() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& kinematicSingleLayer::Ts() const"
+    )   << "Ts field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::Tw() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& kinematicSingleLayer::Tw() const"
+    )   << "Tw field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::cp() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& kinematicSingleLayer::cp() const"
+    )   << "cp field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::kappa() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& kinematicSingleLayer::kappa() const"
+    )   << "kappa field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::massForPrimary() const
+{
+    return massForPrimary_;
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::diametersForPrimary() const
+{
+    return diametersForPrimary_;
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::massPhaseChangeForPrimary() const
+{
+    return massPhaseChangeForPrimary_;
+}
+
+
+void Foam::surfaceFilmModels::kinematicSingleLayer::info() const
+{
+    Info<< "\nSurface film: " << type() << endl;
+
+    // Output Courant number for info only - does not change time step
+    CourantNumber();
+
+    Info<< indent << "added mass        = "
+        << returnReduce<scalar>(addedMass_, sumOp<scalar>()) << nl
+        << indent << "current mass      = "
+        << gSum((deltaRho_*magSf_)()) << nl
+        << indent << "detached mass     = "
+        << returnReduce<scalar>(detachedMass_, sumOp<scalar>()) << nl
+        << indent << "min/max(mag(U))   = " << min(mag(U_)).value() << ", "
+        << max(mag(U_)).value() << nl
+        << indent << "min/max(delta)    = " << min(delta_).value() << ", "
+        << max(delta_).value() << nl;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::surfaceFilmModels::kinematicSingleLayer::Srho() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tSrho
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "kinematicSingleLayer::Srho",
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
+
+    scalarField& Srho = tSrho();
+    const scalarField& V = mesh_.V();
+    const scalar dt = time_.deltaTValue();
+
+    forAll(filmBottomPatchIDs_, i)
+    {
+        const label primaryPatchI = primaryPatchIDs_[i];
+        const directMappedWallPolyPatch& wpp =
+            refCast<const directMappedWallPolyPatch>
+            (
+                 mesh_.boundaryMesh()[primaryPatchI]
+            );
+
+        const mapDistribute& distMap = wpp.map();
+
+        const label filmPatchI = filmBottomPatchIDs_[i];
+
+        scalarField patchMass =
+            massPhaseChangeForPrimary().boundaryField()[filmPatchI];
+
+        distMap.distribute(patchMass);
+
+        const unallocLabelList& cells = wpp.faceCells();
+
+        forAll(patchMass, j)
+        {
+            Srho[cells[j]] = patchMass[j]/(V[cells[j]]*dt);
+        }
+    }
+
+    return tSrho;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::surfaceFilmModels::kinematicSingleLayer::Srho(const label) const
+{
+    return tmp<DimensionedField<scalar, volMesh> >
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "kinematicSingleLayer::Srho(i)",
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::surfaceFilmModels::kinematicSingleLayer::Sh() const
+{
+    return tmp<DimensionedField<scalar, volMesh> >
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "kinematicSingleLayer::Sh",
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+        )
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.H b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.H
index 095b713b59f..1e45c3defff 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.H
+++ b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayer.H
@@ -1,482 +1,558 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::kinematicSingleLayer
-
-Description
-    Kinematic form of single-cell layer surface film model
-
-SourceFiles
-    kinematicSingleLayer.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef kinematicSingleLayer_H
-#define kinematicSingleLayer_H
-
-#include "surfaceFilmModel.H"
-#include "fvMesh.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-#include "fvMatrices.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace surfaceFilmModels
-{
-
-// Forward declaration of classes
-class injectionModel;
-
-/*---------------------------------------------------------------------------*\
-                      Class kinematicSingleLayer Declaration
-\*---------------------------------------------------------------------------*/
-
-class kinematicSingleLayer
-:
-    public surfaceFilmModel
-{
-private:
-
-    // Private member functions
-
-        //- Disallow default bitwise copy construct
-        kinematicSingleLayer(const kinematicSingleLayer&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const kinematicSingleLayer&);
-
-
-protected:
-
-    // Protected data
-
-        // Mesh databases
-
-            //- Film region mesh database
-            fvMesh filmRegion_;
-
-            //- Patch normal vectors
-            volVectorField nHat_;
-
-            //- Face area magnitudes / [m2]
-            volScalarField magSf_;
-
-            //- List of patch IDs on the primary region coupled with the film
-            //  region
-            labelList primaryPatchIDs_;
-
-            //- List of patch IDs on oppositte side of the film region
-            labelList filmTopPatchIDs_;
-
-            //- List of patch IDs on the film region coupled with the primary
-            //  region
-            labelList filmBottomPatchIDs_;
-
-
-        // Solution parameters
-
-            //- Momentum predictor
-            Switch momentumPredictor_;
-
-            //- Number of outer correctors
-            label nOuterCorr_;
-
-            //- Number of PISO-like correctors
-            label nCorr_;
-
-            //- Number of non-orthogonal correctors
-            label nNonOrthCorr_;
-
-            //- Cumulative continuity error
-            scalar cumulativeContErr_;
-
-
-        // Model parameters
-
-            //- Skin frition coefficient for film/main region interface
-            scalar Cf_;
-
-            //- Stable film thickness - film cannot detach until this is reached
-            dimensionedScalar deltaStable_;
-
-
-        // Thermo properties
-
-            //- Density / [kg/m3]
-            volScalarField rho_;
-
-            //- Dynamic viscosity / [Pa.s]
-            volScalarField mu_;
-
-            //- Surface tension / [m/s2]
-            volScalarField sigma_;
-
-
-        // Fields
-
-            //- Film thickness / [m]
-            volScalarField delta_;
-
-            //- Velocity / [m/s]
-            volVectorField U_;
-
-            //- Film thickness*density (helper field) / [kg/m2]
-            volScalarField deltaRho_;
-
-            //- Mass flux (includes film thickness) / [kg.m/s]
-            surfaceScalarField phi_;
-
-
-            // Transfer fields - to the primary region
-
-                //- Return the film mass available for transfer
-                volScalarField massForPrimary_;
-
-                //- Return the parcel diameters originating from film
-                volScalarField diametersForPrimary_;
-
-
-        // Source term fields
-
-            // Film region - registered to the film region mesh
-            // Note: need boundary value mapped from primary region, and then
-            // pushed into the patch internal field
-
-                //- Momementum / [kg/m/s2]
-                volVectorField USp_;
-
-                //- Pressure / [Pa]
-                volScalarField pSp_;
-
-                //- Mass / [kg/m2/s]
-                volScalarField rhoSp_;
-
-
-            // Primary region - registered to the primary region mesh
-            // Internal use only - not read-in
-
-                //- Momementum / [kg/m/s2]
-                volVectorField USpPrimary_;
-
-                //- Pressure / [Pa]
-                volScalarField pSpPrimary_;
-
-                //- Mass / [kg/m2/s]
-                volScalarField rhoSpPrimary_;
-
-
-        // Fields mapped from primary region - registered to the film region
-        // Note: need both boundary AND patch internal fields to be mapped
-
-            //- Velocity / [m/s]
-            volVectorField UPrimary_;
-
-            //- Pressure / [Pa]
-            volScalarField pPrimary_;
-
-
-        // Sub-models
-
-            //- Injection
-            autoPtr<injectionModel> injection_;
-
-
-       // Checks
-
-           //- Cumulative mass added via sources [kg]
-           scalar addedMass_;
-
-
-       // Detached surface properties
-
-           //- Cumulative mass detached [kg]
-           scalar detachedMass_;
-
-
-    // Protected member functions
-
-        //- Initialise the film model - called on construction
-        void initialise();
-
-        //- Read control parameters from dictionary
-        virtual bool read();
-
-        //- Reset source term fields
-        virtual void resetPrimaryRegionSourceTerms();
-
-        //- Transfer fields from the primary region to the film region
-        virtual void transferPrimaryRegionFields();
-
-        //- Correct the source terms for film that detaches from film region
-        virtual void correctDetachedFilm();
-
-        // Explicit pressure source contribution
-        virtual tmp<volScalarField> pu();
-
-        // Implicit pressure source coefficient
-        virtual tmp<volScalarField> pp();
-
-        //- Update the film sub-models
-        virtual void updateSubmodels();
-
-        //- Courant number evaluation
-        virtual scalar CourantNumber() const;
-
-        //- Continuity check
-        virtual void continuityCheck();
-
-        //- Return the stress term for the momentum equation
-        virtual tmp<fvVectorMatrix> tau(volVectorField& dU) const;
-
-        //- Constrain a film region master/slave boundaries of a field to a
-        //  given value
-        template<class Type>
-        void constrainFilmField
-        (
-            Type& field,
-            const typename Type::cmptType& value
-        );
-
-
-        // Equations
-
-            //- Solve continuity equation
-            virtual void solveContinuity();
-
-            //- Solve for film velocity
-            virtual tmp<fvVectorMatrix> solveMomentum
-            (
-                const volScalarField& pu,
-                const volScalarField& pp
-            );
-
-            //- Solve coupled velocity-thickness equations
-            virtual void solveThickness
-            (
-                const volScalarField& pu,
-                const volScalarField& pp,
-                const fvVectorMatrix& UEqn
-            );
-
-
-public:
-
-    //- Runtime type information
-    TypeName("kinematicSingleLayer");
-
-
-    // Constructors
-
-        //- Construct from components
-        kinematicSingleLayer
-        (
-            const word& modelType,
-            const fvMesh& mesh,
-            const dimensionedVector& g
-        );
-
-
-    //- Destructor
-    virtual ~kinematicSingleLayer();
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return the film mesh database
-            virtual const fvMesh& film() const;
-
-            //- Return the patch normal vectors
-            inline const volVectorField& nHat() const;
-
-            //- Return the face area magnitudes / [m2]
-            inline const volScalarField& magSf() const;
-
-            //- Return the list of coupled patches on the film region
-            virtual const labelList& filmBottomPatchIDs() const;
-
-            //- Return the list of patches oppositte coupled patches
-            virtual const labelList& filmTopPatchIDs() const;
-
-            //- Return the list of coupled patches on the primary region
-            virtual const labelList& primaryPatchIDs() const;
-
-
-        // Solution parameters
-
-            //- Return the momentum predictor
-            inline const Switch& momentumPredictor() const;
-
-            //- Return the number of outer correctors
-            inline label nOuterCorr() const;
-
-            //- Return the number of PISO correctors
-            inline label nCorr() const;
-
-            //- Return the number of non-orthogonal correctors
-            inline label nNonOrthCorr() const;
-
-
-        // Model parameters
-
-            //- Return the skin friction coefficient
-            inline scalar Cf() const;
-
-
-        // Thermo properties
-
-            //- Return const access to the dynamic viscosity / [Pa.s]
-            inline const volScalarField& mu() const;
-
-            //- Return const access to the surface tension / [m/s2]
-            inline const volScalarField& sigma() const;
-
-
-        // Fields
-
-            //- Return const access to the film thickness / [m]
-            inline const volScalarField& delta() const;
-
-            //- Return the film velocity [m/s]
-            virtual const volVectorField& U() const;
-
-            //- Return the film density [kg/m3]
-            virtual const volScalarField& rho() const;
-
-            //- Return the film temperature [K]
-            virtual const volScalarField& T() const;
-
-            //- Return the film specific heat capacity [J/kg/K]
-            virtual const volScalarField& cp() const;
-
-
-           // Transfer fields - to the primary region
-
-               //- Return the film mass available for transfer
-               virtual const volScalarField& massForPrimary() const;
-
-               //- Return the parcel diameters originating from film
-               virtual const volScalarField& diametersForPrimary() const;
-
-
-        // External helper functions
-
-            //- Return true if patchI on the primary region is a coupled patch
-            //  to the film region
-            virtual bool isFilmPatch(const label patchI) const;
-
-            //- External hook to add sources to the film
-            virtual void addSources
-            (
-                const label patchI,            // patchI on primary region
-                const label faceI,             // faceI of patchI
-                const scalar massSource,       // [kg]
-                const vector& momentumSource,  // [kg.m/s] (tangential momentum)
-                const scalar pressureSource,   // [kg.m/s] (normal momentum)
-                const scalar energySource = 0  // [J]
-            );
-
-
-         // Source fields (read/write access)
-
-            // Primary region
-
-                //- Momementum / [kg/m/s2]
-                inline volVectorField& USpPrimary();
-
-                //- Pressure / [Pa]
-                inline volScalarField& pSpPrimary();
-
-                //- Mass / [kg/m2/s]
-                inline volScalarField& rhoSpPrimary();
-
-
-            // Film region
-
-                //- Momementum / [kg/m/s2]
-                inline volVectorField& USp();
-
-                //- Pressure / [Pa]
-                inline volScalarField& pSp();
-
-                //- Mass / [kg/m2/s]
-                inline volScalarField& rhoSp();
-
-
-        // Sub-models
-
-            //- Injection
-            inline injectionModel& injection();
-
-
-        // Helper functions
-
-            //- Return the gravity tangential component contributions
-            inline tmp<volVectorField> gTan() const;
-
-            //- Return the gravity normal-to-patch component contribution
-            inline tmp<volScalarField> gNorm() const;
-
-            //- Return the gravity normal-to-patch component contribution
-            //  Clipped so that only non-zero if g & nHat_ < 0
-            inline tmp<volScalarField> gNormClipped() const;
-
-
-       // Evolution
-
-            //- Evolve the film equations
-            virtual void evolveFilm();
-
-
-       // I-O
-
-            //- Provide some feedback
-            virtual void info() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace surfaceFilmModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "kinematicSingleLayerTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "kinematicSingleLayerI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::kinematicSingleLayer
+
+Description
+    Kinematic form of single-cell layer surface film model
+
+SourceFiles
+    kinematicSingleLayer.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kinematicSingleLayer_H
+#define kinematicSingleLayer_H
+
+#include "surfaceFilmModel.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFilmModels
+{
+
+// Forward declaration of classes
+class injectionModel;
+
+/*---------------------------------------------------------------------------*\
+                   Class kinematicSingleLayer Declaration
+\*---------------------------------------------------------------------------*/
+
+class kinematicSingleLayer
+:
+    public surfaceFilmModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        kinematicSingleLayer(const kinematicSingleLayer&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const kinematicSingleLayer&);
+
+
+protected:
+
+    // Protected data
+
+        // Mesh databases
+
+            //- Film region mesh database
+            fvMesh filmRegion_;
+
+            //- Patch normal vectors
+            volVectorField nHat_;
+
+            //- Face area magnitudes / [m2]
+            volScalarField magSf_;
+
+            //- List of patch IDs on the primary region coupled with the film
+            //  region
+            labelList primaryPatchIDs_;
+
+            //- List of patch IDs on oppositte side of the film region
+            labelList filmTopPatchIDs_;
+
+            //- List of patch IDs on the film region coupled with the primary
+            //  region
+            labelList filmBottomPatchIDs_;
+
+
+        // Solution parameters
+
+            //- Momentum predictor
+            Switch momentumPredictor_;
+
+            //- Number of outer correctors
+            label nOuterCorr_;
+
+            //- Number of PISO-like correctors
+            label nCorr_;
+
+            //- Number of non-orthogonal correctors
+            label nNonOrthCorr_;
+
+            //- Cumulative continuity error
+            scalar cumulativeContErr_;
+
+
+        // Model parameters
+
+            //- Skin frition coefficient for film/primary region interface
+            scalar Cf_;
+
+            //- Stable film thickness - film cannot detach until this is reached
+            dimensionedScalar deltaStable_;
+
+
+        // Thermo properties
+
+            // Fields
+
+                //- Initiliased thermo flag
+                bool initialisedThermo_;
+
+                //- Density / [kg/m3]
+                volScalarField rho_;
+
+                //- Dynamic viscosity / [Pa.s]
+                volScalarField mu_;
+
+                //- Surface tension / [m/s2]
+                volScalarField sigma_;
+
+
+        // Fields
+
+            //- Film thickness / [m]
+            volScalarField delta_;
+
+            //- Velocity - mean / [m/s]
+            volVectorField U_;
+
+            //- Velocity - surface / [m/s]
+            volVectorField Us_;
+
+            //- Velocity - wall / [m/s]
+            volVectorField Uw_;
+
+            //- Film thickness*density (helper field) / [kg/m2]
+            volScalarField deltaRho_;
+
+            //- Mass flux (includes film thickness) / [kg.m/s]
+            surfaceScalarField phi_;
+
+
+            // Transfer fields - to the primary region
+
+                //- Film mass available for transfer
+                volScalarField massForPrimary_;
+
+                //- Parcel diameters originating from film
+                volScalarField diametersForPrimary_;
+
+                //- Film mass evolved via phase change
+                volScalarField massPhaseChangeForPrimary_;
+
+
+        // Source term fields
+
+            // Film region - registered to the film region mesh
+            // Note: need boundary value mapped from primary region, and then
+            // pushed into the patch internal field
+
+                //- Momementum / [kg/m/s2]
+                volVectorField USp_;
+
+                //- Pressure / [Pa]
+                volScalarField pSp_;
+
+                //- Mass / [kg/m2/s]
+                volScalarField rhoSp_;
+
+
+            // Primary region - registered to the primary region mesh
+            // Internal use only - not read-in
+
+                //- Momementum / [kg/m/s2]
+                volVectorField USpPrimary_;
+
+                //- Pressure / [Pa]
+                volScalarField pSpPrimary_;
+
+                //- Mass / [kg/m2/s]
+                volScalarField rhoSpPrimary_;
+
+
+        // Fields mapped from primary region - registered to the film region
+        // Note: need both boundary AND patch internal fields to be mapped
+
+            //- Velocity / [m/s]
+            volVectorField UPrimary_;
+
+            //- Pressure / [Pa]
+            volScalarField pPrimary_;
+
+
+        // Sub-models
+
+            //- Injection
+            autoPtr<injectionModel> injection_;
+
+
+       // Checks
+
+           //- Cumulative mass added via sources [kg]
+           scalar addedMass_;
+
+
+       // Detached surface properties
+
+           //- Cumulative mass detached [kg]
+           scalar detachedMass_;
+
+
+    // Protected member functions
+
+        //- Initialise the film model - called on construction
+        void initialise();
+
+        //- Read control parameters from dictionary
+        virtual bool read();
+
+        //- Correct the thermo fields
+        virtual void correctThermoFields();
+
+        //- Reset source term fields
+        virtual void resetPrimaryRegionSourceTerms();
+
+        //- Transfer fields from the primary region to the film region
+        virtual void transferPrimaryRegionFields();
+
+        //- Correct the source terms for film that detaches from film region
+        virtual void correctDetachedFilm();
+
+        // Explicit pressure source contribution
+        virtual tmp<volScalarField> pu();
+
+        // Implicit pressure source coefficient
+        virtual tmp<volScalarField> pp();
+
+        //- Update the film sub-models
+        virtual void updateSubmodels();
+
+        //- Courant number evaluation
+        virtual scalar CourantNumber() const;
+
+        //- Continuity check
+        virtual void continuityCheck();
+
+        //- Update film surface velocities
+        virtual void updateSurfaceVelocities();
+
+        //- Return the stress term for the momentum equation
+        virtual tmp<fvVectorMatrix> tau(volVectorField& dU) const;
+
+        //- Constrain a film region master/slave boundaries of a field to a
+        //  given value
+        template<class Type>
+        void constrainFilmField
+        (
+            Type& field,
+            const typename Type::cmptType& value
+        );
+
+
+        // Equations
+
+            //- Solve continuity equation
+            virtual void solveContinuity();
+
+            //- Solve for film velocity
+            virtual tmp<fvVectorMatrix> solveMomentum
+            (
+                const volScalarField& pu,
+                const volScalarField& pp
+            );
+
+            //- Solve coupled velocity-thickness equations
+            virtual void solveThickness
+            (
+                const volScalarField& pu,
+                const volScalarField& pp,
+                const fvVectorMatrix& UEqn
+            );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("kinematicSingleLayer");
+
+
+    // Constructors
+
+        //- Construct from components
+        kinematicSingleLayer
+        (
+            const word& modelType,
+            const fvMesh& mesh,
+            const dimensionedVector& g
+        );
+
+
+    //- Destructor
+    virtual ~kinematicSingleLayer();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the film mesh database
+            virtual const fvMesh& film() const;
+
+            //- Return the patch normal vectors
+            inline const volVectorField& nHat() const;
+
+            //- Return the face area magnitudes / [m2]
+            inline const volScalarField& magSf() const;
+
+            //- Return the list of coupled patches on the film region
+            virtual const labelList& filmBottomPatchIDs() const;
+
+            //- Return the list of patches oppositte coupled patches
+            virtual const labelList& filmTopPatchIDs() const;
+
+            //- Return the list of coupled patches on the primary region
+            virtual const labelList& primaryPatchIDs() const;
+
+
+        // Solution parameters
+
+            //- Return the momentum predictor
+            inline const Switch& momentumPredictor() const;
+
+            //- Return the number of outer correctors
+            inline label nOuterCorr() const;
+
+            //- Return the number of PISO correctors
+            inline label nCorr() const;
+
+            //- Return the number of non-orthogonal correctors
+            inline label nNonOrthCorr() const;
+
+
+        // Model parameters
+
+            //- Return the skin friction coefficient
+            inline scalar Cf() const;
+
+
+        // Thermo properties
+
+            //- Return const access to the dynamic viscosity / [Pa.s]
+            inline const volScalarField& mu() const;
+
+            //- Return const access to the surface tension / [m/s2]
+            inline const volScalarField& sigma() const;
+
+
+        // Fields
+
+            //- Return const access to the film thickness / [m]
+            inline const volScalarField& delta() const;
+
+            //- Return the film velocity [m/s]
+            virtual const volVectorField& U() const;
+
+            //- Return the film surface velocity [m/s]
+            virtual const volVectorField& Us() const;
+
+            //- Return the film wall velocity [m/s]
+            virtual const volVectorField& Uw() const;
+
+            //- Return the film density [kg/m3]
+            virtual const volScalarField& rho() const;
+
+            //- Return the film mean temperature [K]
+            virtual const volScalarField& T() const;
+
+            //- Return the film surface temperature [K]
+            virtual const volScalarField& Ts() const;
+
+            //- Return the film wall temperature [K]
+            virtual const volScalarField& Tw() const;
+
+            //- Return the film specific heat capacity [J/kg/K]
+            virtual const volScalarField& cp() const;
+
+            //- Return the film thermal conductivity [W/m/K]
+            virtual const volScalarField& kappa() const;
+
+
+           // Transfer fields - to the primary region
+
+               //- Return the film mass available for transfer
+               virtual const volScalarField& massForPrimary() const;
+
+               //- Return the parcel diameters originating from film
+               virtual const volScalarField& diametersForPrimary() const;
+
+               //- Return the film mass evolved via phase change
+               virtual const volScalarField& massPhaseChangeForPrimary() const;
+
+
+        // External helper functions
+
+            //- Return true if patchI on the primary region is a coupled patch
+            //  to the film region
+            virtual bool isFilmPatch(const label patchI) const;
+
+            //- External hook to add sources to the film
+            virtual void addSources
+            (
+                const label patchI,            // patchI on primary region
+                const label faceI,             // faceI of patchI
+                const scalar massSource,       // [kg]
+                const vector& momentumSource,  // [kg.m/s] (tangential momentum)
+                const scalar pressureSource,   // [kg.m/s] (normal momentum)
+                const scalar energySource = 0  // [J]
+            );
+
+
+         // Source fields (read/write access)
+
+            // Primary region
+
+                //- Momementum / [kg/m/s2]
+                inline volVectorField& USpPrimary();
+
+                //- Pressure / [Pa]
+                inline volScalarField& pSpPrimary();
+
+                //- Mass / [kg/m2/s]
+                inline volScalarField& rhoSpPrimary();
+
+
+            // Film region
+
+                //- Momentum / [kg/m/s2]
+                inline volVectorField& USp();
+
+                //- Pressure / [Pa]
+                inline volScalarField& pSp();
+
+                //- Mass / [kg/m2/s]
+                inline volScalarField& rhoSp();
+
+                //- Momentum / [kg/m/s2]
+                inline const volVectorField& USp() const;
+
+                //- Pressure / [Pa]
+                inline const volScalarField& pSp() const;
+
+                //- Mass / [kg/m2/s]
+                inline const volScalarField& rhoSp() const;
+
+
+        // Fields mapped from primary region
+
+            //- Velocity / [m/s]
+            inline const volVectorField& UPrimary() const;
+
+            //- Pressure / [Pa]
+            inline const volScalarField& pPrimary() const;
+
+
+        // Sub-models
+
+            //- Injection
+            inline injectionModel& injection();
+
+
+        // Helper functions
+
+            //- Return the gravity tangential component contributions
+            inline tmp<volVectorField> gTan() const;
+
+            //- Return the gravity normal-to-patch component contribution
+            inline tmp<volScalarField> gNorm() const;
+
+            //- Return the gravity normal-to-patch component contribution
+            //  Clipped so that only non-zero if g & nHat_ < 0
+            inline tmp<volScalarField> gNormClipped() const;
+
+
+        // Evolution
+
+            //- Pre-evolve film hook
+            virtual void preEvolveFilm();
+
+            //- Evolve the film equations
+            virtual void evolveFilm();
+
+
+        // Source fields
+
+            // Mapped into primary region
+
+                //- Return total mass source - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Srho() const;
+
+                //- Return mass source for specie i - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Srho
+                (
+                    const label i
+                ) const;
+
+                //- Return enthalpy source - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Sh() const;
+
+
+       // I-O
+
+            //- Provide some feedback
+            virtual void info() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "kinematicSingleLayerTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "kinematicSingleLayerI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayerI.H b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayerI.H
index 521d460d44f..c0f4a434678 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayerI.H
+++ b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayerI.H
@@ -136,6 +136,41 @@ Foam::surfaceFilmModels::kinematicSingleLayer::rhoSp()
 }
 
 
+inline const Foam::volVectorField&
+Foam::surfaceFilmModels::kinematicSingleLayer::USp() const
+{
+    return USp_;
+}
+
+
+inline const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::pSp() const
+{
+    return pSp_;
+}
+
+
+inline const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::rhoSp() const
+{
+    return rhoSp_;
+}
+
+
+inline const Foam::volVectorField&
+Foam::surfaceFilmModels::kinematicSingleLayer::UPrimary() const
+{
+    return UPrimary_;
+}
+
+
+inline const Foam::volScalarField&
+Foam::surfaceFilmModels::kinematicSingleLayer::pPrimary() const
+{
+    return pPrimary_;
+}
+
+
 inline Foam::surfaceFilmModels::injectionModel&
 Foam::surfaceFilmModels::kinematicSingleLayer::injection()
 {
diff --git a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayerTemplates.C b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayerTemplates.C
index 91117356485..de8a93b2708 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayerTemplates.C
+++ b/src/surfaceFilmModels/surfaceFilmModel/kinematicSingleLayer/kinematicSingleLayerTemplates.C
@@ -1,62 +1,62 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "kinematicSingleLayer.H"
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class Type>
-void Foam::surfaceFilmModels::kinematicSingleLayer::constrainFilmField
-(
-    Type& field,
-    const typename Type::cmptType& value
-)
-{
-    forAll(filmBottomPatchIDs_, i)
-    {
-        label patchI = filmBottomPatchIDs_[i];
-        field.boundaryField()[patchI] = value;
-        if (debug)
-        {
-            Info<< "Constraining " << field.name()
-                << " boundary " << field.boundaryField()[patchI].patch().name()
-                << " to " << value << endl;
-        }
-    }
-    forAll(filmTopPatchIDs_, i)
-    {
-        label patchI = filmTopPatchIDs_[i];
-        field.boundaryField()[patchI] = value;
-        if (debug)
-        {
-            Info<< "Constraining " << field.name()
-                << " boundary " << field.boundaryField()[patchI].patch().name()
-                << " to " << value << endl;
-        }
-    }
-}
-
-
-// ************************************************************************* //
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kinematicSingleLayer.H"
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::surfaceFilmModels::kinematicSingleLayer::constrainFilmField
+(
+    Type& field,
+    const typename Type::cmptType& value
+)
+{
+    forAll(filmBottomPatchIDs_, i)
+    {
+        label patchI = filmBottomPatchIDs_[i];
+        field.boundaryField()[patchI] = value;
+        if (debug)
+        {
+            Info<< "Constraining " << field.name()
+                << " boundary " << field.boundaryField()[patchI].patch().name()
+                << " to " << value << endl;
+        }
+    }
+    forAll(filmTopPatchIDs_, i)
+    {
+        label patchI = filmTopPatchIDs_[i];
+        field.boundaryField()[patchI] = value;
+        if (debug)
+        {
+            Info<< "Constraining " << field.name()
+                << " boundary " << field.boundaryField()[patchI].patch().name()
+                << " to " << value << endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C
index d8fd74960f9..3565cf7c125 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C
+++ b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.C
@@ -76,6 +76,12 @@ Foam::surfaceFilmModels::noFilm::~noFilm()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+void Foam::surfaceFilmModels::noFilm::preEvolveFilm()
+{
+    // do nothing
+}
+
+
 void Foam::surfaceFilmModels::noFilm::evolveFilm()
 {
     // do nothing
@@ -87,7 +93,7 @@ const Foam::fvMesh& Foam::surfaceFilmModels::noFilm::film() const
     FatalErrorIn("const fvMesh& noFilm::film() const")
         << "Cannot return film for noFilm model" << abort(FatalError);
 
-    return reinterpret_cast<const fvMesh&>(null);
+    return mesh();
 }
 
 
@@ -141,6 +147,28 @@ const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::U() const
 }
 
 
+const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Us() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& noFilm::Us() const"
+    )   << "Us field not available for " << type() << abort(FatalError);
+
+    return volVectorField::null();
+}
+
+
+const Foam::volVectorField& Foam::surfaceFilmModels::noFilm::Uw() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& noFilm::Uw() const"
+    )   << "Uw field not available for " << type() << abort(FatalError);
+
+    return volVectorField::null();
+}
+
+
 const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::rho() const
 {
     FatalErrorIn
@@ -163,6 +191,28 @@ const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::T() const
 }
 
 
+const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::Ts() const
+{
+    FatalErrorIn
+    (
+        "const Foam::volScalarField& Foam::noFilm::Ts() const"
+    )   << "Ts field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
+const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::Tw() const
+{
+    FatalErrorIn
+    (
+        "const Foam::volScalarField& Foam::noFilm::Tw() const"
+    )   << "Tw field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
 const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::cp() const
 {
     FatalErrorIn
@@ -174,6 +224,17 @@ const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::cp() const
 }
 
 
+const Foam::volScalarField& Foam::surfaceFilmModels::noFilm::kappa() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& noFilm::kappa() const"
+    )   << "kappa field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
 const Foam::volScalarField&
 Foam::surfaceFilmModels::noFilm::massForPrimary() const
 {
@@ -200,6 +261,19 @@ Foam::surfaceFilmModels::noFilm::diametersForPrimary() const
 }
 
 
+const Foam::volScalarField&
+Foam::surfaceFilmModels::noFilm::massPhaseChangeForPrimary() const
+{
+    FatalErrorIn
+    (
+        "const volScalarField& noFilm::massPhaseChangeForPrimary() const"
+    )   << "massPhaseChangeForPrimary field not available for " << type()
+        << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
 void Foam::surfaceFilmModels::noFilm::info() const
 {
     // do nothing
diff --git a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.H b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.H
index adc97166bbd..7d5fb40021b 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.H
+++ b/src/surfaceFilmModels/surfaceFilmModel/noFilm/noFilm.H
@@ -132,12 +132,27 @@ public:
            //- Return the film density [kg/m3]
            virtual const volScalarField& rho() const;
 
-           //- Return the film temperature [K]
+           //- Return the film surface velocity [m/s]
+           virtual const volVectorField& Us() const;
+
+           //- Return the film wall velocity [m/s]
+           virtual const volVectorField& Uw() const;
+
+           //- Return the film mean temperature [K]
            virtual const volScalarField& T() const;
 
+           //- Return the film surface temperature [K]
+           virtual const volScalarField& Ts() const;
+
+           //- Return the film wall temperature [K]
+           virtual const volScalarField& Tw() const;
+
            //- Return the film specific heat capacity [J/kg/K]
            virtual const volScalarField& cp() const;
 
+           //- Return the film thermal conductivity [W/m/K]
+           virtual const volScalarField& kappa() const;
+
 
            // Transfer fields - to the primary region
 
@@ -147,8 +162,14 @@ public:
                //- Return the parcel diameters originating from film
                virtual const volScalarField& diametersForPrimary() const;
 
+               //- Return the film mass evolved via phase change
+               virtual const volScalarField& massPhaseChangeForPrimary() const;
+
+
+        // Evolution
 
-       // Evolution
+            //- Pre-evolve film hook
+            virtual void preEvolveFilm();
 
             //- Evolve the film equations
             virtual void evolveFilm();
diff --git a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.C b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.C
index 5dfb27e1a0f..68403073845 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.C
+++ b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.C
@@ -38,6 +38,22 @@ namespace Foam
     }
 }
 
+
+template<>
+const char*
+Foam::NamedEnum<Foam::surfaceFilmModels::surfaceFilmModel::thermoModelType, 2>::
+names[] =
+{
+    "constant",
+    "singleComponent"
+};
+
+
+const
+Foam::NamedEnum<Foam::surfaceFilmModels::surfaceFilmModel::thermoModelType, 2>
+    Foam::surfaceFilmModels::surfaceFilmModel::thermoModelTypeNames_;
+
+
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 bool Foam::surfaceFilmModels::surfaceFilmModel::read()
@@ -74,7 +90,7 @@ Foam::surfaceFilmModels::surfaceFilmModel::surfaceFilmModel
             "surfaceFilmProperties",
             mesh.time().constant(),
             mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
     ),
@@ -83,7 +99,8 @@ Foam::surfaceFilmModels::surfaceFilmModel::surfaceFilmModel
     active_(false),
     g_(g),
     filmRegionName_("none"),
-    coeffs_(dictionary::null)
+    coeffs_(dictionary::null),
+    thermoModel_(tmConstant)
 {}
 
 
@@ -101,7 +118,7 @@ Foam::surfaceFilmModels::surfaceFilmModel::surfaceFilmModel
             "surfaceFilmProperties",
             mesh.time().constant(),
             mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::MUST_READ,
             IOobject::NO_WRITE
         )
     ),
@@ -111,7 +128,8 @@ Foam::surfaceFilmModels::surfaceFilmModel::surfaceFilmModel
     active_(lookup("active")),
     g_(g),
     filmRegionName_(lookup("filmRegionName")),
-    coeffs_(subDict(type + "Coeffs"))
+    coeffs_(subDict(type + "Coeffs")),
+    thermoModel_(thermoModelTypeNames_.read(coeffs_.lookup("thermoModel")))
 {}
 
 
@@ -152,6 +170,9 @@ void Foam::surfaceFilmModels::surfaceFilmModel::evolve()
         // Update any input information
         read();
 
+        // Pre-evolve
+        preEvolveFilm();
+
         // Increment the film equations up to the new time level
         evolveFilm();
 
@@ -163,4 +184,40 @@ void Foam::surfaceFilmModels::surfaceFilmModel::evolve()
 }
 
 
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::surfaceFilmModels::surfaceFilmModel::Srho() const
+{
+    notImplemented
+    (
+        "Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> > "
+        "Foam::surfaceFilmModels::surfaceFilmModel::Srho() const"
+    )
+    return tmp<DimensionedField<scalar, volMesh> >(NULL);
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::surfaceFilmModels::surfaceFilmModel::Srho(const label) const
+{
+    notImplemented
+    (
+        "Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> > "
+        "Foam::surfaceFilmModels::surfaceFilmModel::Srho(const label) const"
+    )
+    return tmp<DimensionedField<scalar, volMesh> >(NULL);
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::surfaceFilmModels::surfaceFilmModel::Sh() const
+{
+    notImplemented
+    (
+        "Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> > "
+        "Foam::surfaceFilmModels::surfaceFilmModel::Sh() const"
+    )
+    return tmp<DimensionedField<scalar, volMesh> >(NULL);
+}
+
+
 // ************************************************************************* //
diff --git a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.H b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.H
index 1be3ad456b8..9e9030fe2bb 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.H
+++ b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModel.H
@@ -40,7 +40,9 @@ SourceFiles
 #include "dimensionedVector.H"
 #include "runTimeSelectionTables.H"
 #include "volFieldsFwd.H"
+#include "DimensionedField.H"
 #include "labelList.H"
+#include "NamedEnum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -62,6 +64,21 @@ class surfaceFilmModel
 :
     public IOdictionary
 {
+public:
+
+    // Data types
+
+        //- Enumeration listing the possible thermo types
+        enum thermoModelType
+        {
+            tmConstant,
+            tmSingleComponent
+        };
+
+        //- Named enumeration for the thermoType
+        static const NamedEnum<thermoModelType, 2> thermoModelTypeNames_;
+
+
 private:
 
     // Private Member Functions
@@ -95,15 +112,15 @@ protected:
         //- Model coefficients dictionary
         dictionary coeffs_;
 
+        //- Thermo type
+        thermoModelType thermoModel_;
+
 
     // Protected member functions
 
         //- Read control parameters from dictionary
         virtual bool read();
 
-        //- Evolve the film
-        virtual void evolveFilm() = 0;
-
 
 public:
 
@@ -170,6 +187,9 @@ public:
             //- Return the model coefficients dictionary
             inline const dictionary& coeffs() const;
 
+            //- Return the thermo type
+            inline const thermoModelType& thermoModel() const;
+
             //- Return the film mesh database
             virtual const fvMesh& film() const = 0;
 
@@ -202,15 +222,30 @@ public:
            //- Return the film velocity [m/s]
            virtual const volVectorField& U() const = 0;
 
+           //- Return the film surface velocity [m/s]
+           virtual const volVectorField& Us() const = 0;
+
+           //- Return the film wall velocity [m/s]
+           virtual const volVectorField& Uw() const = 0;
+
            //- Return the film density [kg/m3]
            virtual const volScalarField& rho() const = 0;
 
-           //- Return the film temperature [K]
+           //- Return the film mean temperature [K]
            virtual const volScalarField& T() const = 0;
 
-           //- Return the film specific heat capacity [J/kg/K]
+           //- Return the film surface temperature [K]
+           virtual const volScalarField& Ts() const = 0;
+
+           //- Return the film wall temperature [K]
+           virtual const volScalarField& Tw() const = 0;
+
+            //- Return the film specific heat capacity [J/kg/K]
            virtual const volScalarField& cp() const = 0;
 
+           //- Return the film thermal conductivity [W/m/K]
+           virtual const volScalarField& kappa() const = 0;
+
 
            // Transfer fields - to the primary region
 
@@ -220,14 +255,41 @@ public:
                //- Return the parcel diameters originating from film
                virtual const volScalarField& diametersForPrimary() const = 0;
 
+               //- Return the film mass evolved via phase change
+               virtual const volScalarField& massPhaseChangeForPrimary()
+                   const = 0;
+
 
        // Evolution
 
             //- Evolve the film
             virtual void evolve();
 
+            //- Pre-evolve film hook
+            virtual void preEvolveFilm() = 0;
+
+            //- Evolve the film
+            virtual void evolveFilm() = 0;
+
+
+        // Source fields
+
+            // Mapped into primary region
+
+                //- Return total mass source - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Srho() const;
+
+                //- Return mass source for specie i - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Srho
+                (
+                    const label i
+                ) const;
+
+                //- Return enthalpy source - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Sh() const;
+
 
-       // I-O
+        // I-O
 
             //- Provide some feedback
             virtual void info() const = 0;
diff --git a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModelI.H b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModelI.H
index 9adf1d72cb2..0585867dc01 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModelI.H
+++ b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModelI.H
@@ -41,4 +41,11 @@ Foam::surfaceFilmModels::surfaceFilmModel::coeffs() const
 }
 
 
+inline const Foam::surfaceFilmModels::surfaceFilmModel::thermoModelType&
+Foam::surfaceFilmModels::surfaceFilmModel::thermoModel() const
+{
+    return thermoModel_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModelNew.C b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModelNew.C
index 57347f91325..fb9d7e3bb6c 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModelNew.C
+++ b/src/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel/surfaceFilmModelNew.C
@@ -36,22 +36,24 @@ Foam::surfaceFilmModels::surfaceFilmModel::New
     const dimensionedVector& g
 )
 {
-    // get model name, but do not register the dictionary
-    const word modelType
-    (
-        IOdictionary
+    word modelType;
+
+    {
+        IOdictionary surfaceFilmPropertiesDict
         (
             IOobject
             (
                 "surfaceFilmProperties",
                 mesh.time().constant(),
                 mesh,
-                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::MUST_READ,
                 IOobject::NO_WRITE,
                 false
             )
-        ).lookup("surfaceFilmModel")
-    );
+        );
+
+        surfaceFilmPropertiesDict.lookup("surfaceFilmModel") >> modelType;
+    }
 
     Info<< "Selecting surfaceFilmModel " << modelType << endl;
 
diff --git a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
index b39c8eb22f1..ac251146e89 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.C
@@ -27,11 +27,13 @@ License
 #include "fvcDiv.H"
 #include "fvcLaplacian.H"
 #include "fvm.H"
-#include "phaseChangeModel.H"
 #include "addToRunTimeSelectionTable.H"
+#include "zeroGradientFvPatchFields.H"
+#include "directMappedWallPolyPatch.H"
 
 // Sub-models
-#include "injectionModel.H"
+#include "heatTransferModel.H"
+#include "phaseChangeModel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -49,38 +51,94 @@ namespace Foam
 
 bool Foam::surfaceFilmModels::thermoSingleLayer::read()
 {
-    if (kinematicSingleLayer::read())
-    {
-        coeffs_.lookup("htcw") >> htcw_;
-        coeffs_.lookup("htcs") >> htcs_;
+    // no additional properties to read
+    return kinematicSingleLayer::read();
+}
 
-        return true;
-    }
-    else
-    {
-        return false;
-    }
+
+void Foam::surfaceFilmModels::thermoSingleLayer::resetPrimaryRegionSourceTerms()
+{
+    kinematicSingleLayer::resetPrimaryRegionSourceTerms();
+
+    hsSpPrimary_ == dimensionedScalar("zero", hsSp_.dimensions(), 0.0);
 }
 
 
-void Foam::surfaceFilmModels::thermoSingleLayer::initialise()
+void Foam::surfaceFilmModels::thermoSingleLayer::correctThermoFields()
 {
-    if (debug)
+    switch (thermoModel_)
     {
-        Pout<< "thermoSingleLayer::initialise()" << endl;
-    }
-
-    kinematicSingleLayer::initialise();
+        case tmConstant:
+        {
+            rho_ == dimensionedScalar(coeffs_.lookup("rho0"));
+            mu_ == dimensionedScalar(coeffs_.lookup("mu0"));
+            sigma_ == dimensionedScalar(coeffs_.lookup("sigma0"));
+            cp_ == dimensionedScalar(coeffs_.lookup("cp0"));
+            kappa_ == dimensionedScalar(coeffs_.lookup("kappa0"));
 
-    hs_ == hs(T_);
+            break;
+        }
+        case tmSingleComponent:
+        {
+            const liquid& liq = thermo_.liquids().properties()[liquidId_];
+            forAll(rho_, cellI)
+            {
+                const scalar T = T_[cellI];
+                const scalar p = pPrimary_[cellI];
+                rho_[cellI] = liq.rho(p, T);
+                mu_[cellI] = liq.mu(p, T);
+                sigma_[cellI] = liq.sigma(p, T);
+                cp_[cellI] = liq.cp(p, T);
+                kappa_[cellI] = liq.K(p, T);
+            }
+
+            rho_.correctBoundaryConditions();
+            mu_.correctBoundaryConditions();
+            sigma_.correctBoundaryConditions();
+            cp_.correctBoundaryConditions();
+            kappa_.correctBoundaryConditions();
+
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "void Foam::surfaceFilmModels::thermoSingleLayer::"
+                "correctThermoFields()"
+            )   << "Unknown thermoType enumeration" << abort(FatalError);
+        }
+    }
 }
 
 
-void Foam::surfaceFilmModels::thermoSingleLayer::resetPrimaryRegionSourceTerms()
+void Foam::surfaceFilmModels::thermoSingleLayer::updateSurfaceTemperatures()
 {
-    kinematicSingleLayer::resetPrimaryRegionSourceTerms();
+    // Push boundary film temperature values into internal field
+    for (label i=0; i<filmBottomPatchIDs_.size(); i++)
+    {
+        label patchI = filmBottomPatchIDs_[i];
+        const polyPatch& pp = filmRegion_.boundaryMesh()[patchI];
+        UIndirectList<scalar>(Tw_, pp.faceCells()) =
+            T_.boundaryField()[patchI];
+    }
+    Tw_.correctBoundaryConditions();
 
-    hsSpPrimary_ == dimensionedScalar("zero", hsSp_.dimensions(), 0.0);
+    // Update film surface temperature
+    dimensionedScalar deltaSmall("SMALL", dimLength, SMALL);
+    volScalarField kappaDeltaBy2 = kappa_/(0.5*delta_ + deltaSmall);
+    Ts_ =
+        (
+            // qRad
+          - massPhaseChangeForPrimary_*hs_/(time_.deltaT()*magSf_)
+          + TPrimary_*htcs_->h()
+          + kappaDeltaBy2*T_
+        )
+       /(
+            htcs_->h()
+          + kappaDeltaBy2
+        );
+    Ts_.correctBoundaryConditions();
 }
 
 
@@ -88,9 +146,13 @@ void Foam::surfaceFilmModels::thermoSingleLayer::transferPrimaryRegionFields()
 {
     kinematicSingleLayer::transferPrimaryRegionFields();
 
-    // Update temperature from primary region via direct mapped (coupled)
+    // Update primary region fileds via direct mapped (coupled)
     // boundary conditions
     TPrimary_.correctBoundaryConditions();
+    forAll(YPrimary_, i)
+    {
+        YPrimary_[i].correctBoundaryConditions();
+    }
 
     // Retrieve the source fields from the primary region via direct mapped
     // (coupled) boundary conditions
@@ -100,20 +162,31 @@ void Foam::surfaceFilmModels::thermoSingleLayer::transferPrimaryRegionFields()
 
     // Convert accummulated source terms into per unit area per unit time
     // Note: boundary values will still have original (neat) values
-    const scalar deltaT = filmRegion_.time().deltaTValue();
+    const scalar deltaT = time_.deltaTValue();
     hsSp_.field() /= magSf_*deltaT;
 }
 
 
 void Foam::surfaceFilmModels::thermoSingleLayer::updateSubmodels()
 {
-    kinematicSingleLayer::updateSubmodels();
+    // Update heat transfer coefficient sub-models
+    htcs_->correct();
+    htcw_->correct();
 
-    const dimensionedScalar deltaT = filmRegion_.time().deltaT();
-    hsSp_ -= massForPrimary_*hs_/magSf_/deltaT;
+    // Update phase change
+    massPhaseChangeForPrimary_ == dimensionedScalar("zero", dimMass, 0.0);
+    phaseChange_->correct(time_.deltaTValue(), massPhaseChangeForPrimary_);
+    massPhaseChangeForPrimary_.correctBoundaryConditions();
+    totalMassPhaseChange_ += sum(massPhaseChangeForPrimary_).value();
 
-    // Update the sub-models
-    phaseChange_->correct();
+    // Update kinematic sub-models
+    kinematicSingleLayer::updateSubmodels();
+
+    // Update source fields
+    const dimensionedScalar deltaT = time_.deltaT();
+    rhoSp_ -= massPhaseChangeForPrimary_/magSf_/deltaT;
+    USp_ -= massPhaseChangeForPrimary_*U_/magSf_/deltaT;
+    hsSp_ -= (massForPrimary_ + massPhaseChangeForPrimary_)*hs_/magSf_/deltaT;
 }
 
 
@@ -122,39 +195,12 @@ Foam::tmp<Foam::fvScalarMatrix> Foam::surfaceFilmModels::thermoSingleLayer::q
     volScalarField& hs
 ) const
 {
-    DimensionedField<scalar, volMesh> Tw
-    (
-        IOobject
-        (
-            "Tw",
-            time_.timeName(),
-            filmRegion_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        filmRegion_,
-        dimensionedScalar("zero", dimTemperature, 0.0)
-    );
-
-    forAll(filmBottomPatchIDs_, i)
-    {
-        label patchI = filmBottomPatchIDs_[i];
-        const polyPatch& pp = filmRegion_.boundaryMesh()[patchI];
-        UIndirectList<scalar>(Tw, pp.faceCells()) =
-            TPrimary_.boundaryField()[patchI];
-    }
-
-    // TODO: Use T at film thickness instead of cell value
-    DimensionedField<scalar, volMesh> Ts =
-        TPrimary_.dimensionedInternalField();
-
     return
     (
-      - fvm::Sp(htcs_/cp_, hs)
-      + htcs_*(dimensionedScalar("Tstd", dimTemperature, 298.15) - Ts)
-      - fvm::Sp(htcw_/cp_, hs)
-      + htcw_*(dimensionedScalar("Tstd", dimTemperature, 298.15) - Tw)
+      - fvm::Sp(htcs_->h()/cp_, hs)
+      - htcs_->h()*(dimensionedScalar("Tstd", dimTemperature, 298.15) - Ts_)
+      - fvm::Sp(htcw_->h()/cp_, hs)
+      - htcw_->h()*(dimensionedScalar("Tstd", dimTemperature, 298.15) - Tw_)
     );
 }
 
@@ -166,6 +212,8 @@ void Foam::surfaceFilmModels::thermoSingleLayer::solveEnergy()
         Info<< "thermoSingleLayer::solveEnergy()" << endl;
     }
 
+    updateSurfaceTemperatures();
+
     solve
     (
         fvm::ddt(deltaRho_, hs_)
@@ -174,6 +222,8 @@ void Foam::surfaceFilmModels::thermoSingleLayer::solveEnergy()
         hsSp_
       + q(hs_)
     );
+
+    correctThermoFields();
 }
 
 
@@ -187,6 +237,8 @@ Foam::surfaceFilmModels::thermoSingleLayer::thermoSingleLayer
 )
 :
     kinematicSingleLayer(modelType, mesh, g),
+    thermo_(mesh.lookupObject<SLGThermo>("SLGThermo")),
+    liquidId_(thermo_.liquidId(coeffs_.lookup("liquid"))),
     cp_
     (
         IOobject
@@ -195,19 +247,30 @@ Foam::surfaceFilmModels::thermoSingleLayer::thermoSingleLayer
             time_.timeName(),
             filmRegion_,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         filmRegion_,
-        dimensionedScalar(coeffs_.lookup("cp"))
+        dimensionedScalar("cp", dimEnergy/dimMass/dimTemperature, 0.0),
+        zeroGradientFvPatchScalarField::typeName
     ),
-
-    htcw_
-    (
-        dimensionedScalar("htc", dimEnergy/dimTime/dimArea/dimTemperature, 0.0)
-    ),
-    htcs_
+    kappa_
     (
-        dimensionedScalar("htc", dimEnergy/dimTime/dimArea/dimTemperature, 0.0)
+        IOobject
+        (
+            "kappa",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        filmRegion_,
+        dimensionedScalar
+        (
+            "kappa",
+            dimEnergy/dimTime/dimLength/dimTemperature,
+            0.0
+        ),
+        zeroGradientFvPatchScalarField::typeName
     ),
 
     T_
@@ -222,6 +285,32 @@ Foam::surfaceFilmModels::thermoSingleLayer::thermoSingleLayer
         ),
         filmRegion_
     ),
+    Ts_
+    (
+        IOobject
+        (
+            "Tsf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        T_,
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    Tw_
+    (
+        IOobject
+        (
+            "Twf",
+            time_.timeName(),
+            filmRegion_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        T_,
+        zeroGradientFvPatchScalarField::typeName
+    ),
     hs_
     (
         IOobject
@@ -279,13 +368,46 @@ Foam::surfaceFilmModels::thermoSingleLayer::thermoSingleLayer
         filmRegion_
     ),
 
-    phaseChange_(phaseChangeModel::New(*this, coeffs_)),
+    YPrimary_(),
 
-    hsSpDetach_(filmRegion_.nCells(), 0.0)
+    htcs_
+    (
+        heatTransferModel::New(*this, coeffs_.subDict("upperSurfaceModels"))
+    ),
+    htcw_
+    (
+        heatTransferModel::New(*this, coeffs_.subDict("lowerSurfaceModels"))
+    ),
+    phaseChange_(phaseChangeModel::New(*this, coeffs_)),
+    totalMassPhaseChange_(0.0)
 {
-    read();
+    if (thermo_.hasMultiComponentCarrier())
+    {
+        YPrimary_.setSize(thermo_.carrier().species().size());
+
+        forAll(thermo_.carrier().species(), i)
+        {
+            YPrimary_.set
+            (
+                i,
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        thermo_.carrier().species()[i],
+                        time_.timeName(),
+                        filmRegion_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    filmRegion_,
+                    dimensionedScalar("zero", dimless, 0.0),
+                    pSp_.boundaryField().types()
+                )
+            );
+        }
+    }
 
-    initialise();
 }
 
 
@@ -326,6 +448,25 @@ void Foam::surfaceFilmModels::thermoSingleLayer::addSources
 }
 
 
+void Foam::surfaceFilmModels::thermoSingleLayer::preEvolveFilm()
+{
+    if (!initialisedThermo_)
+    {
+        // Retreive pressure from primary region
+        pPrimary_.correctBoundaryConditions();
+
+        // Correct (temperature dependent) thermo fields
+        correctThermoFields();
+
+        // Update derived fields
+        hs_ == hs(T_);
+        deltaRho_ == delta_*rho_;
+        phi_ = fvc::interpolate(deltaRho_*U_) & filmRegion_.Sf();
+        initialisedThermo_ = true;
+    }
+}
+
+
 void Foam::surfaceFilmModels::thermoSingleLayer::evolveFilm()
 {
     transferPrimaryRegionFields();
@@ -338,22 +479,22 @@ void Foam::surfaceFilmModels::thermoSingleLayer::evolveFilm()
     for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
     {
         // Explicit pressure source contribution
-        tmp<volScalarField> pu = this->pu();
+        tmp<volScalarField> tpu = this->pu();
 
         // Implicit pressure source coefficient
-        tmp<volScalarField> pp = this->pp();
+        tmp<volScalarField> tpp = this->pp();
 
         // Solve for momentum for U_
-        tmp<fvVectorMatrix> UEqn = solveMomentum(pu(), pp());
+        tmp<fvVectorMatrix> UEqn = solveMomentum(tpu(), tpp());
+
+        // Solve energy for hs_
+        solveEnergy();
 
         // Film thickness correction loop
         for (int corr=1; corr<=nCorr_; corr++)
         {
-            // Solve energy for hs_
-            solveEnergy();
-
             // Solve thickness for delta_
-            solveThickness(pu(), pp(), UEqn());
+            solveThickness(tpu(), tpp(), UEqn());
         }
     }
 
@@ -363,11 +504,31 @@ void Foam::surfaceFilmModels::thermoSingleLayer::evolveFilm()
     // Update temperature using latest hs_
     T_ == T(hs_);
 
+    // Update film wall and surface velocities
+    updateSurfaceVelocities();
+
+    // Update film wall and surface temperatures
+    updateSurfaceTemperatures();
+
     // Reset source terms for next time integration
     resetPrimaryRegionSourceTerms();
 }
 
 
+const Foam::volScalarField&
+Foam::surfaceFilmModels::thermoSingleLayer::cp() const
+{
+    return cp_;
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::thermoSingleLayer::kappa() const
+{
+    return kappa_;
+}
+
+
 const Foam::volScalarField&
 Foam::surfaceFilmModels::thermoSingleLayer::T() const
 {
@@ -376,9 +537,23 @@ Foam::surfaceFilmModels::thermoSingleLayer::T() const
 
 
 const Foam::volScalarField&
-Foam::surfaceFilmModels::thermoSingleLayer::cp() const
+Foam::surfaceFilmModels::thermoSingleLayer::Ts() const
 {
-    return cp_;
+    return Ts_;
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::thermoSingleLayer::Tw() const
+{
+    return Tw_;
+}
+
+
+const Foam::volScalarField&
+Foam::surfaceFilmModels::thermoSingleLayer::hs() const
+{
+    return hs_;
 }
 
 
@@ -386,8 +561,128 @@ void Foam::surfaceFilmModels::thermoSingleLayer::info() const
 {
     kinematicSingleLayer::info();
 
-    Info<< indent<< "min/max(T)      = " << min(T_).value() << ", "
-        << max(T_).value() << nl;
+    Info<< indent << "min/max(T)        = " << min(T_).value() << ", "
+        << max(T_).value() << nl
+        << indent << "mass phase change = "
+        << returnReduce(totalMassPhaseChange_, sumOp<scalar>()) << nl;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::surfaceFilmModels::thermoSingleLayer::Srho(const label i) const
+{
+    const label vapId =
+        thermo_.carrierId(thermo_.liquids().components()[liquidId_]);
+
+    tmp<DimensionedField<scalar, volMesh> > tSrho
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "thermoSingleLayer::Srho(i)",
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
+
+    if (vapId == i)
+    {
+        scalarField& Srho = tSrho();
+        const scalarField& V = mesh_.V();
+        const scalar dt = time_.deltaTValue();
+
+        forAll(filmBottomPatchIDs_, i)
+        {
+            const label primaryPatchI = primaryPatchIDs_[i];
+            const directMappedWallPolyPatch& wpp =
+                refCast<const directMappedWallPolyPatch>
+                (
+                     mesh_.boundaryMesh()[primaryPatchI]
+                );
+
+            const mapDistribute& distMap = wpp.map();
+
+            const label filmPatchI = filmBottomPatchIDs_[i];
+
+            scalarField patchMass =
+                massPhaseChangeForPrimary().boundaryField()[filmPatchI];
+
+            distMap.distribute(patchMass);
+
+            const unallocLabelList& cells = wpp.faceCells();
+
+            forAll(patchMass, j)
+            {
+                Srho[cells[j]] = patchMass[j]/(V[cells[j]]*dt);
+            }
+        }
+    }
+
+    return tSrho;
+}
+
+
+Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::surfaceFilmModels::thermoSingleLayer::Sh() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tSh
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                "thermoSingleLayer::Sh",
+                time_.timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+        )
+    );
+
+    scalarField& Sh = tSh();
+    const scalarField& V = mesh_.V();
+    const scalar dt = time_.deltaTValue();
+
+    forAll(filmBottomPatchIDs_, i)
+    {
+        const label primaryPatchI = primaryPatchIDs_[i];
+        const directMappedWallPolyPatch& wpp =
+            refCast<const directMappedWallPolyPatch>
+            (
+                 mesh_.boundaryMesh()[primaryPatchI]
+            );
+
+        const mapDistribute& distMap = wpp.map();
+
+        const label filmPatchI = filmBottomPatchIDs_[i];
+
+        scalarField patchMass =
+            massPhaseChangeForPrimary().boundaryField()[filmPatchI];
+        distMap.distribute(patchMass);
+
+        scalarField patchEnthalpy = hs_.boundaryField()[filmPatchI];
+        distMap.distribute(patchEnthalpy);
+
+        const unallocLabelList& cells = wpp.faceCells();
+
+        forAll(patchMass, j)
+        {
+            Sh[cells[j]] += patchMass[j]*patchEnthalpy[j]/(V[cells[j]]*dt);
+        }
+    }
+
+    return tSh;
 }
 
 
diff --git a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.H b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.H
index eb7fb7800a1..a7255a99944 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.H
+++ b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayer.H
@@ -40,6 +40,7 @@ SourceFiles
 #define thermoSingleLayer_H
 
 #include "kinematicSingleLayer.H"
+#include "SLGThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,6 +50,7 @@ namespace surfaceFilmModels
 {
 
 // Forward declaration of classes
+class heatTransferModel;
 class phaseChangeModel;
 
 /*---------------------------------------------------------------------------*\
@@ -74,29 +76,36 @@ protected:
 
     // Protected data
 
-        // Thermo properties  TODO: currently read from coeffs dictionary
+        // Thermo properties
 
-            //- Specific heat capacity / [J/kg/K]
-            volScalarField cp_;
+            //- Reference to the SLGThermo
+            const SLGThermo& thermo_;
 
+            // Single component
 
-        // Model parameters
+                //- Id of component in thermo database
+                label liquidId_;
 
-            //- Heat transfer coefficient bewteen wall and film [W/m2/K]
-            dimensionedScalar htcw_;
 
-            //- Heat transfer coefficient bewteen film surface and primary
-            //  region [W/m2/K]
-            dimensionedScalar htcs_;
+            // Fields
 
+                //- Specific heat capacity / [J/kg/K]
+                volScalarField cp_;
 
-        // Fields
+                //- Thermal conductivity / [W/m/K]
+                volScalarField kappa_;
 
-            //- Temperature / [K]
-            volScalarField T_;
+                //- Temperature - mean / [K]
+                volScalarField T_;
+
+                //- Temperature - surface / [K]
+                volScalarField Ts_;
 
-            //- Sensible)enthalpy / [J/kg]
-            volScalarField hs_;
+                //- Temperature - wall / [K]
+                volScalarField Tw_;
+
+                //- Sensible enthalpy / [J/kg]
+                volScalarField hs_;
 
 
         // Source term fields
@@ -122,27 +131,37 @@ protected:
             //- Temperature / [K]
             volScalarField TPrimary_;
 
+            //- List of specie mass fractions / [0-1]
+            PtrList<volScalarField> YPrimary_;
+
 
         // Sub-models
 
-            //- Phase change
-            autoPtr<phaseChangeModel> phaseChange_;
+            //- Heat transfer coefficient bewteen film surface and primary
+            //  region [W/m2/K]
+            autoPtr<heatTransferModel> htcs_;
 
+            //- Heat transfer coefficient bewteen wall and film [W/m2/K]
+            autoPtr<heatTransferModel> htcw_;
 
-       // Detached surface properties
+            //- Phase change
+            autoPtr<phaseChangeModel> phaseChange_;
 
-           //- Energy sink from film to carrier phase [kg]
-           scalarField hsSpDetach_;
+            //- Total mass transferred to primary region [kg]
+            scalar totalMassPhaseChange_;
 
 
     // Protected member functions
 
-        //- Initialise the film model - called on construction
-        void initialise();
-
         //- Read control parameters from dictionary
         virtual bool read();
 
+        //- Correct the thermo fields
+        virtual void correctThermoFields();
+
+        //- Correct the film surface and wall temperatures
+        virtual void updateSurfaceTemperatures();
+
         //- Reset source term fields
         virtual void resetPrimaryRegionSourceTerms();
 
@@ -185,38 +204,51 @@ public:
 
     // Member Functions
 
-        // Model parameters
+        // Thermo properties
 
-            //- Return the heat transfer coefficient bewteen wall and film
-            inline const dimensionedScalar& htcw() const;
+            //- Return const reference to the SLGThermo object
+            inline const SLGThermo& thermo() const;
 
-            //- Return the Heat transfer coefficient bewteen film surface and
-            //  primary region
-            inline const dimensionedScalar& htcs() const;
+            // Single component
 
+                //- Return the Id of component in thermo database
+                inline label liquidId() const;
 
-        // Thermo properties
 
-            //- Return sensible enthalpy as a function of temperature
-            inline tmp<Foam::volScalarField> hs
-            (
-                const volScalarField& T
-            ) const;
+            // Fields
 
-            //- Return temperature as a function of sensible enthalpy
-            inline tmp<Foam::volScalarField> T
-            (
-                const volScalarField& hs
-            ) const;
+                //- Return the film specific heat capacity [J/kg/K]
+                virtual const volScalarField& cp() const;
+
+                //- Return the film thermal conductivity [W/m/K]
+                virtual const volScalarField& kappa() const;
+
+                //- Return the film mean temperature [K]
+                virtual const volScalarField& T() const;
 
+                //- Return the film surface temperature [K]
+                virtual const volScalarField& Ts() const;
 
-        // Fields
+                //- Return the film wall temperature [K]
+                virtual const volScalarField& Tw() const;
 
-            //- Return the film temperature [K]
-            virtual const volScalarField& T() const;
+                //- Return the film sensible enthalpy [J/kg]
+                virtual const volScalarField& hs() const;
 
-            //- Return the film specific heat capacity [J/kg/K]
-            virtual const volScalarField& cp() const;
+
+            // Helper functions
+
+                //- Return sensible enthalpy as a function of temperature
+                inline tmp<Foam::volScalarField> hs
+                (
+                    const volScalarField& T
+                ) const;
+
+                //- Return temperature as a function of sensible enthalpy
+                inline tmp<Foam::volScalarField> T
+                (
+                    const volScalarField& hs
+                ) const;
 
 
          // Source fields (read/write access)
@@ -233,30 +265,65 @@ public:
             );
 
 
-            // Primary region
+        // Source term fields
+
+            // Film region
 
                 //- Energy / [J/m2/s]
-                inline volScalarField& hsSpPrimary();
+                inline const volScalarField& hsSp() const;
 
 
-            // Film region
+            // Primary region
 
                 //- Energy / [J/m2/s]
-                inline volScalarField& hsSp();
+                inline const volScalarField& hsSpPrimary() const;
+
+
+        // Fields mapped from the primary region
+
+            //- Temperature / [K]
+            inline const volScalarField& TPrimary() const;
+
+            //- Specie mass fractions / [0-1]
+            inline const PtrList<volScalarField>& YPrimary() const;
+
 
 
         // Sub-models
 
-            //- Phase change
-            inline phaseChangeModel& phaseChange();
+            //- Return const access to the (surface) heat transfer model
+            inline const heatTransferModel& htcs() const;
+
+            //- Return const access to the (wall) heat transfer model
+            inline const heatTransferModel& htcw() const;
+
+            //- Return const access to the phawse change model
+            inline const phaseChangeModel& phaseChange() const;
 
 
         // Evolution
 
+            //- Pre-evolve film hook
+            virtual void preEvolveFilm();
+
             //- Evolve the film equations
             virtual void evolveFilm();
 
 
+        // Source fields
+
+            // Mapped into primary region
+
+                //- Return mass source for specie i - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Srho
+                (
+                    const label i
+                ) const;
+
+                //- Return enthalpy source - Eulerian phase only
+                virtual tmp<DimensionedField<scalar, volMesh> > Sh() const;
+
+
        // I-O
 
             //- Provide some feedback
diff --git a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayerI.H b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayerI.H
index 987816cd30a..a6c979eca8a 100644
--- a/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayerI.H
+++ b/src/surfaceFilmModels/surfaceFilmModel/thermoSingleLayer/thermoSingleLayerI.H
@@ -27,17 +27,16 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-inline const Foam::dimensionedScalar&
-Foam::surfaceFilmModels::thermoSingleLayer::htcw() const
+inline const Foam::SLGThermo&
+Foam::surfaceFilmModels::thermoSingleLayer::thermo() const
 {
-    return htcw_;
+    return thermo_;
 }
 
 
-inline const Foam::dimensionedScalar&
-Foam::surfaceFilmModels::thermoSingleLayer::htcs() const
+inline Foam::label Foam::surfaceFilmModels::thermoSingleLayer::liquidId() const
 {
-    return htcs_;
+    return liquidId_;
 }
 
 
@@ -59,7 +58,8 @@ Foam::surfaceFilmModels::thermoSingleLayer::hs
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            cp_*(T - (dimensionedScalar("Tstd", dimTemperature, 298.15)))
+            cp_*(T - (dimensionedScalar("Tstd", dimTemperature, 298.15))),
+            zeroGradientFvPatchScalarField::typeName
         )
     );
 }
@@ -83,28 +83,57 @@ Foam::surfaceFilmModels::thermoSingleLayer::T
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            hs/cp_ + dimensionedScalar("Tstd", dimTemperature, 298.15)
+            hs/cp_ + dimensionedScalar("Tstd", dimTemperature, 298.15),
+            zeroGradientFvPatchScalarField::typeName
         )
     );
 }
 
 
-inline Foam::volScalarField&
-Foam::surfaceFilmModels::thermoSingleLayer::hsSpPrimary()
+inline const Foam::volScalarField&
+Foam::surfaceFilmModels::thermoSingleLayer::hsSp() const
+{
+    return hsSp_;
+}
+
+
+inline const Foam::volScalarField&
+Foam::surfaceFilmModels::thermoSingleLayer::hsSpPrimary() const
 {
     return hsSpPrimary_;
 }
 
 
-inline Foam::volScalarField&
-Foam::surfaceFilmModels::thermoSingleLayer::hsSp()
+inline const Foam::volScalarField&
+Foam::surfaceFilmModels::thermoSingleLayer::TPrimary() const
 {
-    return hsSp_;
+    return TPrimary_;
+}
+
+
+inline const Foam::PtrList<Foam::volScalarField>&
+Foam::surfaceFilmModels::thermoSingleLayer::YPrimary() const
+{
+    return YPrimary_;
+}
+
+
+inline const Foam::surfaceFilmModels::heatTransferModel&
+Foam::surfaceFilmModels::thermoSingleLayer::htcs() const
+{
+    return htcs_();
+}
+
+
+inline const Foam::surfaceFilmModels::heatTransferModel&
+Foam::surfaceFilmModels::thermoSingleLayer::htcw() const
+{
+    return htcw_();
 }
 
 
-inline Foam::surfaceFilmModels::phaseChangeModel&
-Foam::surfaceFilmModels::thermoSingleLayer::phaseChange()
+inline const Foam::surfaceFilmModels::phaseChangeModel&
+Foam::surfaceFilmModels::thermoSingleLayer::phaseChange() const
 {
     return phaseChange_();
 }
-- 
GitLab