diff --git a/src/lagrangian/coalCombustion/Make/options b/src/lagrangian/coalCombustion/Make/options
index 25ee31bf7dfc939af4f6adaeb43280e93cc02261..594770817691aa072cb71de3748c523f137b9385 100644
--- a/src/lagrangian/coalCombustion/Make/options
+++ b/src/lagrangian/coalCombustion/Make/options
@@ -18,7 +18,10 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -41,4 +44,6 @@ LIB_LIBS = \
     -lregionModels \
     -lsurfaceFilmModels \
     -ldynamicMesh \
-    -ldynamicFvMesh
+    -ldynamicFvMesh \
+    -lregionFaModels \
+    -lfiniteArea
diff --git a/src/lagrangian/intermediate/Make/options b/src/lagrangian/intermediate/Make/options
index d972d5b0345590520cdb943aa98de8f2c2e29f10..74eb95ca4c0e8adf79c358849f15e9692f3d53f0 100644
--- a/src/lagrangian/intermediate/Make/options
+++ b/src/lagrangian/intermediate/Make/options
@@ -16,7 +16,10 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -36,4 +39,6 @@ LIB_LIBS = \
     -lregionModels \
     -lsurfaceFilmModels \
     -ldynamicMesh \
-    -ldynamicFvMesh
+    -ldynamicFvMesh \
+    -lregionFaModels \
+    -lfiniteArea
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
index 5ef693076de29917888cc2252304077a92ddd5bd..4696ac7021fe051c7615157b9d06926291d39996 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
@@ -27,8 +27,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SurfaceFilmModel.H"
-#include "surfaceFilmRegionModel.H"
 #include "mathematicalConstants.H"
+#include "surfaceFilmRegionModel.H"
+#include "liquidFilmBase.H"
 
 using namespace Foam::constant;
 
@@ -102,6 +103,64 @@ Foam::SurfaceFilmModel<CloudType>::~SurfaceFilmModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+template<class CloudTrackType>
+void Foam::SurfaceFilmModel<CloudType>::injectParticles
+(
+    const label primaryPatchi,
+    const labelList& injectorCellsPatch,
+    CloudTrackType& cloud
+)
+{
+    const fvMesh& mesh = this->owner().mesh();
+    const vectorField& Cf = mesh.C().boundaryField()[primaryPatchi];
+    const vectorField& Sf = mesh.Sf().boundaryField()[primaryPatchi];
+    const scalarField& magSf =
+        mesh.magSf().boundaryField()[primaryPatchi];
+
+    forAll(injectorCellsPatch, j)
+    {
+        if (diameterParcelPatch_[j] > 0)
+        {
+            const label celli = injectorCellsPatch[j];
+
+            const scalar offset =
+                max
+                (
+                    diameterParcelPatch_[j],
+                    deltaFilmPatch_[primaryPatchi][j]
+                );
+            const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j];
+
+            // Create a new parcel
+            parcelType* pPtr =
+                new parcelType(this->owner().pMesh(), pos, celli);
+
+            // Check/set new parcel thermo properties
+            cloud.setParcelThermoProperties(*pPtr, 0.0);
+
+            setParcelProperties(*pPtr, j);
+
+            if (pPtr->nParticle() > 0.001)
+            {
+                // Check new parcel properties
+//                cloud.checkParcelProperties(*pPtr, 0.0, true);
+                cloud.checkParcelProperties(*pPtr, 0.0, false);
+
+                // Add the new parcel to the cloud
+                cloud.addParticle(pPtr);
+
+                nParcelsInjected_++;
+            }
+            else
+            {
+                // TODO: cache mass and re-distribute?
+                delete pPtr;
+            }
+        }
+    }
+}
+
 template<class CloudType>
 template<class TrackCloudType>
 void Foam::SurfaceFilmModel<CloudType>::inject(TrackCloudType& cloud)
@@ -111,78 +170,63 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackCloudType& cloud)
         return;
     }
 
+    const fvMesh& mesh = this->owner().mesh();
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
     // Retrieve the film model from the owner database
-    const regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel =
-        this->owner().mesh().time().objectRegistry::template lookupObject
+    const regionModels::surfaceFilmModels::surfaceFilmRegionModel* filmModel =
+        mesh.time().objectRegistry::template findObject
         <regionModels::surfaceFilmModels::surfaceFilmRegionModel>
         (
             "surfaceFilmProperties"
         );
 
-    if (!filmModel.active())
+    // Check the singleLayer type of films
+    if (filmModel && filmModel->active())
     {
-        return;
-    }
 
-    const labelList& filmPatches = filmModel.intCoupledPatchIDs();
-    const labelList& primaryPatches = filmModel.primaryPatchIDs();
+        const labelList& filmPatches = filmModel->intCoupledPatchIDs();
+        const labelList& primaryPatches = filmModel->primaryPatchIDs();
 
-    const fvMesh& mesh = this->owner().mesh();
-    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+        forAll(filmPatches, i)
+        {
+            const label filmPatchi = filmPatches[i];
+            const label primaryPatchi = primaryPatches[i];
 
-    forAll(filmPatches, i)
-    {
-        const label filmPatchi = filmPatches[i];
-        const label primaryPatchi = primaryPatches[i];
+            const labelList& injectorCellsPatch = pbm[primaryPatchi].faceCells();
 
-        const labelList& injectorCellsPatch = pbm[primaryPatchi].faceCells();
+            cacheFilmFields(filmPatchi, primaryPatchi, *filmModel);
 
-        cacheFilmFields(filmPatchi, primaryPatchi, filmModel);
+            injectParticles(primaryPatchi, injectorCellsPatch, cloud);
+        }
+    }
 
-        const vectorField& Cf = mesh.C().boundaryField()[primaryPatchi];
-        const vectorField& Sf = mesh.Sf().boundaryField()[primaryPatchi];
-        const scalarField& magSf = mesh.magSf().boundaryField()[primaryPatchi];
+    // Check finite area films
+    wordList names =
+        mesh.time().objectRegistry::template
+            sortedNames<regionModels::regionFaModel>();
 
-        forAll(injectorCellsPatch, j)
+    forAll (names, i)
+    {
+        const regionModels::regionFaModel* regionFa =
+            mesh.time().objectRegistry::template cfindObject
+            <
+                regionModels::regionFaModel
+            >(names[i]);
+
+        // Check that it is a type areaFilm
+        if (regionFa && isA<areaFilm>(*regionFa))
         {
-            if (diameterParcelPatch_[j] > 0)
-            {
-                const label celli = injectorCellsPatch[j];
-
-                const scalar offset =
-                    max
-                    (
-                        diameterParcelPatch_[j],
-                        deltaFilmPatch_[primaryPatchi][j]
-                    );
-                const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j];
-
-                // Create a new parcel
-                parcelType* pPtr =
-                    new parcelType(this->owner().pMesh(), pos, celli);
-
-                // Check/set new parcel thermo properties
-                cloud.setParcelThermoProperties(*pPtr, 0.0);
-
-                setParcelProperties(*pPtr, j);
-
-                if (pPtr->nParticle() > 0.001)
-                {
-                    // Check new parcel properties
-    //                cloud.checkParcelProperties(*pPtr, 0.0, true);
-                    cloud.checkParcelProperties(*pPtr, 0.0, false);
-
-                    // Add the new parcel to the cloud
-                    cloud.addParticle(pPtr);
-
-                    nParcelsInjected_++;
-                }
-                else
-                {
-                    // TODO: cache mass and re-distribute?
-                    delete pPtr;
-                }
-            }
+            areaFilm& film =
+                const_cast<areaFilm&>(refCast<const areaFilm>(*regionFa));
+
+            const label patchId = regionFa->patchID();
+
+            const labelList& injectorCellsPatch = pbm[patchId].faceCells();
+
+            cacheFilmFields(patchId, film);
+
+            injectParticles(patchId, injectorCellsPatch, cloud);
         }
     }
 }
@@ -215,6 +259,35 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
 }
 
 
+template<class CloudType>
+void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
+(
+    const label filmPatchi,
+    const regionModels::areaSurfaceFilmModels::liquidFilmBase& filmModel
+)
+{
+    // WIP: the finite area film does not support splashing so far
+    //massParcelPatch_ = filmModel.cloudMassTrans().boundaryField()[filmPatchi];
+    //filmModel.toPrimary(filmPatchi, massParcelPatch_);
+    massParcelPatch_.setSize(filmModel.Uf().size(), Zero);
+    diameterParcelPatch_.setSize(filmModel.Uf().size(), Zero);
+    //    filmModel.cloudDiameterTrans().boundaryField()[filmPatchi];
+    //filmModel.toPrimary(filmPatchi, diameterParcelPatch_, maxEqOp<scalar>());
+    const volSurfaceMapping& map = filmModel.region().vsm();
+
+    UFilmPatch_.setSize(filmModel.Uf().size(), Zero);
+    map.mapToField(filmModel.Uf(), UFilmPatch_);
+    //filmModel.toPrimary(filmPatchi, UFilmPatch_);
+
+    rhoFilmPatch_.setSize(UFilmPatch_.size(), Zero);
+    map.mapToField(filmModel.rho(), rhoFilmPatch_);
+    //filmModel.toPrimary(filmPatchi, rhoFilmPatch_);
+
+    deltaFilmPatch_[filmPatchi].setSize(UFilmPatch_.size(), Zero);
+    map.mapToField(filmModel.h(), deltaFilmPatch_[filmPatchi]);
+}
+
+
 template<class CloudType>
 void Foam::SurfaceFilmModel<CloudType>::setParcelProperties
 (
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
index ed232999065c5470add13bed5ac797f4797df9e5..0e5a7a4cef20afb87512b37f4184e75f86fcfec5 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
@@ -60,6 +60,15 @@ namespace regionModels
     }
 }
 
