diff --git a/src/lagrangian/dsmc/boundaryMeasurements/boundaryMeasurements.C b/src/lagrangian/dsmc/boundaryMeasurements/boundaryMeasurements.C
index b296e631..6a649a69 100644
--- a/src/lagrangian/dsmc/boundaryMeasurements/boundaryMeasurements.C
+++ b/src/lagrangian/dsmc/boundaryMeasurements/boundaryMeasurements.C
@@ -325,20 +325,20 @@ void boundaryMeasurements::outputResults()
 {
     if(mesh_.time().outputTime())
     {
-        if(cloud_.boundaries().isAStickingPatch())
-        {
-            writenStuckParticles();
-        }
-
-        if(cloud_.boundaries().isAAbsorbingPatch())
-        {
-            writenAbsorbedParticles();
-        }
-
-        if(cloud_.boundaries().isAFieldPatch())
-        {
-            writePatchFields();
-        }
+//        if(cloud_.boundaries().isAStickingPatch())
+//        {
+//            writenStuckParticles();
+//        }
+
+//        if(cloud_.boundaries().isAAbsorbingPatch())
+//        {
+//            writenAbsorbedParticles();
+//        }
+
+//        if(cloud_.boundaries().isAFieldPatch())
+//        {
+//            writePatchFields();
+//        }
     }
 }
 
diff --git a/src/lagrangian/dsmc/clouds/dsmcCloud.C b/src/lagrangian/dsmc/clouds/dsmcCloud.C
index 48f6d513..2feea632 100644
--- a/src/lagrangian/dsmc/clouds/dsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/dsmcCloud.C
@@ -631,7 +631,14 @@ Foam::dsmcCloud::dsmcCloud
     ),
     collisionSelectionRemainder_(),
     constProps_(),
-    rndGen_(label(clock::getTime()) + 7183*Pstream::myProcNo()), // different seed every time simulation is started - needed for ensemble averaging!
+    rndGen_
+    (
+        particleProperties_.lookupOrDefault<label>
+        (
+            "seedNumber",
+            label(clock::getTime()) + 7183*Pstream::myProcNo()
+        )
+    ), 
     controllers_(t, mesh, *this),
     dynamicLoadBalancing_(t, mesh, *this),
     boundaryMeas_(mesh, *this, true),
@@ -744,7 +751,14 @@ Foam::dsmcCloud::dsmcCloud
     ),
     collisionSelectionRemainder_(),
     constProps_(),
-    rndGen_(label(clock::getTime()) + 1526*Pstream::myProcNo()),
+    rndGen_
+    (
+        particleProperties_.lookupOrDefault<label>
+        (
+            "seedNumber",
+            label(clock::getTime()) + 1526*Pstream::myProcNo()
+        )
+    ),
     controllers_(t, mesh),
     dynamicLoadBalancing_(t, mesh, *this),
     boundaryMeas_(mesh, *this),
@@ -1004,8 +1018,8 @@ void Foam::dsmcCloud::autoMap(const mapPolyMesh& mapper)
 
 Foam::label Foam::dsmcCloud::randomLabel
 (
-            const label valOne,
-            const label valTwo
+    const label valOne,
+    const label valTwo
 )
 {
     if (valOne == valTwo)
diff --git a/src/lagrangian/dsmc/macroscopicProperties/basic/dsmcField/dsmcField.C b/src/lagrangian/dsmc/macroscopicProperties/basic/dsmcField/dsmcField.C
index b91e2126..a5385b95 100644
--- a/src/lagrangian/dsmc/macroscopicProperties/basic/dsmcField/dsmcField.C
+++ b/src/lagrangian/dsmc/macroscopicProperties/basic/dsmcField/dsmcField.C
@@ -103,6 +103,7 @@ autoPtr<dsmcField> dsmcField::New
 dsmcField::~dsmcField()
 {}
 
+
 void dsmcField::updateTime()
 {
     time_++;
@@ -152,11 +153,13 @@ void dsmcField::updateBasicFieldProperties
     }
 }
 
+
 const fileName& dsmcField::casePath() const
 {
     return casePath_;
 }
 
+
 fileName& dsmcField::casePath()
 {
     return casePath_;
@@ -168,13 +171,13 @@ const fileName& dsmcField::timePath() const
     return timePath_;
 }
 
+
 fileName& dsmcField::timePath()
 {
     return timePath_;
 }
 
 
-
 } // End namespace Foam
 
 // ************************************************************************* //
diff --git a/src/lagrangian/dsmc/macroscopicProperties/derived/combined/dsmcVolFields/dsmcVolFields.C b/src/lagrangian/dsmc/macroscopicProperties/derived/combined/dsmcVolFields/dsmcVolFields.C
index 0349056b..ed3071d3 100644
--- a/src/lagrangian/dsmc/macroscopicProperties/derived/combined/dsmcVolFields/dsmcVolFields.C
+++ b/src/lagrangian/dsmc/macroscopicProperties/derived/combined/dsmcVolFields/dsmcVolFields.C
@@ -23,12 +23,14 @@ License
 
 Description
 
-Measures overall temperature, including vibrational temperature, for a single species gas
-or a gas mixture and writes the results to a volume scalar field that can be viewed in Paraview.
+Measures DSMC macroscopic fields for a single species or a gas mixture and
+writes the results to a volume field that can be viewed in Paraview.
 
-Translational, rotatational and vibrational temperature field will also be written automatically.
+Translational, rotatational and vibrational temperature fields will also be
+written automatically.
 
-Boundary fields are measured in conjunction with the boundaryMeasurements class and are also written.
+Boundary fields are measured in conjunction with the boundaryMeasurements class
+and are also written.
 
 \*---------------------------------------------------------------------------*/
 
@@ -64,7 +66,8 @@ void dsmcVolFields::calculateWallUnitVectors()
 
                 //- Wall tangential unit vector. Use the direction between the
                 // face centre and the first vertex in the list
-                t1_[patchi][facei] = fC[facei] - mesh_.points()[mesh_.faces()[pPatch.start() + facei][0]];
+                t1_[patchi][facei] = fC[facei]
+                    - mesh_.points()[mesh_.faces()[pPatch.start() + facei][0]];
                 t1_[patchi][facei] /= mag(t1_[patchi][facei]);
 
                 //- Other tangential unit vector.  Rescaling in case face is not
@@ -72,10 +75,6 @@ void dsmcVolFields::calculateWallUnitVectors()
                 t2_[patchi][facei] = n_[patchi][facei]^t1_[patchi][facei];
                 t2_[patchi][facei] /= mag(t2_[patchi][facei]);
             }
-
-            //Info << "n_[patchi]: " << n_[patchi] << endl;
-            //Info << "t1_[patchi]: " << t1_[patchi] << endl;
-            //Info << "t2_[patchi]: " << t2_[patchi] << endl;
         }
     }
 }
@@ -972,8 +971,8 @@ void dsmcVolFields::createField()
     mfpReferenceTemperature_ =
         propsDict_.lookupOrDefault<scalar>("mfpReferenceTemperature", 273.0);
 
-    averagingAcrossManyRuns_ =
-        propsDict_.lookupOrDefault<bool>("averagingAcrossManyRuns", false);
+    averagingAcrossManyRuns_ = false;
+//        propsDict_.lookupOrDefault<bool>("averagingAcrossManyRuns", false); TODO
 
     //- read in stored data from dictionary
     if (averagingAcrossManyRuns_)
@@ -981,7 +980,7 @@ void dsmcVolFields::createField()
         if (!time_.resetFieldsAtOutput())
         {
             Info<< "Averaging across many runs for field " << fieldName_
-                << " is enabled. Averaging data will be read from file."
+                << " is enabled. Sampled data will be read from file."
                 << endl;
             readIn();
         }
@@ -1036,8 +1035,6 @@ void dsmcVolFields::calculateField()
         }
         else
         {
-            //scalar timer = mesh_.time().elapsedCpuTime(); // TODO VINCENT
-
             forAllConstIter(dsmcCloud, cloud_, iter)
             {
                 const dsmcParcel& p = iter();
@@ -1045,12 +1042,14 @@ void dsmcVolFields::calculateField()
 
                 if (iD != -1 && p.isFree())
                 {
-                    const dsmcParcel::constantProperties& cP = cloud_.constProps(p.typeId());
+                    const dsmcParcel::constantProperties& cP =
+                        cloud_.constProps(p.typeId());
+                    const vector& Up = p.U();  
 
                     const label cell = p.cell();
                     const scalar nParticles = cloud_.nParticles(cell);
                     const scalar mass = cP.mass();
-                    const scalar massBySqMagU = mass*(p.U() & p.U());
+                    const scalar massBySqMagU = mass*(Up & Up);
                     const scalar rotationalDof = cP.rotationalDegreesOfFreedom();
                     const scalarList& electronicEnergies = cP.electronicEnergyList();
 
@@ -1071,7 +1070,7 @@ void dsmcVolFields::calculateField()
                     // cumulative linear kinetic energy of the DSMC parcels
                     dsmcLinKinEnCum_[cell] += massBySqMagU;
                     // cumulative momentum of the DSMC parcels
-                    dsmcMomentumCum_[cell] += mass*p.U();
+                    dsmcMomentumCum_[cell] += mass*Up;
                     // cumulative rotational energy of the DSMC parcels
                     dsmcRotEnCum_[cell] += p.ERot();
                     // cumulative rotational degrees of freedom of the DSMC
@@ -1092,31 +1091,31 @@ void dsmcVolFields::calculateField()
                     // cumulative mass
                     mCum_[cell] += mass*nParticles;
                     // cumulative momentum of the simulated real particles
-                    momentumCum_[cell] += mass*(p.U())*nParticles;
+                    momentumCum_[cell] += mass*(Up)*nParticles;
                     // cumulative linear kinetic energy of the simulated
                     // particles
                     linKinEnCum_[cell] += massBySqMagU*nParticles;
 
                     // cumulative counters for mass times squared velocity
                     // components of the DSMC parcels
-                    // notation: p.U = c = (u v w)
-                    dsmcMuuCum_[cell] += mass*sqr(p.U().x());
-                    dsmcMuvCum_[cell] += mass*p.U().x()*p.U().y();
-                    dsmcMuwCum_[cell] += mass*p.U().x()*p.U().z();
-                    dsmcMvvCum_[cell] += mass*sqr(p.U().y());
-                    dsmcMvwCum_[cell] += mass*p.U().y()*p.U().z();
-                    dsmcMwwCum_[cell] += mass*sqr(p.U().z());
+                    // notation: Up = c = (u v w)
+                    dsmcMuuCum_[cell] += mass*sqr(Up.x());
+                    dsmcMuvCum_[cell] += mass*Up.x()*Up.y();
+                    dsmcMuwCum_[cell] += mass*Up.x()*Up.z();
+                    dsmcMvvCum_[cell] += mass*sqr(Up.y());
+                    dsmcMvwCum_[cell] += mass*Up.y()*Up.z();
+                    dsmcMwwCum_[cell] += mass*sqr(Up.z());
 
                     dsmcMccCum_[cell] += massBySqMagU;
-                    dsmcMccuCum_[cell] += massBySqMagU*(p.U().x());
-                    dsmcMccvCum_[cell] += massBySqMagU*(p.U().y());
-                    dsmcMccwCum_[cell] += massBySqMagU*(p.U().z());
+                    dsmcMccuCum_[cell] += massBySqMagU*(Up.x());
+                    dsmcMccvCum_[cell] += massBySqMagU*(Up.y());
+                    dsmcMccwCum_[cell] += massBySqMagU*(Up.z());
 
                     // cumulative energy times velocity components of the DSMC
                     // parcels
-                    dsmcEuCum_[cell] += (p.ERot() + vibEn)*p.U().x(); // TODO missing Eel
-                    dsmcEvCum_[cell] += (p.ERot() + vibEn)*p.U().y(); // TODO missing Eel
-                    dsmcEwCum_[cell] += (p.ERot() + vibEn)*p.U().z(); // TODO missing Eel
+                    dsmcEuCum_[cell] += (p.ERot() + vibEn)*Up.x(); // TODO missing Eel
+                    dsmcEvCum_[cell] += (p.ERot() + vibEn)*Up.y(); // TODO missing Eel
+                    dsmcEwCum_[cell] += (p.ERot() + vibEn)*Up.z(); // TODO missing Eel
                     dsmcECum_[cell] += (p.ERot() + vibEn); // TODO missing Eel
 
                     if (rotationalDof > 0)
@@ -1162,8 +1161,6 @@ void dsmcVolFields::calculateField()
                 }
             }
 
-            //Info<< "fields myCal" << tab << mesh_.time().elapsedCpuTime() - timer << " s" << endl; // TODO VINCENT
-
             //- Obtain collision quality measurements and mixture translational
             //  temperature
             forAll(cloud_.cellPropMeasurements().collisionSeparation(), cell)
@@ -1185,9 +1182,8 @@ void dsmcVolFields::calculateField()
 
                     rhoN_[cell] = rhoNMean;
                     rhoM_[cell] = rhoMMean;
-
-                    UMean_[cell] = momentumCum_[cell]
-                        /mCum_[cell];
+                    
+                    UMean_[cell] = momentumCum_[cell]/mCum_[cell];
 
                     const scalar linearKEMean = 0.5
                         *linKinEnCum_[cell]
@@ -1450,7 +1446,6 @@ void dsmcVolFields::calculateField()
                     if
                     (
                          dsmcNWithRotDofCum_[cell] > VSMALL
-                      && dsmcNCum_[cell] > VSMALL
                       && dsmcNSpeciesCum_[i][cell] > SMALL
                     )
                     {