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..0ca24de9 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){
@@ -29,7 +31,7 @@ void MPMesh::calculateStrain(){
         MPsStrainRate(mp, 2) =  0.0;
         return;
       }
-
+     
       int numVtx = elm2VtxConn(elm,0);
 
       double v11 = 0.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");
@@ -171,7 +178,7 @@ void MPMesh::calculateStressDivergence(){
   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);
     //Transfer the correct values at owned vertices to halo vertices
     //communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 1, 1);
   }
diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp
index d3aac1b7..094289de 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;
+        int numToRecv = 0;
 
-        int numToSend = 0, numToRecv = 0; 
-        if(mode == 0) {
-          //gather (halos send to owners)
+        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,663 @@ 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){
+
+      int self, numProcsTot;
+
+      MPI_Comm comm = p_MPs->getMPIComm();
+
+      MPI_Comm_rank(comm, &self);
+      MPI_Comm_size(comm, &numProcsTot);
+
+      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;
+
+        pumipic::RecordTime(
+            "SD: CUDA-aware MPI Cache Build m" + std::to_string(mode) +
+                " e" + std::to_string(numEntries) + "-" + 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();
+
+      pumipic::RecordTime(
+          "SD: CUDA-aware MPI Pack m" + std::to_string(mode) +
+              " e" + std::to_string(numEntries) + "-" + std::to_string(self),
+          timer.seconds());
+
+      timer.reset();
+
+      Kokkos::Timer mpiTotalTimer;
+      std::vector<MPI_Request> requests;
+      requests.reserve(2 * numProcsTot);
+      int mpiError = MPI_SUCCESS;
+
+      for(int proc = 0; proc < numProcsTot; proc++){
+        if(proc == self) continue;
+
+        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);
+        }
+
+        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);
+        }
+      }
+
+      pumipic::RecordTime(
+          "SD: CUDA-aware MPI Post m" + std::to_string(mode) +
+              " e" + std::to_string(numEntries) + "-" + std::to_string(self),
+          timer.seconds());
+
+      timer.reset();
+
+      if(mpiError == MPI_SUCCESS && !requests.empty()){
+        mpiError = MPI_Waitall(
+            static_cast<int>(requests.size()),
+            requests.data(),
+            MPI_STATUSES_IGNORE);
+      }
+
+      pumipic::RecordTime(
+          "SD: CUDA-aware MPI Wait m" + std::to_string(mode) +
+              " e" + std::to_string(numEntries) + "-" + std::to_string(self),
+          timer.seconds());
+
+      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;
+        }
+
+        timer.reset();
+
+        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;
+      }
+
+      pumipic::RecordTime(
+          "SD: CUDA-aware MPI Comm m" + std::to_string(mode) +
+              " e" + std::to_string(numEntries) + "-" + std::to_string(self),
+          mpiTotalTimer.seconds());
+
+      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();
+
+      pumipic::RecordTime(
+          "SD: CUDA-aware MPI Contribution m" + std::to_string(mode) +
+              " e" + std::to_string(numEntries) + "-" + 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){
+
+      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..9685bbb5 100644
--- a/src/pmpo_MPMesh_assembly.hpp
+++ b/src/pmpo_MPMesh_assembly.hpp
@@ -130,6 +130,7 @@ void MPMesh::reconstruct_coeff_full(){
     radius=p_mesh->getSphereRadius();
 
   //Assemble matrix for each vertex
+  timer.reset();
   auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
       int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element
@@ -152,7 +153,7 @@ 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());
+  pumipic::RecordTime("VR Assemble Matrix Per Process" + std::to_string(self), timer.seconds());
   //Mode 0 is Gather:  Halos Send to Owners
   //Mode 1 is Scatter: Owners Send to Halos
   //Op 0 is addition
@@ -161,21 +162,27 @@ void MPMesh::reconstruct_coeff_full(){
   int mode = 0;
   int op = 0;
   if (numProcsTot >1){
-    communicate_and_take_halo_contributions1(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
+    communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
+    pumipic::RecordTime("VR Matrix Gather Halo/MPI " + std::to_string(self), timer.seconds());
+
+    timer.reset();
     mode=1; 
     op=1;
-    communicate_and_take_halo_contributions1(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
+    communicate_and_take_halo_contributions1_improved(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
+    pumipic::RecordTime("VR Matrix Scatter Halo/MPI " + std::to_string(self), timer.seconds());
   }
-  pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds());
- 
-  //Stroe the 1st matrix element
+  
+  //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;
 
+  timer.reset();
   invertMatrix(vtxMatrices, radius);
+  pumipic::RecordTime("VR Invert Matrix " + std::to_string(self), timer.seconds());
 }
 
 void MPMesh::invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius){
@@ -380,7 +387,7 @@ void MPMesh::assemblyVtx1(){
 
   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);
   }
   pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds());
 }
diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp
index 53c330e1..73528f2f 100644
--- a/src/pmpo_c.cpp
+++ b/src/pmpo_c.cpp
@@ -728,7 +728,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 +753,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 +806,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);
@@ -1448,6 +1538,20 @@ 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);
@@ -1570,6 +1674,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){
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  p_mesh->calcOceanStressCoeff();
+}
+
 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 +1714,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); 
 }
 
 //Advection Calcualtions
diff --git a/src/pmpo_c.h b/src/pmpo_c.h
index 23d78b70..99ba04f9 100644
--- a/src/pmpo_c.h
+++ b/src/pmpo_c.h
@@ -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
@@ -106,6 +110,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 +119,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);
 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..522f3b22 100644
--- a/src/pmpo_fortran.f90
+++ b/src/pmpo_fortran.f90
@@ -419,6 +419,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 +434,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
@@ -1011,6 +1041,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
@@ -1073,6 +1111,12 @@ subroutine polympo_set_oceanStressCoefficient(mpMesh, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
 
+  subroutine polympo_calculate_oceanStressCoefficient(mpMesh) &
+             bind(C, NAME='polympo_calculate_oceanStressCoefficient_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end subroutine
+
   subroutine  polympo_velocity_grid_solve(mpMesh) &
               bind(C, NAME='polympo_velocity_grid_solve_f')
     use :: iso_c_binding
diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp
index 0e7849ca..5acb9454 100644
--- a/src/pmpo_mesh.cpp
+++ b/src/pmpo_mesh.cpp
@@ -78,6 +78,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 +160,20 @@ namespace polyMPO{
     });
   }
 
+  void Mesh::calcOceanStressCoeff(){
+    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>();
+
+    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) = 0.00536 * 1026.0 * iceAreaVtx(vtx, 0) * sqrt(relVelSq);
+    });
+  }
 
   void Mesh::gridSolveGPU(){
     //Mesh Fields
diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp
index a851f114..40b7e03a 100644
--- a/src/pmpo_mesh.hpp
+++ b/src/pmpo_mesh.hpp
@@ -41,7 +41,8 @@ enum MeshFieldIndex{
     MeshF_SurfaceTilt,
     MeshF_TotalMassFVtx,
     MeshF_OceanStress,
-    MeshF_OceanStressCoeff
+    MeshF_OceanStressCoeff,
+    MeshF_OceanVelocity
 };
 
 enum MeshFieldType{
@@ -77,6 +78,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;
@@ -108,7 +110,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,
@@ -167,6 +170,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 +275,8 @@ class Mesh {
     double getDynamicTimeStep(){
       return dynamicTimeStep_;
     }
+
+    void calcOceanStressCoeff(); 
     void gridSolveGPU();
     void aggregateDeluDyn();
     void applyFreeSlipBC();
@@ -361,6 +367,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);
 }