diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
index 2e4dff7a..cc6ad959 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
@@ -44,6 +44,7 @@ heatTransferGunn::heatTransferGunn
     energyModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
+    verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
     QPartFluidName_(propsDict_.lookupOrDefault<word>("QPartFluidName","QPartFluid")),
     QPartFluid_
     (   IOobject
@@ -81,6 +82,30 @@ heatTransferGunn::heatTransferGunn
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
     ),
+    ReField_
+    (   IOobject
+        (
+            "ReField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+    ),
+    NuField_
+    (   IOobject
+        (
+            "NuField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+    ),
     partRefTemp_("partRefTemp", dimensionSet(0,0,0,1,0,0,0), 0.0),
     calcPartTempField_(propsDict_.lookupOrDefault<bool>("calcPartTempField",false)),
     tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
@@ -95,7 +120,9 @@ heatTransferGunn::heatTransferGunn
     partTempName_(propsDict_.lookup("partTempName")),
     partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookup("partHeatFluxName")),
-    partHeatFlux_(NULL)
+    partHeatFlux_(NULL),
+    partRe_(NULL),
+    partNu_(NULL)
 {
      allocateMyArrays();
 
@@ -114,6 +141,13 @@ heatTransferGunn::heatTransferGunn
 	partRelTempField_.write();
 	Info <<  "Particle temperature field activated." << endl;
     }
+    if (verbose_)
+    {
+        ReField_.writeOpt() = IOobject::AUTO_WRITE;
+	NuField_.writeOpt() = IOobject::AUTO_WRITE;
+	ReField_.write();
+	NuField_.write();
+    }
 }
 
 
@@ -123,6 +157,8 @@ heatTransferGunn::~heatTransferGunn()
 {
     delete partTemp_;
     delete partHeatFlux_;
+    delete partRe_;
+    delete partNu_;
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -132,6 +168,12 @@ void heatTransferGunn::allocateMyArrays() const
     double initVal=0.0;
     particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
     particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
+    
+    if(verbose_)
+    {
+        particleCloud_.dataExchangeM().allocateArray(partRe_,initVal,1);
+        particleCloud_.dataExchangeM().allocateArray(partNu_,initVal,1);
+    }
 }
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
@@ -209,6 +251,9 @@ void heatTransferGunn::calcEnergyContribution()
                     Ufluid = U_[cellI];
                     Tfluid = tempField_[cellI];
                 }
+                
+                if (voidfraction < 0.01)
+		    voidfraction = 0.01;
 
                 // calc relative velocity
                 Us = particleCloud_.velocity(index);
@@ -228,6 +273,12 @@ void heatTransferGunn::calcEnergyContribution()
                 heatFlux(index, h, As, Tfluid);
                 heatFluxCoeff(index, h, As);
 		
+		if(verbose_)
+		{
+		    partRe_[index][0] = Rep;
+		    partNu_[index][0] = Nup;
+		}
+		
 		if(particleCloud_.verbose() && index >=0 && index <2)
                 {
                     Info << "partHeatFlux = " << partHeatFlux_[index][0] << endl;
@@ -254,6 +305,32 @@ void heatTransferGunn::calcEnergyContribution()
 
     QPartFluid_.primitiveFieldRef() /= -QPartFluid_.mesh().V();
 
+    if(verbose_)
+    {
+        ReField_.primitiveFieldRef() = 0.0;
+	NuField_.primitiveFieldRef() = 0.0;
+	particleCloud_.averagingM().resetWeightFields();
+        particleCloud_.averagingM().setScalarAverage
+        (
+            ReField_,
+            partRe_,
+            particleCloud_.particleWeights(),
+	    particleCloud_.averagingM().UsWeightField(),
+            NULL
+        );
+	particleCloud_.averagingM().resetWeightFields();
+        particleCloud_.averagingM().setScalarAverage
+        (
+            NuField_,
+            partNu_,
+            particleCloud_.particleWeights(),
+	    particleCloud_.averagingM().UsWeightField(),
+            NULL
+        );
+	ReField_.primitiveFieldRef() /= ReField_.mesh().V();
+	NuField_.primitiveFieldRef() /= NuField_.mesh().V();
+    }
+    
     // limit source term
     forAll(QPartFluid_,cellI)
     {
@@ -261,6 +338,7 @@ void heatTransferGunn::calcEnergyContribution()
 
         if(mag(EuFieldInCell) > maxSource_ )
         {
+	     Pout << "limiting source term\n"  << endl  ;
              QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
         }
     }
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
index 3437d125..4862e5ed 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
@@ -47,6 +47,8 @@ protected:
     
     bool interpolation_;
     
+    bool verbose_;
+    
     word QPartFluidName_;
         
     volScalarField QPartFluid_;
@@ -55,6 +57,10 @@ protected:
     
     volScalarField partRelTempField_;
     
+    volScalarField ReField_;
+    
+    volScalarField NuField_;
+    
     dimensionedScalar partRefTemp_;
     
     bool calcPartTempField_;
@@ -84,6 +90,10 @@ protected:
     word partHeatFluxName_;
 
     mutable double **partHeatFlux_;              // Lagrangian array
+    
+    mutable double **partRe_;
+    
+    mutable double **partNu_;
 
     void allocateMyArrays() const;