diff --git a/CMakeLists.txt b/CMakeLists.txt
index 080a4186..3f627daa 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -50,4 +50,4 @@ if(IS_TESTING)
   add_subdirectory (test)
 endif()
 
-bob_end_package()
\ No newline at end of file
+bob_end_package()
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 76059965..5eee8a31 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -21,6 +21,25 @@ set(SOURCES
 )
 
 add_library(polyMPO-core ${SOURCES})
+
+if(DEFINED ENV{PE_MPICH_GTL_DIR_nvidia80}
+   AND DEFINED ENV{PE_MPICH_GTL_LIBS_nvidia80})
+
+  target_link_options(polyMPO-core PUBLIC
+    "$ENV{PE_MPICH_GTL_DIR_nvidia80}"
+    "$ENV{PE_MPICH_GTL_LIBS_nvidia80}"
+  )
+
+  message(STATUS
+    "polyMPO: enabling Cray MPICH CUDA GTL linkage")
+else()
+  message(FATAL_ERROR
+    "Cray MPICH CUDA GTL environment variables are unavailable. "
+    "Load cudatoolkit and craype-accel-nvidia80, or set "
+    "CRAY_ACCEL_TARGET=nvidia80 before configuring.")
+endif()
+
+target_compile_definitions(polyMPO-core PUBLIC CUDA_AWARE_MPI)
 set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD "17")
 set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD_REQUIRED ON)
 set_property(TARGET polyMPO-core PROPERTY CXX_EXTENSIONS OFF)
diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp
index 3801dcd9..10aa52aa 100644
--- a/src/pmpo_MPMesh.cpp
+++ b/src/pmpo_MPMesh.cpp
@@ -14,11 +14,13 @@ void MPMesh::calculateStrain(){
   auto MPsBasisGrads = p_MPs->getData<MPF_Basis_Grad_Vals>();
   auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
   auto MPsStrainRate = p_MPs->getData<MPF_Strain_Rate>();
+  auto MPsArea       = p_MPs->getData<polyMPO::MPF_Area>();
   //Mesh Fields
   auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField<MeshF_TanLatVertexRotatedOverRadius>();
   auto elm2VtxConn = p_mesh->getElm2VtxConn();
   auto velField = p_mesh->getMeshField<MeshF_Vel>();
   auto solveStress = p_mesh->getMeshField<polyMPO::MeshF_SolveStress>();
+  auto elasticTimeStep = p_mesh->getElasticTimeStep();
 
   auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
     if(mask){
@@ -38,7 +40,7 @@ void MPMesh::calculateStrain(){
       double v22 = 0.0;
       double uTanOverR = 0.0;
       double vTanOverR = 0.0;
- 
+
       for (int i = 0; i < numVtx; i++){
         int iVertex = elm2VtxConn(elm, i+1)-1;
         v11 = v11 + MPsBasisGrads(mp, i*2 + 0) * velField(iVertex, 0);
@@ -52,6 +54,8 @@ void MPMesh::calculateStrain(){
       MPsStrainRate(mp, 0) =  v11 - vTanOverR;
       MPsStrainRate(mp, 1) =  v22;
       MPsStrainRate(mp, 2) =  0.5*(v12 + v21 + uTanOverR);
+
+      MPsArea(mp, 0) = MPsArea(mp, 0) * exp((v11+v22-vTanOverR)*elasticTimeStep);
     }
   };
   p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate");
@@ -75,13 +79,16 @@ void MPMesh::calculateStress(const int constitutive_relation){
     if(mask){
       Vec3d strain_rate (MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
       Vec3d stress(MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
+      double rep_pressure=MPsRepPressure(mp,0);
 
       if (constitutive_relation == 1)
-        constitutive_evp(strain_rate, stress, MPsIcePressure(mp,0), MPsRepPressure(mp,0), MPsArea(mp,0), elasticTimeStep, dampingTimescale);
+        constitutive_evp(strain_rate, stress, MPsIcePressure(mp,0), rep_pressure, MPsArea(mp,0), elasticTimeStep, dampingTimescale);
       else if(constitutive_relation == 3)
         constitutive_linear(strain_rate, stress);
+
       for (int m=0 ; m<3; m++)
         MPsStress(mp, m) = stress[m]*solveStress(elm);
+      MPsRepPressure(mp,0)=rep_pressure;
     }
   };
   p_MPs->parallel_for(setMPStress, "setMPStress");
@@ -89,11 +96,21 @@ void MPMesh::calculateStress(const int constitutive_relation){
 
 void MPMesh::calculateStressDivergence(){
 
-  Kokkos::Timer timer;
   int self, numProcsTot;
   MPI_Comm comm = p_MPs->getMPIComm();
   MPI_Comm_rank(comm, &self);
   MPI_Comm_size(comm, &numProcsTot);
+  
+  //Kokkos::Timer b0Timer;
+  //Kokkos::fence();
+
+  //const double b0DrainTime = b0Timer.seconds();                        //B0: drain any device work left from the previous phase
+  //MPI_Barrier(comm);
+  //const double b0Sync = b0Timer.seconds();                             //B0: align all ranks so this timed region starts at the same instant
+  //const double b0BarrierWait = b0Sync - b0DrainTime; 
+
+  Kokkos::Timer totalTimer;
+  Kokkos::Timer computeTimer;
 
   //Mesh Information
   auto elm2VtxConn = p_mesh->getElm2VtxConn();
@@ -145,7 +162,7 @@ void MPMesh::calculateStressDivergence(){
                                                                 (1.0 - ramp) * invM * w_vtx;
 
         factor = factor * tanLatVertexRotatedOverRadius(vID, 0);
-      
+
         auto factor1 = ramp * (w_vtx/radius) * (VtxCoeffs_new(vID, 1, 0) + VtxCoeffs_new(vID, 1, 1)*CoordDiffs[1]  +
                                                                            VtxCoeffs_new(vID, 1, 2)*CoordDiffs[2]  +
                                                                            VtxCoeffs_new(vID, 1, 3)*CoordDiffs[3]) -
@@ -165,18 +182,44 @@ void MPMesh::calculateStressDivergence(){
     }
   };
   p_MPs->parallel_for(stress_div, " stress_div_assembly");
-  Kokkos::fence();
-  pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); 
+  Kokkos::fence();                                       //drain device work -> compute really is done on this rank
+
+  const double computeTime = computeTimer.seconds();     //T_before: pure local compute time, no waiting
+
+  //MPI_Barrier(comm);                                      //B1: fast ranks wait here for the slowest rank
+
+  //const double computeTimeSync = computeTimer.seconds();  //time until every rank reached the barrier
+  //const double b1BarrierWait   = computeTimeSync - computeTime; //this rank's wait time = compute load imbalance
+
+  Kokkos::Timer communicationTimer;
 
-  timer.reset();
   if(numProcsTot>1){ 
     //Takes contribution of halo vertices and adds it in owner procs
-    communicate_and_take_halo_contributions1(stress_divUV, numVertices, 2, 0, 0);
+    communicate_and_take_halo_contributions1_improved(stress_divUV, numVertices, 2, 0, 0, "Stress_Divergence");
     //Transfer the correct values at owned vertices to halo vertices
     //communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 1, 1);
   }
-  Kokkos::fence();
-  pumipic::RecordTime("Stress_Divergence Communication" + std::to_string(self), timer.seconds());  
+  Kokkos::fence();                                        //drain device work from the communication step
+
+  const double communicationTime = communicationTimer.seconds(); //pure local communication time, no waiting
+
+  //MPI_Barrier(comm);                                       //B2: fast ranks wait here for the slowest rank
+
+  //const double communicationTimeSync = communicationTimer.seconds();
+  //const double b2BarrierWait         = communicationTimeSync - communicationTime; //communication load imbalance
+
+  const double totalTime = totalTimer.seconds();
+
+  //pumipic::RecordTime("Stress_Divergence_B0_Drain_" + std::to_string(self), b0DrainTime);
+  //pumipic::RecordTime("Stress_Divergence_B0_BarrierWait_" + std::to_string(self), b0BarrierWait);
+
+  pumipic::RecordTime("Stress_Divergence_Compute_" + std::to_string(self),computeTime);
+  //pumipic::RecordTime("Stress_Divergence_Compute_B1_Barrier_Wait_Time_" + std::to_string(self),b1BarrierWait);
+
+  pumipic::RecordTime("Stress_Divergence_Communication_" + std::to_string(self),communicationTime);
+  //pumipic::RecordTime("Stress_Divergence_Communication_B2_Wait_Time_" + std::to_string(self), b2BarrierWait);
+
+  pumipic::RecordTime("Stress_Divergence_Total_" + std::to_string(self),totalTime);
 }
 
 void MPMesh::calcBasis() {
@@ -280,7 +323,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
       Vec3d dx = MPnew-MP;
       while(true){
         int numConnElms = elm2ElmConn(iElm,0);
-                
+
         Vec3d center(elmCenter(iElm, 0), elmCenter(iElm, 1), elmCenter(iElm, 2));
         Vec3d delta = MPnew - center;
 
@@ -441,7 +484,13 @@ void MPMesh::startCommunication(){
   int self, numProcsTot;
   MPI_Comm comm = p_MPs->getMPIComm();
   MPI_Comm_rank(comm, &self);
-  MPI_Comm_size(comm, &numProcsTot); 
+  MPI_Comm_size(comm, &numProcsTot);
+
+  std::cout << "[RankSummary] Rank=" << self
+            << " Vertices(total)=" << p_mesh->getNumVertices()
+            << " Vertices(owned)=" << p_mesh->getNumVerticesOwned()
+            << " Elements=" << p_mesh->getNumElements()
+            << std::endl; 
 
   //The routine should work for elements too, although currently the communication 
   //is done for vertices. For elements, the follwoing three variables should correspond 
@@ -742,8 +791,8 @@ void MPMesh::T2LTracking(Vec2dView dx){
         Vec2d MP(mpPositions(mp,0),mpPositions(mp,1));//XXX:the input is XYZ, but we only support 2d vector
         if(mask){
             int iElm = elm;
-            Vec2d MPnew = MP + dx(mp);    
-            
+            Vec2d MPnew = MP + dx(mp);
+
             while(true){
                 int numVtx = elm2VtxConn(iElm,0);
                 bool goToNeighbour = false;
@@ -753,7 +802,7 @@ void MPMesh::T2LTracking(Vec2dView dx){
                     v[i] = elm2VtxConn(iElm,i+1)-1;
                 //get edges and perpendiculardx
                 Vec2d e[maxVtxsPerElm];
-                double pdx[maxVtxsPerElm];                    
+                double pdx[maxVtxsPerElm];
                 for(int i=0; i< numVtx; i++){
                     int idx_ip1 = (i+1)%numVtx;
                     Vec2d v_i(vtxCoords(v[i],0),vtxCoords(v[i],1));
@@ -761,17 +810,17 @@ void MPMesh::T2LTracking(Vec2dView dx){
                     e[i] = v_ip1 - v_i;
                     pdx[i] = (v_i - MP).cross(dx(mp));
                 }
-                
+
                 for(int i=0; i<numVtx; i++){
                     int ip1 = (i+1)%numVtx;
-                    //pdx*pdx<0 and edge is acrossed 
+                    //pdx*pdx<0 and edge is acrossed
                     if(pdx[i]*pdx[ip1] <0 && e[i].cross(Vec2d(MPnew[0]-vtxCoords(v[i],0),
                                                               MPnew[1]-vtxCoords(v[i],1)))<0){
                         //go to the next elm
                         iElm = elm2ElmConn(iElm,i+1);
                         goToNeighbour = true;
                         if(iElm <0){
-                            mpStatus(mp) = 0;                  
+                            mpStatus(mp) = 0;
                             MPs2Elm(mp) = -1;
                             goToNeighbour = false;
                         }
diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp
index d3aac1b7..3ddb1032 100644
--- a/src/pmpo_MPMesh.hpp
+++ b/src/pmpo_MPMesh.hpp
@@ -4,6 +4,14 @@
 #include "pmpo_utils.hpp"
 #include "pmpo_mesh.hpp"
 #include "pmpo_materialPoints.hpp"
+#include <cstdlib>
+#include <iostream>
+#include <utility>
+#include <stdexcept>
+#include <string>
+#ifdef KOKKOS_ENABLE_CUDA
+#include <cuda_runtime.h>
+#endif
 
 namespace polyMPO{
 
@@ -28,22 +36,27 @@ class MPMesh{
     int numOwnersTot, numHalosTot;
     std::vector<int> numOwnersOnOtherProcs;
     std::vector<int> numHalosOnOtherProcs;
-    std::vector<int>haloOwnerProcs;
+    std::vector<int> haloOwnerProcs;
     std::vector<std::vector<int>> haloOwnerLocalIDs;
     std::vector<std::vector<int>> ownerOwnerLocalIDs;
     std::vector<std::vector<int>> ownerHaloLocalIDs;
 
     void startCommunication();
     
-    void communicate_and_take_halo_contributions(const Kokkos::View<double**>& meshField, int nEntities, int numEntries, int mode, int op);
+    void communicate_and_take_halo_contributions(
+        const Kokkos::View<double**>& meshField,
+        int nEntities,
+        int numEntries,
+        int mode,
+        int op);
     
-    //Now Kokkos views are made 1D
+    // Original CPU-staging function
     template <typename ViewType>
     void communicate_and_take_halo_contributions1(
         const ViewType& meshField,
         int nEntities,
         int numEntries,
-        int mode ,
+        int mode,
         int op){
 
       int self;
@@ -51,95 +64,155 @@ class MPMesh{
       MPI_Comm_rank(comm, &self);
 
       Kokkos::Timer timer;
-      auto reconVals_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), meshField);
+      auto reconVals_host =
+          Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), meshField);
+
       pumipic::RecordTime("SD: GPU-CPU copy-" + std::to_string(self), timer.seconds());
 
       timer.reset();
-      std::vector<std::vector<int>>    recvIDVec;
+      std::vector<std::vector<int>> recvIDVec;
       std::vector<std::vector<double>> recvDataVec;
+
       pumipic::RecordTime("SD: Recv Vec Allocation-" + std::to_string(self), timer.seconds());
       
       timer.reset();
-      //communicateFields1(fieldData1, nEntities, numEntries, mode, recvIDVec, recvDataVec);
-      communicateFields1(reconVals_host, nEntities, numEntries, mode, recvIDVec, recvDataVec);
+
+      communicateFields1(
+          reconVals_host,
+          nEntities,
+          numEntries,
+          mode,
+          recvIDVec,
+          recvDataVec);
+
       pumipic::RecordTime("SD: IP Comm-" + std::to_string(self), timer.seconds());
 
       timer.reset();
+
       int numProcsTot = recvIDVec.size();
-      //Flatten IDs
+
       int totalSize = 0;
-      std::vector<int> offsets(numProcsTot, 0); 
-      for(int i=0; i<numProcsTot; i++) {
+      std::vector<int> offsets(numProcsTot, 0);
+
+      for(int i = 0; i < numProcsTot; i++){
         offsets[i] = totalSize;
         totalSize += recvIDVec[i].size();
       }
-      std::vector<int> flatIDVec(totalSize, 0);
-      for(int i=0; i<numProcsTot; i++){
-        std::copy(recvIDVec[i].begin(), recvIDVec[i].end(), flatIDVec.begin() + offsets[i]);
+
+      Kokkos::View<int*> recvIDGPU("recvIDGPU", totalSize);
+      auto hostView =
+          Kokkos::View<int*, Kokkos::HostSpace>("recvIDCPU", totalSize);
+
+      for(int i = 0; i < numProcsTot; i++){
+        std::copy(
+            recvIDVec[i].begin(),
+            recvIDVec[i].end(),
+            hostView.data() + offsets[i]);
       }
+
       pumipic::RecordTime("SD: Flatten IDs-" + std::to_string(self), timer.seconds());
 
       timer.reset();
-      Kokkos::View<int*> recvIDGPU("recvIDGPU", totalSize);
-      auto hostView = Kokkos::View<int*, Kokkos::HostSpace>("recvIDCPU", totalSize);
-      std::copy(flatIDVec.begin(), flatIDVec.end(), hostView.data());
+
       Kokkos::deep_copy(recvIDGPU, hostView);
       Kokkos::fence();
+
       pumipic::RecordTime("SD: Copy CPU-GPU-" + std::to_string(self), timer.seconds());
 
-      //Flatten Data
       timer.reset();
-      int totalSize_data=0;
+
+      int totalSize_data = 0;
       std::vector<int> offsets_data(numProcsTot, 0);
-      for(int i=0; i<numProcsTot; i++){
+
+      for(int i = 0; i < numProcsTot; i++){
         offsets_data[i] = totalSize_data;
         totalSize_data += recvDataVec[i].size();
       }
-      std::vector<double> flatDataVec(totalSize_data, 0);
-      for(int i=0; i<numProcsTot; i++) {
-        std::copy(recvDataVec[i].begin(), recvDataVec[i].end(), flatDataVec.begin() + offsets_data[i]);
+
+      Kokkos::View<double*> recvDataGPU("recvDataGPU", totalSize_data);
+      auto hostView_data =
+          Kokkos::View<double*, Kokkos::HostSpace>("recvDataCPU", totalSize_data);
+
+      for(int i = 0; i < numProcsTot; i++){
+        std::copy(
+            recvDataVec[i].begin(),
+            recvDataVec[i].end(),
+            hostView_data.data() + offsets_data[i]);
       }
+
       pumipic::RecordTime("SD: Flatten Data-" + std::to_string(self), timer.seconds());
 
       timer.reset();
-      Kokkos::View<double*> recvDataGPU("recvDataGPU", totalSize_data);
-      auto hostView_data= Kokkos::View<double*, Kokkos::HostSpace>("recvDataCPU", totalSize_data);
-      std::copy(flatDataVec.begin(), flatDataVec.end(), hostView_data.data()); 
+
       Kokkos::deep_copy(recvDataGPU, hostView_data);
       Kokkos::fence();
-      assert(totalSize_data == totalSize*numEntries);
-      for (int i=0; i<numProcsTot; i++){
+
+      assert(totalSize_data == totalSize * numEntries);
+
+      for(int i = 0; i < numProcsTot; i++){
         assert(recvDataVec[i].size() == recvIDVec[i].size() * numEntries);
       }
+
       pumipic::RecordTime("SD: Copy CPU-GPU2-" + std::to_string(self), timer.seconds());
  
-      //Take contributions from other procs
       timer.reset();
-      Kokkos::parallel_for("halo contribution", recvIDGPU.size(), KOKKOS_LAMBDA(const int i){
-        int vertex = recvIDGPU(i);
-        for(int k=0; k<numEntries; k++){
-          if(op==0) Kokkos::atomic_add(&meshField(vertex,k), recvDataGPU(i*numEntries+k));
-          if(op==1) meshField(vertex, k) = recvDataGPU(i * numEntries + k);
-        }
-      });
+
+      if(op == 0){
+        Kokkos::parallel_for(
+            "halo contribution add",
+            recvIDGPU.size(),
+            KOKKOS_LAMBDA(const int i){
+              const int vertex = recvIDGPU(i);
+
+              for(int k = 0; k < numEntries; k++){
+                Kokkos::atomic_add(
+                    &meshField(vertex,k),
+                    recvDataGPU(i * numEntries + k));
+              }
+            });
+      }
+      else{
+        Kokkos::parallel_for(
+            "halo contribution assign",
+            recvIDGPU.size(),
+            KOKKOS_LAMBDA(const int i){
+              const int vertex = recvIDGPU(i);
+
+              for(int k = 0; k < numEntries; k++){
+                meshField(vertex,k) =
+                    recvDataGPU(i * numEntries + k);
+              }
+            });
+      }
+
       Kokkos::fence();
+
       pumipic::RecordTime("SD: Contribution" + std::to_string(self), timer.seconds());
     }
 
 
-    void communicateFields(const std::vector<std::vector<double>>& fieldData, const int numEntities, const int numEntries, int mode,
-                              std::vector<std::vector<int>>& recvIDVec,  std::vector<std::vector<double>>& recvDataVec);
+    void communicateFields(
+        const std::vector<std::vector<double>>& fieldData,
+        const int numEntities,
+        const int numEntries,
+        int mode,
+        std::vector<std::vector<int>>& recvIDVec,
+        std::vector<std::vector<double>>& recvDataVec);
    
 
     template <class ViewType> 
     void communicateFields1(
-        const ViewType& fieldData, 
-        const int numEntities, const int numEntries, int mode,
+        const ViewType& fieldData,
+        const int numEntities,
+        const int numEntries,
+        int mode,
         std::vector<std::vector<int>>& recvIDVec,
         std::vector<std::vector<double>>& recvDataVec){
     
       int self, numProcsTot;
+
       MPI_Comm comm = p_MPs->getMPIComm();
+
       MPI_Comm_rank(comm, &self);
       MPI_Comm_size(comm, &numProcsTot);
 
@@ -151,98 +224,191 @@ class MPMesh{
       recvDataVec.resize(numProcsTot);
 
       for(int i = 0; i < numProcsTot; i++){
-        if(i==self) continue;
+        if(i == self) continue;
 
-        int numToSend = 0, numToRecv = 0; 
-        if(mode == 0) {
-          //gather (halos send to owners)
+        int numToSend = 0;
+        int numToRecv = 0;
+
+        if(mode == 0){
           numToSend = numOwnersOnOtherProcs[i];
           numToRecv = numHalosOnOtherProcs[i];
         }
-        else{ 
-          //scatter (owners send to halos)
+        else{
           numToSend = numHalosOnOtherProcs[i];
           numToRecv = numOwnersOnOtherProcs[i];
         }
 
         if(numToSend > 0){
-          sendDataVec[i].reserve(numToSend*numEntries);
+          sendDataVec[i].reserve(numToSend * numEntries);
         }
+
         if(numToRecv > 0){
-          recvDataVec[i].resize(numToRecv*numEntries);
+          recvDataVec[i].resize(numToRecv * numEntries);
           recvIDVec[i].resize(numToRecv);
         }
       }
 
       if(mode == 0){
-        // Halos sends to owners
-        for (int iEnt = 0; iEnt < numHalosTot; iEnt++){
+        for(int iEnt = 0; iEnt < numHalosTot; iEnt++){
           auto ownerProc = haloOwnerProcs[iEnt];
-          for (int iDouble = 0; iDouble < numEntries; iDouble++)
-            sendDataVec[ownerProc].push_back(fieldData(numOwnersTot+iEnt, iDouble));
+
+          for(int iDouble = 0; iDouble < numEntries; iDouble++){
+            sendDataVec[ownerProc].push_back(
+                fieldData(numOwnersTot + iEnt, iDouble));
+          }
         }
       }
       else if(mode == 1){
-        // Owner sends to halos
-        for (size_t iProc=0; iProc<ownerOwnerLocalIDs.size(); iProc++) {
-          for (auto& ownerID : ownerOwnerLocalIDs[iProc]) {
-            for (int iDouble = 0; iDouble < numEntries; iDouble++)
-              sendDataVec[iProc].push_back(fieldData(ownerID, iDouble));
+        for(size_t iProc = 0; iProc < ownerOwnerLocalIDs.size(); iProc++){
+          for(auto& ownerID : ownerOwnerLocalIDs[iProc]){
+            for(int iDouble = 0; iDouble < numEntries; iDouble++){
+              sendDataVec[iProc].push_back(
+                  fieldData(ownerID, iDouble));
+            }
           }
         }
       }
 
       std::vector<MPI_Request> requests;
-      requests.reserve(4*numProcsTot);
+      requests.reserve(4 * numProcsTot);
+
       for(int proc = 0; proc < numProcsTot; proc++){
-        if(proc == self) continue;  
+        if(proc == self) continue;
+
         if(mode == 0 && numHalosOnOtherProcs[proc]){
-          assert(recvIDVec[proc].size() == (size_t)numHalosOnOtherProcs[proc]);
-          assert(recvDataVec[proc].size() == recvIDVec[proc].size() * (size_t)numEntries);
-          MPI_Request req3, req4;
-          MPI_Irecv(recvIDVec[proc].data(), recvIDVec[proc].size(), MPI_INT, proc, 1, comm, &req3);
-          MPI_Irecv(recvDataVec[proc].data(), recvDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req4);
+          assert(recvIDVec[proc].size() ==
+                 static_cast<size_t>(numHalosOnOtherProcs[proc]));
+
+          assert(recvDataVec[proc].size() ==
+                 recvIDVec[proc].size() * static_cast<size_t>(numEntries));
+
+          MPI_Request req3;
+          MPI_Request req4;
+
+          MPI_Irecv(
+              recvIDVec[proc].data(),
+              recvIDVec[proc].size(),
+              MPI_INT,
+              proc,
+              1,
+              comm,
+              &req3);
+
+          MPI_Irecv(
+              recvDataVec[proc].data(),
+              recvDataVec[proc].size(),
+              MPI_DOUBLE,
+              proc,
+              2,
+              comm,
+              &req4);
+
           requests.push_back(req3);
           requests.push_back(req4);
         }
-        if(mode == 0 && numOwnersOnOtherProcs[proc]) {
-          assert(haloOwnerLocalIDs[proc].size() == (size_t)numOwnersOnOtherProcs[proc]);
-          assert(sendDataVec[proc].size() == haloOwnerLocalIDs[proc].size() * (size_t)numEntries);
-          MPI_Request req1, req2;
-          MPI_Isend(haloOwnerLocalIDs[proc].data(), haloOwnerLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req1);
-          MPI_Isend(sendDataVec[proc].data(), sendDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req2);
+
+        if(mode == 0 && numOwnersOnOtherProcs[proc]){
+          assert(haloOwnerLocalIDs[proc].size() ==
+                 static_cast<size_t>(numOwnersOnOtherProcs[proc]));
+
+          assert(sendDataVec[proc].size() ==
+                 haloOwnerLocalIDs[proc].size() * static_cast<size_t>(numEntries));
+
+          MPI_Request req1;
+          MPI_Request req2;
+
+          MPI_Isend(
+              haloOwnerLocalIDs[proc].data(),
+              haloOwnerLocalIDs[proc].size(),
+              MPI_INT,
+              proc,
+              1,
+              comm,
+              &req1);
+
+          MPI_Isend(
+              sendDataVec[proc].data(),
+              sendDataVec[proc].size(),
+              MPI_DOUBLE,
+              proc,
+              2,
+              comm,
+              &req2);
+
           requests.push_back(req1);
           requests.push_back(req2);
         }
 
         if(mode == 1 && numOwnersOnOtherProcs[proc]){
-          MPI_Request req3, req4;
-          MPI_Irecv(recvIDVec[proc].data(), recvIDVec[proc].size(), MPI_INT, proc, 1, comm, &req3);
-          MPI_Irecv(recvDataVec[proc].data(), recvDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req4);
+          MPI_Request req3;
+          MPI_Request req4;
+
+          MPI_Irecv(
+              recvIDVec[proc].data(),
+              recvIDVec[proc].size(),
+              MPI_INT,
+              proc,
+              1,
+              comm,
+              &req3);
+
+          MPI_Irecv(
+              recvDataVec[proc].data(),
+              recvDataVec[proc].size(),
+              MPI_DOUBLE,
+              proc,
+              2,
+              comm,
+              &req4);
+
           requests.push_back(req3);
           requests.push_back(req4);
         }
-        if(mode == 1 && numHalosOnOtherProcs[proc]) {
-          MPI_Request req1, req2;
-          MPI_Isend(ownerHaloLocalIDs[proc].data(), ownerHaloLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req1);
-          MPI_Isend(sendDataVec[proc].data(), sendDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req2);
+
+        if(mode == 1 && numHalosOnOtherProcs[proc]){
+          MPI_Request req1;
+          MPI_Request req2;
+
+          MPI_Isend(
+              ownerHaloLocalIDs[proc].data(),
+              ownerHaloLocalIDs[proc].size(),
+              MPI_INT,
+              proc,
+              1,
+              comm,
+              &req1);
+
+          MPI_Isend(
+              sendDataVec[proc].data(),
+              sendDataVec[proc].size(),
+              MPI_DOUBLE,
+              proc,
+              2,
+              comm,
+              &req2);
+
           requests.push_back(req1);
           requests.push_back(req2);
         }
       }
+
       MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
     }
     
+
     MPMesh(Mesh* inMesh, MaterialPoints* inMPs):
-      p_mesh(inMesh), p_MPs(inMPs) {
+      p_mesh(inMesh),
+      p_MPs(inMPs) {
     };
 
-    ~MPMesh() {
+
+    ~MPMesh(){
       delete p_mesh;
       delete p_MPs;
     }
 
-    //MP advection and tracking
+
+    // MP advection and tracking
     void CVTTrackingEdgeCenterBased(Vec2dView dx);
     void CVTTrackingElmCenterBased(const int printVTPIndex = -1);
     void T2LTracking(Vec2dView dx);
@@ -252,47 +418,705 @@ class MPMesh{
     void push_swap_pos();
     void push();
 
-    //Used before advection to interpolate fields from mesh to MPs
-    //And also before reconstruction
+
+    // Used before advection to interpolate fields from mesh to MPs
+    // And also before reconstruction
     void calcBasis();
 
-    //Reconstruction
+
+    // Reconstruction
     DoubleView assemblyV0();
+
     template <MeshFieldIndex meshFieldIndex>
     void assemblyVtx0();
+
     template <MeshFieldIndex meshFieldIndex>
     void assemblyElm0();
+
     template <MeshFieldIndex meshFieldIndex>
     void assemblyVtx1();
+
     void reconstruct_coeff_full();
-    void invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius);
+
+    void invertMatrix(
+        const Kokkos::View<double**>& vtxMatrices,
+        const double& radius);
+
     Kokkos::View<double*[vec3d_nEntries][vec4d_nEntries]> precomputedVtxCoeffs_new;
-    Kokkos::View<double*> nearAnEdge;   
+    Kokkos::View<double*> nearAnEdge;
     Kokkos::View<double*> vtxMatrixMass;
 
-    //Not used currently
-    std::map<MeshFieldIndex, std::function<void()>> reconstructSlice = std::map<MeshFieldIndex, std::function<void()>>();
+
+    // Not used currently
+    std::map<MeshFieldIndex, std::function<void()>> reconstructSlice =
+        std::map<MeshFieldIndex, std::function<void()>>();
+
     template <MaterialPointSlice index>
     DoubleView wtScaAssembly();
+
     template <MaterialPointSlice index>
     Vec2dView wtVec2Assembly();
+
     template <MeshFieldIndex meshFieldIndex>
-    void assembly(int order, MeshFieldType type, bool basisWeightFlag, bool massWeightFlag);
+    void assembly(
+        int order,
+        MeshFieldType type,
+        bool basisWeightFlag,
+        bool massWeightFlag);
+
     template<MeshFieldIndex meshFieldIndex>
-    void setReconstructSlice(int order, MeshFieldType type);
+    void setReconstructSlice(
+        int order,
+        MeshFieldType type);
+
     void reconstructSlices();
 
     void printVTP_mesh(int printVTPIndex);
-    void writeMPTrackingVTP(int printVTPIndex, int numMPs, const Vec3dView& history, const Vec3dView& resultLeft,
-                            const Vec3dView& resultRight, const Vec3dView& mpTgtPosArray);
 
+    void writeMPTrackingVTP(
+        int printVTPIndex,
+        int numMPs,
+        const Vec3dView& history,
+        const Vec3dView& resultLeft,
+        const Vec3dView& resultRight,
+        const Vec3dView& mpTgtPosArray);
 
     void calculateStrain();
     void calculateStress(const int constitutive_relation);
     void calculateStressDivergence();
+
+
+
+#ifdef CUDA_AWARE_MPI
+
+    // Use explicit CUDA memory space for MPI device buffers when CUDA is
+    // enabled. This avoids ambiguity in the default Kokkos::View memory space
+    // and gives Cray MPICH/GTL plain CUDA allocations to register/export.
+#ifdef KOKKOS_ENABLE_CUDA
+    // Plain cudaMalloc'd device buffer, exposed as an unmanaged Kokkos::View
+    // via .view(). Used only for the 4 buffers below that get handed
+    // directly to MPI_Isend/Irecv under CUDA-aware MPI.
+    //
+    // Why not just a Kokkos::View<T*, Kokkos::CudaSpace>: when Kokkos is
+    // built with Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC=ON (the default since
+    // Kokkos 4.2), View allocations use cudaMallocAsync/memory pools.
+    // cuIpcGetMemHandle (which Cray MPICH/GTL uses for intra-node GPU-to-GPU
+    // sends) rejects pool allocations with CUDA_ERROR_INVALID_VALUE. Since
+    // polyMPO isn't allowed to touch the Kokkos build config, these 4
+    // buffers bypass Kokkos's allocator entirely via a direct cudaMalloc,
+    // which cuIpcGetMemHandle always accepts, regardless of how the rest of
+    // Kokkos (or the rest of the app's Views) is configured.
+    template <typename T>
+    struct RawCudaMPIBuffer{
+      T* ptr = nullptr;
+      size_t count = 0;
+
+      void allocate(size_t n){
+        free();
+        count = n;
+        if(n > 0){
+          cudaError_t err = cudaMalloc(&ptr, n * sizeof(T));
+          if(err != cudaSuccess){
+            throw std::runtime_error(
+                std::string("RawCudaMPIBuffer: cudaMalloc failed: ") +
+                cudaGetErrorString(err));
+          }
+        }
+      }
+
+      void free(){
+        if(ptr != nullptr){ cudaFree(ptr); ptr = nullptr; }
+        count = 0;
+      }
+
+      T* data() const{ return ptr; }
+      size_t size() const{ return count; }
+
+      Kokkos::View<T*, Kokkos::CudaSpace, Kokkos::MemoryUnmanaged> view() const{
+        return Kokkos::View<T*, Kokkos::CudaSpace, Kokkos::MemoryUnmanaged>(
+            ptr, count);
+      }
+
+      RawCudaMPIBuffer() = default;
+      ~RawCudaMPIBuffer(){ free(); }
+
+      RawCudaMPIBuffer(const RawCudaMPIBuffer&) = delete;
+      RawCudaMPIBuffer& operator=(const RawCudaMPIBuffer&) = delete;
+
+      RawCudaMPIBuffer(RawCudaMPIBuffer&& other) noexcept{
+        ptr = other.ptr; count = other.count;
+        other.ptr = nullptr; other.count = 0;
+      }
+
+      RawCudaMPIBuffer& operator=(RawCudaMPIBuffer&& other) noexcept{
+        if(this != &other){
+          free();
+          ptr = other.ptr; count = other.count;
+          other.ptr = nullptr; other.count = 0;
+        }
+        return *this;
+      }
+    };
+
+    using CudaAwareMPIIntBuffer = RawCudaMPIBuffer<int>;
+    using CudaAwareMPIDoubleBuffer = RawCudaMPIBuffer<double>;
+#else
+    // No CUDA backend: no CUDA IPC/pool-allocation concern, so just wrap a
+    // normal Kokkos::View with the same .allocate()/.data()/.view()
+    // interface as RawCudaMPIBuffer above, so the cache struct and its call
+    // sites below don't need to branch on KOKKOS_ENABLE_CUDA.
+    template <typename T>
+    struct KokkosMPIBuffer{
+      Kokkos::View<T*> v;
+
+      void allocate(size_t n){
+        v = Kokkos::View<T*>("cudaAwareMPIBuffer_batched", n);
+      }
+
+      T* data() const{ return v.data(); }
+      size_t size() const{ return v.extent(0); }
+      Kokkos::View<T*> view() const{ return v; }
+    };
+
+    using CudaAwareMPIIntBuffer = KokkosMPIBuffer<int>;
+    using CudaAwareMPIDoubleBuffer = KokkosMPIBuffer<double>;
+#endif
+
+    // Cached CUDA-aware MPI communication metadata and batched GPU buffers.
+    // Every neighbor's data lives in one shared allocation (see
+    // CudaAwareMPIFieldCache below) and MPI is given "base pointer + byte
+    // offset" per neighbor rather than a separate allocation per neighbor.
+    // If you ever need to fall back to one allocation per neighbor (e.g. an
+    // MPI/GPU stack that mishandles offset device pointers for CUDA IPC),
+    // restore the per-proc-buffer version from version control.
+    bool cudaAwareMPICacheValid = true;
+    bool cudaAwareMPIDisabled = false;
+    bool cudaAwareMPIEnvChecked = false;
+    bool cudaAwareMPIForceCPU = false;
+    bool cudaAwareMPILogged = false;
+
+    struct CudaAwareMPIFieldCache{
+      bool valid = false;
+      int cachedNumProcs = -1;
+
+      std::vector<int> sendCounts;
+      std::vector<int> recvCounts;
+      std::vector<int> sendOffsets; // prefix sum of sendCounts, in entities
+      std::vector<int> recvOffsets; // prefix sum of recvCounts, in entities
+
+      int totalSendCount = 0;
+      int totalRecvCount = 0;
+
+      // Single batched GPU buffers (one allocation each, instead of one
+      // Kokkos::View per neighbor proc). Per-proc slices are
+      // [offset, offset + count) for the ID buffers, and
+      // [offset * numEntries, (offset + count) * numEntries) for the data
+      // buffers. MPI is given "buffer base pointer + offset", not a
+      // separate allocation per proc.
+      CudaAwareMPIIntBuffer sendEntityGPU;
+      CudaAwareMPIIntBuffer recvIDGPU;
+
+      CudaAwareMPIDoubleBuffer sendDataGPU;
+      CudaAwareMPIDoubleBuffer recvDataGPU;
+    };
+
+    std::map<std::pair<int, int>, CudaAwareMPIFieldCache> cudaAwareMPICaches;
+
+    bool cudaAwareMPIForceDisabled(){
+      if(!cudaAwareMPIEnvChecked){
+        const char* value = std::getenv("POLYMPO_DISABLE_CUDA_AWARE_MPI");
+        cudaAwareMPIForceCPU =
+            value != nullptr && value[0] != '\0' && value[0] != '0';
+        cudaAwareMPIEnvChecked = true;
+      }
+
+      return cudaAwareMPIForceCPU;
+    }
+
+    // Fully CUDA-aware MPI version, batched buffer variant:
+    // Field data is sent/received using GPU pointers. Receive IDs are cached
+    // once from the fixed halo/owner mapping and are not sent every call.
+    //
+    // Every neighbor's send/recv entity-ID list and data live in ONE big
+    // GPU buffer each (laid out back-to-back in proc order), instead of one
+    // Kokkos::View allocation per neighbor. Packing/unpacking is a single
+    // kernel launch over all neighbors' entities at once instead of one
+    // launch per neighbor, and MPI_Isend/Irecv use "buffer base pointer +
+    // offset" into that single buffer per proc. This is what actually
+    // shrinks MPI_Wait time: fewer, larger, more uniform in-flight
+    // transfers instead of many small independent ones.
+    //
+    // Note: an earlier version of this cache used one Kokkos::View
+    // allocation per neighbor specifically to avoid handing MPI a
+    // "base pointer + offset" GPU address, out of concern for CUDA IPC
+    // issues on Cray MPICH/GTL. That failure mode was root-caused to
+    // Kokkos allocating device Views via cudaMallocAsync (invalid for
+    // cuIpcGetMemHandle), not to offset pointers themselves, and is fixed
+    // by building Kokkos with -DKokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC=OFF.
+    // If you ever do hit IPC trouble that tracks back to offset pointers
+    // specifically, the per-proc-buffer version can be restored from
+    // version control.
+    //
+    // Important:
+    // This function caches communication metadata and GPU buffers per
+    // (mode, numEntries). If the communication pattern changes, clear
+    // cudaAwareMPICaches before the next call.
+    template <typename ViewType>
+    void communicate_and_take_halo_contributions1_improved(
+        const ViewType& meshField,
+        int nEntities,
+        int numEntries,
+        int mode,
+        int op,
+        const std::string& label){
+
+      int self, numProcsTot;
+
+      MPI_Comm comm = p_MPs->getMPIComm();
+
+      MPI_Comm_rank(comm, &self);
+      MPI_Comm_size(comm, &numProcsTot);
+
+      const char* diagnosticsEnv =
+      
+      std::getenv("POLYMPO_MPI_DIAGNOSTICS");
+
+      const bool mpiDiagnostics =
+      diagnosticsEnv != nullptr &&
+      std::atoi(diagnosticsEnv) != 0;
+
+      assert(mode == 0 || mode == 1);
+      assert(op == 0 || op == 1);
+      assert(nEntities == numOwnersTot + numHalosTot);
+
+      if(cudaAwareMPIDisabled || cudaAwareMPIForceDisabled()){
+        communicate_and_take_halo_contributions1(
+            meshField,
+            nEntities,
+            numEntries,
+            mode,
+            op);
+        return;
+      }
+
+#ifdef POLYMPO_VERBOSE_MPI
+      if(self == 0 && !cudaAwareMPILogged){
+        std::cout
+            << "[CUDA_AWARE_MPI] Using batched single-buffer GPU-aware MPI path in communicate_and_take_halo_contributions1_improved()"
+            << "\n";
+        cudaAwareMPILogged = true;
+      }
+#endif
+
+      Kokkos::Timer timer;
+
+      if(!cudaAwareMPICacheValid){
+        cudaAwareMPICaches.clear();
+        cudaAwareMPICacheValid = true;
+      }
+
+      auto& cudaAwareCache =
+          cudaAwareMPICaches[std::make_pair(mode, numEntries)];
+
+      const bool needRebuild =
+          (!cudaAwareCache.valid) ||
+          (cudaAwareCache.cachedNumProcs != numProcsTot);
+
+      if(needRebuild){
+
+        cudaAwareCache.cachedNumProcs = numProcsTot;
+
+        cudaAwareCache.sendCounts.assign(numProcsTot, 0);
+        cudaAwareCache.recvCounts.assign(numProcsTot, 0);
+        cudaAwareCache.sendOffsets.assign(numProcsTot, 0);
+        cudaAwareCache.recvOffsets.assign(numProcsTot, 0);
+
+        for(int proc = 0; proc < numProcsTot; proc++){
+          if(proc == self) continue;
+
+          if(mode == 0){
+            cudaAwareCache.sendCounts[proc] = numOwnersOnOtherProcs[proc];
+            cudaAwareCache.recvCounts[proc] = numHalosOnOtherProcs[proc];
+          }
+          else{
+            cudaAwareCache.sendCounts[proc] = numHalosOnOtherProcs[proc];
+            cudaAwareCache.recvCounts[proc] = numOwnersOnOtherProcs[proc];
+          }
+        }
+
+        int totalSend = 0;
+        int totalRecv = 0;
+
+        for(int proc = 0; proc < numProcsTot; proc++){
+          cudaAwareCache.sendOffsets[proc] = totalSend;
+          totalSend += cudaAwareCache.sendCounts[proc];
+
+          cudaAwareCache.recvOffsets[proc] = totalRecv;
+          totalRecv += cudaAwareCache.recvCounts[proc];
+        }
+
+        cudaAwareCache.totalSendCount = totalSend;
+        cudaAwareCache.totalRecvCount = totalRecv;
+
+        cudaAwareCache.sendEntityGPU.allocate(totalSend);
+        cudaAwareCache.sendDataGPU.allocate(totalSend * numEntries);
+        cudaAwareCache.recvIDGPU.allocate(totalRecv);
+        cudaAwareCache.recvDataGPU.allocate(totalRecv * numEntries);
+
+        // ---- Build the flattened send-entity list (host, then one deep_copy) ----
+        if(totalSend > 0){
+          auto sendEntityCPU =
+              Kokkos::View<int*, Kokkos::HostSpace>(
+                  "sendEntityCPU_batched", totalSend);
+
+          if(mode == 0){
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+              if(cudaAwareCache.sendCounts[proc] <= 0) continue;
+
+              assert(haloOwnerLocalIDs[proc].size() ==
+                     static_cast<size_t>(cudaAwareCache.sendCounts[proc]));
+            }
+
+            std::vector<int> cursor(cudaAwareCache.sendOffsets);
+
+            for(int iEnt = 0; iEnt < numHalosTot; iEnt++){
+              int ownerProc = haloOwnerProcs[iEnt];
+              if(ownerProc == self) continue;
+
+              sendEntityCPU(cursor[ownerProc]) = numOwnersTot + iEnt;
+              cursor[ownerProc]++;
+            }
+
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+
+              assert(cursor[proc] ==
+                     cudaAwareCache.sendOffsets[proc] +
+                         cudaAwareCache.sendCounts[proc]);
+            }
+          }
+          else{
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+
+              int sendCount = cudaAwareCache.sendCounts[proc];
+              if(sendCount <= 0) continue;
+
+              assert(ownerOwnerLocalIDs[proc].size() ==
+                     static_cast<size_t>(sendCount));
+
+              int base = cudaAwareCache.sendOffsets[proc];
+
+              for(int i = 0; i < sendCount; i++){
+                sendEntityCPU(base + i) = ownerOwnerLocalIDs[proc][i];
+              }
+            }
+          }
+
+          Kokkos::deep_copy(cudaAwareCache.sendEntityGPU.view(), sendEntityCPU);
+        }
+
+        // ---- Build the flattened recv-ID list (host, then one deep_copy) ----
+        if(totalRecv > 0){
+          auto recvIDCPU =
+              Kokkos::View<int*, Kokkos::HostSpace>(
+                  "recvIDCPU_batched", totalRecv);
+
+          if(mode == 0){
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+
+              int recvCount = cudaAwareCache.recvCounts[proc];
+              if(recvCount <= 0) continue;
+
+              assert(ownerOwnerLocalIDs[proc].size() ==
+                     static_cast<size_t>(recvCount));
+
+              int base = cudaAwareCache.recvOffsets[proc];
+
+              for(int i = 0; i < recvCount; i++){
+                recvIDCPU(base + i) = ownerOwnerLocalIDs[proc][i];
+              }
+            }
+          }
+          else{
+            std::vector<int> cursor(cudaAwareCache.recvOffsets);
+
+            for(int iEnt = 0; iEnt < numHalosTot; iEnt++){
+              int ownerProc = haloOwnerProcs[iEnt];
+              if(ownerProc == self) continue;
+
+              recvIDCPU(cursor[ownerProc]) = numOwnersTot + iEnt;
+              cursor[ownerProc]++;
+            }
+
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+
+              assert(cursor[proc] ==
+                     cudaAwareCache.recvOffsets[proc] +
+                         cudaAwareCache.recvCounts[proc]);
+            }
+          }
+
+          Kokkos::deep_copy(cudaAwareCache.recvIDGPU.view(), recvIDCPU);
+        }
+
+        cudaAwareCache.valid = true;
+
+        if(mpiDiagnostics){
+            pumipic::RecordTime(label + "_MPI_Diagnostics_CacheBuild_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) +
+            "_rank" + std::to_string(self), timer.seconds());
+      }
+
+        timer.reset();
+      }
+
+      // ---- Pack: ONE kernel over all neighbors' send entities at once ----
+      if(cudaAwareCache.totalSendCount > 0){
+        auto sendEntityGPU = cudaAwareCache.sendEntityGPU.view();
+        auto sendDataGPU = cudaAwareCache.sendDataGPU.view();
+
+        Kokkos::parallel_for(
+            "pack cached cuda-aware mpi send buffer batched",
+            cudaAwareCache.totalSendCount,
+            KOKKOS_LAMBDA(const int i){
+              int entity = sendEntityGPU(i);
+
+              for(int k = 0; k < numEntries; k++){
+                sendDataGPU(i * numEntries + k) =
+                    meshField(entity, k);
+              }
+            });
+      }
+
+      Kokkos::fence();
+
+      if(mpiDiagnostics){
+         pumipic::RecordTime(
+         label + "_MPI_Diagnostics_Pack_m" + std::to_string(mode) + "_e" + std::to_string(numEntries) + "_rank" + std::to_string(self), timer.seconds());
+    }
+
+      timer.reset();
+
+      double postTime    = 0.0;
+      double waitallTime = 0.0;
+
+      std::vector<MPI_Request> requests;
+      requests.reserve(2 * numProcsTot);
+      int mpiError = MPI_SUCCESS;
+
+      // Data volume exchanged by this rank in this call (recorded once per
+      // call, independent of the post/wait timing below), tagged by the
+      // caller (label) so SD, VR, and Reconstruction can be told apart.
+      const double bytesSent =
+          static_cast<double>(cudaAwareCache.totalSendCount) * numEntries * sizeof(double);
+      const double bytesRecv =
+          static_cast<double>(cudaAwareCache.totalRecvCount) * numEntries * sizeof(double);
+
+      pumipic::RecordTime(label + "_MPI_BytesSent_" + std::to_string(self), bytesSent);
+      pumipic::RecordTime(label + "_MPI_BytesRecv_" + std::to_string(self), bytesRecv);
+
+      //Post both the Irecv and the matching Isend for a proc together, in
+      //the same loop iteration and under their own counts (recvCounts for
+      //Irecv, sendCounts for Isend). This replaces the previous two-pass
+      //version (a first loop that posted Irecv only, plus a second loop
+      //that posted Isend only) which existed only because of a leftover,
+      //commented-out duplicate of this same block.
+      int numNeighbors = 0;
+      for(int proc = 0; proc < numProcsTot; proc++){
+        if(proc == self) continue;
+        bool hasComm = false;
+
+        if(cudaAwareCache.recvCounts[proc] > 0){
+          MPI_Request reqData;
+
+          double* recvPtr =
+              cudaAwareCache.recvDataGPU.data() +
+              static_cast<size_t>(cudaAwareCache.recvOffsets[proc]) *
+                  numEntries;
+
+          mpiError = MPI_Irecv(
+              recvPtr,
+              cudaAwareCache.recvCounts[proc] * numEntries,
+              MPI_DOUBLE,
+              proc,
+              2,
+              comm,
+              &reqData);
+          if(mpiError != MPI_SUCCESS) break;
+          requests.push_back(reqData);
+          hasComm = true; 
+        }
+
+        if(cudaAwareCache.sendCounts[proc] > 0){
+          MPI_Request reqData;
+
+          double* sendPtr =
+              cudaAwareCache.sendDataGPU.data() +
+              static_cast<size_t>(cudaAwareCache.sendOffsets[proc]) *
+                  numEntries;
+
+          mpiError = MPI_Isend(
+              sendPtr,
+              cudaAwareCache.sendCounts[proc] * numEntries,
+              MPI_DOUBLE,
+              proc,
+              2,
+              comm,
+              &reqData);
+          if(mpiError != MPI_SUCCESS) break;
+          requests.push_back(reqData);
+          hasComm = true;
+        }
+         if(hasComm) numNeighbors++;
+      }
+      pumipic::RecordTime(label + "_MPI_NumNeighbors_" + std::to_string(self), static_cast<double>(numNeighbors));
+
+      postTime = timer.seconds();
+
+      pumipic::RecordTime(label + "_MPI_Post_" + std::to_string(self), postTime);
+
+      //Barrier here, not before posting: Isend/Irecv are non-blocking and
+      //their cost is local (looping + building MPI_Request objects), so a
+      //barrier before posting would only be measuring how skewed ranks
+      //were on entry to this function, which the calling function's own
+      //barriers already capture. Placed here, right before Waitall, it
+      //makes every rank enter Waitall at the same instant, so the
+      //waitallTime below reflects real message-arrival/network imbalance
+      //instead of being contaminated by skew left over from posting.
+      //Kokkos::Timer barrierTimer;
+      //MPI_Barrier(comm);
+      //const double InternalBarrierWait = barrierTimer.seconds();
+      //pumipic::RecordTime(label + "_MPI_BarrierWait_" + std::to_string(self), InternalBarrierWait);
+
+      timer.reset();
+
+      if(mpiError == MPI_SUCCESS && !requests.empty()){
+        mpiError = MPI_Waitall(
+            static_cast<int>(requests.size()),
+            requests.data(),
+            MPI_STATUSES_IGNORE);
+      }
+
+      waitallTime = timer.seconds();
+
+      pumipic::RecordTime(label + "_MPI_Waitall_" + std::to_string(self), waitallTime);
+
+      if(mpiError != MPI_SUCCESS){
+        cudaAwareMPIDisabled = true;
+
+        if(self == 0){
+          std::cout
+              << "[CUDA_AWARE_MPI] Batched device-pointer MPI failed."
+              << std::endl;
+        }
+
+        if(requests.empty()){
+          if(self == 0){
+            std::cout
+                << "[CUDA_AWARE_MPI] Falling back to CPU-staged communication."
+                << std::endl;
+          }
+
+          communicate_and_take_halo_contributions1(
+              meshField,
+              nEntities,
+              numEntries,
+              mode,
+              op);
+          return;
+        }
+
+        if(self == 0){
+          std::cout
+              << "[CUDA_AWARE_MPI] Failure happened after MPI requests were posted. "
+              << "Set POLYMPO_DISABLE_CUDA_AWARE_MPI=1 before running to force the CPU-staged path."
+              << std::endl;
+        }
+
+        MPI_Abort(comm, mpiError);
+        return;
+      }
+
+      timer.reset();
+
+      // ---- Unpack: ONE kernel over all neighbors' recv entities at once ----
+      if(cudaAwareCache.totalRecvCount > 0){
+        auto recvIDGPU = cudaAwareCache.recvIDGPU.view();
+        auto recvDataGPU = cudaAwareCache.recvDataGPU.view();
+
+        if(op == 0){
+          Kokkos::parallel_for(
+              "halo add cached cuda-aware mpi batched",
+              cudaAwareCache.totalRecvCount,
+              KOKKOS_LAMBDA(const int i){
+                const int vertex = recvIDGPU(i);
+
+                for(int k = 0; k < numEntries; k++){
+#ifdef POLYMPO_ASSUME_UNIQUE_HALO_CONTRIBS
+                  meshField(vertex, k) +=
+                      recvDataGPU(i * numEntries + k);
+#else
+                  Kokkos::atomic_add(
+                      &meshField(vertex, k),
+                      recvDataGPU(i * numEntries + k));
+#endif
+                }
+              });
+        }
+        else{
+          Kokkos::parallel_for(
+              "halo assign cached cuda-aware mpi batched",
+              cudaAwareCache.totalRecvCount,
+              KOKKOS_LAMBDA(const int i){
+                const int vertex = recvIDGPU(i);
+
+                for(int k = 0; k < numEntries; k++){
+                  meshField(vertex, k) =
+                      recvDataGPU(i * numEntries + k);
+                }
+              });
+        }
+      }
+
+      Kokkos::fence();
+
+      if(mpiDiagnostics){
+          pumipic::RecordTime(label + "_MPI_Diagnostics_Contribution_m" + std::to_string(mode) + "_e" +
+          std::to_string(numEntries) + "_rank" + std::to_string(self), timer.seconds());
+    }
+  
+
+    }
+
+#else
+
+    // Fallback path:
+    // if CUDA_AWARE_MPI is not defined, use the original GPU-CPU staging function.
+    template <typename ViewType>
+    void communicate_and_take_halo_contributions1_improved(
+        const ViewType& meshField,
+        int nEntities,
+        int numEntries,
+        int mode,
+        int op,
+        const std::string& label){
+      (void)label; // no per-call diagnostics on the CPU-staged fallback path
+
+      communicate_and_take_halo_contributions1(
+          meshField,
+          nEntities,
+          numEntries,
+          mode,
+          op);
+    }
+
+#endif
+
 };
 
 }//namespace polyMPO end
 
 #endif
-
diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp
index 9f184e05..d57bf4c7 100644
--- a/src/pmpo_MPMesh_assembly.hpp
+++ b/src/pmpo_MPMesh_assembly.hpp
@@ -96,12 +96,27 @@ void MPMesh::assemblyElm0() {
 }
 
 void MPMesh::reconstruct_coeff_full(){
-  Kokkos::Timer timer;
+
   int self, numProcsTot;
   MPI_Comm comm = p_MPs->getMPIComm();
   MPI_Comm_rank(comm, &self);
   MPI_Comm_size(comm, &numProcsTot);
-  
+
+  Kokkos::fence();      //B0: drain any device work left from the previous phase
+  //MPI_Barrier(comm);    //B0: align all ranks so this timed region starts at the same instant
+
+  Kokkos::Timer totalTimer;
+  Kokkos::Timer phaseTimer;
+
+  double preCommunicationComputeTime       = 0.0;
+  //double preCommunicationComputeImbalance  = 0.0;
+  double gatherCommunicationTime           = 0.0;
+  //double gatherCommunicationImbalance      = 0.0;
+  double scatterCommunicationTime          = 0.0;
+  //double scatterCommunicationImbalance     = 0.0;
+  double postCommunicationComputeTime      = 0.0;
+  //double postCommunicationComputeImbalance = 0.0;
+
   static int coeff_count=0;
   if(!self) std::cout<<"===="<<__FUNCTION__<<" "<<coeff_count<<"===="<<std::endl;
   coeff_count++;
@@ -114,6 +129,8 @@ void MPMesh::reconstruct_coeff_full(){
   auto dual_triangle_area=p_mesh->getMeshField<MeshF_DualTriangleArea>();
 
   //Material Points
+  phaseTimer.reset();
+
   calcBasis();
 
   auto weight = p_MPs->getData<MPF_Basis_Vals>();
@@ -152,30 +169,69 @@ void MPMesh::reconstruct_coeff_full(){
   };
   p_MPs->parallel_for(assemble, "assembly");
   Kokkos::fence();
-  pumipic::RecordTime("Assemble Matrix Per Process" + std::to_string(self), timer.seconds());
+
+  preCommunicationComputeTime = phaseTimer.seconds();          //T_before: pure local compute time, no waiting
+  //MPI_Barrier(comm);                                            //B1: fast ranks wait here for the slowest rank
+  //preCommunicationComputeImbalance = phaseTimer.seconds() - preCommunicationComputeTime;
+
   //Mode 0 is Gather:  Halos Send to Owners
   //Mode 1 is Scatter: Owners Send to Halos
   //Op 0 is addition
   //Op 1 is replacement
-  timer.reset();
+
   int mode = 0;
   int op = 0;
   if (numProcsTot >1){
-    communicate_and_take_halo_contributions1(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
+
+    phaseTimer.reset();
+    communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op, "Reconstruction_Gather");
+    
+    Kokkos::fence();
+    gatherCommunicationTime = phaseTimer.seconds();             //T_before
+    //MPI_Barrier(comm);                                            //B2: fast ranks wait here for the slowest rank
+    //gatherCommunicationImbalance = phaseTimer.seconds() - gatherCommunicationTime;
+
     mode=1; 
     op=1;
-    communicate_and_take_halo_contributions1(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
+    
+    phaseTimer.reset();
+    communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op, "Reconstruction_Scatter");
+    Kokkos::fence();
+    scatterCommunicationTime = phaseTimer.seconds();            //T_before
+    //MPI_Barrier(comm);                                            //B3: fast ranks wait here for the slowest rank
+    //scatterCommunicationImbalance = phaseTimer.seconds() - scatterCommunicationTime;
   }
-  pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds());
- 
-  //Stroe the 1st matrix element
+  
+  phaseTimer.reset();
+  //Store the 1st matrix element
   Kokkos::View<double*>vtxMatrixMass_l("vtxMass", numVertices);
   Kokkos::parallel_for("storeMatrixMass", numVertices, KOKKOS_LAMBDA(const int vtx){
     vtxMatrixMass_l(vtx) = vtxMatrices(vtx, 0);
   }); 
+  Kokkos::fence();
   this->vtxMatrixMass = vtxMatrixMass_l;
 
   invertMatrix(vtxMatrices, radius);
+  Kokkos::fence();
+  postCommunicationComputeTime = phaseTimer.seconds();          //T_before
+  //MPI_Barrier(comm);                                              //B4: fast ranks wait here for the slowest rank
+  //postCommunicationComputeImbalance = phaseTimer.seconds() - postCommunicationComputeTime;
+
+  const double computeTime = preCommunicationComputeTime + postCommunicationComputeTime;
+  const double communicationTime = gatherCommunicationTime + scatterCommunicationTime;
+  const double totalTime = totalTimer.seconds();
+
+  pumipic::RecordTime("Reconstruction_PreComm_Compute_" + std::to_string(self), preCommunicationComputeTime);
+  //pumipic::RecordTime("Reconstruction_PreComm_Compute_Imbalance_" + std::to_string(self), preCommunicationComputeImbalance);
+  pumipic::RecordTime("Reconstruction_Gather_Communication_" + std::to_string(self),gatherCommunicationTime);
+  //pumipic::RecordTime("Reconstruction_Gather_Communication_Imbalance_" + std::to_string(self),gatherCommunicationImbalance);
+  pumipic::RecordTime("Reconstruction_Scatter_Communication_" + std::to_string(self), scatterCommunicationTime);
+  //pumipic::RecordTime("Reconstruction_Scatter_Communication_Imbalance_" + std::to_string(self), scatterCommunicationImbalance);
+  pumipic::RecordTime("Reconstruction_PostComm_Compute_" + std::to_string(self), postCommunicationComputeTime);
+  //pumipic::RecordTime("Reconstruction_PostComm_Compute_Imbalance_" + std::to_string(self), postCommunicationComputeImbalance);
+  pumipic::RecordTime("Reconstruction_Compute_Total_" + std::to_string(self), computeTime);
+  pumipic::RecordTime("Reconstruction_Communication_Total_" + std::to_string(self), communicationTime);
+  pumipic::RecordTime("Reconstruction_Total_" + std::to_string(self), totalTime);
 }
 
 void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius){
@@ -321,13 +377,17 @@ void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const doubl
 
 template <MeshFieldIndex meshFieldIndex>
 void MPMesh::assemblyVtx1(){
-  Kokkos::Timer timer;
-
   int self, numProcsTot;
   MPI_Comm comm = p_MPs->getMPIComm();
   MPI_Comm_rank(comm, &self);
   MPI_Comm_size(comm, &numProcsTot);
 
+  Kokkos::fence();      //B0: drain any device work left from the previous phase
+  //MPI_Barrier(comm);    //B0: align all ranks so this timed region starts at the same instant
+
+  Kokkos::Timer totalTimer;
+  Kokkos::Timer computeTimer;
+
   auto VtxCoeffs_new=this->precomputedVtxCoeffs_new;
 
   //Mesh Information
@@ -342,16 +402,24 @@ void MPMesh::assemblyVtx1(){
   p_mesh->fillMeshField<meshFieldIndex>(numVtx, numEntries, 0.0);
   auto meshField = p_mesh->getMeshField<meshFieldIndex>();
 
+  auto vtxRotLon = p_mesh->getMeshField<MeshF_VtxRotLon>();
+
   //Material Points
   auto mpData = p_MPs->getData<mpfIndex>();
   auto weight = p_MPs->getData<MPF_Basis_Vals>();
   auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
-
+  auto curPosRotLatLon = p_MPs->getData<MPF_Cur_Pos_Rot_Lat_Lon>();
+  auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
   //Earth Radius
   double radius = 1.0;
   if(p_mesh->getGeomType() == geom_spherical_surf)
     radius=p_mesh->getSphereRadius();
 
+  bool use_correction_term = false;
+  if constexpr (meshFieldIndex == MeshF_Vel) {
+    use_correction_term = true;
+  }
+
   //Reconstruct
   auto reconstruct = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
@@ -367,22 +435,58 @@ void MPMesh::assemblyVtx1(){
                                                        VtxCoeffs_new(vID,0, 2)*CoordDiffs[2] +
                                                        VtxCoeffs_new(vID,0, 3)*CoordDiffs[3]);
 
-        for (int k=0; k<numEntries; k++){
-          auto val = factor*mpData(mp,k);
-          Kokkos::atomic_add(&meshField(vID,k), val);
+        if (use_correction_term){
+          double transport[2] = {0.0};
+          seaice_mpm_coord_parallel_transport(curPosRotLatLon(mp, 1), vtxRotLon(vID), curPosRotLatLon(mp, 0), transport);
+          const double u =  transport[0] * mpData(mp,0) + transport[1] * mpData(mp,1);
+          const double v = -transport[1] * mpData(mp,0) + transport[0] * mpData(mp,1);
+
+          Kokkos::atomic_add(&meshField(vID,0), factor * u);
+          Kokkos::atomic_add(&meshField(vID,1), factor * v);
+        }
+        else{
+          for (int k=0; k<numEntries; k++){
+            auto val = factor*mpData(mp,k);
+            Kokkos::atomic_add(&meshField(vID,k), val);
+          }
         }
       }
     }
   };
   p_MPs->parallel_for(reconstruct, "reconstruct");
   Kokkos::fence();
-  pumipic::RecordTime("Assemble Field per process" + std::to_string(self), timer.seconds());
+  const double computeTime = computeTimer.seconds();          //T_before: pure local compute time, no waiting
+
+  //MPI_Barrier(comm);                                            //B1: fast ranks wait here for the slowest rank
+  //const double computeTimeSync = computeTimer.seconds();
+  //const double computeImbalance = computeTimeSync - computeTime;
+
+  Kokkos::Timer communicationTimer;
 
-  timer.reset();
   if(numProcsTot>1){ 
-    communicate_and_take_halo_contributions1(meshField, numVertices, numEntries, 0, 0);
+    communicate_and_take_halo_contributions1_improved(meshField, numVertices, numEntries, 0, 0, "Velocity_Reconstruction");
   }
-  pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds());
+  Kokkos::fence();
+
+  const double communicationTime = communicationTimer.seconds();  //T_before
+
+  //MPI_Barrier(comm);                                            //B2: fast ranks wait here for the slowest rank
+  //const double communicationTimeSync = communicationTimer.seconds();
+  //const double communicationImbalance = communicationTimeSync - communicationTime;
+
+  const double totalTime = totalTimer.seconds();
+  
+  if constexpr (meshFieldIndex == MeshF_Vel) {
+    pumipic::RecordTime("Velocity_Reconstruction_Compute_" + std::to_string(self), computeTime);
+    //pumipic::RecordTime("Velocity_Reconstruction_Compute_Imbalance_" + std::to_string(self),computeImbalance);
+
+    pumipic::RecordTime("Velocity_Reconstruction_Communication_" + std::to_string(self),communicationTime);
+    //pumipic::RecordTime("Velocity_Reconstruction_Communication_Imbalance_" + std::to_string(self),communicationImbalance);
+
+    pumipic::RecordTime("Velocity_Reconstruction_Total_" + std::to_string(self),totalTime);
+
+}
+
 }
 
 template <MeshFieldIndex meshFieldIndex>
diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp
index 53c330e1..1bf5edf6 100644
--- a/src/pmpo_c.cpp
+++ b/src/pmpo_c.cpp
@@ -138,7 +138,13 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh,
      new polyMPO::MaterialPoints(numElms, numActiveMPs, mpsPerElm_d, active_mp2Elm_d, active_mpIDs_d, elm2global);
 
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
-  p_MPs->setElmIDoffset(offset);  
+  p_MPs->setElmIDoffset(offset);
+
+  int self;
+  MPI_Comm_rank(p_MPs->getMPIComm(), &self);
+  std::cout << "[RankSummary] Rank=" << self
+            << " MaterialPoints=" << p_MPs->getCount()
+            << std::endl;  
 }
 
 void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh,
@@ -567,10 +573,13 @@ void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs
   pumipic::RecordTime("PolyMPO_getMPMass", timer.seconds());
 }
 
-void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn) {
+void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn, const int callSiteId) {
   Kokkos::Timer timer;
   checkMPMeshValid(p_mpmesh);
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  int self;
+  MPI_Comm comm = p_MPs->getMPIComm();
+  MPI_Comm_rank(comm, &self);
   PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries);
   PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
   //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
@@ -587,7 +596,7 @@ void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs,
     }
   };
   p_MPs->parallel_for(setMPVel, "setMPVel");
-  pumipic::RecordTime("PolyMPO_setMPVel", timer.seconds());
+  pumipic::RecordTime("PolyMPO_setMPVel_site" + std::to_string(callSiteId) + "_" + std::to_string(self), timer.seconds());
 }
 
 void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost) {
@@ -728,7 +737,6 @@ void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numM
 void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost){
   Kokkos::Timer timer;
   checkMPMeshValid(p_mpmesh);
-  
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   //Rank information
   int self;
@@ -754,10 +762,34 @@ void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs
   pumipic::RecordTime("PolyMPO_setMPArea" + std::to_string(self), timer.seconds());
 }
 
+void polympo_getAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost) {
+  Kokkos::Timer timer;
+  checkMPMeshValid(p_mpmesh);
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+
+  PMT_ALWAYS_ASSERT(nComps == 1);
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+
+  auto mpArea = p_MPs->getData<polyMPO::MPF_Area>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+
+  Kokkos::View<double**> mpAreaCopy("mpAreaCopy", nComps, numMPs);
+  auto getMPArea = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
+    if(mask){
+      mpAreaCopy(0,mpAppID(mp)) = mpArea(mp,0);
+    }
+  };
+  p_MPs->parallel_for(getMPArea, "getMPArea");
+  kkDbl2dViewHostU arrayHost(areaMPHost, nComps, numMPs);
+  Kokkos::deep_copy(arrayHost, mpAreaCopy);
+  pumipic::RecordTime("PolyMPO_getMPArea", timer.seconds());
+}
+
+
 void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* icePressureMPHost){
   Kokkos::Timer timer;
   checkMPMeshValid(p_mpmesh);
-  
+
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   //Rank information
   int self;
@@ -783,6 +815,73 @@ void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int
   pumipic::RecordTime("PolyMPO_setIcePressure" + std::to_string(self), timer.seconds());
 }
 
+void polympo_setOceanVelocity_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* uArray, const double* vArray){
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries);
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices);
+  //copy the host array to the device
+  auto oceanVelocity = p_mesh->getMeshField<polyMPO::MeshF_OceanVelocity>();
+  auto h_oceanVelocity = Kokkos::create_mirror_view(oceanVelocity);
+  for(int i=0; i<nVertices; i++){
+    h_oceanVelocity(i, 0) = uArray[i];
+    h_oceanVelocity(i, 1) = vArray[i];
+  }
+  Kokkos::deep_copy(oceanVelocity, h_oceanVelocity);
+}
+
+void polympo_subcycle_prep_arrays_f(MPMesh_ptr p_mpmesh){
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  int numVerticesOwned = p_mesh->getNumVerticesOwned();
+
+  auto solveVelocity = p_mesh->getMeshField<polyMPO::MeshF_SolveVelocity>();
+  auto velField = p_mesh->getMeshField<polyMPO::MeshF_Vel>();
+  auto stress_divUV = p_mesh->getMeshField<polyMPO::MeshF_StressDivergence>();
+  auto oceanStress = p_mesh->getMeshField<polyMPO::MeshF_OceanStress>();
+
+  Kokkos::parallel_for("prep_arrays", numVerticesOwned, KOKKOS_LAMBDA(const int vtx){
+    if(solveVelocity(vtx)==0){
+      for (int k=0; k<2; k ++){
+        velField(vtx, k) = 0.0;
+        stress_divUV(vtx, k) = 0.0;
+        oceanStress(vtx, k) = 0.0;
+      }
+    }
+  });
+}
+
+void polympo_setReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost){
+  Kokkos::Timer timer;
+  checkMPMeshValid(p_mpmesh);
+
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  //Rank information
+  int self;
+  MPI_Comm comm = p_MPs->getMPIComm();
+  MPI_Comm_rank(comm, &self);
+  //Asserts
+  PMT_ALWAYS_ASSERT(nComps == 1);
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+  //MP Data
+  auto mpReplacementPressure = p_MPs->getData<polyMPO::MPF_ReplacementPressure>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+
+  //Copy to device
+  kkViewHostU<const double**> mpRepPressure_h(replacementPressureMPHost, nComps, numMPs);
+  Kokkos::View<double**> mpRepPressure_d("mpRepPressureDevice", nComps, numMPs);
+  Kokkos::deep_copy(mpRepPressure_d, mpRepPressure_h);
+  //Set in PS
+  auto setMPRepPressure = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
+    if(mask){
+      mpReplacementPressure(mp,0) = mpRepPressure_d(0, mpAppID(mp));
+    }
+  };
+  p_MPs->parallel_for(setMPRepPressure, "setMPRepPressure");
+  pumipic::RecordTime("PolyMPO_setReplacementPressure" + std::to_string(self), timer.seconds());
+}
+
 void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost){
   Kokkos::Timer timer;
   checkMPMeshValid(p_mpmesh);
@@ -1117,6 +1216,25 @@ void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double
   }
 }
 
+void polympo_setMeshVtxRotLon_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* longitude){
+  Kokkos::Timer timer;
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  //check the size
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
+
+  //copy the host array to the device
+  auto coordsArray = p_mesh->getMeshField<polyMPO::MeshF_VtxRotLon>();
+  auto h_coordsArray = Kokkos::create_mirror_view(coordsArray);
+  for(int i=0; i<nVertices; i++){
+    h_coordsArray(i) = longitude[i];
+  }
+  Kokkos::deep_copy(coordsArray, h_coordsArray);
+  pumipic::RecordTime("PolyMPO_setMeshVtxRotLon", timer.seconds());
+}
+
 void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelIn, const double* vVelIn){
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
@@ -1135,11 +1253,15 @@ void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const dou
   Kokkos::deep_copy(coordsArray, h_coordsArray);
 }
 
-void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelOut, double* vVelOut){
+void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelOut, double* vVelOut, const int callSiteId){
   Kokkos::Timer timer;
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  int self;
+  MPI_Comm comm = p_MPs->getMPIComm();
+  MPI_Comm_rank(comm, &self);
 
   //check the size
   PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); 
@@ -1151,7 +1273,7 @@ void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* u
     uVelOut[i] = h_coordsArray(i,0);
     vVelOut[i] = h_coordsArray(i,1);
   }
-  pumipic::RecordTime("PolyMPO_getMeshVtxVel", timer.seconds());
+  pumipic::RecordTime("PolyMPO_getMeshVtxVel_site" + std::to_string(callSiteId) + "_" + std::to_string(self), timer.seconds());
 }
 
 void polympo_setMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* vtxMass){
@@ -1448,11 +1570,33 @@ void polympo_setSolveVelocityMesh_f(MPMesh_ptr p_mpmesh, const int nVertices, in
   Kokkos::deep_copy(solveVelocity, h_solveVelocity);
 }
 
+void polympo_setIceAreaVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
+  //copy the host array to the device
+  auto iceArea = p_mesh->getMeshField<polyMPO::MeshF_VtxMass>();
+  auto h_iceArea = Kokkos::create_mirror_view(iceArea);
+  for(int i=0; i<nVertices; i++)
+    h_iceArea(i, 0) = array[i];
+  Kokkos::deep_copy(iceArea, h_iceArea);
+}
+
 void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh){
   //chech validity
   checkMPMeshValid(p_mpmesh);
+  Kokkos::Timer wrapperTimer;
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  int self;
+  MPI_Comm comm = p_MPs->getMPIComm();
+  MPI_Comm_rank(comm, &self);
+
   ((polyMPO::MPMesh*)p_mpmesh) -> calculateStressDivergence();
 
+  pumipic::RecordTime("Wrapper_StressDivergence_Total_" + std::to_string(self), wrapperTimer.seconds());
+
  }
 
 void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray){
@@ -1570,6 +1714,11 @@ void polympo_set_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const int nVertic
   Kokkos::deep_copy(oceanStressCoeff, h_oceanStressCoeff);
 }
 
+void polympo_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const double configIceOceanDragCoeff){
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  p_mesh->calcOceanStressCoeff(configIceOceanDragCoeff);
+}
+
 void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh){
   //Temporary grid Solve after calculateDivergence
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
@@ -1605,11 +1754,19 @@ void polympo_set_free_slip_bc_f(MPMesh_ptr p_mpmesh){
 }
 
 void polympo_set_halo_vel_from_owner_f(MPMesh_ptr p_mpmesh){
+  
+  int numProcsTot;
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  MPI_Comm comm = p_MPs->getMPIComm();
+  MPI_Comm_size(comm, &numProcsTot);
+  if(numProcsTot == 1) return;
+
   auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh);
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
   int numVertices = p_mesh->getNumVertices();
   auto vtxFieldVel = p_mesh->getMeshField<polyMPO::MeshF_Vel>();
-  mpMesh->communicate_and_take_halo_contributions1(vtxFieldVel, numVertices, 2, 1, 1); 
+
+  mpMesh->communicate_and_take_halo_contributions1_improved(vtxFieldVel, numVertices, 2, 1, 1, "halo_vel_from_owner"); 
 }
 
 //Advection Calcualtions
@@ -1688,7 +1845,14 @@ void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh){
 void polympo_reconstruct_coeff_with_MPI_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
   auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
+  Kokkos::Timer wrapperTimer;
+  int self;
+  MPI_Comm comm = mpmesh->p_MPs->getMPIComm();
+  MPI_Comm_rank(comm, &self);
+
   mpmesh->reconstruct_coeff_full();
+
+  pumipic::RecordTime("Wrapper_ReconstructCoeff_Total_" + std::to_string(self), wrapperTimer.seconds());
 }
 
 void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh){
@@ -1700,7 +1864,14 @@ void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh){
 void polympo_reconstruct_velocity_with_MPI_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
   auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
+  Kokkos::Timer wrapperTimer;
+  int self;
+  MPI_Comm comm = mpmesh->p_MPs->getMPIComm();
+  MPI_Comm_rank(comm, &self);
+
   mpmesh->assemblyVtx1<polyMPO::MeshF_Vel>();
+
+ pumipic::RecordTime("Wrapper_ReconstructVelocity_Total_" + std::to_string(self), wrapperTimer.seconds());
 }
 
 void polympo_init_deludelvDyn_f(MPMesh_ptr p_mpmesh){
@@ -1726,7 +1897,7 @@ void polympo_aggregate_deluDyn_f(MPMesh_ptr p_mpmesh){
 }
 
 void polympo_finalize_deludelvDyn_f(MPMesh_ptr p_mpmesh){
-  
+
   checkMPMeshValid(p_mpmesh);
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
 
@@ -1737,7 +1908,7 @@ void polympo_finalize_deludelvDyn_f(MPMesh_ptr p_mpmesh){
   auto vtxField = p_mesh->getMeshField<polyMPO::MeshF_OnSurfDispIncr>();
   auto vtxFieldVel      = p_mesh->getMeshField<polyMPO::MeshF_Vel>();
   auto vtxFieldVel_incr = p_mesh->getMeshField<polyMPO::MeshF_OnSurfVeloIncr>();
- 
+
   Kokkos::parallel_for("Finalize_increments", nVertices, KOKKOS_LAMBDA(const int vtx){
     vtxField(vtx, 0) = vtxField(vtx, 0) * elasticTimeStep;
     vtxField(vtx, 1) = vtxField(vtx, 1) * elasticTimeStep;
diff --git a/src/pmpo_c.h b/src/pmpo_c.h
index 23d78b70..b938570f 100644
--- a/src/pmpo_c.h
+++ b/src/pmpo_c.h
@@ -44,7 +44,7 @@ void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const in
 //MP fields
 void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn);
 void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpMassHost);
-void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn);
+void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn, const int callSiteId);
 void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost);
 void polympo_calculateMPStrainRate_f(MPMesh_ptr p_mpmesh);
 void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStrainRateIn);
@@ -53,7 +53,11 @@ void polympo_calculateMPStress_f(MPMesh_ptr p_mpmesh, const int constitutive_mod
 void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn);
 void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost);
 void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost);
+void polympo_getAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost);
 void polympo_setIcePressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* icePressureMPHost);
+void polympo_setOceanVelocity_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* uArray, const double* vArray);
+void polympo_subcycle_prep_arrays_f(MPMesh_ptr p_mpmesh);
+void polympo_setReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost);
 void polympo_getReplacementPressureMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* replacementPressureMPHost);
 
 //Mesh info
@@ -85,9 +89,10 @@ int polympo_getMeshFElmType_f();
 void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* xArray, const double* yArray, const double* zArray);
 void polympo_getMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, double* xArray, double* yArray, double* zArray);
 void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* latitude);
+void polympo_setMeshVtxRotLon_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* longitude);
 void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double* latitude);
 void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelocity, const double* vVelocity);
-void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelocity, double* vVelocity);
+void polympo_getMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uVelocity, double* vVelocity, const int callSiteId);
 void polympo_setMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* vtxMass);
 void polympo_getMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, double* vtxMass);
 void polympo_setMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, const double* elmMass);
@@ -106,6 +111,7 @@ void polympo_setElasticTimeStep_f(MPMesh_ptr p_mpmesh, const double elasticTimeS
 void polympo_setDynamicTimeStep_f(MPMesh_ptr p_mpmesh, const double dynamicTimeStep);
 void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* array);
 void polympo_setSolveVelocityMesh_f(MPMesh_ptr p_mpmesh, const int nVertices, int* array);
+void polympo_setIceAreaVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
 void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh);
 void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray);
 void polympo_setTotalMassVtx_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
@@ -114,6 +120,7 @@ void polympo_set_surfaceTiltForce_f(MPMesh_ptr p_mpmesh, const int nVertices, do
 void polympo_set_totalMassVertexfVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
 void polympo_set_oceanStress_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray);
 void polympo_set_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
+void polympo_calculate_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const double configIceOceanDragCoeff);
 void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh);
 void polympo_set_boundary_normal_vertex_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* uArray, double* vArray);
 void polympo_set_free_slip_bc_f(MPMesh_ptr p_mpmesh);
diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90
index d0bb903c..e07db516 100644
--- a/src/pmpo_fortran.f90
+++ b/src/pmpo_fortran.f90
@@ -340,12 +340,13 @@ subroutine polympo_getMPMass(mpMesh, nComps, numMPs, array) &
   !> @param numMPs(in) number of the MPs
   !> @param array(in) input MP velocity 1D array (numMPs*2)
   !---------------------------------------------------------------------------
-  subroutine polympo_setMPVel(mpMesh, nComps, numMPs, array) &
+  subroutine polympo_setMPVel(mpMesh, nComps, numMPs, array, callSiteId) &
              bind(C, NAME='polympo_setMPVel_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
     integer(c_int), value :: nComps, numMPs
     type(c_ptr), intent(in), value :: array
+    integer(c_int), value :: callSiteId
   end subroutine
 
   !---------------------------------------------------------------------------
@@ -369,7 +370,7 @@ subroutine polympo_calculateMPStrainRate(mpMesh) &
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
   end subroutine
- 
+
   subroutine polympo_setMPStrainRate(mpMesh, nComps, numMPs, array) &
              bind(C, NAME='polympo_setMPStrainRate_f')
     use :: iso_c_binding
@@ -386,7 +387,7 @@ subroutine polympo_getMPStrainRate(mpMesh, nComps, numMPs, array) &
     type(c_ptr), value :: array
   end subroutine
 
- 
+
   !MP Stress
   subroutine polympo_calculateMPStress(mpMesh, constitutive_model) &
              bind(C, NAME='polympo_calculateMPStress_f')
@@ -419,6 +420,14 @@ subroutine polympo_setAreaMP(mpMesh, nComps, numMPs, array) &
     type(c_ptr), value :: array
   end subroutine
 
+  subroutine polympo_getAreaMP(mpMesh, nComps, numMPs, array) &
+             bind(C, NAME='polympo_getAreaMP_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nComps, numMPs
+    type(c_ptr), value :: array
+  end subroutine
+
   subroutine polympo_setIcePressureMP(mpMesh, nComps, numMPs, array) &
              bind(C, NAME='polympo_setIcePressureMP_f')
     use :: iso_c_binding
@@ -426,7 +435,29 @@ subroutine polympo_setIcePressureMP(mpMesh, nComps, numMPs, array) &
     integer(c_int), value :: nComps, numMPs
     type(c_ptr), value :: array
   end subroutine
-  
+
+  subroutine polympo_setOceanVelocity(mpMesh, nComps, nVertices, uArray, vArray) &
+             bind(C, NAME='polympo_setOceanVelocity_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nComps, nVertices
+    type(c_ptr), value :: uArray, vArray
+  end subroutine
+
+  subroutine polympo_subcycle_prep_arrays(mpMesh) &
+             bind(C, NAME='polympo_subcycle_prep_arrays_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end subroutine
+
+  subroutine polympo_setReplacementPressureMP(mpMesh, nComps, numMPs, array) &
+             bind(C, NAME='polympo_setReplacementPressureMP_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nComps, numMPs
+    type(c_ptr), value :: array
+  end subroutine
+
   subroutine polympo_getReplacementPressureMP(mpMesh, nComps, numMPs, array) &
              bind(C, NAME='polympo_getReplacementPressureMP_f')
     use :: iso_c_binding
@@ -756,6 +787,14 @@ subroutine polympo_getMeshVtxRotLat(mpMesh, nVertices, latitude) &
     type(c_ptr), value :: latitude
   end subroutine
 
+  subroutine polympo_setMeshVtxRotLon(mpMesh, nVertices, longitude) &
+             bind(C, NAME='polympo_setMeshVtxRotLon_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), intent(in), value :: longitude
+  end subroutine
+
   !---------------------------------------------------------------------------
   !> @brief set the vertices velocity from a host array
   !> @param mpmesh(in/out) MPMesh object
@@ -780,12 +819,13 @@ subroutine polympo_setMeshVtxVel(mpMesh, nVertices, uVel, vVel) &
   !> @param vVel(in/out) output vertices v-component velocity
   !>        1D array (numVtx), allocated by user
   !---------------------------------------------------------------------------
-  subroutine polympo_getMeshVtxVel(mpMesh, nVertices, uVel, vVel) &
+  subroutine polympo_getMeshVtxVel(mpMesh, nVertices, uVel, vVel, callSiteId) &
              bind(C, NAME='polympo_getMeshVtxVel_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
     integer(c_int), value :: nVertices
     type(c_ptr), value :: uVel, vVel
+    integer(c_int), value :: callSiteId
   end subroutine
 
   !---------------------------------------------------------------------------
@@ -994,7 +1034,7 @@ subroutine polympo_setDynamicTimeStep(mpMesh, dynamicTimeStep) &
     type(c_ptr), value :: mpMesh
     real(c_double), value :: dynamicTimeStep
   end subroutine
-  
+
   subroutine polympo_setSolveStressMesh(mpMesh, nCells, array) &
              bind(C, NAME='polympo_setSolveStressMesh_f')
     use :: iso_c_binding
@@ -1011,6 +1051,14 @@ subroutine polympo_setSolveVelocityMesh(mpMesh, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
 
+ subroutine polympo_setIceAreaVertex(mpMesh, nVertices, array) &
+             bind(C, NAME='polympo_setIceAreaVertex_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), value :: array
+  end subroutine
+
   subroutine polympo_calculateStressDivergence(mpMesh) &
              bind(C, NAME='polympo_calculateStressDivergence_f')
     use :: iso_c_binding
@@ -1040,7 +1088,7 @@ subroutine polympo_set_airStress(mpMesh, nVertices, uArray, vArray) &
     integer(c_int), value :: nVertices
     type(c_ptr), value :: uArray, vArray
   end subroutine
-  
+
   subroutine polympo_set_surfaceTiltForce(mpMesh, nVertices, uArray, vArray) &
              bind(C, NAME='polympo_set_surfaceTiltForce_f')
     use :: iso_c_binding
@@ -1073,6 +1121,13 @@ subroutine polympo_set_oceanStressCoefficient(mpMesh, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
 
+  subroutine polympo_calculate_oceanStressCoefficient(mpMesh, configIceOceanDragCoeff) &
+             bind(C, NAME='polympo_calculate_oceanStressCoefficient_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    real(c_double), value::configIceOceanDragCoeff
+  end subroutine
+
   subroutine  polympo_velocity_grid_solve(mpMesh) &
               bind(C, NAME='polympo_velocity_grid_solve_f')
     use :: iso_c_binding
@@ -1091,13 +1146,13 @@ subroutine polympo_set_free_slip_bc(mpMesh) &
              bind(C, NAME='polympo_set_free_slip_bc_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
-  end subroutine 
+  end subroutine
 
   subroutine polympo_set_halo_vel_from_owner(mpMesh) &
              bind(C, NAME='polympo_set_halo_vel_from_owner_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
-  end subroutine 
+  end subroutine
 
   !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
@@ -1110,13 +1165,13 @@ subroutine polympo_push(mpMesh) &
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
   end subroutine
- 
+
   subroutine polympo_push_ahead(mpMesh) &
              bind(C, NAME='polympo_push_ahead_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
   end subroutine
- 
+
   !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
   !---------------------------------------------------------------------------
diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp
index 6a48615f..c47f911b 100644
--- a/src/pmpo_materialPoints.hpp
+++ b/src/pmpo_materialPoints.hpp
@@ -148,7 +148,7 @@ class MaterialPoints {
     void rebuild(IntView addedMP2elm, IntView addedMPAppID);
     void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID, Kokkos::View<const int*> addedMPMask);
     void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID);
-    
+
     void finishRebuild();
     bool rebuildOngoing();
 
diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp
index 0e7849ca..4d7b2eae 100644
--- a/src/pmpo_mesh.cpp
+++ b/src/pmpo_mesh.cpp
@@ -20,6 +20,10 @@ namespace polyMPO{
     PMT_ALWAYS_ASSERT(vtxRotLatMapEntry.first == MeshFType_VtxBased);
     vtxRotLat_ = MeshFView<MeshF_VtxRotLat>(vtxRotLatMapEntry.second,numVtxs_);
 
+    auto vtxRotLonMapEntry = meshFields2TypeAndString.at(MeshF_VtxRotLon);
+    PMT_ALWAYS_ASSERT(vtxRotLonMapEntry.first == MeshFType_VtxBased);
+    vtxRotLon_ = MeshFView<MeshF_VtxRotLon>(vtxRotLonMapEntry.second,numVtxs_);
+
     auto vtxVelMapEntry = meshFields2TypeAndString.at(MeshF_Vel);
     PMT_ALWAYS_ASSERT(vtxVelMapEntry.first == MeshFType_VtxBased);
     vtxVel_ = MeshFView<MeshF_Vel>(vtxVelMapEntry.second,numVtxs_);
@@ -78,6 +82,9 @@ namespace polyMPO{
     oceanStress_ = MeshFView<MeshF_OceanStress>(meshFields2TypeAndString.at(MeshF_OceanStress).second, numVtxs_);
     
     oceanStressCoeff_ = MeshFView<MeshF_OceanStressCoeff>(meshFields2TypeAndString.at(MeshF_OceanStressCoeff).second, numVtxs_);
+
+    oceanVelocity_ = MeshFView<MeshF_OceanVelocity>(meshFields2TypeAndString.at(MeshF_OceanVelocity).second, numVtxs_);
+
   }
 
   void Mesh::setMeshElmBasedFieldSize(){
@@ -157,6 +164,21 @@ namespace polyMPO{
     });
   }
 
+  void Mesh::calcOceanStressCoeff(const double configIceOceanDragCoeff){
+    int numVerticesOwned = getNumVerticesOwned();
+    auto iceAreaVtx = getMeshField<polyMPO::MeshF_VtxMass>();
+    auto oceanStressCoeff = getMeshField<MeshF_OceanStressCoeff>();
+    auto velocity = getMeshField<MeshF_Vel>();
+    auto solve_velocity = getMeshField<MeshF_SolveVelocity>();
+    auto oceanVelocity = getMeshField<MeshF_OceanVelocity>();
+    auto seaiceDensitySeaWater_ = polyMPO::seaiceDensitySeaWater;
+
+    Kokkos::parallel_for("calcOceanStressCoeff", numVerticesOwned, KOKKOS_LAMBDA(const int vtx){
+      if(solve_velocity(vtx) == 0) return;
+      auto relVelSq = pow(oceanVelocity(vtx, 0) - velocity(vtx, 0), 2) +  pow(oceanVelocity(vtx, 1) - velocity(vtx, 1), 2);
+      oceanStressCoeff(vtx, 0) = configIceOceanDragCoeff * seaiceDensitySeaWater_ * iceAreaVtx(vtx, 0) * sqrt(relVelSq);
+    });
+  }
 
   void Mesh::gridSolveGPU(){
     //Mesh Fields
diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp
index a851f114..c0e3ffa4 100644
--- a/src/pmpo_mesh.hpp
+++ b/src/pmpo_mesh.hpp
@@ -19,6 +19,7 @@ enum MeshFieldIndex{
     MeshF_Unsupported,
     MeshF_VtxCoords,
     MeshF_VtxRotLat,
+    MeshF_VtxRotLon,
     MeshF_ElmCenterXYZ,
     MeshF_DualTriangleArea,
     MeshF_Vel,
@@ -41,7 +42,8 @@ enum MeshFieldIndex{
     MeshF_SurfaceTilt,
     MeshF_TotalMassFVtx,
     MeshF_OceanStress,
-    MeshF_OceanStressCoeff
+    MeshF_OceanStressCoeff,
+    MeshF_OceanVelocity
 };
 
 enum MeshFieldType{
@@ -54,6 +56,7 @@ enum MeshFieldType{
 template <MeshFieldIndex> struct meshFieldToType;
 template <> struct meshFieldToType < MeshF_VtxCoords         > { using type = Kokkos::View<vec3d_t*>; };
 template <> struct meshFieldToType < MeshF_VtxRotLat         > { using type = DoubleView; };
+template <> struct meshFieldToType < MeshF_VtxRotLon         > { using type = DoubleView; };
 template <> struct meshFieldToType < MeshF_ElmCenterXYZ      > { using type = Kokkos::View<vec3d_t*>; };
 template <> struct meshFieldToType < MeshF_DualTriangleArea  > { using type = Kokkos::View<doubleSclr_t*>; };
 template <> struct meshFieldToType < MeshF_Vel               > { using type = Kokkos::View<vec2d_t*>; };
@@ -77,6 +80,7 @@ template <> struct meshFieldToType < MeshF_SurfaceTilt       > { using type = Ko
 template <> struct meshFieldToType < MeshF_TotalMassFVtx     > { using type = Kokkos::View<doubleSclr_t*>; };
 template <> struct meshFieldToType < MeshF_OceanStress       > { using type = Kokkos::View<vec2d_t*>; };
 template <> struct meshFieldToType < MeshF_OceanStressCoeff  > { using type = Kokkos::View<doubleSclr_t*>; };
+template <> struct meshFieldToType < MeshF_OceanVelocity     > { using type = Kokkos::View<vec2d_t*>; };
 
 template <MeshFieldIndex index>
 using MeshFView = typename meshFieldToType<index>::type;
@@ -86,6 +90,7 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType, std::string>> meshFields
         {MeshF_Unsupported,      {MeshFType_Unsupported,"MeshField_Unsupported"}},
         {MeshF_VtxCoords,        {MeshFType_VtxBased,"MeshField_VerticesCoords"}},
         {MeshF_VtxRotLat,        {MeshFType_VtxBased,"MeshField_VerticesLatitude"}},
+        {MeshF_VtxRotLon,        {MeshFType_VtxBased,"MeshField_VerticesLongitude"}},
         {MeshF_ElmCenterXYZ,     {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}},
         {MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}},
         {MeshF_Vel,              {MeshFType_VtxBased,"MeshField_Velocity"}},
@@ -108,7 +113,8 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType, std::string>> meshFields
         {MeshF_SurfaceTilt,      {MeshFType_VtxBased,"MeshField_SurfaceTilt"}},
         {MeshF_TotalMassFVtx,    {MeshFType_VtxBased,"MeshField_TotalMassFVtx"}},
         {MeshF_OceanStress,      {MeshFType_VtxBased,"MeshField_OceanStress"}},
-        {MeshF_OceanStressCoeff, {MeshFType_VtxBased,"MeshField_OceanStressCoeff"}}
+        {MeshF_OceanStressCoeff, {MeshFType_VtxBased,"MeshField_OceanStressCoeff"}},
+        {MeshF_OceanVelocity,    {MeshFType_VtxBased,"MeshField_OceanVelocity"}}
 };
 
 enum mesh_type {mesh_unrecognized_lower = -1,
@@ -141,6 +147,7 @@ class Mesh {
     //start of meshFields
     MeshFView<MeshF_VtxCoords> vtxCoords_;
     MeshFView<MeshF_VtxRotLat> vtxRotLat_;
+    MeshFView<MeshF_VtxRotLon> vtxRotLon_;
     MeshFView<MeshF_ElmCenterXYZ> elmCenterXYZ_;
     MeshFView<MeshF_DualTriangleArea> dualTriangleArea_;
 
@@ -156,7 +163,7 @@ class Mesh {
     //GnomonicProjection
     MeshFView<MeshF_VtxGnomProj> vtxGnomProj_;
     MeshFView<MeshF_ElmCenterGnomProj> elmCenterGnomProj_;
-    
+
     MeshFView<MeshF_TanLatVertexRotatedOverRadius> tanLatVertexRotatedOverRadius_;
     MeshFView<MeshF_SolveStress> solveStress_;
     MeshFView<MeshF_SolveVelocity> solveVelocity_;
@@ -167,6 +174,7 @@ class Mesh {
     MeshFView<MeshF_TotalMassFVtx> totalMassFVtx_;
     MeshFView<MeshF_OceanStress> oceanStress_;
     MeshFView<MeshF_OceanStressCoeff> oceanStressCoeff_;
+    MeshFView<MeshF_OceanVelocity> oceanVelocity_;
 
     bool isRotatedFlag = false;
     double elasticTimeStep_;
@@ -271,6 +279,8 @@ class Mesh {
     double getDynamicTimeStep(){
       return dynamicTimeStep_;
     }
+
+    void calcOceanStressCoeff(const double configIceOceanDragCoeff);
     void gridSolveGPU();
     void aggregateDeluDyn();
     void applyFreeSlipBC();
@@ -292,6 +302,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_VtxRotLat){
         return vtxRotLat_;
     }
+    else if constexpr (index==MeshF_VtxRotLon){
+        return vtxRotLon_;
+    }
     else if constexpr (index==MeshF_ElmCenterXYZ){
         return elmCenterXYZ_;
     }
@@ -361,6 +374,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_OceanStressCoeff){
         return oceanStressCoeff_;
     }
+    else if constexpr (index==MeshF_OceanVelocity){
+        return oceanVelocity_;
+    }
     fprintf(stderr,"Mesh Field Index error!\n");
     exit(1);
 }
diff --git a/src/pmpo_utils.hpp b/src/pmpo_utils.hpp
index 84fa98e8..c7ff689d 100644
--- a/src/pmpo_utils.hpp
+++ b/src/pmpo_utils.hpp
@@ -45,6 +45,9 @@ using DoubleView = Kokkos::View<double*>;
 using IntView = Kokkos::View<int*>;
 using BoolView = Kokkos::View<bool*>;
 
+//CONSTANTS
+inline constexpr double seaiceDensitySeaWater = 1026.0;
+
 class Vec2d {
   private:
     vec2d_t coords_;
@@ -554,6 +557,24 @@ void lat_lon_from_xyz(double& lat, double& lon, Vec3d& xyz, double r){
   lat = Kokkos::asin(xyz[2]/r);
 }
 
+KOKKOS_INLINE_FUNCTION
+double seaice_mpm_wrap_longitude(const double longitude){
+  const auto PI = 3.141592653589;
+  const auto TAU = 2*PI;
+  const double wrapped = longitude - Kokkos::floor((longitude + PI) / TAU) * TAU;
+  return wrapped;
+}
+
+KOKKOS_INLINE_FUNCTION
+void seaice_mpm_coord_parallel_transport(const double long_start, const double long_end,
+                                         const double latitude, double transport[2]){
+  double dLon = long_end - long_start;
+  dLon = seaice_mpm_wrap_longitude(dLon);
+  const double psi = dLon * Kokkos::sin(latitude);
+  transport[0] = Kokkos::cos(psi);
+  transport[1] = Kokkos::sin(psi);
+}
+
 }//namespace polyMPO end
 
 #endif
diff --git a/src/pmpo_wachspressBasis.hpp b/src/pmpo_wachspressBasis.hpp
index 051537df..fc37f7a4 100644
--- a/src/pmpo_wachspressBasis.hpp
+++ b/src/pmpo_wachspressBasis.hpp
@@ -15,10 +15,13 @@ void sphericalInterpolation(MPMesh& mpMesh){
   auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
   int numVtxs = p_mesh->getNumVertices();
   auto elm2VtxConn = p_mesh->getElm2VtxConn();
+  auto vtxRotLon = p_mesh->getMeshField<MeshF_VtxRotLon>();
 
   auto p_MPs = mpMesh.p_MPs;
   auto MPsPosition = p_MPs->getPositions();
   auto MPsBasis = p_MPs->getData<MPF_Basis_Vals>();
+  auto curPosRotLatLon = p_MPs->getData<MPF_Cur_Pos_Rot_Lat_Lon>();
+  auto MPsAppID = p_MPs->getData<MPF_MP_APP_ID>();
 
   constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
   auto mpField = p_MPs->getData<mpfIndex>();
@@ -26,14 +29,34 @@ void sphericalInterpolation(MPMesh& mpMesh){
   const int numEntries = mpSliceToNumEntries<mpfIndex>();
   auto meshField = p_mesh->getMeshField<meshFieldIndex>(); 
 
+  bool use_correction_term = false;
+  if constexpr (meshFieldIndex == MeshF_OnSurfVeloIncr) {
+    use_correction_term = true;
+  }
+
   auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
       int numVtx = elm2VtxConn(elm,0);
-      for(int entry=0; entry<numEntries; entry++){
-        double mpValue = 0.0;
-        for(int i=1; i<= numVtx; i++)
-          mpValue += meshField(elm2VtxConn(elm,i)-1,entry)*MPsBasis(mp,i-1);
-        mpField(mp,entry) = mpValue;
+
+      if(use_correction_term){
+        double transport[2] = {0.0};
+        double mpField_l[2] = {0.0};
+        for(int i=1; i<=numVtx; i++){
+          int vID = elm2VtxConn(elm,i)-1;
+          seaice_mpm_coord_parallel_transport(vtxRotLon(vID), curPosRotLatLon(mp, 1), curPosRotLatLon(mp, 0), transport);
+          mpField_l[0] += ( transport[0] * meshField(vID, 0) + transport[1] * meshField(vID, 1)) * MPsBasis(mp,i-1);
+          mpField_l[1] += (-transport[1] * meshField(vID, 0) + transport[0] * meshField(vID, 1)) * MPsBasis(mp,i-1);
+        }
+        mpField(mp, 0) = mpField_l[0];
+        mpField(mp, 1) = mpField_l[1];
+      }
+      else{
+        for(int entry=0; entry<numEntries; entry++){
+          double mpValue = 0.0;
+          for(int i=1; i<= numVtx; i++)
+            mpValue += meshField(elm2VtxConn(elm,i)-1,entry)*MPsBasis(mp,i-1);
+          mpField(mp,entry) = mpValue;
+        }
       }
     }
   };
diff --git a/test/testFortran.f90 b/test/testFortran.f90
index 2eff4d9f..5d8538dd 100644
--- a/test/testFortran.f90
+++ b/test/testFortran.f90
@@ -76,7 +76,7 @@ program main
         MParray(i,j) = (i-1)*numMPs + j
     end do
   end do
-  call polympo_setMPVel(mpMesh, numCompsVel, numMPs, c_loc(MParray))
+  call polympo_setMPVel(mpMesh, numCompsVel, numMPs, c_loc(MParray), 5)
   
   ! check MP Fields
   MParray = -1
@@ -134,7 +134,7 @@ program main
   call polympo_setMeshVtxVel(mpMesh, nverts, c_loc(xArray),c_loc(yArray))
   xArray = -1
   yArray = -1
-  call polympo_getMeshVtxVel(mpMesh, nverts, c_loc(xArray),c_loc(yArray))
+  call polympo_getMeshVtxVel(mpMesh, nverts, c_loc(xArray),c_loc(yArray), 4)
   do i = 1, nverts
     call assert((xArray(i) .eq. i+value1), "Assert MeshVel u-component Velocity Fail")
     call assert((yArray(i) .eq. value2-i), "Assert MeshVel v-component Velocity Fail")
diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90
index 4ac1c6c9..1f48729c 100644
--- a/test/testFortranMPAdvection.f90
+++ b/test/testFortranMPAdvection.f90
@@ -105,7 +105,7 @@ subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2
     mpVel = TEST_VAL
 
     call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass))
-    call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel))
+    call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel), 7)
 
     ! Although this test just does 0th order reconstruction testing, and just needs the BasisSlice,
     ! calculating the coefficeints too as that will involve calculating the Basis Slice
@@ -169,7 +169,7 @@ subroutine runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPo
       call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition))
       call polympo_setMeshVtxOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr))
       call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass))
-      call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel))
+      call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel), 8)
       call polympo_setMeshVtxCoords(mpMesh, nVertices, c_loc(xVertex), c_loc(yVertex), c_loc(zVertex))
       call polympo_setMeshVtxRotLat(mpMesh,nVertices,c_loc(latVertex))
 
@@ -177,7 +177,7 @@ subroutine runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPo
       call polympo_getMPVel(mpMesh, 2, numMPs, c_loc(mpVel))
       call polympo_getMeshElmMass(mpMesh,nCells,c_loc(meshElmMass))
       call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass))
-      call polympo_getMeshVtxVel(mpMesh,nVertices, c_loc(xVertex),c_loc(yVertex))
+      call polympo_getMeshVtxVel(mpMesh,nVertices, c_loc(xVertex),c_loc(yVertex), 6)
     end do
 
     deallocate(dispIncr)
diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90
index c98f1b88..86f88238 100644
--- a/test/testFortranMPReconstruction.f90
+++ b/test/testFortranMPReconstruction.f90
@@ -120,21 +120,21 @@ program main
   call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition))
 
   call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass))
-  call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel))
+  call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel), 6)
 
   !First Order reconstruction done before 0th order reconstruction as calculation of coefficients will
   !fill the MpBasis slice
   !Test vtx order 1 reconstruction
   call polympo_reconstruct_coeff_with_MPI(mpmesh)
   call polympo_setReconstructionOfMass(mpMesh,1,polympo_getMeshFVtxType())
-  call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType())
+  !call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType())
   call polympo_applyReconstruction(mpMesh)
   call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass1))
-  call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv))
+  call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv), 5)
   do i = 1, nVertices
     call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass order 1")
-    call assert(meshVtxVelu(i)  < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i)  > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1")
-    call assert(meshVtxVelv(i)  < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i)  > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1")
+    !call assert(meshVtxVelu(i)  < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i)  > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1")
+    !call assert(meshVtxVelv(i)  < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i)  > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1")
   end do
 
   ! Test vtx order 0  reconstruction