+namespace regionModels
+{
+    namespace areaSurfaceFilmModels
+    {
+        class liquidFilmBase;
+    }
+}
+
+
 /*---------------------------------------------------------------------------*\
                       Class SurfaceFilmModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -76,6 +85,9 @@ protected:
         //- Convenience typedef to the cloud's parcel type
         typedef typename CloudType::parcelType parcelType;
 
+        typedef typename
+            regionModels::areaSurfaceFilmModels::liquidFilmBase areaFilm;
+
         //- Gravitational acceleration constant
         const dimensionedVector& g_;
 
@@ -94,13 +106,13 @@ protected:
             scalarList diameterParcelPatch_;
 
             //- Film velocity / patch face
-            List<vector> UFilmPatch_;
+            Field<vector> UFilmPatch_;
 
             //- Film density / patch face
-            scalarList rhoFilmPatch_;
+            scalarField rhoFilmPatch_;
 
             //- Film height of all film patches / patch face
-            scalarListList deltaFilmPatch_;
+            Field<scalarField> deltaFilmPatch_;
 
 
         // Counters
@@ -122,6 +134,22 @@ protected:
             const regionModels::surfaceFilmModels::surfaceFilmRegionModel&
         );
 
+        //- Cache the finite area film fields in preparation for injection
+        virtual void cacheFilmFields
+        (
+            const label primaryPatchi,
+            const regionModels::areaSurfaceFilmModels::liquidFilmBase&
+        );
+
+        //- Inject particles in cloud
+        template<class TrackCloudType>
+        void injectParticles
+        (
+            const label primaryPatchi,
+            const labelList& injectorCellsPatch,
+            TrackCloudType& cloud
+        );
+
         //- Set the individual parcel properties
         virtual void setParcelProperties
         (
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
index d45be7e0427a3851c074e7f28ab6c540394ce4dd..70594de693fe43b4c837154f83824b7cdda31abf 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
@@ -27,6 +27,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ThermoSurfaceFilm.H"
+#include "surfaceFilmRegionModel.H"
+#include "liquidFilmBase.H"
 #include "addToRunTimeSelectionTable.H"
 #include "unitConversion.H"
 #include "Pstream.H"
@@ -131,7 +133,6 @@ Foam::vector Foam::ThermoSurfaceFilm<CloudType>::splashDirection
 template<class CloudType>
 void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
 (
-    regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -159,6 +160,56 @@ void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
     // Parcel tangential velocity
     const vector Ut = Urel - Un;
 
+    if (filmModel_)
+    {
+        filmModel_->addSources
+        (
+            pp.index(),
+            facei,
+            mass,                           // mass
+            mass*Ut,                        // tangential momentum
+            mass*mag(Un),                   // impingement pressure
+            mass*p.hs()                     // energy
+        );
+
+        this->nParcelsTransferred()++;
+
+        keepParticle = false;
+    }
+}
+
+template<class CloudType>
+void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
+(
+    areaFilm& filmModel,
+    const parcelType& p,
+    const polyPatch& pp,
+    const label facei,
+    const scalar mass,
+    bool& keepParticle
+)
+{
+
+    if (debug)
+    {
+        Info<< "Parcel " << p.origId() << " absorbInteraction" << endl;
+    }
+
+    // Patch face normal
+    const vector& nf = pp.faceNormals()[facei];
+
+    // Patch velocity
+    const vector& Up = this->owner().U().boundaryField()[pp.index()][facei];
+
+    // Relative parcel velocity
+    const vector Urel = p.U() - Up;
+
+    // Parcel normal velocity
+    const vector Un = nf*(Urel & nf);
+
+    // Parcel tangential velocity
+    const vector Ut = Urel - Un;
+DebugVar(mass)
     filmModel.addSources
     (
         pp.index(),
@@ -172,6 +223,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
     this->nParcelsTransferred()++;
 
     keepParticle = false;
+
 }
 
 
@@ -208,7 +260,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::bounceInteraction
 template<class CloudType>
 void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
 (
-    regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -249,14 +301,14 @@ void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
 
     if (We < Wec) // Adhesion - assume absorb
     {
-        absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
+        absorbInteraction(p, pp, facei, m, keepParticle);
     }
     else // Splash
     {
         // Ratio of incident mass to splashing mass
         const scalar mRatio = 0.2 + 0.6*rndGen_.sample01<scalar>();
         splashInteraction
-            (filmModel, p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
+            (p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
     }
 }
 
@@ -264,7 +316,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
 template<class CloudType>
 void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
 (
-    regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
     parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -307,7 +359,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
 
     if (We < 2) // Adhesion - assume absorb
     {
-        absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
+        absorbInteraction(p, pp, facei, m, keepParticle);
     }
     else if ((We >= 2) && (We < 20)) // Bounce
     {
@@ -325,7 +377,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
     }
     else if ((We >= 20) && (We < Wec)) // Spread - assume absorb
     {
-        absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
+        absorbInteraction(p, pp, facei, m, keepParticle);
     }
     else    // Splash
     {
@@ -333,7 +385,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
         // splash mass can be > incident mass due to film entrainment
         const scalar mRatio = 0.2 + 0.9*rndGen_.sample01<scalar>();
         splashInteraction
-            (filmModel, p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
+            (p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
     }
 }
 
@@ -341,7 +393,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
 template<class CloudType>
 void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
 (
-    regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -414,7 +466,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     // Switch to absorb if insufficient energy for splash
     if (EKs <= 0)
     {
-        absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
+        absorbInteraction(p, pp, facei, m, keepParticle);
         return;
     }
 
@@ -469,7 +521,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     // Transfer remaining part of parcel to film 0 - splashMass can be -ve
     // if entraining from the film
     const scalar mDash = m - mSplash;
-    absorbInteraction(filmModel, p, pp, facei, mDash, keepParticle);
+    absorbInteraction(p, pp, facei, mDash, keepParticle);
 }
 
 
@@ -488,6 +540,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
     (
         owner.db().objectRegistry::template lookupObject<SLGThermo>("SLGThermo")
     ),
+    filmModel_(nullptr),
     TFilmPatch_(0),
     CpFilmPatch_(0),
     interactionType_
@@ -516,6 +569,18 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
         this->coeffDict().readEntry("Awet", Awet_);
         this->coeffDict().readEntry("Cf", Cf_);
     }
+
+    filmModel_ =
+        const_cast<regionModels::surfaceFilmModels::surfaceFilmRegionModel*>
+        (
+            this->owner().mesh().time().objectRegistry::template findObject
+            <
+                regionModels::surfaceFilmModels::surfaceFilmRegionModel
+            >
+            (
+                "surfaceFilmProperties"
+            )
+        );
 }
 
 
@@ -528,6 +593,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
     SurfaceFilmModel<CloudType>(sfm),
     rndGen_(sfm.rndGen_),
     thermo_(sfm.thermo_),
+    filmModel_(nullptr),
     TFilmPatch_(sfm.TFilmPatch_),
     CpFilmPatch_(sfm.CpFilmPatch_),
     interactionType_(sfm.interactionType_),
@@ -538,7 +604,19 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
     Awet_(sfm.Awet_),
     Cf_(sfm.Cf_),
     nParcelsSplashed_(sfm.nParcelsSplashed_)
-{}
+{
+     filmModel_ =
+        const_cast<regionModels::surfaceFilmModels::surfaceFilmRegionModel*>
+        (
+            this->owner().mesh().time().objectRegistry::template findObject
+            <
+                regionModels::surfaceFilmModels::surfaceFilmRegionModel
+            >
+            (
+                "surfaceFilmProperties"
+            )
+        );
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -558,21 +636,12 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
     bool& keepParticle
 )
 {
-    // Retrieve the film model from the owner database
-    regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel =
-        const_cast<regionModels::surfaceFilmModels::surfaceFilmRegionModel&>
-        (
-            this->owner().db().time().objectRegistry::template
-                lookupObject
-                <regionModels::surfaceFilmModels::surfaceFilmRegionModel>
-                (
-                    "surfaceFilmProperties"
-                )
-        );
+    const fvMesh& mesh = this->owner().mesh();
 
     const label patchi = pp.index();
 
-    if (filmModel.isRegionPatch(patchi))
+    // Check the singleLayer film models
+    if (filmModel_ && filmModel_->isRegionPatch(patchi))
     {
         const label facei = pp.whichFace(p.face());
 
@@ -587,7 +656,7 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
             case itAbsorb:
             {
                 const scalar m = p.nParticle()*p.mass();
-                absorbInteraction(filmModel, p, pp, facei, m, keepParticle);
+                absorbInteraction(p, pp, facei, m, keepParticle);
 
                 break;
             }
@@ -597,11 +666,11 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
 
                 if (dry)
                 {
-                    drySplashInteraction(filmModel, p, pp, facei, keepParticle);
+                    drySplashInteraction(p, pp, facei, keepParticle);
                 }
                 else
                 {
-                    wetSplashInteraction(filmModel, p, pp, facei, keepParticle);
+                    wetSplashInteraction( p, pp, facei, keepParticle);
                 }
 
                 break;
@@ -617,6 +686,54 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
         // Transfer parcel/parcel interactions complete
         return true;
     }
+    else
+    {
+
+        // Check the finite area films
+        wordList names =
+            mesh.time().objectRegistry::template
+                sortedNames<regionModels::regionFaModel>();
+
+        forAll(names, i)
+        {
+            const regionModels::regionFaModel* regionFa =
+                mesh.time().objectRegistry::template findObject
+                <
+                    regionModels::regionFaModel
+                >(names[i]);
+
+            // Check that regionFa is a areaFilm
+            if (regionFa && isA<areaFilm>(*regionFa))
+            {
+                areaFilm& film =
+                    const_cast<areaFilm&>(refCast<const areaFilm>(*regionFa));
+
+                const label facei = pp.whichFace(p.face());
+
+                switch (interactionType_)
+                {
+                    // It only supports absorp model
+                    case itAbsorb:
+                    {
+                        const scalar m = p.nParticle()*p.mass();
+                        absorbInteraction
+                        (
+                            film, p, pp, facei, m, keepParticle
+                        );
+                        break;
+                    }
+                    default:
+                    {
+                        FatalErrorInFunction
+                            << "Unknown interaction type enumeration"
+                            << abort(FatalError);
+                    }
+                }
+                // Transfer parcel/parcel interactions complete
+                return true;
+            }
+        }
+    }
 
     // Parcel not interacting with film
     return false;
@@ -646,6 +763,26 @@ void Foam::ThermoSurfaceFilm<CloudType>::cacheFilmFields
 }
 
 
+template<class CloudType>
+void Foam::ThermoSurfaceFilm<CloudType>::cacheFilmFields
+(
+    const label filmPatchi,
+    const areaFilm& filmModel
+)
+{
+    SurfaceFilmModel<CloudType>::cacheFilmFields
+    (
+        filmPatchi,
+        filmModel
+    );
+    const volSurfaceMapping& map = filmModel.region().vsm();
+
+    TFilmPatch_.setSize(filmModel.Tf().size(), Zero);
+    map.mapToField(filmModel.Tf(), TFilmPatch_);
+    //map.mapToField(filmModel.Cp(), CpFilmPatch_);
+}
+
+
 template<class CloudType>
 void Foam::ThermoSurfaceFilm<CloudType>::setParcelProperties
 (
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
index 186e951dbb32190560873aff17158b149be85e98..b1b821a2944cba994ae1591ce5816bc2bf7968c9 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
@@ -61,6 +61,22 @@ SourceFiles
 namespace Foam
 {
 
+//Forward declaration of classes
+namespace regionModels
+{
+    namespace surfaceFilmModels
+    {
+        class surfaceFilmRegionModel;
+    }
+}
+namespace regionModels
+{
+    namespace areaSurfaceFilmModels
+    {
+        class liquidFilmBase;
+    }
+}
+
 /*---------------------------------------------------------------------------*\
                       Class ThermoSurfaceFilm Declaration
 \*---------------------------------------------------------------------------*/
