diff --git a/source/Applications/bChannelVA/Averaging.cpp b/source/Applications/bChannelVA/Averaging.cpp
index 3d99f3711e1eea26cfb1d14c3b73a25437887cd2..92c34aad4f7149ec27d4a1bdd4e670edfe90fab9 100644
--- a/source/Applications/bChannelVA/Averaging.cpp
+++ b/source/Applications/bChannelVA/Averaging.cpp
@@ -359,6 +359,11 @@ void Averaging::writeVaSumMatrixToImageFile(std::string output)
    array < CbArray3D<double>, 4 > matrix = { sumVaVxMatrix, sumVaVyMatrix, sumVaVzMatrix, sumVaPrMatrix };
    writeMatrixToImageFile(output, matrix);
 }
+void Averaging::writeMeanMatrixToImageFile(std::string output)
+{
+   array < CbArray3D<double>, 4 > matrix = { meanVxMatrix, meanVyMatrix, meanVzMatrix, meanPrMatrix };
+   writeMatrixToImageFile(output, matrix);
+}
 //////////////////////////////////////////////////////////////////////////
 void Averaging::volumeAveragingWithMPI(double l_real)
 {
@@ -451,7 +456,7 @@ void Averaging::volumeAveragingWithMPI(double l_real)
                         if (yy < 0)   yy = dimensions[1] + yy;
                         if (yy >= dimensions[1]) yy = yy - dimensions[1];
 
-                        if (zz < 0)   zz = 0;
+                        if (zz < 0)   zz = dimensions[2] + zz; 
                         if (zz >= dimensions[2]) zz = dimensions[2] - 1;
 
                         double mm = (G((double)x, l)*G((double)y, l)*G((double)z, l)) / lNorm;
@@ -517,9 +522,9 @@ void Averaging::volumeAveragingWithMPI(double l_real)
             for (int x1 = startX1; x1 < stopX1; x1++)
             {
                sendBuffer.push_back(vaVxMatrix(x1, x2, x3));
-               sendBuffer.push_back(vaVxMatrix(x1, x2, x3));
-               sendBuffer.push_back(vaVxMatrix(x1, x2, x3));
-               sendBuffer.push_back(vaVxMatrix(x1, x2, x3));
+               sendBuffer.push_back(vaVyMatrix(x1, x2, x3));
+               sendBuffer.push_back(vaVzMatrix(x1, x2, x3));
+               sendBuffer.push_back(vaPrMatrix(x1, x2, x3));
             }
       int count = (int)sendBuffer.size();
       MPI_Send(&count, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
@@ -581,11 +586,25 @@ void Averaging::readGeoMatrix(string dataNameG)
 }
 void Averaging::writeGeoMatrixToBinaryFiles(std::string fname)
 {
-   writeMatrixToBinaryFiles<int>(geoMatrix.getDataVector(), fname);
+   writeMatrixToBinaryFiles<int>(geoMatrix, fname);
 }
 void Averaging::readGeoMatrixFromBinaryFiles(std::string fname)
 {
-   readMatrixFromBinaryFiles<int>(fname, geoMatrix.getDataVector());
+   readMatrixFromBinaryFiles<int>(fname, geoMatrix);
+}
+void Averaging::writeMqMatrixToBinaryFiles(std::string fname, int timeStep)
+{
+   writeMatrixToBinaryFiles<double>(vxMatrix, fname + "Vx" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vyMatrix, fname + "Vy" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vzMatrix, fname + "Vz" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(prMatrix, fname + "Pr" + UbSystem::toString(timeStep) + ".bin");
+}
+void Averaging::readMqMatrixFromBinaryFiles(std::string fname, int timeStep)
+{
+   readMatrixFromBinaryFiles<double>(fname + "Vx" + UbSystem::toString(timeStep) + ".bin", vxMatrix);
+   readMatrixFromBinaryFiles<double>(fname + "Vy" + UbSystem::toString(timeStep) + ".bin", vyMatrix);
+   readMatrixFromBinaryFiles<double>(fname + "Vz" + UbSystem::toString(timeStep) + ".bin", vzMatrix);
+   readMatrixFromBinaryFiles<double>(fname + "Pr" + UbSystem::toString(timeStep) + ".bin", prMatrix);
 }
 void Averaging::initVolumeAveragingValues()
 {
@@ -708,10 +727,10 @@ void Averaging::sumOfVolumeAveragingValues()
 }
 void Averaging::writeVolumeAveragingValuesToBinaryFiles(std::string ffname, int timeStep)
 {
-   writeMatrixToBinaryFiles<double>(vaVxMatrix.getDataVector(), ffname + "Vx" + UbSystem::toString(timeStep)+".bin");
-   writeMatrixToBinaryFiles<double>(vaVyMatrix.getDataVector(), ffname + "Vy" + UbSystem::toString(timeStep)+".bin");
-   writeMatrixToBinaryFiles<double>(vaVzMatrix.getDataVector(), ffname + "Vz" + UbSystem::toString(timeStep)+".bin");
-   writeMatrixToBinaryFiles<double>(vaPrMatrix.getDataVector(), ffname + "Pr" + UbSystem::toString(timeStep)+".bin");
+   writeMatrixToBinaryFiles<double>(vaVxMatrix, ffname + "Vx" + UbSystem::toString(timeStep)+".bin");
+   writeMatrixToBinaryFiles<double>(vaVyMatrix, ffname + "Vy" + UbSystem::toString(timeStep)+".bin");
+   writeMatrixToBinaryFiles<double>(vaVzMatrix, ffname + "Vz" + UbSystem::toString(timeStep)+".bin");
+   writeMatrixToBinaryFiles<double>(vaPrMatrix, ffname + "Pr" + UbSystem::toString(timeStep)+".bin");
 }
 void Averaging::meanOfVolumeAveragingValues(int numberOfTimeSteps)
 {
@@ -737,10 +756,10 @@ void Averaging::meanOfVolumeAveragingValues(int numberOfTimeSteps)
 }
 void Averaging::writeMeanVolumeAveragingValuesToBinaryFiles(std::string ffname)
 {
-   writeMatrixToBinaryFiles<double>(vaVxMatrix.getDataVector(), ffname + "Vx" + ".bin");
-   writeMatrixToBinaryFiles<double>(vaVyMatrix.getDataVector(), ffname + "Vy" + ".bin");
-   writeMatrixToBinaryFiles<double>(vaVzMatrix.getDataVector(), ffname + "Vz" + ".bin");
-   writeMatrixToBinaryFiles<double>(vaPrMatrix.getDataVector(), ffname + "Pr" + ".bin");
+   writeMatrixToBinaryFiles<double>(vaVxMatrix, ffname + "Vx" + ".bin");
+   writeMatrixToBinaryFiles<double>(vaVyMatrix, ffname + "Vy" + ".bin");
+   writeMatrixToBinaryFiles<double>(vaVzMatrix, ffname + "Vz" + ".bin");
+   writeMatrixToBinaryFiles<double>(vaPrMatrix, ffname + "Pr" + ".bin");
 }
 void Averaging::fluctuationsOfVolumeAveragingValue()
 {
@@ -795,22 +814,22 @@ void Averaging::fluctuationsOfVolumeAveragingValue()
 }
 void Averaging::writeFluctuationsToBinaryFiles(std::string fname, int timeStep)
 {
-   writeMatrixToBinaryFiles<double>(FlucVxMatrix.getDataVector(), fname + "Vx" + UbSystem::toString(timeStep) + ".bin");
-   writeMatrixToBinaryFiles<double>(FlucVyMatrix.getDataVector(), fname + "Vy" + UbSystem::toString(timeStep) + ".bin");
-   writeMatrixToBinaryFiles<double>(FlucVzMatrix.getDataVector(), fname + "Vz" + UbSystem::toString(timeStep) + ".bin");
-   writeMatrixToBinaryFiles<double>(FlucPrMatrix.getDataVector(), fname + "Pr" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(FlucVxMatrix, fname + "Vx" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(FlucVyMatrix, fname + "Vy" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(FlucVzMatrix, fname + "Vz" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(FlucPrMatrix, fname + "Pr" + UbSystem::toString(timeStep) + ".bin");
 }
 void Averaging::writeStressesToBinaryFiles(std::string fname, int timeStep)
 {
-   writeMatrixToBinaryFiles<double>(StressXX.getDataVector(), fname + UbSystem::toString(timeStep) + "XX" + ".bin");
-   writeMatrixToBinaryFiles<double>(StressXY.getDataVector(), fname + UbSystem::toString(timeStep) + "XY" + ".bin");
-   writeMatrixToBinaryFiles<double>(StressXZ.getDataVector(), fname + UbSystem::toString(timeStep) + "XZ" + ".bin");
-   writeMatrixToBinaryFiles<double>(StressYX.getDataVector(), fname + UbSystem::toString(timeStep) + "YX" + ".bin");
-   writeMatrixToBinaryFiles<double>(StressYY.getDataVector(), fname + UbSystem::toString(timeStep) + "YY" + ".bin");
-   writeMatrixToBinaryFiles<double>(StressYZ.getDataVector(), fname + UbSystem::toString(timeStep) + "YZ" + ".bin");
-   writeMatrixToBinaryFiles<double>(StressZX.getDataVector(), fname + UbSystem::toString(timeStep) + "ZX" + ".bin");
-   writeMatrixToBinaryFiles<double>(StressZY.getDataVector(), fname + UbSystem::toString(timeStep) + "ZY" + ".bin");
-   writeMatrixToBinaryFiles<double>(StressZZ.getDataVector(), fname + UbSystem::toString(timeStep) + "ZZ" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressXX, fname + UbSystem::toString(timeStep) + "XX" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressXY, fname + UbSystem::toString(timeStep) + "XY" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressXZ, fname + UbSystem::toString(timeStep) + "XZ" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressYX, fname + UbSystem::toString(timeStep) + "YX" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressYY, fname + UbSystem::toString(timeStep) + "YY" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressYZ, fname + UbSystem::toString(timeStep) + "YZ" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressZX, fname + UbSystem::toString(timeStep) + "ZX" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressZY, fname + UbSystem::toString(timeStep) + "ZY" + ".bin");
+   writeMatrixToBinaryFiles<double>(StressZZ, fname + UbSystem::toString(timeStep) + "ZZ" + ".bin");
 }
 void Averaging::sumOfFluctuations()
 {
@@ -1023,10 +1042,10 @@ void Averaging::planarAveragingMQ(std::array<int, 3> dimensions)
    }
 void Averaging::readVolumeAveragingValuesFromBinaryFiles(std::string fname, int timeStep)
 {
-   readMatrixFromBinaryFiles<double>(fname + "Vx" + UbSystem::toString(timeStep) + ".bin", vaVxMatrix.getDataVector());
-   readMatrixFromBinaryFiles<double>(fname + "Vy" + UbSystem::toString(timeStep) + ".bin", vaVyMatrix.getDataVector());
-   readMatrixFromBinaryFiles<double>(fname + "Vz" + UbSystem::toString(timeStep) + ".bin", vaVzMatrix.getDataVector());
-   readMatrixFromBinaryFiles<double>(fname + "Pr" + UbSystem::toString(timeStep) + ".bin", vaPrMatrix.getDataVector());
+   readMatrixFromBinaryFiles<double>(fname + "Vx" + UbSystem::toString(timeStep) + ".bin", vaVxMatrix);
+   readMatrixFromBinaryFiles<double>(fname + "Vy" + UbSystem::toString(timeStep) + ".bin", vaVyMatrix);
+   readMatrixFromBinaryFiles<double>(fname + "Vz" + UbSystem::toString(timeStep) + ".bin", vaVzMatrix);
+   readMatrixFromBinaryFiles<double>(fname + "Pr" + UbSystem::toString(timeStep) + ".bin", vaPrMatrix);
 }
 //////////////////////////////////////////////////////////////////////////
 double Averaging::G(double x, double l)
diff --git a/source/Applications/bChannelVA/Averaging.h b/source/Applications/bChannelVA/Averaging.h
index 46720264056b48742e519747dd112349696657bf..c8b8d692263b9d4e6931510d44fdbc2e6c3eb368 100644
--- a/source/Applications/bChannelVA/Averaging.h
+++ b/source/Applications/bChannelVA/Averaging.h
@@ -16,10 +16,13 @@ public:
    void writeMqMatrixToImageFile(std::string output);
    void writeVaMatrixToImageFile(std::string output);
    void writeVaSumMatrixToImageFile(std::string output);
+   void writeMeanMatrixToImageFile(std::string output);
    void volumeAveragingWithMPI(double l_real);
    void readGeoMatrix(std::string dataNameG);
    void writeGeoMatrixToBinaryFiles(std::string fname);
    void readGeoMatrixFromBinaryFiles(std::string fname);
+   void writeMqMatrixToBinaryFiles(std::string fname, int timeStep);
+   void readMqMatrixFromBinaryFiles(std::string fname, int timeStep);
    void initVolumeAveragingValues();
    void initMeanVolumeAveragingValues();
    void initFluctuationsofVolumeAveragingValues();
@@ -55,9 +58,9 @@ protected:
    double G(double x, double l);
    
    template <class T>
-   void writeMatrixToBinaryFiles(std::vector<T> &matrix, std::string fname);
+   void writeMatrixToBinaryFiles(CbArray3D<T>& matrix, std::string fname);
    template <class T>
-   void readMatrixFromBinaryFiles(std::string fname, std::vector<T> &matrix);
+   void readMatrixFromBinaryFiles(std::string fname, CbArray3D<T>& matrix);
 private:
    std::array<int, 3> dimensions;
    std::array<int, 6> geo_extent;
@@ -151,7 +154,7 @@ private:
 };
 
 //////////////////////////////////////////////////////////////////////////
-template<class T> void Averaging::writeMatrixToBinaryFiles(std::vector<T> &matrix, std::string fname)
+template<class T> void Averaging::writeMatrixToBinaryFiles(CbArray3D<T>& matrix, std::string fname)
  {
    vtkSmartPointer<vtkTimerLog> timer_write = vtkSmartPointer<vtkTimerLog>::New();
 
@@ -169,7 +172,9 @@ template<class T> void Averaging::writeMatrixToBinaryFiles(std::vector<T> &matri
       if (!ostr) throw UbException(UB_EXARGS, "couldn't open file " + fname);
    }
 
-   ostr.write((char*)&matrix[0], sizeof(T)*matrix.size());
+   std::vector<T>& vec = matrix.getDataVector();
+
+   ostr.write((char*)& vec[0], sizeof(T)*vec.size());
    ostr.close();
 
    UBLOG(logINFO,"write matrix: end");
@@ -177,7 +182,7 @@ template<class T> void Averaging::writeMatrixToBinaryFiles(std::vector<T> &matri
    UBLOG(logINFO,"write matrix time: " + UbSystem::toString(timer_write->GetElapsedTime()) + " s");
 }
 //////////////////////////////////////////////////////////////////////////
-template<class T> void Averaging::readMatrixFromBinaryFiles(std::string fname, std::vector<T> &matrix)
+template<class T> void Averaging::readMatrixFromBinaryFiles(std::string fname, CbArray3D<T>& matrix)
 {
    vtkSmartPointer<vtkTimerLog> timer_write = vtkSmartPointer<vtkTimerLog>::New();
 
@@ -195,11 +200,14 @@ template<class T> void Averaging::readMatrixFromBinaryFiles(std::string fname, s
    rewind(file);
 
    // allocate memory to contain the whole file:
-   matrix.resize(lSize);
-   if (matrix.size() == 0) { fputs("Memory error", stderr); exit(2); }
+   //matrix.resize(lSize);
+   matrix.resize(dimensions[0], dimensions[1], dimensions[2]);
+   std::vector<T>& vec = matrix.getDataVector();
+
+   if (vec.size() == 0) { fputs("Memory error", stderr); exit(2); }
 
    // copy the file into the buffer:
-   size_t result = fread(&matrix[0], sizeof(T), lSize, file);
+   size_t result = fread(&vec[0], sizeof(T), lSize, file);
    if (result != lSize) { fputs("Reading error", stderr); exit(3); }
 
    fclose(file);
diff --git a/source/Applications/bChannelVA/bChannelVA.cpp b/source/Applications/bChannelVA/bChannelVA.cpp
index 81d484cd624a501978347b4def35f276392c14bd..60be26dec239d621b30d0bb24d2f3d2b7eb17e4f 100644
--- a/source/Applications/bChannelVA/bChannelVA.cpp
+++ b/source/Applications/bChannelVA/bChannelVA.cpp
@@ -35,7 +35,7 @@ int main(int argc, char* argv[])
         int startTimeStep = 60000;
         int timeStep = 1000;
         int numberOfTimeSteps = 61000;
-        int numberOfSamples = numberOfTimeSteps / startTimeStep;
+        int numberOfSamples = (numberOfTimeSteps - startTimeStep)/timeStep + 1;
         int numberOfGridPoints = dimensions[0] * dimensions[1] * dimensions[2];
 
         Averaging av;
@@ -53,19 +53,38 @@ int main(int argc, char* argv[])
         for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
         {
             av.createMQMatrix(pathIn + "/mq/mq" + UbSystem::toString(t) + ".pvtu");
-            av.volumeAveragingWithMPI(real_l);
-            if (myid == 0) av.sumOfVolumeAveragingValues();
-            if (myid == 0) av.writeVolumeAveragingValuesToBinaryFiles(pathOut + "/va/vav/vav", t);
+            if (myid == 0) av.writeMqMatrixToBinaryFiles(pathOut + "/va/mq/mq" + UbSystem::toString(t) + "/mq", t);
+            //av.volumeAveragingWithMPI(real_l);
+            //if (myid == 0)
+            //{
+            //   //av.sumOfVolumeAveragingValues();
+            //   av.writeVolumeAveragingValuesToBinaryFiles(pathOut + "/va/vav/vav", t);
+            //}
         }
 
         //av.writeMqMatrixToImageFile(pathOut + "/va/vav/mq");
-        //av.writeVaMatrixToImageFile(pathOut + "/va/vav/vav");
-        if (myid == 0) av.writeVaSumMatrixToImageFile(pathOut + "/va/vav/vavSum");
-        return 0;
+        //if (myid == 0) av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", 60000);
+        //if (myid == 0) av.writeVaMatrixToImageFile(pathOut + "/va/vav/vav");
+
+        //if (myid == 0)
+        //{
+        //   av.initMeanVolumeAveragingValues();
+        //   av.meanOfVolumeAveragingValues(numberOfSamples);
+        //   av.writeMeanVolumeAveragingValuesToBinaryFiles(pathOut + "/va/mean/mean");
+        //   av.writeMeanMatrixToImageFile(pathOut + "/va/vav/vavMean");
+        //}
+
+        //for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
+        //{
+        //   av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", t);
+        //   av.sumOfVolumeAveragingValues();
+        //}
 
         //av.initMeanVolumeAveragingValues();
         //av.meanOfVolumeAveragingValues(numberOfSamples);
         //av.writeMeanVolumeAveragingValuesToBinaryFiles(pathOut + "/va/mean/mean");
+        //av.writeMeanMatrixToImageFile(pathOut + "/va/vav/vavMean");
+
         //av.initFluctuationsofVolumeAveragingValues();
         //av.initSumOfFluctuations();
         //av.initMeanOfFluctuations();
@@ -74,6 +93,7 @@ int main(int argc, char* argv[])
         //av.initMeanOfStresses();
         //av.initPlanarAveragingMQ();
 
+
         //for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
         //{
         //    av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", t);