diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/Make/files b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/Make/files
new file mode 100644
index 00000000..57c616af
--- /dev/null
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/Make/files
@@ -0,0 +1,3 @@
+cfdemSolverPisoTemp.C
+
+EXE = $(CFDEM_APP_DIR)/cfdemSolverPisoTemp
diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/Make/options b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/Make/options
new file mode 100644
index 00000000..68aecd17
--- /dev/null
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/Make/options
@@ -0,0 +1,19 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/TEqn.H b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/TEqn.H
new file mode 100644
index 00000000..4ee4e08e
--- /dev/null
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/TEqn.H
@@ -0,0 +1,15 @@
+        // get scalar source from DEM        
+        particleCloud.energyContributions(Qsource);
+        Qsource.correctBoundaryConditions();
+
+        // solve temperature transport equation assuming const. density and heat capacity
+        fvScalarMatrix TEqn
+        (
+           fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
+         + fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
+         - fvm::laplacian(voidfraction*particleCloud.thermCondM().thermDiff(), T)
+         ==
+           Qsource/(rho*particleCloud.energyM(0).Cp())
+        );
+        TEqn.relax();
+        TEqn.solve();
diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/cfdemSolverPisoTemp.C b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/cfdemSolverPisoTemp.C
new file mode 100644
index 00000000..0096290d
--- /dev/null
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/cfdemSolverPisoTemp.C
@@ -0,0 +1,201 @@
+/*---------------------------------------------------------------------------*\
+    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
+    cfdemSolverPisoScalar
+
+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 pisoFoam in OpenFOAM(R) 1.6,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+
+#include "cfdemCloudEnergy.H"
+#include "implicitCouple.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+#include "energyModel.H"
+#include "thermCondModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudEnergy particleCloud(mesh);
+    #include "checkModelType.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "readPISOControls.H"
+        #include "CourantNo.H"
+
+        // do particle stuff
+        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();
+
+
+        #include "solverDebugInfo.H"
+
+        #include "TEqn.H"
+
+        if(particleCloud.solveFlow())
+        {
+            // Pressure-velocity PISO corrector
+            {
+                // Momentum predictor
+                fvVectorMatrix UEqn
+                (
+                    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+                  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
+//                + turbulence->divDevReff(U)
+                  + particleCloud.divVoidfractionTau(U, voidfraction)
+                  ==
+                  - fvm::Sp(Ksl/rho,U)
+                );
+
+                UEqn.relax();
+                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
+                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
+                else if (momentumPredictor)
+                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+
+                // --- PISO loop
+
+                //for (int corr=0; corr<nCorr; corr++)
+                int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
+
+                for (int corr=0; corr<nCorrSoph; corr++)
+                {
+                    volScalarField rUA = 1.0/UEqn.A();
+
+                    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
+                    volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+                    surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
+
+                    U = rUA*UEqn.H();
+
+                    #ifdef version23
+                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+                        + rUAfvoidfraction*fvc::ddtCorr(U, phi);
+                    #else
+                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
+                    #endif
+                    surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
+                    surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
+
+                    if (modelType=="A")
+                        rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
+
+                    // Non-orthogonal pressure corrector loop
+                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                    {
+                        // Pressure corrector
+                        fvScalarMatrix pEqn
+                        (
+                            fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
+                        );
+                        pEqn.setReference(pRefCell, pRefValue);
+
+                        if
+                        (
+                            corr == nCorr-1
+                         && nonOrth == nNonOrthCorr
+                        )
+                        {
+                            pEqn.solve(mesh.solver("pFinal"));
+                        }
+                        else
+                        {
+                            pEqn.solve();
+                        }
+
+                        if (nonOrth == nNonOrthCorr)
+                        {
+                            phiGes -= pEqn.flux();
+                            phi = phiGes;
+                        }
+
+                    } // end non-orthogonal corrector loop
+
+                    #include "continuityErrorPhiPU.H"
+
+                    if (modelType=="B" || modelType=="Bfull")
+                        U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
+                    else
+                        U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
+
+                    U.correctBoundaryConditions();
+
+                } // end piso loop
+            }
+
+            turbulence->correct();
+        }// end solveFlow
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/createFields.H b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/createFields.H
new file mode 100644
index 00000000..9b7ba85d
--- /dev/null
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/createFields.H
@@ -0,0 +1,156 @@
+    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
+    );
+
+//========================
+// drag law modelling
+//========================
+
+    Info<< "\nReading momentum exchange field Ksl\n" << endl;
+    volScalarField Ksl
+    (
+        IOobject
+        (
+            "Ksl",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+        //dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.0)
+    );
+
+    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 density field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// scalar field modelling
+//========================
+
+    Info<< "\nCreating temperature field\n" << endl;
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating fluid-particle heat flux field\n" << endl;
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+
+//#   include "createPhi.H"
+#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
+
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );