diff --git a/src/regionFaModels/liquidFilm/liquidFilmBase.C b/src/regionFaModels/liquidFilm/liquidFilmBase.C
index 5ae39554858c9408d6e33e6e98e1db7145ed9722..469d97bb5ea3f0f5ec3ca143c6885e1d63450062 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmBase.C
+++ b/src/regionFaModels/liquidFilm/liquidFilmBase.C
@@ -551,6 +551,12 @@ scalar liquidFilmBase::pRef()
     return pRef_;
 }
 
+
+word liquidFilmBase::UName() const
+{
+    return UName_;
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace areaSurfaceFilmModels
diff --git a/src/regionFaModels/liquidFilm/liquidFilmBase.H b/src/regionFaModels/liquidFilm/liquidFilmBase.H
index 8bbce9f9b5213acb2d24a4c7b13fe89613ba4eb1..16470f119617f568fb5f1266af211cd7cd6e76e0 100644
--- a/src/regionFaModels/liquidFilm/liquidFilmBase.H
+++ b/src/regionFaModels/liquidFilm/liquidFilmBase.H
@@ -261,6 +261,9 @@ public:
             //- Access to pRef
             scalar pRef();
 
+            //- Name of the U field
+            word UName() const;
+
 
             // Transfer fields - to the primary region (lagragian injection)
 
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
index d8d18116fe96e5feaab1131481cf2e9b9c138090..c7fd810867b92b61d1c858af18ef98fada489efd 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +27,8 @@ License
 
 #include "filmTurbulenceModel.H"
 #include "gravityMeshObject.H"
+#include "turbulentTransportModel.H"
+#include "turbulentFluidThermoModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,6 +56,16 @@ filmTurbulenceModel::frictionMethodTypeNames_
     { frictionMethodType::mManningStrickler, "ManningStrickler" }
 };
 
+const Enum
+<
+    filmTurbulenceModel::shearMethodType
+>
+filmTurbulenceModel::shearMethodTypeNames_
+{
+    { shearMethodType::msimple, "simple" },
+    { shearMethodType::mwallFunction, "wallFunction" }
+};
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 filmTurbulenceModel::filmTurbulenceModel
@@ -65,8 +77,16 @@ filmTurbulenceModel::filmTurbulenceModel
 :
     film_(film),
     dict_(dict.subDict(modelType + "Coeffs")),
-    method_(frictionMethodTypeNames_.get("friction", dict_))
-{}
+    method_(frictionMethodTypeNames_.get("friction", dict_)),
+    shearMethod_(shearMethodTypeNames_.get("shearStress", dict_)),
+    rhoName_(dict_.getOrDefault<word>("rho", "rho")),
+    rhoRef_(VGREAT)
+{
+    if (rhoName_ == "rhoInf")
+    {
+        rhoRef_ = dict_.get<scalar>("rhoInf");
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -163,6 +183,173 @@ tmp<areaScalarField> filmTurbulenceModel::Cw() const
 }
 
 
+tmp<faVectorMatrix> filmTurbulenceModel::primaryRegionFriction
+(
+    areaVectorField& U
+) const
+{
+    tmp<faVectorMatrix> tshearStress
+    (
+        new faVectorMatrix(U, sqr(U.dimensions())*sqr(dimLength))
+    );
+
+    switch (shearMethod_)
+    {
+        case msimple:
+        {
+            tmp<areaVectorField> Up = film_.Up();
+
+            const dimensionedScalar Cf
+            (
+                "Cf",
+                dimVelocity,
+                dict_.get<scalar>("Cf")
+            );
+
+            tshearStress.ref() += - fam::Sp(Cf, U) + Cf*Up();
+
+            break;
+        }
+        case mwallFunction:
+        {
+            tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
+
+            const volSymmTensorField::Boundary& devRhoReffb
+                = tdevRhoReff().boundaryField();
+
+            label patchi = film_.patchID();
+
+            const surfaceVectorField::Boundary& Sfb =
+                film_.primaryMesh().Sf().boundaryField();
+
+            vectorField fT(Sfb[patchi] & devRhoReffb[patchi]);
+
+            const vectorField& nHat =
+                film_.regionMesh().faceAreaNormals().internalField();
+
+            // Substract normal component
+            fT -= nHat*(fT & nHat);
+
+            tmp<areaVectorField> taForce
+            (
+                new areaVectorField
+                (
+                    IOobject
+                    (
+                        "taForce",
+                        film_.primaryMesh().time().timeName(),
+                        film_.primaryMesh()
+                    ),
+                    film_.regionMesh(),
+                    dimensionedVector(sqr(dimVelocity), Zero)
+                )
+            );
+
+            vectorField& aForce = taForce.ref().primitiveFieldRef();
+
+            // Map ft to surface
+            vectorField afT(film_.vsm().mapToSurface(fT));
+
+            const DimensionedField<scalar, areaMesh>& magSf =
+                film_.regionMesh().S();
+
+            aForce = afT/film_.rho().primitiveField()/magSf;
+
+            tshearStress.ref() += taForce();
+
+            if (film_.regionMesh().time().writeTime())
+            {
+                taForce().write();
+            }
+
+            break;
+        }
+
+    }
+
+    return tshearStress;
+
+}
+
+
+tmp<Foam::volSymmTensorField> filmTurbulenceModel::devRhoReff() const
+{
+    typedef compressible::turbulenceModel cmpTurbModel;
+    typedef incompressible::turbulenceModel icoTurbModel;
+
+    const fvMesh& m = film_.primaryMesh();
+
+    const auto& U = m.lookupObject<volVectorField>(film_.UName());
+
+    if (m.foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
+    {
+        const cmpTurbModel& turb =
+            m.lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);
+
+        return turb.devRhoReff();
+    }
+    else if (m.foundObject<icoTurbModel>(icoTurbModel::propertiesName))
+    {
+        const incompressible::turbulenceModel& turb =
+            m.lookupObject<icoTurbModel>(icoTurbModel::propertiesName);
+
+        return rho()*turb.devReff();
+    }
+    else if (m.foundObject<fluidThermo>(fluidThermo::dictName))
+    {
+        const fluidThermo& thermo =
+            m.lookupObject<fluidThermo>(fluidThermo::dictName);
+
+        return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
+    }
+    else if (m.foundObject<transportModel>("transportProperties"))
+    {
+        const transportModel& laminarT =
+            m.lookupObject<transportModel>("transportProperties");
+
+        return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
+    }
+    else if (m.foundObject<dictionary>("transportProperties"))
+    {
+        const dictionary& transportProperties =
+            m.lookupObject<dictionary>("transportProperties");
+
+        dimensionedScalar nu("nu", dimViscosity, transportProperties);
+
+        return -rho()*nu*dev(twoSymm(fvc::grad(U)));
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "No valid model for viscous stress calculation"
+            << exit(FatalError);
+
+        return volSymmTensorField::null();
+    }
+}
+
+
+tmp<Foam::volScalarField> filmTurbulenceModel::rho() const
+{
+    const fvMesh& m = film_.primaryMesh();
+    if (rhoName_ == "rhoInf")
+    {
+        return tmp<volScalarField>::New
+        (
+            IOobject
+            (
+                "rho",
+                m.time().timeName(),
+                m
+            ),
+            m,
+            dimensionedScalar("rho", dimDensity, rhoRef_)
+        );
+    }
+
+    return(m.lookupObject<volScalarField>(rhoName_));
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace areaSurfaceFilmModels
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
index 5277983eecc508f2507272905a8c714c7cd8020e..6234880139c094175dfc14d96e8aaad3f1a3923f 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -82,6 +82,12 @@ public:
             mManningStrickler
         };
 
+        enum shearMethodType
+        {
+            msimple,
+            mwallFunction
+        };
+
 
 protected:
 
@@ -93,12 +99,24 @@ protected:
         //- Names for friction models
         static const Enum<frictionMethodType> frictionMethodTypeNames_;
 
+        //- Names for shear stress models
+        static const Enum<shearMethodType> shearMethodTypeNames_;
+
         //- Model dictionary
         const dictionary dict_;
 
         //- Method used
         const frictionMethodType method_;
 
+        //- Shear method used
+        const shearMethodType shearMethod_;
+
+        //- Name of density field (optional)
+        word rhoName_;
+
+        //- Reference density needed for incompressible calculations
+        scalar rhoRef_;
+
 
 public:
 
@@ -154,7 +172,19 @@ public:
             const liquidFilmBase& film() const;
 
 
-        // Evolution
+        // Turbulence
+
+             //- Return the effective viscous stress (laminar + turbulent).
+            tmp<volSymmTensorField> devRhoReff() const;
+
+            //- Return primary region frintion
+            tmp<faVectorMatrix> primaryRegionFriction
+            (
+                areaVectorField& U
+            ) const;
+
+            //- Return rho if specified otherwise rhoRef
+            tmp<volScalarField> rho() const;
 
             //- Return the wall film surface friction
             virtual tmp<areaScalarField> Cw() const;
@@ -162,6 +192,10 @@ public:
             //- Return the film turbulence viscosity
             virtual tmp<areaScalarField> mut() const = 0;
 
+
+        // Evolution
+
+
             //- Correct/update the model
             virtual void correct() = 0;
 
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
index 7630e254661109ce28688d7d664dd581ad1a3489..b00b5b38b0dbdb818bac70c997ac1b305dd15fd5 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,8 +51,8 @@ laminar::laminar
     const dictionary& dict
 )
 :
-    filmTurbulenceModel(type(), film, dict),
-    Cf_(dict_.get<scalar>("Cf"))
+    filmTurbulenceModel(type(), film, dict)
+//     Cf_(dict_.get<scalar>("Cf"))
 {}
 
 
@@ -91,23 +91,40 @@ void laminar::correct()
 
 tmp<faVectorMatrix> laminar::Su(areaVectorField& U) const
 {
-    // local references to film fields
-    tmp<areaVectorField> Uw = film_.Uw();
-    tmp<areaVectorField> Up = film_.Up();
+    return
+    (
+        primaryRegionFriction(U) + wallFriction(U)
+    );
+}
 
-    // employ simple coeff-based model
-    const dimensionedScalar Cf("Cf", dimVelocity, Cf_);
 
+tmp<faVectorMatrix> laminar::wallFriction(areaVectorField& U) const
+{
+     // local references to film fields
+    tmp<areaVectorField> Uw = film_.Uw();
     tmp<areaScalarField> wf = Cw();
 
     return
     (
-       - fam::Sp(Cf, U) + Cf*Up()     // surface contribution
        - fam::Sp(wf(), U) + wf()*Uw() // wall contribution
     );
 }
 
 
+// tmp<faVectorMatrix> laminar::primaryRegionFriction(areaVectorField& U) const
+// {
+//     tmp<areaVectorField> Up = film_.Up();
+//
+//     // employ simple coeff-based model
+//     const dimensionedScalar Cf("Cf", dimVelocity, Cf_);
+//
+//     return
+//     (
+//        - fam::Sp(Cf, U) + Cf*Up()     // surface contribution
+//     );
+// }
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace areaSurfaceFilmModels
diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
index b0e6b7054136866f3adf18ce0953ac63d5f09663..ba627b7040d9015abbc32270e5cc97bdd0ea2527 100644
--- a/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
+++ b/src/regionFaModels/liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,7 +60,7 @@ class laminar
     // Private Data
 
         //- Surface roughness coefficient
-        scalar Cf_;
+        //scalar Cf_;
 
 
     // Private Member Functions
@@ -90,16 +90,22 @@ public:
 
     // Member Functions
 
-        // Evolution
+        // Turbulence
+
+            //- Wall friction
+            tmp<faVectorMatrix> wallFriction(areaVectorField& U) const;
 
             //- Return the film turbulence viscosity
             virtual tmp<areaScalarField> mut() const;
 
-            //- Correct/update the model
-            virtual void correct();
-
             //- Return the source for the film momentum equation
             virtual tmp<faVectorMatrix> Su(areaVectorField& U) const;
+
+
+        // Evolution
+
+            //- Correct/update the model
+            virtual void correct();
 };