@@ -102,20 +118,26 @@ protected:
         //- Convenience typedef to the cloud's parcel type
         typedef typename CloudType::parcelType parcelType;
 
+        typedef typename
+            regionModels::areaSurfaceFilmModels::liquidFilmBase areaFilm;
+
         //- Reference to the cloud random number generator
         Random& rndGen_;
 
         //- Reference to the cloud thermo package
         const SLGThermo& thermo_;
 
+        //- Pointer to filmModel
+        regionModels::surfaceFilmModels::surfaceFilmRegionModel* filmModel_;
+
 
         // Cached injector fields per film patch
 
             //- Film temperature / patch face
-            scalarList TFilmPatch_;
+            scalarField TFilmPatch_;
 
             //- Film specific heat capacity / patch face
-            scalarList CpFilmPatch_;
+            scalarField CpFilmPatch_;
 
 
         // Interaction model data
@@ -173,7 +195,18 @@ protected:
             //- Absorb parcel into film
             void absorbInteraction
             (
-                regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                const parcelType& p,
+                const polyPatch& pp,
+                const label facei,
+                const scalar mass,
+                bool& keepParticle
+            );
+
+            //- Absorb parcel into area film
+            void absorbInteraction
+            (
+                areaFilm&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -193,7 +226,7 @@ protected:
             //- Parcel interaction with dry surface
             void drySplashInteraction
             (
-                regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -203,7 +236,7 @@ protected:
             //- Parcel interaction with wetted surface
             void wetSplashInteraction
             (
-                regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
                 parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -213,7 +246,7 @@ protected:
             //- Bai parcel splash interaction model
             void splashInteraction
             (
-                regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -224,18 +257,24 @@ protected:
                 bool& keepParticle
             );
 
+            virtual void cacheFilmFields
+            (
+                const label primaryPatchi,
+                const areaFilm&
+            );
 
         // Injection from sheet (ejection) helper functions
 
-            //- Cache the film fields in preparation for injection
+             //- Cache the film fields in preparation for injection
             virtual void cacheFilmFields
             (
                 const label filmPatchi,
                 const label primaryPatchi,
                 const regionModels::surfaceFilmModels::surfaceFilmRegionModel&
-                    filmModel
             );
 
+
+
             //- Set the individual parcel properties
             virtual void setParcelProperties
             (
diff --git a/src/lagrangian/spray/Make/options b/src/lagrangian/spray/Make/options
index d595ada329b1a9965c39d7994b43dc7260be42f6..443d427398f16ebecef895d4549b1801342f586a 100644
--- a/src/lagrangian/spray/Make/options
+++ b/src/lagrangian/spray/Make/options
@@ -20,7 +20,10 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -44,4 +47,6 @@ LIB_LIBS = \
     -lregionModels \
     -lsurfaceFilmModels \
     -ldynamicMesh \
-    -ldynamicFvMesh
+    -ldynamicFvMesh \
+    -lregionFaModels \
+    -lfiniteArea
diff --git a/src/lagrangian/turbulence/Make/options b/src/lagrangian/turbulence/Make/options
index 7be9fe24ee3e03f04493320be429481c72dc17c9..e2fc84c325a18d7cc42ef1efe3dd87bd74d1def7 100644
--- a/src/lagrangian/turbulence/Make/options
+++ b/src/lagrangian/turbulence/Make/options
@@ -18,7 +18,10 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/regionFaModels/lnInclude \
+    -I$(LIB_SRC)/finiteArea/lnInclude \
+    -I$(LIB_SRC)/faOptions/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -40,4 +43,6 @@ LIB_LIBS = \
     -lincompressibleTransportModels \
     -lregionModels \
     -lsurfaceFilmModels \
-    -ldynamicFvMesh
+    -ldynamicFvMesh \
+    -lregionFaModels \
+    -lfiniteArea
diff --git a/src/regionFaModels/Make/files b/src/regionFaModels/Make/files
index 3a3b9a4adddcd3f2bc9b9f3dd700bd6740cb2f91..70521f391ad92c25a7ef89456845bd38ac426b55 100644
--- a/src/regionFaModels/Make/files
+++ b/src/regionFaModels/Make/files
@@ -13,4 +13,15 @@ KirchhoffShell/KirchhoffShell.C
 derivedFvPatchFields/thermalShell/thermalShellFvPatchScalarField.C
 derivedFvPatchFields/vibrationShell/vibrationShellFvPatchScalarField.C
 
+liquidFilm/liquidFilmBase.C
+liquidFilm/liquidFilmBaseNew.C
+liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
+liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
+liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
+liquidFilm/liquidFilmModel/liquidFilmModel.C
+liquidFilm/kinematicThinFilm/kinematicThinFilm.C
+derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.C
+
+functionObjects/setTimeStep/setTimeStepFaRegionsFunctionObject.C
+
 LIB = $(FOAM_LIBBIN)/libregionFaModels
diff --git a/src/regionFaModels/Make/options b/src/regionFaModels/Make/options
index cd4795b38d01292372b3c1514abe9745dcfa97ad..2af18fc4d144ec9149cc70d765f78d4c39882cf4 100644
--- a/src/regionFaModels/Make/options
+++ b/src/regionFaModels/Make/options
@@ -5,7 +5,12 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/transportModels
 
 LIB_LIBS = \
     -lfiniteVolume \
diff --git a/src/regionFaModels/derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.C b/src/regionFaModels/derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.C
new file mode 100644
index 0000000000000000000000000000000000000000..e2fe02f80fe3d55ad89af703000828bc7565eba4
--- /dev/null
+++ b/src/regionFaModels/derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.C
@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "velocityFilmShellFvPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    mixedFvPatchField<vector>(p, iF),
+    baffle_(),
+    dict_(dictionary::null),
+    curTimeIndex_(-1),
+    zeroWallVelocity_(true)
+{
+    refValue() = 0;
+    refGrad() = 0;
+    valueFraction() = 1;
+}
+
+
+velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
+(
+    const velocityFilmShellFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchField<vector>
+    (
+        ptf,
+        p,
+        iF,
+        mapper
+    ),
+    baffle_(),
+    dict_(ptf.dict_),
+    curTimeIndex_(-1),
+    zeroWallVelocity_(true)
+{}
+
+
+velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchField<vector>(p, iF),
+    baffle_(nullptr),
+    dict_(dict),
+    curTimeIndex_(-1),
+    zeroWallVelocity_(dict.getOrDefault<bool>("zeroWallVelocity", true))
+{
+    fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
+
+    typedef regionModels::areaSurfaceFilmModels::liquidFilmBase baffle;
+
+    if (dict.found("refValue"))
+    {
+        // Full restart
+        refValue() = vectorField("refValue", dict, p.size());
+        refGrad() = vectorField("refGradient", dict, p.size());
+        valueFraction() = scalarField("valueFraction", dict, p.size());
+    }
+    else
+    {
+        // Start from user entered data. Assume fixedValue.
+        refValue() = *this;
+        refGrad() = vector::zero;
+        valueFraction() = 1;
+    }
+
+    if (!baffle_)
+    {
+        baffle_.reset(baffle::New(p, dict).ptr());
+    }
+}
+
+
+velocityFilmShellFvPatchVectorField::velocityFilmShellFvPatchVectorField
+(
+    const velocityFilmShellFvPatchVectorField& ptf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    mixedFvPatchField<vector>(ptf, iF),
+    baffle_(),
+    dict_(ptf.dict_),
+    curTimeIndex_(-1),
+    zeroWallVelocity_(true)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void velocityFilmShellFvPatchVectorField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    // Execute the change to the openFraction only once per time-step
+    if (curTimeIndex_ != this->db().time().timeIndex())
+    {
+        baffle_->evolve();
+
+        volVectorField::Boundary& vfb =
+            db().lookupObjectRef<volVectorField>
+            (
+                this->internalField().name()
+            ).boundaryFieldRef();
+
+        baffle_->vsm().mapToVolume(baffle_->Us(), vfb);
+
+        refGrad() = Zero;
+        valueFraction() = 1;
+
+        if (zeroWallVelocity_)
+        {
+            refValue() = Zero;
+        }
+        else
+        {
+            refValue() = vfb[patch().index()];
+        }
+        curTimeIndex_ = this->db().time().timeIndex();
+    }
+
+    mixedFvPatchField<vector>::updateCoeffs();
+}
+
+
+void velocityFilmShellFvPatchVectorField::write(Ostream& os) const
+{
+    mixedFvPatchField<vector>::write(os);
+
+    // Remove value and type already written by mixedFvPatchField
+    dict_.remove("value");
+    dict_.remove("type");
+    dict_.remove("refValue");
+    dict_.remove("refGradient");
+    dict_.remove("valueFraction");
+    dict_.write(os, false);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchVectorField,
+    velocityFilmShellFvPatchVectorField
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.H b/src/regionFaModels/derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.H
new file mode 100644
index 0000000000000000000000000000000000000000..66ba8205ee435516c2b14b927747b405cc398480
--- /dev/null
+++ b/src/regionFaModels/derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.H
@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::compressible::velocityFilmShellFvPatchVectorField
+
+Group
+    grpLiquidFilmBoundaryConditions
+
+Description
+
+Usage
+    Example of the boundary condition specification:
+    \verbatim
+    <patchName>
+    {
+        type                velocityFilmShell;
+
+        active              true;
+        infoOutput          true;
+
+        U                   U;
+        pRef                1e5;
+        T0                  300;
+
+        deltaWet            1e-4;
+        h0                  1e-8;
+
+        zeroWallVelocity    true;
+
+        thermo
+        {
+            H2O;
+        }
+
+        turbulence          laminar;
+
+        laminarCoeffs
+        {
+            friction        ManningStrickler;   // Wall friction model
+            n               0.005;              // Manning number
+            Cf              0.9;                // Gas friction
+        }
+
+        injectionModels
+        (
+            curvatureSeparation
+        );
+
+        forces ();
+
+        curvatureSeparationCoeffs
+        {
+            definedPatchRadii  0;
+        }
+
+        region              film;
+        liquidFilmModel     kinematicThinFilm;
+
+        value               uniform (0 0 0);
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property    | Description                        | Type | Reqd | Deflt
+      type        | Type name: velocityFilmShell       | word | yes  | -
+      U           | Name of the primary U              | word | yes  | -
+      pRef        | Reference pressure for thermo    | scalar | yes  | -
+      T0          | Film initial temperature         | scalar | no   | READ
+      thermo      | Flow thermo                     | wordRes | yes  | -
+      zeroWallVelocity    | Flag to fix zero U for primary flow <!--
+                  -->                                  | bool | no   | true
+      turbulence  | Type of film turbulence model      | word | yes  | -
+      injectionModels   | Lagrangian injection         |      | no   | -
+      forces      | Film force models               | wordRes | no   | -
+      deltaWet    | Wet film thickness               | scalar | no   | 1e-4
+      h0          | Numerical minimum thickness      | scalar | no   | 1e-7
+      region      | Name of the 2D region            | word   | yes  | -
+      liquidFilmModel     | Film model               | word   | yes  | -
+    \endtable
+
+SourceFiles
+    velocityFilmShellFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef velocityFilmShellFvPatchVectorField_H
+#define velocityFilmShellFvPatchVectorField_H
+
+#include "autoPtr.H"
+#include "liquidFilmBase.H"
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+             Class velocityFilmShellFvPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class velocityFilmShellFvPatchVectorField
+:
+    public mixedFvPatchField<vector>
+{
+    // Private Data
+
+        //- Thermal baffle
+        autoPtr<regionModels::areaSurfaceFilmModels::liquidFilmBase> baffle_;
+
+        //- Dictionary
+        mutable dictionary dict_;
+
+        //- Time index to evolve the film
+        label curTimeIndex_;
+
+        //- Zero wall velocity. Fix U to zero or to film U
+        bool zeroWallVelocity_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("velocityFilmShell");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        velocityFilmShellFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        velocityFilmShellFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //- velocityFilmShellFvPatchVectorField onto a new patch
+        velocityFilmShellFvPatchVectorField
+        (
+            const velocityFilmShellFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new velocityFilmShellFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        velocityFilmShellFvPatchVectorField
+        (
+            const velocityFilmShellFvPatchVectorField&,
+            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 velocityFilmShellFvPatchVectorField(*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/regionFaModels/functionObjects/setTimeStep/setTimeStepFaRegionsFunctionObject.C b/src/regionFaModels/functionObjects/setTimeStep/setTimeStepFaRegionsFunctionObject.C
new file mode 100644
index 0000000000000000000000000000000000000000..3b5e95ffd1fae6229d36057b019e06bcf4b264e3
--- /dev/null
+++ b/src/regionFaModels/functionObjects/setTimeStep/setTimeStepFaRegionsFunctionObject.C
@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "setTimeStepFaRegionsFunctionObject.H"
+#include "addToRunTimeSelectionTable.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(setTimeStepFaRegionsFunctionObject, 0);
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        setTimeStepFaRegionsFunctionObject,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::
+setTimeStepFaRegionsFunctionObject::
+setTimeStepFaRegionsFunctionObject
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    timeFunctionObject(name, runTime)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::setTimeStepFaRegionsFunctionObject::adjustTimeStep()
+{
+    // Wanted timestep
+    scalar newDeltaT = regionDeltaT();
+
+    static label index = -1;
+
+    if ((time_.timeIndex() != index) && (newDeltaT < time_.deltaTValue()))
+    {
+        // Store current time so we don't get infinite recursion (since
+        // setDeltaT calls adjustTimeStep() again)
+        index = time_.timeIndex();
+
+        // Set time, allow deltaT to be adjusted for writeInterval purposes
+        const_cast<Time&>(time_).setDeltaT(newDeltaT, false);
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::setTimeStepFaRegionsFunctionObject::read
+(
+    const dictionary& dict
+)
+{
+    if (timeFunctionObject::read(dict))
+    {
+        // Ensure that adjustTimeStep is active
+        if (!time_.controlDict().lookupOrDefault<bool>("adjustTimeStep", false))
+        {
+            FatalIOErrorInFunction(dict)
+                << "Need to set 'adjustTimeStep' true to allow timestep control"
+                << nl
+                << exit(FatalIOError);
+        }
+
+        return true;
+    }
+
+    return false;
+}
+
+
+Foam::scalar Foam::functionObjects::setTimeStepFaRegionsFunctionObject::
+regionDeltaT() const
+{
+    const wordList names(time_.sortedNames<regionFaModel>());
+
+    scalar Co = 0.0;
+
+    forAll (names, i)
+    {
+        const auto* regionFa = time_.cfindObject<regionFaModel>(names[i]);
+
+        if (regionFa)
+        {
+            const scalar regionCo = regionFa->CourantNumber();
+            if (regionCo > Co)
+            {
+                Co = regionCo;
+            }
+        }
+    }
+
+    if (names.size() > 0)
+    {
+        const scalar regionFaMaxCo =
+            time_.controlDict().get<scalar>("regionFaMaxCo");
+
+        const scalar maxDeltaTFact = regionFaMaxCo/(Co + SMALL);
+        const scalar deltaTFact =
+            min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+        return deltaTFact*time_.deltaTValue();
+    }
+
+    return time_.deltaTValue();
+}
+
+
+bool Foam::functionObjects::setTimeStepFaRegionsFunctionObject::execute()
+{
+    return true;
+}
+
+
+bool Foam::functionObjects::setTimeStepFaRegionsFunctionObject::write()
+{
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/functionObjects/setTimeStep/setTimeStepFaRegionsFunctionObject.H b/src/regionFaModels/functionObjects/setTimeStep/setTimeStepFaRegionsFunctionObject.H
new file mode 100644
index 0000000000000000000000000000000000000000..6c61535720b9f129e55f0e1de44d85db5ad5bc7c
--- /dev/null
+++ b/src/regionFaModels/functionObjects/setTimeStep/setTimeStepFaRegionsFunctionObject.H
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::functionObjects::setTimeStepFaRegionsFunctionObject
+
+Group
+    grpUtilitiesFunctionObjects
+
+Description
+    This function object controls the time step for classes of the type
+    \c regionFaModel. It reads \c regionFaMaxCo entry from \c controlDict
+    and evaluate the time step based on the finite area Courant Number.
+
+    Can only be used with solvers using \c adjustTimeStep control (e.g.
+    \c pimpleFoam). It makes no attempt to co-operate with other time step
+    'controllers', e.g. \c maxCo, other functionObjects. Supports \c enabled
+    flag but none of the other options \c timeStart, \c timeEnd, \c writeControl
+    etc.
+
+Usage
+    Example of function object specification to manipulate the time step:
+    \verbatim
+    setTimeStep1
+    {
+        // Mandatory entries
+        type        setTimeStepFaRegion;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                    | Type | Reqd | Deflt
+      type         | Type name: setTimeStepFaRegion | word | yes  | -
+      enabled      | On/off switch                  | bool | no   | yes
+    \endtable
+
+    The inherited entries are elaborated in:
+     - \link timeFunctionObject.H \endlink
+     - \link regionFaModel.H \endlink
+
+SourceFiles
+    setTimeStepFaRegionsFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_setTimeStepFaRegionsFunctionObject_H
+#define functionObjects_setTimeStepFaRegionsFunctionObject_H
+
+#include "timeFunctionObject.H"
+#include "regionFaModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+using namespace Foam::regionModels;
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                Class setTimeStepFaRegionsFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class setTimeStepFaRegionsFunctionObject
+:
+    public functionObjects::timeFunctionObject
+{
+    // Private Member Functions
+
+        //- No copy construct
+        setTimeStepFaRegionsFunctionObject
+        (
+            const setTimeStepFaRegionsFunctionObject&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const setTimeStepFaRegionsFunctionObject&) = delete;
+
+        //- Return minimum deltaT from fa regions
+        scalar regionDeltaT() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("setTimeStepFaRegion");
+
+
+    // Constructors
+
+        //- Construct from components
+        setTimeStepFaRegionsFunctionObject
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+    virtual ~setTimeStepFaRegionsFunctionObject() = default;
+
+
+    // Member Functions
+
+        //- Called at the end of Time::adjustDeltaT() if adjustTime is true
+        virtual bool adjustTimeStep();
+
+        //- Read and set the function object if its data have changed
+        virtual bool read(const dictionary& dict);
+
+        //- Execute does nothing
+        virtual bool execute();
+
+        //- Write does nothing
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/kinematicThinFilm/kinematicThinFilm.C b/src/regionFaModels/liquidFilm/kinematicThinFilm/kinematicThinFilm.C
new file mode 100644
index 0000000000000000000000000000000000000000..1fc4b873ab7170332fb2aa0e5bbe1aea9d8a79c6
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/kinematicThinFilm/kinematicThinFilm.C
@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "kinematicThinFilm.H"
+#include "addToRunTimeSelectionTable.H"
+#include "uniformDimensionedFields.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(kinematicThinFilm, 0);
+addToRunTimeSelectionTable(liquidFilmBase, kinematicThinFilm, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+kinematicThinFilm::kinematicThinFilm
+(
+    const word& modelType,
+    const fvPatch& patch,
+    const dictionary& dict
+)
+:
+    liquidFilmModel(modelType, patch, dict)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void kinematicThinFilm::preEvolveRegion()
+{
+    rhoSp_.storePrevIter();
+    USp_.storePrevIter();
+    pnSp_.storePrevIter();
+
+    // Update mass exchange sources
+    liquidFilmModel::preEvolveRegion();
+
+    // gas pressure map from primary region
+    ppf_ = pg();
+}
+
+
+void kinematicThinFilm::evolveRegion()
+{
+    if (debug)
+    {
+        InfoInFunction << endl;
+    }
+
+    const areaVectorField& ns = regionMesh().faceAreaNormals();
+
+    const areaVectorField gs(g_ - ns*(ns & g_));
+
+    phi2s_ = fac::interpolate(h_)*phif_;
+
+    for (int oCorr=1; oCorr<=nOuterCorr_; oCorr++)
+    {
+        pf_.storePrevIter();
+
+        faVectorMatrix UsEqn
+        (
+            fam::ddt(h_, Uf_)
+          + fam::div(phi2s_, Uf_)
+          ==
+            gs*h_
+          + turbulence_->Su(Uf_)
+          + faOptions()(h_, Uf_, sqr(dimVelocity))
+          + forces_.correct(Uf_)
+          + USp_
+        );
+
+        UsEqn.relax();
+
+        faOptions().constrain(UsEqn);
+
+        if (momentumPredictor_)
+        {
+            solve(UsEqn == - fac::grad(pf_*h_)/rho_ + pf_*fac::grad(h_)/rho_);
+        }
+
+        for (int corr=1; corr<=nCorr_; corr++)
+        {
+            areaScalarField UsA(UsEqn.A());
+
+            Uf_ = UsEqn.H()/UsA;
+            Uf_.correctBoundaryConditions();
+            faOptions().correct(Uf_);
+
+            phif_ =
+                (fac::interpolate(Uf_) & regionMesh().Le())
+                - fac::interpolate(1.0/(rho_*UsA))
+                * fac::lnGrad(pf_*h_)*regionMesh().magLe()
+                + fac::interpolate(pf_/(rho_*UsA))
+                * fac::lnGrad(h_)*regionMesh().magLe();
+
+            for (int nFilm=1; nFilm<=nFilmCorr_; nFilm++)
+            {
+                faScalarMatrix hEqn
+                (
+                    fam::ddt(h_)
+                  + fam::div(phif_, h_)
+                 ==
+                    faOptions()(rho_, h_, dimVelocity)
+                  + rhoSp_
+                );
+
+                hEqn.relax();
+                faOptions().constrain(hEqn);
+                hEqn.solve();
+                faOptions().correct(h_);
+
+                if (nFilm == nFilmCorr_)
+                {
+                    phi2s_ = hEqn.flux();
+                }
+            }
+
+            // Bound h_
+            h_ = max(h_, h0_);
+
+            pf_ = rho_*gn_*h_ - sigma_*fac::laplacian(h_) + pnSp_ + ppf_;
+            pf_.correctBoundaryConditions();
+            pf_.relax();
+
+            Uf_ -= (1.0/(rho_*UsA))*fac::grad(pf_*h_)
+                - (pf_/(rho_*UsA))*fac::grad(h_);
+            Uf_.correctBoundaryConditions();
+            faOptions().correct(Uf_);
+        }
+    }
+
+     Info<< "Film h min/max   = " << min(h_).value() << ", "
+        << max(h_).value() << endl;
+
+     Info<< "Film U min/max   = " << max(mag(Uf_)).value() << endl;
+}
+
+
+void kinematicThinFilm::postEvolveRegion()
+{
+    // Reset sources
+    liquidFilmModel::postEvolveRegion();
+
+    // Correct thermo
+    correctThermoFields();
+
+    // Correct turbulence
+    turbulence_->correct();
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/kinematicThinFilm/kinematicThinFilm.H b/src/regionFaModels/liquidFilm/kinematicThinFilm/kinematicThinFilm.H
new file mode 100644
index 0000000000000000000000000000000000000000..bad60d49548336a46e13ce4f2140c635624b2707
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/kinematicThinFilm/kinematicThinFilm.H
@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::regionFaModels::kinematicThinFilm
+
+Description
+    Thin film model.
+
+SourceFiles
+    kinematicThinFilm.C
+    kinematicThinFilmI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kinematicThinFilm_H
+#define kinematicThinFilm_H
+
+#include "volFieldsFwd.H"
+#include "liquidFilmModel.H"
+#include "faMesh.H"
+#include "filmTurbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class kinematicThinFilm Declaration
+\*---------------------------------------------------------------------------*/
+
+class kinematicThinFilm
+:
+    public liquidFilmModel
+{
+public:
+
+    //- Runtime type information
+    TypeName("kinematicThinFilm");
+
+
+    // Constructors
+
+        //- Construct from components and dict
+        kinematicThinFilm
+        (
+            const word& modelType,
+            const fvPatch& patch,
+            const dictionary& dict
+        );
+
+        //- No copy construct
+        kinematicThinFilm(const kinematicThinFilm&) = delete;
+
+        //- No copy assignment
+        void operator=(const kinematicThinFilm&) = delete;
+
+
+    //- Destructor
+    virtual ~kinematicThinFilm() = default;
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Pre-evolve film
+            virtual void preEvolveRegion();
+
+            //- Evolve the film
+            virtual void evolveRegion();
+
+            //- Post-evolve film
+            virtual void postEvolveRegion();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/liquidFilmBase.C b/src/regionFaModels/liquidFilm/liquidFilmBase.C
new file mode 100644
index 0000000000000000000000000000000000000000..540bbe9da95c1d1e14807c34c75ce4840031020c
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/liquidFilmBase.C
@@ -0,0 +1,566 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "liquidFilmBase.H"
+#include "faMesh.H"
+#include "faCFD.H"
+#include "uniformDimensionedFields.H"
+#include "gravityMeshObject.H"
+#include "movingWallVelocityFvPatchVectorField.H"
+#include "turbulentFluidThermoModel.H"
+#include "turbulentTransportModel.H"
+#include "calculatedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(liquidFilmBase, 0);
+
+defineRunTimeSelectionTable(liquidFilmBase, dictionary);
+
+const Foam::word liquidFilmName("liquidFilm");
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+
+bool liquidFilmBase::read(const dictionary& dict)
+{
+    regionFaModel::read(dict);
+    if (active_)
+    {
+        const dictionary& solution = this->solution().subDict("PISO");
+        solution.readEntry("momentumPredictor", momentumPredictor_);
+        solution.readIfPresent("nOuterCorr", nOuterCorr_);
+        solution.readEntry("nCorr", nCorr_);
+        solution.readEntry("nNonOrthCorr", nNonOrthCorr_);
+    }
+    return true;
+}
+
+
+scalar liquidFilmBase::CourantNumber() const
+{
+    scalar CoNum = 0.0;
+    scalar meanCoNum = 0.0;
+    scalar velMag = 0.0;
+
+    if (regionMesh().nInternalEdges())
+    {
+        edgeScalarField SfUfbyDelta
+        (
+            regionMesh().edgeInterpolation::deltaCoeffs()*mag(phif_)
+        );
+
+        CoNum = max(SfUfbyDelta/regionMesh().magLe())
+            .value()*time().deltaT().value();
+
+        meanCoNum = (sum(SfUfbyDelta)/sum(regionMesh().magLe()))
+            .value()*time().deltaT().value();
+
+        velMag = max(mag(phif_)/regionMesh().magLe()).value();
+    }
+
+    Info<< "Film Courant Number mean: " << meanCoNum
+        << " max: " << CoNum
+        << " Film velocity magnitude: " << velMag << endl;
+
+    return CoNum;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+liquidFilmBase::liquidFilmBase
+(
+    const word& modelType,
+    const fvPatch& p,
+    const dictionary& dict
+)
+:
+    regionFaModel(p, liquidFilmName, modelType, dict, true),
+
+    momentumPredictor_
+    (
+        this->solution().subDict("PISO").lookup("momentumPredictor")
+    ),
+    nOuterCorr_
+    (
+        this->solution().subDict("PISO").lookupOrDefault("nOuterCorr", 1)
+    ),
+    nCorr_(this->solution().subDict("PISO").get<label>("nCorr")),
+    nNonOrthCorr_
+    (
+        this->solution().subDict("PISO").get<label>("nNonOrthCorr")
+    ),
+
+    h0_("h0", dimLength, SMALL),
+
+    UName_(dict.get<word>("U")),
+
+    pName_(dict.lookupOrDefault<word>("p",  word::null)),
+
+    pRef_(dict.get<scalar>("pRef")),
+
+    h_
+    (
+        IOobject
+        (
+            "hf_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh()
+    ),
+
+    Uf_
+    (
+        IOobject
+        (
+            "Uf_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh()
+    ),
+    pf_
+    (
+        IOobject
+        (
+            "pf_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimPressure, Zero)
+    ),
+    ppf_
+    (
+        IOobject
+        (
+            "ppf_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimPressure, Zero)
+    ),
+    phif_
+    (
+        IOobject
+        (
+            "phif_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fac::interpolate(Uf_) & regionMesh().Le()
+    ),
+
+    phi2s_
+    (
+        IOobject
+        (
+            "phi2s_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        fac::interpolate(h_*Uf_) & regionMesh().Le()
+    ),
+
+    Tf_
+    (
+        IOobject
+        (
+            "Tf_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh()
+    ),
+
+    gn_
+    (
+        IOobject
+        (
+            "gn",
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimAcceleration, Zero)
+    ),
+
+    massSource_
+    (
+        IOobject
+        (
+            "massSource",
+            primaryMesh().time().timeName(),
+            primaryMesh()
+        ),
+        primaryMesh(),
+        dimensionedScalar(dimMass, Zero),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+
+    momentumSource_
+    (
+        IOobject
+        (
+            "momentumSource",
+            primaryMesh().time().timeName(),
+            primaryMesh()
+        ),
+        primaryMesh(),
+        dimensionedVector(dimPressure, Zero),
+        calculatedFvPatchField<vector>::typeName
+    ),
+
+    pnSource_
+    (
+        IOobject
+        (
+            "pnSource",
+            primaryMesh().time().timeName(),
+            primaryMesh()
+        ),
+        primaryMesh(),
+        dimensionedScalar(dimPressure, Zero),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+
+    energySource_
+    (
+        IOobject
+        (
+            "energySource",
+            primaryMesh().time().timeName(),
+            primaryMesh()
+        ),
+        primaryMesh(),
+        dimensionedScalar(dimEnergy, Zero),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+
+    addedMassTotal_(0),
+
+    faOptions_(Foam::fa::options::New(p))
+{
+    const uniformDimensionedVectorField& g =meshObjects::gravity::New(time());
+
+    gn_ = (g & regionMesh().faceAreaNormals());
+
+
+    if (!faOptions_.optionList::size())
+    {
+        Info << "No finite area options present" << endl;
+    }
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+liquidFilmBase::~liquidFilmBase()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::areaVectorField> liquidFilmBase::Uw() const
+{
+    tmp<areaVectorField> tUw
+    (
+        new areaVectorField
+        (
+            IOobject
+            (
+                "tUw",
+                primaryMesh().time().timeName(),
+                primaryMesh()
+            ),
+            regionMesh(),
+            dimensionedVector(dimVelocity, Zero)
+        )
+    );
+
+    areaVectorField& Uw = tUw.ref();
+
+    const polyPatch& pp = primaryMesh().boundaryMesh()[patch_.index()];
+
+    if
+    (
+        primaryMesh().moving()
+     && isA<movingWallVelocityFvPatchVectorField>(pp)
+    )
+    {
+        const movingWallVelocityFvPatchVectorField& wpp =
+            refCast<const movingWallVelocityFvPatchVectorField>(pp);
+
+        tmp<vectorField> tUwall = wpp.Uwall();
+
+         // Map Up to area
+        tmp<vectorField> UsWall = vsmPtr_->mapToSurface(tUwall());
+
+        const vectorField& nHat =
+            regionMesh().faceAreaNormals().internalField();
+
+        Uw.primitiveFieldRef() = UsWall() - nHat*(UsWall() & nHat);
+    }
+
+    return tUw;
+}
+
+
+Foam::tmp<Foam::areaVectorField> liquidFilmBase::Us() const
+{
+    tmp<areaVectorField> tUs
+    (
+        new areaVectorField
+        (
+            IOobject
+            (
+                "tUs",
+                primaryMesh().time().timeName(),
+                primaryMesh()
+            ),
+            regionMesh(),
+            dimensionedVector(dimVelocity, Zero)
+        )
+    );
+
+    // Uf quadratic profile
+    tUs.ref() = Foam::sqrt(2.0)*Uf_;
+
+    return tUs;
+}
+
+
+Foam::tmp<Foam::areaVectorField> liquidFilmBase::Up() const
+{
+    const label patchi = patch_.index();
+
+    const volVectorField& Uprimary =
+        primaryMesh().lookupObject<volVectorField>(UName_);
+
+    const fvPatchVectorField& Uw = Uprimary.boundaryField()[patchi];
+
+    tmp<areaVectorField> tUp
+    (
+        new areaVectorField
+        (
+            IOobject
+            (
+                "tUp",
+                primaryMesh().time().timeName(),
+                primaryMesh()
+            ),
+            regionMesh(),
+            dimensionedVector(dimVelocity, Zero)
+        )
+    );
+
+    areaVectorField& Up = tUp.ref();
+/*
+    typedef compressible::turbulenceModel cmpTurbModelType;
+    typedef incompressible::turbulenceModel incmpTurbModelType;
+
+    word turbName
+    (
+        IOobject::groupName
+        (
+            turbulenceModel::propertiesName,
+            Up_.internalField().group()
+        )
+    );
+
+    scalarField nu(patch_.size(), Zero);
+    if (primaryMesh().foundObject<cmpTurbModelType>(turbName))
+    {
+        const cmpTurbModelType& turbModel =
+            primaryMesh().lookupObject<cmpTurbModelType>(turbName);
+
+        //const basicThermo& thermo = turbModel.transport();
+
+        nu = turbModel.nu(patchi);
+    }
+    if (primaryMesh().foundObject<incmpTurbModelType>(turbName))
+    {
+        const incmpTurbModelType& turbModel =
+            primaryMesh().lookupObject<incmpTurbModelType>(turbName);
+
+        nu = turbModel.nu(patchi);
+    }
+*/
+    scalarField hp(patch_.size(), Zero);
+
+    // map areas h to hp on primary
+    vsmPtr_->mapToField(h_, hp);
+
+    const vectorField& nHat = regionMesh().faceAreaNormals().internalField();
+
+    // U tangential on primary
+    const vectorField Ust(-Uw.snGrad()*hp);
+
+    // Map U tang to surface
+    Up.primitiveFieldRef() = vsmPtr_->mapToSurface(Ust);
+
+    // U tangent on surface
+    Up.primitiveFieldRef() -= nHat*(Up.primitiveField() & nHat);
+
+    return tUp;
+}
+
+tmp<areaScalarField> liquidFilmBase::pg() const
+{
+    tmp<areaScalarField> tpg
+    (
+        new areaScalarField
+        (
+            IOobject
+            (
+                "tpg",
+                primaryMesh().time().timeName(),
+                primaryMesh()
+            ),
+            regionMesh(),
+            dimensionedScalar(dimPressure, Zero)
+        )
+    );
+
+    areaScalarField& pfg = tpg.ref();
+
+    if (pName_ != word::null)
+    {
+        const volScalarField& pp =
+            primaryMesh().lookupObject<volScalarField>(pName_);
+
+        const volScalarField::Boundary& pw = pp.boundaryField();
+
+        pfg.primitiveFieldRef() = vsmPtr_->mapInternalToSurface<scalar>(pw)();
+    }
+
+    return tpg;
+}
+
+
+void liquidFilmBase::addSources
+(
+    const label patchi,
+    const label facei,
+    const scalar massSource,
+    const vector& momentumSource,
+    const scalar pressureSource,
+    const scalar energySource
+)
+{
+
+    massSource_.boundaryFieldRef()[patchi][facei] += massSource;
+
+    pnSource_.boundaryFieldRef()[patchi][facei] += pressureSource;
+
+    momentumSource_.boundaryFieldRef()[patchi][facei] += momentumSource;
+
+    addedMassTotal_ += massSource;
+
+    if (debug)
+    {
+        InfoInFunction
+            << "\nSurface film: " << type() << ": adding to film source:" << nl
+            << "    mass     = " << massSource << nl
+            << "    momentum = " << momentumSource << nl
+            << "    pressure = " << pressureSource << endl;
+    }
+}
+
+
+Foam::fa::options& liquidFilmBase::faOptions()
+{
+     return faOptions_;
+}
+
+
+const areaVectorField& liquidFilmBase::Uf() const
+{
+     return Uf_;
+}
+
+
+const areaScalarField& liquidFilmBase::Tf() const
+{
+     return Tf_;
+}
+
+const areaScalarField& liquidFilmBase::gn() const
+{
+     return gn_;
+}
+
+const areaScalarField& liquidFilmBase::h() const
+{
+     return h_;
+}
+
+const dimensionedScalar& liquidFilmBase::h0() const
+{
+     return h0_;
+}
+
+const regionFaModel& liquidFilmBase::region() const
+{
+    return *this;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/liquidFilmBase.H b/src/regionFaModels/liquidFilm/liquidFilmBase.H
new file mode 100644
index 0000000000000000000000000000000000000000..e0f9522fee9e0c45a0181d099160066ba3022a33
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/liquidFilmBase.H
@@ -0,0 +1,302 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::regionModels::liquidFilmBase
+
+Description
+    Base class for thermal 2D shells
+
+SourceFiles
+    liquidFilmBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef liquidFilmBase_H
+#define liquidFilmBase_H
+
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "areaFieldsFwd.H"
+#include "volFieldsFwd.H"
+#include "regionFaModel.H"
+#include "faOptions.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+
+/*---------------------------------------------------------------------------*\
+                         Class liquidFilmBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class liquidFilmBase
+:
+    public regionFaModel
+{
+private:
+
+    // Private Member Functions
+
+        //- No copy construct
+        liquidFilmBase(const liquidFilmBase&) = delete;
+
+        //- No copy assignment
+        void operator=(const liquidFilmBase&) = delete;
+
+        //- Initialize thermal Baffle
+        void init();
+
+
+protected:
+
+    // Protected Data
+
+        // Solution parameters
+
+            //- Momentum predictor
+            Switch momentumPredictor_;
+
+            //- Number of outer correctors
+            label nOuterCorr_;
+
+            //- Number of PISO-like correctors
+            label nCorr_;
+
+            //- Number of non-orthogonal correctors
+            label nNonOrthCorr_;
+
+            //- Cumulative continuity error
+            scalar cumulativeContErr_;
+
+            //- Small delta
+            dimensionedScalar h0_;
+
+        //- Name of the velocity field
+        word UName_;
+
+        //- Name of the pressure field
+        word pName_;
+
+        //- Reference absolute pressure
+        scalar pRef_;
+
+
+        // Fields
+
+            //- Film hight
+            areaScalarField h_;
+
+            //- Film velocity
+            areaVectorField Uf_;
+
+            //- Film pressure
+            areaScalarField pf_;
+
+             //- Primary region pressure
+            areaScalarField ppf_;
+
+            //- Film momentum flux
+            edgeScalarField phif_;
+
+            //- Film height flux
+            edgeScalarField phi2s_;
+
+            //- Film temperature
+            areaScalarField Tf_;
+
+            //- Normal gravity field
+            areaScalarField gn_;
+
+
+            // Mass exchanage fields
+
+                //- Mass
+                volScalarField massSource_;
+
+                //- Momentum
+                volVectorField momentumSource_;
+
+                //- Normal pressure by particles
+                volScalarField pnSource_;
+
+                //- Energy
+                volScalarField energySource_;
+
+                //- Total mass added
+                scalar addedMassTotal_;
+
+
+        //- Pointer to faOptions
+        Foam::fa::options& faOptions_;
+
+
+    // Protected Member Functions
+
+        //- Read control parameters from dictionary
+        virtual bool read(const dictionary&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("liquidFilmBase");
+
+
+    // Declare runtime constructor selection tables
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             liquidFilmBase,
+             dictionary,
+             (
+                const word& modelType,
+                const fvPatch& patch,
+                const dictionary& dict
+             ),
+             (modelType, patch, dict)
+         );
+
+
+    // Constructors
+
+
+        //- Construct from type name and mesh and dict
+        liquidFilmBase
+        (
+            const word& modelType,
+            const fvPatch& patch,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected model using dictionary
+        static autoPtr<liquidFilmBase> New
+        (
+            const fvPatch& patch,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~liquidFilmBase();
+
+
+    // Member Functions
+
+        //- Courant number evaluation
+        virtual scalar CourantNumber() const;
+
+
+        // Helper functions
+
+            //- Wall velocity
+            tmp<areaVectorField> Uw() const;
+
+            //- Film surface film velocity
+            tmp<areaVectorField> Us() const;
+
+            //- Primary region velocity at film hight. Assume the film to be
+            // in the viscous sub-layer
+            tmp<areaVectorField> Up() const;
+
+            //- Map primary static pressure
+            tmp<areaScalarField> pg() const;
+
+        // Access functions
+
+            //- Return faOptions
+            Foam::fa::options& faOptions();
+
+            //- Access const reference Uf
+            const areaVectorField& Uf() const;
+
+            //- Access const reference Tf
+            const areaScalarField& Tf() const;
+
+            //- Access const reference gn
+            const areaScalarField& gn() const;
+
+            //- Access const reference h
+            const areaScalarField& h() const;
+
+            //- Return h0
+            const dimensionedScalar& h0() const;
+
+            //- Access to this region
+            const regionFaModel& region() const;
+
+
+        // Thermo variables
+
+             //- Access const reference mu
+            virtual const areaScalarField& mu() const = 0;
+
+            //- Access const reference rho
+            virtual const areaScalarField& rho() const = 0;
+
+            //- Access const reference sigma
+            virtual const areaScalarField& sigma() const = 0;
+
+
+        // External hook to add sources and mass exchange variables
+
+
+
+            //- Add sources
+            virtual void addSources
+            (
+                const label patchi,            // patchi on primary region
+                const label facei,             // facei of patchi
+                const scalar massSource,       // [kg]
+                const vector& momentumSource,  // [kg.m/s] (tang'l momentum)
+                const scalar pressureSource,   // [kg.m/s] (normal momentum)
+                const scalar energySource = 0  // [J]
+            );
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/liquidFilmBaseNew.C b/src/regionFaModels/liquidFilm/liquidFilmBaseNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..738345c437be65e664641889dc2ef68f7f498043
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/liquidFilmBaseNew.C
@@ -0,0 +1,71 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "liquidFilmBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<liquidFilmBase> liquidFilmBase::New
+(
+    const fvPatch& p,
+    const dictionary& dict
+)
+{
+    word modelType = dict.get<word>("liquidFilmModel");
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalErrorInFunction
+            << "Unknown liquidFilmBase type "
+            << modelType << nl << nl
+            << "Valid liquidFilmBase types :" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<liquidFilmBase>(cstrIter()(modelType, p, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..b4a07ededc9b5ce39950e94998180af53abd977e
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C
@@ -0,0 +1,245 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "liquidFilmModel.H"
+#include "addToRunTimeSelectionTable.H"
+#include "uniformDimensionedFields.H"
+#include "gravityMeshObject.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(liquidFilmModel, 0);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+bool liquidFilmModel::read(const dictionary& dict)
+{
+    liquidFilmBase::read(dict);
+    return true;
+}
+
+void liquidFilmModel::correctThermoFields()
+{
+    scalarField X(thermo_.size(), 1);
+
+    forAll (rho_, faceI)
+    {
+        rho_[faceI] = thermo_.rho(pRef_, Tf_[faceI], X);
+        mu_[faceI] = thermo_.mu(pRef_, Tf_[faceI], X);
+        sigma_[faceI] = thermo_.sigma(pRef_, Tf_[faceI], X);
+    }
+
+    forAll (regionMesh().boundary(), patchI)
+    {
+        const scalarField& patchTf = Tf_.boundaryFieldRef()[patchI];
+
+        scalarField& patchRho = rho_.boundaryFieldRef()[patchI];
+        scalarField& patchmu = mu_.boundaryFieldRef()[patchI];
+        scalarField& patchsigma = sigma_.boundaryFieldRef()[patchI];
+
+        forAll(patchRho, edgeI)
+        {
+            patchRho[edgeI] = thermo_.rho(pRef_, patchTf[edgeI], X);
+            patchmu[edgeI] = thermo_.mu(pRef_, patchTf[edgeI], X);
+            patchsigma[edgeI] = thermo_.sigma(pRef_, patchTf[edgeI], X);
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+liquidFilmModel::liquidFilmModel
+(
+    const word& modelType,
+    const fvPatch& patch,
+    const dictionary& dict
+)
+:
+    liquidFilmBase(modelType, patch, dict),
+
+    thermo_(dict.subDict("thermo")),
+
+    rho_
+    (
+        IOobject
+        (
+            "rhof",
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimDensity, Zero)
+    ),
+    mu_
+    (
+        IOobject
+        (
+            "muf",
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimViscosity, Zero)
+    ),
+    sigma_
+    (
+        IOobject
+        (
+            "sigmaf",
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimMass/sqr(dimTime), Zero)
+    ),
+    hRho_
+    (
+        IOobject
+        (
+            h_.name() + "*" + rho_.name(),
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(h_.dimensions()*rho_.dimensions(), Zero)
+    ),
+    rhoSp_
+    (
+        IOobject
+        (
+            "rhoSp",
+            primaryMesh().time().timeName(),
+            primaryMesh()
+        ),
+        regionMesh(),
+        dimensionedScalar(dimVelocity, Zero)
+    ),
+    USp_
+    (
+        IOobject
+        (
+            "USp",
+            primaryMesh().time().timeName(),
+            primaryMesh()
+        ),
+        regionMesh(),
+        dimensionedVector(sqr(dimVelocity), Zero)
+    ),
+    pnSp_
+    (
+        IOobject
+        (
+            "pnSp",
+            primaryMesh().time().timeName(),
+            primaryMesh()
+        ),
+        regionMesh(),
+        dimensionedScalar(dimPressure, Zero)
+    ),
+    turbulence_(filmTurbulenceModel::New(*this, dict))
+{
+    correctThermoFields();
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+liquidFilmModel::~liquidFilmModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const areaScalarField& liquidFilmModel::mu() const
+{
+     return mu_;
+}
+
+
+const areaScalarField& liquidFilmModel::rho() const
+{
+     return rho_;
+}
+
+
+const areaScalarField& liquidFilmModel::sigma() const
+{
+     return sigma_;
+}
+
+void liquidFilmModel::preEvolveRegion()
+{
+    const scalar deltaT = primaryMesh().time().deltaTValue();
+    const scalarField rAreaDeltaT = 1/deltaT/regionMesh().S().field();
+
+    // Map the total mass, mom and pnSource from particles
+    rhoSp_.primitiveFieldRef() =
+        vsm().mapToSurface(massSource_.boundaryField()[patchID()]);
+    // [kg.m/s]
+    USp_.primitiveFieldRef() =
+        vsm().mapToSurface(momentumSource_.boundaryField()[patchID()]);
+    pnSp_.primitiveFieldRef() =
+        vsm().mapToSurface(pnSource_.boundaryField()[patchID()]);
+
+
+    // Calculate rate per area
+    rhoSp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
+    USp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
+    pnSp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
+}
+
+void liquidFilmModel::postEvolveRegion()
+{
+    massSource_ == dimensionedScalar(massSource_.dimensions(), Zero);
+    momentumSource_ == dimensionedVector(momentumSource_.dimensions(), Zero);
+    pnSource_ == dimensionedScalar(pnSource_.dimensions(), Zero);
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..0aeebf61189b9bf4a6dcc11da2e243203f14879f
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H
@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::regionFaModels::liquidFilmModel
+
+Description
+    Thin Model Film model.
+
+
+SourceFiles
+    liquidFilmModel.C
+    kinematicThinFilmI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef liquidFilmModel_H
+#define liquidFilmModel_H
+
+#include "volFieldsFwd.H"
+#include "liquidFilmBase.H"
+#include "faMesh.H"
+#include "filmTurbulenceModel.H"
+#include "liquidMixtureProperties.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class liquidFilmModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class liquidFilmModel
+:
+    public liquidFilmBase
+{
+private:
+
+    // Private member functions
+
+        //- No copy construct
+        liquidFilmModel(const liquidFilmModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const liquidFilmModel&) = delete;
+
+        //- Initialize liquidFilmModel
+        void init();
+
+
+protected:
+
+
+    // Thermo properties
+
+
+            //- Liquid thermo
+            liquidMixtureProperties thermo_;
+
+
+            // Fields
+
+                //- Density [kg/m3]
+                areaScalarField rho_;
+
+                //- Dynamic viscosity [Pa.s]
+                areaScalarField mu_;
+
+                //- Surface tension [m/s2]
+                areaScalarField sigma_;
+
+                //- Film rho*height
+                areaScalarField hRho_;
+
+                // Mass exchange sources
+
+                    //- Mass source
+                    areaScalarField rhoSp_;
+
+                    //- Momentum source
+                    areaVectorField USp_;
+
+                    //- Normal pressure by particles
+                    areaScalarField pnSp_;
+
+
+    // General properties
+
+
+
+            //- Turbulence model
+            autoPtr<filmTurbulenceModel> turbulence_;
+
+
+    // Protected member functions
+
+
+        //- Read control parameters from dictionary
+        virtual bool read(const dictionary& dict);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("liquidFilmModel");
+
+
+    // Constructors
+
+
+        //- Construct from components and dict
+        liquidFilmModel
+        (
+            const word& modelType,
+            const fvPatch& patch,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~liquidFilmModel();
+
+
+    // Member Functions
+
+        // Help function
+
+
+            //- Correct thermo
+            void correctThermoFields();
+
+
+        // Access functions
+
+            //- Access const reference mu
+            const areaScalarField& mu() const;
+
+            //- Access const reference rho
+            const areaScalarField& rho() const;
+
+            //- Access const reference sigma
+            const areaScalarField& sigma() const;
+
+
+        // Evolution
+
+            //- Pre-evolve film
+            virtual void preEvolveRegion();
+
+            //- Evolve the film
+            //virtual void evolveRegion();
+
+            //- Post-evolve film
+            virtual void postEvolveRegion();
+
+
+       // I-O
+
+            //- Provide some feedback
+            //virtual void info();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..6c0d5a2731ddad9d78cfa7e147bcc6ac7a9788bb
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "filmTurbulenceModel.H"
+#include "gravityMeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(filmTurbulenceModel, 0);
+defineRunTimeSelectionTable(filmTurbulenceModel, dictionary);
+
+const Enum
+<
+    filmTurbulenceModel::frictionMethodType
+>
+filmTurbulenceModel::frictionMethodTypeNames_
+{
+    { frictionMethodType::mquadraticProfile, "quadraticProfile" },
+    { frictionMethodType::mlinearProfile, "linearProfile" },
+    { frictionMethodType::mDarcyWeisbach, "DarcyWeisbach" },
+    { frictionMethodType::mManningStrickler, "ManningStrickler" }
+};
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+filmTurbulenceModel::filmTurbulenceModel
+(
+    const word& modelType,
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    film_(film),
+    dict_(dict.subDict(modelType + "Coeffs")),
+    method_(frictionMethodTypeNames_.get("friction", dict_))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+filmTurbulenceModel::~filmTurbulenceModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const liquidFilmBase& filmTurbulenceModel::film() const
+{
+    return film_;
+}
+
+
+tmp<areaScalarField> filmTurbulenceModel::Cw() const
+{
+    auto tCw =
+        tmp<areaScalarField>::New
+        (
+            IOobject
+            (
+                "tCw",
+                film_.primaryMesh().time().timeName(),
+                film_.primaryMesh()
+            ),
+            film_.regionMesh(),
+            dimensionedScalar(dimVelocity)
+        );
+    auto& Cw = tCw.ref();
+
+    switch (method_)
+    {
+        case mquadraticProfile:
+        {
+            const scalarField& mu = film_.mu().primitiveField();
+            const scalarField& h = film_.h().primitiveField();
+            const scalar h0 = film_.h0().value();
+            Cw.primitiveFieldRef() = mu/((1.0/3.0)*(h + h0));
+            Cw.min(5000.0);
+            break;
+        }
+        case mlinearProfile:
+        {
+            const scalarField& mu = film_.mu().primitiveField();
+            const scalarField& h = film_.h().primitiveField();
+            const scalar h0 = film_.h0().value();
+            Cw.primitiveFieldRef() = 2*mu/(h + h0);
+            break;
+        }
+        case mDarcyWeisbach:
+        {
+            const uniformDimensionedVectorField& g =
+                meshObjects::gravity::New(film_.primaryMesh().time());
+
+            const vectorField& Uf = film_.Uf().primitiveField();
+
+            scalar Cf = dict_.get<scalar>("Cf");
+            Cw.primitiveFieldRef() = Cf*mag(g.value())*mag(Uf);
+            break;
+        }
+        case mManningStrickler:
+        {
+            const uniformDimensionedVectorField& g =
+                meshObjects::gravity::New(film_.primaryMesh().time());
+
+            const scalar n = dict_.get<scalar>("n");
+
+            const vectorField& Uf = film_.Uf().primitiveField();
+            const scalarField& h = film_.h().primitiveField();
+            const scalar h0 = film_.h0().value();
+
+            Cw.primitiveFieldRef() =
+                sqr(n)*mag(g.value())*mag(Uf)/cbrt(h+ h0);
+            break;
+        }
+        default:
+        {
+            FatalErrorInFunction
+                << "Unimplemented method "
+                << frictionMethodTypeNames_[method_] << nl
+                << "Please set 'frictionMethod' to one of "
+                << flatOutput(frictionMethodTypeNames_.sortedToc()) << nl
+                << exit(FatalError);
+
+            break;
+        }
+    }
+
+    return tCw;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..ce2b620af9cb7e1a46055159e220b1ce8b37b147
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
@@ -0,0 +1,183 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::regionModels::areaSurfaceFilmModels::filmTurbulenceModel
+
+Description
+    Base class for film turbulence models
+
+SourceFiles
+    filmTurbulenceModel.C
+    filmTurbulenceModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filmTurbulenceModel_H
+#define filmTurbulenceModel_H
+
+#include "areaFieldsFwd.H"
+#include "runTimeSelectionTables.H"
+#include "faMatrices.H"
+#include "liquidFilmBase.H"
+#include "fvMesh.H"
+#include "Enum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class filmTurbulenceModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class filmTurbulenceModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        filmTurbulenceModel(const filmTurbulenceModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const filmTurbulenceModel&) = delete;
+
+
+public:
+
+    // Public Enumerations
+
+        //- Options for the friction models
+        enum frictionMethodType
+        {
+            mquadraticProfile,
+            mlinearProfile,
+            mDarcyWeisbach,
+            mManningStrickler
+        };
+
+
+protected:
+
+    // Protected Data
+
+        //- Reference to liquidFilmBase
+        const liquidFilmBase& film_;
+
+        //- Names for friction models
+        static const Enum<frictionMethodType> frictionMethodTypeNames_;
+
+        //- Model dictionary
+        const dictionary dict_;
+
+        //- Method used
+        const frictionMethodType method_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("filmTurbulenceModel");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             filmTurbulenceModel,
+             dictionary,
+             (
+                liquidFilmBase& film,
+                const dictionary& dict
+             ),
+             (film, dict)
+         );
+
+
+    // Constructors
+
+        //- Construct from type name, dictionary and surface film model
+        filmTurbulenceModel
+        (
+            const word& modelType,
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected injection model
+        static autoPtr<filmTurbulenceModel> New
+        (
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~filmTurbulenceModel();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return film
+            const liquidFilmBase& film() const;
+
+
+        // Evolution
+
+            //- Return the wall film surface friction
+            virtual tmp<areaScalarField> Cw() const;
+
+            //- Return the film turbulence viscosity
+            virtual tmp<areaScalarField> mut() const = 0;
+
+            //- Correct/update the model
+            virtual void correct() = 0;
+
+            //- Return the source for the film momentum equation
+            virtual tmp<faVectorMatrix> Su(areaVectorField& U) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..446e59250abddfd41a49ab8b2d35f4d867ab9f7b
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "filmTurbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<filmTurbulenceModel> filmTurbulenceModel::New
+(
+    liquidFilmBase& model,
+    const dictionary& dict
+)
+{
+    const word modelType(dict.get<word>("turbulence"));
+
+    Info<< "    Selecting filmTurbulenceModel " << modelType << endl;
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "filmTurbulenceModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return autoPtr<filmTurbulenceModel>(cstrIter()(model, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
new file mode 100644
index 0000000000000000000000000000000000000000..472f80dee221818236005869bb6217960f17ea4a
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "laminar.H"
+#include "addToRunTimeSelectionTable.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(laminar, 0);
+addToRunTimeSelectionTable(filmTurbulenceModel, laminar, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+laminar::laminar
+(
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    filmTurbulenceModel(type(), film, dict),
+    Cf_(dict_.get<scalar>("Cf"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+laminar::~laminar()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+tmp<areaScalarField> laminar::mut() const
+{
+    auto tmut =
+        tmp<areaScalarField>::New
+        (
+            IOobject
+            (
+                "mut",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            film().regionMesh(),
+            dimensionedScalar(dimMass/dimLength/dimTime)
+        );
+
+    return tmut;
+}
+
+
+void laminar::correct()
+{}
+
+
+tmp<faVectorMatrix> laminar::Su(areaVectorField& U) const
+{
+    // local references to film fields
+    tmp<areaVectorField> Uw = film_.Uw();
+    tmp<areaVectorField> Up = film_.Up();
+
+    // employ simple coeff-based model
+    areaScalarField Cs("Cs", Cf_*mag(Up() - U));
+
+    tmp<areaScalarField> wf = Cw();
+
+    return
+    (
+       - fam::Sp(Cs, U) + Cs*Up()     // surface contribution
+       - fam::Sp(wf(), U) + wf()*Uw() // wall contribution
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
new file mode 100644
index 0000000000000000000000000000000000000000..a405929e2343f555a0eb7c2e24583f9d5a0934b1
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::regionModels::areaSurfaceFilmModels::laminar
+
+Description
+    Film laminar turbulence model.
+
+SourceFiles
+    laminar.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef laminar_H
+#define laminar_H
+
+#include "filmTurbulenceModel.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class laminar Declaration
+\*---------------------------------------------------------------------------*/
+
+class laminar
+:
+    public filmTurbulenceModel
+{
+    // Private Data
+
+        //- Surface roughness coefficient
+        scalar Cf_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        laminar(const laminar&) = delete;
+
+        //- No copy assignment
+        void operator=(const laminar&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("laminar");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        laminar(liquidFilmBase& film, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~laminar();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Return the film turbulence viscosity
+            virtual tmp<areaScalarField> mut() const;
+
+            //- Correct/update the model
+            virtual void correct();
+
+            //- Return the source for the film momentum equation
+            virtual tmp<faVectorMatrix> Su(areaVectorField& U) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //