diff --git a/applications/solvers/rcfdemSolverHeattransfer/TEqImp.H b/applications/solvers/rcfdemSolverHeattransfer/TEqImp.H
new file mode 100644
index 00000000..216c4735
--- /dev/null
+++ b/applications/solvers/rcfdemSolverHeattransfer/TEqImp.H
@@ -0,0 +1,49 @@
+	volScalarField rhoeps = rhoRec*voidfractionRec;
+	
+	particleCloud.energyContributions(Qsource);
+
+        thDiff=particleCloud.thermCondM().thermDiff();
+	
+	dimensionedScalar Cp = particleCloud.thermCondM().Cp()*718/1005; // need cV, hence the factor
+
+	K = 0.5*magSqr(URec);
+
+	addSource = fvc::ddt(rhoeps, K) + fvc::div(phiRec, K) +
+	fvc::div
+	(
+            fvc::absolute(phiRec/fvc::interpolate(rhoRec), voidfractionRec*URec), pRec
+	);
+
+	particleCloud.energyCoefficients(QCoeff);
+	
+	fvScalarMatrix TEqn = 
+	(
+		fvm::ddt(rhoeps, T)
+		+ fvm::div(phiRec, T)
+		+ addSource/Cp
+		- fvm::laplacian(rhoeps*thDiff, T)
+		- Qsource/Cp
+		- fvm::Sp(QCoeff/Cp, T)
+		==
+		fvOptions(rhoeps, T)	// no fvOptions support yet
+	);
+	
+	//TEqn.relax(relaxCoeff);
+	
+	fvOptions.constrain(TEqn);	// no fvOptions support yet
+	
+	TEqn.solve();
+
+        particleCloud.clockM().start(31,"postFlow");
+	counter++;
+	if( (counter - couplingSubStep) % 10 == 0)
+	    particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");
+
+	
+	// change couplingSubStep during course of simulation
+//	if( counter == 16000 )
+//	    couplingSubStep++;
+	
+//	if (doCouple)
+//	    particleCloud.postFlow();
diff --git a/applications/solvers/rcfdemSolverHeattransfer/updateFields.H b/applications/solvers/rcfdemSolverHeattransfer/updateFields.H
new file mode 100644
index 00000000..5e9a8ab1
--- /dev/null
+++ b/applications/solvers/rcfdemSolverHeattransfer/updateFields.H
@@ -0,0 +1,14 @@
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolScalarField("rho",rhoRec);
+recurrenceBase.recM().exportVolScalarField("p",pRec);
+recurrenceBase.recM().exportVolVectorField("Us",UsRec);
+recurrenceBase.recM().exportVolVectorField("U",URec);
+recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);
+
+{
+    volScalarField& NuField(const_cast<volScalarField&>(mesh.lookupObject<volScalarField> ("NuField")));
+    recurrenceBase.recM().exportVolScalarField("NuField",NuField);
+}
+
+
+#include "updateRho.H"
\ No newline at end of file
diff --git a/applications/solvers/rcfdemSolverHeattransfer/updateRho.H b/applications/solvers/rcfdemSolverHeattransfer/updateRho.H
new file mode 100644
index 00000000..9830513f
--- /dev/null
+++ b/applications/solvers/rcfdemSolverHeattransfer/updateRho.H
@@ -0,0 +1,29 @@
+// case 1
+
+forAll(rhoRec,cellI)
+{
+    if (mesh.C()[cellI].z() < 0.00228)
+        rhoRec[cellI] = 1.18+(1.085-1.18)*Foam::exp(-0.065*runTime.timeOutputValue());
+    else if (mesh.C()[cellI].z() < 0.00456)
+        rhoRec[cellI] = 1.18+(1.01-1.18)*Foam::exp(-0.05*runTime.timeOutputValue());
+    else if (mesh.C()[cellI].z() < 0.00684)
+        rhoRec[cellI] = 1.18+(0.98-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());
+    else
+        rhoRec[cellI] = 1.18+(0.955-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());   
+}
+
+
+// case 2
+/*
+forAll(rhoRec,cellI)
+{
+    if (mesh.C()[cellI].z() < 0.00228)
+        rhoRec[cellI] = 1.18+(1.115-1.18)*Foam::exp(-0.065*runTime.timeOutputValue());
+    else if (mesh.C()[cellI].z() < 0.00456)
+        rhoRec[cellI] = 1.18+(1.04-1.18)*Foam::exp(-0.05*runTime.timeOutputValue());
+    else if (mesh.C()[cellI].z() < 0.00684)
+        rhoRec[cellI] = 1.18+(1.005-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());
+    else
+        rhoRec[cellI] = 1.18+(0.96-1.18)*Foam::exp(-0.0425*runTime.timeOutputValue());   
+}
+*/
\ No newline at end of file