diff --git a/src/Allwmake b/src/Allwmake
index e488f93978..0af702324e 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -83,14 +83,18 @@ wmake $targetType ODE
 thermophysicalModels/Allwmake $targetType $*
 TurbulenceModels/Allwmake $targetType $*
 wmake $targetType combustionModels
+
+wmakeLnInclude -u regionFaModels
+wmakeLnInclude -u faOptions
 regionModels/Allwmake $targetType $*
-lagrangian/Allwmake $targetType $*
 
 wmake $targetType fvOptions
 wmake $targetType faOptions
 wmake $targetType fvMotionSolver
 wmake $targetType regionFaModels
 
+lagrangian/Allwmake $targetType $*
+
 wmake $targetType overset
 
 # snappyHexMesh uses overset voxelMesh
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
index 923af2c002..94e4961ca8 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
index 4696ac7021..d9e474202b 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.C
@@ -41,13 +41,16 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel(CloudType& owner)
     CloudSubModelBase<CloudType>(owner),
     g_(owner.g()),
     ejectedParcelType_(0),
+    injectionOffset_(1.1),
+    minDiameter_(0),
     massParcelPatch_(0),
     diameterParcelPatch_(0),
     UFilmPatch_(0),
     rhoFilmPatch_(0),
     deltaFilmPatch_(0),
     nParcelsTransferred_(0),
-    nParcelsInjected_(0)
+    nParcelsInjected_(0),
+    totalMassTransferred_(0)
 {}
 
 
@@ -63,7 +66,15 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel
     g_(owner.g()),
     ejectedParcelType_
     (
-        this->coeffDict().getOrDefault("ejectedParcelType", -1)
+        this->coeffDict().template getOrDefault<label>("ejectedParcelType", -1)
+    ),
+    injectionOffset_
+    (
+        this->coeffDict().template getOrDefault<scalar>("injectionOffset", 1.1)
+    ),
+    minDiameter_
+    (
+        this->coeffDict().template getOrDefault<scalar>("minDiameter", -1)
     ),
     massParcelPatch_(0),
     diameterParcelPatch_(0),
@@ -71,7 +82,8 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel
     rhoFilmPatch_(0),
     deltaFilmPatch_(owner.mesh().boundary().size()),
     nParcelsTransferred_(0),
-    nParcelsInjected_(0)
+    nParcelsInjected_(0),
+    totalMassTransferred_()
 {}
 
 
@@ -84,13 +96,16 @@ Foam::SurfaceFilmModel<CloudType>::SurfaceFilmModel
     CloudSubModelBase<CloudType>(sfm),
     g_(sfm.g_),
     ejectedParcelType_(sfm.ejectedParcelType_),
+    injectionOffset_(sfm.injectionOffset_),
+    minDiameter_(sfm.minDiameter_),
     massParcelPatch_(sfm.massParcelPatch_),
     diameterParcelPatch_(sfm.diameterParcelPatch_),
     UFilmPatch_(sfm.UFilmPatch_),
     rhoFilmPatch_(sfm.rhoFilmPatch_),
     deltaFilmPatch_(sfm.deltaFilmPatch_),
     nParcelsTransferred_(sfm.nParcelsTransferred_),
-    nParcelsInjected_(sfm.nParcelsInjected_)
+    nParcelsInjected_(sfm.nParcelsInjected_),
+    totalMassTransferred_(sfm.totalMassTransferred_)
 {}
 
 
@@ -130,7 +145,7 @@ void Foam::SurfaceFilmModel<CloudType>::injectParticles
                     diameterParcelPatch_[j],
                     deltaFilmPatch_[primaryPatchi][j]
                 );
-            const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j];
+            const point pos = Cf[j] - injectionOffset_*offset*Sf[j]/magSf[j];
 
             // Create a new parcel
             parcelType* pPtr =
@@ -144,7 +159,6 @@ void Foam::SurfaceFilmModel<CloudType>::injectParticles
             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
@@ -215,7 +229,7 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackCloudType& cloud)
             >(names[i]);
 
         // Check that it is a type areaFilm
-        if (regionFa && isA<areaFilm>(*regionFa))
+        if (regionFa && isA<areaFilm>(*regionFa) && regionFa->active())
         {
             areaFilm& film =
                 const_cast<areaFilm&>(refCast<const areaFilm>(*regionFa));
@@ -227,6 +241,22 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackCloudType& cloud)
             cacheFilmFields(patchId, film);
 
             injectParticles(patchId, injectorCellsPatch, cloud);
+
+            forAll(injectorCellsPatch, facei)
+            {
+                if (diameterParcelPatch_[facei] > 0)
+                {
+                    film.addSources
+                    (
+                        patchId,
+                        facei,
+                       -massParcelPatch_[facei],    // mass
+                        Zero,                       // tangential momentum
+                        Zero,                       // impingement
+                        Zero                        // energy
+                    );
+                }
+            }
         }
     }
 }
@@ -266,22 +296,26 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
     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);
+    massParcelPatch_.setSize(filmModel.Uf().size(), Zero);
+
+    const scalarField& massParcelPatch =
+        filmModel.cloudMassTrans().boundaryField()[filmPatchi];
+
+    map.mapToField(massParcelPatch, massParcelPatch_);
+
+    const scalarField&  diameterParcelPatch =
+        filmModel.cloudDiameterTrans().boundaryField()[filmPatchi];
+
+    diameterParcelPatch_.setSize(filmModel.Uf().size(), Zero);
+    map.mapToField(diameterParcelPatch, diameterParcelPatch_);
+
+    UFilmPatch_.setSize(filmModel.Uf().size(), vector::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]);
@@ -303,6 +337,14 @@ void Foam::SurfaceFilmModel<CloudType>::setParcelProperties
 
     p.nParticle() = massParcelPatch_[filmFacei]/p.rho()/vol;
 
+    if (minDiameter_ != -1)
+    {
+        if (p.d() < minDiameter_)
+        {
+            p.nParticle() = 0;
+        }
+    }
+
     if (ejectedParcelType_ >= 0)
     {
         p.typeId() = ejectedParcelType_;
@@ -319,22 +361,33 @@ void Foam::SurfaceFilmModel<CloudType>::info(Ostream& os)
     label nInject0 =
         this->template getModelProperty<label>("nParcelsInjected");
 
+    scalar massTransferred0 =
+        this->template getModelProperty<scalar>("massTransferred");
+
     label nTransTotal =
         nTrans0 + returnReduce(nParcelsTransferred_, sumOp<label>());
 
     label nInjectTotal =
         nInject0 + returnReduce(nParcelsInjected_, sumOp<label>());
 
+    scalar massTransferredTotal =
+        massTransferred0 + returnReduce(totalMassTransferred_, sumOp<scalar>());
+
+
     os  << "    Surface film:" << nl
         << "      - parcels absorbed            = " << nTransTotal << nl
+        << "      - mass absorbed               = " << massTransferredTotal << nl
         << "      - parcels ejected             = " << nInjectTotal << endl;
 
     if (this->writeTime())
     {
         this->setModelProperty("nParcelsTransferred", nTransTotal);
         this->setModelProperty("nParcelsInjected", nInjectTotal);
+        this->setModelProperty("massTransferred", massTransferredTotal);
+
         nParcelsTransferred_ = 0;
         nParcelsInjected_ = 0;
+        totalMassTransferred_ = 0;
     }
 }
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
index 0e5a7a4cef..410dd63e23 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModel.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +47,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "CloudSubModelBase.H"
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -96,14 +98,20 @@ protected:
         //  type
         label ejectedParcelType_;
 
+        //- Injection offset position
+        scalar injectionOffset_;
+
+        //- Minimum diameter particle injection
+        scalar minDiameter_;
+
 
         // Cached injector fields per film patch
 
             //- Parcel mass / patch face
-            scalarList massParcelPatch_;
+            scalarField massParcelPatch_;
 
             //- Parcel diameter / patch face
-            scalarList diameterParcelPatch_;
+            scalarField diameterParcelPatch_;
 
             //- Film velocity / patch face
             Field<vector> UFilmPatch_;
@@ -124,6 +132,12 @@ protected:
             label nParcelsInjected_;
 
 
+        // Total mass info
+
+            //- Total mass transferred to the film
+            scalar totalMassTransferred_;
+
+
     // Protected functions
 
         //- Cache the film fields in preparation for injection
@@ -138,7 +152,7 @@ protected:
         virtual void cacheFilmFields
         (
             const label primaryPatchi,
-            const regionModels::areaSurfaceFilmModels::liquidFilmBase&
+            const areaFilm&
         );
 
         //- Inject particles in cloud
@@ -232,6 +246,12 @@ public:
             //  the film model
             inline label& nParcelsInjected();
 
+            //- Return non-const total mass transferred
+            inline scalar& totalMassTransferred();
+
+            //- Return consr access to mass transferred
+            inline scalar totalMassTransferred() const;
+
 
         // Member Functions
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelI.H b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelI.H
index 67420350c6..4240d00ea3 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelI.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/SurfaceFilmModel/SurfaceFilmModel/SurfaceFilmModelI.H
@@ -63,5 +63,18 @@ Foam::label Foam::SurfaceFilmModel<CloudType>::nParcelsInjected() const
     return nParcelsInjected_;
 }
 
+template<class CloudType>
+Foam::scalar& Foam::SurfaceFilmModel<CloudType>::totalMassTransferred()
+{
+    return totalMassTransferred_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::SurfaceFilmModel<CloudType>::totalMassTransferred() const
+{
+    return totalMassTransferred_;
+}
+
 
 // ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
index 70594de693..bfe1f3fa9b 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C
@@ -131,57 +131,51 @@ Foam::vector Foam::ThermoSurfaceFilm<CloudType>::splashDirection
 
 
 template<class CloudType>
-void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
-(
-    const parcelType& p,
-    const polyPatch& pp,
-    const label facei,
-    const scalar mass,
-    bool& keepParticle
-)
+void Foam::ThermoSurfaceFilm<CloudType>::init()
 {
-    if (debug)
-    {
-        Info<< "Parcel " << p.origId() << " absorbInteraction" << endl;
-    }
+    const fvMesh& mesh = this->owner().mesh();
 
-    // 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;
-
-    if (filmModel_)
-    {
-        filmModel_->addSources
+    // set up filmModel pointer
+    filmModel_ =
+        const_cast<regionFilm*>
         (
-            pp.index(),
-            facei,
-            mass,                           // mass
-            mass*Ut,                        // tangential momentum
-            mass*mag(Un),                   // impingement pressure
-            mass*p.hs()                     // energy
+            mesh.time().objectRegistry::template findObject
+            <
+                regionFilm
+            >
+            (
+                "surfaceFilmProperties"
+            )
         );
 
-        this->nParcelsTransferred()++;
+    // set up areaFilms
+    const wordList names =
+        mesh.time().objectRegistry::template
+            sortedNames<regionModels::regionFaModel>();
 
-        keepParticle = false;
+    forAll(names, i)
+    {
+        const regionModels::regionFaModel* regionFa =
+            mesh.time().objectRegistry::template findObject
+            <
+                regionModels::regionFaModel
+            >(names[i]);
+
+        if (regionFa && isA<areaFilm>(*regionFa))
+        {
+            areaFilm& film =
+                const_cast<areaFilm&>(refCast<const areaFilm>(*regionFa));
+            areaFilms_.append(&film);
+        }
     }
 }
 
+
 template<class CloudType>
+template<class filmType>
 void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
 (
-    areaFilm& filmModel,
+    filmType& film,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -189,7 +183,6 @@ void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
     bool& keepParticle
 )
 {
-
     if (debug)
     {
         Info<< "Parcel " << p.origId() << " absorbInteraction" << endl;
@@ -209,21 +202,22 @@ void Foam::ThermoSurfaceFilm<CloudType>::absorbInteraction
 
     // Parcel tangential velocity
     const vector Ut = Urel - Un;
-DebugVar(mass)
-    filmModel.addSources
+
+    film.addSources
     (
         pp.index(),
         facei,
-        mass,                           // mass
-        mass*Ut,                        // tangential momentum
-        mass*mag(Un),                   // impingement pressure
-        mass*p.hs()                     // energy
+        mass,                             // mass
+        mass*Ut,                          // tangential momentum
+        mass*mag(Un),                     // impingement pressure
+        mass*p.hs()                       // energy
     );
 
     this->nParcelsTransferred()++;
 
-    keepParticle = false;
+    this->totalMassTransferred() += mass;
 
+    keepParticle = false;
 }
 
 
@@ -258,9 +252,10 @@ void Foam::ThermoSurfaceFilm<CloudType>::bounceInteraction
 
 
 template<class CloudType>
+template<class filmType>
 void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
 (
-    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    filmType& filmModel,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -301,22 +296,24 @@ void Foam::ThermoSurfaceFilm<CloudType>::drySplashInteraction
 
     if (We < Wec) // Adhesion - assume absorb
     {
-        absorbInteraction(p, pp, facei, m, keepParticle);
+        absorbInteraction<filmType>
+            (filmModel, 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
-            (p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
+        splashInteraction<filmType>
+            (filmModel, p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
     }
 }
 
 
 template<class CloudType>
+template<class filmType>
 void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
 (
-    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    filmType& filmModel,
     parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -359,7 +356,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
 
     if (We < 2) // Adhesion - assume absorb
     {
-        absorbInteraction(p, pp, facei, m, keepParticle);
+        absorbInteraction<filmType>
+            (filmModel, p, pp, facei, m, keepParticle);
     }
     else if ((We >= 2) && (We < 20)) // Bounce
     {
@@ -377,23 +375,25 @@ void Foam::ThermoSurfaceFilm<CloudType>::wetSplashInteraction
     }
     else if ((We >= 20) && (We < Wec)) // Spread - assume absorb
     {
-        absorbInteraction(p, pp, facei, m, keepParticle);
+        absorbInteraction<filmType>
+            (filmModel, p, pp, facei, m, keepParticle);
     }
     else    // Splash
     {
         // Ratio of incident mass to splashing mass
         // splash mass can be > incident mass due to film entrainment
         const scalar mRatio = 0.2 + 0.9*rndGen_.sample01<scalar>();
-        splashInteraction
-            (p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
+        splashInteraction<filmType>
+            (filmModel, p, pp, facei, mRatio, We, Wec, sigma, keepParticle);
     }
 }
 
 
 template<class CloudType>
+template<class filmType>
 void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
 (
-    //regionModels::surfaceFilmModels::surfaceFilmRegionModel& filmModel,
+    filmType& filmModel,
     const parcelType& p,
     const polyPatch& pp,
     const label facei,
@@ -466,7 +466,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     // Switch to absorb if insufficient energy for splash
     if (EKs <= 0)
     {
-        absorbInteraction(p, pp, facei, m, keepParticle);
+        absorbInteraction<filmType>
+            (filmModel, p, pp, facei, m, keepParticle);
         return;
     }
 
@@ -521,7 +522,8 @@ 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(p, pp, facei, mDash, keepParticle);
+    absorbInteraction<filmType>
+        (filmModel, p, pp, facei, mDash, keepParticle);
 }
 
 
@@ -541,6 +543,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
         owner.db().objectRegistry::template lookupObject<SLGThermo>("SLGThermo")
     ),
     filmModel_(nullptr),
+    areaFilms_(0),
     TFilmPatch_(0),
     CpFilmPatch_(0),
     interactionType_
@@ -570,17 +573,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
         this->coeffDict().readEntry("Cf", Cf_);
     }
 
-    filmModel_ =
-        const_cast<regionModels::surfaceFilmModels::surfaceFilmRegionModel*>
-        (
-            this->owner().mesh().time().objectRegistry::template findObject
-            <
-                regionModels::surfaceFilmModels::surfaceFilmRegionModel
-            >
-            (
-                "surfaceFilmProperties"
-            )
-        );
+    init();
 }
 
 
@@ -594,6 +587,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
     rndGen_(sfm.rndGen_),
     thermo_(sfm.thermo_),
     filmModel_(nullptr),
+    areaFilms_(0),
     TFilmPatch_(sfm.TFilmPatch_),
     CpFilmPatch_(sfm.CpFilmPatch_),
     interactionType_(sfm.interactionType_),
@@ -605,17 +599,7 @@ Foam::ThermoSurfaceFilm<CloudType>::ThermoSurfaceFilm
     Cf_(sfm.Cf_),
     nParcelsSplashed_(sfm.nParcelsSplashed_)
 {
-     filmModel_ =
-        const_cast<regionModels::surfaceFilmModels::surfaceFilmRegionModel*>
-        (
-            this->owner().mesh().time().objectRegistry::template findObject
-            <
-                regionModels::surfaceFilmModels::surfaceFilmRegionModel
-            >
-            (
-                "surfaceFilmProperties"
-            )
-        );
+    init();
 }
 
 
@@ -636,10 +620,10 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
     bool& keepParticle
 )
 {
-    const fvMesh& mesh = this->owner().mesh();
-
     const label patchi = pp.index();
 
+    bool bInteraction(false);
+
     // Check the singleLayer film models
     if (filmModel_ && filmModel_->isRegionPatch(patchi))
     {
@@ -656,7 +640,9 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
             case itAbsorb:
             {
                 const scalar m = p.nParticle()*p.mass();
-                absorbInteraction(p, pp, facei, m, keepParticle);
+
+                absorbInteraction<regionFilm>
+                    (*filmModel_, p, pp, facei, m, keepParticle);
 
                 break;
             }
@@ -666,11 +652,13 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
 
                 if (dry)
                 {
-                    drySplashInteraction(p, pp, facei, keepParticle);
+                    drySplashInteraction<regionFilm>
+                        (*filmModel_, p, pp, facei, keepParticle);
                 }
                 else
                 {
-                    wetSplashInteraction( p, pp, facei, keepParticle);
+                    wetSplashInteraction<regionFilm>
+                        (*filmModel_, p, pp, facei, keepParticle);
                 }
 
                 break;
@@ -684,59 +672,66 @@ bool Foam::ThermoSurfaceFilm<CloudType>::transferParcel
         }
 
         // Transfer parcel/parcel interactions complete
-        return true;
+        bInteraction = true;
     }
-    else
+
+
+    for (areaFilm& film : areaFilms_)
     {
-
-        // Check the finite area films
-        wordList names =
-            mesh.time().objectRegistry::template
-                sortedNames<regionModels::regionFaModel>();
-
-        forAll(names, i)
+        if (patchi == film.patchID())
         {
-            const regionModels::regionFaModel* regionFa =
-                mesh.time().objectRegistry::template findObject
-                <
-                    regionModels::regionFaModel
-                >(names[i]);
+            const label facei = pp.whichFace(p.face());
 
-            // Check that regionFa is a areaFilm
-            if (regionFa && isA<areaFilm>(*regionFa))
+            switch (interactionType_)
             {
-                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:
                 {
-                    // 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);
-                    }
+                    const scalar m = p.nParticle()*p.mass();
+
+                    absorbInteraction<areaFilm>
+                    (
+                        film, p, pp, facei, m, keepParticle
+                    );
+                    break;
+                }
+                case itBounce:
+                {
+                    bounceInteraction(p, pp, facei, keepParticle);
+
+                    break;
+                }
+                case itSplashBai:
+                {
+                    bool dry = film.h()[facei] < deltaWet_;
+
+                    if (dry)
+                    {
+                        drySplashInteraction<areaFilm>
+                            (film, p, pp, facei, keepParticle);
+                    }
+                    else
+                    {
+                        wetSplashInteraction<areaFilm>
+                            (film, p, pp, facei, keepParticle);
+                    }
+
+                    break;
+                }
+                default:
+                {
+                    FatalErrorInFunction
+                        << "Unknown interaction type enumeration"
+                        << abort(FatalError);
                 }
-                // Transfer parcel/parcel interactions complete
-                return true;
             }
+            // Transfer parcel/parcel interactions complete
+            bInteraction = true;
         }
     }
 
     // Parcel not interacting with film
-    return false;
+    return bInteraction;
 }
 
 
@@ -779,7 +774,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::cacheFilmFields
 
     TFilmPatch_.setSize(filmModel.Tf().size(), Zero);
     map.mapToField(filmModel.Tf(), TFilmPatch_);
-    //map.mapToField(filmModel.Cp(), CpFilmPatch_);
+    CpFilmPatch_.setSize(filmModel.Tf().size(), Zero);
+    map.mapToField(filmModel.Cp(), CpFilmPatch_);
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
index b1b821a294..5e5221062b 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.H
@@ -55,6 +55,7 @@ SourceFiles
 #define ThermoSurfaceFilm_H
 
 #include "SurfaceFilmModel.H"
+#include "UPtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -121,6 +122,9 @@ protected:
         typedef typename
             regionModels::areaSurfaceFilmModels::liquidFilmBase areaFilm;
 
+        typedef typename
+            regionModels::surfaceFilmModels::surfaceFilmRegionModel regionFilm;
+
         //- Reference to the cloud random number generator
         Random& rndGen_;
 
@@ -128,8 +132,10 @@ protected:
         const SLGThermo& thermo_;
 
         //- Pointer to filmModel
-        regionModels::surfaceFilmModels::surfaceFilmRegionModel* filmModel_;
+        regionFilm* filmModel_;
 
+        //- UPointers to area films
+        UPtrList<areaFilm> areaFilms_;
 
         // Cached injector fields per film patch
 
@@ -189,24 +195,17 @@ protected:
             const vector& nf
         ) const;
 
+        //- Init films pointers
+        void init();
+
 
         // Interaction models
 
             //- Absorb parcel into film
+            template<class filmType>
             void absorbInteraction
             (
-                //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&,
+                filmType&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -224,9 +223,10 @@ protected:
             ) const;
 
             //- Parcel interaction with dry surface
+            template<class filmType>
             void drySplashInteraction
             (
-                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                filmType&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -234,9 +234,10 @@ protected:
             );
 
             //- Parcel interaction with wetted surface
+            template<class filmType>
             void wetSplashInteraction
             (
-                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                filmType&,
                 parcelType& p,
                 const polyPatch& pp,
                 const label facei,
@@ -244,9 +245,10 @@ protected:
             );
 
             //- Bai parcel splash interaction model
+            template<class filmType>
             void splashInteraction
             (
-                //regionModels::surfaceFilmModels::surfaceFilmRegionModel&,
+                filmType&,
                 const parcelType& p,
                 const polyPatch& pp,
                 const label facei,
diff --git a/src/regionFaModels/Make/files b/src/regionFaModels/Make/files
index 70521f391a..85634c54ee 100644
--- a/src/regionFaModels/Make/files
+++ b/src/regionFaModels/Make/files
@@ -15,9 +15,25 @@ derivedFvPatchFields/vibrationShell/vibrationShellFvPatchScalarField.C
 
 liquidFilm/liquidFilmBase.C
 liquidFilm/liquidFilmBaseNew.C
+
+/* Sub-Model */
+
 liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
 liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
 liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
+
+liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
+
+liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
+liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
+
+liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
+
+liquidFilm/subModels/filmSubModelBase.C
+
+liquidFilm/liquidFilmBase.C
+liquidFilm/liquidFilmBaseNew.C
+
 liquidFilm/liquidFilmModel/liquidFilmModel.C
 liquidFilm/kinematicThinFilm/kinematicThinFilm.C
 derivedFvPatchFields/filmShell/velocityFilmShellFvPatchVectorField.C
diff --git a/src/regionFaModels/liquidFilm/liquidFilmBase.C b/src/regionFaModels/liquidFilm/liquidFilmBase.C
index 540bbe9da9..a1d36f1f55 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmBase.C
+++ b/src/regionFaModels/liquidFilm/liquidFilmBase.C
@@ -28,7 +28,6 @@ License
 #include "liquidFilmBase.H"
 #include "faMesh.H"
 #include "faCFD.H"
-#include "uniformDimensionedFields.H"
 #include "gravityMeshObject.H"
 #include "movingWallVelocityFvPatchVectorField.H"
 #include "turbulentFluidThermoModel.H"
@@ -60,11 +59,11 @@ bool liquidFilmBase::read(const dictionary& dict)
     regionFaModel::read(dict);
     if (active_)
     {
-        const dictionary& solution = this->solution().subDict("PISO");
+        const dictionary& solution = this->solution().subDict("PIMPLE");
         solution.readEntry("momentumPredictor", momentumPredictor_);
-        solution.readIfPresent("nOuterCorr", nOuterCorr_);
+        solution.readEntry("nOuterCorr", nOuterCorr_);
         solution.readEntry("nCorr", nCorr_);
-        solution.readEntry("nNonOrthCorr", nNonOrthCorr_);
+        solution.readEntry("nFilmCorr", nFilmCorr_);
     }
     return true;
 }
@@ -73,28 +72,24 @@ bool liquidFilmBase::read(const dictionary& dict)
 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_)
-        );
+    edgeScalarField SfUfbyDelta
+    (
+        regionMesh().edgeInterpolation::deltaCoeffs()*mag(phif_)
+    );
 
-        CoNum = max(SfUfbyDelta/regionMesh().magLe())
-            .value()*time().deltaT().value();
+    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();
 
-        velMag = max(mag(phif_)/regionMesh().magLe()).value();
-    }
+    reduce(CoNum, maxOp<scalar>());
+    reduce(velMag, maxOp<scalar>());
 
-    Info<< "Film Courant Number mean: " << meanCoNum
+    Info<< "Film Courant Number: "
         << " max: " << CoNum
-        << " Film velocity magnitude: " << velMag << endl;
+        << " Film velocity magnitude: (h)" << velMag << endl;
 
     return CoNum;
 }
@@ -114,19 +109,21 @@ liquidFilmBase::liquidFilmBase
 
     momentumPredictor_
     (
-        this->solution().subDict("PISO").lookup("momentumPredictor")
+        this->solution().subDict("PIMPLE").get<bool>("momentumPredictor")
     ),
     nOuterCorr_
     (
-        this->solution().subDict("PISO").lookupOrDefault("nOuterCorr", 1)
+        this->solution().subDict("PIMPLE").get<label>("nOuterCorr")
     ),
-    nCorr_(this->solution().subDict("PISO").get<label>("nCorr")),
-    nNonOrthCorr_
+    nCorr_(this->solution().subDict("PIMPLE").get<label>("nCorr")),
+    nFilmCorr_
     (
-        this->solution().subDict("PISO").get<label>("nNonOrthCorr")
+        this->solution().subDict("PIMPLE").get<label>("nFilmCorr")
     ),
 
-    h0_("h0", dimLength, SMALL),
+    h0_("h0", dimLength, 1e-7, dict),
+
+    deltaWet_("deltaWet", dimLength, 1e-4, dict),
 
     UName_(dict.get<word>("U")),
 
@@ -211,18 +208,6 @@ liquidFilmBase::liquidFilmBase
         fac::interpolate(h_*Uf_) & regionMesh().Le()
     ),
 
-    Tf_
-    (
-        IOobject
-        (
-            "Tf_" + regionName_,
-            primaryMesh().time().timeName(),
-            primaryMesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        regionMesh()
-    ),
 
     gn_
     (
@@ -238,6 +223,8 @@ liquidFilmBase::liquidFilmBase
         dimensionedScalar(dimAcceleration, Zero)
     ),
 
+    g_(meshObjects::gravity::New(primaryMesh().time())),
+
     massSource_
     (
         IOobject
@@ -294,10 +281,9 @@ liquidFilmBase::liquidFilmBase
 
     faOptions_(Foam::fa::options::New(p))
 {
-    const uniformDimensionedVectorField& g =meshObjects::gravity::New(time());
-
-    gn_ = (g & regionMesh().faceAreaNormals());
+    const areaVectorField& ns = regionMesh().faceAreaNormals();
 
+    gn_ = g_ & ns;
 
     if (!faOptions_.optionList::size())
     {
@@ -407,37 +393,7 @@ Foam::tmp<Foam::areaVectorField> liquidFilmBase::Up() const
     );
 
     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
@@ -481,7 +437,10 @@ tmp<areaScalarField> liquidFilmBase::pg() const
         const volScalarField& pp =
             primaryMesh().lookupObject<volScalarField>(pName_);
 
-        const volScalarField::Boundary& pw = pp.boundaryField();
+        volScalarField::Boundary& pw =
+            const_cast<volScalarField::Boundary&>(pp.boundaryField());
+
+        //pw -= pRef_;
 
         pfg.primitiveFieldRef() = vsmPtr_->mapInternalToSurface<scalar>(pw)();
     }
@@ -490,6 +449,30 @@ tmp<areaScalarField> liquidFilmBase::pg() const
 }
 
 
+tmp<areaScalarField> liquidFilmBase::alpha() const
+{
+    tmp<areaScalarField> talpha
+    (
+        new areaScalarField
+        (
+            IOobject
+            (
+                "talpha",
+                primaryMesh().time().timeName(),
+                primaryMesh()
+            ),
+            regionMesh(),
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    areaScalarField& alpha = talpha.ref();
+
+    alpha = pos0(h_ - deltaWet_);
+
+    return talpha;
+}
+
+
 void liquidFilmBase::addSources
 (
     const label patchi,
@@ -507,46 +490,65 @@ void liquidFilmBase::addSources
 
     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;
-    }
 }
 
 
+void liquidFilmBase::preEvolveRegion()
+{
+    regionFaModel::preEvolveRegion();
+}
+
+void liquidFilmBase::postEvolveRegion()
+{
+    if (debug && primaryMesh().time().writeTime())
+    {
+        massSource_.write();
+        pnSource_.write();
+        momentumSource_.write();
+    }
+
+    massSource_.boundaryFieldRef() = Zero;
+    pnSource_.boundaryFieldRef() = Zero;
+    momentumSource_.boundaryFieldRef() = Zero;
+
+    regionFaModel::postEvolveRegion();
+}
+
 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 uniformDimensionedVectorField& liquidFilmBase::g() const
+{
+    return g_;
+}
+
 const areaScalarField& liquidFilmBase::h() const
 {
      return h_;
 }
 
+const edgeScalarField& liquidFilmBase::phif() const
+{
+    return phif_;
+}
+
+const edgeScalarField& liquidFilmBase::phi2s() const
+{
+    return phi2s_;
+}
+
 const dimensionedScalar& liquidFilmBase::h0() const
 {
      return h0_;
@@ -557,6 +559,7 @@ const regionFaModel& liquidFilmBase::region() const
     return *this;
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace areaSurfaceFilmModels
diff --git a/src/regionFaModels/liquidFilm/liquidFilmBase.H b/src/regionFaModels/liquidFilm/liquidFilmBase.H
index e0f9522fee..172f15dbfc 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmBase.H
+++ b/src/regionFaModels/liquidFilm/liquidFilmBase.H
@@ -41,10 +41,10 @@ SourceFiles
 #include "autoPtr.H"
 #include "areaFieldsFwd.H"
 #include "volFieldsFwd.H"
+#include "uniformDimensionedFields.H"
 #include "regionFaModel.H"
 #include "faOptions.H"
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -54,7 +54,6 @@ namespace regionModels
 namespace areaSurfaceFilmModels
 {
 
-
 /*---------------------------------------------------------------------------*\
                          Class liquidFilmBase Declaration
 \*---------------------------------------------------------------------------*/
@@ -92,15 +91,19 @@ protected:
             //- Number of PISO-like correctors
             label nCorr_;
 
-            //- Number of non-orthogonal correctors
-            label nNonOrthCorr_;
+            //- Number of film thickness correctors
+            label nFilmCorr_;
 
             //- Cumulative continuity error
             scalar cumulativeContErr_;
 
-            //- Small delta
+            //- Smallest numerical thickness
             dimensionedScalar h0_;
 
+            //- Delta wet for sub-models
+            dimensionedScalar deltaWet_;
+
+
         //- Name of the velocity field
         word UName_;
 
@@ -131,14 +134,14 @@ protected:
             //- Film height flux
             edgeScalarField phi2s_;
 
-            //- Film temperature
-            areaScalarField Tf_;
-
             //- Normal gravity field
             areaScalarField gn_;
 
+            //- Gravity
+            uniformDimensionedVectorField g_;
 
-            // Mass exchanage fields
+
+            // Mass exchanage fields from the primary region (lagragian)
 
                 //- Mass
                 volScalarField massSource_;
@@ -156,10 +159,12 @@ protected:
                 scalar addedMassTotal_;
 
 
+
         //- Pointer to faOptions
         Foam::fa::options& faOptions_;
 
 
+
     // Protected Member Functions
 
         //- Read control parameters from dictionary
@@ -235,6 +240,10 @@ public:
             //- Map primary static pressure
             tmp<areaScalarField> pg() const;
 
+            //- Wet indicator using h0
+            tmp<areaScalarField> alpha() const;
+
+
         // Access functions
 
             //- Return faOptions
@@ -243,15 +252,22 @@ public:
             //- Access const reference Uf
             const areaVectorField& Uf() const;
 
-            //- Access const reference Tf
-            const areaScalarField& Tf() const;
-
             //- Access const reference gn
             const areaScalarField& gn() const;
 
+            //- Gravity
+            const uniformDimensionedVectorField& g() const;
+
             //- Access const reference h
             const areaScalarField& h() const;
 
+            //- Access to momentum flux
+            const edgeScalarField& phif() const;
+
+            //- Access continuity flux
+            const edgeScalarField& phi2s() const;
+
+
             //- Return h0
             const dimensionedScalar& h0() const;
 
@@ -259,6 +275,28 @@ public:
             const regionFaModel& region() const;
 
 
+            // Transfer fields - to the primary region (lagragian injection)
+
+                //- Return mass transfer source - Eulerian phase only
+                //virtual tmp<volScalarField> primaryMassTrans() const = 0;
+
+                //- Return the film mass available for transfer to cloud
+                virtual const volScalarField& cloudMassTrans() const = 0;
+
+                //- Return the parcel diameters originating from film to cloud
+                virtual const volScalarField& cloudDiameterTrans() const = 0;
+
+
+
+        // Evolution
+
+            //- Pre-evolve film
+            virtual void preEvolveRegion();
+
+            //- Post-evolve film
+            virtual void postEvolveRegion();
+
+
         // Thermo variables
 
              //- Access const reference mu
@@ -270,11 +308,16 @@ public:
             //- Access const reference sigma
             virtual const areaScalarField& sigma() const = 0;
 
+            //- Access const reference Tf
+            virtual const areaScalarField& Tf() const = 0;
+
+            //- Access const reference Cp
+            virtual const areaScalarField& Cp() const = 0;
+
 
         // External hook to add sources and mass exchange variables
 
 
-
             //- Add sources
             virtual void addSources
             (
@@ -285,7 +328,6 @@ public:
                 const scalar pressureSource,   // [kg.m/s] (normal momentum)
                 const scalar energySource = 0  // [J]
             );
-
 };
 
 
diff --git a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C
index b4a07ededc..306a09b895 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C
+++ b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C
@@ -31,6 +31,7 @@ License
 #include "gravityMeshObject.H"
 #include "volFields.H"
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -61,6 +62,7 @@ void liquidFilmModel::correctThermoFields()
         rho_[faceI] = thermo_.rho(pRef_, Tf_[faceI], X);
         mu_[faceI] = thermo_.mu(pRef_, Tf_[faceI], X);
         sigma_[faceI] = thermo_.sigma(pRef_, Tf_[faceI], X);
+        Cp_[faceI] = thermo_.Cp(pRef_, Tf_[faceI], X);
     }
 
     forAll (regionMesh().boundary(), patchI)
@@ -70,14 +72,19 @@ void liquidFilmModel::correctThermoFields()
         scalarField& patchRho = rho_.boundaryFieldRef()[patchI];
         scalarField& patchmu = mu_.boundaryFieldRef()[patchI];
         scalarField& patchsigma = sigma_.boundaryFieldRef()[patchI];
+        scalarField& patchCp = Cp_.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);
+            patchCp[edgeI] = thermo_.Cp(pRef_, patchTf[edgeI], X);
         }
     }
