diff --git a/CMakeLists.txt b/CMakeLists.txt
index 359ebc7..8c9586d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,7 +20,7 @@ endif()
 
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/")
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/findFFTW")
-set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/findFFTMPI")
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/findHEFFTE")
 set(LIBS)
 
 set(CMAKE_CXX_STANDARD 20) # C++20...
@@ -38,7 +38,7 @@ include_directories(${Boost_INCLUDE_DIRS})
 
 option(USE_MPI "Use MPI" ON)
 option(USE_OMP "Use OpenMP" OFF)
-option(USE_FFTMPI "Use FFTMPI" ON)
+option(USE_HEFFTE "Use Heffte" ON)
 
 ###mpi
 if(USE_MPI)
@@ -64,13 +64,13 @@ if(USE_OMP)
   message(SEND_ERROR "OMP not yet supported / experimental.")
 endif()
 
-if(USE_FFTMPI)
+if(USE_HEFFTE)
+  add_definitions(-DHAVE_HEFFTE)
+  find_package(Heffte REQUIRED)
   find_package(FFTW REQUIRED)
-  find_package(FFTMPI REQUIRED)
-  add_definitions(-DHAVE_FFTMPI)
+  include_directories(${HEFFTE_INCLUDES})
   include_directories(${FFTW_INCLUDE_DIRS})
-  include_directories(${FFTMPI_INCLUDES})
-  set(LIBS ${LIBS} "${FFTMPI_LIBRARIES}" "${FFTW_LIBRARIES}" )
+  set(LIBS ${LIBS} "${HEFFTE_LIBRARIES}" "${FFTW_LIBRARIES}")
 endif()
 
 option(PROFILE "Use gperftools" OFF)
@@ -81,6 +81,12 @@ if(PROFILE)
   set(LIBS ${LIBS} "${GPERFTOOLS_LIBRARIES}")
 endif()
 
+find_package(HDF5 COMPONENTS C CXX HL REQUIRED)
+include_directories( ${HDF5_INCLUDES} )
+include_directories (${HDF5_INCLUDE_DIR})
+set (LIBS ${LIBS} "${HDF5_CXX_LIBRARIES}")
+
+
 include_directories(src)
 
 add_subdirectory (src)
diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 3fa2ced..663d8bb 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -1,15 +1,24 @@
 add_executable (TG ./TG/TG.cpp)
-target_link_libraries (TG  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) 
+target_link_libraries (TG  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${HEFFTE_LIBRARIES} ${LIBS}) 
 install(TARGETS TG DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
 
 add_executable (TG_Poisson ./TG_Piosson/TG.cpp)
-target_link_libraries (TG_Poisson  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS}) 
+target_link_libraries (TG_Poisson  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${LIBS}) 
 install(TARGETS TG_Poisson DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
 
 add_executable (OT ./OT/OT.cpp)
-target_link_libraries (OT  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS})
+target_link_libraries (OT  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${HEFFTE_LIBRARIES} ${LIBS})
 install(TARGETS OT DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
 
 add_executable (OT_Poisson ./OT_Poisson/OT.cpp)
-target_link_libraries (OT_Poisson  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${FFTMPI_LIBRARIES} ${LIBS})
+target_link_libraries (OT_Poisson  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${HEFFTE_LIBRARIES} ${LIBS})
 install(TARGETS OT DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
+
+add_executable (shearLayer ./shearLayer/shearLayer.cpp)
+target_link_libraries (shearLayer  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${LIBS}) 
+install(TARGETS shearLayer DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
+
+add_executable (mhdjet ./mhdjet/mhdjet.cpp)
+target_link_libraries (mhdjet  NOVA ${Boost_LIBRARIES} ${MPI_LIBRARIES} ${LIBS}) 
+install(TARGETS mhdjet DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
+
diff --git a/app/OT/BEqn.H b/app/OT/BEqn.H
index 9f6224d..f988aae 100644
--- a/app/OT/BEqn.H
+++ b/app/OT/BEqn.H
@@ -29,12 +29,12 @@
 
     }
 
-    Bstar.correctBoundaryConditions();
-    pBEqn.rhs(  ex::div(Bstar) );  // /time.dt()
-    pBEqn.solve();
-    pB.correctBoundaryConditions();
-    B = Bstar - ex::grad(pB); //*time.dt()
-    B.correctBoundaryConditions();
+//    Bstar.correctBoundaryConditions();
+//    pBEqn.rhs(  ex::div(Bstar) );  // /time.dt()
+//    pBEqn.solve();
+//    pB.correctBoundaryConditions();
+//    B = Bstar - ex::grad(pB); //*time.dt()
+//    B.correctBoundaryConditions();
 
     if( time.writelog() )
     {   
diff --git a/app/OT_Poisson/OT.cpp b/app/OT_Poisson/OT.cpp
index 1fe3e67..3bd3cd7 100644
--- a/app/OT_Poisson/OT.cpp
+++ b/app/OT_Poisson/OT.cpp
@@ -23,7 +23,7 @@
 int main(int argc, char* argv[])
 {
     settings::process( argc, argv ); 
-    Time time( 0.001, 8, 50 ); //args: dt, endT, write interval / steps
+    Time time( 0.001, 2.01, 10 ); //args: dt, endT, write interval / steps
 
     const scalar pi = tools::pi;
     parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); 
@@ -45,12 +45,35 @@ int main(int argc, char* argv[])
     std::shared_ptr<Field<scalar> > pB_ptr( std::make_shared<Field<scalar> >( mesh, 0, "pB" ) );
     auto& pB = (*pB_ptr);
 
-    scalar mu = 0.01;
-    scalar eta = 0.01;
+    scalar mu = 0.057;
+    scalar eta = 0.057;
 
     poisson pEqn(p_ptr);
     poisson pBEqn(pB_ptr);
 
+    std::default_random_engine eng( parallelCom::myProcNo() );
+    std::uniform_real_distribution<double> dist(-1.0,1.0);
+
+/*
+    //initial conditions
+    for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2; i++ )
+    {   
+        for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2; j++ )
+        {   
+            for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2; k++ )
+            {   
+                scalar x = -pi+(i-settings::m()/2)*mesh.dx()+mesh.origin().x();
+                scalar y = -pi+(j-settings::m()/2)*mesh.dy()+mesh.origin().y();
+                scalar z = -pi+(k-settings::m()/2)*mesh.dz()+mesh.origin().z();
+
+                U(i, j, k) = vector(std::sin(x)*std::cos(y)*std::cos(z),-std::cos(x)*std::sin(y)*std::cos(z),0);
+		B(i, j, k) = vector(0.2, 0, 0);
+
+            }   
+        }   
+    }
+*/
+
     //initial conditions
     for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2; i++ )
     {
@@ -58,13 +81,15 @@ int main(int argc, char* argv[])
         {
             for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2; k++ )
             {
+		//OT Vortex initial conditions
                 scalar x = (i-settings::m()/2)*mesh.dx()+mesh.origin().x();
                 scalar y = (j-settings::m()/2)*mesh.dy()+mesh.origin().y();
                 scalar z = (k-settings::m()/2)*mesh.dz()+mesh.origin().z();
-		
+
 		U(i, j, k) = vector((-2.0)*std::sin(y),2.0*std::sin(x),0.0);
-                B(i, j, k) = 0.818*vector( ((-2.0)*std::sin(y*2.0))+std::sin(z),((2.0*std::sin(x))+std::sin(z)),( std::sin(x)+std::sin(y) ));
-            }
+		B(i, j, k) = 0.8*vector( ((-2.0)*std::sin(2.0*y))+std::sin(z),((2.0*std::sin(x))+std::sin(z)),( std::sin(x)+std::sin(y) ));
+
+	    }
         }
     }
 
@@ -75,6 +100,7 @@ int main(int argc, char* argv[])
 
     std::ofstream data( settings::zoneName()+"/"+"dat.dat");
     scalar EkOld = 0.0;
+    scalar EmOld = 0.0;
 
     boost::timer::cpu_timer timer;
 
@@ -84,7 +110,7 @@ int main(int argc, char* argv[])
 
         time++;
 
-        #include "RKLoop.H" 
+        #include "UEqn.H" 
 	#include "BEqn.H"	
 
         mesh.write(settings::zoneName()+"/data");
@@ -94,8 +120,8 @@ int main(int argc, char* argv[])
 	J.write(settings::zoneName()+"/data", "J");
 	pB.write(settings::zoneName()+"/data", "pB");
 
-        if( time.writelog() )
-        {
+//        if( time.writelog() )
+//        {
             if( parallelCom::master() )
             { 
                 std::cout << "Step: " << time.timeStep() << ".              Time: " << time.curTime() << std::endl;
@@ -103,11 +129,12 @@ int main(int argc, char* argv[])
             }
 
             tools::CFL( U, mesh );
-        }
+//        }
 
         scalar Ek=0.0;
         scalar Em=0.0;
         scalar Jmax=0.0;
+	scalar Umax=0.0;
 	omegav = ex::curl(U);
         scalar Jav=0.0;
         scalar omega=0.0;
@@ -122,6 +149,7 @@ int main(int argc, char* argv[])
                     Ek += 0.5 * (U(i, j, k).x() * U(i, j, k).x() + U(i, j, k).y() * U(i, j, k).y() + U(i, j, k).z() * U(i, j, k).z() );
                     Em += 0.5 * (B(i, j, k).x() * B(i, j, k).x() + B(i, j, k).y() * B(i, j, k).y() + B(i, j, k).z() * B(i, j, k).z() );
 		    Jmax = std::max( Jmax, sqrt(J(i, j, k).x() * J(i, j, k).x() + J(i, j, k).y() * J(i, j, k).y() + J(i, j, k).z() * J(i, j, k).z())); 
+		    Umax = std::max( Umax, sqrt(U(i, j, k).x() * U(i, j, k).x() + U(i, j, k).y() * U(i, j, k).y() + U(i, j, k).z() * U(i, j, k).z()));
 		    omega += sqrt(omegav(i, j, k).x() * omegav(i, j, k).x() + omegav(i, j, k).y() * omegav(i, j, k).y() + omegav(i, j, k).z() * omegav(i, j, k).z());
                     Jav += sqrt(J(i, j, k).x() * J(i, j, k).x() + J(i, j, k).y() * J(i, j, k).y() + J(i, j, k).z() * J(i, j, k).z());
 		    n++;
@@ -131,6 +159,7 @@ int main(int argc, char* argv[])
     
         reduce( Ek, plusOp<scalar>() );
 	reduce( Jmax, maxOp<scalar>() );
+	reduce( Umax, maxOp<scalar>() );
         reduce( Em, plusOp<scalar>() );
         reduce( Jav, plusOp<scalar>() );
         reduce( omega, plusOp<scalar>() );
@@ -141,15 +170,18 @@ int main(int argc, char* argv[])
         Jav /= n;
         omega /= n;
 
+//	Jmax=Jmax*17.54;
+
 	scalar epsilon=0.0;
         epsilon = -mu*(omega*omega) - eta*(Jav*Jav);
 
         if( time.curTime() > time.dt() && parallelCom::master() )
         {
-            data<<time.curTime()<<" "<<std::setprecision(15)<<Ek<<" "<<std::setprecision(15)<<Em<<" "<<std::setprecision(15)<<Jav<<" "<<std::setprecision(15)<<Jmax<<" "<<std::setprecision(15)<<omega<<" "<<std::setprecision(15)<<epsilon<<" "<<std::setprecision(15)<<-(Ek-EkOld)/time.dt()<<std::endl;
+            data<<time.curTime()<<" "<<std::setprecision(15)<<Ek<<" "<<std::setprecision(15)<<Em<<" "<<std::setprecision(15)<<Jav<<" "<<std::setprecision(15)<<Jmax<<" "<<std::setprecision(15)<<omega<<" "<<std::setprecision(15)<<epsilon<<" "<<std::setprecision(15)<<-(Ek-EkOld)/time.dt()<<std::setprecision(15)<<" "<<(Em-EmOld)/time.dt()<<std::setprecision(15)<<" "<<Umax<<std::endl;
 	}
 
         EkOld = Ek;
+	EmOld = Em;
  
     }
 
diff --git a/app/OT_Poisson/RKLoop.H b/app/OT_Poisson/UEqn.H
similarity index 88%
rename from app/OT_Poisson/RKLoop.H
rename to app/OT_Poisson/UEqn.H
index 6a64632..06cd3c8 100644
--- a/app/OT_Poisson/RKLoop.H
+++ b/app/OT_Poisson/UEqn.H
@@ -28,14 +28,15 @@
             Ustar = Ustar+settings::RK4a()[rk][j]*dU[j];
         }
 
+    }
+
         Ustar.correctBoundaryConditions();
-        pEqn.rhs(ex::div(Ustar) / time.dt() / settings::RK4c()[rk]);
+        pEqn.rhs(ex::div(Ustar) / time.dt());
         pEqn.solve();
         p.correctBoundaryConditions();
-        U = Ustar - settings::RK4c()[rk]*time.dt() * ex::grad(p);
+        U = Ustar - time.dt() * ex::grad(p);
 
         U.correctBoundaryConditions();
-    }
 
 
     if( time.writelog() )
diff --git a/app/TG/TG.cpp b/app/TG/TG.cpp
index c9b225e..eaa0bfb 100644
--- a/app/TG/TG.cpp
+++ b/app/TG/TG.cpp
@@ -21,7 +21,7 @@
 int main(int argc, char* argv[])
 {
     settings::process( argc, argv ); 
-    Time time( 0.001, 20, 1000 ); //args: dt, endT, write interval / steps
+    Time time( 0.01, 6.01, 100 ); //args: dt, endT, write interval / steps
 
     const scalar pi = 3.1415926536;
     parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); 
diff --git a/app/TG_Piosson/RKLoop.H b/app/TG_Piosson/RKLoop.H
index 3ae261b..38e8fe2 100644
--- a/app/TG_Piosson/RKLoop.H
+++ b/app/TG_Piosson/RKLoop.H
@@ -1,4 +1,5 @@
 {
+
     std::array<std::shared_ptr<Field<vector> >, 4> dU;
     for( int i=0; i<4; i++ )
     {
@@ -25,16 +26,18 @@
         {
             Ustar = Ustar+settings::RK4a()[rk][j]*dU[j];
         }
-
+  }
         Ustar.correctBoundaryConditions();
-        pEqn.rhs(ex::div(Ustar) / time.dt() / settings::RK4c()[rk]);
+
+        pEqn.rhs(ex::div(Ustar) / time.dt() );
+
 	pEqn.solve();
         p.correctBoundaryConditions();
-        U = Ustar - settings::RK4c()[rk]*time.dt() * ex::grad(p);
+        U = Ustar - time.dt() * ex::grad(p);
 
         U.correctBoundaryConditions();
-    }
-
+  
+  
 
     if( time.writelog() )
     {
diff --git a/app/TG_Piosson/TG.cpp b/app/TG_Piosson/TG.cpp
index 9ff8dc2..8c85fb6 100644
--- a/app/TG_Piosson/TG.cpp
+++ b/app/TG_Piosson/TG.cpp
@@ -22,7 +22,7 @@
 int main(int argc, char* argv[])
 {
     settings::process( argc, argv ); 
-    Time time( 0.01, 20, 1.0 ); //args: dt, endT, write interval / steps
+    Time time( 0.01, 20, 1 ); //args: dt, endT, write interval / steps
 
     const scalar pi = tools::pi;
     parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); 
@@ -53,6 +53,7 @@ int main(int argc, char* argv[])
                 scalar z = -pi+(k-settings::m()/2)*mesh.dz()+mesh.origin().z();
 
                 U(i, j, k) = vector(std::sin(x)*std::cos(y)*std::cos(z),-std::cos(x)*std::sin(y)*std::cos(z),0);
+//		U(i, j, k) = vector(0.0, 0.0, 0.0);
 	    }
         }
     }
@@ -79,8 +80,8 @@ int main(int argc, char* argv[])
         p.write(settings::zoneName()+"/data", "p"); 
 
 
-        if( time.writelog() )
-        {
+//        if( time.writelog() )
+//        {
             if( parallelCom::master() )
             { 
                 std::cout << "Step: " << time.timeStep() << ".              Time: " << time.curTime() << std::endl;
@@ -88,7 +89,7 @@ int main(int argc, char* argv[])
             }
 
             tools::CFL( U, mesh );
-        }
+//        }
 
         scalar Ek=0.0;
         int n=0;
diff --git a/app/mhdjet/BEqn.H b/app/mhdjet/BEqn.H
new file mode 100644
index 0000000..f7861ab
--- /dev/null
+++ b/app/mhdjet/BEqn.H
@@ -0,0 +1,56 @@
+{
+    std::array<std::shared_ptr<Field<vector> >, 4> dB; 
+    for( int i=0; i<4; i++ )
+    {   
+        reuseTmp<vector> dBt( mesh );
+        dB[i]= dBt();
+    }
+
+
+    reuseTmp<vector> Boldt( mesh );
+    std::shared_ptr<Field<vector> > Bold( Boldt() );
+
+    *Bold = B;
+
+    for( int rk=0; rk<4; rk++ )
+    {
+
+	dB[rk] =
+        time.dt() *
+        (
+           -ex::div( U, B )
+	   +ex::div( B, U ) 
+           +ex::laplacian( eta, B ) 
+        );
+	
+	Bstar = Bold;
+	for( int j=0; j<=rk; j++ )
+        {   
+	      Bstar = Bstar + settings::RK4a()[rk][j]*dB[j]; 
+	} 
+
+	B = Bstar;
+	B.correctBoundaryConditions();
+
+    }
+
+    pBEqn.rhs(ex::div(B));
+    pBEqn.solve();
+    pB.correctBoundaryConditions();
+    B -= ex::grad(pB);
+    B.correctBoundaryConditions();
+
+    J = ex::curl(B);
+    J.correctBoundaryConditions();
+
+    if( time.writelog() )
+    {   
+       auto divB = ex::div(B);
+       auto maxdivB = tools::maximum(divB);
+        if( parallelCom::master() )
+        {   
+            std::cout<<"Magnetic field divergence: " << maxdivB<<std::endl;
+        }   
+    }
+
+}
diff --git a/app/mhdjet/UEqn.H b/app/mhdjet/UEqn.H
new file mode 100644
index 0000000..745a66c
--- /dev/null
+++ b/app/mhdjet/UEqn.H
@@ -0,0 +1,59 @@
+{
+    std::array<std::shared_ptr<Field<vector> >, 4> dU;
+    for( int i=0; i<4; i++ )
+    {
+        reuseTmp<vector> dUt( mesh );
+        dU[i]= dUt();
+    }
+
+    reuseTmp<vector> UOldt( mesh );
+    std::shared_ptr<Field<vector> > UOld( UOldt() );
+
+    reuseTmp<vector> Ustartt( mesh );
+    std::shared_ptr<Field<vector> > Ustart( Ustartt() ); 
+    *Ustart = Ustart;
+
+    *UOld = U;
+
+    for( int rk=0; rk<4; rk++ )
+    {
+        dU[rk] =
+        time.dt() *
+        (
+           -ex::div( U, U )
+           +ex::laplacian( mu, U ) 
+           +ex::div( B, B )
+           -ex::grad( 0.5*(B&B) )
+	);
+
+        Ustar = UOld;
+        for( int j=0; j<=rk; j++ )
+        {
+            Ustar = Ustar+settings::RK4a()[rk][j]*dU[j];
+        }
+
+    }
+
+        Ustar.correctBoundaryConditions();
+        pEqn.rhs(ex::div(Ustar) / time.dt());
+        pEqn.solve();
+        p.correctBoundaryConditions();
+        U = Ustar - time.dt() * ex::grad(p);
+
+        U.correctBoundaryConditions();
+
+    Umag = (U&U);
+
+    if( time.writelog() )
+    {
+        auto divU = ex::div(U);
+        auto maxdivU = tools::maximum(divU);
+
+        if( parallelCom::master() )
+        {
+            std::cout<<"Velocity divergence: " << maxdivU<<std::endl;
+        }
+    }
+
+
+}
diff --git a/app/mhdjet/mhdjet.cpp b/app/mhdjet/mhdjet.cpp
new file mode 100644
index 0000000..b9949a2
--- /dev/null
+++ b/app/mhdjet/mhdjet.cpp
@@ -0,0 +1,470 @@
+#include "Field/Field.h"
+#include <iostream>
+#include "Types/vector/vector.h"
+#include "ex/grad/grad.h"
+#include "ex/div/div.h"
+#include "ex/curl/curl.h"
+#include "ex/laplacian/laplacian.h"
+#include "Time/Time.h"
+#include "parallelCom/parallelCom.h"
+#include "settings/settings.h"
+#include "BC/BC.h"
+#include <functional>
+#include "Tools/Tools.h"
+#include "poisson/poisson.h"
+#include <boost/timer/timer.hpp>
+
+#include "H5cpp.h"
+
+#include <cmath>
+#include <memory>
+#include <random>
+
+    //H5 Create
+    //H5 Write dataset using MPI
+    void writeVectorFieldToHDF5( MPI_Comm comm, const std::string &fileName, const std::vector<std::vector<std::vector<double>>> &vectorField) {
+    
+    MPI_Comm world = parallelCom::world();
+    int rank, size;
+    MPI_Comm_rank(comm, &rank);
+    MPI_Comm_size(comm, &size);
+
+    // Dimensions of the vector field
+    hsize_t hni = vectorField.size();
+    hsize_t hnj = vectorField[0].size();
+    hsize_t hnk = vectorField[0][0].size();
+
+    //prepare data
+    std::vector<double> data(hni * hnj * hnk );  
+    int n=0;
+    for (int i = 0; i < hni; ++i) {
+        for (int j = 0; j < hnj; ++j) {
+            for (int k = 0; k < hnk; ++k) {
+            data[n] = vectorField[i][j][k];
+            n++;
+            }
+        }
+    }
+
+    // Create an HDF5 file
+    hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
+    H5Pset_fapl_mpio(plist_id, parallelCom::world(), MPI_INFO_NULL);
+    hid_t file_id = H5Fcreate(fileName.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
+
+    // Create a data space
+    hsize_t dims[] = {   static_cast<hsize_t>((hni)*parallelCom::ni() ),
+			 static_cast<hsize_t>((hnj)*parallelCom::nj() ),
+			 static_cast<hsize_t>((hnk)*parallelCom::nk() )  };  // 3 for the vector components (x, y, z)
+    hid_t dataspace_id = H5Screate_simple(3, dims, NULL);
+
+    //Create a dataset
+    hid_t dataset_id = H5Dcreate2(file_id, "/3DArray", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+    // Create the local dataspace for each MPI rank
+    hsize_t local_dims[3] = {hni, hnj, hnk};
+    hsize_t start[3] = {static_cast<hsize_t>((hni*parallelCom::i()) ) , static_cast<hsize_t>((hnj*parallelCom::j()) ) , static_cast<hsize_t>((hnk*parallelCom::k()) ) };
+    hsize_t count[3] = {hni, hnj, hnk};
+    hid_t memspace_id = H5Screate_simple(3, local_dims, NULL);
+
+
+    //Sanity Checks for data indexing
+//    std::cout << "  ni:  " << parallelCom::ni() << "  nj:  " << parallelCom::nj() << "  nk:  " << parallelCom::nk() << "  i:  " <<parallelCom::i()  << "  j:  " <<parallelCom::j()<< "  k:  " <<parallelCom::k() << "  hni  " << hni << "  hnj  " << hnj << "  hnk  " << hnk << "  dims_x  " << dims[0] << "  dims_y  " << dims[1] << "  dims_z  " << dims[2] << "  LocalStart_x  " << start[0] << "  LocalStart_y  " << start[1]<< "  LocalStart_z  " << start[2] << std::endl;
+//    std::cout << " myProNo_x " << parallelCom::myProcNo() << std::endl;
+
+
+    // Select the hyperslab to write
+    H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, start, NULL, count, NULL);
+
+    // Write the data from each processor into the HDF5 dataset
+    H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, memspace_id, dataspace_id, H5P_DEFAULT, data.data());
+
+    // Close the dataset and file
+    // Close resources
+    H5Pclose(plist_id);
+    H5Sclose(dataspace_id);
+    H5Sclose(memspace_id);
+    H5Dclose(dataset_id);
+    H5Fclose(file_id);
+
+}
+
+
+/*
+//H5 Write dataset in serial
+    void writeVectorFieldToHDF5(const std::string &fileName, const Field<vector> &vectorField) {
+    // Dimensions of the vector field
+    hsize_t hni = vectorField[0]-settings::m();
+    hsize_t hnj = vectorField[1]-settings::m();
+    hsize_t hnk = vectorField[2]-settings::m();
+
+    // Create an HDF5 file
+    H5::H5File file(fileName, H5F_ACC_TRUNC);
+
+    // Create a data space
+    hsize_t dims[3] = {hni, hnj, hnk };  // 3 for the vector components (x, y, z)
+    H5::DataSpace dataspace(3, dims);
+
+   // Create a dataset
+    H5::DataSet dataset = file.createDataSet("VectorField", H5::PredType::IEEE_F64LE, dataspace);
+
+    // Prepare data to write
+    std::vector<double> data(hni * hnj * hnk * 3); 
+    for (int i = 0; i < hni; ++i) {
+        for (int j = 0; j < hnj; ++j) {
+            for (int k = 0; k < hnk; ++k) {
+                data[(i * hnj * hnk + j * hnk + k)  + 0] = vectorField[0];  // x component
+                data[(i * hnj * hnk + j * hnk + k)  + 1] = vectorField[1];  // y component
+                data[(i * hnj * hnk + j * hnk + k)  + 2] = vectorField[2];  // z component
+                
+            }
+        }
+    }
+
+    // Write data to the dataset
+    dataset.write(data.data(), H5::PredType::IEEE_F64LE);
+
+    // Close the dataset and file
+    dataset.close();
+    file.close();
+}
+*/
+
+
+
+int main(int argc, char* argv[])
+{
+    settings::restart() = true;
+
+    settings::process( argc, argv ); 
+    Time time( 0.02, 160.02, 50 ); //args: dt, endT, write interval / steps
+
+    const scalar pi = tools::pi;
+    parallelCom::decompose( settings::zoneName()+"/"+"mesh" ); 
+
+    Mesh mesh( settings::zoneName()+"/"+"mesh", time );
+    
+    mesh.write(settings::zoneName()+"/data");
+
+
+    Field<vector> U( mesh, vector(0,0,0), "U" );
+    Field<vector> Ustar( mesh, vector(0,0,0), "U" );
+    Field<vector> Ustart( mesh, vector(0,0,0), "U" );
+    Field<vector> B( mesh, vector(0,0,0), "B" );
+    Field<vector> Bstar( mesh, vector(0,0,0), "B" );
+    Field<vector> J( mesh, vector(0,0,0), "B" );
+    Field<vector> omegav( mesh, vector(0,0,0), "U" );
+
+    std::shared_ptr<Field<scalar> > Umag_ptr(std::make_shared<Field<scalar> >(mesh, 0, "U") );
+    auto& Umag = (*Umag_ptr);
+
+    std::shared_ptr<Field<scalar> > p_ptr( std::make_shared<Field<scalar> >( mesh, 0, "p" ) );
+    auto& p = (*p_ptr);
+
+    std::shared_ptr<Field<scalar> > pB_ptr( std::make_shared<Field<scalar> >( mesh, 0, "pB" ) );
+    auto& pB = (*pB_ptr);
+
+    std::shared_ptr<Field<scalar> > fftx_ptr( std::make_shared<Field<scalar> >( mesh, 0, "U" ) );
+    auto& fftx = (*fftx_ptr);
+
+    scalar mu = 0.000285714;
+    scalar eta = 0.000285714;
+
+    poisson fftxDo(fftx_ptr);
+    poisson pEqn(p_ptr);
+    poisson pBEqn(pB_ptr);
+
+    std::default_random_engine eng( parallelCom::myProcNo() );
+    std::uniform_real_distribution<double> dist(-1.0,1.0);
+
+    //initial conditions
+    for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2; i++ )
+    {
+        for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2; j++ )
+        {
+            for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2; k++ )
+            {
+                scalar x = (i-settings::m()/2)*mesh.dx()+mesh.origin().x();
+                scalar y = -10+(j-settings::m()/2)*mesh.dy()+mesh.origin().y();
+                scalar z = (k-settings::m()/2)*mesh.dz()+mesh.origin().z();
+		
+		int perU=0.0;
+                for( int px=0; px<2; px++ )
+                {   
+                    perU += dist(eng)*sin((px*x)-(dist(eng)*13.96));
+                }
+
+		for( int px=3; px<mesh.ni()/3; px++ )
+		{
+		    perU += dist(eng)*sin((px*x)-(dist(eng)*13.96));
+		}
+
+                int perZ=0.0;
+                for( int pz=0; pz<2; pz++ )
+                {   
+                    perZ += dist(eng)*sin((pz*z)-(dist(eng)*13.96));
+                }   
+
+                for( int pz=3; pz<mesh.nk()/3; pz++ )
+                {   
+                    perZ += dist(eng)*sin((pz*z)-(dist(eng)*13.96));
+                } 
+
+		U(i, j, k) = vector(1.0/(std::cosh(y)*std::cosh(y)), 0.0, 0.0);
+		U(i, j, k).y() += (0.0001*dist(eng))*std::exp(-(y*y));
+		U(i, j, k).x() += (0.001*std::sin(x*0.9) + 0.0001*dist(eng))*std::exp(-(y*y));
+		U(i, j, k).z() += (0.0001*dist(eng))*std::exp(-(y*y));		
+
+		B(i, j, k) = vector(0.015, 0.0, 0.0);
+
+	    }
+        }
+    }
+
+    U.correctBoundaryConditions();
+    p.correctBoundaryConditions();
+    B.correctBoundaryConditions();
+    pB.correctBoundaryConditions();
+
+    std::ofstream data( settings::zoneName()+"/"+"dat.dat");
+    scalar EkOld = 0.0;
+    scalar EmOld = 0.0;
+
+    boost::timer::cpu_timer timer;
+
+
+    while( time.run() )
+    {
+        double start = timer.elapsed().wall / 1e9;
+
+
+        mesh.write(settings::zoneName()+"/data");
+        U.write(settings::zoneName()+"/data", "U");
+        p.write(settings::zoneName()+"/data", "p"); 
+	B.write(settings::zoneName()+"/data", "B");
+	J.write(settings::zoneName()+"/data", "J");
+	pB.write(settings::zoneName()+"/data", "pB");
+
+        time++;
+        #include "UEqn.H" 
+        #include "BEqn.H"
+
+        if( time.writelog() )
+        {
+            if( parallelCom::master() )
+            { 
+                std::cout << "Step: " << time.timeStep() << ".              Time: " << time.curTime() << std::endl;
+                std::cout << "Elapsed wall time for timestep: " <<  timer.elapsed().wall / 1e9 - start <<std::endl;
+            }
+
+            tools::CFL( U, mesh );
+        }
+
+
+//	HDF5 data writing
+	if( time.write() )
+        {   
+        int hni = mesh.glob_n()[0]/parallelCom::ni();
+        int hnj = mesh.glob_n()[1]/parallelCom::nj();
+        int hnk = mesh.glob_n()[2]/parallelCom::nk();
+
+	std::vector<std::vector<std::vector<double>>> vectorField(hni, std::vector<std::vector<double>>(hnj, std::vector<double>(hnk)));
+	//write Ux
+	int n=0;
+	for (int k = 0; k < hnk; ++k) {
+        	for (int j = 0; j < hnj; ++j) {
+            	    for (int i = 0; i < hni; ++i) {
+                	vectorField[i][j][k] = U( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).x();  // x component
+            		n++;
+			}
+        	}
+    	}
+	std::string H5title = "./h5data/Ux" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+	writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+	//write Uy
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = U( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).y();  // x component
+                        n++;
+                        }   
+                }   
+        }	
+        H5title = "./h5data/Uy_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+        //write Uz
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = U( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).z();  // x component
+                        n++;
+                        }   
+                }   
+        }           
+        H5title = "./h5data/Uz_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+        //write p
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = p( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) );  // x component
+                        n++;
+                        }   
+                }   
+        }           
+        H5title = "./h5data/p_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+        //write Bx
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = B( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).x();  // x component
+                        n++;
+                        }   
+                }   
+        }    
+        H5title = "./h5data/Bx_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+        //write By
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = B( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).y();  // x component
+                        n++;
+                        }
+                }
+        }
+        H5title = "./h5data/By_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+        //write Bz
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = B( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).z();  // x component
+                        n++;
+                        }
+                }
+        }
+        H5title = "./h5data/Bz_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+	//write Jx
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = J( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).x();  // x component
+                        n++;
+                        }   
+                }   
+        }   
+        H5title = "./h5data/Jx_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+        //write Jy
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = J( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).y();  // x component
+                        n++;
+                        }
+                }
+        }
+        H5title = "./h5data/Jy_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+        //write Jz
+        n=0;
+        for (int k = 0; k < hnk; ++k) {
+                for (int j = 0; j < hnj; ++j) {
+                    for (int i = 0; i < hni; ++i) {
+                        vectorField[i][j][k] = J( (i+settings::m()/2), (j+settings::m()/2), (k+settings::m()/2) ).z();  // x component
+                        n++;
+                        }
+                }
+        }
+        H5title = "./h5data/Jz_" +  std::to_string( (time.curTime()-time.dt()) ) + ".h5";
+        writeVectorFieldToHDF5(parallelCom::world(), H5title, vectorField);
+
+
+
+
+
+	if( parallelCom::master() )
+	  {
+          std::cout << "Vector field written to vector_field.h5" << std::endl;
+	  }
+	}
+
+
+        scalar Ek=0.0;
+        scalar Em=0.0;
+        scalar Jmax=0.0;
+	scalar Umax=0.0;
+	omegav = ex::curl(U);
+        scalar Jav=0.0;
+        scalar omega=0.0;
+        int n=0;
+    
+        for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2-1; i++ )
+        {
+            for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2-1; j++ )
+            {
+                for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2-1; k++ )
+                {
+                    Ek += 0.5 * (U(i, j, k).x() * U(i, j, k).x() + U(i, j, k).y() * U(i, j, k).y() + U(i, j, k).z() * U(i, j, k).z() );
+                    Em += 0.5 * (B(i, j, k).x() * B(i, j, k).x() + B(i, j, k).y() * B(i, j, k).y() + B(i, j, k).z() * B(i, j, k).z() );
+		    Jmax = std::max( Jmax, sqrt(J(i, j, k).x() * J(i, j, k).x() + J(i, j, k).y() * J(i, j, k).y() + J(i, j, k).z() * J(i, j, k).z())); 
+		    Umax = std::max( Umax, sqrt(U(i, j, k).x() * U(i, j, k).x() + U(i, j, k).y() * U(i, j, k).y() + U(i, j, k).z() * U(i, j, k).z()));
+		    omega += sqrt(omegav(i, j, k).x() * omegav(i, j, k).x() + omegav(i, j, k).y() * omegav(i, j, k).y() + omegav(i, j, k).z() * omegav(i, j, k).z());
+                    Jav += sqrt(J(i, j, k).x() * J(i, j, k).x() + J(i, j, k).y() * J(i, j, k).y() + J(i, j, k).z() * J(i, j, k).z());
+		    n++;
+                }
+            }
+        }
+    
+        reduce( Ek, plusOp<scalar>() );
+	reduce( Jmax, maxOp<scalar>() );
+	reduce( Umax, maxOp<scalar>() );
+        reduce( Em, plusOp<scalar>() );
+        reduce( Jav, plusOp<scalar>() );
+        reduce( omega, plusOp<scalar>() );
+	reduce( n, plusOp<int>() );
+
+        Ek /= n;
+	Em /= n;
+        Jav /= n;
+        omega /= n;
+
+	scalar epsilon=0.0;
+        epsilon = -mu*(omega*omega) - eta*(Jav*Jav);
+
+        if( time.curTime() > time.dt() && parallelCom::master() )
+        {
+            data<<time.curTime()<<" "<<std::setprecision(15)<<Ek<<" "<<std::setprecision(15)<<Em<<" "<<std::setprecision(15)<<Jav<<" "<<std::setprecision(15)<<Jmax<<" "<<std::setprecision(15)<<omega<<" "<<std::setprecision(15)<<epsilon<<" "<<std::setprecision(15)<<-(Ek-EkOld)/time.dt()<<std::setprecision(15)<<" "<<(Em-EmOld)/time.dt()<<std::setprecision(15)<<" "<<Umax<<std::endl;
+	}
+
+        EkOld = Ek;
+	EmOld = Em;
+ 
+    }
+
+
+    std::cout<< timer.elapsed().wall / 1e9 <<std::endl;
+    
+    mesh.write(settings::zoneName()+"/data");
+
+    return 0;
+}
diff --git a/app/shearLayer/BEqn.H b/app/shearLayer/BEqn.H
new file mode 100644
index 0000000..42dd404
--- /dev/null
+++ b/app/shearLayer/BEqn.H
@@ -0,0 +1,56 @@
+{
+    std::array<std::shared_ptr<Field<vector> >, 4> dB; 
+    for( int i=0; i<4; i++ )
+    {   
+        reuseTmp<vector> dBt( mesh );
+        dB[i]= dBt();
+    }
+
+
+    reuseTmp<vector> Boldt( mesh );
+    std::shared_ptr<Field<vector> > Bold( Boldt() );
+
+    *Bold = B;
+
+    for( int rk=0; rk<4; rk++ )
+    {
+
+	dB[rk] =
+        time.dt() *
+        (
+           -ex::div( U, B )
+	   +ex::div( B, U ) 
+           +ex::laplacian( eta, B ) 
+        );
+	
+	Bstar = Bold;
+	for( int j=0; j<=rk; j++ )
+        {   
+	      Bstar = Bstar + settings::RK4a()[rk][j]*dB[j]; 
+	} 
+
+	B = Bstar;
+	B.correctBoundaryConditions();
+
+    }
+
+    pBEqn.rhs(ex::div(B));
+    pBEqn.solve();
+    pB.correctBoundaryConditions();
+    B -= ex::grad(pB);
+    B.correctBoundaryConditions();
+
+//    J = ex::curl(B);
+//    J.correctBoundaryConditions();
+
+    if( time.writelog() )
+    {   
+       auto divB = ex::div(B);
+       auto maxdivB = tools::maximum(divB);
+        if( parallelCom::master() )
+        {   
+            std::cout<<"Magnetic field divergence: " << maxdivB<<std::endl;
+        }   
+    }
+
+}
diff --git a/app/shearLayer/TEqn.H b/app/shearLayer/TEqn.H
new file mode 100644
index 0000000..59b7933
--- /dev/null
+++ b/app/shearLayer/TEqn.H
@@ -0,0 +1,38 @@
+{
+    reuseTmp<scalar> dTt( mesh );
+    std::shared_ptr<Field<scalar> > dT( dTt() );
+  
+    reuseTmp<scalar> Tct( mesh );
+    std::shared_ptr<Field<scalar> > Tc( Tct() );
+
+    reuseTmp<scalar> TOldt( mesh );
+    std::shared_ptr<Field<scalar> > TOld( TOldt() );
+
+    *Tc = T;
+    *TOld = T;
+
+    for( int rk=0; rk<4; rk++ )
+    {
+        *dT =
+        time.dt() *
+        (
+           -ex::div( U, T )
+           +ex::laplacian( ka, T )
+        );
+
+//        drhoT->correctBoundaryConditions();
+
+        *Tc = Tc+RK4[rk]*dT;
+
+        if( rk<3 )
+        {
+            T = TOld + (int(rk/2)+1)/2.0 * (dT);
+        }
+        else
+        {
+            T=Tc;
+        }
+    
+        T.correctBoundaryConditions();
+    }
+}
diff --git a/app/shearLayer/UEqn.H b/app/shearLayer/UEqn.H
new file mode 100644
index 0000000..b19fa90
--- /dev/null
+++ b/app/shearLayer/UEqn.H
@@ -0,0 +1,52 @@
+{
+    reuseTmp<vector> dUt( mesh );
+    std::shared_ptr<Field<vector> > dU( dUt() );
+
+    reuseTmp<vector> UOldt( mesh );
+    std::shared_ptr<Field<vector> > UOld( UOldt() );
+
+    reuseTmp<vector> Uct( mesh );
+    std::shared_ptr<Field<vector> > Uc( Uct() );
+
+    *UOld = U;
+
+
+    for( int rk=0; rk<4; rk++ )
+    {
+        *dU =
+        time.dt() *
+        (
+           -ex::div( U, U )
+           +ex::laplacian( mu, U ) 
+//           +mu*ex::div( ( dev2(transpose(gradU) ) ) )
+//           -( (rho0*T)*(g) )  
+           +ex::div( B, B ) 
+           -ex::grad( 0.5*(B&B) )
+        );
+
+
+        Ustar = UOld;
+        Ustar = Ustar+RK4[rk]*dU;
+   
+	Ustar.correctBoundaryConditions();
+        pEqn.rhs(ex::div(Ustar) / time.dt() / settings::RK4c()[rk]);
+        pEqn.solve();
+        p.correctBoundaryConditions();
+        U = Ustar - settings::RK4c()[rk]*time.dt() * ex::grad(p); 
+
+        U.correctBoundaryConditions();
+   } 
+
+//    if( time.writelog() )
+//    {   
+        auto divU = ex::div(U);
+        auto maxdivU = tools::maximum(divU);
+
+        if( parallelCom::master() )
+        {   
+            std::cout<<"Velocity divergence: " << maxdivU<<std::endl;
+        }   
+//    }
+
+
+}
diff --git a/app/shearLayer/shearLayer.cpp b/app/shearLayer/shearLayer.cpp
new file mode 100644
index 0000000..4b34c15
--- /dev/null
+++ b/app/shearLayer/shearLayer.cpp
@@ -0,0 +1,225 @@
+#include "Field/Field.h"
+#include <iostream>
+#include "Types/vector/vector.h"
+#include "ex/grad/grad.h"
+#include "ex/div/div.h"
+#include "ex/laplacian/laplacian.h"
+#include "Time/Time.h"
+#include "parallelCom/parallelCom.h"
+#include "settings/settings.h"
+#include "BC/BC.h"
+#include <functional>
+#include "Tools/Tools.h"
+#include "poisson/poisson.h"
+
+#include <boost/timer/timer.hpp>
+
+#include <cmath>
+#include <memory>
+#include <random>
+
+int main(int argc, char* argv[])
+{ 
+    settings::process( argc, argv );  
+    Time time( 0.001, 420, 1 ); //args: dt, endT, write interval / steps
+
+    const scalar pi = tools::pi;
+    parallelCom::decompose( settings::zoneName()+"/"+"mesh" );  
+
+    Mesh mesh( settings::zoneName()+"/"+"mesh", time );
+
+    mesh.write(settings::zoneName()+"/data"); 
+/*--------------------
+    scalar Mach=0.2;
+    scalar Re0=1967.0;
+    scalar Ri0=0.08;
+    scalar Pr=7;
+
+    scalar rho0=999.57;
+    scalar c=0.0417;
+
+    scalar u0=c*Mach;
+    scalar h0=0.2;//Re0*mu/rho0/u0;
+    scalar t0=0.1;//-Ri0*u0*u0/g[2]/h0;
+    scalar mu = rho0*u0*h0/Re0;
+    scalar ka = mu/Pr;
+    scalar eta = mu;
+//    vector g( 0, 0, -Ri0*u0*u0/t0/h0 );
+    vector g( 0, 0, -0.000011129 );
+
+    Field<vector> U( mesh, vector(1,0,0), "U"  );  
+    Field<scalar> T( mesh, 0, "T" );
+    Field<scalar> pHyd( mesh,0, "pHyd" );
+    Field<vector> Ustar( mesh, vector(0,0,0), "U" );    
+//    Field<vector> B( mesh, vector(0,0,0), "B" );
+//    Field<vector> Bstar( mesh, vector(0,0,0), "B" );
+
+----------------------*/
+
+	Field<vector> U( mesh, vector(0,0,0), "U"  ); 
+	Field<vector> Ustar( mesh, vector(0,0,0), "U" );
+	Field<vector> B( mesh, vector(0,0,0), "B" );
+	Field<vector> Bstar( mesh, vector(0,0,0), "B" );
+
+
+	scalar mu = 1.0/28000.0;
+	scalar u0=1.0;
+	scalar delta = 1.0/28.0;
+	scalar cn = 0.001;
+	scalar eta = mu*10.0;	
+
+
+    std::default_random_engine eng( parallelCom::myProcNo() );
+    std::uniform_real_distribution<double> dist(-1.0,1.0);
+
+    std::shared_ptr<Field<scalar> > p_ptr( std::make_shared<Field<scalar> >( mesh, (0.87), "p" ) );
+    auto& p = (*p_ptr);
+
+    std::shared_ptr<Field<scalar> > pB_ptr( std::make_shared<Field<scalar> >( mesh, 0, "pB" ) );
+    auto& pB = (*pB_ptr);
+
+    poisson pEqn(p_ptr);
+    poisson pBEqn(pB_ptr);
+
+    {
+        //initial conditions
+        for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2; i++ )
+        {
+          for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2; k++ )
+          {
+
+            for( int j=settings::m()/2; j<(mesh.nj()-settings::m()/2); j++ )
+            { 
+
+
+//		    scalar z = (k-settings::m()/2)*mesh.dz()+mesh.origin().z();
+                    scalar x = (i-settings::m()/2)*mesh.dx()+mesh.origin().x();
+		    scalar y = (j-settings::m()/2)*mesh.dy()+mesh.origin().y();
+
+		    scalar psi = u0*std::exp(-(((y-0.5)*(y-0.5))/(delta*delta)))*(std::cos(8.0*pi*x)+std::cos(20.0*pi*x));
+		    scalar psi2 = u0*std::exp(-(((y-1.5)*(y-1.5))/(delta*delta)))*(std::cos(8.0*pi*x)+std::cos(20.0*pi*x));
+//                    scalar k0 = 2.0*3.1415926536 / (14.0*h0/2.0);
+//                    scalar a=0.005;
+//                    scalar b=0.4177;
+
+		if ( y < 1.0 )
+		{
+		    U(i, j, k) = vector((u0*std::tanh(((4.0 * y) - 2.0)/delta)), 0, 0);
+                    U(i, j, k).x() += ( cn * mesh.dy() * psi);
+                    U(i, j, k).y() += -cn * mesh.dx() * psi;
+		    B(i, j, k) = vector((u0*std::tanh(((4.0 * y) - 2.0)/delta)), 0, 0);
+		}
+		else
+		{  
+		    U(i, j, k) = vector(-u0*std::tanh(( (4.0 * y) - 6.0)/delta), 0, 0);
+                    U(i, j, k).x() += ( cn * mesh.dy() * psi2);
+                    U(i, j, k).y() += -cn * mesh.dx() * psi2;
+		    B(i, j, k) = vector(-u0*std::tanh(( (4.0 * y) - 6.0)/delta), 0, 0);
+		}
+	    }
+
+//                    U(i, j, k).x() += ( cn * mesh.dy() * psi);
+//		    U(i, j, k).y() += -cn * mesh.dx() * psi;
+
+//                    U(i, j, k).z() += a*u0/2.0 *dist(eng) / std::cosh(2.0*z/h0) / std::cosh(2.0*z/h0);
+//		    B(i, j, k) = vector( u0/2.0*std::tanh( 2.0*z / h0 ), 0, 0 );
+
+//                    T(i, j, k) = t0/2.0*std::tanh( 2.0*z / h0 ); //  + a*t0/2.0 *dist(eng) / std::cosh(2.0*z/h0) / std::cosh(2.0*z/h0);
+	    }
+        }
+        
+
+
+        U.correctBoundaryConditions();
+        p.correctBoundaryConditions();
+//        T.correctBoundaryConditions();
+	B.correctBoundaryConditions();
+
+    	U.write(settings::zoneName()+"/data", "U");
+  //      T.write(settings::zoneName()+"/data", "T");
+        p.write(settings::zoneName()+"/data", "p");
+	B.write(settings::zoneName()+"/data", "B");
+    }
+
+   
+    scalar RK4[4] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0 };
+
+    std::ofstream data( settings::zoneName()+"/"+"dat.dat");
+
+//    if( parallelCom::master() )
+//    {   
+//        scalar epsPlus = 6.5e-4*std::pow(u0,3)/h0;
+//        scalar LbMinus = std::pow( std::pow( mu/rho0,3 )/epsPlus, 0.25 )*std::pow(Pr,-0.5);
+//
+//        std::cout<<"Suggested mesh spacing: " << 1.5*LbMinus << " . Actual grid spacing: " << std::max(mesh.dx(), std::max( mesh.dy(), mesh.dz() ) )<<std::endl;
+//    } 
+
+    boost::timer::cpu_timer timer;
+
+    while( time.run() )
+    {
+        double start = timer.elapsed().wall / 1e9;
+//        scalar cfl = settings::restart() ? 0.5 : std::min( (0.5-0.0001)/1000.0 * time.timeStep(), 0.5 ); 
+//        time.dt() = cfl*std::min(mesh.dx(), std::min( mesh.dy(), mesh.dz() ) )/(u0/Mach+u0/2.0);
+
+        time++;
+
+        if( parallelCom::master() )
+            {   
+                std::cout << "Step: " << time.timeStep() << ".              Time: " << time.curTime() << std::endl;
+                std::cout << "Elapsed wall time for timestep: " <<  timer.elapsed().wall / 1e9 - start <<std::endl;
+  //              std::cout << "drhoz: " << drhodz << "pBbottom: " << (t0/2.0*g[2]*rho0) <<std::endl;
+	    }   
+
+            tools::CFL( U, mesh );
+        
+        
+        #include "UEqn.H" 
+//        #include "TEqn.H"
+//	#include "BEqn.H"
+
+	mesh.write(settings::zoneName()+"/data");
+        U.write(settings::zoneName()+"/data", "U");
+//	T.write(settings::zoneName()+"/data", "T");
+	p.write(settings::zoneName()+"/data", "p");
+//	B.write(settings::zoneName()+"/data", "B");
+
+//        if( time.writelog() )
+//        {
+//	    reuseTmp<vector> upt( mesh );
+//            std::shared_ptr<Field<vector> > up( upt() );
+//            reuseTmp<scalar> Tpt( mesh );
+//            std::shared_ptr<Field<scalar> > Tp( Tpt() );
+//            #include "layerAverage.H"
+//
+//            reuseTmp<scalar> epst( mesh );
+//            std::shared_ptr<Field<scalar> > eps( epst() );
+//            *eps = 2.0*( ex::grad( up ) && ex::grad( up ) )*mu/rho0;
+//            (*eps).write("e");
+//
+//            reuseTmp<scalar> xit( mesh );
+//            std::shared_ptr<Field<scalar> > xi( xit() );
+//            *xi = 2.0*( ex::grad( Tp ) & ex::grad( Tp ) )*ka/rho0;
+//            (*xi).write("x");
+//
+//            reuseTmp<vector> vortt( mesh );
+//            std::shared_ptr<Field<vector> > vort( vortt() );
+//            *vort = ex::curl( U );
+//            (*vort).write("v");
+//     
+//            reuseTmp<vector> gpt( mesh );
+//            std::shared_ptr<Field<vector> > gp( gpt() );
+//            *gp = ex::grad(p);
+//            (*gp).write("z");
+//
+//            mesh.writeSOS();
+//	}
+    }
+        
+
+    std::cout<< timer.elapsed().wall / 1e9 <<std::endl;
+    
+    mesh.write(settings::zoneName()+"/data");
+
+    return 0;
+}
diff --git a/app/shearLayerWC/TEqn.H b/app/shearLayerWC/TEqn.H
new file mode 100644
index 0000000..19063d4
--- /dev/null
+++ b/app/shearLayerWC/TEqn.H
@@ -0,0 +1,38 @@
+{
+    reuseTmp<scalar> drhoTt( mesh );
+    std::shared_ptr<Field<scalar> > drhoT( drhoTt() );
+  
+    reuseTmp<scalar> rhoTct( mesh );
+    std::shared_ptr<Field<scalar> > rhoTc( rhoTct() );
+
+    reuseTmp<scalar> rhoTOldt( mesh );
+    std::shared_ptr<Field<scalar> > rhoTOld( rhoTOldt() );
+
+    *rhoTc = rho*T;
+    *rhoTOld = rho*T;
+
+    for( int rk=0; rk<4; rk++ )
+    {
+        *drhoT =
+        time.dt() *
+        (
+           -ex::div( rho*U, T )
+           +ex::laplacian( ka, T )
+        );
+
+//        drhoT->correctBoundaryConditions();
+
+        *rhoTc = rhoTc+RK4[rk]*drhoT;
+
+        if( rk<3 )
+        {
+            T = rhoTOld/rho + (int(rk/2)+1)/2.0 * (drhoT/rho);
+        }
+        else
+        {
+            T=rhoTc/rho;
+        }
+    
+        T.correctBoundaryConditions();
+    }
+}
diff --git a/app/shearLayerWC/UEqn.H b/app/shearLayerWC/UEqn.H
new file mode 100644
index 0000000..025635c
--- /dev/null
+++ b/app/shearLayerWC/UEqn.H
@@ -0,0 +1,42 @@
+{
+    reuseTmp<vector> drhoUt( mesh );
+    std::shared_ptr<Field<vector> > drhoU( drhoUt() );
+
+    reuseTmp<vector> rhoUOldt( mesh );
+    std::shared_ptr<Field<vector> > rhoUOld( rhoUOldt() );
+
+    reuseTmp<vector> rhoUct( mesh );
+    std::shared_ptr<Field<vector> > rhoUc( rhoUct() );
+
+    *rhoUc = rho*U;
+    *rhoUOld = rho*U;
+
+//    std::shared_ptr<Field<tensor> > gradU( fdc::grad( U ) );
+//    gradU->correctBoundaryConditions();
+
+    for( int rk=0; rk<4; rk++ )
+    {
+        *drhoU =
+        time.dt() *
+        (
+           -ex::div( rho*U, U )
+           -ex::grad( p )
+           +ex::laplacian( mu, U ) 
+//           +mu*ex::div( ( dev2(transpose(gradU) ) ) )
+           -( (rho*T)*g )  
+        );
+
+        *rhoUc = rhoUc+RK4[rk]*drhoU;
+
+        if( rk<3 )
+        {
+            U = rhoUOld/rho + (int(rk/2)+1)/2.0 * (drhoU/rho);
+        }
+        else
+        {
+            U=rhoUc/rho;
+        }    
+
+        U.correctBoundaryConditions();
+    }
+}
diff --git a/app/shearLayerWC/hydrostaticP.H b/app/shearLayerWC/hydrostaticP.H
new file mode 100644
index 0000000..9f636e3
--- /dev/null
+++ b/app/shearLayerWC/hydrostaticP.H
@@ -0,0 +1,39 @@
+{
+    for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2-1; i++ )
+    {
+        for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2; j++ )
+        {
+            double pHydc=0.0;
+
+            for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2-1; k++ )
+            {
+                pHydc += T(i, j, k)*mesh.dz()*g[2]*rho0;
+            }
+
+            std::vector<scalar> cumTot( parallelCom::worldSize(), scalar(0.0) );
+            cumTot[parallelCom::myProcNo()] = pHydc;
+
+            all_reduce( cumTot, std::plus<scalar>() );
+
+            pHydc = 0.0;
+            double aboveProcs=0.0;
+
+            for( int pk=parallelCom::k()+1; pk<parallelCom::nk(); pk++ )
+            {
+                aboveProcs += cumTot[parallelCom::procNo( parallelCom::i(), parallelCom::j(), pk )];
+            } 
+
+            for( int k=mesh.nk()-2-settings::m()/2; k>=settings::m()/2; k-- )
+            {
+                pHydc += T(i, j, k)*mesh.dz()*g[2]*rho0;
+                pHyd(i, j, k) = pHydc; // + aboveProcs;
+            }
+
+        }
+    }
+//add field and implement null boundary condition
+//    pHyd.correctBoundaryConditions();
+
+    p += pHyd;
+}
+
diff --git a/app/shearLayerWC/layerAverage.H b/app/shearLayerWC/layerAverage.H
new file mode 100644
index 0000000..0e7087e
--- /dev/null
+++ b/app/shearLayerWC/layerAverage.H
@@ -0,0 +1,52 @@
+{
+    for( int k=0; k<mesh.nk(); k++ )
+    {
+        std::vector<vector> proc_avgU( parallelCom::worldSize(), 0.0 );
+        std::vector<scalar> proc_avgT( parallelCom::worldSize(), 0.0 );
+        std::vector<int> proc_N( parallelCom::worldSize(), 0 );
+        
+        for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2-1; j++ )
+        {
+            for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2-1; i++ )
+            {
+                proc_avgU[parallelCom::myProcNo()] += U(i, j, k);
+                proc_avgT[parallelCom::myProcNo()] += T(i,j,k);
+                proc_N[parallelCom::myProcNo()] += 1;
+            }
+        }
+
+        all_reduce(proc_avgU, std::plus<vector>());
+        all_reduce(proc_avgT, std::plus<scalar>());
+        all_reduce(proc_N, std::plus<int>());
+        
+        vector avgU(0,0,0);
+        scalar avgT=0.0;
+        int N=0;
+
+        for( int proci=0; proci<parallelCom::ni(); proci++ )
+        {
+            for( int procj=0; procj<parallelCom::nj(); procj++ )
+            {
+                avgU+=proc_avgU[parallelCom::procNo(proci, procj, parallelCom::k())];
+                avgT+=proc_avgT[parallelCom::procNo(proci, procj, parallelCom::k())];
+                N+=proc_N[parallelCom::procNo(proci, procj, parallelCom::k())];
+            }
+        }
+
+        avgU/=N;
+        avgT/=N;
+ 
+        for( int j=0; j<mesh.nj(); j++ )
+        {   
+            for( int i=0; i<mesh.ni(); i++ )
+            {
+                (*up)(i, j, k) = U(i, j, k)-avgU;
+                (*Tp)(i, j, k) = T(i, j, k)-avgT;
+            }
+        }
+
+        //up->correctBoundaryConditions();
+        //Tp->correctBoundaryConditions();
+    }
+}
+
diff --git a/app/shearLayerWC/pEqn.H b/app/shearLayerWC/pEqn.H
new file mode 100644
index 0000000..cd6ff37
--- /dev/null
+++ b/app/shearLayerWC/pEqn.H
@@ -0,0 +1,4 @@
+{ 
+    p = (c*c/2.0/rho0)*rho*rho;
+}
+
diff --git a/app/shearLayerWC/rhoEqn.H b/app/shearLayerWC/rhoEqn.H
new file mode 100644
index 0000000..4d4334c
--- /dev/null
+++ b/app/shearLayerWC/rhoEqn.H
@@ -0,0 +1,38 @@
+{
+    reuseTmp<scalar> drhot( mesh );
+    std::shared_ptr<Field<scalar> > drho( drhot() );
+  
+    reuseTmp<scalar> rhoct( mesh );
+    std::shared_ptr<Field<scalar> > rhoc( rhoct() );
+
+    reuseTmp<scalar> rhoOldt( mesh );
+    std::shared_ptr<Field<scalar> > rhoOld( rhoOldt() );
+
+    *rhoOld = rho;
+    *rhoc = rho;
+
+    for( int rk=0; rk<4; rk++ )
+    {
+        *drho =
+        time.dt() *
+        (
+           -ex::div( rho*U )
+        );
+
+        //drho->correctBoundaryConditions();
+
+        *rhoc = rhoc+RK4[rk]*drho;
+
+        if( rk<3 )
+        {
+            rho = rhoOld + (int(rk/2)+1)/2.0 * drho;
+        }
+        else
+        {
+            rho=rhoc;
+        }
+
+        rho.correctBoundaryConditions(); 
+    }
+
+}
diff --git a/app/shearLayerWC/shearLayer.cpp b/app/shearLayerWC/shearLayer.cpp
new file mode 100644
index 0000000..3c04145
--- /dev/null
+++ b/app/shearLayerWC/shearLayer.cpp
@@ -0,0 +1,247 @@
+#include "Field/Field.h"
+#include <iostream>
+#include "Types/vector/vector.h"
+#include "ex/grad/grad.h"
+#include "ex/div/div.h"
+#include "ex/curl/curl.h"
+#include "ex/laplacian/laplacian.h"
+#include "Time/Time.h"
+#include "parallelCom/parallelCom.h"
+#include "settings/settings.h"
+#include "BC/BC.h"
+#include <functional>
+#include "Tools/Tools.h"
+
+#include <boost/timer/timer.hpp>
+
+#include <cmath>
+#include <random>
+
+int main(int argc, char* argv[])
+{ 
+    settings::process( argc, argv );  
+    Time time( 0.01, 1000000, 100 ); //args: dt, endT, write interval / steps
+
+    parallelCom::decompose( settings::zoneName()+"/"+"mesh" );  
+
+    Mesh mesh( settings::zoneName()+"/"+"mesh", time );
+
+    mesh.write(settings::zoneName()+"/data"); 
+
+    scalar Mach=0.2;
+    scalar Re0=1967.0;
+    scalar Ri0=0.08;
+    scalar Pr=7;
+
+    scalar rho0=999.57;
+    scalar c=0.0417;
+
+    scalar u0=c*Mach;
+    scalar h0=0.2;//Re0*mu/rho0/u0;
+    scalar t0=0.1;//-Ri0*u0*u0/g[2]/h0;
+    scalar mu = rho0*u0*h0/Re0;
+    scalar ka = mu/Pr;
+    vector g( 0, 0, -Ri0*u0*u0/t0/h0 );
+
+    Field<vector> U( mesh, vector(1,0,0), "U"  );  
+    Field<scalar> T( mesh, 0, "T" );
+    Field<scalar> p( mesh,0, "p" );
+    Field<scalar> pHyd( mesh,0, "pHyd" );
+    Field<scalar> rho( mesh, rho0, "rho" );
+    
+//-------------------------------------------
+//    Time time( 5e-9, 1000000, scalar(6*h0/u0) );
+//
+//    int ni=840*1.5;
+//    int nj=102*1.5;
+//    int nk=900*1.5;
+// 
+//    parallelCom::decompose( ni, nj, nk ); 
+//
+//    Mesh mesh( ni, nj, nk, 14.0*h0, 1.7*h0, 15*h0, time );
+//    
+//
+//    time.dt() = 0.5*std::min(mesh.dx(), std::min( mesh.dy(), mesh.dz() ) )/(u0/Mach);
+//
+//    if( parallelCom::master() )
+//    {
+//        std::cout<<"Speed of sound: "<<u0/Mach<<"\nCFL_c: ";
+//        std::cout<<u0/Mach*time.dt()/std::min(mesh.dx(), std::min( mesh.dy(), mesh.dz() ) );
+//        std::cout<< "\nh0: " << h0 << "\nu0: "<< u0 << "\nt0: " << t0 <<std::endl<<std::endl;
+//    }
+//-------------------------------------------
+
+
+    //set BCs 
+//-------------------------------------------
+//    U.bc()(top,0) = BCtype(neumann,0.0);
+//    U.bc()(bottom,0) = BCtype(neumann,0.0);
+//    U.bc()(top,1) = BCtype(neumann,0.0);
+//    U.bc()(bottom,1) = BCtype(neumann,0.0);
+//    U.bc()(top,2) = BCtype(dirichlet,0.0);
+//    U.bc()(bottom,2) = BCtype(dirichlet,0.0); 
+//
+//    T.bc()(top) = BCtype(neumann,0.0);
+//    T.bc()(bottom) = BCtype(neumann,0.0);
+//
+    scalar drhodz=(std::pow( ( (c*c/2.0/rho0)*rho0*rho0 + (-t0/2.0*g[2]*rho0)*0.1e-8 )  / (c*c/2.0/rho0), 0.5 ) - rho0)/0.1e-8;
+
+//    rho.bc()(top) = BCtype(neumann, drhodz);
+//    rho.bc()(bottom) = BCtype(neumann, -drhodz);
+//
+//    p.bc()(top) = BCtype(neumann,-t0/2.0*g[2]*rho0);
+//    p.bc()(bottom) = BCtype(neumann,t0/2.0*g[2]*rho0);
+//-------------------------------------------
+
+    std::default_random_engine eng( parallelCom::myProcNo() );
+    std::uniform_real_distribution<double> dist(-1.0,1.0);
+
+
+//    if( settings::restart() )
+//    {
+//        U.read("U");
+//        rho.read("r");
+//        T.read("T");
+//        p.read("p");
+//
+//        U.correctBoundaryConditions();
+//        p.correctBoundaryConditions();
+//        rho.correctBoundaryConditions(); 
+//        T.correctBoundaryConditions();
+//    }
+//    else
+    {
+        #include "pEqn.H"
+        //initial conditions
+        for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2; i++ )
+        {
+            for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2; j++ )
+            {
+                for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2; k++ )
+                {
+                    scalar z = -7.5*h0 + (k-settings::m()/2)*mesh.dz()+mesh.origin().z();
+                    scalar x = (i-settings::m()/2)*mesh.dx()+mesh.origin().x();
+
+                    scalar k0 = 2.0*3.1415926536 / (14.0*h0/2.0);
+                    scalar a=0.005;
+                    scalar b=0.4177;
+
+                    U(i, j, k) = vector( u0/2.0*std::tanh( 2.0*z / h0 ), 0, 0 );
+                    U(i, j, k).x() += a*u0/2.0 / k0 * (-std::cos(2*k0*x)+2.0*b*std::cos(k0*x))*std::tanh(2.0*z / h0)/std::cosh(2.0*z/h0);
+                    U(i, j, k).y() += a*u0/2.0 *dist(eng) / std::cosh(2.0*z/h0) / std::cosh(2.0*z/h0);
+                    U(i, j, k).z() += a*u0/2.0 *dist(eng) / std::cosh(2.0*z/h0) / std::cosh(2.0*z/h0);
+
+                    T(i, j, k) = t0/2.0*std::tanh( 2.0*z / h0 )  + a*t0/2.0 *dist(eng) / std::cosh(2.0*z/h0) / std::cosh(2.0*z/h0);
+		}
+            }
+        }
+
+        #include "hydrostaticP.H"
+    
+        for( int i=settings::m()/2; i<mesh.ni()-settings::m()/2; i++ )
+        {   
+            for( int j=settings::m()/2; j<mesh.nj()-settings::m()/2; j++ )
+            {   
+                for( int k=settings::m()/2; k<mesh.nk()-settings::m()/2; k++ )
+                {     
+                    rho(i,j,k) = std::pow( p(i,j,k) / (c*c/2.0/rho0), 0.5 );
+                }
+            }
+        } 
+
+        U.correctBoundaryConditions();
+        p.correctBoundaryConditions();
+        rho.correctBoundaryConditions(); 
+        T.correctBoundaryConditions();
+
+    	U.write(settings::zoneName()+"/data", "U");
+        rho.write(settings::zoneName()+"/data", "rho");
+        T.write(settings::zoneName()+"/data", "T");
+        p.write(settings::zoneName()+"/data", "p");
+    }
+
+   
+    scalar RK4[4] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0 };
+
+    std::ofstream data( settings::zoneName()+"/"+"dat.dat");
+
+    if( parallelCom::master() )
+    {   
+        scalar epsPlus = 6.5e-4*std::pow(u0,3)/h0;
+        scalar LbMinus = std::pow( std::pow( mu/rho0,3 )/epsPlus, 0.25 )*std::pow(Pr,-0.5);
+
+        std::cout<<"Suggested mesh spacing: " << 1.5*LbMinus << " . Actual grid spacing: " << std::max(mesh.dx(), std::max( mesh.dy(), mesh.dz() ) )<<std::endl;
+    } 
+
+    boost::timer::cpu_timer timer;
+
+    while( time.run() )
+    {
+        double start = timer.elapsed().wall / 1e9;
+        scalar cfl = settings::restart() ? 0.5 : std::min( (0.5-0.0001)/1000.0 * time.timeStep(), 0.5 ); 
+        time.dt() = cfl*std::min(mesh.dx(), std::min( mesh.dy(), mesh.dz() ) )/(u0/Mach+u0/2.0);
+
+        time++;
+
+        if( parallelCom::master() )
+        {
+            {   
+                std::cout << "Step: " << time.timeStep() << ".              Time: " << time.curTime() << std::endl;
+                std::cout << "Elapsed wall time for timestep: " <<  timer.elapsed().wall / 1e9 - start <<std::endl;
+                std::cout << "drhoz: " << drhodz << "pBbottom: " << (t0/2.0*g[2]*rho0) <<std::endl;
+	    }   
+
+            tools::CFL( U, mesh );
+        }
+        
+        #include "UEqn.H" 
+        #include "rhoEqn.H"
+        #include "pEqn.H"  
+        #include "TEqn.H"
+
+	mesh.write(settings::zoneName()+"/data");
+        U.write(settings::zoneName()+"/data", "U");
+        rho.write(settings::zoneName()+"/data", "rho");
+	T.write(settings::zoneName()+"/data", "T");
+	p.write(settings::zoneName()+"/data", "p");
+
+
+//        if( time.writelog() )
+//        {
+//	    reuseTmp<vector> upt( mesh );
+//            std::shared_ptr<Field<vector> > up( upt() );
+//            reuseTmp<scalar> Tpt( mesh );
+//            std::shared_ptr<Field<scalar> > Tp( Tpt() );
+//            #include "layerAverage.H"
+//
+//            reuseTmp<scalar> epst( mesh );
+//            std::shared_ptr<Field<scalar> > eps( epst() );
+//            *eps = 2.0*( ex::grad( up ) && ex::grad( up ) )*mu/rho0;
+//            (*eps).write("e");
+//
+//            reuseTmp<scalar> xit( mesh );
+//            std::shared_ptr<Field<scalar> > xi( xit() );
+//            *xi = 2.0*( ex::grad( Tp ) & ex::grad( Tp ) )*ka/rho0;
+//            (*xi).write("x");
+//
+//            reuseTmp<vector> vortt( mesh );
+//            std::shared_ptr<Field<vector> > vort( vortt() );
+//            *vort = ex::curl( U );
+//            (*vort).write("v");
+//     
+//            reuseTmp<vector> gpt( mesh );
+//            std::shared_ptr<Field<vector> > gp( gpt() );
+//            *gp = ex::grad(p);
+//            (*gp).write("z");
+//
+//            mesh.writeSOS();
+//	}
+    }
+        
+
+    std::cout<< timer.elapsed().wall / 1e9 <<std::endl;
+    
+//    mesh.write(settings::zoneName()+"/data");
+
+    return 0;
+}
diff --git a/examples/.DS_Store b/examples/.DS_Store
new file mode 100644
index 0000000..1ff44b8
Binary files /dev/null and b/examples/.DS_Store differ
diff --git a/examples/OT/.DS_Store b/examples/OT/.DS_Store
new file mode 100644
index 0000000..5e783bb
Binary files /dev/null and b/examples/OT/.DS_Store differ
diff --git a/examples/OT/UBC b/examples/OT/UBC
index 30205a2..6eb1e0c 100644
--- a/examples/OT/UBC
+++ b/examples/OT/UBC
@@ -1,18 +1,18 @@
-north periodic 0 0
-south periodic 0 0
-east periodic 0 0
-west periodic 0 0
-top periodic 0 0
-bottom periodic 0 0
-north periodic 0 1
-south periodic 0 1
-east periodic 0 1
-west periodic 0 1
-top periodic 0 1
-bottom periodic 0 1
-north periodic 0 2
-south periodic 0 2
-east periodic 0 2
-west periodic 0 2
-top periodic 0 2
-bottom periodic 0 2
+north periodic 0 0 
+south periodic 0 0 
+east periodic 0 0 
+west periodic 0 0 
+top periodic 0 0 
+bottom periodic 0 0 
+north periodic 0 1 
+south periodic 0 1 
+east periodic 0 1 
+west periodic 0 1 
+top periodic 0 1 
+bottom periodic 0 1 
+north periodic 0 2 
+south periodic 0 2 
+east periodic 0 2 
+west periodic 0 2 
+top periodic 0 2 
+bottom periodic 0 2 
diff --git a/examples/TG/UBC b/examples/TG/UBC
index 30205a2..6eb1e0c 100644
--- a/examples/TG/UBC
+++ b/examples/TG/UBC
@@ -1,18 +1,18 @@
-north periodic 0 0
-south periodic 0 0
-east periodic 0 0
-west periodic 0 0
-top periodic 0 0
-bottom periodic 0 0
-north periodic 0 1
-south periodic 0 1
-east periodic 0 1
-west periodic 0 1
-top periodic 0 1
-bottom periodic 0 1
-north periodic 0 2
-south periodic 0 2
-east periodic 0 2
-west periodic 0 2
-top periodic 0 2
-bottom periodic 0 2
+north periodic 0 0 
+south periodic 0 0 
+east periodic 0 0 
+west periodic 0 0 
+top periodic 0 0 
+bottom periodic 0 0 
+north periodic 0 1 
+south periodic 0 1 
+east periodic 0 1 
+west periodic 0 1 
+top periodic 0 1 
+bottom periodic 0 1 
+north periodic 0 2 
+south periodic 0 2 
+east periodic 0 2 
+west periodic 0 2 
+top periodic 0 2 
+bottom periodic 0 2 
diff --git a/examples/TG/mesh b/examples/TG/mesh
index ff64612..428a7c9 100644
--- a/examples/TG/mesh
+++ b/examples/TG/mesh
@@ -1,3 +1,3 @@
 origin: 0 0 0
-N: 64 64 64
+N: 32 32 32
 L: 6.2831853072 6.2831853072 6.2831853072
diff --git a/examples/TG/pBC b/examples/TG/pBC
index 56388c1..808bc04 100644
--- a/examples/TG/pBC
+++ b/examples/TG/pBC
@@ -1,6 +1,6 @@
-east periodic 0 0
-west periodic 0 0
-bottom periodic 0 0
-north periodic 0 0
-top periodic 0 0
-south periodic 0 0
+east periodic 0 0 
+west periodic 0 0 
+bottom periodic 0 0 
+north periodic 0 0 
+top periodic 0 0 
+south periodic 0 0 
diff --git a/examples/mhdjet/.DS_Store b/examples/mhdjet/.DS_Store
new file mode 100644
index 0000000..5135dfe
Binary files /dev/null and b/examples/mhdjet/.DS_Store differ
diff --git a/examples/mhdjet/BBC b/examples/mhdjet/BBC
new file mode 100644
index 0000000..6eb1e0c
--- /dev/null
+++ b/examples/mhdjet/BBC
@@ -0,0 +1,18 @@
+north periodic 0 0 
+south periodic 0 0 
+east periodic 0 0 
+west periodic 0 0 
+top periodic 0 0 
+bottom periodic 0 0 
+north periodic 0 1 
+south periodic 0 1 
+east periodic 0 1 
+west periodic 0 1 
+top periodic 0 1 
+bottom periodic 0 1 
+north periodic 0 2 
+south periodic 0 2 
+east periodic 0 2 
+west periodic 0 2 
+top periodic 0 2 
+bottom periodic 0 2 
diff --git a/examples/mhdjet/UBC b/examples/mhdjet/UBC
new file mode 100644
index 0000000..6eb1e0c
--- /dev/null
+++ b/examples/mhdjet/UBC
@@ -0,0 +1,18 @@
+north periodic 0 0 
+south periodic 0 0 
+east periodic 0 0 
+west periodic 0 0 
+top periodic 0 0 
+bottom periodic 0 0 
+north periodic 0 1 
+south periodic 0 1 
+east periodic 0 1 
+west periodic 0 1 
+top periodic 0 1 
+bottom periodic 0 1 
+north periodic 0 2 
+south periodic 0 2 
+east periodic 0 2 
+west periodic 0 2 
+top periodic 0 2 
+bottom periodic 0 2 
diff --git a/examples/mhdjet/mesh b/examples/mhdjet/mesh
new file mode 100644
index 0000000..6ff470f
--- /dev/null
+++ b/examples/mhdjet/mesh
@@ -0,0 +1,3 @@
+origin: 0 0 0
+N: 32 64 32
+L: 13.96 20 13.96
diff --git a/examples/mhdjet/pBBC b/examples/mhdjet/pBBC
new file mode 100644
index 0000000..29be9ae
--- /dev/null
+++ b/examples/mhdjet/pBBC
@@ -0,0 +1,6 @@
+east periodic 0 0 
+west periodic 0 0 
+bottom periodic 0 0 
+north periodic 0 0 
+top periodic 0 0 
+south periodic 0 0
diff --git a/examples/mhdjet/pBC b/examples/mhdjet/pBC
new file mode 100644
index 0000000..0422dd6
--- /dev/null
+++ b/examples/mhdjet/pBC
@@ -0,0 +1,6 @@
+east periodic 0 0
+west periodic 0 0
+bottom periodic 0 0
+top periodic 0 0
+north periodic 0 0
+south periodic 0 0
diff --git a/examples/mhdjet/rhoBC b/examples/mhdjet/rhoBC
new file mode 100644
index 0000000..962525f
--- /dev/null
+++ b/examples/mhdjet/rhoBC
@@ -0,0 +1,6 @@
+north periodic 0 0
+south periodic 0 0
+east periodic 0 0
+west periodic 0 0
+top periodic 0 0
+bottom periodic 0 0
diff --git a/examples/shearLayer/BBC b/examples/shearLayer/BBC
new file mode 100644
index 0000000..30205a2
--- /dev/null
+++ b/examples/shearLayer/BBC
@@ -0,0 +1,18 @@
+north periodic 0 0
+south periodic 0 0
+east periodic 0 0
+west periodic 0 0
+top periodic 0 0
+bottom periodic 0 0
+north periodic 0 1
+south periodic 0 1
+east periodic 0 1
+west periodic 0 1
+top periodic 0 1
+bottom periodic 0 1
+north periodic 0 2
+south periodic 0 2
+east periodic 0 2
+west periodic 0 2
+top periodic 0 2
+bottom periodic 0 2
diff --git a/examples/shearLayer/TBC b/examples/shearLayer/TBC
new file mode 100644
index 0000000..8d8e35a
--- /dev/null
+++ b/examples/shearLayer/TBC
@@ -0,0 +1,6 @@
+east periodic 0 0
+west periodic 0 0
+bottom neumann 0 0
+north periodic 0 0
+top neumann 0 0
+south periodic 0 0
diff --git a/examples/shearLayer/UBC b/examples/shearLayer/UBC
new file mode 100644
index 0000000..6e032c5
--- /dev/null
+++ b/examples/shearLayer/UBC
@@ -0,0 +1,18 @@
+north neumann 0 0
+south neumann 0 0
+east  periodic 0 0
+west  periodic 0 0
+top periodic 0 0
+bottom periodic 0 0
+north neumann 0 1
+south neumann 0 1
+east  periodic 0 1
+west  periodic 0 1
+top periodic 0 1
+bottom periodic 0 1
+north dirichlet 0 2
+south dirichlet 0 2
+east  periodic 0 2
+west  periodic 0 2
+top periodic 0 2
+bottom periodic 0 2
diff --git a/examples/shearLayer/mesh b/examples/shearLayer/mesh
new file mode 100644
index 0000000..746cca6
--- /dev/null
+++ b/examples/shearLayer/mesh
@@ -0,0 +1,3 @@
+origin: 0 0 0
+N: 64 128 8 
+L: 1.0 2.0 0.125
diff --git a/examples/shearLayer/pBBC b/examples/shearLayer/pBBC
new file mode 100644
index 0000000..4a894f0
--- /dev/null
+++ b/examples/shearLayer/pBBC
@@ -0,0 +1,6 @@
+east periodic 0 0
+west periodic 0 0
+bottom periodic 0 0 
+north periodic 0 0
+top periodic 0 0
+south periodic 0 0
diff --git a/examples/shearLayer/pBC b/examples/shearLayer/pBC
new file mode 100644
index 0000000..c85943a
--- /dev/null
+++ b/examples/shearLayer/pBC
@@ -0,0 +1,6 @@
+east periodic 0 0
+west periodic 0 0
+bottom periodic 0 0 
+north neumann 0 0
+top periodic 0 0
+south neumann 0 0
diff --git a/findFFTMPI/FindFFTMPI.cmake b/findFFTMPI/FindFFTMPI.cmake
deleted file mode 100644
index 3476e48..0000000
--- a/findFFTMPI/FindFFTMPI.cmake
+++ /dev/null
@@ -1,64 +0,0 @@
-# - Find the PFFT library
-#
-# Usage:
-#   find_package(PFFT [REQUIRED] [QUIET] )
-#     
-# It sets the following variables:
-#   PFFT_FOUND               ... true if PFFT is found on the system
-#   PFFT_LIBRARIES           ... full path to PFFT library
-#   PFFT_INCLUDES            ... PFFT include directory
-#
-# The following variables will be checked by the function
-#   PFFT_USE_STATIC_LIBS    ... if true, only static libraries are found
-#   PFFT_ROOT               ... if set, the libraries are exclusively searched
-#                               under this path
-#   PFFT_LIBRARY            ... PFFT library to use
-#   PFFT_INCLUDE_DIR        ... PFFT include directory
-#
-
-#If environment variable PFFTDIR is specified, it has same effect as PFFT_ROOT
-if( NOT FFTMPI_ROOT AND ENV{FFTMPIDIR} )
-  set( FFTMPI_ROOT $ENV{FFTMPIDIR} )
-endif()
-
-# Check if we can use PkgConfig
-find_package(PkgConfig)
-
-#Determine from PKG
-if( PKG_CONFIG_FOUND AND NOT FFTMPI_ROOT )
-  pkg_check_modules( PKG_FFTMPI "fftmpi" QUIET )
-endif()
-
-#Check whether to search static or dynamic libs
-set( CMAKE_FIND_LIBRARY_SUFFIXES_SAV ${CMAKE_FIND_LIBRARY_SUFFIXES} )
-
-if( ${FFTMPI_USE_STATIC_LIBS} )
-  set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_STATIC_LIBRARY_SUFFIX} )
-else()
-  set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_SHARED_LIBRARY_SUFFIX} )
-endif()
-
-#find libs
-find_library(
-  FFTMPI_LIBRARIES
-  NAMES "libfft3dmpi"
-  PATHS ${FFTMPI_ROOT} ${PKG_FFTMPI_PREFIX} ${PKG_FFTMPI_LIBRARY_DIRS} ${LIB_INSTALL_DIR}
-  PATH_SUFFIXES "lib" "lib64" "so"
-)
-
-#find includes
-find_path(
-  FFTMPI_INCLUDES
-  NAMES "fft3d.h"
-  PATHS ${FFTMPI_ROOT} ${PKG_FFTMPI_PREFIX} ${PKG_FFTMPI_INCLUDE_DIRS} ${INCLUDE_INSTALL_DIR}
-  PATH_SUFFIXES "include"
-)
-
-set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_FIND_LIBRARY_SUFFIXES_SAV} )
-
-include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(FFTMPI DEFAULT_MSG
-                                  FFTMPI_INCLUDES FFTMPI_LIBRARIES)
-
-mark_as_advanced(FFTMPI_INCLUDES FFTMPI_LIBRARIES)
-
diff --git a/findHEFFTE/findHEFFTE.cmake b/findHEFFTE/findHEFFTE.cmake
new file mode 100644
index 0000000..be81c52
--- /dev/null
+++ b/findHEFFTE/findHEFFTE.cmake
@@ -0,0 +1,27 @@
+# - Find heFFTe 
+
+
+#find libs
+find_library(
+  HEFFTE_LIBRARIES
+  NAMES "libheffte"
+  PATHS  ${CMAKE_INSTALL_PREFIX} 
+  PATH_SUFFIXES "lib" "lib64"
+)
+
+#find includes
+find_path(
+  HEFFTE_INCLUDES
+  NAMES "heffte.h"
+  PATHS ${CMAKE_INSTALL_PREFIX} 
+  PATH_SUFFIXES "include"
+)
+
+
+set( CMAKE_FIND_LIBRARY_SUFFIXES ${CMAKE_FIND_LIBRARY_SUFFIXES_SAV} )
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(Heffte DEFAULT_MSG
+                                  HEFFTE_INCLUDES HEFFTE_LIBRARIES)
+
+mark_as_advanced(HEFFTE_INCLUDES HEFFTE_LIBRARIES)
diff --git a/src/BC/BC.h b/src/BC/BC.h
index 9f7e15f..f447c5d 100644
--- a/src/BC/BC.h
+++ b/src/BC/BC.h
@@ -54,7 +54,7 @@ class BC
     }
 
     //access
-    scalar& val()
+    scalar val()
     {
         return val_;
     }
diff --git a/src/BC/dirichletBC.hxx b/src/BC/dirichletBC.hxx
index f9a31da..75aee4a 100644
--- a/src/BC/dirichletBC.hxx
+++ b/src/BC/dirichletBC.hxx
@@ -4,6 +4,7 @@ void dirichletBC<T>::update
     Field<T>* fld_ptr
 )
 {
+    //bottom
     if( this->dir() == bottom && parallelCom::k() == 0 )
     {
         for( int i=0; i<fld_ptr->ni(); i++ )
@@ -11,20 +12,116 @@ void dirichletBC<T>::update
             for( int j=0; j<fld_ptr->nj(); j++ )
             {  
                 int k=settings::m()/2;
-                fld_ptr->operator()(i,j,k) =  this->val();
-            }
+                component(fld_ptr->operator()(i, j, k), this->comp()) = this->val();
+		
+		for( int l=0; l<settings::m()/2; l++ )
+                {
+                  component(fld_ptr->operator()(i, j, k-(l+1)), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i, j, k+(l+1) ), this->comp());
+                }	    
+
+	    }
         }
     }
-    else if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 )
+
+
+    //top
+    if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 )
     { 
         for( int i=0; i<fld_ptr->ni(); i++ )
         {      
             for( int j=0; j<fld_ptr->nj(); j++ )
             {   
                 int k=fld_ptr->nk()-settings::m()/2-1;
-                fld_ptr->operator()(i,j,k) =  this->val();
+		component(fld_ptr->operator()(i, j, k), this->comp()) = this->val();
+
+                for( int l=0; l<settings::m()/2; l++ )
+                {   
+                  component(fld_ptr->operator()(i, j, k+(l+1)), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i, j, k-(l+1) ), this->comp());
+                } 
+
             }
         }
     }
 
+
+    //east
+    if( this->dir() == east && parallelCom::i() == 0 )
+    {   
+        for( int j=0; j<fld_ptr->nj(); j++ )
+        {          
+            for( int k=0; k<fld_ptr->nk(); k++ )
+            {   
+                int i=settings::m()/2;
+                component(fld_ptr->operator()(i, j, k), this->comp()) = this->val();
+
+                for( int l=0; l<settings::m()/2; l++ )
+		{   
+                  component(fld_ptr->operator()(i-(l+1), j, k), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i+(l+1), j, k ), this->comp());
+		}   
+
+            }   
+        }   
+    }
+
+
+    //west
+    if( this->dir() == west && parallelCom::i() == parallelCom::ni()-1 )
+    {   
+        for( int j=0; j<fld_ptr->nj(); j++ )
+        {          
+            for( int k=0; k<fld_ptr->nk(); k++ )
+            {   
+                int i=fld_ptr->ni()-settings::m()/2-1;
+                component(fld_ptr->operator()(i, j, k), this->comp()) = this->val();
+
+                for( int l=0; l<settings::m()/2; l++ )
+                {   
+                  component(fld_ptr->operator()(i+(l+1), j, k), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i-(l+1), j, k ), this->comp());
+                }   
+
+            }   
+        }   
+    } 
+    
+
+    //south
+    if( this->dir() == south && parallelCom::j() == 0 ) 
+    {   
+        for( int i=0; i<fld_ptr->ni(); i++ )
+        {   
+            for( int k=0; k<fld_ptr->nk(); k++ )
+            {   
+                int j=settings::m()/2;
+                component(fld_ptr->operator()(i, j, k), this->comp()) = this->val();
+    
+                for( int l=0; l<settings::m()/2; l++ )
+                {   
+                  component(fld_ptr->operator()(i, j-(l+1), k), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i, j+(l+1), k ), this->comp());
+                }    
+
+            }   
+        }   
+    }
+
+
+    //north
+    if( this->dir() == north && parallelCom::j() == parallelCom::nj()-1 )
+    {   
+        for( int i=0; i<fld_ptr->ni(); i++ )
+        {    
+            for( int k=0; k<fld_ptr->nk(); k++ )
+            {   
+                int j=fld_ptr->nj()-settings::m()/2-1;
+                component(fld_ptr->operator()(i, j, k), this->comp()) = this->val();
+
+                for( int l=0; l<settings::m()/2; l++ )
+                {   
+                  component(fld_ptr->operator()(i, j+(l+1), k), this->comp()) = 2.0*(l+1)*this->val() - component(fld_ptr->operator()(i, j-(l+1), k ), this->comp());
+                }   
+
+            }   
+        }   
+    } 
+
+
 }   
diff --git a/src/BC/neumannBC.hxx b/src/BC/neumannBC.hxx
index 947b45f..43386ed 100644
--- a/src/BC/neumannBC.hxx
+++ b/src/BC/neumannBC.hxx
@@ -4,27 +4,131 @@ void neumannBC<T>::update
     Field<T>* fld_ptr
 )
 {
-    if( this->dir() == bottom && parallelCom::k() == 0 )
-    {
-        for( int i=0; i<fld_ptr->ni(); i++ )
+	//bottom
+        if( this->dir() == bottom && parallelCom::k() == 0 ) 
+        {   
+            for( int i=0; i<fld_ptr->ni(); i++ )
+            {   
+                for( int j=0; j<fld_ptr->nj(); j++ )
+                {   
+                    int k=settings::m()/2; 
+
+                    component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i, j, k+1), this->comp())
+                    -this->val()*fld_ptr->mesh().dz();
+
+                    for( int l=0; l<settings::m()/2; l++ )
+                    {
+                        component(fld_ptr->operator()(i, j, k-(l+1)), this->comp()) = component(fld_ptr->operator()(i, j, k+(l+1)), this->comp())
+                       -2.0*(l+1)*this->val()*fld_ptr->mesh().dz();
+		    }
+                }   
+            }   
+        }   
+
+	//top
+        if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 )
+        {   
+            for( int i=0; i<fld_ptr->ni(); i++ )
+            {   
+                for( int j=0; j<fld_ptr->nj(); j++ )
+                {   
+                    int k=fld_ptr->nk()-settings::m()/2-1;
+                
+		    component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i, j, k-1), this->comp()) 
+                    + this->val()*fld_ptr->mesh().dz();
+
+                    for( int l=0; l<settings::m()/2; l++ )
+                    {   
+                        component(fld_ptr->operator()(i, j, k+(l+1)), this->comp()) = component(fld_ptr->operator()(i, j, k-(l+1)), this->comp())
+                       -2.0*(l+1)*this->val()*fld_ptr->mesh().dz();
+                     } 
+                }   
+            }   
+        }      
+
+
+	//east
+        if( this->dir() == east && parallelCom::i() == 0 ) 
         {   
             for( int j=0; j<fld_ptr->nj(); j++ )
-            {  
-                int k=settings::m()/2;
-                fld_ptr->operator()(i,j,k) =  this->val();
-            }
+            {   
+                for( int k=0; k<fld_ptr->nk(); k++ )
+                {   
+                    int i=settings::m()/2; 
+
+                    component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i+1, j, k), this->comp())
+                    -this->val()*fld_ptr->mesh().dx();
+
+                    for( int l=0; l<settings::m()/2; l++ )
+                    {   
+                        component(fld_ptr->operator()(i-(l+1), j, k), this->comp()) = component(fld_ptr->operator()(i+(l+1), j, k), this->comp())
+                       -2.0*(l+1)*this->val()*fld_ptr->mesh().dx();
+                    }   
+                }   
+            }   
         }
-    }
-    else if( this->dir() == top && parallelCom::k() == parallelCom::nk()-1 )
-    {
-        for( int i=0; i<fld_ptr->ni(); i++ )
+
+        //west
+        if( this->dir() == west && parallelCom::i() == parallelCom::ni()-1 )
         {   
             for( int j=0; j<fld_ptr->nj(); j++ )
-            {  
-                int k=fld_ptr->nk()-settings::m()/2-1;
-                fld_ptr->operator()(i,j,k) =  this->val(); 
+            {   
+                for( int k=0; k<fld_ptr->nk(); k++ )
+                {   
+                    int i=fld_ptr->ni()-settings::m()/2-1;
+    
+                    component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i-1, j, k), this->comp()) 
+                    + this->val()*fld_ptr->mesh().dx();
+
+                    for( int l=0; l<settings::m()/2; l++ )
+                    {   
+                        component(fld_ptr->operator()(i+(l+1), j, k), this->comp()) = component(fld_ptr->operator()(i-(l+1), j, k), this->comp())
+                       -2.0*(l+1)*this->val()*fld_ptr->mesh().dx();
+                    }   
+                }   
+            }   
+        }
+
+        //south
+        if( this->dir() == south && parallelCom::j() == 0 ) 
+        {   
+            for( int i=0; i<fld_ptr->ni(); i++ )
+            {   
+                for( int k=0; k<fld_ptr->nk(); k++ )
+                {   
+                    int j=settings::m()/2; 
+
+                    component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i, j+1, k), this->comp())
+                    -this->val()*fld_ptr->mesh().dy();
+
+                    for( int l=0; l<settings::m()/2; l++ )
+                    {
+                        component(fld_ptr->operator()(i, j-(l+1), k), this->comp()) = component(fld_ptr->operator()(i, j+(l+1), k), this->comp())
+                       -2.0*(l+1)*this->val()*fld_ptr->mesh().dy();
+                    }
+                }
             }
         }
-    }
-}
 
+	//north
+        if( this->dir() == north && parallelCom::j() == parallelCom::nj()-1 )
+        {
+            for( int i=0; i<fld_ptr->ni(); i++ )
+            {
+                for( int k=0; k<fld_ptr->nk(); k++ )
+                {
+                    int j=fld_ptr->nj()-settings::m()/2-1;
+
+                    component(fld_ptr->operator()(i, j, k), this->comp()) = component(fld_ptr->operator()(i, j-1, k), this->comp())
+                    + this->val()*fld_ptr->mesh().dy();
+
+                    for( int l=0; l<settings::m()/2; l++ )
+                    {
+                        component(fld_ptr->operator()(i, j+(l+1), k), this->comp()) = component(fld_ptr->operator()(i, j-(l+1), k), this->comp())
+                       -2.0*(l+1)*this->val()*fld_ptr->mesh().dy();
+                    }
+                }
+            }
+        }
+
+}
diff --git a/src/Field/Field.h b/src/Field/Field.h
index 57edb6c..223b8b8 100644
--- a/src/Field/Field.h
+++ b/src/Field/Field.h
@@ -18,7 +18,6 @@
 */
 #ifndef FIELD_H
 #define FIELD_H
-
 #if defined(_OPENMP)
 #include <omp.h>
 #endif
@@ -129,6 +128,7 @@ class Field
 
     void correctBoundaryConditions();  
     void write(std::string, std::string);
+//    void H5write(std::string, std::string);
     void read(std::string, std::string);
     void read(std::string, std::string, int);
     void setBC( std::string );
diff --git a/src/Mesh/Mesh.cpp b/src/Mesh/Mesh.cpp
index 4c1dd5e..043deca 100644
--- a/src/Mesh/Mesh.cpp
+++ b/src/Mesh/Mesh.cpp
@@ -18,7 +18,8 @@
 */
 #include "Mesh/Mesh.h"
 #include <sstream>
-
+#include <memory>
+#include "H5Cpp.h"
 
 Mesh::Mesh( std::string fileName, Time& runTime )
 :
@@ -103,7 +104,7 @@ Mesh::Mesh( std::string fileName, Time& runTime )
     }
 
     origin_ = vector( (ni_-1)*dx_*parallelCom::i()+ox, (nj_-1)*dy_*parallelCom::j()+oy, (nk_-1)*dz_*parallelCom::k()+oz );
-#ifdef HAVE_FFTMPI
+#ifdef HAVE_HEFFTE
     localStart_[0] = origin_[0];
     localStart_[1] = origin_[1];
     localStart_[2] = origin_[2];
@@ -324,8 +325,12 @@ void Mesh::writeCase( std::string path )
 
     fout << "VARIABLE" << std::endl;
     fout << "vector per node: velocity " << boost::format("U.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ;
+    fout << "vector per node: Uprime " << boost::format("Uprime.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ;
+    fout << "vector per node: Udash " << boost::format("Udash.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ;
     fout << "scalar per node: p  " << boost::format("p.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ;
     fout << "scalar per node: T  " << boost::format("T.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ;
+    fout << "vector per node: B " << boost::format("B.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ;
+    fout << "vector per node: J " << boost::format("J.%|04|.*****.dat") % parallelCom::myProcNo() << std::endl ;
 
 	fout << "TIME" << std::endl ;
 	fout << "time set: 1" << std::endl ;
@@ -343,7 +348,7 @@ void Mesh::writeCase( std::string path )
 
 	fout.close() ;
 
-    this->writeRestartInfo( path );
+   this->writeRestartInfo( path );
 }
 
 void Mesh::writeRestartInfo( std::string path )
@@ -364,13 +369,14 @@ void Mesh::writeRestartInfo( std::string path )
 }
 
 
+
 scalar Mesh::z( int k ) const
 {
     scalar z;
     if( hyperbollic_ )
     {
         scalar A = std::pow( s1_/s2_, 0.5 );
-#ifdef HAVE_FFTMPI
+#ifdef HAVE_HEFFTE
         int globK = k-settings::m()/2 + localStart_[2];
 #else 
         int globK = k-settings::m()/2 + (parallelCom::k())*(nk_-settings::m()-1);
@@ -438,3 +444,4 @@ void Mesh::sendMUIPoints()
     }
 }
 #endif
+
diff --git a/src/Mesh/Mesh.h b/src/Mesh/Mesh.h
index 17794a3..c655d10 100644
--- a/src/Mesh/Mesh.h
+++ b/src/Mesh/Mesh.h
@@ -28,6 +28,7 @@
 #include <fstream>
 #include <iostream>
 #include <iomanip>
+#include "H5Cpp.h"
 
 
 class Time;
@@ -57,7 +58,7 @@ class Mesh
 
     vector origin_;
 
-    Time& runTime_;
+//    Time& runTime_;
 
     Box bbox_;
     Box bboxHalo_;
@@ -67,13 +68,14 @@ class Mesh
 
     bool hyperbollic_;
 
-#ifdef HAVE_FFTMPI
+#ifdef HAVE_HEFFTE
     ptrdiff_t localStart_[3]; 
 #endif
 
     public:
     Mesh( int, int, int, scalar, scalar, scalar, Time& );
     Mesh( std::string, Time& );
+    Time& runTime_;
 
     //access
     int ni() const
@@ -156,7 +158,7 @@ class Mesh
         return bboxHalo_;
     } 
 
-#ifdef HAVE_FFTMPI
+#ifdef HAVE_HEFFTE
     ptrdiff_t* glob_n()
     {
         return glob_n_;
@@ -172,6 +174,7 @@ class Mesh
     void writeCase( std::string );
     void writeSOS();
     void writeRestartInfo( std::string );
+//    void H5write( std::string );
 };
 
 #endif
diff --git a/src/parallelCom/parallelCom.cpp b/src/parallelCom/parallelCom.cpp
index 7911b28..4f16883 100644
--- a/src/parallelCom/parallelCom.cpp
+++ b/src/parallelCom/parallelCom.cpp
@@ -60,13 +60,13 @@ void parallelCom::decompose(int ni, int nj, int nk)
 {
     int np=worldSize_;
     
-    double lsq=1e10;
+    double lsq=1e10;    
 
     for( int ii=1; ii<=np; ii++ )
     {
         for( int jj=1; jj<=np/ii; jj++ )
         {
-            if( std::fabs( double(np)/double(ii)/double(jj) - np/ii/jj ) < 1e-5  ) 
+            if( std::fabs( double(np)/double(ii)/double(jj) - np/ii/jj ) < 1e-5  )
             {
                 int kk=np/ii/jj;
 
@@ -82,6 +82,28 @@ void parallelCom::decompose(int ni, int nj, int nk)
             }
         }
     }
+//    pencil decomposition template
+
+//    int ii=1;
+//      for( int jj=1; jj<=np; jj++ )
+//        {   
+//            if( std::fabs( double(np)/double(jj) - np/jj ) < 1e-5  )   
+//            {   
+//                int kk=np/jj;
+//
+//                double lsqc = std::pow( double(nj)/double(jj), 2 ) + std::pow( double(nk)/double(kk), 2 );
+//
+//                if( lsqc < lsq )
+//               {   
+//                    lsq=lsqc;
+//                    ni_=ii;
+//                    nj_=jj;
+//                    nk_=kk;
+//                }   
+//            }   
+//        }      
+      
+    
             
     if( master() ) 
     {
diff --git a/src/poisson/poisson.cpp b/src/poisson/poisson.cpp
index 4acbe98..4833fa7 100644
--- a/src/poisson/poisson.cpp
+++ b/src/poisson/poisson.cpp
@@ -12,22 +12,17 @@ poisson::poisson( std::shared_ptr<Field<scalar> > pptr )
     nk_(pptr->mesh().nk()-settings::m()-1),
     localStart_(pptr->mesh().localStart())
 {
-fft = std::make_unique<FFT3d>(parallelCom::world(), 2);
-  //fft options
-  //tflag=0;
-  permute = 0;
+  //Define FFT grid. This can be trimmed now HeFFTe is the prefered library.
+  int NFAST = n_[0];
+  int NMID = n_[1];
+  int NSLOW = n_[2];
 
-  NFAST = n_[0];
-  NMID = n_[1];
-  NSLOW = n_[2];
-
-//  int MPI_Init( int *argc, char ***argv );
-      MPI_Comm world = parallelCom::world();
-  
   nfast = NFAST;
   nmid = NMID;
   nslow = NSLOW;
 
+  MPI_Comm world = parallelCom::world();
+
   npfast  = parallelCom::ni();
   npmid   = parallelCom::nj();
   npslow  = parallelCom::nk();
@@ -42,7 +37,7 @@ fft = std::make_unique<FFT3d>(parallelCom::world(), 2);
   klo = (int) 1.0*(ipslow)*nslow/npslow;
   khi = (int) 1.0*(ipslow+1)*nslow/npslow - 1;
 
-  //define local coords for wavenumbers
+  //Define local coords for wavenumber calculations
   inxlo = static_cast<int> (1.0 * ipfast * nfast / npfast);
   inxhi = static_cast<int> (1.0 * (ipfast+1) * nfast / npfast) - 1;
 
@@ -54,6 +49,7 @@ fft = std::make_unique<FFT3d>(parallelCom::world(), 2);
 
   nfft_in = (inxhi-inxlo+1) * (inyhi-inylo+1) * (inzhi-inzlo+1);
 
+//Indexing for MPI running
 {
   for (int iproc = 0; iproc < parallelCom::worldSize(); iproc++) {
     if (parallelCom::myProcNo() != iproc) continue;
@@ -69,9 +65,9 @@ fft = std::make_unique<FFT3d>(parallelCom::world(), 2);
           iglobal = inxlo + ilocal;
           jglobal = inylo + jlocal;
           kglobal = inzlo + klocal;
-          printf("Value (%d,%d,%d) on proc %d \n",
-               iglobal,jglobal,kglobal,
-               parallelCom::myProcNo());
+//          printf("Value (%d,%d,%d) on proc %d \n",
+//               iglobal,jglobal,kglobal,
+//               parallelCom::myProcNo());
 	if (m == 0)
 	{
 	localStart_[0] = iglobal;
@@ -83,29 +79,88 @@ fft = std::make_unique<FFT3d>(parallelCom::world(), 2);
   if (parallelCom::myProcNo() < parallelCom::worldSize()) MPI_Send(&tmp,0,MPI_INT,parallelCom::myProcNo()+1,0,world);
   } 
 }
+	//Sanity check for fft coordinates. Uncomment for debugging    
+    //printf("Assigned local start coordinates (%d,%d,%d) on proc %d \n",
+    //            localStart_[0],localStart_[1],localStart_[2],parallelCom::myProcNo());
+    //printf("local range ijk (%d,%d,%d,%d,%d,%d)\n", inxlo,inxhi,inylo,inyhi,inzlo,inzhi);
 
-std::cout<<ni_<<" "<<nj_<<" "<<nk_<<" ";
-std::cout<<ilo<<" "<<ihi<<" "<<jlo<<" "<<jhi<<" "<<klo<<" "<<khi<<" "<<npfast<<" "<<npmid<<" "<<npslow<<" "<<nfast<<" "<<nmid<<" "<<nslow<<" ";
 
-    //fft plan
-    fft->setup(nfast, nmid, nslow,             // 3d verion
-          ilo, ihi, jlo, 
-          jhi, klo, khi,
-          ilo, ihi, jlo, 
-          jhi, klo, khi,
-          permute, fftsize, sendsize, recvsize);
+ 	//initialize local fft domain
 
+    up_N = std::make_unique<heffte::box3d<>>(std::array<int, 3> ({inxlo, inylo, inzlo})  ,std::array<int, 3> ({inxhi, inyhi, inzhi}));
+    N = up_N.get();
+
+    up_NDCT = std::make_unique<heffte::box3d<>>(std::array<int, 3> ({inxlo, inylo, inzlo})  ,std::array<int, 3> ({inxhi, inyhi, inzhi}));
+    NDCT = up_NDCT.get();
+
+	//Construct and allocate memory for fft and dct
+	//HeFFTe::fftw_cos1 selects FFTW DCT TYPE REDFT00
+    dct = std::make_unique<heffte::fft3d<heffte::backend::fftw_cos1>>(*NDCT, *NDCT, world);
+    fft2d = std::make_unique<heffte::fft3d<heffte::backend::fftw>>(*NDCT, *NDCT, world);
+    dct->nullify_executor(0); 
+    dct->nullify_executor(2);
+    fft2d->nullify_executor(1);
+
+
+    fftx = std::make_unique<heffte::fft3d<heffte::backend::fftw>>(*N, *N, world);
+    fftx->nullify_executor(1);
+    fftx->nullify_executor(2);
+
+    fft = std::make_unique<heffte::fft3d<heffte::backend::fftw>>(*N, *N, world);
+
+	//plan 3d DCT-FFT-FFT
+    //perform 1d r2r_DCT
+    //translate domain from std::double to std::complex<double>
+    //perform 2d r2c_DFT
+    //convolute and solve
+    //perform 2d c2r_IDFT
+    //translate domain from std::complex<double> to std::double
+    //perform 1d r2r_DCT
+    //std::double cosphi
+    //std::complex<double> dctphi
+    //std::complex<double> dctphihat
+    //heffte::fft1d<heffte::backend::fftw_cos>>
+    //heffte::fft2d<heffte::backend::fftw>>
 
-    //fft allocate memory
-    std::cout<<fftsize<<" ";
-        
-    up_phi_ = std::make_unique_for_overwrite<FFT_SCALAR[]>(2*fftsize);
-    phi_ = up_phi_.get();
 
-    up_k2_ = std::make_unique<double[]>(2*fftsize);
+    //Create and allocate memory for FFT/DCT inputs & outputs.
+    up_k2_ = std::make_unique<double[]>((fft->size_inbox()));
     k2_ = up_k2_.get();
+
+    up_k2c_ = std::make_unique<double[]>((2*dct->size_inbox()));
+    k2c_ = up_k2c_.get();
+
+    up_k2cy_ = std::make_unique<double[]>((2*dct->size_inbox()));
+    k2cy_ = up_k2cy_.get(); 
+
+    up_dctworkspace_ = std::make_unique<std::vector<double>>(2*dct->size_workspace());  
+    dctworkspace_ = up_dctworkspace_.get();
+
+    up_cosphihat_ = std::make_unique<std::vector<double>>((2*dct->size_outbox()));
+    cosphihat_ = up_cosphihat_.get();
+
+    up_cosphi_ = std::make_unique<std::vector<double>>((2*dct->size_inbox()));
+    cosphi_ = up_cosphi_.get();
+
+    up_cosphicomplex_ = std::make_unique<std::vector<std::complex<double>>>((2*fft2d->size_outbox()));
+    cosphicomplex_ = up_cosphicomplex_.get();
     
-    wavenumbers_();
+    up_phi_ = std::make_unique<std::vector<double>>((fft->size_inbox()));
+    phi_ = up_phi_.get();
+   
+    up_phihat_ = std::make_unique<std::vector<std::complex<double>>>((fft->size_outbox()));
+    phihat_ = up_phihat_.get();
+   
+    up_psi_ = std::make_unique<std::vector<double>>((fftx->size_inbox()));
+    psi_ = up_psi_.get();
+   
+    up_psihat_ = std::make_unique<std::vector<std::complex<double>>>((fftx->size_outbox()));
+    psihat_ = up_psihat_.get();
+
+
+
+    wavenumbersFFT_();
+    wavenumbersDCT_();    
 }
 
 
@@ -118,51 +173,81 @@ void poisson::rhs( std::shared_ptr<Field<scalar> > f )
         {
             for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ )
 	    {
-                phi_[l] = (double) f->operator()(i, j, k); 
-		l++;
-		phi_[l] = 0.0;
+                phi_[0][l] = f->operator()(i, j, k) ; 
 		l++;
             }   
+        }
+    }   
+}
+
+void poisson::rhsDCT( std::shared_ptr<Field<scalar> > f )
+{ 
+
+    int l=0;
+    for( int k=(settings::m()/2); k<(pptr_->nk()-1-settings::m()/2); k++ )
+    {   
+        for( int j=(settings::m()/2); j<(pptr_->nj()-1-settings::m()/2); j++ )
+        {   
+            for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ )
+            {                  
+		cosphi_[0][l] = f->operator()(i, j, k) ; 
+                l++;
+            }   
         }   
-    }
-    fft->compute(phi_,phi_,1);
+    } 
+
 }
 
+void poisson::x_FFT( std::shared_ptr<Field<scalar> > f )
+{  
+    int l=0;
+    for( int k=(settings::m()/2)+1; k<(settings::m()/2)+2; k++ )
+    {
+        for( int j=(settings::m()/2); j<(pptr_->nj()-1-settings::m()/2); j++ )
+        {
+            for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ )
+            {
+                psi_[0][l] = f->operator()(i, j, k) ;
+                l++;
+            }  
+        }
+    } 
+    fftx->forward(psi_->data(), psihat_->data());
+    
+    std::ofstream data("fftxdat.dat");
+    for(size_t ll=0; ll<((fftx->size_outbox())/nk_); ll++) {
+    data<<std::real(psihat_[0][ll])<<" "<<std::setprecision(15)<<std::imag(psihat_[0][ll])<<" "<<std::setprecision(15)<<std::endl;    
+    }
+}
 
 void poisson::solve()
 {
-   int l=0;
-    for( int k=0; k<nk_; k++ )
-    {   
-        for( int j=0; j<nj_; j++ )
+    fft->forward(phi_->data(), phihat_->data());
+    
+    //dealiasing
+//    scalar kk = (n_[0]*n_[1]*n_[2]) * 2.0/3.0;
+//
+    int l=0;
+        for( int k=0; k<nk_; k++ )
         {   
-	    for( int i=0; i<ni_; i++ )
-	    {
+            for( int j=0; j<nj_; j++ )
+            {   
+                for( int i=0; i<ni_; i++ )
+                {
 		if( k2_[l] < tools::eps )
                 {
-                    phi_[l] = 0.0;
+		    phihat_[0][l] = std::complex<double>(0.0, 0.0);
 		}
                 else
-                {
-                    phi_[l] = -phi_[l] / k2_[l];
-		}
-		l++;	
-           
-		if( k2_[l] < tools::eps )
-                {   
-                    phi_[l] = 0.0;
-                }   
-                else
-                {   
-                    phi_[l] = -phi_[l] / k2_[l];
-                }   
-                l++; 
+                { 
+                     phihat_[0][l] = std::complex<double>(( -std::real(phihat_[0][l]) / k2_[l]) , -std::imag(phihat_[0][l]) / k2_[l]);
+                } 
 
+		l++;	
 	    }
-        }
+       }
     }
-
-    fft->compute(phi_,phi_,-1);
+    fft->backward(phihat_->data(), phi_->data());
     
     l=0;
     for( int k=(settings::m()/2); k<(pptr_->nk()-1-settings::m()/2); k++ )
@@ -171,32 +256,98 @@ void poisson::solve()
         {
             for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ )
 	    {
-		pptr_->operator()(i,j,k) = phi_[l];
+		pptr_->operator()(i,j,k) = phi_[0][l] / (n_[0]*n_[1]*n_[2]);
     		l++;
-		l++;
 	
 	    }
         }
     }
+}
+
+
+void poisson::solveDCT()
+{
+    dct->forward(cosphi_->data(), cosphihat_->data());
+
+       int l=0;
+        for( int k=0; k<nk_; k++ )
+        {   
+            for( int j=0; j<nj_; j++ )
+            {   
+                for( int i=0; i<ni_; i++ )
+                {   
+                {   
+                     cosphihat_[0][l] =  (cosphihat_[0][l]) * k2cy_[l];
+                }   
+                l++;
+                }   
+            }   
+        }
+ 
+    fft2d->forward(cosphihat_->data(), cosphicomplex_->data());
+
+	l=0;
+        for( int k=0; k<nk_; k++ )
+        {
+            for( int j=0; j<nj_; j++ )
+            {
+                for( int i=0; i<ni_; i++ )
+                {
+		if( k2c_[l] < tools::eps )
+                {
+		    cosphicomplex_[0][l] = std::complex<double>(0.0 , 0.0);
+		}
+                else
+                {
+                     cosphicomplex_[0][l] =( std::complex<double>(( -std::real(cosphicomplex_[0][l]) / k2c_[l]) , -std::imag(cosphicomplex_[0][l]) / k2c_[l]  ) );
+                }   
+		l++;
+                }   
+            }   
+        } 
+  
+    fft2d->backward(cosphicomplex_->data(), cosphihat_->data());
+
+    dct->backward(cosphihat_->data(), cosphi_->data(), dctworkspace_->data());
+
+    scalar pmax=0.0;
+    l=0;
+    for( int k=(settings::m()/2); k<(pptr_->nk()-1-settings::m()/2); k++ )
+    {
+        for( int j=(settings::m()/2); j<(pptr_->nj()-1-settings::m()/2); j++ )
+        {
+            for( int i=(settings::m()/2); i<(pptr_->ni()-1-settings::m()/2); i++ )
+            {
+//		pmax = std::max( pmax, (cosphi_[0][l] / (2.0*n_[0]*n_[1]*(n_[2]-1) )));
+                pptr_->operator()(i,j,k) = cosphi_[0][l] / (2*(n_[0])*(n_[1]-1)*n_[2]);
+                l++;
+
+            }
+        }
+    }
+//    std::cout<<" "<<dct->size_inbox()<<std::endl;
+//    std::cout<<" "<<dct->size_outbox()<<std::endl;
+//std::cout<<" "<<fft2d->size_outbox()<<std::endl;
+//std::cout<<" "<<fft2d->size_outbox()<<std::endl;
+//std::cout<<" "<<fft->size_outbox()<<std::endl;
+//std::cout<<" "<<fft->size_inbox()<<std::endl;
 
 }
 
-void poisson::wavenumbers_()
+void poisson::wavenumbersFFT_()
 {
     Mesh& m = pptr_->mesh();
     auto pi = tools::pi;
 
-    printf("Assigned local start coordinates (%d,%d,%d) on proc %d \n",
-                localStart_[0],localStart_[1],localStart_[2],parallelCom::myProcNo());
     
     scalar kx, ky, kz; //TODO: will be complex for non-central schemes
     int l=0;
-    for( int k=0; k<nk_; k++ )
-    {   
-        for( int j=0; j<nj_; j++ )
-        {   
-            for( int i=0; i<ni_; i++ )
+	for( int k=0; k<nk_; k++ )
+	{
+	    for( int j=0; j<nj_; j++ )
 	    {
+		for( int i=0; i<ni_; i++ )
+		{
 
 		if (2*( i + localStart_[0] )<n_[0])
                     kx = 2.0*pi*( i + localStart_[0]  )/m.lx();
@@ -239,14 +390,84 @@ void poisson::wavenumbers_()
                 
                 k2_[l] = kpx*kpx + kpy*kpy + kpz*kpz;
                 l++;
-		k2_[l] = kpx*kpx + kpy*kpy + kpz*kpz;
-		l++;
 	    }
         }
     }
 }
 
+void poisson::wavenumbersDCT_()
+{
+    Mesh& m = pptr_->mesh();
+    auto pi = tools::pi;
+
+
+    scalar kx, ky, kz; //TODO: will be complex for non-central schemes
+    int l=0;
+    kx=0;
+    ky=0;
+    kz=0;
+        for( int k=0; k<nk_; k++ )
+        {
+            for( int j=0; j<nj_; j++ )
+            {
+                for( int i=0; i<ni_; i++ )
+                {
+                if (2*( i + localStart_[0] )<n_[0])
+                    kx = 2.0*pi*( i + localStart_[0]  )/m.lx();
+                else 
+                    kx = pi*(n_[0] - i - localStart_[0] )/m.lx();
+
+                if (2*( j + localStart_[1] ) <n_[1])
+                    ky = pi*( j + localStart_[1] )/m.ly();
+                else 
+                    ky = pi*(n_[1] - j - localStart_[1] )/m.ly();
+
+
+                    kz = 2.0*pi*( k + localStart_[2] )/m.lz();
+
+                
+
+		//modified wavenumbers
+                scalar kpx=0.0;
+                scalar kpy=0.0;
+                scalar kpz=0.0;
+                if(kx > tools::eps )
+                {   
+                    for( int ll=0; ll<(settings::m()/2); ll++ )
+                    {   
+                        kpx+=kx*settings::coef()[ll]/(ll+1)*std::sin((ll+1)*kx*m.dx())/(kx*m.dx());
+                    }   
+                }   
+                if(ky > tools::eps )
+                {   
+                    for( int ll=0; ll<(settings::m()/2); ll++ )
+                    {   
+                        kpy+=ky*settings::coef()[ll]/(ll+1)*std::sin((ll+1)*ky*m.dy())/(ky*m.dy());
+                    }   
+                }   
+                if(kz > tools::eps )
+                {   
+                    for( int ll=0; ll<(settings::m()/2); ll++ )
+                    {   
+                        kpz+=kz*settings::coef()[ll]/(ll+1)*std::sin((ll+1)*kz*m.dz())/(kz*m.dz());
+                    }   
+                }    
+
+		k2cy_[l] = (kpy/(n_[1]-1)) * (kpy/(n_[1]-1) ); 
+                k2c_[l] = kpx*kpx + kpz*kpz;
+                l++;
+            }   
+        }   
+    }   
+}
+
+//Destructor no longer needed, can remove
 poisson::~poisson()
 {
-	free(phi_);
+	free(N);
+	fftw_free(phi_);
+	fftw_free(cosphi_);
+	fftw_free(phihat_);
+
 }
+
diff --git a/src/poisson/poisson.h b/src/poisson/poisson.h
index 08ef128..78ed9c5 100644
--- a/src/poisson/poisson.h
+++ b/src/poisson/poisson.h
@@ -21,19 +21,25 @@
 #define POISSON_H
 
 #include <stdint.h>
-#include "fft3d.h"
+#include "heffte.h"
+#include <heffte_common.h>
+
 #include "Mesh/Mesh.h"
 #include "Field/Field.h"
 #include "Tools/Tools.h"
 #include <mpi.h>
 #include <memory>
+#include <complex>
 
-using namespace FFTMPI_NS;
+using namespace heffte;
 
 class poisson
 {
     protected:
-    std::unique_ptr<FFT3d> fft;
+    std::unique_ptr<heffte::fft3d<heffte::backend::fftw>> fft;
+    std::unique_ptr<heffte::fft3d<heffte::backend::fftw>> fft2d;
+    std::unique_ptr<heffte::fft3d<heffte::backend::fftw_cos1>> dct;
+    std::unique_ptr<heffte::fft3d<heffte::backend::fftw>> fftx;
     std::shared_ptr<Field<scalar> > pptr_; 
     int fftsize;
     unsigned ni_, nj_, nk_;
@@ -51,27 +57,63 @@ class poisson
     int nfft_in;
     int nfft_out;
     int ilo,ihi,jlo,jhi,klo,khi;
-    int permute,sendsize,recvsize;
+    int permute,exchange,sendsize,recvsize;
     int flag,niter,tflag;
     double tmax; 
     int memorysize;
+ 
+    std::unique_ptr<std::vector<std::complex<double>>> up_phihat_;
+    std::vector<std::complex<double>>* phihat_; 
+
+    std::unique_ptr<std::vector<double>> up_cosphi_;
+    std::vector<double>* cosphi_;
+
+    std::unique_ptr<std::vector<std::complex<double>>> up_cosphicomplex_;
+    std::vector<std::complex<double>>* cosphicomplex_;
+
+    std::unique_ptr<std::vector<double>> up_cosphihat_;
+    std::vector<double>* cosphihat_;
+
+    std::unique_ptr<std::vector<double>> up_dctworkspace_;
+    std::vector<double>* dctworkspace_;
+
+    std::unique_ptr<std::vector<double>> up_phi_;
+    std::vector<double>* phi_;
+
+    std::unique_ptr<std::vector<std::complex<double>>> up_psihat_;
+    std::vector<std::complex<double>>* psihat_;
+
+    std::unique_ptr<std::vector<double>> up_psi_;
+    std::vector<double>* psi_;
+
+    std::unique_ptr<double[]> up_k2cy_;
+    double* k2cy_;
 
-    std::unique_ptr<FFT_SCALAR[]> up_phi_;
-    FFT_SCALAR* phi_;
     std::unique_ptr<double[]> up_k2_;
     double* k2_;
 
+    std::unique_ptr<double[]> up_k2c_;
+    double* k2c_;
+
+    heffte::box3d<>* N;
+    std::unique_ptr<heffte::box3d<>> up_N;
+
+    heffte::box3d<>* NDCT;
+    std::unique_ptr<heffte::box3d<>> up_NDCT;
 
     public:
+
     poisson( std::shared_ptr<Field<scalar> > );
     ~poisson();
-    void wavenumbers_();
-    
+    void wavenumbersFFT_();
+    void wavenumbersDCT_();
+    void x_FFT( std::shared_ptr<Field<scalar> > f );    
     void rhs(std::shared_ptr<Field<scalar> > f );
+    void rhsDCT(std::shared_ptr<Field<scalar> > f );
     void solve();
+    void solveDCT();
 };
 
 #endif
 
 
-
diff --git a/src/reuseTmp/reuseTmp.hxx b/src/reuseTmp/reuseTmp.hxx
index be178e9..4b8b1aa 100644
--- a/src/reuseTmp/reuseTmp.hxx
+++ b/src/reuseTmp/reuseTmp.hxx
@@ -32,7 +32,8 @@ reuseTmp<T>::reuseTmp
     for( int i=0; i<=20; i++ )
     {
         assert( i!=20 );
-        if( tmp_[i].unique() || tmp_[i] == NULL)
+//        if( tmp_[i].unique() || tmp_[i] == NULL)
+	if( tmp_[i].use_count() == 1 || tmp_[i] == NULL)
         {
             if( tmp_[i] == NULL )
             {