diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp
index 322317f8..0e569250 100644
--- a/src/pmpo_MPMesh.cpp
+++ b/src/pmpo_MPMesh.cpp
@@ -74,6 +74,7 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){
 
 
 void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
+    Kokkos::Timer timer;
     int numVtxs = p_mesh->getNumVertices();
     int numElms = p_mesh->getNumElements();
     auto numMPs = p_MPs->getCount();
@@ -84,7 +85,9 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
 
     auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
     auto mpTgtPos = p_MPs->getData<MPF_Tgt_Pos_XYZ>();
-    auto MPs2Elm = p_MPs->getData<MPF_Tgt_Elm_ID>();;
+    auto MPs2Elm = p_MPs->getData<MPF_Tgt_Elm_ID>();
+    auto MPs2Proc = p_MPs->getData<MPF_Tgt_Proc_ID>();
+    auto elm2Process = p_mesh->getElm2Process();
     
     if(printVTPIndex>=0) {
       printVTP_mesh(printVTPIndex);
@@ -132,6 +135,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
                 }
                 if(closestElm<0){
                     MPs2Elm(mp) = iElm;
+                    if (elm2Process.size() > 0) 
+                        MPs2Proc(mp) = elm2Process(iElm);
                     break;
                 }else{
                     iElm = closestElm;
@@ -193,6 +198,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
         fprintf(pFile,"        </DataArray>\n      </Lines>\n    </Piece>\n  </PolyData>\n</VTKFile>\n");
         fclose(pFile);
     }
+    pumipic::RecordTime("PolyMPO_CVTTrackingElmCenterBased", timer.seconds());
 }
 
 void MPMesh::T2LTracking(Vec2dView dx){
@@ -258,17 +264,37 @@ void MPMesh::T2LTracking(Vec2dView dx){
     p_MPs->parallel_for(T2LCalc,"T2lTrackingCalc");
 }
 
+bool getAnyIsMigrating(bool isMigrating) {
+  Kokkos::Timer timer;
+  int comm_rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
+  int comm_size;
+  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
+
+  bool anyIsMigrating = false;
+  MPI_Allreduce(&isMigrating, &anyIsMigrating, 1, MPI_C_BOOL, MPI_LOR, MPI_COMM_WORLD);
+  pumipic::RecordTime("PolyMPO_getAnyIsMigrating", timer.seconds());
+  return anyIsMigrating;
+}
+
 void MPMesh::push(){
+  Kokkos::Timer timer;
   p_mesh->computeRotLatLonIncr();
   sphericalInterpolation<MeshF_RotLatLonIncr, MPF_Rot_Lat_Lon_Incr>(*this);
   p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ
 
-  CVTTrackingElmCenterBased(); // move to Tgt_XYZ
+  bool anyIsMigrating = false;
+  do {
+    CVTTrackingElmCenterBased(); // move to Tgt_XYZ
+    p_MPs->updateMPSlice<MPF_Cur_Pos_XYZ, MPF_Tgt_Pos_XYZ>(); // Tgt_XYZ becomes Cur_XYZ
+    p_MPs->updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon, MPF_Tgt_Pos_Rot_Lat_Lon>(); // Tgt becomes Cur
+    p_MPs->updateMPElmID(); //update mpElm IDs slices
+    bool isMigrating = p_MPs->migrate();
+    anyIsMigrating = getAnyIsMigrating(isMigrating);
+  } 
+  while (anyIsMigrating);
 
-  p_MPs->updateMPSlice<MPF_Cur_Pos_XYZ, MPF_Tgt_Pos_XYZ>(); // Tgt_XYZ becomes Cur_XYZ
-  p_MPs->updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon, MPF_Tgt_Pos_Rot_Lat_Lon>(); // Tgt becomes Cur
-  p_MPs->rebuild(); //rebuild pumi-pic
-  p_MPs->updateMPElmID(); //update mpElm IDs slices
+  pumipic::RecordTime("PolyMPO_push", timer.seconds());
 }
 
 void MPMesh::printVTP_mesh(int printVTPIndex){
diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp
index df280d48..4991c80e 100644
--- a/src/pmpo_c.cpp
+++ b/src/pmpo_c.cpp
@@ -595,7 +595,7 @@ void polympo_getMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, cons
   Kokkos::deep_copy(arrayHost, vtxField);
 }
 
-void polympo_setMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) {
+void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) {
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
@@ -642,3 +642,29 @@ void polympo_push_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
   ((polyMPO::MPMesh*)p_mpmesh) ->push();
 }
+
+void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array){
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; 
+  PMT_ALWAYS_ASSERT(p_mesh->meshEditable());
+  kkViewHostU<const int*> arrayHost(array,nCells); 
+
+  //check the size
+  PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements());
+
+  Kokkos::View<int*> owningProc("owningProc",nCells);
+  Kokkos::deep_copy(owningProc, arrayHost);
+  p_mesh->setOwningProc(owningProc);
+}
+
+void polympo_enableTiming_f(){
+  pumipic::EnableTiming();
+}
+
+void polympo_summarizeTime_f(){
+  pumipic::SummarizeTime();
+}
+
+void polympo_setTimingVerbosity_f(int v){
+  pumipic::SetTimingVerbosity(v);
+}
\ No newline at end of file
diff --git a/src/pmpo_c.h b/src/pmpo_c.h
index 23200b42..db8e35cd 100644
--- a/src/pmpo_c.h
+++ b/src/pmpo_c.h
@@ -50,6 +50,7 @@ void polympo_getMeshNumElms_f(MPMesh_ptr p_mpmesh, int & numElms);
 void polympo_setMeshNumEdgesPerElm_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array);
 void polympo_setMeshElm2VtxConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int nCells, const int* array);
 void polympo_setMeshElm2ElmConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int nCells, const int* array);
+void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array);
 
 //Mesh fields
 void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* xArray, const double* yArray, const double* zArray);
@@ -58,10 +59,15 @@ void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const
 void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double* latitude);
 void polympo_setMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
 void polympo_getMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
-void polympo_setMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
+void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
 void polympo_getMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
 
 // calculations
 void polympo_push_f(MPMesh_ptr p_mpmesh);
+
+// Timing
+void polympo_enableTiming_f();
+void polympo_summarizeTime_f();
+void polympo_setTimingVerbosity_f(int v);
 }
 #endif
diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90
index d2051ee3..4753e179 100644
--- a/src/pmpo_fortran.f90
+++ b/src/pmpo_fortran.f90
@@ -483,8 +483,8 @@ subroutine polympo_getMeshOnSurfVeloIncr(mpMesh, nComps, nVertices, array) &
   !> @param nVertices(in) numVertices
   !> @param array(in) input mesh velocity 2D array (2,numVtx)
   !---------------------------------------------------------------------------
-  subroutine polympo_setMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
-             bind(C, NAME='polympo_setMeshOnSurfDispIncr_f')
+  subroutine polympo_setMeshVtxOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
+             bind(C, NAME='polympo_setMeshVtxOnSurfDispIncr_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
     integer(c_int), value :: nComps, nVertices
@@ -507,6 +507,19 @@ subroutine polympo_getMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
   !---------------------------------------------------------------------------
+  !> @brief set the owning process array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nCells(in) number of cells
+  !> @param array(in) input mesh cell to process array
+  !---------------------------------------------------------------------------
+  subroutine polympo_setOwningProc(mpMesh, nCells, array) &
+    bind(C, NAME='polympo_setOwningProc_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nCells
+    type(c_ptr), intent(in), value :: array
+  end subroutine
+  !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
   !>        longitude, update the MP slices
   !>        MPs MUST have rotated flag set to True(>0)
@@ -517,6 +530,25 @@ subroutine polympo_push(mpMesh) &
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
   end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief enable timing functions
+  !---------------------------------------------------------------------------
+  subroutine polympo_enableTiming() bind(C, NAME='polympo_enableTiming_f')
+    use :: iso_c_binding
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief print summary of timings collected
+  !---------------------------------------------------------------------------
+  subroutine polympo_summarizeTime() bind(C, NAME='polympo_summarizeTime_f')
+    use :: iso_c_binding
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief set verbosity of timings collected
+  !---------------------------------------------------------------------------
+  subroutine polympo_setTimingVerbosity(v) bind(C, NAME='polympo_setTimingVerbosity_f')
+    use :: iso_c_binding
+    integer(c_int), value :: v
+  end subroutine
   end interface
   contains
   !---------------------------------------------------------------------------
diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp
index c0c53b69..c1570ef8 100644
--- a/src/pmpo_materialPoints.cpp
+++ b/src/pmpo_materialPoints.cpp
@@ -101,6 +101,33 @@ void MaterialPoints::finishRebuild() {
   rebuildFields.ongoing = false;
 }
 
+bool MaterialPoints::migrate() {
+  Kokkos::Timer timer;
+  auto MPs2Elm = getData<MPF_Tgt_Elm_ID>();
+  auto MPs2Proc = getData<MPF_Tgt_Proc_ID>();
+
+  IntView new_elem("new_elem", MPs->capacity());
+  IntView new_process("new_process", MPs->capacity());
+  IntView isMigrating("isMigrating", 1);
+
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  auto setMigrationFields = PS_LAMBDA(const int& e, const int& mp, const bool& mask) {
+    if (mask) {
+      new_elem(mp) = MPs2Elm(mp);
+      new_process(mp) = MPs2Proc(mp);
+      if (new_process(mp) != rank) isMigrating(0) = 1;
+    }
+  };
+  parallel_for(setMigrationFields, "setMigrationFields");
+  MPs->migrate(new_elem, new_process);
+
+  if (getOpMode() == polyMPO::MP_DEBUG)
+    printf("Material point migration: %f\n", timer.seconds());
+  pumipic::RecordTime("PolyMPO_migrate", timer.seconds());
+  return pumipic::getLastValue(isMigrating) > 0;
+}
+
 bool MaterialPoints::rebuildOngoing() { return rebuildFields.ongoing; }
 
 void MaterialPoints::setAppIDFunc(IntFunc getAppIDIn) { getAppID = getAppIDIn; }
diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp
index 33beff2f..ae93a2e4 100644
--- a/src/pmpo_materialPoints.hpp
+++ b/src/pmpo_materialPoints.hpp
@@ -21,6 +21,7 @@ using defaultSpace = Kokkos::DefaultExecutionSpace::memory_space;
 typedef int mp_flag_t;
 typedef int mp_id_t;
 typedef int  mp_elm_id_t;
+typedef int mp_proc_id_t;
 typedef double mp_sclr_t[1];//TODO
 typedef vec2d_t mp_vec2d_t;
 typedef double mp_vec3d_t[3];//TODO
@@ -49,7 +50,8 @@ enum MaterialPointSlice {
   MPF_Stress_Div,
   MPF_Shear_Traction,
   MPF_Constv_Mdl_Param,
-  MPF_MP_APP_ID
+  MPF_MP_APP_ID,
+  MPF_Tgt_Proc_ID
 };
 
 enum Operating_Mode{
@@ -76,7 +78,8 @@ const static std::map<MaterialPointSlice, std::pair<int,MeshFieldIndex>>
                            {MPF_Stress_Div,      {3,MeshF_Unsupported}},
                            {MPF_Shear_Traction,  {3,MeshF_Unsupported}},
                            {MPF_Constv_Mdl_Param,{12,MeshF_Unsupported}},
-                           {MPF_MP_APP_ID,       {1,MeshF_Invalid}}};
+                           {MPF_MP_APP_ID,       {1,MeshF_Invalid}},
+                           {MPF_Tgt_Proc_ID,     {0,MeshF_Invalid}}}; //TODO: What does this integer mean?
 
 const static std::vector<std::pair<MaterialPointSlice, MaterialPointSlice>>
         mpSliceSwap = {{MPF_Cur_Elm_ID, MPF_Tgt_Elm_ID},
@@ -101,7 +104,8 @@ typedef MemberTypes<mp_flag_t,              //MP_Status
                     mp_vec3d_t,             //MP_Stress_Div
                     mp_vec3d_t,             //MP_Shear_Traction
                     mp_constv_mdl_param_t,  //MP_Constv_Mdl_Param
-                    mp_id_t                 //MP_APP_ID
+                    mp_id_t,                //MP_APP_ID
+                    mp_proc_id_t            //MPF_Tgt_Proc_ID
                     >MaterialPointTypes;
 typedef ps::ParticleStructure<MaterialPointTypes> PS;
 
@@ -134,7 +138,9 @@ class MaterialPoints {
     void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID, Kokkos::View<const int*> addedMPMask);
     void finishRebuild();
     bool rebuildOngoing();
-    
+
+    bool migrate();
+
     template<int mpSliceIndex, typename mpSliceData>
     typename std::enable_if<mpSliceData::rank==1>::type
     setRebuildMPSlice(mpSliceData mpSliceIn);
@@ -202,6 +208,7 @@ class MaterialPoints {
         updateMPSlice<MPF_Cur_Pos_XYZ,MPF_Tgt_Pos_XYZ>();
     }
     void updateRotLatLonAndXYZ2Tgt(const double radius){
+        Kokkos::Timer timer;
         auto curPosRotLatLon = MPs->get<MPF_Cur_Pos_Rot_Lat_Lon>();
         auto tgtPosRotLatLon = MPs->get<MPF_Tgt_Pos_Rot_Lat_Lon>();
         auto tgtPosXYZ = MPs->get<MPF_Tgt_Pos_XYZ>();
@@ -227,6 +234,7 @@ class MaterialPoints {
             PMT_ALWAYS_ASSERT(false);
         } 
         ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude"); 
+        pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds()); 
     } 
 
     template <int index>
diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp
index 78b3a7be..269aa4cb 100644
--- a/src/pmpo_mesh.cpp
+++ b/src/pmpo_mesh.cpp
@@ -40,6 +40,7 @@ namespace polyMPO{
     }
     
     void Mesh::computeRotLatLonIncr(){
+        Kokkos::Timer timer;
         PMT_ALWAYS_ASSERT(geomType_ == geom_spherical_surf);
        
         auto dispIncr = getMeshField<MeshF_OnSurfDispIncr>();
@@ -53,6 +54,11 @@ namespace polyMPO{
             rotLatLonIncr(iVtx, 0) = dispIncr(iVtx, 1)/sphereRadius;
             rotLatLonIncr(iVtx, 1) = dispIncr(iVtx, 0)/(sphereRadius * std::cos(lat(iVtx)));
         });
+        pumipic::RecordTime("PolyMPO_computeRotLatLonIncr", timer.seconds());
+    }
+
+    IntView Mesh::getElm2Process() {
+        return owningProc_;
     }
 
 } // namespace polyMPO
diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp
index 3f94256a..9e82c741 100644
--- a/src/pmpo_mesh.hpp
+++ b/src/pmpo_mesh.hpp
@@ -2,6 +2,7 @@
 #define POLYMPO_MESH_H
 
 #include "pmpo_utils.hpp"
+#include <pumipic_kktypes.hpp>
 
 namespace polyMPO{
 
@@ -66,7 +67,9 @@ class Mesh {
     //IntView nEdgesPerElm_;
     IntVtx2ElmView elm2VtxConn_;
     IntElm2ElmView elm2ElmConn_;
-  
+
+    IntView owningProc_;
+
     //start of meshFields
     DoubleVec3dView vtxCoords_;
     DoubleSclrView vtxRotLat_;
@@ -103,6 +106,8 @@ class Mesh {
     bool checkMeshType(int meshType);
     bool checkGeomType(int geomType);
 
+    IntView getElm2Process();
+
     mesh_type getMeshType() { return meshType_; }
     geom_type getGeomType() { return geomType_; }
     double getSphereRadius() { return sphereRadius_; }
@@ -129,6 +134,8 @@ class Mesh {
                                                      elm2VtxConn_ = elm2VtxConn; }
     void setElm2ElmConn(IntElm2ElmView elm2ElmConn) {PMT_ALWAYS_ASSERT(meshEdit_);
                                                      elm2ElmConn_ = elm2ElmConn; }
+    void setOwningProc(IntView owningProc) {PMT_ALWAYS_ASSERT(meshEdit_);
+                                            owningProc_ = owningProc; }
     
     void computeRotLatLonIncr();
 };
diff --git a/src/pmpo_wachspressBasis.hpp b/src/pmpo_wachspressBasis.hpp
index cff2a35a..77e1e68b 100644
--- a/src/pmpo_wachspressBasis.hpp
+++ b/src/pmpo_wachspressBasis.hpp
@@ -306,6 +306,7 @@ void getBasisByAreaGblForm_1(Vec2d MP, int numVtxs, Vec2d* vtxCoords, double* ba
 // spherical interpolation of values from mesh vertices to MPsi
 template <MeshFieldIndex mfIndex, MaterialPointSlice mpfIndex>
 void sphericalInterpolation(MPMesh& mpMesh){
+    Kokkos::Timer timer;
     auto p_mesh = mpMesh.p_mesh;
     auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
     int numVtxs = p_mesh->getNumVertices();
@@ -355,6 +356,7 @@ void sphericalInterpolation(MPMesh& mpMesh){
         }
     };
     p_MPs->parallel_for(interpolation, "interpolation");
+    pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds());
 }
 
 } //namespace polyMPO end
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 637a5b9f..500fb58e 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -48,6 +48,7 @@ pmpo_add_exe(testReconstruction testReconstruction.cpp)
 pmpo_add_exe(testTracking testTracking.cpp)
 pmpo_add_exe(testRebuild testRebuild.cpp)
 pmpo_add_exe(timeAssmblyWachspress testTiming.cpp)
+pmpo_add_exe(testMigration testMigration.cpp)
 pmpo_add_fortran_exe(testFortranInit testFortranInit.f90)
 pmpo_add_fortran_exe(testFortran testFortran.f90)
 pmpo_add_fortran_exe(testFortranReadMPAS testFortranReadMPAS.f90)
@@ -57,26 +58,26 @@ pmpo_add_fortran_exe(testFortranMPAdvection testFortranMPAdvection.f90)
 pmpo_add_fortran_exe(testFortranCreateRebuildMPs testFortranCreateRebuildMPs.f90)
 
 #add test
-find_program(MPIRUN_EXECUTABLE NAMES mpirun) 
+find_program(MPIRUN NAMES mpirun mpiexec srun) 
 
-function(pmpo_add_test testname command)
-  if(MPIRUN_EXECUTABLE)
+function(pmpo_add_test testname procs command)
+  if(MPIRUN)
     add_test(NAME ${testname} 
-             COMMAND mpirun --bind-to core -np 1 ${command} ${ARGN})
+             COMMAND ${MPIRUN} -n ${procs} ${command} ${ARGN})
   else()
     add_test(NAME ${testname} COMMAND ${command} ${ARGN})
   endif()
 endfunction()
 
-pmpo_add_test(unit_test ./unitTest)
-pmpo_add_test(mp_unit_test ./mpUnitTest)
-pmpo_add_test(test_reconstruction ./testReconstruction)
-pmpo_add_test(test_tracking ./testTracking)
-pmpo_add_test(test_rebuild ./testRebuild)
-pmpo_add_test(testFortranInit ./testFortranInit)
-pmpo_add_test(testFortran ./testFortran)
-pmpo_add_test(testFortranMPMeshModule ./testFortranMPMeshModule)
-pmpo_add_test(testFortranMPAppIDs ./testFortranMPAppIDs)
+pmpo_add_test(unit_test 1 ./unitTest)
+pmpo_add_test(mp_unit_test 1 ./mpUnitTest)
+pmpo_add_test(test_reconstruction 1 ./testReconstruction)
+pmpo_add_test(test_tracking 1 ./testTracking)
+pmpo_add_test(test_rebuild 1 ./testRebuild)
+pmpo_add_test(testFortranInit 1 ./testFortranInit)
+pmpo_add_test(testFortran 1 ./testFortran)
+pmpo_add_test(testFortranMPMeshModule 1 ./testFortranMPMeshModule)
+pmpo_add_test(testFortranMPAppIDs 1 ./testFortranMPAppIDs)
 
 #set NC file for test
 set(TEST_NC_FILE_PLANAR "${CMAKE_SOURCE_DIR}/test/sample_mpas_meshes/planar_nonuniform_cvt_for_square_673elms.nc")
@@ -86,12 +87,15 @@ if(EXISTS ${TEST_NC_FILE_SPHERICAL})
     message(STATUS "Found Spherical NC File: ${TEST_NC_FILE_SPHERICAL}") 
     set(TEST_NC_FILE ${TEST_NC_FILE_SPHERICAL})
     # this test only works for spherical NC files!
-    pmpo_add_test(testFortranMPAdvection ./testFortranMPAdvection ${TEST_NC_FILE}) 
+    pmpo_add_test(testFortranMPAdvection 1 ./testFortranMPAdvection 5 5 ${TEST_NC_FILE}) 
 else()
     set(TEST_NC_FILE ${TEST_NC_FILE_PLANAR})
 endif()
-pmpo_add_test(testFortranReadMPAS ./testFortranReadMPAS ${TEST_NC_FILE})
-pmpo_add_test(testFortranCreateRebuildMPs ./testFortranCreateRebuildMPs ${TEST_NC_FILE})
+pmpo_add_test(testFortranReadMPAS 1 ./testFortranReadMPAS ${TEST_NC_FILE})
+pmpo_add_test(testFortranCreateRebuildMPs 1 ./testFortranCreateRebuildMPs ${TEST_NC_FILE})
 
-pmpo_add_test(test_timing ./timeAssmblyWachspress 1 ${TEST_NC_FILE})
-pmpo_add_test(test_wachspress ./testWachspress ${TEST_NC_FILE})
+pmpo_add_test(test_timing 1 ./timeAssmblyWachspress 1 ${TEST_NC_FILE})
+pmpo_add_test(test_wachspress 1 ./testWachspress ${TEST_NC_FILE})
+if(MPIRUN)
+  pmpo_add_test(test_migration 8 ./testMigration ${TEST_NC_FILE})
+endif()
\ No newline at end of file
diff --git a/test/calculateDisplacement.f90 b/test/calculateDisplacement.f90
new file mode 100644
index 00000000..c499008b
--- /dev/null
+++ b/test/calculateDisplacement.f90
@@ -0,0 +1,40 @@
+subroutine calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+  use :: polympo
+  use :: readMPAS
+  use :: iso_c_binding
+  implicit none
+
+  real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
+  real(kind=MPAS_RKIND) :: maxlon, minlon, deltaLon, sphereRadius
+  integer, dimension(:), pointer :: nEdgesOnCell
+  integer, dimension(:,:), pointer :: verticesOnCell
+  integer :: i, j, nVertices, nCompsDisp
+  real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr
+  type(c_ptr) :: mpMesh
+
+  nCompsDisp = 2
+  allocate(dispIncr(nCompsDisp,nVertices))
+
+  ! check first element/cell for delta
+  maxlon = minval(lonVertex)
+  minlon = maxval(lonVertex)
+  do i = 1, nEdgesOnCell(1)
+    j = verticesOnCell(i,1)
+    if(maxlon .lt. lonVertex(j)) then
+      maxlon = lonVertex(j)
+    endif
+    if(minlon .gt. lonVertex(j)) then
+      minlon = lonVertex(j)
+    endif
+  end do
+
+  deltaLon = maxlon - minlon
+
+  do i = 1,nVertices
+    dispIncr(1,i) = sphereRadius*cos(latVertex(i))*deltaLon
+    dispIncr(2,i) = 0.0_MPAS_RKIND
+  end do
+  call polympo_setMeshVtxOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr))
+
+  deallocate(dispIncr)
+end subroutine
\ No newline at end of file
diff --git a/test/testFortran.f90 b/test/testFortran.f90
index dd953ace..92ecb2ea 100644
--- a/test/testFortran.f90
+++ b/test/testFortran.f90
@@ -75,7 +75,7 @@ program main
     end do
   end do
   call polympo_setMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
-  call polympo_setMeshOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
+  call polympo_setMeshVtxOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
 
   Mesharray = 1
   call polympo_getMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90
index a1b3fb88..d0a013f5 100644
--- a/test/testFortranMPAdvection.f90
+++ b/test/testFortranMPAdvection.f90
@@ -1,3 +1,60 @@
+module advectionTests
+  contains
+  include "calculateDisplacement.f90"
+
+  subroutine setProcWedges(mpMesh, nCells, comm_size, nEdgesOnCell, verticesOnCell, lonCell)
+    use :: polympo
+    use :: readMPAS
+    use :: iso_c_binding
+    implicit none
+
+    type(c_ptr) :: mpMesh
+    integer :: i, j, k, nCells, comm_size
+    integer, dimension(:), pointer :: owningProc, nEdgesOnCell
+    integer, dimension(:,:), pointer :: verticesOnCell
+    real(kind=MPAS_RKIND), dimension(:), pointer :: lonCell
+    real(kind=MPAS_RKIND) :: normalizedLat, min, max
+    
+    allocate(owningProc(nCells))
+    
+    min = 1000
+    max = -1000
+    do i = 1, nCells
+      if (lonCell(i) < min) min = lonCell(i)
+      if (lonCell(i) > max) max = lonCell(i)
+    end do
+
+    do i = 1, nCells
+      normalizedLat = (lonCell(i) - min) / (max - min) * .99
+      owningProc(i) = normalizedLat * comm_size
+    end do
+
+    call polympo_startMeshFill(mpMesh)
+    call polympo_setOwningProc(mpMesh, nCells, c_loc(owningProc))
+  end subroutine
+
+  subroutine runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+    use :: polympo
+    use :: readMPAS
+    use :: iso_c_binding
+    implicit none
+
+    type(c_ptr) :: mpMesh
+    integer :: i, numPush, nVertices
+    real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
+    integer, dimension(:), pointer :: nEdgesOnCell
+    integer, dimension(:,:), pointer :: verticesOnCell
+    real(kind=MPAS_RKIND) :: sphereRadius
+
+
+    do i = 1, numPush
+      call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+      call polympo_push(mpMesh)
+    end do
+
+  end subroutine
+
+end module
 !---------------------------------------------------------------------------
 !> todo add a discription
 !---------------------------------------------------------------------------
@@ -5,47 +62,52 @@ program main
   use :: polympo
   use :: readMPAS
   use :: iso_c_binding
+  use :: advectionTests
   implicit none
   include 'mpif.h'
 
   !integer, parameter :: APP_RKIND = selected_real_kind(15)
   type(c_ptr) :: mpMesh
-  integer :: ierr, self
-  integer :: argc, i, j, arglen, k
+  integer :: ierr, self, comm_size
+  integer :: argc, i, j, arglen, k, m, mpsScaleFactorPerVtx, localNumMPs
   integer :: setMeshOption, setMPOption
   integer :: maxEdges, vertexDegree, nCells, nVertices
-  integer :: nCompsDisp
   integer :: mpi_comm_handle = MPI_COMM_WORLD
-  real(kind=MPAS_RKIND) :: xc, yc, zc, radius, maxlon, minlon, deltaLon, lon
+  real(kind=MPAS_RKIND) :: xc, yc, zc, xMP, yMP, zMP, radius, lon
   real(kind=MPAS_RKIND) :: pi = 4.0_MPAS_RKIND*atan(1.0_MPAS_RKIND)
-  character (len=2048) :: filename
-  real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr
+  character (len=2048) :: filename, input
   character (len=64) :: onSphere
-  real(kind=MPAS_RKIND) :: sphereRadius, xComputed, yComputed, zComputed, latComputed, lonComputed
+  real(kind=MPAS_RKIND) :: sphereRadius
   integer, dimension(:), pointer :: nEdgesOnCell
   real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
-  real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
+  real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex, lonCell
   integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
-  integer :: numMPs 
+  integer :: numMPs, numMPsCount, numPush
   integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition, mpLatLon
-  logical :: inBound
   integer, parameter :: MP_ACTIVE = 1
   integer, parameter :: MP_INACTIVE = 0
   integer, parameter :: INVALID_ELM_ID = -1
 
   call mpi_init(ierr)
   call mpi_comm_rank(mpi_comm_handle, self, ierr)
+  call mpi_comm_size(mpi_comm_handle, comm_size, ierr)
 
   call polympo_setMPICommunicator(mpi_comm_handle)
   call polympo_initialize()
+  call polympo_enableTiming()
 
   call polympo_checkPrecisionForRealKind(MPAS_RKIND)
   argc = command_argument_count()
-  if(argc == 1) then
-    call get_command_argument(1, filename)
+  if(argc == 3) then
+    call get_command_argument(1, input)
+    read(input, '(I7)') mpsScaleFactorPerVtx
+    call get_command_argument(2, input)
+    read(input, '(I7)') numPush
+    call get_command_argument(3, filename)
   else
-    write(0, *) "Usage: ./testFortranInterpolatePush <path to the nc file>"
+    write(0, *) "Usage: ./testFortranMPAdvection <mpsScaleFactorPerVtx> <numPush> <path to the nc file>"
+    call exit(1)
   end if
 
   call readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
@@ -68,11 +130,17 @@ program main
                         xVertex, yVertex, zVertex, &
                         latVertex, &
                         verticesOnCell, cellsOnCell)
- 
-  nCompsDisp = 2
-  allocate(dispIncr(nCompsDisp,nVertices))
+
   !createMPs
-  numMPs = nCells
+  numMPs = 0
+  do i = 1, nCells
+    numMPs = numMPs + nEdgesOnCell(i) * mpsScaleFactorPerVtx
+  end do
+
+  print *, "Scale Factor", mpsScaleFactorPerVtx
+  print *, "NUM MPs", numMPs
+
+  allocate(lonCell(nCells))
   allocate(mpsPerElm(nCells))
   allocate(mp2Elm(numMPs))
   allocate(isMPActive(numMPs))
@@ -80,92 +148,78 @@ program main
   allocate(mpLatLon(2,numMPs))
 
   isMPActive = MP_ACTIVE !all active MPs and some changed below
-  mpsPerElm = 1 !all elements have 1 MP and some changed below
-  do i = 1,numMPs
-    mp2Elm(i) = i
+
+  numMPsCount = 0
+  do i = 1, nCells
+    localNumMPs = nEdgesOnCell(i) * mpsScaleFactorPerVtx
+    mp2Elm(numMPsCount+1:numMPsCount+localNumMPs) = i
+    mpsPerElm(i) = localNumMPs
+    numMPsCount = numMPsCount + localNumMPs
   end do
-  do i = 1, numMPs
-    inBound = .true.
+
+  call assert(numMPsCount == numMPs, "num mps miscounted")
+
+  numMPsCount = 0
+  do i = 1, nCells
+    xc = 0.0_MPAS_RKIND
+    yc = 0.0_MPAS_RKIND
+    zc = 0.0_MPAS_RKIND
     do k = 1, nEdgesOnCell(i)
       j = verticesOnCell(k,i)
-      if ((latVertex(j) .gt. 0.4*pi) .or. (latVertex(j) .lt. -0.4*pi)) then
-        inBound = .false.
-        isMPActive(i) = MP_INACTIVE
-        mpsPerElm(i) = 0
-        mp2Elm(i) = INVALID_ELM_ID
-        EXIT  
-      endif
+      xc = xc + xVertex(j) 
+      yc = yc + yVertex(j) 
+      zc = zc + zVertex(j) 
     end do
+    xc = xc/nEdgesOnCell(i)
+    yc = yc/nEdgesOnCell(i)
+    zc = zc/nEdgesOnCell(i)
 
-    if (inBound) then
-      xc = 0.0_MPAS_RKIND
-      yc = 0.0_MPAS_RKIND
-      zc = 0.0_MPAS_RKIND
-      do k = 1, nEdgesOnCell(i)
-        j = verticesOnCell(k,i)
-        xc = xc + xVertex(j) 
-        yc = yc + yVertex(j) 
-        zc = zc + zVertex(j) 
-        xComputed = sphereRadius*cos(latVertex(j))*cos(lonVertex(j))
-        yComputed = sphereRadius*cos(latVertex(j))*sin(lonVertex(j))
-        zComputed = sphereRadius*sin(latVertex(j))
-        latComputed = asin(zVertex(j)/sphereRadius)
-        lonComputed = atan2(yVertex(j),xVertex(j))
-        if (lonComputed .le. 0.0_MPAS_RKIND) then ! lon[0,2pi]
-          lonComputed = lonComputed + 2.0_MPAS_RKIND*pi
-        endif
+    lonCell(i) = atan2(yc,xc)
+    if (lonCell(i) .le. 0.0_MPAS_RKIND) then ! lon[0,2pi]
+      lonCell(i) = lonCell(i) + 2.0_MPAS_RKIND*pi
+    endif 
 
+    do k = 1, nEdgesOnCell(i)
+      j = verticesOnCell(k,i)
+      
+      ! note: m=0 not included but should lead to x_mp=xc and m=mpsScaleFactorPerVtx+1 is also not include but should lead to x_mp=xVertex(j)
+      ! so'mpsScaleFactorPerVtx+1' segments and  'mpsScaleFactorPerVtx+2' points along line from 'xc' to 'xVertex(j)'
+      ! taking only "inner" or interior 'mpsScaleFactorPerVtx' points (i.e., exclude end points of 'xc' and 'xVertex(j)') and same applies to y- and z-coordinates
+      do m = 1, mpsScaleFactorPerVtx
+        numMPsCount = numMPsCount + 1
+        xMP = (mpsScaleFactorPerVtx+1 - m) * xc + m*xVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation
+        yMP = (mpsScaleFactorPerVtx+1 - m) * yc + m*yVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation
+        zMP = (mpsScaleFactorPerVtx+1 - m) * zc + m*zVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation
+        
+        ! normalize to project each MP to be on sphere of radius 'sphereRadius'
+        radius = sqrt(xMP*xMP + yMP*yMP + zMP*zMP) ! assuming sphere center to be at origin
+        xMP = xMP/radius * sphereRadius
+        yMP = yMP/radius * sphereRadius
+        zMP = zMP/radius * sphereRadius
+        mpPosition(1,numMPsCount) = xMP
+        mpPosition(2,numMPsCount) = yMP
+        mpPosition(3,numMPsCount) = zMP
+        mpLatLon(1,numMPsCount) = asin(zMP/sphereRadius)
+        lon = atan2(yMP,xMP)
+        if (lon .le. 0.0_MPAS_RKIND) then ! lon[0,2pi]
+          lon = lon + 2.0_MPAS_RKIND*pi
+        endif 
+        mpLatLon(2,numMPsCount) = lon
       end do
-      xc = xc/nEdgesOnCell(i)
-      yc = yc/nEdgesOnCell(i)
-      zc = zc/nEdgesOnCell(i)
-      ! normalize
-      radius = sqrt(xc*xc + yc*yc + zc*zc)! assuming sphere center to be at origin
-      xc = xc/radius * sphereRadius
-      yc = yc/radius * sphereRadius
-      zc = zc/radius * sphereRadius
-      mpPosition(1,i) = xc
-      mpPosition(2,i) = yc
-      mpPosition(3,i) = zc
-      mpLatLon(1,i) = asin(zc/sphereRadius)
-      lon = atan2(yc,xc)
-      if (lon .le. 0.0_MPAS_RKIND) then ! lon[0,2pi]
-        lon = lon + 2.0_MPAS_RKIND*pi
-      endif 
-      mpLatLon(2,i) = lon
-    endif
-  end do
-  ! check first element/cell for delta
-  maxlon = minval(lonVertex)
-  minlon = maxval(lonVertex)
-  do i = 1, nEdgesOnCell(1)
-    j = verticesOnCell(i,1)
-    if(maxlon .lt. lonVertex(j)) then
-      maxlon = lonVertex(j)
-    endif
-    if(minlon .gt. lonVertex(j)) then
-      minlon = lonVertex(j)
-    endif
+    end do
   end do
+
+  call assert(numMPsCount == numMPs, "num mps miscounted")
+
   call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
   call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon))
   call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition))
 
-  
-  deltaLon = maxlon - minlon
+  call setProcWedges(mpMesh, nCells, comm_size, nEdgesOnCell, verticesOnCell, lonCell)
+  call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
+
+  call polympo_summarizeTime();
 
-  do i = 1,nVertices
-    dispIncr(1,i) = sphereRadius*cos(latVertex(i))*deltaLon
-    dispIncr(2,i) = 0.0_MPAS_RKIND
-  end do
-  call polympo_setMeshOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr))
-  call polympo_push(mpMesh)
-  do i = 1,nVertices
-    dispIncr(1,i) = sphereRadius*cos(latVertex(i))*2*deltaLon
-    dispIncr(2,i) = 0.0_MPAS_RKIND
-  end do
-  call polympo_setMeshOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr))
-  call polympo_push(mpMesh)
   call polympo_deleteMPMesh(mpMesh)
   call polympo_finalize()
 
@@ -179,7 +233,6 @@ program main
   deallocate(lonVertex)
   deallocate(verticesOnCell)
   deallocate(cellsOnCell)
-  deallocate(dispIncr)
   deallocate(mpsPerElm)
   deallocate(mp2Elm)
   deallocate(isMPActive)
@@ -187,4 +240,5 @@ program main
   deallocate(mpLatLon)
 
   stop
+
 end program
diff --git a/test/testFortranMPMeshModule.f90 b/test/testFortranMPMeshModule.f90
index 6c90fdf0..ee02098b 100644
--- a/test/testFortranMPMeshModule.f90
+++ b/test/testFortranMPMeshModule.f90
@@ -76,7 +76,7 @@ program main
   Mesharray = value1
   call polympo_setMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
   Mesharray = value2
-  call polympo_setMeshOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
+  call polympo_setMeshVtxOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
 
   Mesharray = 1
   call polympo_getMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
diff --git a/test/testMigration.cpp b/test/testMigration.cpp
new file mode 100644
index 00000000..cd12ad71
--- /dev/null
+++ b/test/testMigration.cpp
@@ -0,0 +1,84 @@
+#include "pmpo_MPMesh.hpp"
+#include "testUtils.hpp"
+
+using namespace polyMPO;
+
+IntView procWedges(Mesh* mesh, int numElms, int comm_size) {
+    auto elm2VtxConn = mesh->getElm2VtxConn();
+    auto vtxLat = mesh->getMeshField<polyMPO::MeshF_VtxRotLat>(); 
+    
+    DoubleView elmLat("elmLat",numElms);
+    DoubleView max("max", 1);
+    DoubleView min("min", 1);
+    Kokkos::parallel_for("calcElementLat", numElms, KOKKOS_LAMBDA(const int elm){  
+        int numVtx = elm2VtxConn(elm,0);
+        double sum = 0.0;
+        for(int i=1; i<= numVtx; i++){
+            sum += vtxLat(elm2VtxConn(elm,i)-1);
+        }
+        elmLat(elm) = sum/numVtx;
+        Kokkos::atomic_max(&max(0), elmLat(elm));
+        Kokkos::atomic_min(&min(0), elmLat(elm));
+    });
+
+    IntView owningProc("owningProc", numElms);
+    Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){
+        double normalizedLat = (elmLat(elm) - min(0)) / (max(0) - min(0)) * .9;
+        owningProc(elm) = normalizedLat * comm_size;
+    });
+    return owningProc;
+}
+
+int main(int argc, char* argv[] ) {
+    MPI_Init(&argc, &argv);
+    Kokkos::initialize(argc,argv);
+
+    int comm_rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank);
+    int comm_size;
+    MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
+
+    {
+        MPMesh* mpmesh = NULL;
+        void* p;
+        setWithMPASMeshByFortran(&p, argv[1], (int)strlen(argv[1]));
+        mpmesh = (MPMesh*)p;
+        Mesh* mesh = mpmesh->p_mesh;
+        int numElms = mesh->getNumElements();
+        int numMPs = numElms;
+
+        IntView mpsPerElm("mpsPerElm", numElms);
+        IntView activeMP2Elm("activeMP2Elm", numMPs);
+        IntView activeMPIDs("activeMPIDs", numMPs);
+
+        Kokkos::parallel_for("setMPsPerElm", numElms, KOKKOS_LAMBDA(const int elm){
+            mpsPerElm(elm) = 1;
+        });
+        Kokkos::parallel_for("setMPs", numMPs, KOKKOS_LAMBDA(const int mp){
+            activeMP2Elm(mp) = mp;
+            activeMPIDs(mp) = mp;
+        });
+
+        MaterialPoints p_MPs(numElms, numMPs, mpsPerElm, activeMP2Elm, activeMPIDs);
+        mpmesh->p_MPs = &p_MPs;
+
+        IntView owningProc = procWedges(mesh, numElms, comm_size);
+        mesh->setMeshEdit(true);
+        mesh->setOwningProc(owningProc);
+        for (int i=0; i<5; i++)
+            mpmesh->push();
+
+        auto MPs2Proc = p_MPs.getData<MPF_Tgt_Proc_ID>();
+        auto checkPostBackMigrate = PS_LAMBDA(const int& e, const int& mp, const bool& mask) {
+            if (mask) {
+                assert(MPs2Proc(mp) == comm_rank);
+                assert(owningProc(e) == comm_rank);
+            }
+        };
+        p_MPs.parallel_for(checkPostBackMigrate, "checkPostBackMigrate");
+    } 
+
+    Kokkos::finalize();
+    MPI_Finalize();
+    return 0;
+}
\ No newline at end of file