+
+    //Initialize pf_
+    pf_ = rho_*gn_*h_ - sigma_*fac::laplacian(h_);
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -119,6 +126,32 @@ liquidFilmModel::liquidFilmModel
         regionMesh(),
         dimensionedScalar(dimViscosity, Zero)
     ),
+    Tf_
+    (
+        IOobject
+        (
+            "Tf_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimTemperature, Zero)
+    ),
+    Cp_
+    (
+        IOobject
+        (
+            "Cp_" + regionName_,
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        regionMesh(),
+        dimensionedScalar(dimEnergy/dimTemperature, Zero)
+    ),
     sigma_
     (
         IOobject
@@ -151,7 +184,9 @@ liquidFilmModel::liquidFilmModel
         (
             "rhoSp",
             primaryMesh().time().timeName(),
-            primaryMesh()
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
         ),
         regionMesh(),
         dimensionedScalar(dimVelocity, Zero)
@@ -178,8 +213,49 @@ liquidFilmModel::liquidFilmModel
         regionMesh(),
         dimensionedScalar(dimPressure, Zero)
     ),
-    turbulence_(filmTurbulenceModel::New(*this, dict))
+    cloudMassTrans_
+    (
+        IOobject
+        (
+            "cloudMassTrans",
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        primaryMesh(),
+        dimensionedScalar(dimMass, Zero),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+
+    cloudDiameterTrans_
+    (
+        IOobject
+        (
+            "cloudDiameterTrans",
+            primaryMesh().time().timeName(),
+            primaryMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        primaryMesh(),
+        dimensionedScalar(dimLength, Zero),
+        calculatedFvPatchField<scalar>::typeName
+    ),
+
+    turbulence_(filmTurbulenceModel::New(*this, dict)),
+
+    availableMass_(regionMesh().faces().size(), Zero),
+
+    injection_(*this, dict),
+
+    forces_(*this, dict)
 {
+
+    if (dict.found("T0"))
+    {
+        Tf_ = dimensionedScalar("T0", dimTemperature, dict);
+    }
     correctThermoFields();
 }
 
@@ -208,8 +284,39 @@ const areaScalarField& liquidFilmModel::sigma() const
      return sigma_;
 }
 
+
+const areaScalarField& liquidFilmModel::Tf() const
+{
+     return Tf_;
+}
+
+
+const areaScalarField& liquidFilmModel::Cp() const
+{
+     return Cp_;
+}
+
+const volScalarField& liquidFilmModel::cloudMassTrans() const
+{
+    return cloudMassTrans_;
+}
+
+const volScalarField& liquidFilmModel::cloudDiameterTrans() const
+{
+    return cloudDiameterTrans_;
+}
+
+
+
 void liquidFilmModel::preEvolveRegion()
 {
+    liquidFilmBase::preEvolveRegion();
+
+
+
+    cloudMassTrans_ == dimensionedScalar(dimMass, Zero);
+    cloudDiameterTrans_ == dimensionedScalar(dimLength, Zero);
+
     const scalar deltaT = primaryMesh().time().deltaTValue();
     const scalarField rAreaDeltaT = 1/deltaT/regionMesh().S().field();
 
@@ -219,6 +326,7 @@ void liquidFilmModel::preEvolveRegion()
     // [kg.m/s]
     USp_.primitiveFieldRef() =
         vsm().mapToSurface(momentumSource_.boundaryField()[patchID()]);
+
     pnSp_.primitiveFieldRef() =
         vsm().mapToSurface(pnSource_.boundaryField()[patchID()]);
 
@@ -227,13 +335,38 @@ void liquidFilmModel::preEvolveRegion()
     rhoSp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
     USp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
     pnSp_.primitiveFieldRef() *= rAreaDeltaT/rho_;
+
+    rhoSp_.relax();
+    pnSp_.relax();
+    USp_.relax();
 }
 
 void liquidFilmModel::postEvolveRegion()
 {
-    massSource_ == dimensionedScalar(massSource_.dimensions(), Zero);
-    momentumSource_ == dimensionedVector(momentumSource_.dimensions(), Zero);
-    pnSource_ == dimensionedScalar(pnSource_.dimensions(), Zero);
+    availableMass_ = (h() - h0_)*rho()*regionMesh().S();
+    injection_.correct(availableMass_, cloudMassTrans_, cloudDiameterTrans_);
+    liquidFilmBase::postEvolveRegion();
+}
+
+
+void liquidFilmModel::info()
+{
+    Info<< "\nSurface film: " << type() << " on patch: " << patchID() << endl;
+
+    const DimensionedField<scalar, areaMesh>& sf = regionMesh().S();
+
+    Info<< indent << "min/max(mag(Uf))    = " << gMin(mag(Uf_.field())) << ", "
+        << gMax(mag(Uf_.field())) << nl
+        << indent << "min/max(delta)     = " << gMin(h_.field()) << ", " << gMax(h_.field()) << nl
+        << indent << "coverage           = "
+        << gSum(alpha()().field()*mag(sf.field()))/gSum(mag(sf.field())) <<  nl
+        << indent << "total mass         = "
+        << gSum(availableMass_) << nl;
+
+
+    Info<< indent << CourantNumber() << endl;
+
+    injection_.info(Info);
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H
index 0aeebf6118..89225b70ca 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H
+++ b/src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.H
@@ -42,8 +42,16 @@ SourceFiles
 #include "volFieldsFwd.H"
 #include "liquidFilmBase.H"
 #include "faMesh.H"
+
 #include "filmTurbulenceModel.H"
+
 #include "liquidMixtureProperties.H"
+#include "injectionModelList.H"
+
+#include "faCFD.H"
+
+//#include "transferModelList.H"
+#include "forceList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,6 +62,7 @@ namespace regionModels
 namespace areaSurfaceFilmModels
 {
 
+
 /*---------------------------------------------------------------------------*\
                         Class liquidFilmModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -94,6 +103,12 @@ protected:
                 //- Dynamic viscosity [Pa.s]
                 areaScalarField mu_;
 
+                //- Film temperature
+                areaScalarField Tf_;
+
+                 //- Film Heat capacity
+                areaScalarField Cp_;
+
                 //- Surface tension [m/s2]
                 areaScalarField sigma_;
 
@@ -112,6 +127,15 @@ protected:
                     areaScalarField pnSp_;
 
 
+                 // Transfer fields
+
+                    //- Film mass  for transfer to cloud
+                    volScalarField cloudMassTrans_;
+
+                    //- Parcel diameters originating from film to cloud
+                    volScalarField cloudDiameterTrans_;
+
+
     // General properties
 
 
@@ -120,6 +144,21 @@ protected:
             autoPtr<filmTurbulenceModel> turbulence_;
 
 
+            // Sub-models
+
+                //- Available mass for transfer via sub-models
+                scalarField availableMass_;
+
+                //- Cloud injection
+                injectionModelList injection_;
+
+                //- Transfer with the continuous phase
+                //transferModelList transfer_;
+
+                //- List of film forces
+                forceList forces_;
+
+
     // Protected member functions
 
 
@@ -169,6 +208,20 @@ public:
             //- Access const reference sigma
             const areaScalarField& sigma() const;
 
+            //- Access const reference Tf
+            const areaScalarField& Tf() const;
+
+            //- Access const reference Cp
+            const areaScalarField& Cp() const;
+
+
+            // Transfer fields - to the primary region (lagragian injection)
+
+                //- Return the film mass available for transfer to cloud
+                virtual const volScalarField& cloudMassTrans() const;
+
+                //- Return the parcel diameters originating from film to cloud
+                virtual const volScalarField& cloudDiameterTrans() const;
 
         // Evolution
 
@@ -185,7 +238,7 @@ public:
        // I-O
 
             //- Provide some feedback
-            //virtual void info();
+            virtual void info();
 };
 
 
diff --git a/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.C b/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.C
new file mode 100644
index 0000000000..9b072df6db
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.C
@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 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 "filmSubModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+filmSubModelBase::filmSubModelBase(liquidFilmBase& film)
+:
+    subModelBase(film.outputProperties()),
+    filmModel_(film)
+{}
+
+
+filmSubModelBase::filmSubModelBase
+(
+    liquidFilmBase& film,
+    const dictionary& dict,
+    const word& baseName,
+    const word& modelType,
+    const word& dictExt
+)
+:
+    subModelBase
+    (
+        film.outputProperties(),
+        dict,
+        baseName,
+        modelType,
+        dictExt
+    ),
+    filmModel_(film)
+{}
+
+
+filmSubModelBase::filmSubModelBase
+(
+    const word& modelName,
+    liquidFilmBase& film,
+    const dictionary& dict,
+    const word& baseName,
+    const word& modelType
+)
+:
+    subModelBase
+    (
+        modelName,
+        film.outputProperties(),
+        dict,
+        baseName,
+        modelType
+    ),
+    filmModel_(film)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+filmSubModelBase::~filmSubModelBase()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool filmSubModelBase::writeTime() const
+{
+    return active() && filmModel_.time().writeTime();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.H b/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.H
new file mode 100644
index 0000000000..12afa909d5
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/filmSubModelBase.H
@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 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::surfaceFilmModels::filmSubModelBase
+
+Description
+    Base class for surface film sub-models
+
+SourceFiles
+    filmSubModelBaseI.H
+    filmSubModelBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filmSubModelBase_H
+#define filmSubModelBase_H
+
+#include "liquidFilmBase.H"
+#include "subModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class filmSubModelBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class filmSubModelBase
+:
+    public subModelBase
+{
+protected:
+
+    // Protected Data
+
+        //- Reference to the film surface film model
+        liquidFilmBase& filmModel_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        filmSubModelBase(liquidFilmBase& film);
+
+        //- Construct from film film without name
+        filmSubModelBase
+        (
+            liquidFilmBase& film,
+            const dictionary& dict,
+            const word& baseName,
+            const word& modelType,
+            const word& dictExt = "Coeffs"
+        );
+
+        //- Construct from film film with name
+        filmSubModelBase
+        (
+            const word& modelName,
+            liquidFilmBase& film,
+            const dictionary& dict,
+            const word& baseName,
+            const word& modelType
+        );
+
+
+    //- Destructor
+    virtual ~filmSubModelBase();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Flag to indicate when to write a property
+            virtual bool writeTime() const;
+
+            //- Return const access to the film surface film model
+            inline const liquidFilmBase& film() const;
+
+            //- Return the reference to the film surface film model
+            inline liquidFilmBase& film();
+
+            template<class FilmType>
+            inline const FilmType& filmType() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "filmSubModelBaseI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "filmSubModelBaseTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseI.H b/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseI.H
new file mode 100644
index 0000000000..70f07eb03f
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseI.H
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline const liquidFilmBase& filmSubModelBase::film() const
+{
+    return filmModel_;
+}
+
+
+inline liquidFilmBase& filmSubModelBase::film()
+{
+    return filmModel_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseTemplates.C b/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseTemplates.C
new file mode 100644
index 0000000000..355f1a8851
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/filmSubModelBaseTemplates.C
@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class FilmType>
+const FilmType& filmSubModelBase::filmType() const
+{
+    if (!isA<FilmType>(filmModel_))
+    {
+        FatalErrorInFunction
+            << "Model " << this->modelType() << " requested film type "
+            << FilmType::typeName << " but film is type " << filmModel_.type()
+            << abort(FatalError);
+    }
+
+    return refCast<const FilmType>(filmModel_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
index 6c0d5a2731..cec1b7acb9 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
@@ -106,7 +106,7 @@ tmp<areaScalarField> filmTurbulenceModel::Cw() const
             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.primitiveFieldRef() = 3*mu/(h + h0);
             Cw.min(5000.0);
             break;
         }
@@ -141,7 +141,7 @@ tmp<areaScalarField> filmTurbulenceModel::Cw() const
             const scalar h0 = film_.h0().value();
 
             Cw.primitiveFieldRef() =
-                sqr(n)*mag(g.value())*mag(Uf)/cbrt(h+ h0);
+                sqr(n)*mag(g.value())*mag(Uf)/cbrt(h+h0);
             break;
         }
         default:
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
index ce2b620af9..5277983eec 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef filmTurbulenceModel_H
-#define filmTurbulenceModel_H
+#ifndef regionFaModels_filmTurbulenceModel_H
+#define regionFaModels_filmTurbulenceModel_H
 
 #include "areaFieldsFwd.H"
 #include "runTimeSelectionTables.H"
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
index 472f80dee2..7630e25466 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
@@ -96,13 +96,13 @@ tmp<faVectorMatrix> laminar::Su(areaVectorField& U) const
     tmp<areaVectorField> Up = film_.Up();
 
     // employ simple coeff-based model
-    areaScalarField Cs("Cs", Cf_*mag(Up() - U));
+    const dimensionedScalar Cf("Cf", dimVelocity, Cf_);
 
     tmp<areaScalarField> wf = Cw();
 
     return
     (
-       - fam::Sp(Cs, U) + Cs*Up()     // surface contribution
+       - fam::Sp(Cf, U) + Cf*Up()     // surface contribution
        - fam::Sp(wf(), U) + wf()*Uw() // wall contribution
     );
 }
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
index a405929e23..b0e6b70541 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
@@ -34,10 +34,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef laminar_H
-#define laminar_H
+#ifndef regionFaModels_surfaceFilmModels_laminar_H
+#define regionFaModels_surfaceFilmModels_laminar_H
 
-#include "filmTurbulenceModel.H"
+#include "liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H"
 #include "faCFD.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
new file mode 100644
index 0000000000..ee40e48424
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
@@ -0,0 +1,324 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    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 "curvatureSeparation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "stringListOps.H"
+#include "cyclicPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(curvatureSeparation, 0);
+addToRunTimeSelectionTable
+(
+    injectionModel,
+    curvatureSeparation,
+    dictionary
+);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<areaScalarField> curvatureSeparation::calcInvR1
+(
+    const areaVectorField& U,
+    const scalarField& calcCosAngle
+) const
+{
+    const dimensionedScalar smallU(dimVelocity, ROOTVSMALL);
+    const areaVectorField UHat(U/(mag(U) + smallU));
+    tmp<areaScalarField> tinvR1
+    (
+        new areaScalarField("invR1", UHat & (UHat & -gradNHat_))
+    );
+
+    scalarField& invR1 = tinvR1.ref().primitiveFieldRef();
+
+    // apply defined patch radii
+    const scalar rMin = 1e-6;
+    const scalar definedInvR1 = 1.0/max(rMin, definedPatchRadii_);
+
+    if (definedPatchRadii_ > 0)
+    {
+        invR1 = definedInvR1;
+    }
+
+    // filter out large radii
+    const scalar rMax = 1e6;
+    forAll(invR1, i)
+    {
+        if ((mag(invR1[i]) < 1/rMax))
+        {
+            invR1[i] = -1.0;
+        }
+    }
+
+    return tinvR1;
+}
+
+
+tmp<scalarField> curvatureSeparation::calcCosAngle
+(
+    const edgeScalarField& phi
+) const
+{
+    const areaVectorField& U = film().Uf();
+    const dimensionedScalar smallU(dimVelocity, ROOTVSMALL);
+    const areaVectorField UHat(U/(mag(U) + smallU));
+
+    const faMesh& mesh = film().regionMesh();
+    const labelUList& own = mesh.edgeOwner();
+    const labelUList& nbr = mesh.edgeNeighbour();
+
+    scalarField phiMax(mesh.nFaces(), -GREAT);
+    scalarField cosAngle(UHat.size(), Zero);
+
+    const scalarField invR1(calcInvR1(U, cosAngle));
+
+    forAll(nbr, edgei)
+    {
+        const label cellO = own[edgei];
+        const label cellN = nbr[edgei];
+
+        if (phi[edgei] > phiMax[cellO])
+        {
+            phiMax[cellO] = phi[edgei];
+            cosAngle[cellO] = -gHat_ & UHat[cellN];
+        }
+        if (-phi[edgei] > phiMax[cellN])
+        {
+            phiMax[cellN] = -phi[edgei];
+            cosAngle[cellN] = -gHat_ & UHat[cellO];
+        }
+    }
+
+    cosAngle *= pos(invR1);
+
+    // checks
+    if (debug && mesh.time().writeTime())
+    {
+        areaScalarField volCosAngle
+        (
+            IOobject
+            (
+                "cosAngle",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            film().regionMesh(),
+            dimensionedScalar(dimless, Zero)
+        );
+        volCosAngle.primitiveFieldRef() = cosAngle;
+        volCosAngle.correctBoundaryConditions();
+        volCosAngle.write();
+    }
+
+    return max(min(cosAngle, scalar(1)), scalar(-1));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+curvatureSeparation::curvatureSeparation
+(
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    injectionModel(type(), film, dict),
+    gradNHat_(fac::grad(film.regionMesh().faceAreaNormals())),
+    deltaByR1Min_(coeffDict_.getOrDefault<scalar>("deltaByR1Min", 0)),
+    definedPatchRadii_
+    (
+        coeffDict_.getOrDefault<scalar>("definedPatchRadii", 0)
+    ),
+    magG_(mag(film.g().value())),
+    gHat_(Zero),
+    fThreshold_
+    (
+        coeffDict_.getOrDefault<scalar>("fThreshold", 1e-8)
+    ),
+    minInvR1_
+    (
+        coeffDict_.getOrDefault<scalar>("minInvR1", 5)
+    )
+{
+    if (magG_ < ROOTVSMALL)
+    {
+        FatalErrorInFunction
+            << "Acceleration due to gravity must be non-zero"
+            << exit(FatalError);
+    }
+
+    gHat_ = film.g().value()/magG_;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void curvatureSeparation::correct
+(
+    scalarField& availableMass,
+    scalarField& massToInject,
+    scalarField& diameterToInject
+)
+{
+    const faMesh& mesh = film().regionMesh();
+
+    const areaScalarField& delta = film().h();
+    const areaVectorField& U = film().Uf();
+    const edgeScalarField& phi = film().phi2s();
+    const areaScalarField& rho = film().rho();
+    const scalarField magSqrU(magSqr(film().Uf()));
+    const areaScalarField& sigma = film().sigma();
+
+    const scalarField cosAngle(calcCosAngle(phi));
+    const scalarField invR1(calcInvR1(U, cosAngle));
+
+
+    // calculate force balance
+    scalarField Fnet(mesh.nFaces(), Zero);
+    scalarField separated(mesh.nFaces(), Zero);
+
+    forAll(invR1, i)
+    {
+        if ((invR1[i] > minInvR1_) && (delta[i]*invR1[i] > deltaByR1Min_))
+        {
+            const scalar R1 = 1.0/(invR1[i] + ROOTVSMALL);
+            const scalar R2 = R1 + delta[i];
+
+            // inertial force
+            const scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i];
+
+            // body force
+            const scalar Fb =
+               - 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];
+
+            // surface force
+            const scalar Fs = sigma[i]/R2;
+
+            Fnet[i] = Fi + Fb + Fs;
+
+            if (Fnet[i] + fThreshold_ < 0)
+            {
+                separated[i] = 1.0;
+            }
+        }
+    }
+
+    // inject all available mass
+    massToInject = separated*availableMass;
+    diameterToInject = separated*delta;
+    availableMass -= separated*availableMass;
+
+    addToInjectedMass(sum(massToInject));
+
+    if (debug && mesh.time().writeTime())
+    {
+        areaScalarField volFnet
+        (
+            IOobject
+            (
+                "Fnet",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar(dimForce, Zero)
+        );
+        volFnet.primitiveFieldRef() = Fnet;
+        volFnet.write();
+
+        areaScalarField areaSeparated
+        (
+            IOobject
+            (
+                "separated",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar(dimMass, Zero)
+        );
+        areaSeparated.primitiveFieldRef() = separated;
+        areaSeparated.write();
+
+        areaScalarField areaMassToInject
+        (
+            IOobject
+            (
+                "massToInject",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar(dimMass, Zero)
+        );
+        areaMassToInject.primitiveFieldRef() = massToInject;
+        areaMassToInject.write();
+
+        areaScalarField areaInvR1
+        (
+            IOobject
+            (
+                "InvR1",
+                film().primaryMesh().time().timeName(),
+                film().primaryMesh(),
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar(inv(dimLength), Zero)
+        );
+        areaInvR1.primitiveFieldRef() = invR1;
+        areaInvR1.write();
+    }
+
+    injectionModel::correct();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H
new file mode 100644
index 0000000000..0ec037e432
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H
@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2021 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::curvatureSeparation
+
+Description
+    Curvature film separation model
+
+    Assesses film curvature via the mesh geometry and calculates a force
+    balance of the form:
+
+        F_sum = F_inertial + F_body + F_surface
+
+    If F_sum < 0, the film separates. Similarly, if F_sum > 0 the film will
+    remain attached.
+
+    Based on description given by
+        Owen and D. J. Ryley. The flow of thin liquid films around corners.
+        International Journal of Multiphase Flow, 11(1):51-62, 1985.
+
+SourceFiles
+    curvatureSeparation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef curvatureSeparation_H
+#define curvatureSeparation_H
+
+#include "injectionModel.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class curvatureSeparation Declaration
+\*---------------------------------------------------------------------------*/
+
+class curvatureSeparation
+:
+    public injectionModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        curvatureSeparation(const curvatureSeparation&) = delete;
+
+        //- No copy assignment
+        void operator=(const curvatureSeparation&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        //- Gradient of surface normals
+        areaTensorField gradNHat_;
+
+        //- Minimum gravity driven film thickness (non-dimensionalised delta/R1)
+        scalar deltaByR1Min_;
+
+        //- List of radii for patches - if patch not defined, radius
+        //- calculated based on mesh geometry
+        scalar definedPatchRadii_;
+
+        //- Magnitude of gravity vector
+        scalar magG_;
+
+        //- Direction of gravity vector
+        vector gHat_;
+
+        //- Threshold force for separation
+        scalar fThreshold_;
+
+        //- Minimum inv R1 for separation
+        scalar minInvR1_;
+
+
+    // Protected Member Functions
+
+        //- Calculate local (inverse) radius of curvature
+        tmp<areaScalarField> calcInvR1
+        (
+            const areaVectorField& U,
+            const scalarField& calcCosAngle
+        ) const;
+
+        //- Calculate the cosine of the angle between gravity vector and
+        //- cell out flow direction
+        tmp<scalarField> calcCosAngle(const edgeScalarField& phi) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("curvatureSeparation");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        curvatureSeparation
+        (
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~curvatureSeparation() = default;
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                scalarField& availableMass,
+                scalarField& massToInject,
+                scalarField& diameterToInject
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
new file mode 100644
index 0000000000..b1e5475a85
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 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 "injectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(injectionModel, 0);
+defineRunTimeSelectionTable(injectionModel, dictionary);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void injectionModel::addToInjectedMass(const scalar dMass)
+{
+    injectedMass_ += dMass;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+injectionModel::injectionModel(liquidFilmBase& film)
+:
+    filmSubModelBase(film),
+    injectedMass_(0.0)
+{}
+
+
+injectionModel::injectionModel
+(
+    const word& modelType,
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    filmSubModelBase(film, dict, typeName, modelType),
+    injectedMass_(0.0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+injectionModel::~injectionModel()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void injectionModel::correct()
+{
+    if (writeTime())
+    {
+        scalar injectedMass0 = getModelProperty<scalar>("injectedMass");
+        injectedMass0 += returnReduce(injectedMass_, sumOp<scalar>());
+        setModelProperty<scalar>("injectedMass", injectedMass0);
+        injectedMass_ = 0.0;
+    }
+}
+
+
+scalar injectionModel::injectedMassTotal() const
+{
+    const scalar injectedMass0 = getModelProperty<scalar>("injectedMass");
+    return injectedMass0 + returnReduce(injectedMass_, sumOp<scalar>());
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.H
new file mode 100644
index 0000000000..5bf117f7a9
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.H
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 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::injectionModel
+
+Description
+    Base class for film injection models, handling mass transfer from the
+    film.
+
+SourceFiles
+    injectionModel.C
+    injectionModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef injectionModel_H
+#define injectionModel_H
+
+#include "filmSubModelBase.H"
+#include "runTimeSelectionTables.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class injectionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class injectionModel
+:
+    public filmSubModelBase
+{
+    // Private Data
+
+        //- Injected mass
+        scalar injectedMass_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        injectionModel(const injectionModel&) = delete;
+
+        //- No copy assignment
+        void operator=(const injectionModel&) = delete;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Add to injected mass
+        void addToInjectedMass(const scalar dMass);
+
+        //- Correct
+        void correct();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("injectionModel");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             injectionModel,
+             dictionary,
+             (
+                liquidFilmBase& film,
+                const dictionary& dict
+             ),
+             (film, dict)
+         );
+
+
+    // Constructors
+
+        //- Construct null
+        injectionModel(liquidFilmBase& film);
+
+        //- Construct from type name, dictionary and surface film model
+        injectionModel
+        (
+            const word& modelType,
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected injection model
+        static autoPtr<injectionModel> New
+        (
+            liquidFilmBase& film,
+            const dictionary& dict,
+            const word& mdoelType
+        );
+
+
+    //- Destructor
+    virtual ~injectionModel();
+
+
+    // Member Functions
+
+        //- Correct
+        virtual void correct
+        (
+            scalarField& availableMass,
+            scalarField& massToInject,
+            scalarField& diameterToInject
+        ) = 0;
+
+        //- Return the total mass injected
+        virtual scalar injectedMassTotal() const;
+
+        //- Accumulate the total mass injected for the patches into the
+        //- scalarField provided
+        virtual void patchInjectedMassTotals(scalar& patchMasses) const
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
new file mode 100644
index 0000000000..e449744693
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 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 "injectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<injectionModel> injectionModel::New
+(
+    liquidFilmBase& model,
+    const dictionary& dict,
+    const word& modelType
+)
+{
+    Info<< "        " << modelType << endl;
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "injectionModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return autoPtr<injectionModel>(cstrIter()(model, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
new file mode 100644
index 0000000000..01900da387
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "injectionModelList.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+injectionModelList::injectionModelList(liquidFilmBase& film)
+:
+    PtrList<injectionModel>(),
+    filmSubModelBase(film)
+{}
+
+
+injectionModelList::injectionModelList
+(
+    liquidFilmBase& film,
+    const dictionary& dict
+)
+:
+    PtrList<injectionModel>(),
+    filmSubModelBase
+    (
+        "injectionModelList",
+        film,
+        dict,
+        "injectionModelList",
+        "injectionModelList"
+    ),
+    massInjected_(Zero)
+{
+    const wordList activeModels(dict.lookup("injectionModels"));
+
+    wordHashSet models(activeModels);
+
+    Info<< "    Selecting film injection models" << endl;
+    if (models.size())
+    {
+        this->setSize(models.size());
+
+        label i = 0;
+        for (const word& model : models)
+        {
+            set(i, injectionModel::New(film, dict, model));
+            i++;
+        }
+    }
+    else
+    {
+        Info<< "        none" << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+injectionModelList::~injectionModelList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void injectionModelList::correct
+(
+    scalarField& availableMass,
+    volScalarField& massToInject,
+    volScalarField& diameterToInject
+)
+{
+    const label patchi = film().patchID();
+
+    // Correct models that accumulate mass and diameter transfers
+    forAll(*this, i)
+    {
+        injectionModel& im = operator[](i);
+        im.correct
+        (
+            availableMass,
+            massToInject.boundaryFieldRef()[patchi],
+            diameterToInject.boundaryFieldRef()[patchi]
+        );
+    }
+
+    massInjected_ += gSum(massToInject.boundaryField()[patchi]);
+}
+
+
+void injectionModelList::info(Ostream& os)
+{
+    const polyBoundaryMesh& pbm = film().primaryMesh().boundaryMesh();
+
+    scalar injectedMass = 0;
+    scalar patchInjectedMasses = 0;
+
+    forAll(*this, i)
+    {
+        const injectionModel& im = operator[](i);
+        injectedMass += im.injectedMassTotal();
+        im.patchInjectedMassTotals(patchInjectedMasses);
+    }
+
+    os  << indent << "injected mass      = " << injectedMass << nl;
+
+    const label patchi = film().patchID();
+
+    if (mag(patchInjectedMasses) > VSMALL)
+    {
+        os  << indent << indent << "from patch " << pbm[patchi].name()
+            << " = " << patchInjectedMasses << nl;
+    }
+
+    scalar mass0(Zero);
+    this->getBaseProperty("massInjected", mass0);
+
+    scalar mass(massInjected_);
+    mass += mass0;
+
+    Info<< indent << "  - patch: " << pbm[patchi].name() << "  "
+        << mass << endl;
+
+    if (film().primaryMesh().time().writeTime())
+    {
+        setBaseProperty("massInjected", mass);
+        massInjected_ = 0.0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.H
new file mode 100644
index 0000000000..b159f9486a
--- /dev/null
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.H
@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::surfaceFilmModels::injectionModelList
+
+Description
+    List container for film injection models
+
+SourceFiles
+    injectionModelList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef injectionModelList_H
+#define injectionModelList_H
+
+#include "PtrList.H"
+#include "injectionModel.H"
+#include "filmSubModelBase.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace areaSurfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class injectionModelList Declaration
+\*---------------------------------------------------------------------------*/
+
+class injectionModelList
+:
+    public PtrList<injectionModel>,
+    public filmSubModelBase
+{
+    // Private Data
+
+        //- List of mass injected per patch
+        scalar massInjected_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        injectionModelList(const injectionModelList&) = delete;
+
+        //- No copy assignment
+        void operator=(const injectionModelList&) = delete;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        injectionModelList(liquidFilmBase& film);
+
+        //- Construct from type name, dictionary and surface film model
+        injectionModelList
+        (
+            liquidFilmBase& film,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~injectionModelList();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                scalarField& availableMass,
+                volScalarField& massToInject,
+                volScalarField& diameterToInject
+            );
+
+
+        // I-O
+
+            //- Provide some info
+            virtual void info(Ostream& os);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace areaSurfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //