diff --git a/src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.C b/src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.C
index 34b9da0da1..e282ff8c21 100644
--- a/src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.C
+++ b/src/OpenFOAM/db/functionObjects/outputFilterOutputControl/outputFilterOutputControl.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,11 +30,8 @@ License
 namespace Foam
 {
     template<>
-    const char* Foam::NamedEnum
-    <
-        Foam::outputFilterOutputControl::outputControls,
-        2
-    >::names[] =
+    const char* NamedEnum<outputFilterOutputControl::outputControls, 2>::
+    names[] =
     {
         "timeStep",
         "outputTime"
diff --git a/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.C b/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.C
index bef68f3c4e..34db0a7e3e 100644
--- a/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.C
+++ b/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.C
@@ -41,8 +41,18 @@ namespace Foam
         fixedTemperatureSource,
         dictionary
     );
+
+    template<>
+    const char* NamedEnum<fixedTemperatureSource::temperatureMode, 2>::names[] =
+    {
+        "constant",
+        "lookup"
+    };
 }
 
+const Foam::NamedEnum<Foam::fixedTemperatureSource::temperatureMode, 2>
+    Foam::fixedTemperatureSource::temperatureModeNames_;
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -55,8 +65,29 @@ Foam::fixedTemperatureSource::fixedTemperatureSource
 )
 :
     basicSource(name, modelType, dict, mesh),
-    T_(readScalar(coeffs_.lookup("temperature")))
+    mode_(temperatureModeNames_.read(coeffs_.lookup("mode"))),
+    Tconstant_(0.0),
+    TName_("T")
 {
+    switch (mode_)
+    {
+        case tmConstant:
+        {
+            coeffs_.lookup("temperature") >> Tconstant_;
+            break;
+        }
+        case tmLookup:
+        {
+            TName_ = coeffs_.lookupOrDefault<word>("TName", "T");
+            break;
+        }
+        default:
+        {
+            // error handling done by NamedEnum
+        }
+    }
+
+
     fieldNames_.setSize(1, "energy");
     applied_.setSize(1, false);
 }
@@ -81,8 +112,31 @@ void Foam::fixedTemperatureSource::setValue
 
     if (eqn.psi().name() == thermo.he().name())
     {
-        const scalarField Tfield(cells_.size(), T_);
-        eqn.setValues(cells_, thermo.he(thermo.p(), Tfield, cells_));
+        switch (mode_)
+        {
+            case tmConstant:
+            {
+                scalarField Tconst(cells_.size(), Tconstant_);
+                eqn.setValues(cells_, thermo.he(thermo.p(), Tconst, cells_));
+
+                break;
+            }
+            case tmLookup:
+            {
+                const volScalarField& T =
+                    mesh().lookupObject<volScalarField>(TName_);
+
+                scalarField Tlookup(T, cells_);
+                eqn.setValues(cells_, thermo.he(thermo.p(), Tlookup, cells_));
+
+                break;
+            }
+            default:
+            {
+                // error handling done by NamedEnum
+            }
+        }
+
     }
 }
 
@@ -98,7 +152,8 @@ bool Foam::fixedTemperatureSource::read(const dictionary& dict)
 {
     if (basicSource::read(dict))
     {
-        coeffs_.readIfPresent("T", T_);
+        coeffs_.readIfPresent("temperature", Tconstant_);
+        coeffs_.readIfPresent("TName", TName_);
 
         return true;
     }
diff --git a/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.H b/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.H
index 4281a38e2c..3aefdb5645 100644
--- a/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.H
+++ b/src/fieldSources/derived/fixedTemperatureSource/fixedTemperatureSource.H
@@ -33,7 +33,13 @@ Description
         fixedTemperatureSourceCoeffs
         {
             fieldNames      (h e hs);   // names of fields to apply source
+            mode            constant;   // constant or lookup
+
+            // constant option
             temperature     500;        // fixed temperature [K]
+
+            // loolup option
+            // TName        T;          // optional temperature field name
         }
 
 
@@ -46,6 +52,7 @@ SourceFiles
 #define fixedTemperatureSource_H
 
 #include "basicSource.H"
+#include "NamedEnum.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,12 +68,32 @@ class fixedTemperatureSource
     public basicSource
 {
 
+public:
+
+    //- Temperature mode
+    enum temperatureMode
+    {
+        tmConstant,
+        tmLookup
+    };
+
+
+    //- String representation of temperatureMode enums
+    static const NamedEnum<temperatureMode, 2> temperatureModeNames_;
+
+
 protected:
 
     // Protected data
 
-        //- Fixed temperature [K]
-        scalar T_;
+        //- Operation mode
+        temperatureMode mode_;
+
+        //- Constant temperature [K]
+        scalar Tconstant_;
+
+        //- Temperature field name
+        word TName_;
 
 
 private:
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C
index 9ef4be21a4..f3a7615b4a 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedJump/fixedJumpFvPatchField.C
@@ -63,7 +63,19 @@ Foam::fixedJumpFvPatchField<Type>::fixedJumpFvPatchField
 :
     jumpCyclicFvPatchField<Type>(p, iF),
     jump_("jump", dict, p.size())
-{}
+{
+    if (dict.found("value"))
+    {
+        fvPatchField<Type>::operator=
+        (
+            Field<Type>("value", dict, p.size())
+        );
+    }
+    else
+    {
+        this->evaluate(Pstream::blocking);
+    }
+}
 
 
 template<class Type>
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C
index 4c3290cf2d..ad6fb925ae 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fixedJumpAMI/fixedJumpAMIFvPatchField.C
@@ -63,7 +63,19 @@ Foam::fixedJumpAMIFvPatchField<Type>::fixedJumpAMIFvPatchField
 :
     jumpCyclicAMIFvPatchField<Type>(p, iF),
     jump_("jump", dict, p.size())
-{}
+{
+    if (dict.found("value"))
+    {
+        fvPatchField<Type>::operator=
+        (
+            Field<Type>("value", dict, p.size())
+        );
+    }
+    else
+    {
+        this->evaluate(Pstream::blocking);
+    }
+}
 
 
 template<class Type>
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
index fbfc5eb96a..29590cc051 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
@@ -53,7 +53,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(-1),
     endSampledValues_(0),
     endAverage_(pTraits<Type>::zero),
-    offSet_()
+    offset_()
 {}
 
 
@@ -79,7 +79,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(-1),
     endSampledValues_(0),
     endAverage_(pTraits<Type>::zero),
-    offSet_()
+    offset_()
 {}
 
 
@@ -104,7 +104,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(-1),
     endSampledValues_(0),
     endAverage_(pTraits<Type>::zero),
-    offSet_(DataEntry<Type>::New("offSet", dict))
+    offset_(DataEntry<Type>::New("offset", dict))
 {
     dict.readIfPresent("fieldTableName", fieldTableName_);
 
@@ -138,7 +138,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(ptf.endSampleTime_),
     endSampledValues_(ptf.endSampledValues_),
     endAverage_(ptf.endAverage_),
-    offSet_(ptf.offSet_().clone().ptr())
+    offset_(ptf.offset_().clone().ptr())
 {}
 
 
@@ -163,7 +163,7 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(ptf.endSampleTime_),
     endSampledValues_(ptf.endSampledValues_),
     endAverage_(ptf.endAverage_),
-    offSet_(ptf.offSet_().clone().ptr())
+    offset_(ptf.offset_().clone().ptr())
 {}
 
 
@@ -487,7 +487,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::updateCoeffs()
 
     // apply offset to mapped values
     const scalar t = this->db().time().timeOutputValue();
-    this->operator==(*this + offSet_->value(t));
+    this->operator==(*this + offset_->value(t));
 
     if (debug)
     {
@@ -512,7 +512,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::write(Ostream& os) const
             << token::END_STATEMENT << nl;
     }
 
-    offSet_->writeData(os);
+    offset_->writeData(os);
 
     this->writeEntry("value", os);
 }
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
index 01e61301eb..aec6c322b9 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
@@ -132,7 +132,7 @@ class timeVaryingMappedFixedValueFvPatchField
         Type endAverage_;
 
         //- Time varying offset values to interpolated data
-        autoPtr<DataEntry<Type> > offSet_;
+        autoPtr<DataEntry<Type> > offset_;
 
 
 public:
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C
index 06f4e3c134..0670119ab0 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/uniformJumpAMI/uniformJumpAMIFvPatchField.C
@@ -73,6 +73,10 @@ Foam::uniformJumpAMIFvPatchField<Type>::uniformJumpAMIFvPatchField
     {
         fvPatchField<Type>::operator=(Field<Type>("value", dict, p.size()));
     }
+    else
+    {
+        this->evaluate(Pstream::blocking);
+    }
 }
 
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 58217adbf5..7ec52b1ee8 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -326,7 +326,7 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         p.age() += dt;
 
-        td.cloud().functions().postMove(p, cellI, dt);
+        td.cloud().functions().postMove(p, cellI, dt, td.keepParticle);
     }
 
     return td.keepParticle;
@@ -340,7 +340,7 @@ void Foam::KinematicParcel<ParcelType>::hitFace(TrackData& td)
     typename TrackData::cloudType::parcelType& p =
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
-    td.cloud().functions().postFace(p, p.face());
+    td.cloud().functions().postFace(p, p.face(), td.keepParticle);
 }
 
 
@@ -364,7 +364,14 @@ bool Foam::KinematicParcel<ParcelType>::hitPatch
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
     // Invoke post-processing model
-    td.cloud().functions().postPatch(p, pp, trackFraction, tetIs);
+    td.cloud().functions().postPatch
+    (
+        p,
+        pp,
+        trackFraction,
+        tetIs,
+        td.keepParticle
+    );
 
     // Invoke surface film model
     if (td.cloud().surfaceFilm().transferParcel(p, pp, td.keepParticle))
diff --git a/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H b/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
index 9e19beae34..783878f480 100644
--- a/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
+++ b/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ License
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "FacePostProcessing.H"
+#include "ParticleCollector.H"
 #include "ParticleErosion.H"
 #include "ParticleTracks.H"
 #include "ParticleTrap.H"
@@ -42,6 +43,7 @@ License
     makeCloudFunctionObject(CloudType);                                       \
                                                                               \
     makeCloudFunctionObjectType(FacePostProcessing, CloudType);               \
+    makeCloudFunctionObjectType(ParticleCollector, CloudType);                \
     makeCloudFunctionObjectType(ParticleErosion, CloudType);                  \
     makeCloudFunctionObjectType(ParticleTracks, CloudType);                   \
     makeCloudFunctionObjectType(ParticleTrap, CloudType);                     \
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C
index dbc359ae9a..5a208a4ca8 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.C
@@ -96,7 +96,8 @@ void Foam::CloudFunctionObject<CloudType>::postMove
 (
     const typename CloudType::parcelType&,
     const label,
-    const scalar
+    const scalar,
+    bool&
 )
 {
     // do nothing
@@ -109,7 +110,8 @@ void Foam::CloudFunctionObject<CloudType>::postPatch
     const typename CloudType::parcelType&,
     const polyPatch&,
     const scalar,
-    const tetIndices&
+    const tetIndices&,
+    bool&
 )
 {
     // do nothing
@@ -120,7 +122,8 @@ template<class CloudType>
 void Foam::CloudFunctionObject<CloudType>::postFace
 (
     const typename CloudType::parcelType&,
-    const label
+    const label,
+    bool&
 )
 {
     // do nothing
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H
index 5baf0f7ae9..976fd3867f 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObject/CloudFunctionObject.H
@@ -46,6 +46,9 @@ SourceFiles
 namespace Foam
 {
 
+class polyPatch;
+class tetIndices;
+
 /*---------------------------------------------------------------------------*\
                     Class CloudFunctionObject Declaration
 \*---------------------------------------------------------------------------*/
@@ -134,7 +137,8 @@ public:
             (
                 const typename CloudType::parcelType& p,
                 const label cellI,
-                const scalar dt
+                const scalar dt,
+                bool& keepParticle
             );
 
             //- Post-patch hook
@@ -143,14 +147,16 @@ public:
                 const typename CloudType::parcelType& p,
                 const polyPatch& pp,
                 const scalar trackFraction,
-                const tetIndices& testIs
+                const tetIndices& testIs,
+                bool& keepParticle
             );
 
             //- Post-face hook
             virtual void postFace
             (
                 const typename CloudType::parcelType& p,
-                const label faceI
+                const label faceI,
+                bool& keepParticle
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C
index 2f97468e36..e0d57b6057 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.C
@@ -132,12 +132,18 @@ void Foam::CloudFunctionObjectList<CloudType>::postMove
 (
     const typename CloudType::parcelType& p,
     const label cellI,
-    const scalar dt
+    const scalar dt,
+    bool& keepParticle
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postMove(p, cellI, dt);
+        this->operator[](i).postMove(p, cellI, dt, keepParticle);
+
+        if (!keepParticle)
+        {
+            return;
+        }
     }
 }
 
@@ -148,12 +154,25 @@ void Foam::CloudFunctionObjectList<CloudType>::postPatch
     const typename CloudType::parcelType& p,
     const polyPatch& pp,
     const scalar trackFraction,
-    const tetIndices& tetIs
+    const tetIndices& tetIs,
+    bool& keepParticle
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postPatch(p, pp, trackFraction, tetIs);
+        this->operator[](i).postPatch
+        (
+            p,
+            pp,
+            trackFraction,
+            tetIs,
+            keepParticle
+        );
+
+        if (!keepParticle)
+        {
+            return;
+        }
     }
 }
 
@@ -162,12 +181,18 @@ template<class CloudType>
 void Foam::CloudFunctionObjectList<CloudType>::postFace
 (
     const typename CloudType::parcelType& p,
-    const label faceI
+    const label faceI,
+    bool& keepParticle
 )
 {
     forAll(*this, i)
     {
-        this->operator[](i).postFace(p, faceI);
+        this->operator[](i).postFace(p, faceI, keepParticle);
+
+        if (!keepParticle)
+        {
+            return;
+        }
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H
index 5daeab9f29..de9c3afc02 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/CloudFunctionObjectList/CloudFunctionObjectList.H
@@ -114,7 +114,8 @@ public:
             (
                 const typename CloudType::parcelType& p,
                 const label cellI,
-                const scalar dt
+                const scalar dt,
+                bool& keepParticle
             );
 
             //- Post-patch hook
@@ -123,14 +124,16 @@ public:
                 const typename CloudType::parcelType& p,
                 const polyPatch& pp,
                 const scalar trackFraction,
-                const tetIndices& tetIs
+                const tetIndices& tetIs,
+                bool& keepParticle
             );
 
             //- Post-face hook
             virtual void postFace
             (
                 const typename CloudType::parcelType& p,
-                const label faceI
+                const label faceI,
+                bool& keepParticle
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C
index ae02927bbd..8c64e1b62f 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -133,7 +133,7 @@ void Foam::FacePostProcessing<CloudType>::write()
         {
             OFstream& os = outputFilePtr_[zoneI];
             os  << time.timeName() << token::TAB << sumMassTotal << token::TAB
-                <<  sumMassFlowRate<< endl;
+                << sumMassFlowRate<< endl;
         }
     }
 
@@ -376,7 +376,8 @@ template<class CloudType>
 void Foam::FacePostProcessing<CloudType>::postFace
 (
     const parcelType& p,
-    const label faceI
+    const label faceI,
+    bool&
 )
 {
     if
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H
index 06387d234b..3195b9fc57 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H
@@ -160,7 +160,8 @@ public:
             virtual void postFace
             (
                 const parcelType& p,
-                const label faceI
+                const label faceI,
+                bool& keepParticle
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
new file mode 100644
index 0000000000..06f836a3d5
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -0,0 +1,677 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "ParticleCollector.H"
+#include "Pstream.H"
+#include "surfaceWriter.H"
+#include "unitConversion.H"
+#include "Random.H"
+#include "triangle.H"
+#include "cloud.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::makeLogFile
+(
+    const faceList& faces,
+    const Field<point>& points,
+    const Field<scalar>& area
+)
+{
+    // Create the output file if not already created
+    if (log_)
+    {
+        if (debug)
+        {
+            Info<< "Creating output file" << endl;
+        }
+
+        if (Pstream::master())
+        {
+            const fileName logDir = outputDir_/this->owner().time().timeName();
+
+            // Create directory if does not exist
+            mkDir(logDir);
+
+            // Open new file at start up
+            outputFilePtr_.reset
+            (
+                new OFstream(logDir/(type() + ".dat"))
+            );
+
+            outputFilePtr_()
+                << "# Source     : " << type() << nl
+                << "# Total area : " << sum(area) << nl
+                << "# Time";
+
+            forAll(faces, i)
+            {
+                word id = Foam::name(i);
+
+                outputFilePtr_()
+                    << tab << "area[" << id << "]"
+                    << tab << "mass[" << id << "]"
+                    << tab << "massFlowRate[" << id << "]"
+                    << endl;
+            }
+        }
+    }
+}
+
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::initPolygons()
+{
+    mode_ = mtPolygon;
+
+    List<Field<point> > polygons(this->coeffDict().lookup("polygons"));
+    label nPoints = 0;
+    forAll(polygons, polyI)
+    {
+        label np = polygons[polyI].size();
+        if (np < 3)
+        {
+            FatalIOErrorIn
+            (
+                "Foam::ParticleCollector<CloudType>::initPolygons()",
+                this->coeffDict()
+            )
+                << "polygons must consist of at least 3 points"
+                << exit(FatalIOError);
+        }
+
+        nPoints += np;
+    }
+
+    label pointOffset = 0;
+    points_.setSize(nPoints);
+    faces_.setSize(polygons.size());
+    faceTris_.setSize(polygons.size());
+    area_.setSize(polygons.size());
+    forAll(faces_, faceI)
+    {
+        const Field<point>& polyPoints = polygons[faceI];
+        face f(identity(polyPoints.size()) + pointOffset);
+        UIndirectList<point>(points_, f) = polyPoints;
+        area_[faceI] = f.mag(points_);
+
+        DynamicList<face> tris;
+        f.triangles(points_, tris);
+        faceTris_[faceI].transfer(tris);
+
+        faces_[faceI].transfer(f);
+
+        pointOffset += polyPoints.size();
+    }
+}
+
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::initConcentricCircles()
+{
+    mode_ = mtConcentricCircle;
+
+    vector origin(this->coeffDict().lookup("origin"));
+
+    radius_ = this->coeffDict().lookup("radius");
+    nSector_ = readLabel(this->coeffDict().lookup("nSector"));
+ 
+    label nS = nSector_;
+
+    vector refDir;
+    if (nSector_ > 1)
+    {
+        refDir = this->coeffDict().lookup("refDir");
+        refDir -= normal_*(normal_ & refDir);
+        refDir /= mag(refDir);
+    }
+    else
+    {
+        // set 4 quadrants for single sector cases
+        nS = 4;
+
+        vector tangent = vector::zero;
+        scalar magTangent = 0.0;
+
+        Random rnd(1234);
+        while (magTangent < SMALL)
+        {
+            vector v = rnd.vector01();
+
+            tangent = v - (v & normal_)*normal_;
+            magTangent = mag(tangent);
+        }
+
+        refDir = tangent/magTangent;
+    }
+
+    scalar dTheta = 5.0;
+    scalar dThetaSector = 360.0/scalar(nS);
+    label intervalPerSector = max(1, ceil(dThetaSector/dTheta));
+    dTheta = dThetaSector/scalar(intervalPerSector);
+
+    label nPointPerSector = intervalPerSector + 1;
+
+    label nPointPerRadius = nS*(nPointPerSector - 1);
+    label nPoint = radius_.size()*nPointPerRadius;
+    label nFace = radius_.size()*nS;
+
+    // add origin
+    nPoint++;
+
+    points_.setSize(nPoint);
+    faces_.setSize(nFace);
+    area_.setSize(nFace);
+
+    coordSys_ = cylindricalCS("coordSys", origin, normal_, refDir, false);
+
+    List<label> ptIDs(identity(nPointPerRadius));
+
+    points_[0] = origin;
+
+    // points
+    forAll(radius_, radI)
+    {
+        label pointOffset = radI*nPointPerRadius + 1;
+
+        for (label i = 0; i < nPointPerRadius; i++)
+        {
+            label pI = i + pointOffset;
+            point pCyl(radius_[radI], degToRad(i*dTheta), 0.0);
+            points_[pI] = coordSys_.globalPosition(pCyl);
+        }
+    }
+
+    // faces
+    DynamicList<label> facePts(2*nPointPerSector);
+    forAll(radius_, radI)
+    {
+        if (radI == 0)
+        {
+            for (label secI = 0; secI < nS; secI++)
+            {
+                facePts.clear();
+
+                // append origin point
+                facePts.append(0);
+
+                for (label ptI = 0; ptI < nPointPerSector; ptI++)
+                {
+                    label i = ptI + secI*(nPointPerSector - 1);
+                    label id = ptIDs.fcIndex(i - 1) + 1;
+                    facePts.append(id);
+                }
+
+                label faceI = secI + radI*nS;
+
+                faces_[faceI] = face(facePts);
+                area_[faceI] = faces_[faceI].mag(points_);
+            }
+        }
+        else
+        {
+            for (label secI = 0; secI < nS; secI++)
+            {
+                facePts.clear();
+
+                label offset = (radI - 1)*nPointPerRadius + 1;
+
+                for (label ptI = 0; ptI < nPointPerSector; ptI++)
+                {
+                    label i = ptI + secI*(nPointPerSector - 1);
+                    label id = offset + ptIDs.fcIndex(i - 1);
+                    facePts.append(id);
+                }
+                for (label ptI = nPointPerSector-1; ptI >= 0; ptI--)
+                {
+                    label i = ptI + secI*(nPointPerSector - 1);
+                    label id = offset + nPointPerRadius + ptIDs.fcIndex(i - 1);
+                    facePts.append(id);
+                }
+
+                label faceI = secI + radI*nS;
+
+                faces_[faceI] = face(facePts);
+                area_[faceI] = faces_[faceI].mag(points_);
+            }
+        }
+    }
+}
+
+
+template<class CloudType>
+Foam::label Foam::ParticleCollector<CloudType>::collectParcelPolygon
+(
+    const point& position,
+    const vector& U
+) const
+{
+    scalar dt = this->owner().db().time().deltaTValue();
+
+    point end(position + dt*U);
+
+    label dummyNearType = -1;
+    label dummyNearLabel = -1;
+
+    forAll(faces_, faceI)
+    {
+        const label facePoint0 = faces_[faceI][0];
+
+        const point p0 = points_[facePoint0];
+
+        const scalar d1 = normal_ & (position - p0);
+        const scalar d2 = normal_ & (end - p0);
+
+        if (sign(d1) == sign(d2))
+        {
+            // did not cross polygon plane
+            continue;
+        }
+
+        // intersection point
+        point pIntersect = position + (d1/(d1 - d2))*dt*U;
+
+        const List<face>& tris = faceTris_[faceI];
+
+        // identify if point is within poly bounds
+        forAll(tris, triI)
+        {
+            const face& tri = tris[triI];
+            triPointRef t
+            (
+                points_[tri[0]],
+                points_[tri[1]],
+                points_[tri[2]]
+            );
+
+            if (t.classify(pIntersect, dummyNearType, dummyNearLabel))
+            {
+                return faceI;
+            }
+        }
+    }
+
+    return -1;
+}
+
+
+template<class CloudType>
+Foam::label Foam::ParticleCollector<CloudType>::collectParcelConcentricCircles
+(
+    const point& position,
+    const vector& U
+) const
+{
+    label secI = -1;
+
+    scalar dt = this->owner().db().time().deltaTValue();
+
+    point end(position + dt*U);
+
+    const scalar d1 = normal_ & (position - coordSys_.origin());
+    const scalar d2 = normal_ & (end - coordSys_.origin());
+
+    if (sign(d1) == sign(d2))
+    {
+        // did not cross plane
+        return secI;
+    }
+
+    // intersection point in cylindrical co-ordinate system
+    point pCyl = coordSys_.localPosition(position + (d1/(d1 - d2))*dt*U);
+
+    scalar r = pCyl[0];
+
+    if (r < radius_.last())
+    {
+        label radI = 0;
+        while (r > radius_[radI])
+        {
+            radI++;
+        }
+
+        if (nSector_ == 1)
+        {
+            secI = 4*radI;
+        }
+        else
+        {
+            scalar theta = pCyl[1] + constant::mathematical::pi;
+
+            secI =
+                nSector_*radI
+              + floor
+                (
+                    scalar(nSector_)*theta/constant::mathematical::twoPi
+                );
+        }
+    }
+
+    return secI;
+}
+
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::write()
+{
+    const fvMesh& mesh = this->owner().mesh();
+    const Time& time = mesh.time();
+    scalar timeNew = time.value();
+    scalar timeElapsed = timeNew - timeOld_;
+
+    totalTime_ += timeElapsed;
+
+    const scalar alpha = (totalTime_ - timeElapsed)/totalTime_;
+    const scalar beta = timeElapsed/totalTime_;
+
+    forAll(faces_, faceI)
+    {
+        massFlowRate_[faceI] =
+            alpha*massFlowRate_[faceI] + beta*mass_[faceI]/timeElapsed;
+        massTotal_[faceI] += mass_[faceI];
+    }
+
+    const label procI = Pstream::myProcNo();
+
+    Info<< type() << " output:" << nl;
+
+    if (outputFilePtr_.valid())
+    {
+        outputFilePtr_() << time.timeName();
+    }
+
+
+    Field<scalar> faceMassTotal(mass_.size());
+    Field<scalar> faceMassFlowRate(massFlowRate_.size());
+    
+    forAll(faces_, faceI)
+    {
+        scalarList allProcMass(Pstream::nProcs());
+        allProcMass[procI] = massTotal_[faceI];
+        Pstream::gatherList(allProcMass);
+        faceMassTotal[faceI] = sum(allProcMass);
+
+        scalarList allProcMassFlowRate(Pstream::nProcs());
+        allProcMassFlowRate[procI] = massFlowRate_[faceI];
+        Pstream::gatherList(allProcMassFlowRate);
+        faceMassFlowRate[faceI] = sum(allProcMassFlowRate);
+
+        Info<< "    face " << faceI
+            << ": total mass = " << faceMassTotal[faceI]
+            << "; average mass flow rate = " << faceMassFlowRate[faceI]
+            << nl;
+
+        if (outputFilePtr_.valid())
+        {
+            outputFilePtr_()
+                << tab << area_[faceI]
+                << tab << faceMassTotal[faceI]
+                << tab << faceMassFlowRate[faceI]
+                << endl;
+        }
+    }
+
+    Info<< endl;
+
+
+    if (surfaceFormat_ != "none")
+    {
+        if (Pstream::master())
+        {
+            autoPtr<surfaceWriter> writer(surfaceWriter::New(surfaceFormat_));
+
+            writer->write
+            (
+                outputDir_/time.timeName(),
+                "collector",
+                points_,
+                faces_,
+                "massTotal",
+                faceMassTotal,
+                false
+            );
+
+            writer->write
+            (
+                outputDir_/time.timeName(),
+                "collector",
+                points_,
+                faces_,
+                "massFlowRate",
+                faceMassFlowRate,
+                false
+            );
+        }
+    }
+
+
+    if (resetOnWrite_)
+    {
+        forAll(faces_, faceI)
+        {
+            massFlowRate_[faceI] = 0.0;
+        }
+        timeOld_ = timeNew;
+        totalTime_ = 0.0;
+    }
+
+    forAll(faces_, faceI)
+    {
+        mass_[faceI] = 0.0;
+    }
+
+    // writeProperties();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ParticleCollector<CloudType>::ParticleCollector
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    CloudFunctionObject<CloudType>(dict, owner, typeName),
+    mode_(mtUnknown),
+    parcelType_(this->coeffDict().lookupOrDefault("parcelType", -1)),
+    removeCollected_(this->coeffDict().lookup("removeCollected")),
+    points_(),
+    faces_(),
+    faceTris_(),
+    nSector_(0),
+    radius_(),
+    coordSys_(false),
+    normal_(this->coeffDict().lookup("normal")),
+    negateParcelsOppositeNormal_
+    (
+        readBool(this->coeffDict().lookup("negateParcelsOppositeNormal"))
+    ),
+    surfaceFormat_(this->coeffDict().lookup("surfaceFormat")),
+    resetOnWrite_(this->coeffDict().lookup("resetOnWrite")),
+    totalTime_(0.0),
+    mass_(),
+    massTotal_(),
+    massFlowRate_(),
+    log_(this->coeffDict().lookup("log")),
+    outputFilePtr_(),
+    outputDir_(owner.mesh().time().path()),
+    timeOld_(owner.mesh().time().value())
+{
+    if (Pstream::parRun())
+    {
+        // Put in undecomposed case (Note: gives problems for
+        // distributed data running)
+        outputDir_ =
+            outputDir_/".."/"postProcessing"/cloud::prefix/owner.name();
+    }
+    else
+    {
+        outputDir_ =
+            outputDir_/"postProcessing"/cloud::prefix/owner.name();
+    }
+
+    normal_ /= mag(normal_);
+
+    word mode(this->coeffDict().lookup("mode"));
+    if (mode == "polygon")
+    {
+        initPolygons();
+    }
+    else if (mode == "concentricCircle")
+    {
+        initConcentricCircles();
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::ParticleCollector<CloudType>::ParticleCollector"
+            "("
+                "const dictionary& dict,"
+                "CloudType& owner"
+            ")"
+        )
+            << "Unknown mode " << mode << ".  Available options are "
+            << "polygon and concentricCircle" << exit(FatalError);
+    }
+
+    mass_.setSize(faces_.size(), 0.0);
+    massTotal_.setSize(faces_.size(), 0.0);
+    massFlowRate_.setSize(faces_.size(), 0.0);
+
+    makeLogFile(faces_, points_, area_);
+
+    // readProperties(); AND initialise mass... fields
+}
+
+
+template<class CloudType>
+Foam::ParticleCollector<CloudType>::ParticleCollector
+(
+    const ParticleCollector<CloudType>& pc
+)
+:
+    CloudFunctionObject<CloudType>(pc),
+    mode_(pc.mode_),
+    parcelType_(pc.parcelType_),
+    points_(pc.points_),
+    faces_(pc.faces_),
+    faceTris_(pc.faceTris_),
+    nSector_(pc.nSector_),
+    radius_(pc.radius_),
+    coordSys_(pc.coordSys_),
+    normal_(pc.normal_),
+    negateParcelsOppositeNormal_(pc.negateParcelsOppositeNormal_),
+    surfaceFormat_(pc.surfaceFormat_),
+    resetOnWrite_(pc.resetOnWrite_),
+    totalTime_(pc.totalTime_),
+    mass_(pc.mass_),
+    massTotal_(pc.massTotal_),
+    massFlowRate_(pc.massFlowRate_),
+    log_(pc.log_),
+    outputFilePtr_(),
+    outputDir_(pc.outputDir_),
+    timeOld_(0.0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ParticleCollector<CloudType>::~ParticleCollector()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ParticleCollector<CloudType>::postMove
+(
+    const parcelType& p,
+    const label cellI,
+    const scalar dt,
+    bool& keepParticle
+)
+{
+    if ((parcelType_ != -1) && (parcelType_ != p.typeId()))
+    {
+        return;
+    }
+
+    label faceI = -1;
+
+    switch (mode_)
+    {
+        case mtPolygon:
+        {
+            faceI = collectParcelPolygon(p.position(), p.U());
+            break;
+        }
+        case mtConcentricCircle:
+        {
+            faceI = collectParcelConcentricCircles(p.position(), p.U());
+            break;
+        }
+        default:
+        {
+        }
+    }
+
+    if (faceI != -1)
+    {
+        scalar m = p.nParticle()*p.mass();
+
+        if (negateParcelsOppositeNormal_)
+        {
+            vector Uhat = p.U();
+            Uhat /= mag(Uhat) + ROOTVSMALL;
+            if ((Uhat & normal_) < 0)
+            {
+                m *= -1.0;
+            }
+        }
+
+        // add mass contribution
+        mass_[faceI] += m;
+
+        if (nSector_ == 1)
+        {
+            mass_[faceI + 1] += m;
+            mass_[faceI + 2] += m;
+            mass_[faceI + 3] += m;
+        }
+
+        if (removeCollected_)
+        {
+            keepParticle = false;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
new file mode 100644
index 0000000000..283246cfec
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
@@ -0,0 +1,259 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::ParticleCollector
+
+Description
+    Function object to collect the parcel mass- and mass flow rate over a
+    set of polygons.  The polygons are defined as lists of points.  If a
+    parcel is 'collected', it is subsequently flagged to be removed from the
+    domain.
+
+SourceFiles
+    ParticleCollector.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ParticleCollector_H
+#define ParticleCollector_H
+
+#include "CloudFunctionObject.H"
+#include "cylindricalCS.H"
+#include "face.H"
+#include "Switch.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class ParticleCollector Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ParticleCollector
+:
+    public CloudFunctionObject<CloudType>
+{
+public:
+
+    enum modeType
+    {
+        mtPolygon,
+        mtConcentricCircle,
+        mtUnknown
+    };
+
+
+private:
+
+    // Private Data
+
+        // Typedefs
+
+            //- Convenience typedef for parcel type
+            typedef typename CloudType::parcelType parcelType;
+
+        //- Collector mode type
+        modeType mode_;
+
+        //- Index of parcel types to collect (-1 by default = all particles)
+        const label parcelType_;
+
+        //- Flag to remove collected particles
+        Switch removeCollected_;
+
+        //- List of points
+        Field<point> points_;
+
+        //- List of faces
+        List<face> faces_;
+
+
+        // Polygon collector
+
+            //- Triangulation of faces
+            List<List<face> > faceTris_;
+
+        // Concentric circles collector
+
+            //- Number of sectors per circle
+            label nSector_;
+
+            //- List of radii
+            List<scalar> radius_;
+
+            //- Cylindrical co-ordinate system
+            cylindricalCS coordSys_;
+
+
+        //- Face areas
+        Field<scalar> area_;
+
+        //- Polygon normal vector
+        vector normal_;
+
+        //- Remove mass of parcel travelling in opposite direction to normal_
+        bool negateParcelsOppositeNormal_;
+
+        //- Surface output format
+        const word surfaceFormat_;
+
+        //- Flag to indicate whether data should be reset/cleared on writing
+        Switch resetOnWrite_;
+
+        //- Total time
+        scalar totalTime_;
+
+        //- Mass storage
+        List<scalar> mass_;
+
+        //- Mass total storage
+        List<scalar> massTotal_;
+
+        //- Mass flow rate storage
+        List<scalar> massFlowRate_;
+
+        //- Flag to indicate whether data should be written to file
+        Switch log_;
+
+        //- Output file pointer
+        autoPtr<OFstream> outputFilePtr_;
+
+        //- Output directory
+        fileName outputDir_;
+
+        //- Last calculation time
+        scalar timeOld_;
+
+
+    // Private Member Functions
+
+        //- Helper function to create log files
+        void makeLogFile
+        (
+            const faceList& faces,
+            const Field<point>& points,
+            const Field<scalar>& area
+        );
+
+        //- Initialise polygon collectors
+        void initPolygons();
+
+        //- Initialise concentric circle collectors
+        void initConcentricCircles();
+
+        //- Collect parcels in polygon collectors
+        label collectParcelPolygon
+        (
+            const point& position,  
+            const vector& U
+        ) const;
+
+        //- Collect parcels in concentric circle collectors
+        label collectParcelConcentricCircles
+        (
+            const point& position,  
+            const vector& U
+        ) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Write post-processing info
+        void write();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("particleCollector");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        ParticleCollector(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        ParticleCollector(const ParticleCollector<CloudType>& pc);
+
+        //- Construct and return a clone
+        virtual autoPtr<CloudFunctionObject<CloudType> > clone() const
+        {
+            return autoPtr<CloudFunctionObject<CloudType> >
+            (
+                new ParticleCollector<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ParticleCollector();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to the reset on write flag
+            inline const Switch& resetOnWrite() const;
+
+
+        // Evaluation
+
+            //- Post-move hook
+            virtual void postMove
+            (
+                const parcelType& p,
+                const label cellI,
+                const scalar dt,
+                bool& keepParticle
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ParticleCollectorI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "ParticleCollector.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollectorI.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollectorI.H
new file mode 100644
index 0000000000..019deba51a
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollectorI.H
@@ -0,0 +1,34 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+inline const Foam::Switch&
+Foam::ParticleCollector<CloudType>::resetOnWrite() const
+{
+    return resetOnWrite_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
index 81772debfe..c0683d0dbf 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.C
@@ -168,7 +168,8 @@ void Foam::ParticleErosion<CloudType>::postPatch
     const parcelType& p,
     const polyPatch& pp,
     const scalar trackFraction,
-    const tetIndices& tetIs
+    const tetIndices& tetIs,
+    bool&
 )
 {
     const label patchI = pp.index();
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H
index c741bd4004..bedc447188 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleErosion/ParticleErosion.H
@@ -128,7 +128,8 @@ public:
                 const parcelType& p,
                 const polyPatch& pp,
                 const scalar trackFraction,
-                const tetIndices& tetIs
+                const tetIndices& tetIs,
+                bool& keepParticle
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTracks/ParticleTracks.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTracks/ParticleTracks.C
index 9721a5a55d..69f48bf5e0 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTracks/ParticleTracks.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTracks/ParticleTracks.C
@@ -112,7 +112,8 @@ template<class CloudType>
 void Foam::ParticleTracks<CloudType>::postFace
 (
     const parcelType& p,
-    const label
+    const label,
+    bool&
 )
 {
     if
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTracks/ParticleTracks.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTracks/ParticleTracks.H
index 9a9939dbb2..1eeda7bcb8 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTracks/ParticleTracks.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTracks/ParticleTracks.H
@@ -143,7 +143,8 @@ public:
             virtual void postFace
             (
                 const parcelType& p,
-                const label faceI
+                const label faceI,
+                bool& keepParticle
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTrap/ParticleTrap.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTrap/ParticleTrap.C
index f1b19bde59..0a5ff09123 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTrap/ParticleTrap.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTrap/ParticleTrap.C
@@ -104,7 +104,8 @@ void Foam::ParticleTrap<CloudType>::postMove
 (
     parcelType& p,
     const label cellI,
-    const scalar
+    const scalar,
+    bool&
 )
 {
     if (alphaPtr_->internalField()[cellI] < threshold_)
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTrap/ParticleTrap.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTrap/ParticleTrap.H
index 39321af3c0..43fac16c58 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTrap/ParticleTrap.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleTrap/ParticleTrap.H
@@ -125,7 +125,8 @@ public:
             (
                 typename CloudType::parcelType& p,
                 const label cellI,
-                const scalar dt
+                const scalar dt,
+                bool& keepParticle
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C
index f66a872cfc..a267b86266 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C
@@ -220,7 +220,8 @@ void Foam::PatchPostProcessing<CloudType>::postPatch
     const parcelType& p,
     const polyPatch& pp,
     const scalar,
-    const tetIndices& tetIs
+    const tetIndices& tetIs,
+    bool&
 )
 {
     const label patchI = pp.index();
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H
index f7d618e6ce..cb4605e666 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H
@@ -129,7 +129,8 @@ public:
                 const parcelType& p,
                 const polyPatch& pp,
                 const scalar trackFraction,
-                const tetIndices& tetIs
+                const tetIndices& tetIs,
+                bool& keepParticle
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
index dfef9f4835..e5b64f65ca 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -125,7 +125,8 @@ void Foam::VoidFraction<CloudType>::postMove
 (
     const parcelType& p,
     const label cellI,
-    const scalar dt
+    const scalar dt,
+    bool&
 )
 {
     volScalarField& theta = thetaPtr_();
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H
index 26f0f65018..e8ed507f7d 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/VoidFraction/VoidFraction.H
@@ -115,7 +115,8 @@ public:
             (
                 const parcelType& p,
                 const label cellI,
-                const scalar dt
+                const scalar dt,
+                bool& keepParticle
             );
 };
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
index 6315eac827..8804fb6759 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
@@ -302,13 +302,15 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                     if (particleHit)
                     {
-                        this->owner().functions().postFace(p, realFaceI);
+                        bool keep = true;
+                        this->owner().functions().postFace(p, realFaceI, keep);
                         this->owner().functions().postPatch
                         (
                             p,
                             mesh.boundaryMesh()[patchI],
                             1.0,
-                            p.currentTetIndices()
+                            p.currentTetIndices(),
+                            keep
                         );
                      }
                 }
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
index 4eee7d990c..3e1a29e988 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
@@ -398,13 +398,15 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     // surface energy of secondary parcels [J]
     scalar ESigmaSec = 0;
 
-    // sample splash distribution to detrmine secondary parcel diameters
+    // sample splash distribution to determine secondary parcel diameters
     scalarList dNew(parcelsPerSplash_);
+    scalarList npNew(parcelsPerSplash_);
     forAll(dNew, i)
     {
         const scalar y = rndGen_.sample01<scalar>();
         dNew[i] = -dBarSplash*log(exp(-dMin/dBarSplash) - y*K);
-        ESigmaSec += sigma*p.areaS(dNew[i]);
+        npNew[i] = mRatio*np*pow3(d)/pow3(dNew[i])/parcelsPerSplash_;
+        ESigmaSec += npNew[i]*sigma*p.areaS(dNew[i]);
     }
 
     // incident kinetic energy [J]
@@ -459,7 +461,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
         // perturb new parcels towards the owner cell centre
         pPtr->position() += 0.5*rndGen_.sample01<scalar>()*(posC - posCf);
 
-        pPtr->nParticle() = mRatio*np*pow3(d)/pow3(dNew[i])/parcelsPerSplash_;
+        pPtr->nParticle() = npNew[i];
 
         pPtr->d() = dNew[i];
 
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index 5c36e4ad53..3588060c72 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -90,58 +90,6 @@ bool reactingOneDim::read(const dictionary& dict)
 }
 
 
-void reactingOneDim::updateQr()
-{
-    // Retrieve field from coupled region using mapped boundary conditions
-    QrCoupled_.correctBoundaryConditions();
-
-    // Update local Qr from coupled Qr field
-    Qr_ == dimensionedScalar("zero", Qr_.dimensions(), 0.0);
-    forAll(intCoupledPatchIDs_, i)
-    {
-        const label patchI = intCoupledPatchIDs_[i];
-
-        scalarField& Qrp = Qr_.boundaryField()[patchI];
-
-        // Qr is negative going out the solid
-        // If the surface is emitting the radiative flux is set to zero
-        Qrp = max(Qrp, scalar(0.0));
-    }
-
-    const volScalarField kappaRad_(kappaRad());
-
-    // Propagate Qr through 1-D regions
-    label totalFaceId = 0;
-    forAll(intCoupledPatchIDs_, i)
-    {
-        const label patchI = intCoupledPatchIDs_[i];
-
-        const scalarField& Qrp = Qr_.boundaryField()[patchI];
-        const vectorField& Cf = regionMesh().Cf().boundaryField()[patchI];
-
-        forAll(Qrp, faceI)
-        {
-            const scalar Qr0 = Qrp[faceI];
-            point Cf0 = Cf[faceI];
-            const labelList& cells = boundaryFaceCells_[totalFaceId];
-            scalar kappaInt = 0.0;
-            forAll(cells, k)
-            {
-                const label cellI = cells[k];
-                const point& Cf1 = regionMesh().cellCentres()[cellI];
-                const scalar delta = mag(Cf1 - Cf0);
-                kappaInt += kappaRad_[cellI]*delta;
-                Qr_[cellI] = Qr0*exp(-kappaInt);
-                Cf0 = Cf1;
-            }
-            totalFaceId ++;
-        }
-    }
-
-    Qr_.correctBoundaryConditions();
-}
-
-
 void reactingOneDim::updatePhiGas()
 {
     phiHsGas_ ==  dimensionedScalar("zero", phiHsGas_.dimensions(), 0.0);
@@ -198,8 +146,6 @@ void reactingOneDim::updatePhiGas()
 
 void reactingOneDim::updateFields()
 {
-    updateQr();
-
     updatePhiGas();
 }
 
@@ -305,8 +251,6 @@ void reactingOneDim::solveEnergy()
 
     tmp<volScalarField> alpha(solidThermo_.alpha());
 
-    const surfaceScalarField phiQr(fvc::interpolate(Qr_)*nMagSf());
-
     const surfaceScalarField phiGas(fvc::interpolate(phiHsGas_));
 
     fvScalarMatrix hEqn
@@ -315,7 +259,6 @@ void reactingOneDim::solveEnergy()
       - fvm::laplacian(alpha, h_)
      ==
         chemistrySh_
-      + fvc::div(phiQr)
       + fvc::div(phiGas)
     );
 
@@ -380,7 +323,6 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
     ),
     Ys_(solidThermo_.composition().Y()),
     h_(solidThermo_.he()),
-    primaryRadFluxName_(coeffs().lookupOrDefault<word>("radFluxName", "Qr")),
     nNonOrthCorr_(-1),
     maxDiff_(10),
     minimumDelta_(1e-4),
@@ -427,34 +369,6 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
         dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
     ),
 
-    QrCoupled_
-    (
-        IOobject
-        (
-            primaryRadFluxName_,
-            time().timeName(),
-            regionMesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh()
-    ),
-
-    Qr_
-    (
-        IOobject
-        (
-            "QrPyr",
-            time().timeName(),
-            regionMesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh(),
-        dimensionedScalar("zero", dimEnergy/dimArea/dimTime, 0.0),
-        zeroGradientFvPatchVectorField::typeName
-    ),
-
     lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
     addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
     totalGasMassFlux_(0.0),
@@ -492,7 +406,6 @@ reactingOneDim::reactingOneDim
     ),
     Ys_(solidThermo_.composition().Y()),
     h_(solidThermo_.he()),
-    primaryRadFluxName_(dict.lookupOrDefault<word>("radFluxName", "Qr")),
     nNonOrthCorr_(-1),
     maxDiff_(10),
     minimumDelta_(1e-4),
@@ -539,34 +452,6 @@ reactingOneDim::reactingOneDim
         dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
     ),
 
-    QrCoupled_
-    (
-        IOobject
-        (
-            primaryRadFluxName_,
-            time().timeName(),
-            regionMesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh()
-    ),
-
-    Qr_
-    (
-        IOobject
-        (
-            "QrPyr",
-            time().timeName(),
-            regionMesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh(),
-        dimensionedScalar("zero", dimEnergy/dimArea/dimTime, 0.0),
-        zeroGradientFvPatchVectorField::typeName
-    ),
-
     lostSolidMass_(dimensionedScalar("zero", dimMass, 0.0)),
     addedGasMass_(dimensionedScalar("zero", dimMass, 0.0)),
     totalGasMassFlux_(0.0),
@@ -684,36 +569,6 @@ void reactingOneDim::preEvolveRegion()
     {
         solidChemistry_->setCellReacting(cellI, true);
     }
-
-    // De-activate reactions if pyrolysis region coupled to (valid) film
-    if (filmCoupled_)
-    {
-        const volScalarField& filmDelta = filmDeltaPtr_();
-
-        forAll(intCoupledPatchIDs_, i)
-        {
-            const label patchI = intCoupledPatchIDs_[i];
-            const scalarField& filmDeltap = filmDelta.boundaryField()[patchI];
-
-            forAll(filmDeltap, faceI)
-            {
-                const scalar filmDelta0 = filmDeltap[faceI];
-                if (filmDelta0 > reactionDeltaMin_)
-                {
-                    const labelList& cells = boundaryFaceCells_[faceI];
-
-                    // TODO: only limit cell adjacent to film?
-                    //solidChemistry_->setCellNoReacting(cells[0])
-
-                    // Propagate flag through 1-D region
-                    forAll(cells, k)
-                    {
-                        solidChemistry_->setCellReacting(cells[k], false);
-                    }
-                }
-            }
-        }
-    }
 }
 
 
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
index 76cf66a60f..104e56c34e 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -125,16 +125,6 @@ protected:
             volScalarField chemistrySh_;
 
 
-        // Source term fields
-
-            //- Coupled region radiative heat flux [W/m2]
-            //  Requires user to input mapping info for coupled patches
-            volScalarField QrCoupled_;
-
-            //- In depth radiative heat flux [W/m2]
-            volScalarField Qr_;
-
-
         // Checks
 
             //- Cumulative lost mass of the condensed phase [kg]
@@ -164,9 +154,6 @@ protected:
         //- Update/move mesh based on change in mass
         void updateMesh(const scalarField& mass0);
 
-        //- Update radiative flux in pyrolysis region
-        void updateQr();
-
         //- Update enthalpy flux for pyrolysis gases
         void updatePhiGas();
 
@@ -259,12 +246,6 @@ public:
             virtual scalar solidRegionDiffNo() const;
 
 
-        // Source fields (read/write access)
-
-            //- In depth radiative heat flux
-            inline const volScalarField& Qr() const;
-
-
        // Evolution
 
             //- Pre-evolve region
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDimI.H b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDimI.H
index eec1e7fa0a..375c378080 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDimI.H
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDimI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,11 +34,4 @@ Foam::regionModels::pyrolysisModels::reactingOneDim::nNonOrthCorr() const
 }
 
 
-inline const Foam::volScalarField&
-Foam::regionModels::pyrolysisModels::reactingOneDim::Qr() const
-{
-    return Qr_;
-}
-
-
 // ************************************************************************* //
diff --git a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.C b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.C
index b0f8d1629a..1ffe2ae82f 100644
--- a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.C
+++ b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,8 +40,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(p, iF),
     phiName_("phi"),
-    rhoName_("rho"),
-    deltaWet_(1e-6)
+    rhoName_("rho")
 {}
 
 
@@ -56,8 +55,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
-    deltaWet_(ptf.deltaWet_)
+    rhoName_(ptf.rhoName_)
 {}
 
 
@@ -71,8 +69,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(p, iF),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
-    deltaWet_(dict.lookupOrDefault<scalar>("deltaWet", 1e-6))
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
 {
     fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
 }
@@ -86,8 +83,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(fptpsf),
     phiName_(fptpsf.phiName_),
-    rhoName_(fptpsf.rhoName_),
-    deltaWet_(fptpsf.deltaWet_)
+    rhoName_(fptpsf.rhoName_)
 {}
 
 
@@ -100,8 +96,7 @@ filmPyrolysisTemperatureCoupledFvPatchScalarField
 :
     fixedValueFvPatchScalarField(fptpsf, iF),
     phiName_(fptpsf.phiName_),
-    rhoName_(fptpsf.rhoName_),
-    deltaWet_(fptpsf.deltaWet_)
+    rhoName_(fptpsf.rhoName_)
 {}
 
 
@@ -151,13 +146,12 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
-    scalarField deltaFilm = filmModel.delta().boundaryField()[filmPatchI];
-    filmModel.toPrimary(filmPatchI, deltaFilm);
+    scalarField alphaFilm = filmModel.alpha().boundaryField()[filmPatchI];
+    filmModel.toPrimary(filmPatchI, alphaFilm);
 
     scalarField TFilm = filmModel.Ts().boundaryField()[filmPatchI];
     filmModel.toPrimary(filmPatchI, TFilm);
 
-
     // Retrieve pyrolysis model
     const pyrModelType& pyrModel =
         db().lookupObject<pyrModelType>("pyrolysisProperties");
@@ -168,19 +162,8 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::updateCoeffs()
     pyrModel.toPrimary(pyrPatchI, TPyr);
 
 
-    forAll(deltaFilm, i)
-    {
-        if (deltaFilm[i] > deltaWet_)
-        {
-            // temperature set by film
-            Tp[i] = TFilm[i];
-        }
-        else
-        {
-            // temperature set by pyrolysis model
-            Tp[i] = TPyr[i];
-        }
-    }
+    // Evaluate temperature
+    Tp = alphaFilm*TFilm + (1.0 - alphaFilm)*TPyr;
 
     // Restore tag
     UPstream::msgType() = oldTag;
@@ -197,7 +180,6 @@ void Foam::filmPyrolysisTemperatureCoupledFvPatchScalarField::write
     fvPatchScalarField::write(os);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    os.writeKeyword("deltaWet") << deltaWet_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 
diff --git a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.H b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.H
index 503b023f95..bc4ace8abd 100644
--- a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.H
+++ b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisTemperatureCoupled/filmPyrolysisTemperatureCoupledFvPatchScalarField.H
@@ -28,15 +28,7 @@ Description
     This boundary condition is designed to be used in conjunction with surface
     film and pyrolysis modelling.  It provides a temperature boundary condition
     for patches on the primary region based on whether the patch is seen to
-    be 'wet', specified by:
-
-        \f[
-            delta > delta_wet
-        \f]
-
-    where
-        \var delta     = film height [m]
-        \var delta_wet = film height above which the surface is considered wet
+    be 'wet', retrieved from the film alpha field.
 
     \li if the patch is wet, the temperature is set using the film temperature
     \li otherwise, it is set using pyrolysis temperature
@@ -84,9 +76,6 @@ class filmPyrolysisTemperatureCoupledFvPatchScalarField
         //- Name of density field
         word rhoName_;
 
-        //- Film height threshold beyond which it is considered 'wet' [m]
-        scalar deltaWet_;
-
 
 public:
 
diff --git a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C
index 29b3e34586..9210d30e57 100644
--- a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C
+++ b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,8 +40,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(p, iF),
     phiName_("phi"),
-    rhoName_("rho"),
-    deltaWet_(1e-6)
+    rhoName_("rho")
 {}
 
 
@@ -56,8 +55,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(ptf, p, iF, mapper),
     phiName_(ptf.phiName_),
-    rhoName_(ptf.rhoName_),
-    deltaWet_(ptf.deltaWet_)
+    rhoName_(ptf.rhoName_)
 {}
 
 
@@ -71,8 +69,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(p, iF),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
-    deltaWet_(dict.lookupOrDefault<scalar>("deltaWet", 1e-6))
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
 {
     fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
 }
@@ -86,8 +83,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(fpvpvf),
     phiName_(fpvpvf.phiName_),
-    rhoName_(fpvpvf.rhoName_),
-    deltaWet_(fpvpvf.deltaWet_)
+    rhoName_(fpvpvf.rhoName_)
 {}
 
 
@@ -100,8 +96,7 @@ filmPyrolysisVelocityCoupledFvPatchVectorField
 :
     fixedValueFvPatchVectorField(fpvpvf, iF),
     phiName_(fpvpvf.phiName_),
-    rhoName_(fpvpvf.rhoName_),
-    deltaWet_(fpvpvf.deltaWet_)
+    rhoName_(fpvpvf.rhoName_)
 {}
 
 
@@ -154,13 +149,12 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
 
     const label filmPatchI = filmModel.regionPatchID(patchI);
 
-    scalarField deltaFilm = filmModel.delta().boundaryField()[filmPatchI];
-    filmModel.toPrimary(filmPatchI, deltaFilm);
+    scalarField alphaFilm = filmModel.alpha().boundaryField()[filmPatchI];
+    filmModel.toPrimary(filmPatchI, alphaFilm);
 
     vectorField UFilm = filmModel.Us().boundaryField()[filmPatchI];
     filmModel.toPrimary(filmPatchI, UFilm);
 
-
     // Retrieve pyrolysis model
     const pyrModelType& pyrModel =
         db().objectRegistry::lookupObject<pyrModelType>
@@ -203,19 +197,9 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::updateCoeffs()
     const scalarField UAvePyr(-phiPyr/patch().magSf());
     const vectorField& nf = patch().nf();
 
-    forAll(deltaFilm, i)
-    {
-        if (deltaFilm[i] > deltaWet_)
-        {
-            // velocity set by film
-            Up[i] = UFilm[i];
-        }
-        else
-        {
-            // velocity set by pyrolysis model
-            Up[i] = UAvePyr[i]*nf[i];
-        }
-    }
+
+    // Evaluate velocity
+    Up = alphaFilm*UFilm + (1.0 - alphaFilm)*UAvePyr*nf;
 
     // Restore tag
     UPstream::msgType() = oldTag;
@@ -232,7 +216,6 @@ void Foam::filmPyrolysisVelocityCoupledFvPatchVectorField::write
     fvPatchVectorField::write(os);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    os.writeKeyword("deltaWet") << deltaWet_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 
diff --git a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.H b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.H
index c7c4c5d24b..b37157ddf0 100644
--- a/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.H
+++ b/src/regionModels/regionCoupling/derivedFvPatchFields/filmPyrolysisVelocityCoupled/filmPyrolysisVelocityCoupledFvPatchVectorField.H
@@ -28,15 +28,7 @@ Description
     This boundary condition is designed to be used in conjunction with surface
     film and pyrolysis modelling.  It provides a velocity boundary condition
     for patches on the primary region based on whether the patch is seen to
-    be 'wet', specified by:
-
-        \f[
-            delta > delta_wet
-        \f]
-
-    where
-        \var delta     = film height [m]
-        \var delta_wet = film height above which the surface is considered wet
+    be 'wet', retrieved from the film alpha field.
 
     \li if the patch is wet, the velocity is set using the film velocity
     \li otherwise, it is set using pyrolysis out-gassing velocity
@@ -84,9 +76,6 @@ class filmPyrolysisVelocityCoupledFvPatchVectorField
         //- Name of density field
         word rhoName_;
 
-        //- Film height threshold beyond which it is considered 'wet'
-        scalar deltaWet_;
-
 
 public:
 
diff --git a/src/regionModels/regionModel/Make/files b/src/regionModels/regionModel/Make/files
index ab18530890..11002e0a29 100644
--- a/src/regionModels/regionModel/Make/files
+++ b/src/regionModels/regionModel/Make/files
@@ -8,4 +8,8 @@ derivedFvPatches/mappedVariableThicknessWall/mappedVariableThicknessWallFvPatch.
 
 regionProperties/regionProperties.C
 
+regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.C
+regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectNew.C
+regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.C
+
 LIB = $(FOAM_LIBBIN)/libregionModels
diff --git a/src/regionModels/regionModel/regionModel/regionModel.C b/src/regionModels/regionModel/regionModel/regionModel.C
index 98e152c3f1..227af0adb4 100644
--- a/src/regionModels/regionModel/regionModel/regionModel.C
+++ b/src/regionModels/regionModel/regionModel/regionModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -190,6 +190,66 @@ bool Foam::regionModels::regionModel::read(const dictionary& dict)
 }
 
 
+Foam::label Foam::regionModels::regionModel::nbrCoupledPatchID
+(
+    const regionModel& nbrRegion,
+    const label regionPatchI
+) const
+{
+    label nbrPatchI = -1;
+
+    // region
+    const fvMesh& nbrRegionMesh = nbrRegion.regionMesh();
+
+    // boundary mesh
+    const polyBoundaryMesh& nbrPbm = nbrRegionMesh.boundaryMesh();
+
+    const mappedPatchBase& mpb =
+        refCast<const mappedPatchBase>
+        (
+            regionMesh().boundaryMesh()[regionPatchI]
+        );
+
+    // sample patch name on the primary region
+    const word& primaryPatchName = mpb.samplePatch();
+
+    // find patch on nbr region that has the same sample patch name
+    forAll(nbrRegion.intCoupledPatchIDs(), j)
+    {
+        const label nbrRegionPatchI = nbrRegion.intCoupledPatchIDs()[j];
+
+        const mappedPatchBase& mpb =
+            refCast<const mappedPatchBase>(nbrPbm[nbrRegionPatchI]);
+
+        if (mpb.samplePatch() == primaryPatchName)
+        {
+            nbrPatchI = nbrRegionPatchI;
+            break;
+        }
+    }
+
+    if (nbrPatchI == -1)
+    {
+        const polyPatch& p = regionMesh().boundaryMesh()[regionPatchI];
+
+        FatalErrorIn
+        (
+            "Foam::tmp<Foam::Field<Type> > "
+            "Foam::regionModels::regionModel::nbrCoupledPatchID"
+            "("
+                "const regionModel& , "
+                "const label"
+            ") const"
+        )
+            << "Unable to find patch pair for local patch "
+            << p.name() << " and region " << nbrRegion.name()
+            << abort(FatalError);
+    }
+
+    return nbrPatchI;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::regionModels::regionModel::regionModel(const fvMesh& mesh)
@@ -214,7 +274,8 @@ Foam::regionModels::regionModel::regionModel(const fvMesh& mesh)
     coeffs_(dictionary::null),
     primaryPatchIDs_(),
     intCoupledPatchIDs_(),
-    regionName_("none")
+    regionName_("none"),
+    functions_(*this)
 {}
 
 
@@ -246,7 +307,8 @@ Foam::regionModels::regionModel::regionModel
     coeffs_(subOrEmptyDict(modelName + "Coeffs")),
     primaryPatchIDs_(),
     intCoupledPatchIDs_(),
-    regionName_(lookup("regionName"))
+    regionName_(lookup("regionName")),
+    functions_(*this, subOrEmptyDict("functions"))
 {
     if (active_)
     {
@@ -274,7 +336,7 @@ Foam::regionModels::regionModel::regionModel
     (
         IOobject
         (
-            regionType,
+            regionType + "Properties",
             mesh.time().constant(),
             mesh,
             IOobject::NO_READ,
@@ -292,7 +354,8 @@ Foam::regionModels::regionModel::regionModel
     coeffs_(dict.subOrEmptyDict(modelName + "Coeffs")),
     primaryPatchIDs_(),
     intCoupledPatchIDs_(),
-    regionName_(dict.lookup("regionName"))
+    regionName_(dict.lookup("regionName")),
+    functions_(*this, subOrEmptyDict("functions"))
 {
     if (active_)
     {
@@ -315,18 +378,6 @@ Foam::regionModels::regionModel::~regionModel()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::regionModels::regionModel::preEvolveRegion()
-{
-    // do nothing
-}
-
-
-void Foam::regionModels::regionModel::evolveRegion()
-{
-    // do nothing
-}
-
-
 void Foam::regionModels::regionModel::evolve()
 {
     if (active_)
@@ -334,15 +385,14 @@ void Foam::regionModels::regionModel::evolve()
         Info<< "\nEvolving " << modelName_ << " for region "
             << regionMesh().name() << endl;
 
-        // Update any input information
         //read();
 
-        // Pre-evolve
         preEvolveRegion();
 
-        // Increment the region equations up to the new time level
         evolveRegion();
 
+        postEvolveRegion();
+
         // Provide some feedback
         if (infoOutput_)
         {
@@ -354,6 +404,24 @@ void Foam::regionModels::regionModel::evolve()
 }
 
 
+void Foam::regionModels::regionModel::preEvolveRegion()
+{
+    functions_.preEvolveRegion();
+}
+
+
+void Foam::regionModels::regionModel::evolveRegion()
+{
+    // do nothing
+}
+
+
+void Foam::regionModels::regionModel::postEvolveRegion()
+{
+    functions_.postEvolveRegion();
+}
+
+
 void Foam::regionModels::regionModel::info() const
 {
     // do nothing
diff --git a/src/regionModels/regionModel/regionModel/regionModel.H b/src/regionModels/regionModel/regionModel/regionModel.H
index 6fbad3d372..9f2ec23e67 100644
--- a/src/regionModels/regionModel/regionModel/regionModel.H
+++ b/src/regionModels/regionModel/regionModel/regionModel.H
@@ -41,16 +41,13 @@ SourceFiles
 #include "labelList.H"
 #include "volFields.H"
 #include "mappedPatchBase.H"
+#include "regionModelFunctionObjectList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward declaration of classes
-//class fvMesh;
-//class Time;
-
 namespace regionModels
 {
 
@@ -121,6 +118,9 @@ protected:
         //- Region name
         word regionName_;
 
+        //- Region model function objects
+        regionModelFunctionObjectList functions_;
+
 
     // Protected member functions
 
@@ -219,6 +219,46 @@ public:
 
         // Helper
 
+            //- Return the coupled patch ID paired with coupled patch
+            //  regionPatchI
+            label nbrCoupledPatchID
+            (
+                const regionModel& nbrRegion,
+                const label regionPatchI
+            ) const;
+
+            //- Map patch field from another region model to local patch
+            template<class Type>
+            tmp<Foam::Field<Type> > mapRegionPatchField
+            (
+                const regionModel& nbrRegion,
+                const label regionPatchI,
+                const label nbrPatchI,
+                const Field<Type>& nbrField,
+                const bool flip = false
+            ) const;
+
+            //- Map patch field from another region model to local patch
+            template<class Type>
+            tmp<Field<Type> > mapRegionPatchField
+            (
+                const word& regionModelName,
+                const word& fieldName,
+                const label regionPatchI,
+                const bool flip = false
+            ) const;
+
+            //- Map patch internal field from another region model to local
+            //  patch
+            template<class Type>
+            tmp<Field<Type> > mapRegionPatchInternalField
+            (
+                const word& regionModelName,
+                const word& fieldName,
+                const label regionPatchI,
+                const bool flip = false
+            ) const;
+
             //- Convert a local region field to the primary region
             template<class Type>
             void toPrimary
@@ -256,14 +296,17 @@ public:
 
         // Evolution
 
+            //- Main driver routing to evolve the region - calls other evolves
+            virtual void evolve();
+
             //- Pre-evolve region
             virtual void preEvolveRegion();
 
             //- Evolve the region
             virtual void evolveRegion();
 
-            //- Evolve the region
-            virtual void evolve();
+            //- Post-evolve region
+            virtual void postEvolveRegion();
 
 
         // I-O
diff --git a/src/regionModels/regionModel/regionModel/regionModelTemplates.C b/src/regionModels/regionModel/regionModel/regionModelTemplates.C
index b78772de01..4e98b2e2fb 100644
--- a/src/regionModels/regionModel/regionModel/regionModelTemplates.C
+++ b/src/regionModels/regionModel/regionModel/regionModelTemplates.C
@@ -23,6 +23,164 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::regionModels::regionModel::mapRegionPatchField
+(
+    const regionModel& nbrRegion,
+    const label regionPatchI,
+    const label nbrPatchI,
+    const Field<Type>& nbrField,
+    const bool flip
+) const
+{
+    const fvMesh& nbrRegionMesh = nbrRegion.regionMesh();
+
+    const polyPatch& p = regionMesh().boundaryMesh()[regionPatchI];
+
+    const polyPatch& nbrP = nbrRegionMesh.boundaryMesh()[nbrPatchI];
+
+    int oldTag = UPstream::msgType();
+    UPstream::msgType() = oldTag + 1;
+
+    AMIPatchToPatchInterpolation ami(p, nbrP, faceAreaIntersect::tmMesh, flip);
+
+    tmp<Field<Type> > tresult(ami.interpolateToSource(nbrField));
+
+    UPstream::msgType() = oldTag;
+
+    return tresult;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::regionModels::regionModel::mapRegionPatchField
+(
+    const word& regionModelName,
+    const word& fieldName,
+    const label regionPatchI,
+    const bool flip
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    const regionModel& nbrRegion =
+        this->primaryMesh_.lookupObject<regionModel>(regionModelName);
+
+    const fvMesh& nbrRegionMesh = nbrRegion.regionMesh();
+
+    const polyPatch& p = regionMesh().boundaryMesh()[regionPatchI];
+
+    if (nbrRegionMesh.foundObject<fieldType>(fieldName))
+    {
+        const fieldType& nbrField =
+            nbrRegionMesh.lookupObject<fieldType>(fieldName);
+
+        const label nbrPatchI = nbrCoupledPatchID(nbrRegion, regionPatchI);
+
+        const polyPatch& nbrP = nbrRegionMesh.boundaryMesh()[nbrPatchI];
+
+        int oldTag = UPstream::msgType();
+        UPstream::msgType() = oldTag + 1;
+
+        AMIPatchToPatchInterpolation ami
+        (
+            p,
+            nbrP,
+            faceAreaIntersect::tmMesh,
+            flip
+        );
+
+        const Field<Type>& nbrFieldp = nbrField.boundaryField()[nbrPatchI];
+
+        tmp<Field<Type> > tresult(ami.interpolateToSource(nbrFieldp));
+
+        UPstream::msgType() = oldTag;
+
+        return tresult;
+    }
+    else
+    {
+        return
+            tmp<Field<Type> >
+            (
+                new Field<Type>
+                (
+                    p.size(),
+                    pTraits<Type>::zero
+                )
+            );
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::regionModels::regionModel::mapRegionPatchInternalField
+(
+    const word& regionModelName,
+    const word& fieldName,
+    const label regionPatchI,
+    const bool flip
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    const regionModel& nbrRegion =
+        this->primaryMesh_.lookupObject<regionModel>(regionModelName);
+
+    const fvMesh& nbrRegionMesh = nbrRegion.regionMesh();
+
+    const polyPatch& p = regionMesh().boundaryMesh()[regionPatchI];
+
+    if (nbrRegionMesh.foundObject<fieldType>(fieldName))
+    {
+        const fieldType& nbrField =
+            nbrRegionMesh.lookupObject<fieldType>(fieldName);
+
+        const label nbrPatchI = nbrCoupledPatchID(nbrRegion, regionPatchI);
+
+        const polyPatch& nbrP = nbrRegionMesh.boundaryMesh()[nbrPatchI];
+
+        int oldTag = UPstream::msgType();
+        UPstream::msgType() = oldTag + 1;
+
+        AMIPatchToPatchInterpolation ami
+        (
+            p,
+            nbrP,
+            faceAreaIntersect::tmMesh,
+            flip
+        );
+
+        const fvPatchField<Type>& nbrFieldp =
+            nbrField.boundaryField()[nbrPatchI];
+
+        tmp<Field<Type> > tresult
+        (
+            ami.interpolateToSource(nbrFieldp.patchInternalField())
+        );
+
+        UPstream::msgType() = oldTag;
+
+        return tresult;
+    }
+    else
+    {
+        return
+            tmp<Field<Type> >
+            (
+                new Field<Type>
+                (
+                    p.size(),
+                    pTraits<Type>::zero
+                )
+            );
+    }
+}
+
+
 template<class Type>
 void Foam::regionModels::regionModel::toPrimary
 (
diff --git a/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.C b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.C
new file mode 100644
index 0000000000..182258b637
--- /dev/null
+++ b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.C
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "regionModelFunctionObject.H"
+#include "regionModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+    defineTypeNameAndDebug(regionModelFunctionObject, 0);
+    defineRunTimeSelectionTable(regionModelFunctionObject, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::regionModels::regionModelFunctionObject::regionModelFunctionObject
+(
+    regionModel& owner
+)
+:
+    dict_(dictionary::null),
+    owner_(owner),
+    modelType_("modelType")
+{}
+
+
+Foam::regionModels::regionModelFunctionObject::regionModelFunctionObject
+(
+    const dictionary& dict,
+    regionModel& owner,
+    const word& type
+)
+:
+    dict_(dict),
+    owner_(owner),
+    modelType_(type)
+{}
+
+
+Foam::regionModels::regionModelFunctionObject::regionModelFunctionObject
+(
+    const regionModelFunctionObject& rmfo
+)
+:
+    dict_(rmfo.dict_),
+    owner_(rmfo.owner_),
+    modelType_(rmfo.modelType_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::regionModels::regionModelFunctionObject::~regionModelFunctionObject()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::regionModels::regionModelFunctionObject::preEvolveRegion()
+{
+    // do nothing
+}
+
+
+void Foam::regionModels::regionModelFunctionObject::postEvolveRegion()
+{
+    if (owner_.regionMesh().time().outputTime())
+    {
+        write();
+    }
+}
+
+
+void Foam::regionModels::regionModelFunctionObject::write() const
+{
+    // do nothing
+}
+
+
+// ************************************************************************* //
diff --git a/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.H b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.H
new file mode 100644
index 0000000000..77e56b14d1
--- /dev/null
+++ b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObject.H
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::regionModelFunctionObject
+
+Description
+    Region model function object base class
+
+SourceFiles
+    regionModelFunctionObject.C
+    regionModelFunctionObjectNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionModelFunctionObject_H
+#define regionModelFunctionObject_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace regionModels
+{
+
+class regionModel;
+
+/*---------------------------------------------------------------------------*\
+                  Class regionModelFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionModelFunctionObject
+{
+protected:
+
+    // Protected data
+
+        //- Dictionary
+        dictionary dict_;
+
+        //- Reference to the region model
+        regionModel& owner_;
+
+        //- Model type name
+        word modelType_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("regionModelFunctionObject");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        regionModelFunctionObject,
+        dictionary,
+        (
+            const dictionary& dict,
+            regionModel& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        regionModelFunctionObject(regionModel& owner);
+
+        //- Construct from dictionary
+        regionModelFunctionObject
+        (
+            const dictionary& dict,
+            regionModel& owner,
+            const word& modelType
+        );
+
+        //- Construct copy
+        regionModelFunctionObject(const regionModelFunctionObject& ppm);
+
+        //- Construct and return a clone
+        virtual autoPtr<regionModelFunctionObject> clone() const
+        {
+            return autoPtr<regionModelFunctionObject>
+            (
+                new regionModelFunctionObject(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~regionModelFunctionObject();
+
+
+    //- Selector
+    static autoPtr<regionModelFunctionObject> New
+    (
+        const dictionary& dict,
+        regionModel& owner,
+        const word& modelType
+    );
+
+
+    // Member Functions
+
+        // Evaluation
+
+            //- Pre-evolve region hook
+            virtual void preEvolveRegion();
+
+            //- Post-evolve region hook
+            virtual void postEvolveRegion();
+
+    // I-O
+
+        //- write
+        virtual void write() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.C b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.C
new file mode 100644
index 0000000000..366ec496ba
--- /dev/null
+++ b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "regionModelFunctionObjectList.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::regionModels::regionModelFunctionObjectList::regionModelFunctionObjectList
+(
+    regionModel& owner
+)
+:
+    PtrList<regionModelFunctionObject>(),
+    owner_(owner),
+    dict_(dictionary::null)
+{}
+
+
+Foam::regionModels::regionModelFunctionObjectList::regionModelFunctionObjectList
+(
+    regionModel& owner,
+    const dictionary& dict,
+    const bool readFields
+)
+:
+    PtrList<regionModelFunctionObject>(),
+    owner_(owner),
+    dict_(dict)
+{
+    if (readFields)
+    {
+        wordList modelNames(dict.toc());
+
+        Info<< "    Selecting region model functions" << endl;
+
+        if (modelNames.size() > 0)
+        {
+            this->setSize(modelNames.size());
+
+            forAll(modelNames, i)
+            {
+                const word& modelName = modelNames[i];
+
+                this->set
+                (
+                    i,
+                    regionModelFunctionObject::New
+                    (
+                        dict,
+                        owner,
+                        modelName
+                    )
+                );
+            }
+        }
+        else
+        {
+            Info<< "    none" << endl;
+        }
+    }
+}
+
+
+Foam::regionModels::regionModelFunctionObjectList::regionModelFunctionObjectList
+(
+    const regionModelFunctionObjectList& cfol
+)
+:
+    PtrList<regionModelFunctionObject>(cfol),
+    owner_(cfol.owner_),
+    dict_(cfol.dict_)
+{}
+
+
+// * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
+
+Foam::regionModels::regionModelFunctionObjectList::
+~regionModelFunctionObjectList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::regionModels::regionModelFunctionObjectList::preEvolveRegion()
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).preEvolveRegion();
+    }
+}
+
+
+void Foam::regionModels::regionModelFunctionObjectList::postEvolveRegion()
+{
+    forAll(*this, i)
+    {
+        this->operator[](i).postEvolveRegion();
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.H b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.H
new file mode 100644
index 0000000000..22488f18c8
--- /dev/null
+++ b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectList.H
@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::regionModelFunctionObjectList
+
+Description
+    List of cloud function objects
+
+SourceFiles
+    regionModelFunctionObjectListI.H
+    regionModelFunctionObjectList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionModelFunctionObjectList_H
+#define regionModelFunctionObjectList_H
+
+#include "PtrList.H"
+#include "regionModelFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+
+class regionModel;
+
+/*---------------------------------------------------------------------------*\
+                Class regionModelFunctionObjectList Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionModelFunctionObjectList
+:
+    public PtrList<regionModelFunctionObject>
+{
+protected:
+
+    // Protected Data
+
+        //- Reference to the owner region model
+        regionModel& owner_;
+
+        //- Dictionary
+        const dictionary dict_;
+
+
+public:
+
+    // Constructors
+
+        //- Null constructor
+        regionModelFunctionObjectList(regionModel& owner);
+
+        //- Construct from mesh
+        regionModelFunctionObjectList
+        (
+            regionModel& owner,
+            const dictionary& dict,
+            const bool readFields = true
+        );
+
+        //- Construct copy
+        regionModelFunctionObjectList
+        (
+            const regionModelFunctionObjectList& rmfol
+        );
+
+
+    //- Destructor
+    virtual ~regionModelFunctionObjectList();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to the cloud owner
+            inline const regionModel& owner() const;
+
+            //- Return refernce to the cloud owner
+            inline regionModel& owner();
+
+            //- Return the forces dictionary
+            inline const dictionary& dict() const;
+
+
+        // Evaluation
+
+            //- Pre-evolve hook
+            virtual void preEvolveRegion();
+
+            //- Post-evolve hook
+            virtual void postEvolveRegion();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "regionModelFunctionObjectListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectListI.H b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectListI.H
new file mode 100644
index 0000000000..a0a90d0f67
--- /dev/null
+++ b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectListI.H
@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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/>.
+
+\*---------------------------------------------------------------------------*/
+
+inline const Foam::regionModels::regionModel&
+Foam::regionModels::regionModelFunctionObjectList::owner() const
+{
+    return owner_;
+}
+
+
+inline Foam::regionModels::regionModel&
+Foam::regionModels::regionModelFunctionObjectList::owner()
+{
+    return owner_;
+}
+
+
+inline const Foam::dictionary&
+Foam::regionModels::regionModelFunctionObjectList::dict() const
+{
+    return dict_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectNew.C b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectNew.C
new file mode 100644
index 0000000000..11371b1f79
--- /dev/null
+++ b/src/regionModels/regionModel/regionModelFunctionObject/regionModelFunctionObject/regionModelFunctionObjectNew.C
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "regionModelFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::regionModels::regionModelFunctionObject>
+Foam::regionModels::regionModelFunctionObject::New
+(
+    const dictionary& dict,
+    regionModel& owner,
+    const word& modelName
+)
+{
+    const word modelType = dict.subDict(modelName).lookup("type");
+
+    Info<< "        " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "regionModelFunctionObject::New"
+            "("
+                "const dictionary&, "
+                "regionModel&, "
+                "const word&"
+            ")"
+        )   << "Unknown region model function type "
+            << modelType << nl << nl
+            << "Valid region model function types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return
+        autoPtr<regionModelFunctionObject>
+        (
+            cstrIter()
+            (
+                dict.subDict(modelName),
+                owner
+            )
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/Make/files b/src/regionModels/surfaceFilmModels/Make/files
index 12f80a7ce7..1151e2f7fd 100644
--- a/src/regionModels/surfaceFilmModels/Make/files
+++ b/src/regionModels/surfaceFilmModels/Make/files
@@ -38,12 +38,16 @@ $(THERMOMODELS)/heatTransferModel/mappedConvectiveHeatTransfer/mappedConvectiveH
 $(THERMOMODELS)/filmRadiationModel/filmRadiationModel/filmRadiationModel.C
 $(THERMOMODELS)/filmRadiationModel/filmRadiationModel/filmRadiationModelNew.C
 $(THERMOMODELS)/filmRadiationModel/noRadiation/noRadiation.C
+$(THERMOMODELS)/filmRadiationModel/constantRadiation/constantRadiation.C
+$(THERMOMODELS)/filmRadiationModel/primaryRadiation/primaryRadiation.C
 $(THERMOMODELS)/filmRadiationModel/standardRadiation/standardRadiation.C
 
 
 /* Boundary conditions */
 PATCHFIELDS=derivedFvPatchFields
 $(PATCHFIELDS)/filmHeightInletVelocity/filmHeightInletVelocityFvPatchVectorField.C
+$(PATCHFIELDS)/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.C
+$(PATCHFIELDS)/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.C
 
 /* Wall functions for primary region */
 WALLFUNCS=$(PATCHFIELDS)/wallFunctions
diff --git a/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.C b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.C
new file mode 100644
index 0000000000..b36832a1d3
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.C
@@ -0,0 +1,213 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "inclinedFilmNusseltHeightFvPatchScalarField.H"
+#include "volFields.H"
+#include "kinematicSingleLayer.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    GammaMean_(),
+    a_(),
+    omega_()
+{}
+
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const inclinedFilmNusseltHeightFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    GammaMean_(ptf.GammaMean_().clone().ptr()),
+    a_(ptf.a_().clone().ptr()),
+    omega_(ptf.omega_().clone().ptr())
+{}
+
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    GammaMean_(DataEntry<scalar>::New("GammaMean", dict)),
+    a_(DataEntry<scalar>::New("a", dict)),
+    omega_(DataEntry<scalar>::New("omega", dict))
+{
+    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
+}
+
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const inclinedFilmNusseltHeightFvPatchScalarField& wmfrhpsf
+)
+:
+    fixedValueFvPatchScalarField(wmfrhpsf),
+    GammaMean_(wmfrhpsf.GammaMean_().clone().ptr()),
+    a_(wmfrhpsf.a_().clone().ptr()),
+    omega_(wmfrhpsf.omega_().clone().ptr())
+{}
+
+
+Foam::inclinedFilmNusseltHeightFvPatchScalarField::
+inclinedFilmNusseltHeightFvPatchScalarField
+(
+    const inclinedFilmNusseltHeightFvPatchScalarField& wmfrhpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(wmfrhpsf, iF),
+    GammaMean_(wmfrhpsf.GammaMean_().clone().ptr()),
+    a_(wmfrhpsf.a_().clone().ptr()),
+    omega_(wmfrhpsf.omega_().clone().ptr())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::inclinedFilmNusseltHeightFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const label patchI = patch().index();
+
+    const scalar t = db().time().timeOutputValue();
+
+    // retrieve the film region from the database
+
+    const regionModels::regionModel& region =
+        db().lookupObject<regionModels::regionModel>("surfaceFilmProperties");
+
+    const regionModels::surfaceFilmModels::kinematicSingleLayer& film =
+        dynamic_cast
+        <
+            const regionModels::surfaceFilmModels::kinematicSingleLayer&
+        >(region);
+
+    // calculate the vector tangential to the patch
+
+    const vectorField n(patch().nf());
+
+    const volVectorField& nHat = film.nHat();
+
+    const vectorField nHatp(nHat.boundaryField()[patchI].patchInternalField());
+
+    vectorField nTan(nHatp ^ n);
+    nTan /= mag(nTan) + ROOTVSMALL;
+
+    // calculate distance in patch tangential direction
+
+    const vectorField& Cf = patch().Cf();
+    scalarField d(nTan & Cf);
+
+    // calculate the wavy film height
+
+    const scalar GMean = GammaMean_->value(t);
+    const scalar a = a_->value(t);
+    const scalar omega = omega_->value(t);
+
+    const scalarField G(GMean + a*sin(omega*constant::mathematical::twoPi*d));
+
+    const volScalarField& mu = film.mu();
+    const scalarField mup(mu.boundaryField()[patchI].patchInternalField());
+
+    const volScalarField& rho = film.rho();
+    const scalarField rhop(rho.boundaryField()[patchI].patchInternalField());
+
+    const scalarField Re(max(G, 0.0)/mup);
+
+    // TODO: currently re-evaluating the entire gTan field to return this patch
+    const scalarField gTan(film.gTan()().boundaryField()[patchI] & n);
+
+    if (max(mag(gTan)) < SMALL)
+    {
+        WarningIn
+        (
+            "void Foam::inclinedFilmNusseltHeightFvPatchScalarField::"
+            "updateCoeffs()"
+        )
+            << "Tangential gravity component is zero.  This boundary condition "
+            << "is designed to operate on patches inclined with respect to "
+            << "gravity"
+            << endl;
+    }
+
+    operator==
+    (
+        pow(3.0*sqr(mup/rhop)/(gTan + ROOTVSMALL), 0.333)*pow(Re, 0.333)
+    );
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::inclinedFilmNusseltHeightFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fixedValueFvPatchScalarField::write(os);
+    GammaMean_->writeData(os);
+    a_->writeData(os);
+    omega_->writeData(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        inclinedFilmNusseltHeightFvPatchScalarField
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.H b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.H
new file mode 100644
index 0000000000..cb9f999237
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltHeight/inclinedFilmNusseltHeightFvPatchScalarField.H
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::inclinedFilmNusseltHeightFvPatchScalarField
+
+Description
+    Film height boundary condition for inclined films that imposes a
+    sinusoidal perturbation on top of a mean flow rate, where the height is
+    calculated using the Nusselt solution.
+
+SourceFiles
+    inclinedFilmNusseltHeightFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef inclinedFilmNusseltHeightFvPatchScalarField_H
+#define inclinedFilmNusseltHeightFvPatchScalarField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+#include "DataEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+         Class inclinedFilmNusseltHeightFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class inclinedFilmNusseltHeightFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Mean mass flow rate per unit length [kg/s/m]
+        autoPtr<DataEntry<scalar> > GammaMean_;
+
+        //- Perturbation amplitude [m]
+        autoPtr<DataEntry<scalar> > a_;
+
+        //- Perturbation frequency [rad/s/m]
+        autoPtr<DataEntry<scalar> > omega_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("inclinedFilmNusseltHeight");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given 
+        // inclinedFilmNusseltHeightFvPatchScalarField onto a new patch
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const inclinedFilmNusseltHeightFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const inclinedFilmNusseltHeightFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new inclinedFilmNusseltHeightFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        inclinedFilmNusseltHeightFvPatchScalarField
+        (
+            const inclinedFilmNusseltHeightFvPatchScalarField&,
+            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 inclinedFilmNusseltHeightFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.C b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.C
new file mode 100644
index 0000000000..7f21d9a2e4
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.C
@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "inclinedFilmNusseltInletVelocityFvPatchVectorField.H"
+#include "volFields.H"
+#include "kinematicSingleLayer.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    GammaMean_(),
+    a_(),
+    omega_()
+{}
+
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const inclinedFilmNusseltInletVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
+    GammaMean_(ptf.GammaMean_().clone().ptr()),
+    a_(ptf.a_().clone().ptr()),
+    omega_(ptf.omega_().clone().ptr())
+{}
+
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchVectorField(p, iF),
+    GammaMean_(DataEntry<scalar>::New("GammaMean", dict)),
+    a_(DataEntry<scalar>::New("a", dict)),
+    omega_(DataEntry<scalar>::New("omega", dict))
+{
+    fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
+}
+
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const inclinedFilmNusseltInletVelocityFvPatchVectorField& fmfrpvf
+)
+:
+    fixedValueFvPatchVectorField(fmfrpvf),
+    GammaMean_(fmfrpvf.GammaMean_().clone().ptr()),
+    a_(fmfrpvf.a_().clone().ptr()),
+    omega_(fmfrpvf.omega_().clone().ptr())
+{}
+
+
+Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::
+inclinedFilmNusseltInletVelocityFvPatchVectorField
+(
+    const inclinedFilmNusseltInletVelocityFvPatchVectorField& fmfrpvf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchVectorField(fmfrpvf, iF),
+    GammaMean_(fmfrpvf.GammaMean_().clone().ptr()),
+    a_(fmfrpvf.a_().clone().ptr()),
+    omega_(fmfrpvf.omega_().clone().ptr())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const label patchI = patch().index();
+
+    const scalar t = db().time().timeOutputValue();
+
+    // retrieve the film region from the database
+
+    const regionModels::regionModel& region =
+        db().lookupObject<regionModels::regionModel>("surfaceFilmProperties");
+
+    const regionModels::surfaceFilmModels::kinematicSingleLayer& film =
+        dynamic_cast
+        <
+            const regionModels::surfaceFilmModels::kinematicSingleLayer&
+        >(region);
+
+    // calculate the vector tangential to the patch
+
+    const vectorField n(patch().nf());
+
+    const volVectorField& nHat = film.nHat();
+
+    const vectorField nHatp(nHat.boundaryField()[patchI].patchInternalField());
+
+    vectorField nTan(nHatp ^ n);
+    nTan /= mag(nTan) + ROOTVSMALL;
+
+    // calculate distance in patch tangential direction
+
+    const vectorField& Cf = patch().Cf();
+    scalarField d(nTan & Cf);
+
+    // calculate the wavy film height
+
+    const scalar GMean = GammaMean_->value(t);
+    const scalar a = a_->value(t);
+    const scalar omega = omega_->value(t);
+
+    const scalarField G(GMean + a*sin(omega*constant::mathematical::twoPi*d));
+
+    const volScalarField& mu = film.mu();
+    const scalarField mup(mu.boundaryField()[patchI].patchInternalField());
+
+    const volScalarField& rho = film.rho();
+    const scalarField rhop(rho.boundaryField()[patchI].patchInternalField());
+
+    const scalarField Re(max(G, 0.0)/mup);
+
+    // TODO: currently re-evaluating the entire gTan field to return this patch
+    const scalarField gTan(film.gTan()().boundaryField()[patchI] & n);
+
+    if (max(mag(gTan)) < SMALL)
+    {
+        WarningIn
+        (
+            "void Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::"
+            "updateCoeffs()"
+        )
+            << "Tangential gravity component is zero.  This boundary condition "
+            << "is designed to operate on patches inclined with respect to "
+            << "gravity"
+            << endl;
+    }
+
+    operator==(n*pow(gTan*mup/(3.0*rhop), 0.333)*pow(Re, 0.666));
+
+    fixedValueFvPatchVectorField::updateCoeffs();
+}
+
+
+void Foam::inclinedFilmNusseltInletVelocityFvPatchVectorField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchVectorField::write(os);
+    GammaMean_->writeData(os);
+    a_->writeData(os);
+    omega_->writeData(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchVectorField,
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.H b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.H
new file mode 100644
index 0000000000..211ea77561
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/derivedFvPatchFields/inclinedFilmNusseltInletVelocity/inclinedFilmNusseltInletVelocityFvPatchVectorField.H
@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::inclinedFilmNusseltInletVelocityFvPatchVectorField
+
+Description
+    Film velocity boundary condition for inclined films that imposes a
+    sinusoidal perturbation on top of a mean flow rate, where the velocity is
+    calculated using the Nusselt solution.
+
+SourceFiles
+    inclinedFilmNusseltInletVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef inclinedFilmNusseltInletVelocityFvPatchVectorField_H
+#define inclinedFilmNusseltInletVelocityFvPatchVectorField_H
+
+#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
+#include "DataEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+     Class inclinedFilmNusseltInletVelocityFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class inclinedFilmNusseltInletVelocityFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Mean mass flow rate per unit length [kg/s/m]
+        autoPtr<DataEntry<scalar> > GammaMean_;
+
+        //- Perturbation amplitude [m]
+        autoPtr<DataEntry<scalar> > a_;
+
+        //- Perturbation frequency [rad/s/m]
+        autoPtr<DataEntry<scalar> > omega_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("inclinedFilmNusseltInletVelocity");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // inclinedFilmNusseltInletVelocityFvPatchVectorField onto a new patch
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const inclinedFilmNusseltInletVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const inclinedFilmNusseltInletVelocityFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new inclinedFilmNusseltInletVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        inclinedFilmNusseltInletVelocityFvPatchVectorField
+        (
+            const inclinedFilmNusseltInletVelocityFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+            const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new inclinedFilmNusseltInletVelocityFvPatchVectorField
+                (
+                    *this, iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index 86f5b34b63..941da7e848 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -61,7 +61,7 @@ bool kinematicSingleLayer::read()
     {
         const dictionary& solution = this->solution().subDict("PISO");
         solution.lookup("momentumPredictor") >> momentumPredictor_;
-        solution.lookup("nOuterCorr") >> nOuterCorr_;
+        solution.readIfPresent("nOuterCorr", nOuterCorr_);
         solution.lookup("nCorr") >> nCorr_;
         solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
 
@@ -197,6 +197,12 @@ tmp<volScalarField> kinematicSingleLayer::pp()
 }
 
 
+void kinematicSingleLayer::correctAlpha()
+{
+    alpha_ == pos(delta_ - dimensionedScalar("SMALL", dimLength, SMALL));
+}
+
+
 void kinematicSingleLayer::updateSubmodels()
 {
     if (debug)
@@ -276,8 +282,8 @@ void kinematicSingleLayer::updateSurfaceVelocities()
     Uw_ -= nHat()*(Uw_ & nHat());
     Uw_.correctBoundaryConditions();
 
-    // TODO: apply quadratic profile to determine surface velocity
-    Us_ = U_;
+    // apply quadratic profile to surface velocity
+    Us_ = 2.0*U_;
     Us_.correctBoundaryConditions();
 }
 
@@ -439,7 +445,7 @@ kinematicSingleLayer::kinematicSingleLayer
     surfaceFilmModel(modelType, mesh, g),
 
     momentumPredictor_(solution().subDict("PISO").lookup("momentumPredictor")),
-    nOuterCorr_(readLabel(solution().subDict("PISO").lookup("nOuterCorr"))),
+    nOuterCorr_(solution().subDict("PISO").lookupOrDefault("nOuterCorr", 1)),
     nCorr_(readLabel(solution().subDict("PISO").lookup("nCorr"))),
     nNonOrthCorr_
     (
@@ -503,6 +509,19 @@ kinematicSingleLayer::kinematicSingleLayer
         ),
         regionMesh()
     ),
+    alpha_
+    (
+        IOobject
+        (
+            "alpha",
+            time().timeName(),
+            regionMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar("zero", dimless, 0.0)
+    ),
     U_
     (
         IOobject
@@ -817,6 +836,8 @@ void kinematicSingleLayer::preEvolveRegion()
         Info<< "kinematicSingleLayer::preEvolveRegion()" << endl;
     }
 
+    surfaceFilmModel::preEvolveRegion();
+
     transferPrimaryRegionThermoFields();
 
     correctThermoFields();
@@ -838,6 +859,8 @@ void kinematicSingleLayer::evolveRegion()
         Info<< "kinematicSingleLayer::evolveRegion()" << endl;
     }
 
+    correctAlpha();
+
     updateSubmodels();
 
     // Solve continuity for deltaRho_
@@ -846,7 +869,7 @@ void kinematicSingleLayer::evolveRegion()
     // Implicit pressure source coefficient - constant
     tmp<volScalarField> tpp(this->pp());
 
-    for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
+    for (int oCorr=1; oCorr<=nOuterCorr_; oCorr++)
     {
         // Explicit pressure source contribution - varies with delta_
         tmp<volScalarField> tpu(this->pu());
@@ -1032,7 +1055,9 @@ void kinematicSingleLayer::info() const
         << 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;
+        << max(delta_).value() << nl
+        << indent << "coverage           = "
+        << gSum(alpha_.internalField()*magSf())/gSum(magSf()) <<  nl;
 
     injection_.info(Info);
 }
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H
index b084f885d4..079999d0cd 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,6 +113,9 @@ protected:
             //- Film thickness / [m]
             volScalarField delta_;
 
+            //- Film coverage indicator, 1 = covered, 0 = uncovered / []
+            volScalarField alpha_;
+
             //- Velocity - mean / [m/s]
             volVectorField U_;
 
@@ -221,12 +224,15 @@ protected:
         //- Transfer source fields from the primary region to the film region
         virtual void transferPrimaryRegionSourceFields();
 
-        // Explicit pressure source contribution
+        //- Explicit pressure source contribution
         virtual tmp<volScalarField> pu();
 
-        // Implicit pressure source coefficient
+        //- Implicit pressure source coefficient
         virtual tmp<volScalarField> pp();
 
+        //- Correct film coverage field
+        virtual void correctAlpha();
+
         //- Update the film sub-models
         virtual void updateSubmodels();
 
@@ -323,6 +329,9 @@ public:
             //- Return const access to the film thickness / [m]
             inline const volScalarField& delta() const;
 
+            //- Return the film coverage, 1 = covered, 0 = uncovered / []
+            inline const volScalarField& alpha() const;
+
             //- Return the film velocity [m/s]
             virtual const volVectorField& U() const;
 
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
index e717acc430..6a77b614bc 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,6 +79,12 @@ inline const volScalarField& kinematicSingleLayer::delta() const
 }
 
 
+inline const volScalarField& kinematicSingleLayer::alpha() const
+{
+    return alpha_;
+}
+
+
 inline volVectorField& kinematicSingleLayer::USpPrimary()
 {
     return USpPrimary_;
@@ -173,7 +179,8 @@ inline tmp<volScalarField> kinematicSingleLayer::netMass() const
 {
     dimensionedScalar d0("SMALL", dimLength, ROOTVSMALL);
     return
-        fvc::surfaceSum(phi_/(fvc::interpolate(delta_) + d0))*time().deltaT()
+        fvc::surfaceSum(pos(phi_)*phi_/(fvc::interpolate(delta_) + d0))
+       *time().deltaT()
       + rho_*delta_*magSf();
 }
 
diff --git a/src/regionModels/surfaceFilmModels/noFilm/noFilm.C b/src/regionModels/surfaceFilmModels/noFilm/noFilm.C
index 1c267307ec..82f3d0afc9 100644
--- a/src/regionModels/surfaceFilmModels/noFilm/noFilm.C
+++ b/src/regionModels/surfaceFilmModels/noFilm/noFilm.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,10 +102,10 @@ const volScalarField& noFilm::delta() const
 }
 
 
-const volScalarField& noFilm::sigma() const
+const volScalarField& noFilm::alpha() const
 {
-    FatalErrorIn("const volScalarField& noFilm::sigma() const")
-        << "sigma field not available for " << type() << abort(FatalError);
+    FatalErrorIn("const volScalarField& noFilm::alpha() const")
+        << "alpha field not available for " << type() << abort(FatalError);
 
     return volScalarField::null();
 }
@@ -192,6 +192,15 @@ const volScalarField& noFilm::kappa() const
 }
 
 
+const volScalarField& noFilm::sigma() const
+{
+    FatalErrorIn("const volScalarField& noFilm::sigma() const")
+        << "sigma field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
 tmp<volScalarField> noFilm::primaryMassTrans() const
 {
     return tmp<volScalarField>
diff --git a/src/regionModels/surfaceFilmModels/noFilm/noFilm.H b/src/regionModels/surfaceFilmModels/noFilm/noFilm.H
index 76501a4fd8..17915f7788 100644
--- a/src/regionModels/surfaceFilmModels/noFilm/noFilm.H
+++ b/src/regionModels/surfaceFilmModels/noFilm/noFilm.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,8 +115,8 @@ public:
             //- Return the film thickness [m]
             virtual const volScalarField& delta() const;
 
-            //- Return const access to the surface tension / [m/s2]
-            inline const volScalarField& sigma() const;
+            //- Return the film coverage, 1 = covered, 0 = uncovered / []
+            virtual const volScalarField& alpha() const;
 
             //- Return the film velocity [m/s]
             virtual const volVectorField& U() const;
@@ -145,6 +145,9 @@ public:
             //- Return the film thermal conductivity [W/m/K]
             virtual const volScalarField& kappa() const;
 
+            //- Return const access to the surface tension / [m/s2]
+            inline const volScalarField& sigma() const;
+
 
             // Transfer fields - to the primary region
 
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
index 5f31601474..2ceff1429f 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,7 +52,6 @@ contactAngleForce::contactAngleForce
 )
 :
     force(typeName, owner, dict),
-    deltaWet_(readScalar(coeffs_.lookup("deltaWet"))),
     Ccf_(readScalar(coeffs_.lookup("Ccf"))),
     rndGen_(label(0), -1),
     distribution_
@@ -82,7 +81,7 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         (
             IOobject
             (
-                "contactForce",
+                typeName + ".contactForce",
                 owner_.time().timeName(),
                 owner_.regionMesh(),
                 IOobject::NO_READ,
@@ -100,17 +99,11 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
 
     const scalarField& magSf = owner_.magSf();
 
-    const volScalarField& delta = owner_.delta();
+    const volScalarField& alpha = owner_.alpha();
     const volScalarField& sigma = owner_.sigma();
 
-    volScalarField alpha
-    (
-        "alpha",
-        pos(delta - dimensionedScalar("deltaWet", dimLength, deltaWet_))
-    );
     volVectorField gradAlpha(fvc::grad(alpha));
 
-
     scalarField nHits(owner_.regionMesh().nCells(), 0.0);
 
     forAll(nbr, faceI)
@@ -119,19 +112,17 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         const label cellN = nbr[faceI];
 
         label cellI = -1;
-        if ((delta[cellO] > deltaWet_) && (delta[cellN] < deltaWet_))
+        if ((alpha[cellO] > 0.5) && (alpha[cellN] < 0.5))
         {
             cellI = cellO;
         }
-        else if ((delta[cellO] < deltaWet_) && (delta[cellN] > deltaWet_))
+        else if ((alpha[cellO] < 0.5) && (alpha[cellN] > 0.5))
         {
             cellI = cellN;
         }
 
         if (cellI != -1)
         {
-//            const scalar dx = Foam::sqrt(magSf[cellI]);
-            // bit of a cheat, but ok for regular meshes
             const scalar dx = owner_.regionMesh().deltaCoeffs()[faceI];
             const vector n =
                 gradAlpha[cellI]/(mag(gradAlpha[cellI]) + ROOTVSMALL);
@@ -141,17 +132,17 @@ tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
         }
     }
 
-    forAll(delta.boundaryField(), patchI)
+    forAll(alpha.boundaryField(), patchI)
     {
-        const fvPatchField<scalar>& df = delta.boundaryField()[patchI];
-        const scalarField& dx = df.patch().deltaCoeffs();
-        const labelUList& faceCells = df.patch().faceCells();
+        const fvPatchField<scalar>& alphaf = alpha.boundaryField()[patchI];
+        const scalarField& dx = alphaf.patch().deltaCoeffs();
+        const labelUList& faceCells = alphaf.patch().faceCells();
 
-        forAll(df, faceI)
+        forAll(alphaf, faceI)
         {
             label cellO = faceCells[faceI];
 
-            if ((delta[cellO] > deltaWet_) && (df[faceI] < deltaWet_))
+            if ((alpha[cellO] > 0.5) && (alphaf[faceI] < 0.5))
             {
                 const vector n =
                     gradAlpha[cellO]/(mag(gradAlpha[cellO]) + ROOTVSMALL);
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H
index 2ed08e02ee..8a6dc7a8b3 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/contactAngleForce/contactAngleForce.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,9 +60,6 @@ private:
 
     // Private Data
 
-        //- Threshold film thickness beyon which the film is 'wet'
-        scalar deltaWet_;
-
         //- Coefficient applied to the contact angle force
         scalar Ccf_;
 
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/surfaceShearForce/surfaceShearForce.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/surfaceShearForce/surfaceShearForce.C
index c5419d2909..a2c910c412 100644
--- a/src/regionModels/surfaceFilmModels/submodels/kinematic/force/surfaceShearForce/surfaceShearForce.C
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/force/surfaceShearForce/surfaceShearForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,11 +80,11 @@ tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
     tmp<volScalarField> tCs;
 
     typedef compressible::turbulenceModel turbModel;
-    if (film.primaryMesh().foundObject<turbModel>("turbulenceProperties"))
+    if (film.primaryMesh().foundObject<turbModel>("turbulenceModel"))
     {
         // local reference to turbulence model
         const turbModel& turb =
-            film.primaryMesh().lookupObject<turbModel>("turbulenceProperties");
+            film.primaryMesh().lookupObject<turbModel>("turbulenceModel");
 
         // calculate and store the stress on the primary region
         const volSymmTensorField primaryReff(turb.devRhoReff());
@@ -117,13 +117,23 @@ tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
         Reff.correctBoundaryConditions();
 
         dimensionedScalar U0("SMALL", U.dimensions(), SMALL);
-        tCs = Cf_*mag(-film.nHat() & Reff)/(mag(Up - U) + U0);
+        volVectorField UHat("UHat", (Up - U)/(mag(Up - U) + U0));
+
+        // shear stress tangential to the film
+        volVectorField tauTan
+        (
+            "tauTan",
+            UHat & (Reff + film.nHat()*(-film.nHat() & Reff))
+        );
+
+        // note: Cf_ 'should' be 1 in this case
+        tCs = Cf_*mag(tauTan)/(mag(Up - U) + U0);
     }
     else
     {
         // laminar case - employ simple coeff-based model
-        const volScalarField& rho = film.rho();
-        tCs = Cf_*rho*mag(Up - U);
+        const volScalarField& rhop = film.rhoPrimary();
+        tCs = Cf_*rhop*mag(Up - U);
     }
 
     dimensionedScalar d0("SMALL", delta.dimensions(), SMALL);
@@ -131,7 +141,7 @@ tmp<fvVectorMatrix> surfaceShearForce::correct(volVectorField& U)
     // linear coeffs to apply to velocity
     const volScalarField& Cs = tCs();
     volScalarField Cw("Cw", mu/(0.3333*(delta + d0)));
-    Cw.min(1.0e+06);
+    Cw.min(5000.0);
 
     return
     (
diff --git a/src/regionModels/surfaceFilmModels/submodels/subModelBase.C b/src/regionModels/surfaceFilmModels/submodels/subModelBase.C
index 469b62e855..b2283172b3 100644
--- a/src/regionModels/surfaceFilmModels/submodels/subModelBase.C
+++ b/src/regionModels/surfaceFilmModels/submodels/subModelBase.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,7 +51,7 @@ subModelBase::subModelBase
 )
 :
     owner_(owner),
-    coeffs_(dict.subDict(type + "Coeffs"))
+    coeffs_(dict.subOrEmptyDict(type + "Coeffs"))
 {}
 
 
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C
new file mode 100644
index 0000000000..ac0358f141
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.C
@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "constantRadiation.H"
+#include "volFields.H"
+#include "zeroGradientFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(constantRadiation, 0);
+
+addToRunTimeSelectionTable
+(
+    filmRadiationModel,
+    constantRadiation,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+constantRadiation::constantRadiation
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    filmRadiationModel(typeName, owner, dict),
+    QrConst_
+    (
+        IOobject
+        (
+            typeName + "::QrConst",
+            owner.time().timeName(),
+            owner.regionMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        owner.regionMesh()
+    ),
+    mask_
+    (
+        IOobject
+        (
+            typeName + "::mask",
+            owner.time().timeName(),
+            owner.regionMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        owner.regionMesh(),
+        dimensionedScalar("one", dimless, 1.0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    timeStart_(readScalar(coeffs_.lookup("timeStart"))),
+    duration_(readScalar(coeffs_.lookup("duration")))
+{
+    mask_ = pos(mask_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+constantRadiation::~constantRadiation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void constantRadiation::correct()
+{}
+
+
+tmp<volScalarField> constantRadiation::Shs()
+{
+    tmp<volScalarField> tShs
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                typeName + "::Shs",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    const scalar time = owner().time().value();
+
+    if ((time >= timeStart_) && (time <= timeStart_ + duration_))
+    {
+        scalarField& Shs = tShs();
+        const scalarField& Qr = QrConst_.internalField();
+        const scalarField& alpha = owner_.alpha().internalField();
+
+        Shs = mask_*Qr*alpha;
+    }
+
+    return tShs;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.H b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.H
new file mode 100644
index 0000000000..06e6fe676b
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/constantRadiation/constantRadiation.H
@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::constantRadiation
+
+Description
+    Film constant radiation model.  The constant radiative flux is specified
+    by the user, and operated over a time interval defined by a start time and
+    duration.  In addition, a mask can be applied to shield the film from the
+    radiation.
+
+SourceFiles
+    constantRadiation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantRadiation_H
+#define constantRadiation_H
+
+#include "filmRadiationModel.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class constantRadiation Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantRadiation
+:
+    public filmRadiationModel
+{
+private:
+
+    // Private data
+
+        //- Constant radiative flux [kg/s3]
+        volScalarField QrConst_;
+
+        //- Radiation mask
+        volScalarField mask_;
+
+        //- Time start [s]
+        const scalar timeStart_;
+
+        //- Duration [s]
+        const scalar duration_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        constantRadiation(const constantRadiation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const constantRadiation&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantRadiation");
+
+
+    // Constructors
+
+        //- Construct from surface film model and dictionary
+        constantRadiation
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~constantRadiation();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct();
+
+            //- Return the radiation sensible enthalpy source
+            //  Also updates QrNet
+            virtual tmp<volScalarField> Shs();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
new file mode 100644
index 0000000000..2fdccf63f7
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.C
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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 "primaryRadiation.H"
+#include "volFields.H"
+#include "zeroGradientFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(primaryRadiation, 0);
+
+addToRunTimeSelectionTable
+(
+    filmRadiationModel,
+    primaryRadiation,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+primaryRadiation::primaryRadiation
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    filmRadiationModel(typeName, owner, dict),
+    QinPrimary_
+    (
+        IOobject
+        (
+            "Qin", // same name as Qin on primary region to enable mapping
+            owner.time().timeName(),
+            owner.regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        owner.regionMesh(),
+        dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
+        owner.mappedPushedFieldPatchTypes<scalar>()
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+primaryRadiation::~primaryRadiation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void primaryRadiation::correct()
+{
+    // Transfer Qin from primary region
+    QinPrimary_.correctBoundaryConditions();
+}
+
+
+tmp<volScalarField> primaryRadiation::Shs()
+{
+    tmp<volScalarField> tShs
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                typeName + "::Shs",
+                owner().time().timeName(),
+                owner().regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner().regionMesh(),
+            dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    scalarField& Shs = tShs();
+    const scalarField& QinP = QinPrimary_.internalField();
+    const scalarField& alpha = owner_.alpha().internalField();
+
+    Shs = QinP*alpha;
+
+    return tShs;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H
new file mode 100644
index 0000000000..a8904cb1e2
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/primaryRadiation/primaryRadiation.H
@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::primaryRadiation
+
+Description
+    Radiation model whereby the radiative heat flux is mapped from the primary
+    region
+
+SourceFiles
+    primaryRadiation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef primaryRadiation_H
+#define primaryRadiation_H
+
+#include "filmRadiationModel.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class primaryRadiation Declaration
+\*---------------------------------------------------------------------------*/
+
+class primaryRadiation
+:
+    public filmRadiationModel
+{
+private:
+
+    // Private data
+
+        //- Incident radiative flux mapped from  the primary region / [kg/s3]
+        volScalarField QinPrimary_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        primaryRadiation(const primaryRadiation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const primaryRadiation&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("primaryRadiation");
+
+
+    // Constructors
+
+        //- Construct from surface film model and dictionary
+        primaryRadiation
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~primaryRadiation();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct();
+
+            //- Return the radiation sensible enthalpy source
+            //  Also updates QrNet
+            virtual tmp<volScalarField> Shs();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
index be5460757f..1ff3274716 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C
@@ -85,8 +85,6 @@ standardRadiation::standardRadiation
         dimensionedScalar("zero", dimMass/pow3(dimTime), 0.0),
         zeroGradientFvPatchScalarField::typeName
     ),
-    delta_(owner.delta()),
-    deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
     beta_(readScalar(coeffs_.lookup("beta"))),
     kappaBar_(readScalar(coeffs_.lookup("kappaBar")))
 {}
@@ -129,9 +127,10 @@ tmp<volScalarField> standardRadiation::Shs()
 
     scalarField& Shs = tShs();
     const scalarField& QinP = QinPrimary_.internalField();
-    const scalarField& delta = delta_.internalField();
+    const scalarField& delta = owner_.delta().internalField();
+    const scalarField& alpha = owner_.alpha().internalField();
 
-    Shs = beta_*(QinP*pos(delta - deltaMin_))*(1.0 - exp(-kappaBar_*delta));
+    Shs = beta_*QinP*alpha*(1.0 - exp(-kappaBar_*delta));
 
     // Update net Qr on local region
     QrNet_.internalField() = QinP - Shs;
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H
index 9b96fb60fe..799bb6c649 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.H
@@ -65,15 +65,9 @@ private:
         //- Remaining radiative flux after removing local contribution
         volScalarField QrNet_;
 
-        //- Reference to the film thickness field / [m]
-        const volScalarField& delta_;
-
 
         // Model coefficients
 
-            //- Minimum thickness to apply radiation model
-            scalar deltaMin_;
-
             //- Beta coefficient
             scalar beta_;
 
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
index 552bd0dda4..4a08755ee3 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,7 +77,6 @@ standardPhaseChange::standardPhaseChange
 )
 :
     phaseChangeModel(typeName, owner, dict),
-    Tb_(readScalar(coeffs_.lookup("Tb"))),
     deltaMin_(readScalar(coeffs_.lookup("deltaMin"))),
     L_(readScalar(coeffs_.lookup("L"))),
     TbFactor_(coeffs_.lookupOrDefault<scalar>("TbFactor", 1.1))
@@ -113,14 +112,10 @@ void standardPhaseChange::correctModel
     const scalarField& YInf = film.YPrimary()[vapId];
     const scalarField& pInf = film.pPrimary();
     const scalarField& T = film.T();
-    const scalarField& Tw = film.Tw();
     const scalarField& rho = film.rho();
-    const scalarField& TInf = film.TPrimary();
     const scalarField& rhoInf = film.rhoPrimary();
     const scalarField& muInf = film.muPrimary();
     const scalarField& magSf = film.magSf();
-    const scalarField hInf(film.htcs().h());
-    const scalarField hFilm(film.htcw().h());
     const vectorField dU(film.UPrimary() - film.Us());
     const scalarField limMass
     (
@@ -134,8 +129,11 @@ void standardPhaseChange::correctModel
             // cell pressure [Pa]
             const scalar pc = pInf[cellI];
 
+            // calculate the boiling temperature
+            const scalar Tb = liq.pvInvert(pc);
+
             // local temperature - impose lower limit of 200 K for stability
-            const scalar Tloc = min(TbFactor_*Tb_, max(200.0, T[cellI]));
+            const scalar Tloc = min(TbFactor_*Tb, max(200.0, T[cellI]));
 
             // saturation pressure [Pa]
             const scalar pSat = liq.pv(pc, Tloc);
@@ -147,15 +145,10 @@ void standardPhaseChange::correctModel
             if (pSat >= 0.95*pc)
             {
                 // boiling
-                const scalar qDotInf = hInf[cellI]*(TInf[cellI] - T[cellI]);
-                const scalar qDotFilm = hFilm[cellI]*(T[cellI] - Tw[cellI]);
-
                 const scalar Cp = liq.Cp(pc, Tloc);
-                const scalar Tcorr = max(0.0, T[cellI] - Tb_);
+                const scalar Tcorr = max(0.0, T[cellI] - Tb);
                 const scalar qCorr = limMass[cellI]*Cp*(Tcorr);
-                dMass[cellI] =
-                    dt*magSf[cellI]/hVap*(qDotInf + qDotFilm)
-                  + qCorr/hVap;
+                dMass[cellI] = qCorr/hVap;
             }
             else
             {
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
index 1fff3384ff..ee1b21d6d0 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/standardPhaseChange/standardPhaseChange.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,9 +69,6 @@ protected:
 
     // Protected data
 
-        //- Boiling temperature / [K]
-        const scalar Tb_;
-
         //- Minimum film height for model to be active
         const scalar deltaMin_;
 
diff --git a/src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel.H b/src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel.H
index 32c71b3fed..2aa4423196 100644
--- a/src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel.H
+++ b/src/regionModels/surfaceFilmModels/surfaceFilmModel/surfaceFilmModel.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -186,6 +186,9 @@ public:
             //- Return the film thickness [m]
             virtual const volScalarField& delta() const = 0;
 
+            //- Return the film coverage, 1 = covered, 0 = uncovered / []
+            virtual const volScalarField& alpha() const = 0;
+
             //- Return the film velocity [m/s]
             virtual const volVectorField& U() const = 0;
 
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
index 1cb392eb6a..21953ae94c 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
@@ -227,6 +227,33 @@ void thermoSingleLayer::transferPrimaryRegionSourceFields()
 }
 
 
+void thermoSingleLayer::correctAlpha()
+{
+    if (hydrophilic_)
+    {
+        const scalar hydrophilicDry = hydrophilicDryScale_*deltaWet_;
+        const scalar hydrophilicWet = hydrophilicWetScale_*deltaWet_;
+
+        forAll(alpha_, i)
+        {
+            if ((alpha_[i] < 0.5) && (delta_[i] > hydrophilicDry))
+            {
+                alpha_[i] = 1.0;
+            }
+            else if ((alpha_[i] > 0.5) && (delta_[i] < hydrophilicWet))
+            {
+                alpha_[i] = 0.0;
+            }
+        }
+    }
+    else
+    {
+        alpha_ ==
+            pos(delta_ - dimensionedScalar("deltaWet", dimLength, deltaWet_));
+    }
+}
+
+
 void thermoSingleLayer::updateSubmodels()
 {
     if (debug)
@@ -290,7 +317,8 @@ void thermoSingleLayer::solveEnergy()
       - hsSp_
       + q(hs_)
       + radiation_->Shs()
-      - fvm::SuSp(rhoSp_, hs_)
+//      - fvm::SuSp(rhoSp_, hs_)
+      - rhoSp_*hs_
     );
 
     correctThermoFields();
@@ -426,6 +454,11 @@ thermoSingleLayer::thermoSingleLayer
         zeroGradientFvPatchScalarField::typeName
     ),
 
+    deltaWet_(readScalar(coeffs_.lookup("deltaWet"))),
+    hydrophilic_(readBool(coeffs_.lookup("hydrophilic"))),
+    hydrophilicDryScale_(0.0),
+    hydrophilicWetScale_(0.0),
+
     hsSp_
     (
         IOobject
@@ -481,8 +514,20 @@ thermoSingleLayer::thermoSingleLayer
         heatTransferModel::New(*this, coeffs().subDict("lowerSurfaceModels"))
     ),
     phaseChange_(phaseChangeModel::New(*this, coeffs())),
-    radiation_(filmRadiationModel::New(*this, coeffs()))
+    radiation_(filmRadiationModel::New(*this, coeffs())),
+    Tmin_(-VGREAT),
+    Tmax_(VGREAT)
 {
+    if (coeffs().readIfPresent("Tmin", Tmin_))
+    {
+        Info<< "    limiting minimum temperature to " << Tmin_ << endl;
+    }
+
+    if (coeffs().readIfPresent("Tmax", Tmax_))
+    {
+        Info<< "    limiting maximum temperature to " << Tmax_ << endl;
+    }
+
     if (thermo_.hasMultiComponentCarrier())
     {
         YPrimary_.setSize(thermo_.carrier().species().size());
@@ -510,6 +555,12 @@ thermoSingleLayer::thermoSingleLayer
         }
     }
 
+    if (hydrophilic_)
+    {
+        coeffs_.lookup("hydrophilicDryScale") >> hydrophilicDryScale_;
+        coeffs_.lookup("hydrophilicWetScale") >> hydrophilicWetScale_;
+    }
+
     if (readFields)
     {
         transferPrimaryRegionThermoFields();
@@ -584,12 +635,14 @@ void thermoSingleLayer::evolveRegion()
         Info<< "thermoSingleLayer::evolveRegion()" << endl;
     }
 
+    correctAlpha();
+
     updateSubmodels();
 
     // Solve continuity for deltaRho_
     solveContinuity();
 
-    for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
+    for (int oCorr=1; oCorr<=nOuterCorr_; oCorr++)
     {
         // Explicit pressure source contribution
         tmp<volScalarField> tpu(this->pu());
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H
index 6814e25a3e..696728e7f9 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -123,6 +123,22 @@ protected:
                 volScalarField primaryEnergyPCTrans_;
 
 
+        //- Threshold film thickness beyond which the film is considered 'wet'
+        scalar deltaWet_;
+
+
+        // Hyprophilic/phobic properties
+
+            //- Activation flag
+            bool hydrophilic_;
+
+            //- Length scale applied to deltaWet_ for dry faces, typically 0.5
+            scalar hydrophilicDryScale_;
+
+            //- Length scale applied to deltaWet_ for wet faces, typically 0.001
+            scalar hydrophilicWetScale_;
+
+
         // Source term fields
 
             // Film region - registered to the film region mesh
@@ -166,6 +182,14 @@ protected:
             autoPtr<filmRadiationModel> radiation_;
 
 
+        // Limits
+
+            //- Minimum temperature limit (optional)
+            scalar Tmin_;
+
+            //- Maximum temperature limit (optional)
+            scalar Tmax_;
+
 
     // Protected member functions
 
@@ -190,6 +214,9 @@ protected:
         //- Transfer source fields from the primary region to the film region
         virtual void transferPrimaryRegionSourceFields();
 
+        //- Correct film coverage field
+        virtual void correctAlpha();
+
         //- Update the film sub-models
         virtual void updateSubmodels();
 
@@ -342,6 +369,15 @@ public:
             inline const filmRadiationModel& radiation() const;
 
 
+        // Derived fields (calculated on-the-fly)
+
+            //- Return the convective heat energy from film to wall
+            inline tmp<scalarField> Qconvw(const label patchI) const;
+
+            //- Return the convective heat energy from primary region to film
+            inline tmp<scalarField> Qconvp(const label patchI) const;
+
+
         // Evolution
 
             //- Pre-evolve film hook
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
index 73dc08496b..659a315f49 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "thermoSingleLayer.H"
+#include "heatTransferModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -88,7 +89,7 @@ inline tmp<volScalarField> thermoSingleLayer::T
     const volScalarField& hs
 ) const
 {
-    return tmp<volScalarField>
+    tmp<volScalarField> tT
     (
         new volScalarField
         (
@@ -104,6 +105,11 @@ inline tmp<volScalarField> thermoSingleLayer::T
             zeroGradientFvPatchScalarField::typeName
         )
     );
+
+    tT().min(Tmax_);
+    tT().max(Tmin_);
+
+    return tT;
 }
 
 
@@ -155,6 +161,26 @@ inline const filmRadiationModel& thermoSingleLayer::radiation() const
 }
 
 
+inline tmp<scalarField> thermoSingleLayer::Qconvw(const label patchI) const
+{
+    const scalarField htc(htcw_->h()().boundaryField()[patchI]);
+    const scalarField& Tp = T_.boundaryField()[patchI];
+    const scalarField& Twp = Tw_.boundaryField()[patchI];
+
+    return htc*(Tp - Twp);
+}
+
+
+inline tmp<scalarField> thermoSingleLayer::Qconvp(const label patchI) const
+{
+    const scalarField htc(htcs_->h()().boundaryField()[patchI]);
+    const scalarField& Tp = T_.boundaryField()[patchI];
+    const scalarField& Tpp = TPrimary_.boundaryField()[patchI];
+
+    return htc*(Tp - Tpp);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace surfaceFilmModels
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C
index 9e5578b886..9851f44b6f 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJump/energyJumpFvPatchScalarField.C
@@ -68,6 +68,10 @@ Foam::energyJumpFvPatchScalarField::energyJumpFvPatchScalarField
             scalarField("value", dict, p.size())
         );
     }
+    else
+    {
+        evaluate(Pstream::blocking);
+    }
 }
 
 
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C
index 7f6a09c181..cc8675dad9 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/energyJump/energyJumpAMI/energyJumpAMIFvPatchScalarField.C
@@ -68,6 +68,10 @@ Foam::energyJumpAMIFvPatchScalarField::energyJumpAMIFvPatchScalarField
             scalarField("value", dict, p.size())
         );
     }
+    else
+    {
+        evaluate(Pstream::blocking);
+    }
 }
 
 
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.C b/src/thermophysicalModels/basic/heThermo/heThermo.C
index bba8da1514..8efaf0f5d5 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.C
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.C
@@ -35,6 +35,39 @@ License
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+template<class BasicThermo, class MixtureType>
+Foam::wordList Foam::heThermo<BasicThermo, MixtureType>::heBoundaryBaseTypes()
+{
+    const volScalarField::GeometricBoundaryField& tbf =
+        this->T_.boundaryField();
+
+    wordList hbt(tbf.size(), word::null);
+
+    forAll(tbf, patchi)
+    {
+        if (isA<fixedJumpFvPatchScalarField>(tbf[patchi]))
+        {
+            const fixedJumpFvPatchScalarField& pf =
+                dynamic_cast<const fixedJumpFvPatchScalarField&>(tbf[patchi]);
+
+            hbt[patchi] = pf.interfaceFieldType();
+        }
+        else if (isA<fixedJumpAMIFvPatchScalarField>(tbf[patchi]))
+        {
+            const fixedJumpAMIFvPatchScalarField& pf =
+                dynamic_cast<const fixedJumpAMIFvPatchScalarField&>
+                (
+                    tbf[patchi]
+                );
+
+            hbt[patchi] = pf.interfaceFieldType();
+        }
+    }
+
+    return hbt;
+}
+
+
 template<class BasicThermo, class MixtureType>
 Foam::wordList Foam::heThermo<BasicThermo, MixtureType>::heBoundaryTypes()
 {
@@ -149,7 +182,8 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         ),
         mesh,
         dimEnergy/dimMass,
-        this->heBoundaryTypes()
+        this->heBoundaryTypes(),
+        this->heBoundaryBaseTypes()
     )
 {
     init();
@@ -179,7 +213,8 @@ Foam::heThermo<BasicThermo, MixtureType>::heThermo
         ),
         mesh,
         dimEnergy/dimMass,
-        this->heBoundaryTypes()
+        this->heBoundaryTypes(),
+        this->heBoundaryBaseTypes()
     )
 {
     init();
diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.H b/src/thermophysicalModels/basic/heThermo/heThermo.H
index 2c894eee43..ef9919a99c 100644
--- a/src/thermophysicalModels/basic/heThermo/heThermo.H
+++ b/src/thermophysicalModels/basic/heThermo/heThermo.H
@@ -68,6 +68,10 @@ protected:
             //  by interrogating the temperature field boundary types
             wordList heBoundaryTypes();
 
+            //- Return the enthalpy/internal energy field boundary base types
+            //  by interrogating the temperature field boundary types
+            wordList heBoundaryBaseTypes();
+
             //- Correct the enthalpy/internal energy field boundaries
             void heBoundaryCorrection(volScalarField& he);
 
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/MRFProperties b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/MRFProperties
index 3782a377cd..c4c20f18d0 100644
--- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/MRFProperties
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/MRFProperties
@@ -19,6 +19,8 @@ zone1
 {
     cellZone    rotor;
 
+    active      on;
+
     // Fixed patches (by default they 'move' with the MRF zone)
     nonRotatingPatches ();
 
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicitFixedCoeff/constant/sourcesProperties b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicitFixedCoeff/constant/sourcesProperties
index 4386bc5250..404024507d 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicitFixedCoeff/constant/sourcesProperties
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctExplicitFixedCoeff/constant/sourcesProperties
@@ -26,6 +26,7 @@ source1
 
     fixedTemperatureSourceCoeffs
     {
+        mode            constant;
         temperature     350;
     }
 }
diff --git a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/sourcesProperties b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/sourcesProperties
index 4386bc5250..404024507d 100644
--- a/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/sourcesProperties
+++ b/tutorials/compressible/rhoPorousMRFSimpleFoam/angledDuctImplicit/constant/sourcesProperties
@@ -26,6 +26,7 @@ source1
 
     fixedTemperatureSourceCoeffs
     {
+        mode            constant;
         temperature     350;
     }
 }
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
index 715f0ebdd9..f3eac3d8de 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/sourcesProperties
@@ -26,6 +26,7 @@ source1
 
     fixedTemperatureSourceCoeffs
     {
+        mode            constant;
         temperature     2000;
     }
 }
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/surfaceFilmProperties b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/surfaceFilmProperties
index 98c728fdec..e327a82959 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/surfaceFilmProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/surfaceFilmProperties
@@ -26,6 +26,9 @@ thermoSingleLayerCoeffs
     thermoModel singleComponent;
     liquid      H2O;
 
+    deltaWet    1e-4;
+    hydrophilic no;
+
     forces
     (
         surfaceShear
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/wallFilmRegion/fvSolution b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/wallFilmRegion/fvSolution
index 381a5c4f27..eedbd097d3 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/wallFilmRegion/fvSolution
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/system/wallFilmRegion/fvSolution
@@ -37,8 +37,7 @@ solvers
 PISO
 {
     momentumPredictor true;
-    nOuterCorr      1;
-    nCorr           2;
+    nCorr           1;
     nNonOrthCorr    0;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties
index 6984dcbb1e..96e33c1d68 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties
@@ -28,6 +28,9 @@ thermoSingleLayerCoeffs
 
     liquid      H2O;
 
+    deltaWet    1e-4;
+    hydrophilic no;
+
     radiationModel  none;
 
     upperSurfaceModels
@@ -65,7 +68,6 @@ thermoSingleLayerCoeffs
 
     contactAngleCoeffs
     {
-        deltaWet        1e-4;
         Ccf             0.085;
         contactAngleDistribution
         {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/wallFilmRegion.org/fvSolution b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/wallFilmRegion.org/fvSolution
index cfae8b06ce..9385e6ec40 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/wallFilmRegion.org/fvSolution
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/wallFilmRegion.org/fvSolution
@@ -39,8 +39,7 @@ solvers
 PISO
 {
     momentumPredictor true;
-    nOuterCorr      1;
-    nCorr           2;
+    nCorr           1;
     nNonOrthCorr    0;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/surfaceFilmProperties b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/surfaceFilmProperties
index d9a68e059d..8440f3f479 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/surfaceFilmProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/constant/surfaceFilmProperties
@@ -43,7 +43,6 @@ kinematicSingleLayerCoeffs
 
     contactAngleCoeffs
     {
-        deltaWet        1e-4;
         Ccf             1;
         contactAngleDistribution
         {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/wallFilmRegion/fvSolution b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/wallFilmRegion/fvSolution
index 645d42a456..86cdc388bf 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/wallFilmRegion/fvSolution
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/system/wallFilmRegion/fvSolution
@@ -37,8 +37,7 @@ solvers
 PISO
 {
     momentumPredictor  true;
-    nOuterCorr      1;
-    nCorr           2;
+    nCorr           1;
     nNonOrthCorr    0;
 }
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/surfaceFilmProperties b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/surfaceFilmProperties
index 30a2dd79f3..8335fadfc2 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/surfaceFilmProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/surfaceFilmProperties
@@ -27,6 +27,9 @@ thermoSingleLayerCoeffs
 
     liquid      H2O;
 
+    deltaWet    1e-4;
+    hydrophilic no;
+
     forces
     (
         surfaceShear
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/system/wallFilmRegion/fvSolution b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/system/wallFilmRegion/fvSolution
index 381a5c4f27..eedbd097d3 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/system/wallFilmRegion/fvSolution
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/system/wallFilmRegion/fvSolution
@@ -37,8 +37,7 @@ solvers
 PISO
 {
     momentumPredictor true;
-    nOuterCorr      1;
-    nCorr           2;
+    nCorr           1;
     nNonOrthCorr    0;
 }