diff --git a/applications/solvers/cfdemSolverPimple/Make/files b/applications/solvers/cfdemSolverPimple/Make/files
new file mode 100644
index 00000000..9fd6788a
--- /dev/null
+++ b/applications/solvers/cfdemSolverPimple/Make/files
@@ -0,0 +1,3 @@
+cfdemSolverPimple.C
+
+EXE=$(CFDEM_APP_DIR)/cfdemSolverPimpleRadl
diff --git a/applications/solvers/cfdemSolverPimple/Make/options b/applications/solvers/cfdemSolverPimple/Make/options
new file mode 100644
index 00000000..a6db87b6
--- /dev/null
+++ b/applications/solvers/cfdemSolverPimple/Make/options
@@ -0,0 +1,25 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -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)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
diff --git a/applications/solvers/cfdemSolverPimple/UEqn.H b/applications/solvers/cfdemSolverPimple/UEqn.H
new file mode 100644
index 00000000..5d615d86
--- /dev/null
+++ b/applications/solvers/cfdemSolverPimple/UEqn.H
@@ -0,0 +1,40 @@
+particleCloud.otherForces(fOther);
+
+tmp<fvVectorMatrix> tUEqn
+(
+    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U) 
+  - fOther/rho
+  ==
+    fvOptions(U)
+  - fvm::Sp(Ksl/rho,U)
+);
+
+fvVectorMatrix& UEqn = tUEqn.ref();
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+volScalarField rUA = 1.0/UEqn.A();
+
+surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+
+surfaceScalarField voidfractionf = fvc::interpolate(voidfraction);
+
+surfaceScalarField phicForces
+(
+	fvc::interpolate(rUA*(Ksl*Us)/rho) & mesh.Sf()
+);
+
+if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve(UEqn == fvc::reconstruct(phicForces/rUAf - fvc::snGrad(p)*mesh.magSf()));
+    fvOptions.correct(U);
+}
+else if (pimple.momentumPredictor())
+{
+    solve(UEqn == fvc::reconstruct(phicForces/rUAf - fvc::snGrad(p)*voidfractionf*mesh.magSf()));
+    fvOptions.correct(U);
+}
+
diff --git a/applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C b/applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C
new file mode 100644
index 00000000..91572c5d
--- /dev/null
+++ b/applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C
@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-     DCS Computing GmbH,Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling 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.
+
+    CFDEMcoupling 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 CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverPimple
+
+Description
+    Transient solver for incompressible flow.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver pimpleFoam in OpenFOAM(R) 6.0,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+
+#include "cfdemCloud.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloud particleCloud(mesh);
+    #include "checkModelType.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.loop())
+    {
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "CourantNo.H"
+
+        // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        if(particleCloud.solveFlow())
+        {
+            // Pressure-velocity PIMPLE corrector
+            while (pimple.loop())
+            {
+                // Momentum predictor
+                 #include "UEqn.H"
+
+                // --- Inner PIMPLE loop
+
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+            }
+
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        runTime.write();
+        
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/cfdemSolverPimple/createFields.H b/applications/solvers/cfdemSolverPimple/createFields.H
new file mode 100644
index 00000000..e8df010a
--- /dev/null
+++ b/applications/solvers/cfdemSolverPimple/createFields.H
@@ -0,0 +1,155 @@
+//===============================
+// Fluid Fields
+//===============================      
+   
+   Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading physical velocity field U" << endl;
+    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================
+// particle interaction modelling
+//===============================
+
+    Info<< "\nReading momentum exchange field Ksl\n" << endl;
+    volScalarField Ksl
+    (
+        IOobject
+        (
+            "Ksl",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
+    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating dummy density field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================
+
+#ifndef createPhi_H
+#define createPhi_H
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+         mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(U*voidfraction) & mesh.Sf()
+);
+#endif
+
+ Info<< "Generating interstitial flux field phiByVoidfraction (this is the INTERSTITIAL field!)\n" << endl;
+ surfaceScalarField phiByVoidfraction
+ (
+     IOobject
+     (
+        "phiByVoidfraction",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+     ),
+     linearInterpolate(U) & mesh.Sf()
+ );
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+#include "createMRF.H"
diff --git a/applications/solvers/cfdemSolverPimple/pEqn.H b/applications/solvers/cfdemSolverPimple/pEqn.H
new file mode 100644
index 00000000..3206f449
--- /dev/null
+++ b/applications/solvers/cfdemSolverPimple/pEqn.H
@@ -0,0 +1,59 @@
+volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+
+surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
+
+volVectorField HbyA("HbyA", U);
+
+HbyA = rUA*UEqn.H();
+
+phi = voidfractionf*phiByVoidfraction;
+
+surfaceScalarField phiHbyA
+(
+	"phiHbyA",
+	(
+	    (fvc::interpolate(HbyA) & mesh.Sf() )
+	    + phicForces                               //explicit contribution
+	    + rUAfvoidfraction*fvc::ddtCorr(U, phiByVoidfraction) //correction
+	)
+);
+
+if (modelType=="A")
+    rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
+    
+// Update the fixedFluxPressure BCs to ensure flux consistency
+if (modelType=="A")
+{
+	volScalarField rUsed = rUA*voidfraction;
+	constrainPressure(p, U, phiHbyA, rUsed,MRF);
+}
+else constrainPressure(p, U, phiHbyA, rUA,MRF);
+
+while (pimple.correctNonOrthogonal())
+{
+    // Pressure corrector
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rUAvoidfraction, p) == fvc::div(voidfractionf*phiHbyA) + particleCloud.ddtVoidfraction()
+    );
+    pEqn.setReference(pRefCell, pRefValue);
+    
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+    
+    if (pimple.finalNonOrthogonalIter())
+    {
+	    phiByVoidfraction = phiHbyA - pEqn.flux()/voidfractionf;
+	    phi = voidfractionf*phiByVoidfraction;
+
+	    #include "continuityErrorPhiPU.H"
+
+	    // Explicitly relax pressure for momentum corrector
+	    p.relax();
+
+	    U = fvc::reconstruct(phiHbyA)
+	      - rUA*fvc::reconstruct(pEqn.flux()/voidfractionf/rUAf);
+
+	    U.correctBoundaryConditions();
+	    fvOptions.correct(U);
+    }
+}