diff --git a/applications/utilities/preProcessing/setTurbulenceFields/Make/files b/applications/utilities/preProcessing/setTurbulenceFields/Make/files
new file mode 100644
index 0000000000..d4a1527323
--- /dev/null
+++ b/applications/utilities/preProcessing/setTurbulenceFields/Make/files
@@ -0,0 +1,3 @@
+setTurbulenceFields.C
+
+EXE = $(FOAM_APPBIN)/setTurbulenceFields
diff --git a/applications/utilities/preProcessing/setTurbulenceFields/Make/options b/applications/utilities/preProcessing/setTurbulenceFields/Make/options
new file mode 100644
index 0000000000..a887778b7b
--- /dev/null
+++ b/applications/utilities/preProcessing/setTurbulenceFields/Make/options
@@ -0,0 +1,22 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfluidThermophysicalModels \
+    -lincompressibleTransportModels \
+    -lcompressibleTransportModels \
+    -lgenericPatchFields
diff --git a/applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C b/applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C
new file mode 100644
index 0000000000..6e69f83457
--- /dev/null
+++ b/applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C
@@ -0,0 +1,531 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 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/>.
+
+Application
+    setTurbulenceFields
+
+Group
+    grpPreProcessingUtilities
+
+Description
+    Initialises turbulence fields according to
+    various empirical governing equations.
+
+    References:
+    \verbatim
+        Initialisation method (tag:M):
+            Manceau, R. (n.d.).
+            A two-step automatic initialization
+            procedure for RANS computations.
+            (Unpublished).
+
+        Previous-version of the initialisation model (tag:LM):
+            Lardeau, S., & Manceau, R. (2014).
+            Computations of complex flow configurations using
+            a modified elliptic-blending Reynolds-stress model.
+            10th International ERCOFTAC Symposium on Engineering
+            Turbulence Modelling and Measurements. Marbella, Spain.
+            https://hal.archives-ouvertes.fr/hal-01051799
+    \endverbatim
+
+Usage
+    Minimal example by using \c system/setTurbulenceFieldsDict:
+    \verbatim
+        // Mandatory entries
+        uRef            <scalar>;
+
+        // Optional entries
+        initialiseU     <bool>;
+        initialiseEpsilon <bool>;
+        initialiseK     <bool>;
+        initialiseOmega <bool>;
+        initialiseR     <bool>;
+        writeF          <bool>;
+
+        kappa           <scalar>;
+        Cmu             <scalar>;
+        dPlusRef        <scalar>;
+
+        f               <word>;
+        U               <word>;
+        epsilon         <word>;
+        k               <word>;
+        omega           <word>;
+        R               <word>;
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                    | Type | Reqd | Deflt
+      uRef      | Reference speed              | scalar | yes  | -
+      initialiseU | Flag to initialise U         | bool | no   | false
+      initialiseEpsilon | Flag to initialise epsilon    | bool | no   | false
+      initialiseK | Flag to initialise k         | bool | no   | false
+      initialiseOmega | Flag to initialise omega | bool | no   | false
+      initialiseR | Flag to initialise R         | bool | no   | false
+      writeF | Flag to write elliptic-blending field, f | bool | no   | false
+      kappa  | von Karman constant             | scalar | no   | 0.41
+      Cmu    | Empirical constant              | scalar | no   | 0.09
+      dPlusRef | Reference dPlus               | scalar | no   | 15
+      f      | Name of operand f field         | word   | no   | f
+      U      | Name of operand U field         | word   | no   | U
+      epsilon  | Name of operand epsilon field | word   | no   | epsilon
+      k      | Name of operand k field         | word   | no   | k
+      omega      | Name of operand omega field | word   | no   | omega
+      R      | Name of operand R field         | word   | no   | R
+    \endtable
+
+Note
+  - Time that the utility applies to is determined by the
+    \c controlDict.startFrom and \c controlDict.startTime entries.
+  - The utility modifies near-wall fields, hence
+    can be more effective for low-Re mesh cases.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "turbulentFluidThermoModel.H"
+#include "wallFvPatch.H"
+#include "fixedValueFvPatchFields.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void InfoField(const word& fldName)
+{
+    Info<< "Writing field: " << fldName << nl << endl;
+}
+
+
+IOobject createIOobject
+(
+    const fvMesh& mesh,
+    const word& name,
+    IOobject::readOption rOpt = IOobject::READ_IF_PRESENT
+)
+{
+    return
+        IOobject
+        (
+            name,
+            mesh.time().timeName(),
+            mesh,
+            rOpt,
+            IOobject::NO_WRITE,
+            false // do not register
+        );
+}
+
+
+tmp<volScalarField> nu
+(
+    const fvMesh& mesh,
+    const volVectorField& U
+)
+{
+    const Time& runTime = mesh.time();
+
+    if
+    (
+        IOobject
+        (
+            basicThermo::dictName,
+            runTime.constant(),
+            mesh
+        ).typeHeaderOk<IOdictionary>(true)
+    )
+    {
+        // Compressible
+        autoPtr<fluidThermo> pThermo(fluidThermo::New(mesh));
+        fluidThermo& thermo = pThermo();
+        volScalarField rho(thermo.rho());
+
+        // Create phi
+        #include "compressibleCreatePhi.H"
+
+        autoPtr<compressible::turbulenceModel> turbulence
+        (
+            compressible::turbulenceModel::New
+            (
+                rho,
+                U,
+                phi,
+                thermo
+            )
+        );
+
+        turbulence->validate();
+
+        return tmp<volScalarField>::New(turbulence->nu());
+    }
+    else
+    {
+        // Incompressible
+
+        // Create phi
+        #include "createPhi.H"
+
+        singlePhaseTransportModel laminarTransport(U, phi);
+
+        autoPtr<incompressible::turbulenceModel> turbulence
+        (
+            incompressible::turbulenceModel::New(U, phi, laminarTransport)
+        );
+
+        turbulence->validate();
+
+        return tmp<volScalarField>::New(turbulence->nu());
+    }
+}
+
+
+// Calculate elliptic blending function
+// between near-wall and weakly-inhomogeneous regions
+void calcF
+(
+    const volScalarField& L,
+    volScalarField& f
+)
+{
+    tmp<volScalarField> tinvLsqr = scalar(1)/sqr(L);
+    const volScalarField& invLsqr = tinvLsqr.cref();
+
+    // (M:Eq. 6)
+    tmp<fvScalarMatrix> fEqn
+    (
+        fvm::Sp(invLsqr, f)
+      - fvm::laplacian(f)
+      ==
+        invLsqr
+    );
+
+    tinvLsqr.clear();
+
+    fEqn.ref().relax();
+    solve(fEqn);
+
+    // (M:p. 2)
+    const dimensioned<scalarMinMax> fMinMax
+    (
+        dimless,
+        scalarMinMax(Zero, scalar(1) - Foam::exp(-scalar(400)/scalar(50)))
+    );
+
+    f.clip(fMinMax);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Sets initial turbulence fields based on"
+        " various empirical equations"
+    );
+    argList::noFunctionObjects();
+    argList::addOption("dict", "file", "Alternative setTurbulenceFieldsDict");
+
+    #include "addRegionOption.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    const word dictName("setTurbulenceFieldsDict");
+    #include "setSystemRunTimeDictionaryIO.H"
+    Info<< "Reading " << dictIO.name() << nl << endl;
+    IOdictionary dict(dictIO);
+
+    #include "createNamedMesh.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    IOstream::defaultPrecision(15);
+
+
+    // Dictionary input (M:p. 2)
+
+    const scalar uRef = dict.getCheck<scalar>("uRef", scalarMinMax::ge(SMALL));
+    const dimensionedScalar uTau(dimVelocity, 0.05*uRef);
+    const scalar kappa =
+        dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(SMALL));
+    const scalar Cmu =
+        dict.getCheckOrDefault<scalar>("Cmu", 0.09, scalarMinMax::ge(SMALL));
+    const scalar dPlusRef =
+        dict.getCheckOrDefault<scalar>("dPlusRef", 15, scalarMinMax::ge(SMALL));
+    const word fName = dict.getOrDefault<word>("f", "f");
+    const word UName = dict.getOrDefault<word>("U", "U");
+    const word epsilonName = dict.getOrDefault<word>("epsilon", "epsilon");
+    const word kName = dict.getOrDefault<word>("k", "k");
+    const word omegaName = dict.getOrDefault<word>("omega", "omega");
+    const word RName = dict.getOrDefault<word>("R", "R");
+    const bool initU = dict.getOrDefault<bool>("initialiseU", false);
+    const bool initEpsilon =
+        dict.getOrDefault<bool>("initialiseEpsilon", false);
+    const bool initK = dict.getOrDefault<bool>("initialiseK", false);
+    const bool initOmega = dict.getOrDefault<bool>("initialiseOmega", false);
+    const bool initR = dict.getOrDefault<bool>("initialiseR", false);
+    const bool writeF = dict.getOrDefault<bool>("writeF", false);
+
+
+    // Start initialising the operand fields
+
+    // Read operand fields
+
+    auto tU = tmp<volVectorField>::New
+    (
+        createIOobject(mesh, UName, IOobject::MUST_READ),
+        mesh
+    );
+
+    // Infer the initial BCs from the velocity
+    const wordList bcTypes
+    (
+        tU.cref().boundaryField().size(),
+        fixedValueFvPatchScalarField::typeName
+    );
+
+    tmp<volScalarField> tepsilon;
+    tmp<volScalarField> tk;
+    tmp<volScalarField> tomega;
+    tmp<volSymmTensorField> tR;
+
+    if (initEpsilon)
+    {
+        tepsilon = tmp<volScalarField>::New
+        (
+            createIOobject(mesh, epsilonName),
+            mesh,
+            dimensionedScalar(sqr(dimLength)/pow3(dimTime), SMALL),
+            bcTypes
+        );
+    }
+
+    if (initK)
+    {
+        tk = tmp<volScalarField>::New
+        (
+            createIOobject(mesh, kName),
+            mesh,
+            dimensionedScalar(sqr(dimLength/dimTime), SMALL),
+            bcTypes
+        );
+    }
+
+    if (initOmega)
+    {
+        tomega = tmp<volScalarField>::New
+        (
+            createIOobject(mesh, omegaName),
+            mesh,
+            dimensionedScalar(dimless/dimTime, SMALL),
+            bcTypes
+        );
+    }
+
+    if (initR)
+    {
+        tR = tmp<volSymmTensorField>::New
+        (
+            createIOobject(mesh, RName),
+            mesh,
+            dimensionedSymmTensor(sqr(dimLength/dimTime), Zero),
+            bcTypes
+        );
+    }
+
+
+    // Create elliptic blending factor field
+
+    volScalarField f
+    (
+        IOobject
+        (
+            fName,
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false  // do not register
+        ),
+        mesh,
+        dimensionedScalar(dimless, scalar(1)),
+        fixedValueFvPatchField<scalar>::typeName
+    );
+
+    for (fvPatchScalarField& pfld : f.boundaryFieldRef())
+    {
+        if (isA<wallFvPatch>(pfld.patch()))
+        {
+            pfld == Zero;
+        }
+    }
+
+
+    // Create auxillary fields for the initialisation
+
+    tmp<volScalarField> tnu = nu(mesh, tU.cref());
+    const volScalarField& nu = tnu.cref();
+
+    // (M:p. 2)
+    tmp<volScalarField> tL = scalar(50)*nu/uTau;
+    const volScalarField& L = tL.cref();
+
+    calcF(L, f);
+
+    // (M:Eq. 8)
+    tmp<volScalarField> td = -tL*Foam::log(scalar(1) - f);
+    const volScalarField& d = td.cref();
+
+    // (M:p. 2)
+    const volScalarField dPlus(d*uTau/nu);
+
+    // (M:Eq. 11)
+    const volScalarField epsilon
+    (
+        pow4(uTau)/(kappa*nu*max(dPlus, dPlusRef))
+    );
+
+    // (M:Eq. 13)
+    const volScalarField fk(Foam::exp(-dPlus/scalar(25)));
+
+    // (M:Eq. 12)
+    const volScalarField k
+    (
+        (epsilon*sqr(td)*sqr(fk))/(2*nu)
+        + sqr(uTau)*sqr(scalar(1) - fk)/Foam::sqrt(Cmu)
+    );
+
+
+    // Initialise operand fields
+
+    if (initU)
+    {
+        volVectorField& U = tU.ref();
+
+        // Reichardt’s law (M:Eq. 10)
+        const scalar C = 7.8;
+        const scalar B1 = 11;
+        const scalar B2 = 3;
+        const volScalarField fRei
+        (
+            Foam::log(scalar(1) + kappa*dPlus)/kappa
+          + C*
+            (
+                scalar(1)
+              - Foam::exp(-dPlus/B1)
+              - dPlus/B1*Foam::exp(-dPlus/B2)
+            )
+        );
+
+        // (M:Eq. 9)
+        const dimensionedScalar maxU(dimVelocity, SMALL);
+        U *= min(scalar(1), fRei*uTau/max(mag(U), maxU));
+    }
+
+    if (tepsilon.valid())
+    {
+        tepsilon.ref() = epsilon;
+    }
+
+    if (tk.valid())
+    {
+        tk.ref() = k;
+    }
+
+    if (tomega.valid())
+    {
+        const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
+        tomega.ref() = Cmu*epsilon/(k + k0);
+    }
+
+    if (tR.valid())
+    {
+        volSymmTensorField& R = tR.ref();
+
+        // (M:Eq. 3)
+        const volSphericalTensorField Rdiag(k*twoThirdsI);
+        forAll(R, celli)
+        {
+            R[celli] = Rdiag[celli];
+        }
+    }
+
+
+    // Write operand fields
+
+    Info<< endl;
+
+    if (initU)
+    {
+        InfoField(tU->name());
+        tU->write();
+    }
+
+    if (tepsilon.valid())
+    {
+        InfoField(tepsilon->name());
+        tepsilon->write();
+    }
+
+    if (tk.valid())
+    {
+        InfoField(tk->name());
+        tk->write();
+    }
+
+    if (tomega.valid())
+    {
+        InfoField(tomega->name());
+        tomega->write();
+    }
+
+    if (tR.valid())
+    {
+        InfoField(tR->name());
+        tR->write();
+    }
+
+    if (writeF)
+    {
+        InfoField(f.name());
+        f.write();
+    }
+
+
+    Info<< nl;
+    runTime.printExecutionTime(Info);
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //