diff --git a/src/.gitignore b/src/.gitignore
index 44599544..ce60e8c7 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -5,3 +5,5 @@
 log_*
 log.*
 *~
+
+lnInclude
diff --git a/src/recurrence/Make/files b/src/recurrence/Make/files
index cd69ddc6..cacbc10d 100644
--- a/src/recurrence/Make/files
+++ b/src/recurrence/Make/files
@@ -2,6 +2,7 @@ recBase/recBase.C
 recModel/recModel/recModel.C
 recModel/recModel/newRecModel.C
 recModel/standardRecModel/standardRecModel.C
+recModel/gerhardsRecModel/gerhardsRecModel.C
 recNorm/recNorm/recNorm.C
 recNorm/recNorm/newRecNorm.C
 recNorm/sqrDiffNorm/sqrDiffNorm.C
diff --git a/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModel.C b/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModel.C
new file mode 100644
index 00000000..cb6b8889
--- /dev/null
+++ b/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModel.C
@@ -0,0 +1,1026 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "Random.H"
+#include "gerhardsRecModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(gerhardsRecModel, 0);
+
+addToRunTimeSelectionTable
+(
+    recModel,
+    gerhardsRecModel,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+gerhardsRecModel::gerhardsRecModel
+(
+    const dictionary& dict,
+    recBase& base
+)
+:
+    recModel(dict,base),
+    propsDict_(dict.subDict(typeName + "Props")),
+    dataBaseName_(propsDict_.lookupOrDefault<word>("dataBase", word("dataBase"))),
+    recTime(fileName(dataBaseName_), "", "../system", "../constant", false),
+    timeDirs(recTime.times()),
+    skipZero_(propsDict_.lookupOrDefault<Switch>("skipZero", Switch(false))),
+    checkSkipZero_(checkSkipZero()),
+    numRecFields_(skipZero_ ? label(timeDirs.size())-1 : label(timeDirs.size())),
+    verboseVerbose_(propsDict_.lookupOrDefault<Switch>("verboseVerbose", Switch(false))),
+    verboseVacate_(propsDict_.lookupOrDefault<Switch>("verboseVacate", Switch(false))),
+    recurrenceMatrix_(numRecFields_, scalar(-1.0)),
+    timeIndexList_(numRecFields_-1),
+    revTimeIndexList_(numRecFields_-1),
+    timeValueList_(numRecFields_-1),
+    contTimeIndex(0),
+    lowerSeqLim_(max(1, label(numRecFields_/20))),
+    upperSeqLim_(label(numRecFields_/5)),
+    numDataBaseFields_(propsDict_.lookupOrDefault<label>("numDataBaseFields",numRecFields_)),
+    storageIndex_(numRecFields_),
+    storageUsageList_(numRecFields_),
+    volScalarFieldList_(volScalarFieldNames_.size()),
+    volVectorFieldList_(volVectorFieldNames_.size()),
+    surfaceScalarFieldList_(surfaceScalarFieldNames_.size()),
+    nrOfReadsFromDisk_(0)
+{
+    /*
+        Sanity checks
+        
+        "Sanity is just another weird form of Madness"
+                                                --- Darkwell
+    */
+    if (numDataBaseFields_ < 2)
+    {
+        FatalError << "Number of fields in dataBase specified in "
+            << propsDict_.name() << " is smaller than 2!"
+            << abort(FatalError);
+    }
+    
+    if (numDataBaseFields_ > numRecFields_)
+    {   
+        FatalError << "Number of fields in dataBase specified in "
+            << propsDict_.name() << " is larger than number of snapshots!"
+            << abort(FatalError);
+    }
+    
+    if (dataBaseName_ == "")
+    {
+        FatalError << "Empty dataBase path provided"
+            << abort(FatalError);
+    }
+    
+    if (upperSeqLim_ == 0)
+    {
+        FatalError << "Bad sequence limits! The dataBase is most probably too small."
+            << abort(FatalError);
+    }
+    
+    // initialize data structures
+    forAll(storageUsageList_, i)
+    {
+        storageUsageList_[i] = storageUsageList_.size()+1;
+    }
+    
+    if (verbose_)
+    {
+        // be informative on properties of the "recTime" Time-object
+        Info << "recTime.rootPath() " << recTime.rootPath() << endl;
+        Info << "recTime.caseName() " << recTime.caseName() << endl;
+        Info << "recTime.path() " << recTime.path() << endl;
+        Info << "recTime.timePath() " << recTime.timePath() << endl;
+        Info << "recTime.timeName() " << recTime.timeName() << endl;
+        Info << "timeDirs " << timeDirs << endl;
+        Info << "ignoring 0 directory: " << skipZero_ << endl;
+        Info << "number of snapshots: " << numRecFields_ << endl;
+        Info << "number of dataBase slots: " << numDataBaseFields_ << endl;
+    }
+    readTimeSeries();
+
+    recTimeStep_ = checkTimeStep();
+    totRecSteps_ = 1 + static_cast<label>( (endTime_-startTime_) / recTimeStep_ );
+
+    for(int i=0; i<volScalarFieldNames_.size(); i++)
+    {
+        volScalarFieldList_[i].setSize(numDataBaseFields_);
+    }
+
+    for(int i=0; i<volVectorFieldNames_.size(); i++)
+    {
+        volVectorFieldList_[i].setSize(numDataBaseFields_);
+    }
+
+    for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
+    {
+        surfaceScalarFieldList_[i].setSize(numDataBaseFields_);
+    }
+
+    //   setRecFields();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+gerhardsRecModel::~gerhardsRecModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+scalar gerhardsRecModel::checkTimeStep()
+{
+   // check time step of provided data
+    scalar dtCur(0.0);
+    scalar dtOld(0.0);
+    
+    if (verbose_)
+    {
+    	Info << "timeValueList : " << timeValueList_ << endl;
+    }
+    
+    forAll(timeValueList_, i)
+    {
+    	// skip zero
+    	if (skipZero_ and timeDirs[i].value() == 0)
+    	{
+    	    if (verbose_)
+    	    {
+    	        Info << " ... skipping 0 in checkTimeStep()" << endl;
+    	    }
+    	    continue;
+    	}
+    	
+    	// compute time step
+    	if (timeDirs[i].value() == timeDirs.last().value())
+    	{
+    		if (verbose_)
+    		{
+    			Info << ".. leaving loop at " << timeDirs[i] << endl;
+    		}
+    		// leave loop
+    		break;
+    	}
+    	
+    	if (verbose_)
+    	{
+    		Info << "timeDirs.fcIndex(i)].value(),  timeDirs[i].value() : " 
+    			<< timeDirs[timeDirs.fcIndex(i)].value() << "   " << timeDirs[i].value()
+    			<< endl;
+    	}
+    	
+    	// the documentation is in the code ;-)
+    	//	fcIndex() - return forward circular index, i.e. the next index
+        dtCur = timeDirs[timeDirs.fcIndex(i)].value() - timeDirs[i].value();
+        
+        if (dtOld < SMALL)
+        {
+        	dtOld = dtCur;
+        }
+        
+        if (abs(dtOld - dtCur) > SMALL)
+        {
+        	Info << "dtCur, dtOld = " << dtCur << "   " << dtOld << endl;
+        	FatalError << "    in setting up data" << nl
+				<< "    non-constant time-step of provided simulation data" 
+				<< abort(FatalError);
+        }
+    }
+        
+    // set deltaT
+    recTime.setDeltaT(dtCur, false);
+	
+	if (verbose_)
+    {
+		Info << "Setting deltaRecT to " << dtCur << endl;
+		Info << "Actual recTime.deltaT = " << recTime.deltaTValue() << endl;
+		Info << "Actual runTime.deltaT = " << timeStep_ << endl;
+    }
+    return dtCur;
+}
+
+void gerhardsRecModel::readFieldSeries()
+{
+    Info << "Checking all " << numRecFields() << " fields\n" << endl;
+
+    label size = timeDirs.size();
+    label counter = 0;
+    label percentage = 0;
+
+    for (instantList::iterator it=timeDirs.begin(); it != timeDirs.end(); ++it)
+    {
+        if(counter >= 0.1 * percentage * size)
+        {
+            Info << "\t" << 10 * percentage << " \% done" << endl;
+            percentage++;
+        }
+        counter++;
+
+        // set time
+        recTime.setTime(*it, it->value());
+
+        // skip zero
+        if (skipZero_ and recTime.timeName() == "0")
+        {
+            if (verbose_)
+    	    {
+    	        Info << " ... skipping 0 in readFieldSeries()" << endl;
+    	    }
+    	    
+    	    continue;
+        }
+        
+        // skip constant
+        if (recTime.timeName() == "constant")
+        {
+            continue;
+        }
+
+        if (verbose_)
+        {
+            Info << "Checking fields at t = " << recTime.timeName() << endl;
+        }
+        
+        for (int i=0; i<volScalarFieldNames_.size(); i++)
+        {
+            IOobject header
+            (
+                volScalarFieldNames_[i],
+                recTime.timePath(),
+                base_.mesh(),
+                IOobject::MUST_READ
+            );
+            
+            // Check if volScalarFieldNames_[i] is a valid field
+            /*
+                COMPAT: replace the method typeHeaderOk<volScalarField>(true) with headerOk() 
+                for OpenFOAM versions prior to OpenFOAM-5.0
+                Do this the other way around for OpenFOAM-5.0 and potentially later versions
+            */
+            if (not header.headerOk())
+            {
+                FatalError
+                << "Field " << volScalarFieldNames_[i] << " not found"
+                << exit(FatalError);
+            }
+        }
+        
+        for (int i=0; i<volVectorFieldNames_.size(); i++)
+        {
+            IOobject header
+            (
+                volVectorFieldNames_[i],
+                recTime.timePath(),
+                base_.mesh(),
+                IOobject::MUST_READ
+            );
+            
+            // Check if volVectorFieldNames_[i] is a valid field
+            /*
+                COMPAT: replace the method typeHeaderOk<volScalarField>(true) with headerOk() 
+                for OpenFOAM versions prior to OpenFOAM-5.0
+                Do this the other way around for OpenFOAM-5.0 and potentially later versions
+            */
+            if (not header.headerOk())
+            {
+                FatalError
+                << "Field " << volVectorFieldNames_[i] << " not found"
+                << exit(FatalError);
+            }
+        }
+        
+        for (int i=0; i<surfaceScalarFieldNames_.size(); i++)
+        {
+            IOobject header
+            (
+                surfaceScalarFieldNames_[i],
+                recTime.timePath(),
+                base_.mesh(),
+                IOobject::MUST_READ
+            );
+            
+            // Check if surfaceScalarFieldNames_[i] is a valid field
+            /*
+                COMPAT: replace the method typeHeaderOk<volScalarField>(true) with headerOk() 
+                for OpenFOAM versions prior to OpenFOAM-5.0
+                Do this the other way around for OpenFOAM-5.0 and potentially later versions
+            */
+            if (not header.headerOk())
+            {
+                FatalError
+                << "Field " << surfaceScalarFieldNames_[i] << " not found"
+                << exit(FatalError);
+            }
+        }
+        
+        // reset storageIndex_ to -1
+        storageIndex_[timeIndexList_(recTime.timeName())] = -1;
+    }
+    
+    Info << "Checking fields done" << nl << endl;
+    
+    
+    Info << "Reading fields for " << numDataBaseFields() << " dataBase slots\n" << endl;
+    
+    size = numDataBaseFields();
+    counter = 0;
+    percentage = 0;
+    
+    
+    label fieldCounter(0);
+    for (instantList::iterator it=timeDirs.begin(); it != timeDirs.end(); ++it)
+    {
+        if(counter >= 0.1 * percentage * size)
+        {
+            Info << "\t" << 10 * percentage << " \% done" << endl;
+            percentage++;
+        }
+        counter++;
+
+        // set time
+        recTime.setTime(*it, it->value());
+
+        // skip constant
+        if (recTime.timeName() == "constant")
+        {
+            continue;
+        }
+        
+        // skip zero
+        if (skipZero_ and recTime.timeName() == "0")
+        {
+            if (verbose_)
+    	    {
+    	        Info << " ... skipping 0 in readTimeSeries()" << endl;
+    	    }
+            
+            continue;
+        }
+
+        if (verbose_)
+        {
+            Info << "Reading fields at t = " << recTime.timeName() << endl;
+        }
+        
+        nrOfReadsFromDisk_++;
+        
+        for (int i=0; i<volScalarFieldNames_.size(); i++)
+        {
+            volScalarFieldList_[i].set
+            (
+                fieldCounter,
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        volScalarFieldNames_[i],
+                        recTime.timePath(),
+                        base_.mesh(),
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    base_.mesh()
+                )
+            );
+        }
+        
+        for (int i=0; i<volVectorFieldNames_.size(); i++)
+        {
+            volVectorFieldList_[i].set
+            (
+                fieldCounter,
+                new volVectorField
+                (
+                    IOobject
+                    (
+                        volVectorFieldNames_[i],
+                        recTime.timePath(),
+                        base_.mesh(),
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    base_.mesh()
+                )
+            );
+        }
+        
+        for (int i=0; i<surfaceScalarFieldNames_.size(); i++)
+        {
+            surfaceScalarFieldList_[i].set
+            (
+                fieldCounter,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        surfaceScalarFieldNames_[i],
+                        recTime.timePath(),
+                        base_.mesh(),
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    base_.mesh()
+                )
+            );
+        }
+        
+        // set storageIndex_
+        //  storageIndex_[index] translates the index of the snapshot 
+        //  to an appropriate index of the dataBase
+        storageIndex_[timeIndexList_(recTime.timeName())] = fieldCounter;
+        
+        
+        fieldCounter++;
+        
+        if (fieldCounter >= numDataBaseFields_)
+        {
+            break;
+        }
+    }
+    
+    Info << "Reading fields done" << endl;
+}
+
+
+void gerhardsRecModel::readTimeSeries()
+{
+    bool firsttime = true;
+    // fill the data structure for the time indices
+    for (instantList::iterator it=timeDirs.begin(); it != timeDirs.end(); ++it)
+    {
+    	// set run-time
+    	recTime.setTime(*it, it->value());
+    	
+    	
+    	// skip constant
+    	if (recTime.timeName() == "constant")
+    	{
+        	continue;
+        }
+        
+        // skip zero
+        if (skipZero_ and recTime.timeName() == "0")
+        {
+            if (verbose_)
+    	    {
+    	        Info << " ... skipping 0 in readTimeSeries()" << endl;
+    	    }
+            
+            continue;
+        }
+        
+        if (firsttime)
+        {
+            firsttime = false;
+            recStartTime_ = recTime.value();
+        }
+        recEndTime_ = recTime.value();
+        
+        // insert the time name into the hash-table with a continuous second index
+        timeIndexList_.insert(recTime.timeName(), contTimeIndex);
+        revTimeIndexList_.insert(contTimeIndex, recTime.timeName());
+        
+        
+        if (verbose_)
+    	{
+    		Info << "current time " << recTime.timeName() << endl;
+    		Info << "insert " << recTime.timeName() << " , " << contTimeIndex << endl;
+    	}
+        
+        
+        // insert the time value
+        timeValueList_.insert(contTimeIndex, recTime.timeOutputValue());
+        
+        // increment continuousIndex
+        contTimeIndex++;
+        
+        if (verbose_)
+    	{
+        	Info << "contTimeIndex " << contTimeIndex << endl;
+        }
+    }
+
+    if (verbose_)
+    {
+        Info << endl;
+        Info << "Found " << label(timeDirs.size()) << " time folders" << endl;
+        Info << "Found " << label(timeIndexList_.size()) << " time steps" << endl;
+        Info << "database start time = " << recStartTime_ << endl;
+        Info << "database end time = " << recEndTime_ << endl;
+    }
+}
+
+
+
+void gerhardsRecModel::fetchStateForDataBase(label index)
+{
+    if (verboseVerbose_)
+    {
+        Info << " ... calling fetchStateForDataBase(" << index << ")" << endl;
+        Info << "     storageUsageList_ : " << storageUsageList_ << endl;
+        Info << "     storageIndex_ : " << storageIndex_ << endl;
+    }
+    
+    // decide which state to remove from the data base
+    //  find maximum element in storageUsageList_, as it is the index of the least used state
+    label oldIndex(0);
+    label maxVal(0);
+    
+    // find maximum value in storageUsageList_
+    forAll(storageUsageList_, i)
+    {
+        if (storageUsageList_[i] > maxVal and storageIndex_[i] > -1)
+        {
+            oldIndex = i;
+            maxVal = storageUsageList_[i];
+        }
+    }
+    if (verboseVerbose_ or verboseVacate_)
+    {
+        Info << " --> vacate dataBase  : " << oldIndex << endl;
+    }
+    
+    // read state for index from disk and put it into the data base
+    readNewSnapshot(index, oldIndex);
+    
+    
+    // reset storageIndex
+    //  the slot at oldIndex in the dataBase is replaced by the 
+    //  snapshot for index
+    storageIndex_[index] = storageIndex_[oldIndex];
+    storageIndex_[oldIndex] = -1;
+    
+    // update storageIndex
+    updateStorageUsageList(index);
+    
+    if (verboseVerbose_)
+    {
+        Info << " ... finished fetchStateForDataBase(" << index << ")" << endl;
+        Info << "     storageUsageList_ : " << storageUsageList_ << endl;
+        Info << "     storageIndex_ : " << storageIndex_ << endl;
+    }
+}
+
+
+void gerhardsRecModel::updateStorageUsageList(label index)
+{
+    if (verboseVerbose_)
+    {
+        Info << " ... calling updateStorageUsageList(" << index << ")" << endl;
+        Info << "     storageUsageList_ : " << storageUsageList_ << endl;
+        Info << "   storageIndex_ : " << storageIndex_ << endl;
+    }
+    
+    label oldValue(storageUsageList_[index]);
+    
+    /*
+        increment (make less recent) all elements more recent (smaller) 
+        than the element at slot index
+    */
+    forAll(storageUsageList_, i)
+    {
+        if (storageUsageList_[i] < oldValue)
+        {
+            storageUsageList_[i]++;
+        }
+    }
+    
+    // storageUsageList_[index] is the most recently used element in dataBase
+    storageUsageList_[index] = 1;
+    
+    if (verboseVerbose_)
+    {
+        Info << " ... finished updateStorageUsageList(" << index << ")" << endl;
+        Info << "     storageUsageList_ : " << storageUsageList_ << endl;
+        Info << "     storageIndex_ : " << storageIndex_ << nl << endl;
+    }
+}
+
+
+void gerhardsRecModel::readNewSnapshot(label index, label oldIndex)
+{
+    if (verboseVerbose_)
+    {
+        Info << "      ... calling readNewSnapshot(" << index << ", " << oldIndex << ")" << endl;
+        Info << "          index = " << index << "; oldIndex = " << oldIndex << endl;
+        Info << "          storageUsageList_ : " << storageUsageList_ << endl;
+        Info << "          storageIndex_ : " << storageIndex_ << nl << endl;
+    }
+    
+    /*
+        Replace all fields in the slot "oldIndex" of the dataBase with fields 
+        for the snapshot at "index"
+    */
+    
+    // vacate slot
+    for (int i=0; i<volScalarFieldNames_.size(); i++)
+    {
+        volScalarFieldList_[i][storageIndex_[oldIndex]].clear();
+    }
+    for (int i=0; i<volVectorFieldNames_.size(); i++)
+    {
+        volVectorFieldList_[i][storageIndex_[oldIndex]].clear();
+    }
+    for (int i=0; i<surfaceScalarFieldNames_.size(); i++)
+    {
+        surfaceScalarFieldList_[i][storageIndex_[oldIndex]].clear();
+    }
+    
+    
+    // build path to read from
+    word readPath(dataBaseName_+"/"+timeDirs[index].name());
+    
+    if (dataBaseName_ == ".")
+    {
+        readPath = timeDirs[index].name();
+    }
+    
+    nrOfReadsFromDisk_++;
+    
+    // repopulate slot
+    for (int i=0; i<volScalarFieldNames_.size(); i++)
+    {
+        if (verboseVerbose_)
+        {
+            Info << " .... reading " << volScalarFieldNames_[i] << " at time " 
+                << timeDirs[index] << endl;
+        }
+        
+        volScalarFieldList_[i].set
+        (
+            storageIndex_[oldIndex],
+            new volScalarField
+            (
+                IOobject
+                (
+                    volScalarFieldNames_[i],
+                    readPath,
+                    base_.mesh(),
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                base_.mesh()
+            )
+        );
+    }
+        
+    for (int i=0; i<volVectorFieldNames_.size(); i++)
+    {
+        volVectorFieldList_[i].set
+        (
+            storageIndex_[oldIndex],
+            new volVectorField
+            (
+                IOobject
+                (
+                    volVectorFieldNames_[i],
+                    readPath,
+                    base_.mesh(),
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                base_.mesh()
+            )
+        );
+    }
+    
+    for (int i=0; i<surfaceScalarFieldNames_.size(); i++)
+    {
+        surfaceScalarFieldList_[i].set
+        (
+            storageIndex_[oldIndex],
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    surfaceScalarFieldNames_[i],
+                    readPath,
+                    base_.mesh(),
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                base_.mesh()
+            )
+        );
+    }
+
+}
+
+
+
+
+void gerhardsRecModel::exportVolScalarField(word fieldname, volScalarField& field) 
+{
+    field = exportVolScalarField(fieldname, virtualTimeIndex);
+}
+
+
+void gerhardsRecModel::exportVolVectorField(word fieldname, volVectorField& field) 
+{
+    field = exportVolVectorField(fieldname, virtualTimeIndex);
+}
+
+
+void gerhardsRecModel::exportSurfaceScalarField(word fieldname, surfaceScalarField& field) 
+{
+    field = exportSurfaceScalarField(fieldname, virtualTimeIndex);
+}
+
+
+
+
+const volScalarField& gerhardsRecModel::exportVolScalarField(word fieldname, label index) 
+{
+    const label fieldI = getVolScalarFieldIndex(fieldname, index);
+    
+    if (verboseVerbose_)
+    {
+        Info << nl << " ... calling exportVolScalarField(" << fieldname << ", " << index << ")" << endl;
+    }
+    
+    if (numDataBaseFields_ == numRecFields_)
+    {
+        return volScalarFieldList_[fieldI][index];
+    }
+    
+    // check whether field with index is in the dataBase
+    //  translate index to storageIndex
+    if (storageIndex_[index] < 0)
+    {
+        fetchStateForDataBase(index);
+    }
+    else
+    {
+        updateStorageUsageList(index);
+    }
+
+    return volScalarFieldList_[fieldI][storageIndex_[index]];
+}
+
+
+
+const volVectorField& gerhardsRecModel::exportVolVectorField(word fieldname, label index) 
+{
+    const label fieldI = getVolVectorFieldIndex(fieldname, index);
+    
+    if (verboseVerbose_)
+    {
+        Info << nl << " ... calling exportVolVectorField(" << fieldname << ", " << index << ")" << endl;
+    }
+    
+    if (numDataBaseFields_ == numRecFields_)
+    {
+        return volVectorFieldList_[fieldI][index];
+    }
+    
+    
+    // check whether field with index is in the dataBase
+    //  translate intex to storageIndex
+    if (storageIndex_[index] < 0)
+    {
+        fetchStateForDataBase(index);
+    }
+    else
+    {
+        updateStorageUsageList(index);
+    }
+
+    return volVectorFieldList_[fieldI][storageIndex_[index]];
+}
+
+const surfaceScalarField& gerhardsRecModel::exportSurfaceScalarField(word fieldname, label index) 
+{
+    const label fieldI = getSurfaceScalarFieldIndex(fieldname, index);
+    
+    if (numDataBaseFields_ == numRecFields_)
+    {
+        return surfaceScalarFieldList_[fieldI][index];
+    }
+    
+    // check whether field with index is in the dataBase
+    //  translate intex to storageIndex
+    if (storageIndex_[index] < 0)
+    {
+        fetchStateForDataBase(index);
+    }
+    else
+    {
+        updateStorageUsageList(index);
+    }
+
+    return surfaceScalarFieldList_[fieldI][storageIndex_[index]];
+}
+
+
+
+
+SymmetricSquareMatrix<scalar>& gerhardsRecModel::recurrenceMatrix()
+{
+    return recurrenceMatrix_;
+}
+
+
+const HashTable<label,word>& gerhardsRecModel::timeIndexList() const
+{
+    return timeIndexList_;
+}
+
+
+label gerhardsRecModel::lowerSeqLim() const
+{
+    return lowerSeqLim_; 
+}
+
+label gerhardsRecModel::upperSeqLim() const
+{
+    return upperSeqLim_; 
+}
+
+
+label gerhardsRecModel::numRecFields() const
+{
+    return numRecFields_;
+}
+
+label gerhardsRecModel::numDataBaseFields() const
+{
+    return numDataBaseFields_;
+}
+
+void gerhardsRecModel::updateRecFields()
+{
+    virtualTimeIndex=virtualTimeIndexNext;
+    virtualTimeIndexNext++;
+    if (virtualTimeIndexNext>sequenceEnd)
+    {
+        virtualTimeIndexListPos++;
+        sequenceStart=virtualTimeIndexList_[virtualTimeIndexListPos].first();
+        sequenceEnd=virtualTimeIndexList_[virtualTimeIndexListPos].second();
+        virtualTimeIndexNext=sequenceStart;
+    }
+
+    if (verbose_)
+    {
+        Info << "\nUpdating virtual time index to " << virtualTimeIndex << ".\n" << endl;  
+    }
+}
+
+
+void gerhardsRecModel::writeRecMatrix() const
+{
+    Info << nl << "Writing recurrence matrix to disk." << nl
+        << "Nr. of reads from disk : " << nrOfReadsFromDisk_
+        << "; compared to a theoretical minimum of "
+        << numRecFields() << " reads." << endl;
+    
+    OFstream matrixFile("recurrenceMatrix");
+    matrixFile << recurrenceMatrix_;
+}
+
+
+
+
+Switch gerhardsRecModel::checkSkipZero()
+{
+    if (skipZero_)
+    {
+        bool foundZero(false);
+        
+        forAll(timeDirs, i)
+        {
+            if (timeDirs[i].value() == 0)
+            {
+                foundZero = true;
+            }
+        }
+        
+        if (not foundZero)
+        {
+            skipZero_ = false;
+        }
+    }
+    
+    return skipZero_;
+}
+
+
+
+// tmp<surfaceScalarField> gerhardsRecModel::exportAveragedSurfaceScalarField(word fieldname, scalar threshold, label index)
+// {
+//     label timeIndex;
+//     if (index < 0)
+//     {
+//         timeIndex = virtualTimeIndex;
+//     }
+//     else
+//     {
+//         timeIndex = index; 
+//     }
+//     const label fieldI = getSurfaceScalarFieldIndex(fieldname, timeIndex);
+//     
+//     tmp<surfaceScalarField> tAveragedSurfaceScalarField(surfaceScalarFieldList_[fieldI][timeIndex]);
+//     
+//     label counter = 1;
+//     scalar recErr;
+//     label delay = 10;
+//     label lastMin = -1000;
+//     
+//     for(int runningTimeIndex = 1; runningTimeIndex < numRecFields_-1 ; runningTimeIndex++)
+//     {
+//         recErr = recurrenceMatrix_[timeIndex][runningTimeIndex];
+//         if(recErr > threshold) continue;
+//         if(recErr > recurrenceMatrix_[timeIndex][runningTimeIndex-1]) continue;
+//         if(recErr > recurrenceMatrix_[timeIndex][runningTimeIndex+1]) continue;
+//         if(abs(runningTimeIndex - timeIndex) < delay) continue;
+//         if(abs(runningTimeIndex - lastMin) < delay) continue;
+// 
+//         lastMin = runningTimeIndex;
+//         counter++;
+//         tAveragedSurfaceScalarField += surfaceScalarFieldList_[fieldI][runningTimeIndex];
+//     }
+//     
+//     tAveragedSurfaceScalarField /= counter;
+//     return tAveragedSurfaceScalarField;
+// }
+
+/*void gerhardsRecModel::exportAveragedVolVectorField(volVectorField& smoothfield, word fieldname, scalar threshold, label index) const
+{
+    label timeIndex;
+    if (index < 0)
+    {
+        timeIndex = virtualTimeIndex;
+    }
+    else
+    {
+        timeIndex = index; 
+    }
+    const label fieldI = getVolVectorFieldIndex(fieldname, timeIndex);
+    
+    smoothfield = volVectorFieldList_[fieldI][timeIndex];
+     
+    label counter = 1;
+    scalar recErr;
+    label delay = 1;
+    label lastMin = -1000;
+
+    for(int runningTimeIndex = 0; runningTimeIndex < numRecFields_ ; runningTimeIndex++)
+    {
+        recErr = recurrenceMatrix_[timeIndex][runningTimeIndex];
+        if(recErr > threshold) continue;
+   //     if(recErr > recurrenceMatrix_[timeIndex][runningTimeIndex-1]) continue;
+   //     if(recErr > recurrenceMatrix_[timeIndex][runningTimeIndex+1]) continue;
+        if(abs(runningTimeIndex - timeIndex) < delay) continue;
+        if(abs(runningTimeIndex - lastMin) < delay) continue;
+
+        lastMin = runningTimeIndex;
+        counter++;
+        smoothfield += volVectorFieldList_[fieldI][runningTimeIndex];
+    }
+    Info << "time index = " << index << ", counter = " << counter << endl;
+    smoothfield /= counter;
+}
+*/
+
+
+
+
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModel.H b/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModel.H
new file mode 100644
index 00000000..606be634
--- /dev/null
+++ b/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModel.H
@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+    
+Class
+    Foam::gerhardsRecModel
+
+Description
+    A recurrence model for a dataBase with a smaller, or equal number of slots 
+    than there are snapshots on disk.
+    
+    In the case of less slots than snapshots, we follow the Least Recently Used
+    page replacement algorithm described in computer science literature.
+    A label list is used to keep track of which slot was most recently used,
+    the least recently used slot is then vacated if a new snapshot needs to be 
+    read from disk.
+    
+    Reference:
+    \verbatim
+        "Modern Operating Systems"
+        Andrew S. Tannenbaum,
+        Prentice Hall, 1992
+    \endverbatim
+
+SourceFiles
+    gerhardsRecModelI.H
+    gerhardsRecModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef gerhardsRecModel_H
+#define gerhardsRecModel_H
+
+#include "recModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class gerhardsRecModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class gerhardsRecModel
+:
+    public recModel
+{
+protected:
+  
+    dictionary propsDict_;
+    word dataBaseName_;
+    Foam::Time recTime;
+    instantList timeDirs;
+    Switch skipZero_;
+    Switch checkSkipZero_;
+    label numRecFields_;
+    
+    const Switch verboseVerbose_;
+    const Switch verboseVacate_;
+
+    // matrix that contains the recurrence ERROR
+    SymmetricSquareMatrix<scalar> recurrenceMatrix_;
+    
+    // create a data structure for the time indices
+    //	constant will not be contained
+    //	runTimeIndex -> continuousIndex
+    HashTable<label,word> timeIndexList_;
+    
+    HashTable<word, label> revTimeIndexList_;
+    
+    // create a data structure for the time values
+    //	constant will not be contained
+    //	continuousIndex -> time.value()
+    HashTable<label,scalar> timeValueList_;
+    
+    label contTimeIndex;
+    
+    const label lowerSeqLim_;
+    const label upperSeqLim_;
+    
+    scalar checkTimeStep();
+    
+    inline label getVolScalarFieldIndex(word, label) const;
+    inline label getVolVectorFieldIndex(word, label) const;
+    inline label getSurfaceScalarFieldIndex(word, label) const;
+  
+    void readFieldSeries();
+    
+    void readTimeSeries();
+    
+    void fetchStateForDataBase(label);
+    
+    void updateStorageUsageList(label);
+    
+    void readNewSnapshot(label, label);
+    
+    Switch checkSkipZero();
+
+public:
+
+    //- Runtime type information
+    TypeName("gerhardsRecModel");
+
+
+    // Constructors
+
+        //- Construct from components
+        gerhardsRecModel
+        (
+            const dictionary& dict,
+            recBase& base
+        );
+
+    // Destructor
+
+        ~gerhardsRecModel();
+
+    void exportVolScalarField(word, volScalarField&);
+    void exportVolVectorField(word, volVectorField&);
+    void exportSurfaceScalarField(word, surfaceScalarField&);
+    
+    const volScalarField& exportVolScalarField(word, label);
+    const volVectorField& exportVolVectorField(word, label);
+    const surfaceScalarField& exportSurfaceScalarField(word, label);
+    
+  //  tmp<surfaceScalarField> exportAveragedSurfaceScalarField(word, scalar, label index = -1);
+    void exportAveragedVolVectorField(volVectorField&, word, scalar, label index = -1) const;
+   
+    SymmetricSquareMatrix<scalar>& recurrenceMatrix();
+    
+    const HashTable<label,word>& timeIndexList() const;
+    
+    label lowerSeqLim() const;
+    label upperSeqLim() const;
+    
+    label numRecFields() const;
+    label numDataBaseFields() const;
+    
+    void updateRecFields();
+    
+    void writeRecMatrix() const;
+	
+private:
+    
+    /*
+        Number of fields kept in the dataBase, this number is smaller, or equal, as the 
+        number of snapshots.
+    */
+    const label numDataBaseFields_;
+    
+    /*
+        As the data base holds less states than there are snapshots, we need to translate
+        the index of the snapshot to an appropriate index of the dataBase.
+    */
+    List<label> storageIndex_;
+    
+    /*
+        As we repeatedly need to replace fields in the dataBase with new snapshots from 
+        disk, we need a way to judge, which slot to vacate for the new snapshot.
+        Here, we follow the lead of paging algorithms used by OS'es in memory management.
+        
+        The storageUsageList_ represents the recentness of a slot's use; with 1 for the 
+        most recently used element, and higher values for less recently used ones.
+    */
+    List<label> storageUsageList_;
+    
+    // this class does not read all the fields
+    List<PtrList<volScalarField> > volScalarFieldList_;
+    List<PtrList<volVectorField> > volVectorFieldList_;
+    List<PtrList<surfaceScalarField> > surfaceScalarFieldList_;
+    
+    // book-keeping
+    label nrOfReadsFromDisk_;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "gerhardsRecModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModelI.H b/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModelI.H
new file mode 100644
index 00000000..69081502
--- /dev/null
+++ b/src/recurrence/recModel/gerhardsRecModel/gerhardsRecModelI.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::gerhardsRecModel::getVolScalarFieldIndex(word fieldname, label index) const
+{
+    if(index < 0 || index > numRecFields_ - 1)
+    {
+        FatalError<<"gerhardsRecModel: Index out of bounds for volScalarField with name " 
+            << fieldname << abort(FatalError);
+    }
+    
+    label fieldI(-1);
+    
+    for(int i=0; i<volScalarFieldNames_.size(); i++)
+    {
+        if(volScalarFieldNames_[i].match(fieldname))
+	    {
+	        fieldI = i;
+	    }
+	}
+	
+	if (fieldI == -1)
+	{
+	    FatalError<<"gerhardsRecModel: Could not find volScalarField with name " 
+	        << fieldname << abort(FatalError);
+	}
+	
+	return fieldI;
+}
+
+inline Foam::label Foam::gerhardsRecModel::getVolVectorFieldIndex(word fieldname, label index) const
+{
+    if(index < 0 || index > numRecFields_ - 1)
+    {
+        FatalError<<"gerhardsRecModel: Index out of bounds for volVectorField with name " 
+            << fieldname << abort(FatalError);
+    }
+    
+    label fieldI(-1);
+    
+    for(int i=0; i<volVectorFieldNames_.size(); i++)
+    {
+        if(volVectorFieldNames_[i].match(fieldname))
+	    {
+	        fieldI = i;
+	    }
+	}
+	
+	if (fieldI == -1)
+	{
+	    FatalError<<"gerhardsRecModel: Could not find volVectorField with name " 
+	        << fieldname << abort(FatalError);
+	}
+	
+	return fieldI;
+}
+
+inline Foam::label Foam::gerhardsRecModel::getSurfaceScalarFieldIndex(word fieldname, label index) const
+{
+    if(index < 0 || index > numRecFields_ - 1)
+    {
+        FatalError<<"gerhardsRecModel: Index out of bounds for surfaceScalarField with name " 
+            << fieldname << abort(FatalError);
+    }
+
+    label fieldI(-1);
+
+    for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
+    {
+        if(surfaceScalarFieldNames_[i].match(fieldname))
+        {
+            fieldI = i;
+        }
+    }
+
+    if (fieldI == -1)
+    {
+        FatalError<<"gerhardsRecModel: Could not find surfaceScalarField with name " 
+            << fieldname << abort(FatalError);
+    }
+
+    return fieldI;
+}
diff --git a/src/recurrence/recModel/recModel/recModel.C b/src/recurrence/recModel/recModel/recModel.C
index b5bd17b6..e469b46e 100644
--- a/src/recurrence/recModel/recModel/recModel.C
+++ b/src/recurrence/recModel/recModel/recModel.C
@@ -63,7 +63,7 @@ recModel::recModel
             IOobject::NO_WRITE
         )
     ),
-    verbose_(false),
+    verbose_(dict.lookupOrDefault<Switch>("verbose", false)),
     volScalarFieldNames_(recProperties_.lookup("volScalarFields")),
     volVectorFieldNames_(recProperties_.lookup("volVectorFields")),
     surfaceScalarFieldNames_(recProperties_.lookup("surfaceScalarFields")),
@@ -80,7 +80,6 @@ recModel::recModel
     virtualTimeIndexList_(0),
     virtualTimeIndexListPos(0)
 {
-    if (recProperties_.found("verbose")) verbose_=true;
     recTimeStep_ = -1.0;
     totRecSteps_ = -1;
 }
diff --git a/src/recurrence/recModel/recModel/recModel.H b/src/recurrence/recModel/recModel/recModel.H
index 260ed995..64061da9 100644
--- a/src/recurrence/recModel/recModel/recModel.H
+++ b/src/recurrence/recModel/recModel/recModel.H
@@ -125,16 +125,24 @@ public:
 
     // Member Functions
       
-    virtual void exportVolScalarField(word, volScalarField&) const = 0;
+    /*virtual void exportVolScalarField(word, volScalarField&) const = 0;
     virtual void exportVolVectorField(word, volVectorField&) const = 0;
     virtual void exportSurfaceScalarField(word, surfaceScalarField&) const = 0;
     
     virtual const volScalarField& exportVolScalarField(word, label) const = 0;
     virtual const volVectorField& exportVolVectorField(word, label) const = 0;
-    virtual const surfaceScalarField& exportSurfaceScalarField(word, label) const = 0;
+    virtual const surfaceScalarField& exportSurfaceScalarField(word, label) const = 0;*/
+    
+    virtual void exportVolScalarField(word, volScalarField&)  = 0;
+    virtual void exportVolVectorField(word, volVectorField&)  = 0;
+    virtual void exportSurfaceScalarField(word, surfaceScalarField&)  = 0;
+
+    virtual const volScalarField& exportVolScalarField(word, label)  = 0;
+    virtual const volVectorField& exportVolVectorField(word, label)  = 0;
+    virtual const surfaceScalarField& exportSurfaceScalarField(word, label)  = 0;
     
   //  virtual tmp<surfaceScalarField> exportAveragedSurfaceScalarField(word, scalar, label index = -1) = 0;
-    virtual void exportAveragedVolVectorField(volVectorField&, word, scalar, label index = -1) const = 0;
+    //virtual void exportAveragedVolVectorField(volVectorField&, word, scalar, label index = -1) const = 0;
     
     virtual SymmetricSquareMatrix<scalar>& recurrenceMatrix() = 0;
     
@@ -146,6 +154,7 @@ public:
     virtual label upperSeqLim() const = 0;
     
     virtual label numRecFields() const = 0;
+    virtual label numDataBaseFields() const = 0;
     
     label totRecSteps() const;
     
diff --git a/src/recurrence/recModel/standardRecModel/standardRecModel.C b/src/recurrence/recModel/standardRecModel/standardRecModel.C
index 47a5b5dc..5e3d6795 100644
--- a/src/recurrence/recModel/standardRecModel/standardRecModel.C
+++ b/src/recurrence/recModel/standardRecModel/standardRecModel.C
@@ -26,7 +26,6 @@ License
 #include "Random.H"
 #include "standardRecModel.H"
 #include "addToRunTimeSelectionTable.H"
-#include <mpi.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,12 +55,14 @@ standardRecModel::standardRecModel
 :
     recModel(dict,base),
     propsDict_(dict.subDict(typeName + "Props")),
+    //dataBaseName_("dataBase"),
     dataBaseName_(propsDict_.lookupOrDefault<word>("dataBase", word("dataBase"))),
+    //recTime("dataBase", "", "../system", "../constant", false),
     recTime(fileName(dataBaseName_), "", "../system", "../constant", false),
     timeDirs(recTime.times()),
     skipZero_(propsDict_.lookupOrDefault<Switch>("skipZero", Switch(false))),
     numRecFields_(skipZero_ ? label(timeDirs.size())-1 : label(timeDirs.size())),
-    recurrenceMatrix_(numRecFields_,scalar(0.0)),
+    recurrenceMatrix_(numRecFields_,scalar(-1.0)),
     timeIndexList_(numRecFields_-1),
     timeValueList_(numRecFields_-1),
     contTimeIndex(0),
@@ -71,38 +72,39 @@ standardRecModel::standardRecModel
     volScalarFieldList_(volScalarFieldNames_.size()),
     volVectorFieldList_(volVectorFieldNames_.size()),
     surfaceScalarFieldList_(surfaceScalarFieldNames_.size())
-{
-  
-  if (verbose_)
+{   
+    if (verbose_)
     {
-    	// be informative on properties of the "recTime" Time-object
-    	Info << "recTime.rootPath() " << recTime.rootPath() << endl;
-    	Info << "recTime.caseName() " << recTime.caseName() << endl;
-    	Info << "recTime.path() " << recTime.path() << endl;
-    	Info << "recTime.timePath() " << recTime.timePath() << endl;
-    	Info << "recTime.timeName() " << recTime.timeName() << endl;
-    	Info << "timeDirs " << timeDirs << endl;
-    	Info << "ignoring 0 directory: " << skipZero_ << endl;
+        // be informative on properties of the "recTime" Time-object
+        Info << "recTime.rootPath() " << recTime.rootPath() << endl;
+        Info << "recTime.caseName() " << recTime.caseName() << endl;
+        Info << "recTime.path() " << recTime.path() << endl;
+        Info << "recTime.timePath() " << recTime.timePath() << endl;
+        Info << "recTime.timeName() " << recTime.timeName() << endl;
+        Info << "timeDirs " << timeDirs << endl;
+        Info << "ignoring 0 directory: " << skipZero_ << endl;
     }
     readTimeSeries();
-  
+
     recTimeStep_ = checkTimeStep();
-    totRecSteps_ = 1+static_cast<label> ((endTime_-startTime_) / recTimeStep_);
-  
+    totRecSteps_ = 1 + static_cast<label>( (endTime_-startTime_) / recTimeStep_ );
+
     for(int i=0; i<volScalarFieldNames_.size(); i++)
+    {
         volScalarFieldList_[i].setSize(numRecFields_);
-    
+    }
+
     for(int i=0; i<volVectorFieldNames_.size(); i++)
+    {
         volVectorFieldList_[i].setSize(numRecFields_);
-    
+    }
+
     for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
+    {
         surfaceScalarFieldList_[i].setSize(numRecFields_);
-       
-
-    
- //   setRecFields();
-    
+    }
 
+    //   setRecFields();
 }
 
 
@@ -126,9 +128,9 @@ scalar standardRecModel::checkTimeStep()
     
     if (timeIndexList_.size() == 1)
     {
-        return 1e10;
+        return 1.e10;
     }
-
+    
     forAll(timeValueList_, i)
     {
     	// skip zero
@@ -200,11 +202,11 @@ void standardRecModel::readFieldSeries()
     for (instantList::iterator it=timeDirs.begin(); it != timeDirs.end(); ++it)
     {
         if(counter >= 0.1 * percentage * size)
-	{
-	    Info << "\t" << 10 * percentage << " \% done" << endl;
-	    percentage++;
-	}
-	counter++;
+	    {
+	        Info << "\t" << 10 * percentage << " \% done" << endl;
+	        percentage++;
+	    }
+	    counter++;
       
         // set time
         recTime.setTime(*it, it->value());
@@ -232,9 +234,10 @@ void standardRecModel::readFieldSeries()
         }
         
         for(int i=0; i<volScalarFieldNames_.size(); i++)
+        {
             volScalarFieldList_[i].set
             (
-	        timeIndexList_(recTime.timeName()),
+	            timeIndexList_(recTime.timeName()),
                 new volScalarField
                 (
                     IOobject
@@ -247,12 +250,14 @@ void standardRecModel::readFieldSeries()
                     ),
                     base_.mesh()
                 )
-	    );
+	        );
+	    }
 	    
-	for(int i=0; i<volVectorFieldNames_.size(); i++)
+	    for(int i=0; i<volVectorFieldNames_.size(); i++)
+	    {
             volVectorFieldList_[i].set
             (
-	        timeIndexList_(recTime.timeName()),
+	            timeIndexList_(recTime.timeName()),
                 new volVectorField
                 (
                     IOobject
@@ -265,12 +270,14 @@ void standardRecModel::readFieldSeries()
                     ),
                     base_.mesh()
                 )
-	    );
+	        );
+	    }
 	    
-	for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
+	    for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
+	    {
             surfaceScalarFieldList_[i].set
             (
-	        timeIndexList_(recTime.timeName()),
+                timeIndexList_(recTime.timeName()),
                 new surfaceScalarField
                 (
                     IOobject
@@ -283,7 +290,8 @@ void standardRecModel::readFieldSeries()
                     ),
                     base_.mesh()
                 )
-            );    
+            );
+        } 
     }
     Info << "Reading fields done" << endl;
 }
@@ -292,7 +300,7 @@ void standardRecModel::readFieldSeries()
 void standardRecModel::readTimeSeries()
 {
     bool firsttime = true;
-  // fill the data structure for the time indices
+    // fill the data structure for the time indices
     for (instantList::iterator it=timeDirs.begin(); it != timeDirs.end(); ++it)
     {
     	// set run-time
@@ -317,11 +325,11 @@ void standardRecModel::readTimeSeries()
         }
         
         if (firsttime)
-	{
-	    firsttime = false;
-	    recStartTime_ = recTime.value();
-	}
-	recEndTime_ = recTime.value();
+        {
+            firsttime = false;
+            recStartTime_ = recTime.value();
+        }
+	    recEndTime_ = recTime.value();
         
         // insert the time name into the hash-table with a continuous second index
         timeIndexList_.insert(recTime.timeName(), contTimeIndex);
@@ -348,48 +356,48 @@ void standardRecModel::readTimeSeries()
 
     if (verbose_)
     {
-    	Info << endl;
-    	Info << "Found " << label(timeDirs.size()) << " time folders" << endl;
+        Info << endl;
+        Info << "Found " << label(timeDirs.size()) << " time folders" << endl;
         Info << "Found " << label(timeIndexList_.size()) << " time steps" << endl;
-	Info << "database start time = " << recStartTime_ << endl;
-	Info << "database end time = " << recEndTime_ << endl;
+        Info << "database start time = " << recStartTime_ << endl;
+        Info << "database end time = " << recEndTime_ << endl;
     }
 }
 
 
-void standardRecModel::exportVolScalarField(word fieldname, volScalarField& field) const
+void standardRecModel::exportVolScalarField(word fieldname, volScalarField& field)
 {
     field = exportVolScalarField(fieldname, virtualTimeIndex);
 }
 
 
-void standardRecModel::exportVolVectorField(word fieldname, volVectorField& field) const
+void standardRecModel::exportVolVectorField(word fieldname, volVectorField& field)
 {
     field = exportVolVectorField(fieldname, virtualTimeIndex);
 }
 
 
-void standardRecModel::exportSurfaceScalarField(word fieldname, surfaceScalarField& field) const
+void standardRecModel::exportSurfaceScalarField(word fieldname, surfaceScalarField& field)
 {
     field = exportSurfaceScalarField(fieldname, virtualTimeIndex);
 }
 
 
-const volScalarField& standardRecModel::exportVolScalarField(word fieldname, label index) const
+const volScalarField& standardRecModel::exportVolScalarField(word fieldname, label index)
 {
     const label fieldI = getVolScalarFieldIndex(fieldname, index);
     
     return volScalarFieldList_[fieldI][index];
 }
 
-const volVectorField& standardRecModel::exportVolVectorField(word fieldname, label index) const
+const volVectorField& standardRecModel::exportVolVectorField(word fieldname, label index)
 {
     const label fieldI = getVolVectorFieldIndex(fieldname, index);
     
     return volVectorFieldList_[fieldI][index];
 }
 
-const surfaceScalarField& standardRecModel::exportSurfaceScalarField(word fieldname, label index) const
+const surfaceScalarField& standardRecModel::exportSurfaceScalarField(word fieldname, label index)
 {
     const label fieldI = getSurfaceScalarFieldIndex(fieldname, index);
     
@@ -425,21 +433,36 @@ label standardRecModel::numRecFields() const
     return numRecFields_;
 }
 
+label standardRecModel::numDataBaseFields() const
+{
+    return numRecFields_;
+}
+
 
 void standardRecModel::updateRecFields()
 {
-  virtualTimeIndex=virtualTimeIndexNext;
-  virtualTimeIndexNext++;
-  if (virtualTimeIndexNext>sequenceEnd)
-  {
-    virtualTimeIndexListPos++;
-    sequenceStart=virtualTimeIndexList_[virtualTimeIndexListPos].first();
-    sequenceEnd=virtualTimeIndexList_[virtualTimeIndexListPos].second();
-    virtualTimeIndexNext=sequenceStart;
-  }
-  
-  if (verbose_)
-      Info << "\nUpdating virtual time index to " << virtualTimeIndex << ".\n" << endl;  
+    virtualTimeIndex = virtualTimeIndexNext;
+    virtualTimeIndexNext++;
+    
+    if (virtualTimeIndexNext > sequenceEnd)
+    {
+        virtualTimeIndexListPos++; // this is problematic with noPath
+        
+        sequenceStart = virtualTimeIndexList_[virtualTimeIndexListPos].first();
+        sequenceEnd = virtualTimeIndexList_[virtualTimeIndexListPos].second();
+        
+        virtualTimeIndexNext = sequenceStart;
+        
+        if (verbose_)
+        {
+            Info << " new sequence (start/end) : " << sequenceStart << " / " << sequenceEnd << endl;
+        }
+    }
+
+    if (verbose_)
+    {
+        Info << "\nUpdating virtual time index to " << virtualTimeIndex << ".\n" << endl;
+    }
 }
 
 
diff --git a/src/recurrence/recModel/standardRecModel/standardRecModel.H b/src/recurrence/recModel/standardRecModel/standardRecModel.H
index 3ec05bfa..2ced11f8 100644
--- a/src/recurrence/recModel/standardRecModel/standardRecModel.H
+++ b/src/recurrence/recModel/standardRecModel/standardRecModel.H
@@ -20,6 +20,14 @@ License
 
     You should have received a copy of the GNU General Public License
     along with CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Standard Recurrence Model
+    
+    This class sets limits for the length of a recurrence-sequence. As of now, 
+    such a sequence is between one 20-th and 5-th of the total number of data 
+    base entries.
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef standardRecModel_H
@@ -96,13 +104,13 @@ public:
 
         ~standardRecModel();
 
-    void exportVolScalarField(word, volScalarField&) const;
-    void exportVolVectorField(word, volVectorField&) const;
-    void exportSurfaceScalarField(word, surfaceScalarField&) const;
+    void exportVolScalarField(word, volScalarField&);
+    void exportVolVectorField(word, volVectorField&);
+    void exportSurfaceScalarField(word, surfaceScalarField&);
     
-    const volScalarField& exportVolScalarField(word, label) const;
-    const volVectorField& exportVolVectorField(word, label) const;
-    const surfaceScalarField& exportSurfaceScalarField(word, label) const;
+    const volScalarField& exportVolScalarField(word, label);
+    const volVectorField& exportVolVectorField(word, label);
+    const surfaceScalarField& exportSurfaceScalarField(word, label);
     
   //  tmp<surfaceScalarField> exportAveragedSurfaceScalarField(word, scalar, label index = -1);
     void exportAveragedVolVectorField(volVectorField&, word, scalar, label index = -1) const;
@@ -115,6 +123,7 @@ public:
     label upperSeqLim() const;
     
     label numRecFields() const;
+    label numDataBaseFields() const;
     
     void updateRecFields();
     
@@ -125,6 +134,7 @@ private:
     List<PtrList<volScalarField> > volScalarFieldList_;
     List<PtrList<volVectorField> > volVectorFieldList_;
     List<PtrList<surfaceScalarField> > surfaceScalarFieldList_;
+    
 };
 
 
diff --git a/src/recurrence/recModel/standardRecModel/standardRecModelI.H b/src/recurrence/recModel/standardRecModel/standardRecModelI.H
index 05e42f64..ce5d4cec 100644
--- a/src/recurrence/recModel/standardRecModel/standardRecModelI.H
+++ b/src/recurrence/recModel/standardRecModel/standardRecModelI.H
@@ -85,22 +85,22 @@ inline Foam::label Foam::standardRecModel::getSurfaceScalarFieldIndex(word field
         FatalError<<"standardRecModel: Index out of bounds for surfaceScalarField with name " 
             << fieldname << abort(FatalError);
     }
-    
+
     label fieldI(-1);
-    
+
     for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
     {
         if(surfaceScalarFieldNames_[i].match(fieldname))
-	    {
-	        fieldI = i;
-	    }
-	}
-	
-	if (fieldI == -1)
-	{
-	    FatalError<<"standardRecModel: Could not find surfaceScalarField with name " 
-	        << fieldname << abort(FatalError);
-	}
-	
-	return fieldI;
+        {
+            fieldI = i;
+        }
+    }
+
+    if (fieldI == -1)
+    {
+        FatalError<<"standardRecModel: Could not find surfaceScalarField with name " 
+            << fieldname << abort(FatalError);
+    }
+
+    return fieldI;
 }
diff --git a/src/recurrence/recNorm/recNorm/recNorm.C b/src/recurrence/recNorm/recNorm/recNorm.C
index 61e58d97..01ec2f7c 100644
--- a/src/recurrence/recNorm/recNorm/recNorm.C
+++ b/src/recurrence/recNorm/recNorm/recNorm.C
@@ -24,7 +24,7 @@ License
 
 #include "error.H"
 #include "recNorm.H"
-#include <unistd.h>
+//#include <unistd.h>
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,9 +52,8 @@ recNorm::recNorm
 :
     base_(base),
     recProperties_(dict),
-    verbose_(false)
+    verbose_(dict.lookupOrDefault<Switch>("verbose", false))
 {
-    if (recProperties_.found("verbose")) verbose_=true;
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/recurrence/recNorm/recNorm/recNorm.H b/src/recurrence/recNorm/recNorm/recNorm.H
index c34ef568..e5ac6a83 100644
--- a/src/recurrence/recNorm/recNorm/recNorm.H
+++ b/src/recurrence/recNorm/recNorm/recNorm.H
@@ -26,7 +26,7 @@ License
 #define recNorm_H
 
 #include "recBase.H"
-#include "fvCFD.H"
+//#include "fvCFD.H"
 #include "HashTable.H"
 // #include "labelPair.H"
 #include "OFstream.H"
diff --git a/src/recurrence/recNorm/sqrDiffNorm/sqrDiffNorm.C b/src/recurrence/recNorm/sqrDiffNorm/sqrDiffNorm.C
index 1fadb8e2..057f6174 100644
--- a/src/recurrence/recNorm/sqrDiffNorm/sqrDiffNorm.C
+++ b/src/recurrence/recNorm/sqrDiffNorm/sqrDiffNorm.C
@@ -69,71 +69,141 @@ sqrDiffNorm::~sqrDiffNorm()
 
 void sqrDiffNorm::computeRecMatrix()
 {
-    Info<< "\nsqrDiffNorm: computing recurrence matrix\n" << endl;
-    
+    Info << nl << type() << ": computing recurrence matrix\n" << endl;
+
     const HashTable<label,word>& timeIndexList( base_.recM().timeIndexList() );
     SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );
-    
+
     scalar normIJ(0.0);
     scalar maxNormIJ(0.0);
-    
-    label size = timeIndexList.size();
+
+    /*
+        total number of computed elements: total number of matrix entries
+            minus the main diagonal entries, divided by two, 
+            since the matrix is symmetric
+    */
+    label size = (timeIndexList.size()*(timeIndexList.size()-1))/2;
     label counter = 0;
     label percentage = 0;
-
-    // perform un-normalized calculation for lower half of the recurrence matrix
-    forAll(timeIndexList, ti)
+    
+    
+    label N(this->base_.recM().numRecFields());
+    label M(this->base_.recM().numDataBaseFields());
+    
+    if (verbose_)
     {
-        if(counter >= 0.1 * percentage * size)
-	{
-	    Info << "\t" << 10 * percentage << " \% done" << endl;
-	    percentage++;
-	}
-	counter++;
-    	forAll(timeIndexList, tj)
-    	{
-	        if(verbose_)
-		    Info<<"\n Doing calculation for element " << ti << " " << tj << "\n" << endl;
-		
-    		// skip main diagonal and upper half
-    		if (ti >= tj)
-    		{
-    			recurrenceMatrix[ti][tj] = 0;
-    			continue;
-    		}
-    		
-    		// compute elements
-    		if (fieldType_ == "volScalarField")
-		    normIJ = normVSF(ti,tj);
-		else if (fieldType_ == "volVectorField")
-		    normIJ = normVVF(ti,tj);
-		else if (fieldType_ == "surfaceScalarField")
-		    normIJ = normSSF(ti,tj);
-		else
-		    FatalError<<"sqrDiffNorm: Unknown field type " << fieldType_ << abort(FatalError);
-		
-    		recurrenceMatrix[ti][tj] = normIJ;
-		
-		if (normIJ > maxNormIJ)
-		    maxNormIJ = normIJ;
-    	}
+        Info << " N = " << N << ",  M = " << M << endl;
     }
     
+    label ti(0);
+    label tj(0);
+    label tmp(-1);
+    
+    for (int j=0; j<=N/(M-1); j++)
+    {
+        for (int i=0; i<timeIndexList.size(); i++)
+        {
+            if (verbose_)
+            {
+                Info << " i = " << i << ",  j = " << j << endl;
+            }
+            
+            if(counter >= 0.1 * percentage * size)
+	        {
+	            Info << "\t" << 10 * percentage << " \% done" << endl;
+	            percentage++;
+	        }
+	        
+            
+            for (int k=i; k<i+(M-1); k++)
+            {
+                ti = i;
+                tj = j*(M-1) + k;
+                
+                if (ti > tj)
+                {
+                    tmp = ti;
+                    ti = tj;
+                    tj = tmp;
+                }
+                
+                // skip coordinates outside the recurrence space                
+                if (ti >= N or tj >= N)
+                {
+                    continue;
+                }
+                
+                // start
+                // skip main diagonal and upper half
+                if (ti >= tj)
+                {
+                    recurrenceMatrix[ti][tj] = 0;
+                    continue;
+                }
+                
+                if (verbose_)
+                {
+                    Info << " Doing calculation for element " 
+                        << ti << " " << tj << endl;
+                }
+                
+                counter++;
+
+                // compute elements
+                if (fieldType_ == "volScalarField")
+                {
+                    normIJ = normVSF(ti,tj);
+                }
+                else if (fieldType_ == "volVectorField")
+                {
+                    normIJ = normVVF(ti,tj);
+                }
+                else if (fieldType_ == "surfaceScalarField")
+                {
+                    normIJ = normSSF(ti,tj);
+                }
+                else
+                {
+                    FatalError
+                        << "sqrDiffNorm: Unknown field type " << fieldType_ 
+                        << abort(FatalError);
+                }
+
+                recurrenceMatrix[ti][tj] = normIJ;
+
+                if (normIJ > maxNormIJ)
+                {
+                    maxNormIJ = normIJ;
+                }
+                // end
+            }
+        }
+    }
+    
+
+
     // normalize matrix and copy lower to upper half
     if(normConstant_ > 0.0) maxNormIJ = normConstant_;
+    
     forAll(timeIndexList, ti)
     {
-    	forAll(timeIndexList, tj)
-    	{
-	    if (ti >= tj)
-	        continue;
-	    
-	    recurrenceMatrix[ti][tj] /= maxNormIJ;
-	    recurrenceMatrix[tj][ti] = recurrenceMatrix[ti][tj];
-	}
+        forAll(timeIndexList, tj)
+        {
+            if (ti >= tj) continue;
+            
+            if (recurrenceMatrix[ti][tj] < 0)
+            {
+                FatalErrorInFunction << "Error in computation of recurrence matrix!"
+                    << nl << "Negative elements encountered. This should not happen!"
+                    << abort(FatalError);
+            }
+
+            recurrenceMatrix[ti][tj] /= maxNormIJ;
+            recurrenceMatrix[tj][ti] = recurrenceMatrix[ti][tj];
+        }
     }
-    
-    Info<< "\nComputing recurrence matrix done\n" << endl;
+
+    Info << "\nComputing recurrence matrix done\n" << endl;
 }
 
 scalar sqrDiffNorm::normVSF(label ti, label tj)
@@ -141,6 +211,7 @@ scalar sqrDiffNorm::normVSF(label ti, label tj)
     const volScalarField& t1( base_.recM().exportVolScalarField(fieldName_,ti) );
     const volScalarField& t2( base_.recM().exportVolScalarField(fieldName_,tj) );
     dimensionedScalar tNorm( fvc::domainIntegrate( sqr( t1 - t2 ) ) );
+    
     return tNorm.value();
 }
 
@@ -149,6 +220,7 @@ scalar sqrDiffNorm::normVVF(label ti, label tj)
     const volVectorField& t1( base_.recM().exportVolVectorField(fieldName_,ti) );
     const volVectorField& t2( base_.recM().exportVolVectorField(fieldName_,tj) );
     dimensionedScalar tNorm( fvc::domainIntegrate( magSqr( t1 - t2 ) ) );
+    
     return tNorm.value();
 }
 
@@ -158,6 +230,7 @@ scalar sqrDiffNorm::normSSF(label ti, label tj)
     const surfaceScalarField& t2( base_.recM().exportSurfaceScalarField(fieldName_,tj) );
     volVectorField t12 (fvc::reconstruct( t1-t2 ) );
     dimensionedScalar tNorm( fvc::domainIntegrate( magSqr( t12 ) ) );
+    
     return tNorm.value();
 }
 
diff --git a/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C b/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C
index 210051b0..62227e74 100644
--- a/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C
+++ b/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C
@@ -125,27 +125,26 @@ multiIntervalPath::~multiIntervalPath()
 void multiIntervalPath::getRecPath()
 {
     label numRecIntervals = 0;
-    int root = 0;
-    int rank;
-    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
     
-    if(rank==root)
+    if(Pstream::master())
     {
-      computeRecPath();
-      numRecIntervals=virtualTimeIndexList_.size();
+        computeRecPath();
+        numRecIntervals = virtualTimeIndexList_.size();
     }
     
-    MPI_Bcast(&numRecIntervals,sizeof(numRecIntervals),MPI_BYTE,root,MPI_COMM_WORLD);
+    Pstream::scatter(numRecIntervals);
     
-    if(rank!=root)
-      virtualTimeIndexList_.setSize(numRecIntervals);
+    if(not Pstream::master())
+    {
+        virtualTimeIndexList_.setSize(numRecIntervals);
+    }
     
-    int listSizeBytes = 2*sizeof(numRecIntervals)*numRecIntervals;
-    
-    MPI_Bcast(&virtualTimeIndexList_[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
+    Pstream::scatter(virtualTimeIndexList_);
     
     if(verbose_)
-      Info << "\nRecurrence path communicated to all processors.\n" << endl;
+    {
+        Info << "\nRecurrence path communicated to all processors.\n" << endl;
+    }
 }
 
 
diff --git a/src/recurrence/recPath/noPath/noPath.C b/src/recurrence/recPath/noPath/noPath.C
index 933cc64f..9b025421 100644
--- a/src/recurrence/recPath/noPath/noPath.C
+++ b/src/recurrence/recPath/noPath/noPath.C
@@ -68,27 +68,26 @@ noPath::~noPath()
 void noPath::getRecPath()
 {
     label numRecIntervals = 0;
-    int root = 0;
-    int rank;
-    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
     
-    if(rank==root)
+    if(Pstream::master())
     {
-      computeRecPath();
-      numRecIntervals=virtualTimeIndexList_.size();
+        computeRecPath();
+        numRecIntervals=virtualTimeIndexList_.size();
     }
     
-    MPI_Bcast(&numRecIntervals,sizeof(numRecIntervals),MPI_BYTE,root,MPI_COMM_WORLD);
+    Pstream::scatter(numRecIntervals);
     
-    if(rank!=root)
-      virtualTimeIndexList_.setSize(numRecIntervals);
+    if(not Pstream::master())
+    {
+        virtualTimeIndexList_.setSize(numRecIntervals);
+    }
     
-    int listSizeBytes = 2*sizeof(numRecIntervals)*numRecIntervals;
-    
-    MPI_Bcast(&virtualTimeIndexList_[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
+    Pstream::scatter(virtualTimeIndexList_);
     
     if(verbose_)
-      Info << "\nRecurrence path communicated to all processors.\n" << endl;
+    {
+        Info << "\nRecurrence path communicated to all processors.\n" << endl;
+    }
 }
 
 
@@ -96,11 +95,16 @@ void noPath::computeRecPath()
 {
     labelPair seqStartEnd(0,1);
     virtualTimeIndexList_.append(seqStartEnd);
+    
+    if (verbose_)
+    {
+        Info << " virtualTimeIndexList_ : " << virtualTimeIndexList_ << endl;
+    }
 }
 
 label noPath::seqEnd(label seqStart, label & seqLength)
 {
-  return 0;
+    return 0;
 }
 
 void noPath::computeJumpVector()
diff --git a/src/recurrence/recPath/recPath/recPath.C b/src/recurrence/recPath/recPath/recPath.C
index 4070148a..c752b863 100644
--- a/src/recurrence/recPath/recPath/recPath.C
+++ b/src/recurrence/recPath/recPath/recPath.C
@@ -53,10 +53,9 @@ recPath::recPath
 :
     base_(base),
     recProperties_(dict),
-    verbose_(false),
+    verbose_(dict.lookupOrDefault<Switch>("verbose", false)),
     virtualTimeIndexList_( base_.recM().virtualTimeIndexList() )
 {
-    if (recProperties_.found("verbose")) verbose_=true;
 }
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
diff --git a/src/recurrence/recPath/recPath/recPath.H b/src/recurrence/recPath/recPath/recPath.H
index dfee7869..22f40e0c 100644
--- a/src/recurrence/recPath/recPath/recPath.H
+++ b/src/recurrence/recPath/recPath/recPath.H
@@ -30,7 +30,7 @@ License
 #include "HashTable.H"
 #include "labelPair.H"
 #include "OFstream.H"
-#include <mpi.h>
+//#include <mpi.h>
 
 
 namespace Foam
diff --git a/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C b/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C
index a838546e..8c480934 100644
--- a/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C
+++ b/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C
@@ -68,27 +68,26 @@ simpleRandomPath::~simpleRandomPath()
 void simpleRandomPath::getRecPath()
 {
     label numRecIntervals = 0;
-    int root = 0;
-    int rank;
-    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
     
-    if(rank==root)
+    if(Pstream::master())
     {
-      computeRecPath();
-      numRecIntervals=virtualTimeIndexList_.size();
+        computeRecPath();
+        numRecIntervals = virtualTimeIndexList_.size();
     }
     
-    MPI_Bcast(&numRecIntervals,sizeof(numRecIntervals),MPI_BYTE,root,MPI_COMM_WORLD);
+    Pstream::scatter(numRecIntervals);
     
-    if(rank!=root)
-      virtualTimeIndexList_.setSize(numRecIntervals);
+    if(not Pstream::master())
+    {
+        virtualTimeIndexList_.setSize(numRecIntervals);
+    }
     
-    int listSizeBytes = 2*sizeof(numRecIntervals)*numRecIntervals;
-    
-    MPI_Bcast(&virtualTimeIndexList_[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
+    Pstream::scatter(virtualTimeIndexList_);
     
     if(verbose_)
-      Info << "\nRecurrence path communicated to all processors.\n" << endl;
+    {
+        Info << "\nRecurrence path communicated to all processors.\n" << endl;
+    }
 }
 
 
@@ -99,21 +98,21 @@ void simpleRandomPath::computeRecPath()
     
     Random ranGen(osRandomInteger());
     
-    label virtualTimeIndex=0;
-    label recSteps=0;
-    label seqStart=0;
+    label virtualTimeIndex = 0;
+    label recSteps = 0;
+    label seqStart = 0;
     label lowerSeqLim( base_.recM().lowerSeqLim() );
     label upperSeqLim( base_.recM().upperSeqLim() );
-    label seqLength=ranGen.integer(lowerSeqLim, upperSeqLim);
+    label seqLength = ranGen.integer(lowerSeqLim, upperSeqLim);
     
-    virtualTimeIndex=seqEnd(seqStart,seqLength);
+    virtualTimeIndex = seqEnd(seqStart,seqLength);
     labelPair seqStartEnd(seqStart,virtualTimeIndex);
     virtualTimeIndexList_.append(seqStartEnd);
-    recSteps+=seqLength;
-
+    recSteps += seqLength;
+    
     SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );
     label numRecFields( base_.recM().numRecFields() );
-
+    
     if(base_.recM().totRecSteps() == 1)
     {
         Info<< "\nPrimitive recurrence path with one element.\n" << endl;
@@ -129,9 +128,9 @@ void simpleRandomPath::computeRecPath()
         // the new starting point of the next sequence
         if (virtualTimeIndex < numRecFields/2)
         {
-        	startLoop = numRecFields/2;
-        	endLoop = numRecFields-1 - upperSeqLim/2;  // avoid running into end of recurrence database too often
-		endLoop--; // start of next sequence one snapshot AFTER minimum position
+            startLoop = numRecFields/2;
+            endLoop = numRecFields-1 - upperSeqLim/2;  // avoid running into end of recurrence database too often
+            endLoop--; // start of next sequence one snapshot AFTER minimum position
         }
         else
         {
@@ -151,61 +150,73 @@ void simpleRandomPath::computeRecPath()
         }
         
         seqLength = ranGen.integer(lowerSeqLim, upperSeqLim);  
-        virtualTimeIndex=seqEnd(seqStart,seqLength);
+        virtualTimeIndex = seqEnd(seqStart,seqLength);
         labelPair seqStartEnd(seqStart,virtualTimeIndex);
         virtualTimeIndexList_.append(seqStartEnd);
-        recSteps+=seqLength;
+        recSteps += seqLength;
     }
+    
     Info<< "\nComputing recurrence path done\n" << endl;
+    
+    if (verbose_)
+    {
+        Info << " virtualTimeIndexList_ : " << virtualTimeIndexList_ << endl;
+    }
+    
     computeJumpVector();
 }
 
 label simpleRandomPath::seqEnd(label seqStart, label & seqLength)
 {
-  label numRecFields( base_.recM().numRecFields() );
-  if(seqStart+seqLength > numRecFields-1)
-      seqLength=numRecFields-1-seqStart;
-  return seqStart+seqLength;
+    label numRecFields( base_.recM().numRecFields() );
+    if(seqStart+seqLength > numRecFields-1)
+    {
+        seqLength=numRecFields-1-seqStart;
+    }
+    
+    return seqStart+seqLength;
 }
 
 void simpleRandomPath::computeJumpVector()
 {
     Info << "\nComputing recurrence jump vector\n" << endl;
-    
+
     OFstream jumpvec("rec_jump.dat");
     label numRecFields( base_.recM().numRecFields() );
     label startLoop = 0;
     label endLoop = 0;
     SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );
-    
+
     for (label i = 0; i < numRecFields; i++)
     {
-       if (i < numRecFields/2)
+        if (i < numRecFields/2)
         {
-        	startLoop = numRecFields/2;
-        	endLoop = numRecFields-1;
+            startLoop = numRecFields/2;
+            endLoop = numRecFields-1;
         }
         else
         {
-        	startLoop = 0;
-        	endLoop = numRecFields/2-1;
+            startLoop = 0;
+            endLoop = numRecFields/2-1;
         }
-        
+
         scalar nextMinimum(GREAT);
-	label jumpdest = 0;
+        label jumpdest = 0;
+        
         for (label j = startLoop; j < endLoop; j++)
         {
-        	if (recurrenceMatrix[j][i] < nextMinimum)
-        	{
-        		nextMinimum = recurrenceMatrix[j][i];
-        		jumpdest = j+1;
-        		continue;
-        	}
+            if (recurrenceMatrix[j][i] < nextMinimum)
+            {
+                nextMinimum = recurrenceMatrix[j][i];
+                jumpdest = j+1;
+                
+                continue;
+            }
         }
 
         jumpvec << jumpdest << endl;
     }
-    
+
     Info<< "\nComputing recurrence jump vector done\n" << endl;
 }
 
diff --git a/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.H b/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.H
index 5f42b69e..0c04aa66 100644
--- a/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.H
+++ b/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.H
@@ -20,6 +20,11 @@ License
 
     You should have received a copy of the GNU General Public License
     along with CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Create a path based on the recurrence statistics and add some randomness.
+    The minimum data base size is 5 entries for this model to work.
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef simpleRandomPath_H
diff --git a/src/recurrence/recStatAnalysis/noRecStatAnalysis/noRecStatAnalysis.C b/src/recurrence/recStatAnalysis/noRecStatAnalysis/noRecStatAnalysis.C
index 5c647177..1dde9e0f 100644
--- a/src/recurrence/recStatAnalysis/noRecStatAnalysis/noRecStatAnalysis.C
+++ b/src/recurrence/recStatAnalysis/noRecStatAnalysis/noRecStatAnalysis.C
@@ -27,8 +27,6 @@ License
 #include "noRecStatAnalysis.H"
 #include "recModel.H"
 #include "addToRunTimeSelectionTable.H"
-#include <mpi.h>
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
diff --git a/src/recurrence/recStatAnalysis/standardRecStatAnalysis/standardRecStatAnalysis.C b/src/recurrence/recStatAnalysis/standardRecStatAnalysis/standardRecStatAnalysis.C
index 3ae8b6bc..36ac44c9 100644
--- a/src/recurrence/recStatAnalysis/standardRecStatAnalysis/standardRecStatAnalysis.C
+++ b/src/recurrence/recStatAnalysis/standardRecStatAnalysis/standardRecStatAnalysis.C
@@ -27,7 +27,6 @@ License
 #include "standardRecStatAnalysis.H"
 #include "recModel.H"
 #include "addToRunTimeSelectionTable.H"
-#include <mpi.h>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //