diff --git a/src/core_landice/Registry.xml b/src/core_landice/Registry.xml
index 7062959278..c534626142 100644
--- a/src/core_landice/Registry.xml
+++ b/src/core_landice/Registry.xml
@@ -931,7 +931,7 @@ is the value of that variable from the *previous* time level!
                 <var name="dynamicThickening" type="real" dimensions="nCells Time" units="m a^{-1}"
                      description="diagnostic field of dynamic thickening rate (calculated as negative of flux divergence)"
                 />
-                <var name="cellMask" type="integer" dimensions="nCells Time" units="none"
+                <var name="cellMask" type="integer" dimensions="nCells Time" units="none" time_levs="1"
                      description="bitmask indicating various properties about the ice sheet on cells.  cellMask only needs to be a restart field if config_allow_additional_advance = false (to keep the mask of initial ice extent)"
                 />
                 <var name="edgeMask" type="integer" dimensions="nEdges Time" units="none"
@@ -1127,6 +1127,10 @@ is the value of that variable from the *previous* time level!
                      units="none"  description="flag needed by external velocity solvers that indicates if the region to solve on the block's domain has changed (treated as a logical)"
                         packages="higherOrderVelocity"
                 />
+                <var name="cellIceMaskChanged"              type="integer"  dimensions="Time"
+                     units="none"  description="flag needed by external velocity solvers that indicates if the region to solve on the block's domain has changed (treated as a logical)"
+                        packages="higherOrderVelocity"
+                />
                 <var name="dirichletVelocityMask"            type="integer"     dimensions="nVertInterfaces nCells Time"
                         units="none"  time_levs="2"
                         description="mask of where Dirichlet boundary conditions should be applied to the velocity solution.  1 means apply a Dirichlet boundary condition, 0 means do not. (higher-order dycores only)"
@@ -1139,9 +1143,14 @@ is the value of that variable from the *previous* time level!
                 <var name="floatingEdges"
                         type="integer" dimensions="nEdges Time"
                         units="unitless"
+                        time_levs="2"
                         description="edges which are floating have a value of 1.  non floating edges have a value of 0."
                         packages="higherOrderVelocity"
                 />
+                <var name="floatingEdgesChanged"              type="integer"  dimensions="Time"
+                     units="none"  description="flag needed by external velocity solvers that indicates if the floating edges ids have changed (treated as a logical)"
+                        packages="higherOrderVelocity"
+                />
         </var_struct>
 
 <!-- ================ -->
diff --git a/src/core_landice/analysis_members/mpas_li_global_stats.F b/src/core_landice/analysis_members/mpas_li_global_stats.F
index 41ddc9aebd..b367875979 100644
--- a/src/core_landice/analysis_members/mpas_li_global_stats.F
+++ b/src/core_landice/analysis_members/mpas_li_global_stats.F
@@ -283,7 +283,7 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err)
          call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1)
          call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied)
          call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied)
          call mpas_pool_get_array(geometryPool, 'basalMassBal', basalMassBal)
diff --git a/src/core_landice/analysis_members/mpas_li_regional_stats.F b/src/core_landice/analysis_members/mpas_li_regional_stats.F
index 0e0543df57..4c6df7def9 100644
--- a/src/core_landice/analysis_members/mpas_li_regional_stats.F
+++ b/src/core_landice/analysis_members/mpas_li_regional_stats.F
@@ -275,7 +275,7 @@ subroutine li_compute_regional_stats(domain, memberName, timeLevel, err)
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1)
          call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
          call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied)
          call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied)
diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp
index 6612b7bba0..671ed32bcd 100644
--- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp
+++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp
@@ -14,31 +14,29 @@ distributed with this code, or at http://mpas-dev.github.com/license.html
 #include <algorithm>
 #include <set>
 #include <fstream>
+#include <sstream>
+#include <limits>
 #include "Interface_velocity_solver.hpp"
-//#include <lifev/ice_sheet/interface_with_mpas/Interface.hpp>
-//#include <lifev/ice_sheet/solver/BuildMeshFromBareData.hpp>
-//#include <lifev/ice_sheet/solver/Extrude3DMesh.hpp>
 
 // ===================================================
 //! Namespaces
 // ===================================================
 
-//typedef std::list<exchange> exchangeList_Type;
+#define changeTrianglesOwnership
+
 
 // ice_problem pointer
 
 int Ordering = 0; //ordering ==0 means that the mesh is extruded layerwise, whereas ordering==1 means that the mesh is extruded columnwise.
 MPI_Comm comm, reducedComm;
 bool isDomainEmpty = true;
-bool initialize_velocity = true;
 bool first_time_step = true;
 int nCells_F, nEdges_F, nVertices_F;
 int nCellsSolve_F, nEdgesSolve_F, nVerticesSolve_F;
-int nVertices, nEdges, nTriangles, nGlobalVertices, nGlobalEdges,
-    nGlobalTriangles;
+int nVertices, nEdges, nTriangles, globalVertexStride, globalEdgeStride,globalTriangleStride;
 int maxNEdgesOnCell_F;
 int const *cellsOnEdge_F, *cellsOnVertex_F, *verticesOnCell_F,
-    *verticesOnEdge_F, *edgesOnCell_F, *indexToCellID_F, *nEdgesOnCells_F,
+    *verticesOnEdge_F, *edgesOnCell_F, *indexToCellID_F, *indexToEdgeID_F, *indexToVertexID_F, *nEdgesOnCells_F,
     *verticesMask_F, *cellsMask_F, *dirichletCellsMask_F, *floatingEdgesMask_F;
 std::vector<double> layersRatio, levelsNormalizedThickness;
 int nLayers;
@@ -63,13 +61,13 @@ int original_stderr; // the location of stderr before we captured it
 int Interface_stdout; // the location of stdout as we use it here
 
 //void *phgGrid = 0;
-std::vector<int> edgesToReceive, fCellsToReceive, indexToTriangleID,
-    verticesOnTria, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge,
-    trianglesProcIds,  reduced_ranks;
-std::vector<int> indexToVertexID, vertexToFCell, triangleToFVertex, indexToEdgeID, edgeToFEdge,
-    mask, fVertexToTriangleID, fCellToVertex, floatingEdgesIds, dirichletNodesIDs;
-std::vector<double> temperatureOnTetra, dissipationHeatOnTetra, velocityOnVertices, velocityOnCells,
-    elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, temperatureData, smbData, thicknessOnCells;
+std::vector<int> indexToTriangleID,
+    verticesOnTria, trianglesOnEdge, verticesOnEdge,
+    trianglesProcIds, reduced_ranks;
+std::vector<int> indexToVertexID, vertexToFCell, vertexProcIDs, triangleToFVertex, indexToEdgeID, edgeToFEdge,
+    fVertexToTriangleID, fCellToVertex, floatingEdgesIds, dirichletNodesIDs;
+std::vector<double> dissipationHeatOnPrisms, velocityOnVertices, velocityOnCells,
+    elevationData, thicknessData, betaData, bedTopographyData, stiffnessFactorData, effecPressData, temperatureDataOnPrisms, smbData, thicknessOnCells;
 std::vector<bool> isVertexBoundary, isBoundaryEdge;
 
 // only needed for creating ASCII mesh
@@ -86,8 +84,8 @@ double radius;
 exchangeList_Type const *sendCellsList_F = 0, *recvCellsList_F = 0;
 exchangeList_Type const *sendEdgesList_F = 0, *recvEdgesList_F = 0;
 exchangeList_Type const *sendVerticesList_F = 0, *recvVerticesList_F = 0;
-exchangeList_Type sendCellsListReversed, recvCellsListReversed,
-    sendEdgesListReversed, recvEdgesListReversed;
+exchangeList_Type sendVerticesListReversed, recvVerticesListReversed,
+    sendCellsListReversed, recvCellsListReversed;
 
 exchange::exchange(int _procID, int const* vec_first, int const* vec_last,
     int fieldDim) :
@@ -106,7 +104,7 @@ int velocity_solver_init_mpi(int* fComm) {
   // get MPI_Comm from Fortran
   comm = MPI_Comm_f2c(*fComm);
   reducedComm = MPI_COMM_NULL;  // initialize to null so we can check if set
-  return velocity_solver_init_mpi__(fComm);
+  return velocity_solver_init_mpi__(comm);
 }
 
 
@@ -131,11 +129,6 @@ void velocity_solver_export_2d_data(double const* lowerSurface_F,
     double const* thickness_F, double const* beta_F) {
   if (isDomainEmpty)
     return;
-#ifdef LIFEV
-  import2DFields(lowerSurface_F, thickness_F, beta_F,  minThickneess);
-  velocity_solver_export_2d_data__(reducedComm, elevationData, thicknessData,
-      betaData, indexToVertexID);
-#endif
 }
 
 void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
@@ -146,6 +139,8 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
     int const* _verticesOnCell_F, int const* _verticesOnEdge_F,
     int const* _edgesOnCell_F, int const* _nEdgesOnCells_F,
     int const* _indexToCellID_F,
+    int const* _indexToEdgeID_F,
+    int const* _indexToVertexID_F,
     double const* _xCell_F, double const* _yCell_F, double const* _zCell_F,
     double const* _xVertex_F, double const* _yVertex_F, double const* _zVertex_F,
     double const* _areaTriangle_F,
@@ -168,7 +163,9 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
   verticesOnEdge_F = _verticesOnEdge_F;
   edgesOnCell_F = _edgesOnCell_F;
   nEdgesOnCells_F = _nEdgesOnCells_F;
+  indexToEdgeID_F = _indexToEdgeID_F;
   indexToCellID_F = _indexToCellID_F;
+  indexToVertexID_F = _indexToVertexID_F;
   xCell_F = _xCell_F;
   yCell_F = _yCell_F;
   zCell_F = _zCell_F;
@@ -176,7 +173,6 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
   yVertex_F = _yVertex_F;
   zVertex_F = _zVertex_F;
   areaTriangle_F = _areaTriangle_F;
-  mask.resize(nVertices_F);
 
   thicknessOnCells.resize(nCellsSolve_F);
 
@@ -189,9 +185,6 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
   recvVerticesList_F = new exchangeList_Type(
       unpackMpiArray(recvVerticesArray_F));
 
-  trianglesProcIds.resize(nVertices_F);
-  getProcIds(trianglesProcIds, recvVerticesList_F);
-
   if (radius > 10) {
     xCellProjected.resize(nCells_F);
     yCellProjected.resize(nCells_F);
@@ -210,23 +203,6 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
 }
 
 void velocity_solver_init_l1l2(double const* levelsRatio_F) {
-#ifdef LIFEV
-  velocityOnVertices.resize(2 * nVertices * (nLayers + 1), 0.);
-  velocityOnCells.resize(2 * nCells_F * (nLayers + 1), 0.);
-
-  if (isDomainEmpty)
-    return;
-
-  layersRatio.resize(nLayers);
-  // !!Indexing of layers is reversed
-  for (int i = 0; i < nLayers; i++)
-    layersRatio[i] = levelsRatio_F[nLayers - 1 - i];
-  //std::copy(levelsRatio_F, levelsRatio_F+nLayers, layersRatio.begin());
-  mapCellsToVertices(velocityOnCells, velocityOnVertices, 2, nLayers, Ordering);
-
-  velocity_solver_init_l1l2__(layersRatio, velocityOnVertices, initialize_velocity);
-  initialize_velocity = false;
-#endif
 }
 
 
@@ -236,82 +212,9 @@ void velocity_solver_solve_l1l2(double const* lowerSurface_F,
     double const* thickness_F, double const* beta_F, double const* temperature_F,
     double* const dirichletVelocityXValue, double* const dirichletVelocitYValue,
     double* u_normal_F, double* xVelocityOnCell, double* yVelocityOnCell) {
-
-#ifdef LIFEV
-
-  std::fill(u_normal_F, u_normal_F + nEdges_F * (nLayers+1), 0.);
-
-  double localSum(0), sum(0);
-
-  for (int i = 0; i < nCellsSolve_F; i++) {
-    localSum = std::max(localSum,
-        std::fabs(thickness_F[i] - thicknessOnCells[i]));
-  }
-
-  MPI_Allreduce(&localSum, &sum, 1, MPI_DOUBLE, MPI_MAX, comm);
-
-  std::cout << "Thickness change: " << sum << std::endl;
-  std::copy(thickness_F, thickness_F + nCellsSolve_F, &thicknessOnCells[0]);
-
-  if (!isDomainEmpty) {
-    std::vector<double> temperatureData(nLayers * nVertices);
-
-    import2DFields(lowerSurface_F, thickness_F, beta_F,  minThickness);
-
-    for (int index = 0; index < nVertices; index++) {
-      int iCell = vertexToFCell[index];
-      for (int il = 0; il < nLayers; il++) {
-        temperatureData[index + il * nVertices] = temperature_F[iCell * nLayers
-            + (nLayers - il - 1)];
-      }
-    }
-
-
-    velocity_solver_solve_l1l2__(elevationData, thicknessData, betaData,
-        temperatureData, indexToVertexID, velocityOnVertices);
-  }
-
-
-   mapVerticesToCells (velocityOnVertices, &velocityOnCells[0], 2, nLayers, Ordering);
-
-   //computing x, yVelocityOnCell
-   int sizeVelOnCell = nCells_F * (nLayers + 1);
-   for(int iCell=0; iCell<nCells_F; ++iCell)
-     for(int il=0; il<nLayers + 1; ++il) {
-       int ilReversed = nLayers - il;
-       int indexReversed = iCell * (nLayers+1) + ilReversed;
-       int index = iCell * (nLayers + 1) +il;
-     xVelocityOnCell[indexReversed] = velocityOnCells[index];
-     yVelocityOnCell[indexReversed] = velocityOnCells[index+sizeVelOnCell];
-   }
-
-   if (!isDomainEmpty)
-     get_prism_velocity_on_FEdges(u_normal_F, velocityOnCells, edgeToFEdge);
-
-   allToAll (u_normal_F,  &sendEdgesListReversed, &recvEdgesListReversed, nLayers+1);
-   allToAll (u_normal_F,  sendEdgesList_F, recvEdgesList_F, nLayers+1);
-
-#endif
-
 }
 
-void velocity_solver_export_l1l2_velocity() {
-
-  if (isDomainEmpty)
-     return;
-
-  std::vector<double> regulThk(thicknessData);
-      for (int index = 0; index < nVertices; index++)
-        regulThk[index] = std::max(1e-4, thicknessData[index]);
-
-      std::vector<int> mpasIndexToVertexID(nVertices);
-      for (int i = 0; i < nVertices; i++) {
-        mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]];
-      }
-#ifdef LIFEV
-  velocity_solver_export_l1l2_velocity__(layersRatio, elevationData, regulThk, mpasIndexToVertexID, reducedComm);
-#endif
-}
+void velocity_solver_export_l1l2_velocity() {}
 
 void velocity_solver_init_fo(double const *levelsRatio_F) {
 
@@ -325,16 +228,8 @@ void velocity_solver_init_fo(double const *levelsRatio_F) {
   // !!Indexing of layers is reversed
   for (int i = 0; i < nLayers; i++)
     layersRatio[i] = levelsRatio_F[nLayers - 1 - i];
-  //std::copy(levelsRatio_F, levelsRatio_F+nLayers, layersRatio.begin());
-
 
   mapCellsToVertices(velocityOnCells, velocityOnVertices, 2, nLayers, Ordering);
-
-#ifdef LIFEV
-  velocity_solver_init_fo__(layersRatio, velocityOnVertices, indexToVertexID, initialize_velocity);
-#endif
-  //    iceProblemPtr->initializeSolverFO(layersRatio, velocityOnVertices, thicknessData, elevationData, indexToVertexID, initialize_velocity);
-  initialize_velocity = false;
 }
 
 void velocity_solver_solve_fo(double const* bedTopography_F, double const* lowerSurface_F,
@@ -368,53 +263,30 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower
 
   if (!isDomainEmpty) {
 
-#ifdef LIFEV
-   double localSum(0), sum(0);
-
-   for (int i = 0; i < nCellsSolve_F; i++) {
-   localSum = std::max(localSum,
-   std::fabs(thickness_F[i] - thicknessOnCells[i]));
-   }
-
-   MPI_Allreduce(&localSum, &sum, 1, MPI_DOUBLE, MPI_MAX, comm);
-
-   std::cout << "Thickness change: " << sum << std::endl;
-   std::copy(thickness_F, thickness_F + nCellsSolve_F, &thicknessOnCells[0]);
-#endif
-
-
-
     std::map<int, int> bdExtensionMap;
-    import2DFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, temperature_F, smb_F,  minThickness);
+    importFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, stiffnessFactor_F, effecPress_F, temperature_F, smb_F,  minThickness);
 
     std::vector<double> regulThk(thicknessData);
     for (int index = 0; index < nVertices; index++)
       regulThk[index] = std::max(1e-4, thicknessData[index]);
 
-    importP0Temperature();
-
-    dissipationHeatOnTetra.resize(3 * nLayers * indexToTriangleID.size());
+    dissipationHeatOnPrisms.resize(nLayers * indexToTriangleID.size());
 
     std::cout << "\n\nTimeStep: "<< *deltat << "\n\n"<< std::endl;
 
     double dt = (*deltat)/secondsInAYear;
     int albany_error;
-    velocity_solver_solve_fo__(nLayers, nGlobalVertices, nGlobalTriangles,
+    velocity_solver_solve_fo__(nLayers, globalVertexStride, globalTriangleStride,
         Ordering, first_time_step, indexToVertexID, indexToTriangleID, minBeta,
         regulThk, levelsNormalizedThickness, elevationData, thicknessData,
         betaData, bedTopographyData, smbData,
         stiffnessFactorData, effecPressData,
-        temperatureOnTetra, dissipationHeatOnTetra, velocityOnVertices,
+        temperatureDataOnPrisms, dissipationHeatOnPrisms, velocityOnVertices,
         albany_error, dt);
     *error=albany_error;
   }
-  
-  exportDissipationHeat(dissipation_heat_F);
 
-  std::vector<int> mpasIndexToVertexID(nVertices);
-  for (int i = 0; i < nVertices; i++) {
-      mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]];
-  }
+  exportDissipationHeat(dissipation_heat_F);
 
   mapVerticesToCells(velocityOnVertices, &velocityOnCells[0], 2, nLayers,
       Ordering);
@@ -428,81 +300,15 @@ void velocity_solver_solve_fo(double const* bedTopography_F, double const* lower
     xVelocityOnCell[indexReversed] = velocityOnCells[index];
     yVelocityOnCell[indexReversed] = velocityOnCells[index+sizeVelOnCell];
   }
-
-  if (!isDomainEmpty)
     get_prism_velocity_on_FEdges(u_normal_F, velocityOnCells, edgeToFEdge);
-
+    allToAll(u_normal_F, sendEdgesList_F, recvEdgesList_F, nLayers+1);
   std::vector<double> velOnEdges(nEdges * (nLayers+1));
   for (int i = 0; i < nEdges; i++) {
     for (int il = 0; il < nLayers+1; il++) {
       velOnEdges[i * (nLayers+1) + il] = u_normal_F[edgeToFEdge[i] * (nLayers+1) + il];
     }
   }
-
-  allToAll(u_normal_F, &sendEdgesListReversed, &recvEdgesListReversed, nLayers+1);
-
-  allToAll(u_normal_F, sendEdgesList_F, recvEdgesList_F, nLayers+1);
-
   first_time_step = false;
-
-#ifdef LIFEV
-
-  std::vector<int> edgesProcId(nEdges_F), trianglesProcIds(nVertices_F);
-  getProcIds(edgesProcId, recvEdgesList_F);
-  getProcIds(trianglesProcIds, recvVerticesList_F);
-
-  int localSumInt(0), sumInt(0);
-
-  for (int i = 0; i < nEdges; i++) {
-    for (int il = 0; il < 1; il++) {
-      if (std::fabs(
-          velOnEdges[i * (nLayers+1) + il]
-              - u_normal_F[edgeToFEdge[i] * nLayers + il]) > 1e-9)
-      //  if(edgeToFEdge[i]>nEdgesSolve_F)
-          {
-        localSumInt++;
-        int edge = edgeToFEdge[i];
-        int gEdge = indexToEdgeID[i];
-        ID fVertex0 = verticesOnEdge_F[2 * edge] - 1;
-        ID fVertex1 = verticesOnEdge_F[2 * edge + 1] - 1;
-        ID triaId0 = fVertexToTriangleID[fVertex0];
-        ID triaId1 = fVertexToTriangleID[fVertex1];
-        ID procTria0 = trianglesProcIds[fVertex0];
-        ID procTria1 = trianglesProcIds[fVertex1];
-        std::cout << "vs( " << velOnEdges[i * (nLayers+1) + il] << ", "
-            << u_normal_F[edgeToFEdge[i] * nLayers + il] << ")  ";
-        std::cout << "edge: " << edge << ", gEdge: " << gEdge << ", on proc: "
-            << edgesProcId[edgeToFEdge[i]];
-        if (triaId0 != NotAnId) {
-          std::cout << ". first tria0: " << triaId0 << " on proc: "
-              << procTria0;
-        }
-        if (triaId1 != NotAnId) {
-          std::cout << ".. second tria0:" << std::endl;
-        }
-        if ((triaId0 == NotAnId) || (triaId1 == NotAnId)) {
-          std::cout << ". and to Tria: " << triaId1 << " on proc: " << procTria1
-              << std::endl;
-        }
-
-      }
-
-      //localSum = std::max(localSum, std::fabs(velOnEdges[i*nLayers+il] - u_normal_F[edgeToFEdge[i]*nLayers+il]));
-    }
-  }
-
-  MPI_Allreduce(&localSumInt, &sumInt, 1, MPI_INT, MPI_SUM, comm);
-
-  int localNum(sendEdgesListReversed.size()), num(0);
-
-  MPI_Allreduce(&localNum, &num, 1, MPI_INT, MPI_SUM, comm);
-
-  std::cout << "Edges change: " << sumInt << " " << num << std::endl;
-
-#endif
-
-
-
 }
 
 
@@ -526,7 +332,7 @@ void velocity_solver_finalize() {
 
 /*duality:
  *
- *   mpas(F) |  lifev
+ *   mpas(F) |  Albany LandIce (C++)
  *  ---------|---------
  *   cell    |  vertex
  *   vertex  |  triangle
@@ -536,7 +342,6 @@ void velocity_solver_finalize() {
 
 void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cellsMask_F, int const* _dirichletCellsMask_F, int const* _floatingEdgesMask_F) {
   int numProcs, me;
-
   verticesMask_F = _verticesMask_F;
   cellsMask_F = _cellsMask_F;
   verticesMask_F = _verticesMask_F;
@@ -549,88 +354,156 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel
       numProcs + 1), globalOffsetVertices(numProcs + 1), globalOffsetEdge(
       numProcs + 1);
 
+
+  // First, we compute the FE triangles belonging to this processor.
+  // If changeTrianglesOwnership is not define, the triangles belonging to this
+  // processor will be the subset of the triangles (MPAS vertices) owned by this proc
+  // that contain dynamic ice.
+  // If changeTrianglesOwnership is define, we rearrange the ownership of the triangles
+  // to improve the quality of the FE mesh and avoid corner cases (see below).
+
   triangleToFVertex.clear();
   triangleToFVertex.reserve(nVertices_F);
   std::vector<int> fVertexToTriangle(nVertices_F, NotAnId);
-  bool changed = false;
-  for (int i(0); i < nVerticesSolve_F; i++) {
-    if ((verticesMask_F[i] & dynamic_ice_bit_value) && !isGhostTriangle(i)) {
-      fVertexToTriangle[i] = triangleToFVertex.size();
-      triangleToFVertex.push_back(i);
-    }
-    changed = changed || (verticesMask_F[i] != mask[i]);
-  }
 
-  for (int i(0); i < nVertices_F; i++)
-    mask[i] = verticesMask_F[i];
+  //vector containing proc ranks for owned and shared FE triangles
+  trianglesProcIds.assign(nVertices_F,NotAnId);
+
+  //vector containing proc ranks for owned and shared MPAS cells
+  std::vector<int> fCellsProcIds(nCells_F);
+  getProcIds(fCellsProcIds, recvCellsList_F);
+#ifdef changeTrianglesOwnership
+  std::vector<int> fVerticesProcIds(nVertices_F);
+  getProcIds(fVerticesProcIds, recvVerticesList_F);
+  for (int i(0); i < nVertices_F; i++) {
+    int cellWithMinID=nCellsSolve_F;
+    if ((verticesMask_F[i] & dynamic_ice_bit_value)) {
+      int minCellId = std::numeric_limits<int>::max();
+      int minCellIdProc(0);
+
+      int cellProc[3];
+      bool invalidCell=false;
+      for (int j = 0; j < 3; j++) {
+        int iCell = cellsOnVertex_F[3 * i + j] - 1;
+        if(iCell >= nCells_F) {
+          invalidCell = true;
+          break;
+        }
+        int cellID = indexToCellID_F[iCell];
+        cellProc[j] = fCellsProcIds[iCell];
+        if(cellID < minCellId) {
+          minCellId = cellID;
+          cellWithMinID = iCell;
+          minCellIdProc = cellProc[j];
+        }
+      }
+
+      if(invalidCell) continue;
+
+      // the proc that owns at least 2 nodes of the triangle i. If all nodes belong to different procs, procOwns2Nodes is set to -1
+      int procOwns2Nodes = ((cellProc[0] ==  cellProc[1]) || (cellProc[0] ==  cellProc[2])) ? cellProc[0] :
+                           (cellProc[1] == cellProc[2]) ? cellProc[1] : -1;
 
-  if (changed)
-    std::cout << "mask changed!!" << std::endl;
+      int vertexProc = fVerticesProcIds[i];
+      bool triangleOwnsANode = (cellProc[0] == vertexProc) || (cellProc[1] == vertexProc) || (cellProc[2] == vertexProc);
 
-  if ((me == 0) && (triangleToFVertex.size() == 0))
-    for (int i(0); i < nVerticesSolve_F; i++) {
-      if (!isGhostTriangle(i)) {
+      //A triangle will be owned by a proc if:
+      // 1. the proc owns at least 2 nodes of the triangle associated to that vertex, OR
+      // 2. all the nodes of the triangle belong to three different procs, and the proc owns the fortran vertex  and a node OR
+      // 3. the three nodes of the triangle and the fortran vertex belong to four different procs, and the proc owns the node with the minimum ID
+
+      trianglesProcIds[i] = (procOwns2Nodes != -1) ? procOwns2Nodes :
+                       triangleOwnsANode ? vertexProc :
+                       minCellIdProc;
+
+      if (trianglesProcIds[i] == me) {
         fVertexToTriangle[i] = triangleToFVertex.size();
         triangleToFVertex.push_back(i);
-        break;
-      }
+      } 
     }
+  }
+#else
+  //in this case we just set the proc ranks for owned and shared FE triangles to the be the same as MPAS owned and shared vertices
+  getProcIds(trianglesProcIds, recvVerticesList_F);
+  for (int i(0); i < nVerticesSolve_F; i++) {
+    if (verticesMask_F[i] & dynamic_ice_bit_value) {
+      fVertexToTriangle[i] = triangleToFVertex.size();
+      triangleToFVertex.push_back(i);
+    }
+  }
+#endif
 
   nTriangles = triangleToFVertex.size();
 
+  //Initialize the ice sheet problem with the number of FE triangles on this prov
   initialize_iceProblem(nTriangles);
 
-  //Compute the global number of triangles, and the localOffset on the local processor, such that a globalID = localOffset + index
-  int localOffset(0);
-  nGlobalTriangles = 0;
-  computeLocalOffset(nTriangles, localOffset, nGlobalTriangles);
-
-  //Communicate the globalIDs, computed locally, to the other processors.
-  indexToTriangleID.resize(nTriangles);
-
-  //To make local, not used
+  //Create a list of global IDs for FE triangles, using MPAS vertices IDs
   fVertexToTriangleID.assign(nVertices_F, NotAnId);
-  //   std::vector<int> fVertexToTriangleID(nVertices_F, NotAnId);
-  for (int index(0); index < nTriangles; index++)
-    fVertexToTriangleID[triangleToFVertex[index]] = index + localOffset;
-
+  for (int index(0); index < nTriangles; index++) {
+    fVertexToTriangleID[triangleToFVertex[index]] = indexToVertexID_F[triangleToFVertex[index]];
+  }
+
+#ifdef changeTrianglesOwnership
+  // because we change the ownership of some triangles, we need to first communicate back to the processors that used to own those triangles
+  // the data of the newly owned triangles. We do this by defining "reversed" send and receive lists, communicate back using those lists, and
+  // then communicate "forward" using the usual send and receive lists.
+  // We could join these two step in one communication, but for the moment we do that separately
+  createReverseExchangeLists(sendVerticesListReversed, recvVerticesListReversed, trianglesProcIds, indexToVertexID_F, recvVerticesList_F);
+  allToAll(fVertexToTriangleID, &sendVerticesListReversed, &recvVerticesListReversed);
+  allToAll(fVertexToTriangle, &sendVerticesListReversed, &recvVerticesListReversed);
+  allToAll(trianglesProcIds, sendVerticesList_F, recvVerticesList_F);
+#endif
   allToAll(fVertexToTriangleID, sendVerticesList_F, recvVerticesList_F);
+  allToAll(fVertexToTriangle, sendVerticesList_F, recvVerticesList_F);
 
-  for (int index(0); index < nTriangles; index++)
+  //we define the vector of global triangles Ids and compute the stride between the largest and the smallest Id globally
+  //This will be needed by the velocity solver to create the 3D FE mesh.
+  indexToTriangleID.resize(nTriangles);
+  int maxTriangleID=std::numeric_limits<int>::min(), minTriangleID=std::numeric_limits<int>::max(), maxGlobalTriangleID, minGlobalTriangleID;
+  for (int index(0); index < nTriangles; index++) {
     indexToTriangleID[index] = fVertexToTriangleID[triangleToFVertex[index]];
+    maxTriangleID = (indexToTriangleID[index] > maxTriangleID) ? indexToTriangleID[index] : maxTriangleID;
+    minTriangleID = (indexToTriangleID[index] < minTriangleID) ? indexToTriangleID[index] : minTriangleID;
+  }
 
-  //Compute triangle edges
-  std::vector<int> fEdgeToEdge(nEdges_F), edgesToSend, trianglesProcIds(
-      nVertices_F);
-  getProcIds(trianglesProcIds, recvVerticesList_F);
+  MPI_Allreduce(&maxTriangleID, &maxGlobalTriangleID, 1, MPI_INT, MPI_MAX, comm);
+  MPI_Allreduce(&minTriangleID, &minGlobalTriangleID, 1, MPI_INT, MPI_MIN, comm);
+  globalTriangleStride = maxGlobalTriangleID - minGlobalTriangleID +1;
 
-  int interfaceSize(0);
+  // Second, we compute the FE edges belonging to the FE triangles owned by this processor.
+  // We first compute boundary edges, and then all the other edges.
+  std::vector<int> fEdgeToEdge(nEdges_F);
 
   std::vector<int> fEdgeToEdgeID(nEdges_F, NotAnId);
-  edgesToReceive.clear();
   edgeToFEdge.clear();
   isBoundaryEdge.clear();
   trianglesOnEdge.clear();
 
-  edgesToReceive.reserve(nEdges_F - nEdgesSolve_F);
   edgeToFEdge.reserve(nEdges_F);
   trianglesOnEdge.reserve(nEdges_F * 2);
-  edgesToSend.reserve(nEdgesSolve_F);
   isBoundaryEdge.reserve(nEdges_F);
 
-  //first, we compute boundary edges (boundary edges must be the first edges)
+  //we compute boundary edges (boundary edges must be the first edges)
   for (int i = 0; i < nEdges_F; i++) {
     ID fVertex1(verticesOnEdge_F[2 * i] - 1), fVertex2(
         verticesOnEdge_F[2 * i + 1] - 1);
+
+    // skip the (shared) edge when the associated MPAS vertices are not  valid
+    if((fVertex1>=nVertices_F) || (fVertex2>=nVertices_F))
+      continue;
+
     ID triaId_1 = fVertexToTriangleID[fVertex1];
     ID triaId_2 = fVertexToTriangleID[fVertex2];
     bool isboundary = (triaId_1 == NotAnId) || (triaId_2 == NotAnId);
 
     ID iTria1 = fVertexToTriangle[fVertex1];
     ID iTria2 = fVertexToTriangle[fVertex2];
-    if (iTria1 == NotAnId)
+    if (trianglesProcIds[fVertex1] != me) {
+      std::swap(fVertex1, fVertex2);
       std::swap(iTria1, iTria2);
-    bool belongsToLocalTriangle = (iTria1 != NotAnId) || (iTria2 != NotAnId);
+    }
+    bool belongsToLocalTriangle = (trianglesProcIds[fVertex1] == me);
 
     if (belongsToLocalTriangle) {
       if (isboundary) {
@@ -639,36 +512,36 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel
         trianglesOnEdge.push_back(iTria1);
         trianglesOnEdge.push_back(iTria2);
         isBoundaryEdge.push_back(true);
-      } else
-        interfaceSize += (iTria2 == NotAnId);
+      }
     }
   }
 
   numBoundaryEdges = edgeToFEdge.size();
 
-  //procOnInterfaceEdge contains the pairs <id of interface edge, rank of processor that owns the non local element adjacent to the edge>.
-  std::vector < std::pair<int, int> > procOnInterfaceEdge;
-  procOnInterfaceEdge.reserve(interfaceSize);
-
   //then, we compute the other edges
   for (int i = 0; i < nEdges_F; i++) {
 
     ID fVertex1(verticesOnEdge_F[2 * i] - 1), fVertex2(
         verticesOnEdge_F[2 * i + 1] - 1);
+
+    // skip the (shared) edge when the associated MPAS vertices are not  valid
+    if((fVertex1>=nVertices_F) || (fVertex2>=nVertices_F))
+      continue;
+
     ID iTria1 = fVertexToTriangle[fVertex1];
     ID iTria2 = fVertexToTriangle[fVertex2];
 
     ID triaId_1 = fVertexToTriangleID[fVertex1]; //global Triangle
     ID triaId_2 = fVertexToTriangleID[fVertex2]; //global Triangle
 
-    if (iTria1 == NotAnId) {
-      std::swap(iTria1, iTria2);
+    if (trianglesProcIds[fVertex1] != me) {
       std::swap(fVertex1, fVertex2);
+      std::swap(iTria1, iTria2);
     }
 
     bool belongsToAnyTriangle = (triaId_1 != NotAnId) || (triaId_2 != NotAnId);
     bool isboundary = (triaId_1 == NotAnId) || (triaId_2 == NotAnId);
-    bool belongsToLocalTriangle = (iTria1 != NotAnId);
+    bool belongsToLocalTriangle = (trianglesProcIds[fVertex1] == me);
     bool isMine = i < nEdgesSolve_F;
 
     if (belongsToLocalTriangle && !isboundary) {
@@ -677,42 +550,33 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel
       trianglesOnEdge.push_back(iTria1);
       trianglesOnEdge.push_back(iTria2);
       isBoundaryEdge.push_back(false);
-      if (iTria2 == NotAnId)
-        procOnInterfaceEdge.push_back(
-            std::make_pair(fEdgeToEdge[i], trianglesProcIds[fVertex2]));
-    }
-
-    if (belongsToAnyTriangle && isMine) {
-      edgesToSend.push_back(i);
-      if (!belongsToLocalTriangle)
-        edgesToReceive.push_back(i);
     }
-
   }
 
-  //Compute the global number of edges, and the localOffset on the local processor, such that a globalID = localOffset + index
-  computeLocalOffset(edgesToSend.size(), localOffset, nGlobalEdges);
-
-  //Communicate the globalIDs, computed locally, to the other processors.
-  for (ID index = 0; index < edgesToSend.size(); index++)
-    fEdgeToEdgeID[edgesToSend[index]] = index + localOffset;
-
-  allToAll(fEdgeToEdgeID, sendEdgesList_F, recvEdgesList_F);
+  for (int fEdge = 0; fEdge < nEdges_F; fEdge++)
+    fEdgeToEdgeID[fEdge] = indexToEdgeID_F[fEdge];
 
   nEdges = edgeToFEdge.size();
   indexToEdgeID.resize(nEdges);
   floatingEdgesIds.clear();
   floatingEdgesIds.reserve(nEdges);
+  int maxEdgeID=std::numeric_limits<int>::min(), minEdgeID=std::numeric_limits<int>::max(), maxGlobalEdgeID, minGlobalEdgeID;
   for (int index = 0; index < nEdges; index++) {
     int fEdge = edgeToFEdge[index];
     indexToEdgeID[index] = fEdgeToEdgeID[fEdge];
+    maxEdgeID = (indexToEdgeID[index] > maxEdgeID) ? indexToEdgeID[index] : maxEdgeID;
+    minEdgeID = (indexToEdgeID[index] < minEdgeID) ? indexToEdgeID[index] : minEdgeID;
     if((floatingEdgesMask_F[fEdge]!=0)&&(index<numBoundaryEdges))
       floatingEdgesIds.push_back(indexToEdgeID[index]);
   }
 
-  //Compute vertices:
-  std::vector<int> fCellsToSend;
-  fCellsToSend.reserve(nCellsSolve_F);
+  MPI_Allreduce(&maxEdgeID, &maxGlobalEdgeID, 1, MPI_INT, MPI_MAX, comm);
+  MPI_Allreduce(&minEdgeID, &minGlobalEdgeID, 1, MPI_INT, MPI_MIN, comm);
+  globalEdgeStride = maxGlobalEdgeID - minGlobalEdgeID + 1;
+
+  // Third, we compute the FE vertices belonging to the FE triangles owned by this processor.
+  // We need to make sure that an FE vertex is owned by a proc that owns a FE triangle that contain that vertex
+  // Otherwise we might end up with weird situation where a vertex could belong to a process with no associated triangle.
 
   vertexToFCell.clear();
   vertexToFCell.reserve(nCells_F);
@@ -720,70 +584,97 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel
   fCellToVertex.assign(nCells_F, NotAnId);
   std::vector<int> fCellToVertexID(nCells_F, NotAnId);
 
-  fCellsToReceive.clear();
+  vertexProcIDs.clear();
+
+  std::vector<int> verticesProcIds(nCells_F, NotAnId);
 
-  // if(! isDomainEmpty)
-  // {
-  fCellsToReceive.reserve(nCells_F - nCellsSolve_F);
+  vertexProcIDs.reserve(nCells_F);
   for (int i = 0; i < nCells_F; i++) {
     bool isMine = i < nCellsSolve_F;
     bool belongsToLocalTriangle = false;
     bool belongsToAnyTriangle = false;
     int nEdg = nEdgesOnCells_F[i];
+    bool invalidVertex = false;
+    int minTriangleProcId = std::numeric_limits<int>::max();
+    bool nodeOwnedByATriaProc = false;
     for (int j = 0; j < nEdg; j++) {
       ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F * i + j] - 1);
-      ID iTria = fVertexToTriangle[fVertex];
+      if(fVertex >= nVertices_F) {
+        invalidVertex = true;
+        break;
+      }
+
       ID triaId = fVertexToTriangleID[fVertex];
-      belongsToLocalTriangle = belongsToLocalTriangle || (iTria != NotAnId);
+      belongsToLocalTriangle = belongsToLocalTriangle || (trianglesProcIds[fVertex] == me);
       belongsToAnyTriangle = belongsToAnyTriangle || (triaId != NotAnId);
-    }
 
-    if (belongsToAnyTriangle && isMine) {
-      fCellsToSend.push_back(i);
-      if (!belongsToLocalTriangle)
-        fCellsToReceive.push_back(i);
+      if(triaId != NotAnId) {
+        nodeOwnedByATriaProc = nodeOwnedByATriaProc || (trianglesProcIds[fVertex] == fCellsProcIds[i]);
+        minTriangleProcId = (trianglesProcIds[fVertex] < minTriangleProcId) ? trianglesProcIds[fVertex] : minTriangleProcId;
+      }
     }
 
+    if(invalidVertex) continue;
+
+    if(belongsToAnyTriangle)
+      verticesProcIds[i] = nodeOwnedByATriaProc ? fCellsProcIds[i] : minTriangleProcId;
+
     if (belongsToLocalTriangle) {
       fCellToVertex[i] = vertexToFCell.size();
       vertexToFCell.push_back(i);
+      vertexProcIDs.push_back(reduced_ranks[verticesProcIds[i]]);
     }
   }
-  //        }
-
-  //Compute the global number of vertices, and the localOffset on the local processor, such that a globalID = localOffset + index
-  computeLocalOffset(fCellsToSend.size(), localOffset, nGlobalVertices);
-
-  //Communicate the globalIDs, computed locally, to the other processors.
-  for (int index = 0; index < int(fCellsToSend.size()); index++)
-    fCellToVertexID[fCellsToSend[index]] = index + localOffset;
+  nVertices = vertexToFCell.size();
 
-  allToAll(fCellToVertexID, sendCellsList_F, recvCellsList_F);
+  for (int fcell = 0; fcell < nCells_F; fcell++)
+    fCellToVertexID[fcell] = indexToCellID_F[fcell];
 
-  nVertices = vertexToFCell.size();
-  int vertexColumnShift = (Ordering == 1) ? 1 : nGlobalVertices;
-  int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1;
+  int maxVertexID=std::numeric_limits<int>::min(), minVertexID=std::numeric_limits<int>::max(), maxGlobalVertexID, minGlobalVertexID;
+  indexToVertexID.resize(nVertices);
+  for (int index = 0; index < nVertices; index++) {
+    int fCell = vertexToFCell[index];
+    indexToVertexID[index] = fCellToVertexID[fCell];
+    maxVertexID = (indexToVertexID[index] > maxVertexID) ? indexToVertexID[index] : maxVertexID;
+    minVertexID = (indexToVertexID[index] < minVertexID) ? indexToVertexID[index] : minVertexID;
+  }
 
+  MPI_Allreduce(&maxVertexID, &maxGlobalVertexID, 1, MPI_INT, MPI_MAX, comm);
+  MPI_Allreduce(&minVertexID, &minGlobalVertexID, 1, MPI_INT, MPI_MIN, comm);
+  globalVertexStride = maxGlobalVertexID - minGlobalVertexID + 1;
 
-  std::cout << "\n nvertices: " << nVertices << " " << nGlobalVertices << "\n"
-      << std::endl;
 
-  indexToVertexID.resize(nVertices);
+  int vertexColumnShift = (Ordering == 1) ? 1 : globalVertexStride;
+  int vertexLayerShift = (Ordering == 0) ? 1 : nLayers + 1;
   dirichletNodesIDs.clear();
   dirichletNodesIDs.reserve(nVertices); //need to improve storage efficiency
+  isVertexBoundary.assign(nVertices, false);
   for (int index = 0; index < nVertices; index++) {
     int fCell = vertexToFCell[index];
-    indexToVertexID[index] = fCellToVertexID[fCell];
     for(int il=0; il< nLayers+1; ++il)
     {
       int imask_F = il+(nLayers+1)*fCell;
       if(dirichletCellsMask_F[imask_F]!=0)
         dirichletNodesIDs.push_back((nLayers-il)*vertexColumnShift+indexToVertexID[index]*vertexLayerShift);
     }
+
+    int nEdg = nEdgesOnCells_F[fCell];
+    int j = 0;
+    bool isBoundary;
+    do {
+      int fVertex = verticesOnCell_F[maxNEdgesOnCell_F * fCell + j++] - 1;
+      isBoundary = !(verticesMask_F[fVertex] & dynamic_ice_bit_value);
+    } while ((j < nEdg) && (!isBoundary));
+    isVertexBoundary[index] = isBoundary;
   }
 
-  createReverseCellsExchangeLists(sendCellsListReversed, recvCellsListReversed,
-      fVertexToTriangleID, fCellToVertexID);
+  // because we change the ownership of some vertices, we need to first communicate back to the processors that used to own those vertices
+  // the data of the newly owned vertices. We do this by defining "reversed" send and receive lists, communicate back using that list, and
+  // then communicate forward with the usual send and receive lists.
+  // We could join these two step in one communication, but for the moment we do that separately.
+  // We need to communicate info about the vertices when we get the ice velocity on vertices form the velocity solver/
+
+  createReverseExchangeLists(sendCellsListReversed, recvCellsListReversed, verticesProcIds, indexToCellID_F, recvCellsList_F);
 
   //construct the local vector vertices on triangles making sure the area is positive
   verticesOnTria.resize(nTriangles * 3);
@@ -802,132 +693,24 @@ void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _cel
       std::swap(verticesOnTria[3 * index + 1], verticesOnTria[3 * index + 2]);
   }
 
-  //construct the local vector vertices on edges
-  trianglesPositionsOnEdge.resize(2 * nEdges);
-  isVertexBoundary.assign(nVertices, false);
-
   verticesOnEdge.resize(2 * nEdges);
-
-  //contains the local id of a triangle and the global id of the edges of the triangle.
-  //dataForGhostTria[4*i] contains the triangle id
-  //dataForGhostTria[4*i+1+k] contains the global id of the edge (at position k = 0,1,2) of the triangle.
-  //Possible Optimization: for our purposes it would be enough to store two of the three edges of a triangle.
-  std::vector<int> dataForGhostTria(nVertices_F * 4, NotAnId);
-
-  //*
-  for (int iV = 0; iV < nVertices; iV++) {
-    int fCell = vertexToFCell[iV];
-    int nEdg = nEdgesOnCells_F[fCell];
-    int j = 0;
-    bool isBoundary;
-    do {
-      int fVertex = verticesOnCell_F[maxNEdgesOnCell_F * fCell + j++] - 1;
-      isBoundary = !(verticesMask_F[fVertex] & dynamic_ice_bit_value);
-    } while ((j < nEdg) && (!isBoundary));
-    isVertexBoundary[iV] = isBoundary;
-  }
-  /*/
-   for(int index=0; index<numBoundaryEdges; index++)
-   {
-   int fEdge = edgeToFEdge[index];
-   int v1 = fCellToVertex[cellsOnEdge_F[2*fEdge]-1];
-   int v2 = fCellToVertex[cellsOnEdge_F[2*fEdge+1]-1];
-   isVertexBoundary[ v1 ] = isVertexBoundary[ v2 ] = true;
-   }
-   //*/
-  //computing the position and local ids of the triangles adjacent to edges.
-  //in the case an adjacent triangle is not local, the position and id will be NotAnId.
-  //We will get the needed information about the non local edge leter on,
-  //for this purpose we fill the vector dataForGhostTria.
   for (int index = 0; index < nEdges; index++) {
-    int p1, p2, v1, v2, v3, t, position;
     int fEdge = edgeToFEdge[index];
     int fCell1 = cellsOnEdge_F[2 * fEdge] - 1;
     int fCell2 = cellsOnEdge_F[2 * fEdge + 1] - 1;
-    verticesOnEdge[2 * index] = p1 = fCellToVertex[fCell1];
-    verticesOnEdge[2 * index + 1] = p2 = fCellToVertex[fCell2];
-
-    for (int k = 0; k < 2 && (t = trianglesOnEdge[2 * index + k]) != NotAnId;
-        k++) {
-      v1 = verticesOnTria[3 * t];
-      v2 = verticesOnTria[3 * t + 1];
-      v3 = verticesOnTria[3 * t + 2];
-      position = (((p1 == v2) && (p2 == v3)) || ((p1 == v3) && (p2 == v2)))
-          + 2 * (((p1 == v1) && (p2 == v3)) || ((p1 == v3) && (p2 == v1)));
-      trianglesPositionsOnEdge[2 * index + k] = position;
-      int dataIndex = 4 * triangleToFVertex[t];
-
-      dataForGhostTria[dataIndex] = t;
-      dataForGhostTria[dataIndex + position + 1] = indexToEdgeID[index];
-    }
+    verticesOnEdge[2 * index] = fCellToVertex[fCell1];
+    verticesOnEdge[2 * index + 1] = fCellToVertex[fCell2];
   }
 
-  createReverseEdgesExchangeLists(sendEdgesListReversed, recvEdgesListReversed,
-      fVertexToTriangleID, fEdgeToEdgeID);
-
-  //send the information about ghost elements.
-  allToAll(dataForGhostTria, sendVerticesList_F, recvVerticesList_F, 4);
-
-  //retrieving the position, local id and owner processor of non local triangles adjacent to interface edges.
-  for (int index = 0; index < (int) procOnInterfaceEdge.size(); index++) {
-    int iEdge = procOnInterfaceEdge[index].first;
-    int fEdge = edgeToFEdge[iEdge];
-    ID fVertex1(verticesOnEdge_F[2 * fEdge] - 1), fVertex2(
-        verticesOnEdge_F[2 * fEdge + 1] - 1);
-    if ((ID) fVertexToTriangle[fVertex1] == NotAnId)
-      fVertex2 = fVertex1;
-    int dataIndex = 4 * fVertex2;
-    int edgeID = indexToEdgeID[iEdge];
-    int position = (dataForGhostTria[dataIndex + 2] == edgeID)
-        + 2 * (dataForGhostTria[dataIndex + 3] == edgeID);
-    trianglesOnEdge[2 * iEdge + 1] = dataForGhostTria[dataIndex];
-    trianglesPositionsOnEdge[2 * iEdge + 1] = position;
-  }
 
+  //call the velocity solver only on procs owning some FE triangles
   if (isDomainEmpty)
     return;
 
-#ifdef LIFEV
-    std::vector<double> verticesCoords(3 * nVertices);
-
-    for (int index = 0; index < nVertices; index++) {
-      int iCell = vertexToFCell[index];
-      verticesCoords[index * 3] = xCell_F[iCell] / unit_length;
-      verticesCoords[index * 3 + 1] = yCell_F[iCell] / unit_length;
-      verticesCoords[index * 3 + 2] = zCell_F[iCell] / unit_length;
-    }
-
-    velocity_solver_compute_2d_grid__(nGlobalTriangles,
-        nGlobalVertices, nGlobalEdges, indexToVertexID, verticesCoords,
-        isVertexBoundary, verticesOnTria,isBoundaryEdge, trianglesOnEdge,
-        trianglesPositionsOnEdge, verticesOnEdge, indexToEdgeID,
-        indexToTriangleID, procOnInterfaceEdge );
-#else
-    velocity_solver_compute_2d_grid__(reducedComm);
-#endif
-
-  /*
-
-   //initialize the mesh
-   iceProblemPtr->mesh2DPtr.reset (new RegionMesh<LinearTriangle>() );
-
-   //construct the mesh nodes
-   constructNodes ( * (iceProblemPtr->mesh2DPtr), indexToVertexID, verticesCoords, isVertexBoundary, nGlobalVertices, 3);
-
-   //construct the mesh elements
-   constructElements ( * (iceProblemPtr->mesh2DPtr), indexToTriangleID, verticesOnTria, nGlobalTriangles);
-
-   //construct the mesh facets
-   constructFacets ( * (iceProblemPtr->mesh2DPtr), isBoundaryEdge, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge, indexToEdgeID, procOnInterfaceEdge, nGlobalEdges, 3);
-
-   Switch sw;
-   std::vector<bool> elSign;
-   checkVolumes ( * (iceProblemPtr->mesh2DPtr), elSign, sw );
-   */
+  velocity_solver_compute_2d_grid__(reducedComm);
 }
 
-void velocity_solver_extrude_3d_grid(double const* levelsRatio_F,
-    double const* lowerSurface_F, double const* thickness_F) {
+void velocity_solver_extrude_3d_grid(double const* levelsRatio_F) {
 
   if (isDomainEmpty)
     return;
@@ -936,7 +719,6 @@ void velocity_solver_extrude_3d_grid(double const* levelsRatio_F,
   // !!Indexing of layers is reversed
   for (int i = 0; i < nLayers; i++)
     layersRatio[i] = levelsRatio_F[nLayers - 1 - i];
-  //std::copy(levelsRatio_F, levelsRatio_F+nLayers, layersRatio.begin());
 
   levelsNormalizedThickness.resize(nLayers + 1);
 
@@ -945,10 +727,6 @@ void velocity_solver_extrude_3d_grid(double const* levelsRatio_F,
     levelsNormalizedThickness[i + 1] = levelsNormalizedThickness[i]
         + layersRatio[i];
 
-  std::vector<int> mpasIndexToVertexID(nVertices);
-  for (int i = 0; i < nVertices; i++)
-    mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]];
-
   //construct the local vector of coordinates
   std::vector<double> verticesCoords(3 * nVertices);
 
@@ -959,10 +737,14 @@ void velocity_solver_extrude_3d_grid(double const* levelsRatio_F,
     verticesCoords[index * 3 + 2] = zCell_F[iCell] / unit_length;
   }
 
-  velocity_solver_extrude_3d_grid__(nLayers, nGlobalTriangles, nGlobalVertices,
-      nGlobalEdges, Ordering, reducedComm, indexToVertexID, mpasIndexToVertexID,
-      verticesCoords, isVertexBoundary, verticesOnTria, isBoundaryEdge,
-      trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge, indexToEdgeID,
+  std::vector<std::vector<int>> procsSharingVertices(nVertices);
+  for(int i=0; i< nVertices; i++)
+    procsSharingVertex(i, procsSharingVertices[i]);
+
+  velocity_solver_extrude_3d_grid__(nLayers, globalTriangleStride, globalVertexStride,
+      globalEdgeStride, Ordering, reducedComm, indexToVertexID, vertexProcIDs,
+      verticesCoords, verticesOnTria, procsSharingVertices, isBoundaryEdge,
+      trianglesOnEdge, verticesOnEdge, indexToEdgeID,
       indexToTriangleID, dirichletNodesIDs, floatingEdgesIds);
   }
 
@@ -1004,11 +786,7 @@ void interface_init_log(){
   }
   strcat(oformat, ".out");
 
-  if (me == 0) {
-    sprintf(albany_log_filename, oformat, me);
-  } else {
-    strcpy(albany_log_filename, "/dev/null");
-  }
+  sprintf(albany_log_filename, oformat, me);
 
   Interface_stdout = open(albany_log_filename, O_CREAT|O_WRONLY|O_TRUNC,0644);
   if(Interface_stdout >=  0) {
@@ -1070,14 +848,9 @@ void get_prism_velocity_on_FEdges(double * uNormal,
 
   UInt nPoints3D = nCells_F * (nLayers + 1);
 
-  // Loop over all edges of the triangulation - MPAS will decide which edges it should use.
-  for (int i = 0; i < nEdges; i++) {
-
-    //identifying vertices on the edge
-    ID lId0 = verticesOnEdge[2 * i];
-    ID lId1 = verticesOnEdge[2 * i + 1];
-    int iCell0 = vertexToFCell[lId0];
-    int iCell1 = vertexToFCell[lId1];
+  for (int iEdge = 0; iEdge < nEdgesSolve_F; iEdge++) {
+    int iCell0 = cellsOnEdge_F[2 * iEdge] - 1;
+    int iCell1 = cellsOnEdge_F[2 * iEdge + 1] - 1;
 
     //computing normal to the cell edge (dual of triangular edge)
     double nx = xCell_F[iCell1] - xCell_F[iCell0];
@@ -1087,9 +860,14 @@ void get_prism_velocity_on_FEdges(double * uNormal,
     ny /= n;
 
     //identifying triangles that shares the edge
-    ID iEdge = edgeToFEdge[i];
     ID fVertex0 = verticesOnEdge_F[2 * iEdge] - 1;
     ID fVertex1 = verticesOnEdge_F[2 * iEdge + 1] - 1;
+
+
+    int iTria0 = fVertexToTriangleID[fVertex0];
+    int iTria1 = fVertexToTriangleID[fVertex1];
+    if((iTria0 == NotAnId) && (iTria1 == NotAnId)) continue;
+
     double t0[2*3], t1[2*3]; //t0[0] contains the x-coords of vertices of triangle 0 and t0[1] its y-coords.
     for (int j = 0; j < 3; j++) {
       int iCell = cellsOnVertex_F[3 * fVertex0 + j] - 1;
@@ -1120,11 +898,11 @@ void get_prism_velocity_on_FEdges(double * uNormal,
       for (int j = 0; j < 3; j++)
         iCells[j] = cellsOnVertex_F[3 * fVertex1 + j] - 1;
       }
-   else if(i<numBoundaryEdges) {
+   else if((iTria0 == NotAnId) || (iTria1 == NotAnId)) { //is on Boundary
       //edge is on boundary and wasn't found by previous two cases
       //For boundary edges one of the two triangles sharing the edge won't be part of the velocity solver's mesh and the dynamic_ice_bit_value will be 0.
 
-      //Compute iCells containinig the vertices of the triangle that is part of the mesh and bcoords the corresponding barycentric coordinates
+      //Compute iCells containing the vertices of the triangle that is part of the mesh and bcoords the corresponding barycentric coordinates
       if(verticesMask_F[fVertex0] & dynamic_ice_bit_value) { belongToTria(e_mid, t0, bcoords);
         for (int j = 0; j < 3; j++)
           iCells[j] = cellsOnVertex_F[3 * fVertex0 + j] - 1;
@@ -1144,7 +922,7 @@ void get_prism_velocity_on_FEdges(double * uNormal,
       //Scale the coordinates so that they sum to 1.
       for(int j=0; j<3; j++)
         bcoords[j] /= sum;
-      }    
+      }
    else { //error, edge midpont does not belong to either triangle
       std::cout << "Error, edge midpont does not belong to either triangle" << std::endl;
 
@@ -1157,7 +935,6 @@ void get_prism_velocity_on_FEdges(double * uNormal,
       exit(1);
    }
 
-
     for (int il = 0; il < nLayers+1; il++) { //loop over layers
       int ilReversed = nLayers - il;
       int index = iEdge * (nLayers+1) + ilReversed;
@@ -1210,113 +987,24 @@ void mapVerticesToCells(const std::vector<double>& velocityOnVertices,
   }
 }
 
-void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F,
-    exchangeList_Type& receiveListReverse_F,
-    const std::vector<int>& fVertexToTriangleID,
-    const std::vector<int>& fCellToVertexID) {
-  sendListReverse_F.clear();
-  receiveListReverse_F.clear();
-  //std::map<int, std::vector<int> > sendMap;
-  std::map<int, std::map<int, int> > sendMap, receiveMap;
-  std::vector<int> cellsProcId(nCells_F), trianglesProcIds(nVertices_F);
-  getProcIds(cellsProcId, recvCellsList_F);
-  getProcIds(trianglesProcIds, recvVerticesList_F);
-
-  //std::cout << "SendList " ;
-  for (int i = 0; i < nVertices; i++) {
-    int iCell = vertexToFCell[i];
-    if (iCell < nCellsSolve_F)
-      continue;
-    bool belongToTriaOnSameProc = false;
-    int j(0);
-    int nEdg = nEdgesOnCells_F[iCell];
-    do {
-      ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F * iCell + j] - 1);
-      ID triaId = fVertexToTriangleID[fVertex];
-      belongToTriaOnSameProc = (triaId != NotAnId)
-          && (trianglesProcIds[fVertex] == cellsProcId[iCell]);
-    } while ((belongToTriaOnSameProc == false) && (++j < nEdg));
-    if (!belongToTriaOnSameProc) {
-      sendMap[cellsProcId[iCell]].insert(
-          std::make_pair(fCellToVertexID[iCell], iCell));
-      // std::cout<< "(" << cellsProcId[iCell] << "," << iCell << ") ";
-    }
-
-  }
-  //std::cout <<std::endl;
-
-  for (std::map<int, std::map<int, int> >::const_iterator it = sendMap.begin();
-      it != sendMap.end(); it++) {
-    std::vector<int> sendVec(it->second.size());
-    int i = 0;
-    for (std::map<int, int>::const_iterator iter = it->second.begin();
-        iter != it->second.end(); iter++)
-      sendVec[i++] = iter->second;
-    sendListReverse_F.push_back(
-        exchange(it->first, &sendVec[0], &sendVec[0] + sendVec.size()));
-  }
-
-  //std::cout << "ReceiveList " ;
-  for (UInt i = 0; i < fCellsToReceive.size(); i++) {
-    int iCell = fCellsToReceive[i];
-    int nEdg = nEdgesOnCells_F[iCell];
-    for (int j = 0; j < nEdg; j++) {
-      ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F * iCell + j] - 1);
-      ID triaId = fVertexToTriangleID[fVertex];
-      if (triaId != NotAnId) {
-        receiveMap[trianglesProcIds[fVertex]].insert(
-            std::make_pair(fCellToVertexID[iCell], iCell));
-        // std::cout<< "(" << trianglesProcIds[fVertex] << "," << iCell << ") ";
-      }
-    }
-  }
-  //std::cout <<std::endl;
 
-  for (std::map<int, std::map<int, int> >::const_iterator it =
-      receiveMap.begin(); it != receiveMap.end(); it++) {
-    std::vector<int> receiveVec(it->second.size());
-    int i = 0;
-    for (std::map<int, int>::const_iterator iter = it->second.begin();
-        iter != it->second.end(); iter++)
-      receiveVec[i++] = iter->second;
-    receiveListReverse_F.push_back(
-        exchange(it->first, &receiveVec[0],
-            &receiveVec[0] + receiveVec.size()));
-  }
-}
-
-void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F,
+void createReverseExchangeLists(exchangeList_Type& sendListReverse_F,
     exchangeList_Type& receiveListReverse_F,
-    const std::vector<int>& fVertexToTriangleID,
-    const std::vector<int>& fEdgeToEdgeID) {
+    const std::vector<int>& newProcIds, const int* indexToID_F, exchangeList_Type const * recvList_F) {
   sendListReverse_F.clear();
   receiveListReverse_F.clear();
-  //std::map<int, std::vector<int> > sendMap;
   std::map<int, std::map<int, int> > sendMap, receiveMap;
-  std::vector<int> edgesProcId(nEdges_F), trianglesProcIds(nVertices_F);
-  getProcIds(edgesProcId, recvEdgesList_F);
-  getProcIds(trianglesProcIds, recvVerticesList_F);
-
-  //std::cout << "EdgesSendList " ;
-  for (int i = 0; i < nEdges; i++) {
-    int iEdge = edgeToFEdge[i];
-    if (iEdge < nEdgesSolve_F)
-      continue;
-    bool belongToTriaOnSameProc = false;
-    int j(0);
-    do {
-      ID fVertex(verticesOnEdge_F[2 * iEdge + j] - 1);
-      ID triaId = fVertexToTriangleID[fVertex];
-      belongToTriaOnSameProc = (triaId != NotAnId)
-          && (trianglesProcIds[fVertex] == edgesProcId[iEdge]);
-    } while ((belongToTriaOnSameProc == false) && (++j < 2));
-    if (!belongToTriaOnSameProc) {
-      sendMap[edgesProcId[iEdge]].insert(
-          std::make_pair(fEdgeToEdgeID[iEdge], iEdge));
-      //std::cout<< "(" << edgesProcId[iEdge] << "," << iEdge << ") ";
+  int nFEntities = newProcIds.size();
+  std::vector<int> procIds(nFEntities);
+  getProcIds(procIds, recvList_F);
+  int me;
+  MPI_Comm_rank(comm, &me);
+  for (int fEntity = 0; fEntity < nFEntities; fEntity++) {
+    if ((procIds[fEntity] != me) && (newProcIds[fEntity] == me)) {
+      sendMap[procIds[fEntity]].insert(
+        std::make_pair(indexToID_F[fEntity], fEntity));
     }
   }
-  //std::cout <<std::endl;
 
   for (std::map<int, std::map<int, int> >::const_iterator it = sendMap.begin();
       it != sendMap.end(); it++) {
@@ -1329,20 +1017,12 @@ void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F,
         exchange(it->first, &sendVec[0], &sendVec[0] + sendVec.size()));
   }
 
-  //std::cout << "EdgesReceiveList " ;
-  for (UInt i = 0; i < edgesToReceive.size(); i++) {
-    int iEdge = edgesToReceive[i];
-    for (int j = 0; j < 2; j++) {
-      ID fVertex(verticesOnEdge_F[2 * iEdge + j] - 1);
-      ID triaId = fVertexToTriangleID[fVertex];
-      if (triaId != NotAnId) {
-        receiveMap[trianglesProcIds[fVertex]].insert(
-            std::make_pair(fEdgeToEdgeID[iEdge], iEdge));
-        //  std::cout<< "(" << trianglesProcIds[fVertex] << "," << iEdge << ") ";
-      }
+  for (int fEntity = 0; fEntity < nFEntities; fEntity++) {
+    if((procIds[fEntity] == me) && (newProcIds[fEntity] != NotAnId) && (newProcIds[fEntity] != me)) {
+      receiveMap[newProcIds[fEntity]].insert(
+        std::make_pair(indexToID_F[fEntity], fEntity));
     }
   }
-  // std::cout <<std::endl;
 
   for (std::map<int, std::map<int, int> >::const_iterator it =
       receiveMap.begin(); it != receiveMap.end(); it++) {
@@ -1381,19 +1061,6 @@ void mapCellsToVertices(const std::vector<double>& velocityOnCells,
   }
 }
 
-bool isGhostTriangle(int i, double relTol) {
-  double x[3], y[3], area;
-
-  for (int j = 0; j < 3; j++) {
-    int iCell = cellsOnVertex_F[3 * i + j] - 1;
-    x[j] = xCell_F[iCell];
-    y[j] = yCell_F[iCell];
-  }
-
-  area = std::fabs(signedTriangleArea(x, y));
-  return false; //(std::fabs(areaTriangle_F[i]-area)/areaTriangle_F[i] > relTol);
-}
-
 double signedTriangleArea(const double* x, const double* y) {
   double u[2] = { x[1] - x[0], y[1] - y[0] };
   double v[2] = { x[2] - x[0], y[2] - y[0] };
@@ -1417,27 +1084,8 @@ double signedTriangleAreaOnSphere(const double* x, const double* y,
           > 0) ? area : -area;
 }
 
-//TO BE FIXED, Access To verticesOnCell_F is not correct
-void extendMaskByOneLayer(int const* verticesMask_F,
-    std::vector<int>& extendedFVerticesMask) {
-  extendedFVerticesMask.resize(nVertices_F);
-  extendedFVerticesMask.assign(&verticesMask_F[0],
-      &verticesMask_F[0] + nVertices_F);
-  for (int i = 0; i < nCells_F; i++) {
-    bool belongsToMarkedTriangle = false;
-    int nEdg = nEdgesOnCells_F[i];
-    for (UInt k = 0; k < nEdg && !belongsToMarkedTriangle; k++)
-      belongsToMarkedTriangle = belongsToMarkedTriangle
-          || verticesMask_F[verticesOnCell_F[maxNEdgesOnCell_F * i + k] - 1];
-    if (belongsToMarkedTriangle)
-      for (UInt k = 0; k < nEdg; k++) {
-        ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F * i + k] - 1);
-        extendedFVerticesMask[fVertex] = !isGhostTriangle(fVertex);
-      }
-  }
-}
 
-void import2DFields(std::map<int, int> bdExtensionMap, double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F,
+void importFields(std::map<int, int> bdExtensionMap, double const* bedTopography_F, double const * lowerSurface_F, double const * thickness_F,
     double const * beta_F, double const* stiffnessFactor_F, double const* effecPress_F,
     double const * temperature_F, double const * smb_F, double eps) {
 
@@ -1452,7 +1100,7 @@ void import2DFields(std::map<int, int> bdExtensionMap, double const* bedTopograp
   if (effecPress_F != 0)
     effecPressData.assign(nVertices, 1e10);
   if(temperature_F != 0)
-    temperatureData.assign(nLayers * nTriangles, 1e10);
+    temperatureDataOnPrisms.assign(nLayers * nTriangles, 1e10);
   if (smb_F != 0)
     smbData.assign(nVertices, 1e10);
 
@@ -1473,6 +1121,8 @@ void import2DFields(std::map<int, int> bdExtensionMap, double const* bedTopograp
       effecPressData[index] = effecPress_F[iCell] / unit_length;  
   }
 
+  int lElemColumnShift = (Ordering == 1) ? 1 : nTriangles;
+  int elemLayerShift = (Ordering == 0) ? 1 : nLayers;
   if(temperature_F != 0) {
     for (int index = 0; index < nTriangles; index++) {
       for (int il = 0; il < nLayers; il++) {
@@ -1482,21 +1132,20 @@ void import2DFields(std::map<int, int> bdExtensionMap, double const* bedTopograp
         for (int iVertex = 0; iVertex < 3; iVertex++) {
       int v = verticesOnTria[iVertex + 3 * index];
       int iCell = vertexToFCell[v];
-          //compute temperature by averaging tmeperature values of triangles vertices where ice is present
-      if (cellsMask_F[iCell] & ice_present_bit_value) { 
+      //compute temperature by averaging temperature values of triangles vertices where ice is present
+      if (cellsMask_F[iCell] & ice_present_bit_value) {
         temperature += temperature_F[iCell * nLayers + ilReversed];
         nPoints++;
          }
         }
         if (nPoints == 0)  //if triangle is in an ice-free area, set the temperature to T0
-      temperatureData[index+il*nTriangles] = T0;
+      temperatureDataOnPrisms[index*elemLayerShift + il*lElemColumnShift] = T0;
         else
-      temperatureData[index+il*nTriangles] = temperature / nPoints;
+      temperatureDataOnPrisms[index*elemLayerShift + il*lElemColumnShift] = temperature / nPoints;
       }
     }
   }
 
-
   //extend thickness elevation and basal friction data to the border for floating vertices
   std::set<int>::const_iterator iter;
 
@@ -1519,18 +1168,11 @@ void import2DFields(std::map<int, int> bdExtensionMap, double const* bedTopograp
          bool foundNeighbor = false;
          for (int j = 0; j < nEdg; j++) {
            int fEdge = edgesOnCell_F[maxNEdgesOnCell_F * fCell + j] - 1;
-           //bool keep = (mask[verticesOnEdge_F[2 * fEdge] - 1] & dynamic_ice_bit_value)
-           //    && (mask[verticesOnEdge_F[2 * fEdge + 1] - 1] & dynamic_ice_bit_value);
-           //if (!keep)
-           //  continue;
-
            int c0 = cellsOnEdge_F[2 * fEdge] - 1;
            int c1 = cellsOnEdge_F[2 * fEdge + 1] - 1;
            c = (fCellToVertex[c0] == iV) ? c1 : c0;
-           //if(!(cellsMask_F[c] & ice_present_bit_value)) continue;
            if((cellsMask_F[c] & dynamic_ice_bit_value)) {
            double elev = thickness_F[c] + lowerSurface_F[c]; // - 1e-8*std::sqrt(pow(xCell_F[c0],2)+std::pow(yCell_F[c0],2));
-           //std::cout << "  elev="<<elev<<std::endl;
            if (elev < elevTemp) {
              elevTemp = elev;
              bdExtensionMap[iV] = c;
@@ -1545,7 +1187,7 @@ void import2DFields(std::map<int, int> bdExtensionMap, double const* bedTopograp
             // Check that this margin location is not below sea level!
             if (bedTopographyData[iV] < 0.0) { // This check is probably redundant...
                thicknessData[iV] = eps*3.0; // insert special value here to make identifying these points easier in exo output
-               elevationData[iV] = (rho_ocean / rho_ice - 1.0) * thicknessData[iV];  // floating surface
+               elevationData[iV] = (1.0 - rho_ice / rho_ocean) * thicknessData[iV];  // floating surface
             }
          }
       } else {
@@ -1554,14 +1196,10 @@ void import2DFields(std::map<int, int> bdExtensionMap, double const* bedTopograp
         // Otherwise, it will have a surface elevation below sea level,
         // which is unphysical and can cause large slopes and other issues.
         if (bedTopographyData[iV] < 0.0) {
-          //std::cout<<"non-floating boundary below sea level at iV="<<iV<<std::endl;
           if (std::find(dirichletNodesIDs.begin(), dirichletNodesIDs.end(), indexToVertexID[iV]*vertexLayerShift) != dirichletNodesIDs.end()) { // Don't do this if location is a Dirichlet node!
-            //std::cout<<"  non-floating boundary below sea level node is Dirichlet so skipping. iV="<<iV<<std::endl;
           } else {
-            //for (int i = 0; i < dirichletNodesIDs.size(); i++)
-            //            std::cout << dirichletNodesIDs.at(i) << ' ';  // print entire list of Diri nodes for debugging.
             thicknessData[iV] = eps*2.0; // insert special small value here to make identifying these points easier in exo output
-            elevationData[iV] = (rho_ocean / rho_ice - 1.0) * thicknessData[iV];  // floating surface
+            elevationData[iV] = (1.0 - rho_ice / rho_ocean) * thicknessData[iV];  // floating surface
           }
         } // if below sea level
       }  // floating or not
@@ -1628,7 +1266,7 @@ void import2DFieldsObservations(std::map<int, int> bdExtensionMap,
     indexToCellIDData[index] = indexToCellID_F[iCell];
   }
 
-  //extend to the border for floating vertices (using map created by import2DFields above)
+  //extend to the border for floating vertices (using map created by importFields above)
   for (std::map<int, int>::iterator it = bdExtensionMap.begin();
       it != bdExtensionMap.end(); ++it) {
     int iv = it->first;
@@ -1649,36 +1287,20 @@ void import2DFieldsObservations(std::map<int, int> bdExtensionMap,
   }
 }
 
-
-void importP0Temperature() {
-  int lElemColumnShift = (Ordering == 1) ? 3 : 3 * indexToTriangleID.size();
-  int elemLayerShift = (Ordering == 0) ? 3 : 3 * nLayers;
-  temperatureOnTetra.resize(3 * nLayers * indexToTriangleID.size());
-  for (int index = 0; index < nTriangles; index++) {
-    for (int il = 0; il < nLayers; il++) {
-      for (int k = 0; k < 3; k++)
-        temperatureOnTetra[index * elemLayerShift + il * lElemColumnShift + k] =
-            temperatureData[index+il*nTriangles];
-    }
-  }
-}
-
 void exportDissipationHeat(double * dissipationHeat_F) {
   std::fill(dissipationHeat_F, dissipationHeat_F + nVertices_F * (nLayers), 0.);
-  int lElemColumnShift = (Ordering == 1) ? 3 : 3 * indexToTriangleID.size();
-  int elemLayerShift = (Ordering == 0) ? 3 : 3 * nLayers;
+  int lElemColumnShift = (Ordering == 1) ? 1 : indexToTriangleID.size();
+  int elemLayerShift = (Ordering == 0) ? 1 : nLayers;
   for (int index = 0; index < nTriangles; index++) {
     for (int il = 0; il < nLayers; il++) {
       int ilReversed = nLayers - il - 1;
       int fVertex = triangleToFVertex[index];
-      double dissipationHeat = 0;
-      for (int k = 0; k < 3; k++)
-        dissipationHeat += dissipationHeatOnTetra[index * elemLayerShift + il * lElemColumnShift + k]/3; //TODO: should be weighted average based on tets volumes
-      dissipationHeat_F[fVertex * nLayers + ilReversed] = dissipationHeat;
-
+      dissipationHeat_F[fVertex * nLayers + ilReversed] = dissipationHeatOnPrisms[index * elemLayerShift + il * lElemColumnShift];
     }
   }
-  //allToAll (dissipationHeat_F,  &sendVerticesListReversed, &recvEdgesListReversed, nLayers+1);
+#ifdef  changeTrianglesOwnership
+  allToAll (dissipationHeat_F,  &sendVerticesListReversed, &recvVerticesListReversed, nLayers);
+#endif
   allToAll (dissipationHeat_F,  sendVerticesList_F, recvVerticesList_F, nLayers);
 }
 
@@ -1693,7 +1315,7 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id) {
   int nonEmpty = int(nLocalEntities > 0);
   MPI_Allgather(&nonEmpty, 1, MPI_INT, &haveElements[0], 1, MPI_INT, comm);
   std::vector<int> ranks;
-  reduced_ranks.resize(numProcs,0);
+  reduced_ranks.assign(numProcs,-1);
   for (int i = 0; i < numProcs; i++) {
     if (haveElements[i]) {
       reduced_ranks[i] = ranks.size();
@@ -1706,24 +1328,6 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id) {
   MPI_Comm_create(comm, reduced_group_id, &reduced_comm_id);
 }
 
-void computeLocalOffset(int nLocalEntities, int& localOffset,
-    int& nGlobalEntities) {
-  int numProcs, me;
-  MPI_Comm_size(comm, &numProcs);
-  MPI_Comm_rank(comm, &me);
-  std::vector<int> offsetVec(numProcs);
-
-  MPI_Allgather(&nLocalEntities, 1, MPI_INT, &offsetVec[0], 1, MPI_INT, comm);
-
-  localOffset = 0;
-  for (int i = 0; i < me; i++)
-    localOffset += offsetVec[i];
-
-  nGlobalEntities = localOffset;
-  for (int i = me; i < numProcs; i++)
-    nGlobalEntities += offsetVec[i];
-}
-
 void getProcIds(std::vector<int>& field, int const * recvArray) {
   int me;
   MPI_Comm_rank(comm, &me);
@@ -1913,13 +1517,6 @@ int initialize_iceProblem(int nTriangles) {
 
   isDomainEmpty = !keep_proc;
 
-
-#ifdef LIFEV
-if(!isDomainEmpty) {
-  velocity_solver_initialize_iceProblem__(keep_proc, reducedComm);
-}
-#endif
-
   // initialize ice problem pointer
   if (keep_proc) {
     std::cout << nTriangles
@@ -1954,164 +1551,12 @@ bool belongToTria(double const* x, double const* t, double bcoords[3], double ep
          ( bcoords[2] > -eps );
 }
 
-int prismType(long long int const* prismVertexMpasIds, int& minIndex)
-{
-  int PrismVerticesMap[6][6] = {{0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}, {3, 5, 4, 0, 2, 1}, {4, 3, 5, 1, 0, 2}, {5, 4, 3, 2, 1, 0}};
-  minIndex = std::min_element (prismVertexMpasIds, prismVertexMpasIds + 3) - prismVertexMpasIds;
-
-  int v1 (prismVertexMpasIds[PrismVerticesMap[minIndex][1]]);
-  int v2 (prismVertexMpasIds[PrismVerticesMap[minIndex][2]]);
-
-  return v1  > v2;
-}
-
- void tetrasFromPrismStructured (long long int const* prismVertexMpasIds, long long int const* prismVertexGIds, long long int tetrasIdsOnPrism[][4])
-  {
-      int PrismVerticesMap[6][6] = {{0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}, {3, 5, 4, 0, 2, 1}, {4, 3, 5, 1, 0, 2}, {5, 4, 3, 2, 1, 0}};
-
-      int tetraOfPrism[2][3][4] = {{{0, 1, 2, 5}, {0, 1, 5, 4}, {0, 4, 5, 3}}, {{0, 1, 2, 4}, {0, 4, 2, 5}, {0, 4, 5, 3}}};
-
-      int tetraAdjacentToPrismLateralFace[2][3][2] = {{{1, 2}, {0, 1}, {0, 2}}, {{0, 2}, {0, 1}, {1, 2}}};
-      int tetraFaceIdOnPrismLateralFace[2][3][2] = {{{0, 0}, {1, 1}, {2, 2}}, {{0, 0}, {1, 1}, {2, 2}}};
-      int tetraAdjacentToBottomFace = 0; //does not depend on type;
-      int tetraAdjacentToUpperFace = 2; //does not depend on type;
-      int tetraFaceIdOnBottomFace = 3; //does not depend on type;
-      int tetraFaceIdOnUpperFace = 0; //does not depend on type;
-
-      int minIndex;
-      int prismT = prismType(prismVertexMpasIds, minIndex);
-
-      long long int reorderedPrismLIds[6];
-
-      for (int ii = 0; ii < 6; ii++)
-      {
-          reorderedPrismLIds[ii] = prismVertexGIds[PrismVerticesMap[minIndex][ii]];
-      }
-
-      for (int iTetra = 0; iTetra < 3; iTetra++)
-          for (int iVertex = 0; iVertex < 4; iVertex++)
-          {
-              tetrasIdsOnPrism[iTetra][iVertex] = reorderedPrismLIds[tetraOfPrism[prismT][iTetra][iVertex]];
-          }
-  }
-
-
-  void computeMap()
-  {
-    int PrismVerticesMap[6][6] = {{0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}, {3, 5, 4, 0, 2, 1}, {4, 3, 5, 1, 0, 2}, {5, 4, 3, 2, 1, 0}};
-
-    int tetraOfPrism[2][3][4] = {{{0, 1, 2, 5}, {0, 1, 5, 4}, {0, 4, 5, 3}}, {{0, 1, 2, 4}, {0, 4, 2, 5}, {0, 4, 5, 3}}};
-
-    int TetraFaces[4][3] = {{0 , 1 , 3}, {1 , 2 , 3}, {0 , 3 , 2}, {0 , 2 , 1}};
-
-    int PrismFaces[5][4] = {{0 , 1 , 4 , 3}, {1 , 2 , 5 , 4}, {0 , 3 , 5 , 2}, {0 , 2 , 1 , -1}, {3 , 4 , 5, -1}};
-
-
-    for(int pType=0; pType<2; ++pType){
-      std::cout<< "pType: " << pType <<std::endl;
-      for(int minIndex = 0; minIndex<6; ++minIndex){
-        std::cout<< "mIndex: " << minIndex <<std::endl;
-        for(int pFace = 0; pFace<5; ++pFace){
-          for(int tFace =0; tFace<4; ++tFace){
-            for(int iTetra =0; iTetra<3; ++iTetra){
-              int count=0;
-              for(int in =0; in<3; ++in){
-                int node=PrismVerticesMap[minIndex][tetraOfPrism[pType][iTetra][TetraFaces[tFace][in]]];
-                for(int i=0; i<4; ++i)
-                  count += (node == PrismFaces[pFace][i]);
-              }
-              if(count == 3)
-                std::cout << pFace << " " << tFace << " " << iTetra << std::endl;
-              }
-            }
-          }
-        }
-      std::cout<<std::endl;
-      }
-    std::cout<<std::endl;
-    }
-
-
-  void tetrasFromPrismStructured (int const* prismVertexMpasIds, int const* prismVertexGIds, int tetrasIdsOnPrism[][4])
-  {
-      int PrismVerticesMap[6][6] = {{0, 1, 2, 3, 4, 5}, {1, 2, 0, 4, 5, 3}, {2, 0, 1, 5, 3, 4}, {3, 5, 4, 0, 2, 1}, {4, 3, 5, 1, 0, 2}, {5, 4, 3, 2, 1, 0}};
-
-      int tetraOfPrism[2][3][4] = {{{0, 1, 2, 5}, {0, 1, 5, 4}, {0, 4, 5, 3}}, {{0, 1, 2, 4}, {0, 4, 2, 5}, {0, 4, 5, 3}}};
-
-      int minIndex = std::min_element (prismVertexMpasIds, prismVertexMpasIds + 3) - prismVertexMpasIds;
-
-      int v1 (prismVertexMpasIds[PrismVerticesMap[minIndex][1]]);
-      int v2 (prismVertexMpasIds[PrismVerticesMap[minIndex][2]]);
-
-      int prismType = v1  > v2;
-
-      for (int iTetra = 0; iTetra < 3; iTetra++)
-          for (int iVertex = 0; iVertex < 4; iVertex++)
-          {
-              tetrasIdsOnPrism[iTetra][iVertex] = prismVertexGIds[tetraOfPrism[prismType][iTetra][iVertex]];
-          }
-
-      // return;
-
-      int reorderedPrismLIds[6];
-
-      for (int ii = 0; ii < 6; ii++)
-      {
-          reorderedPrismLIds[ii] = prismVertexGIds[PrismVerticesMap[minIndex][ii]];
-      }
-
-      for (int iTetra = 0; iTetra < 3; iTetra++)
-          for (int iVertex = 0; iVertex < 4; iVertex++)
-          {
-              tetrasIdsOnPrism[iTetra][iVertex] = reorderedPrismLIds[tetraOfPrism[prismType][iTetra][iVertex]];
-          }
-  }
-
-
-
-
-  void setBdFacesOnPrism (const std::vector<std::vector<std::vector<int> > >& prismStruct, const std::vector<int>& prismFaceIds, std::vector<int>& tetraPos, std::vector<int>& facePos)
-  {
-    int numTriaFaces = prismFaceIds.size() - 2;
-    tetraPos.assign(numTriaFaces,-1);
-    facePos.assign(numTriaFaces,-1);
-
-
-    for (int iTetra (0), k (0); (iTetra < 3 && k < numTriaFaces); iTetra++)
-    {
-      bool found;
-      for (int jFaceLocalId = 0; jFaceLocalId < 4; jFaceLocalId++ )
-      {
-        found = true;
-        for (int ip (0); ip < 3 && found; ip++)
-        {
-          int localId = prismStruct[iTetra][jFaceLocalId][ip];
-          int j = 0;
-          found = false;
-          while ( (j < prismFaceIds.size()) && !found )
-          {
-            found = (localId == prismFaceIds[j]);
-            j++;
-          }
-        }
-        if (found)
-        {
-          tetraPos[k] = iTetra;
-          facePos[k] = jFaceLocalId;
-          k += found;
-          break;
-        }
-      }
-    }
-  }
-
   void procsSharingVertex(const int vertex, std::vector<int>& procIds) {
     int fCell = vertexToFCell[vertex];
     procIds.clear();
     int nEdg = nEdgesOnCells_F[fCell];
     int me;
     MPI_Comm_rank(comm, &me);
-    procIds.reserve(nEdg);
     for(int i=0; i<nEdg; ++i) {
       int fVertex = verticesOnCell_F[maxNEdgesOnCell_F * fCell + i]-1;
       if (verticesMask_F[fVertex] & dynamic_ice_bit_value) {
@@ -2135,60 +1580,63 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex)
     double const* observedThicknessTendency_F, double const * observedThicknessTendencyUncertainty_F) {
 
 
-//    std::vector<double> regulThk(thicknessData);
-//    for (int index = 0; index < nVertices; index++)
-//      regulThk[index] = std::max(1e-4, thicknessData[index]);  //TODO Make limit a parameter
-//
-//    importP0Temperature(temperature_F);
-//
-//    std::cout << "\n\nTimeStep: "<< *deltat << "\n\n"<< std::endl;
-//
-//    std::vector<int> mpasIndexToVertexID(nVertices);
-//    for (int i = 0; i < nVertices; i++) {
-//      mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]];
-//    }
-//
-//    mapVerticesToCells(velocityOnVertices, &velocityOnCells[0], 2, nLayers,
-//        Ordering);
-
-
     // Write out ASCII format
 
     std::cout << "Writing mesh to albany.msh." << std::endl;
     // msh file
     std::ofstream outfile;
-    outfile.open ("albany.msh", std::ios::out | std::ios::trunc);
+    outfile.precision(15);
+    std::stringstream name;
+    int me;
+    MPI_Comm_rank(comm, &me);
+    name <<  "albany.msh." <<  me;
+    outfile.open (name.str(), std::ios::out | std::ios::trunc);
     if (outfile.is_open()) {
 
        int nVerticesBoundaryEdge = 0;
-       for (int index = 0; index < nVertices; index++) {
+       for (int index = 0; index < nEdges; index++) {
           if (isBoundaryEdge[index]) nVerticesBoundaryEdge += 1;
        }
-       std::vector<int> verticesOnBoundaryEdge;
-       verticesOnBoundaryEdge.resize(2 * nVerticesBoundaryEdge);
+
+       //creating set from vector so that we can find elements in the set in log time
+       std::set<int> floatingEdgesSet;
+       for (int i = 0; i < floatingEdgesIds.size(); i++)
+         floatingEdgesSet.insert(floatingEdgesIds[i]);
+
+
+       std::vector<int> boundaryEdges; //list of edge vertices and edge label
+       boundaryEdges.resize(3 * nVerticesBoundaryEdge);
        int iVerticesBoundaryEdge = 0;
-       for (int index = 0; index < nVertices; index++) {
+       for (int index = 0; index < nEdges; index++) {
           if (isBoundaryEdge[index]) {
-             verticesOnBoundaryEdge[0 + 2 * iVerticesBoundaryEdge] = verticesOnEdge[0 + 2 * index];
-             verticesOnBoundaryEdge[1 + 2 * iVerticesBoundaryEdge] = verticesOnEdge[1 + 2 * index];
-             iVerticesBoundaryEdge += 1;
+             boundaryEdges[0 + 3 * iVerticesBoundaryEdge] = verticesOnEdge[0 + 2 * index];
+             boundaryEdges[1 + 3 * iVerticesBoundaryEdge] = verticesOnEdge[1 + 2 * index];
+             auto search = floatingEdgesSet.find(indexToEdgeID[index]);  //slow but executed only when printing
+             //boundary edges labels: 2 if floating, 1 otherwise
+             boundaryEdges[2 + 3 * iVerticesBoundaryEdge] = (search != floatingEdgesSet.end()) ? 2 : 1;
+             iVerticesBoundaryEdge ++;
           }
        }
-       //std::cout<<"final count: "<<iVerticesBoundaryEdge << "  nVerticesBoundaryEdge="<<nVerticesBoundaryEdge<<std::endl;
 
-       outfile << "Triangle " << 3 << "\n";  // first line saying it is a mesh of triangles
+       outfile << "Format: " << 1 << "\n";  // first line stating the format (we are going to provide global ids)
+       outfile << "Triangle " << 3 << "\n";  // second line saying it is a mesh of triangles
        outfile << nVertices << " " << nTriangles << " " << nVerticesBoundaryEdge << "\n";  // second line
 
        for (int index = 0; index < nVertices; index++) { //coordinates lines
           int iCell = vertexToFCell[index];
-          outfile << xCell_F[iCell] / unit_length << " " << yCell_F[iCell] / unit_length << " " << isVertexBoundary[index] << "\n"  ;
+          //setting boundary vertices labels, 2 for dirichlet nodes, 1 otherwise
+          bool isDirichletVertex = dirichletCellsMask_F[(nLayers+1)*iCell] != 0;
+          int vertexLabel = (!isVertexBoundary[index]) ? 0 :
+                            isDirichletVertex ? 3 : 1;
+
+          outfile << indexToVertexID[index] << " " << xCell_F[iCell] / unit_length << " " << yCell_F[iCell] / unit_length << " " << vertexLabel << "\n"  ;
        }
 
        for (int index = 0; index < nTriangles; index++) //triangles lines
-        outfile << verticesOnTria[0 + 3 * index] + 1 << " " << verticesOnTria[1 + 3 * index] + 1 << " " << verticesOnTria[2 + 3 * index] + 1 << " " << 1 << "\n"; // last digit can be used to specify a 'material'.  Not used by Albany LandIce, so giving dummy value
+        outfile << indexToTriangleID[index] << " " << verticesOnTria[0 + 3 * index] + 1 << " " << verticesOnTria[1 + 3 * index] + 1 << " " << verticesOnTria[2 + 3 * index] + 1 << " " << 1 << "\n"; // last digit can be used to specify a 'material'.  Not used by Albany LandIce, so giving dummy value
 
        for (int index = 0; index < nVerticesBoundaryEdge; index++) // boundary edges lines
-       outfile <<  verticesOnBoundaryEdge[0 + 2 * index] + 1 << " " << verticesOnBoundaryEdge[1 + 2 * index] + 1 << " " << 1 << "\n"; //last digit can be used to tell whether it's floating or not.. but let's worry about this later.
+       outfile <<  indexToEdgeID[index] << " " << boundaryEdges[0 + 3 * index] + 1 << " " << boundaryEdges[1 + 3 * index] + 1 << " " << boundaryEdges[2 + 3 * index] << "\n"; //last digit can be used to tell whether it's floating or not.. but let's worry about this later.
 
        outfile.close();
        }
@@ -2198,9 +1646,9 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex)
 
     // individual field values
     // Call needed functions to process MPAS fields to Albany units/format
-    
-    std::map<int, int> bdExtensionMap;  // local map to be created by import2DFields
-    import2DFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F, 
+
+    std::map<int, int> bdExtensionMap;  // local map to be created by importFields
+    importFields(bdExtensionMap, bedTopography_F, lowerSurface_F, thickness_F, beta_F,
                    stiffnessFactor_F, effecPress_F, temperature_F, smb_F, minThickness);
 
     import2DFieldsObservations(bdExtensionMap,
@@ -2235,20 +1683,23 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex)
     std::cout << "Writing temperature.ascii." << std::endl;
     outfile.open ("temperature.ascii", std::ios::out | std::ios::trunc);
     if (outfile.is_open()) {
-       outfile << nTriangles << " " << nLayers << "\n";  //number of triangles and number of layers on first line
-       
-       double midLayer = 0;
-       for (int il = 0; il < nLayers; il++) { //sigma coordinates for temperature
-         midLayer += layersRatio[il]/2.0;
-         outfile << midLayer << "\n";
-         midLayer += layersRatio[il]/2.0;
-       }    
-
-       for(int il = 0; il<nLayers; ++il)
-         for(int i = 0; i<nTriangles; ++i)  //temperature values layer by layer
-           outfile << temperatureData[i + il*nTriangles]<<"\n";
-       outfile.close();
-       }
+      outfile << nTriangles << " " << nLayers << "\n";  //number of triangles and number of layers on first line
+
+      double midLayer = 0;
+      for (int il = 0; il < nLayers; il++) { //sigma coordinates for temperature
+        midLayer += layersRatio[il]/2.0;
+        outfile << midLayer << "\n";
+        midLayer += layersRatio[il]/2.0;
+      }
+
+      int lElemColumnShift = (Ordering == 1) ? 1 : nTriangles;
+      int elemLayerShift = (Ordering == 0) ? 1 : nLayers;
+      for(int il = 0; il<nLayers; ++il)
+        for(int i = 0; i<nTriangles; ++i)  //temperature values layer by layer
+          outfile << temperatureDataOnPrisms[i*elemLayerShift + il*lElemColumnShift]<<"\n";
+
+      outfile.close();
+    }
     else {
        std::cout << "Failed to open tempertature ascii file!"<< std::endl;
     }
@@ -2282,6 +1733,7 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex)
     std::string filename = filenamebase+".ascii";
     std::cout << "Writing " << filename << std::endl;
     std::ofstream outfile;
+    outfile.precision(15);
     outfile.open (filename.c_str(), std::ios::out | std::ios::trunc);
     if (outfile.is_open()) {
        outfile << nVertices << "\n";  //number of vertices on first line
@@ -2311,4 +1763,4 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex)
        std::cout << "Error: Failed to open  "+filename << std::endl;
     }
   }
- 
+
diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.hpp b/src/core_landice/mode_forward/Interface_velocity_solver.hpp
index dbcf210bd8..2dccaf6e89 100644
--- a/src/core_landice/mode_forward/Interface_velocity_solver.hpp
+++ b/src/core_landice/mode_forward/Interface_velocity_solver.hpp
@@ -10,7 +10,6 @@ distributed with this code, or at http://mpas-dev.github.com/license.html
 // ===================================================
 //! Includes
 // ===================================================
-//#include <boost/program_options.hpp>
 #include <cstring>
 #include <vector>
 #include <mpi.h>
@@ -45,11 +44,6 @@ distributed with this code, or at http://mpas-dev.github.com/license.html
 #define write_ascii_mesh write_ascii_mesh_
 #endif
 
-//#include <lifev/core/algorithm/PreconditionerIfpack.hpp>
-//#include <lifev/core/algorithm/PreconditionerML.hpp>
-
-//#include <lifev/ice_sheet/solver/IceProblem.hpp>
-
 struct exchange {
   const int procID;
   const std::vector<int> vec;
@@ -112,6 +106,8 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
     int const* _verticesOnCell_F, int const* _verticesOnEdge_F,
     int const* _edgesOnCell_F, int const* _nEdgesOnCells_F,
     int const* _indexToCellID_F,
+    int const* _indexToEdgeID_F,
+    int const* _indexToVertexID_F,
     double const* _xCell_F, double const* _yCell_F, double const* _zCell_F,
     double const* _xVertex_F, double const* _yVertex_F, double const* _zVertex_F,
     double const* _areaTriangle_F,
@@ -119,15 +115,12 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
     int const* sendEdgesArray_F, int const* recvEdgesArray_F,
     int const* sendVerticesArray_F, int const* recvVerticesArray_F);
 
-void velocity_solver_extrude_3d_grid(double const* levelsRatio_F,
-    double const* lowerSurface_F, double const* thickness_F);
+void velocity_solver_extrude_3d_grid(double const* levelsRatio_F);
 
 void velocity_solver_export_l1l2_velocity();
 
 void velocity_solver_export_fo_velocity();
 
-//void velocity_solver_estimate_SS_SMB (const double* u_normal_F, double* sfcMassBal);
-
 void interface_init_log();
 
 void interface_redirect_stdout(int const* iTimestep);
@@ -148,24 +141,9 @@ void write_ascii_mesh(int const* indexToCellID_F,
 
 } // extern "C"
 
-extern int velocity_solver_init_mpi__(int* fComm);
+extern int velocity_solver_init_mpi__(MPI_Comm comm);
 extern void velocity_solver_finalize__();
 
-#ifdef LIFEV
-extern void velocity_solver_init_l1l2__(const std::vector<double>& layersRatio, const std::vector<double>& velocityOnVertices, bool initialize_velocity);
-
-extern void velocity_solver_solve_l1l2__(const std::vector<double>& elevationData,
-    const std::vector<double>& thicknessData, const std::vector<double>& betaData,
-    const std::vector<double>& temperatureData, const std::vector<int>& indexToVertexID,
-    std::vector<double>& velocityOnVertices);
-
-extern void velocity_solver_init_fo__(const std::vector<double>& layersRatio, const std::vector<double>& velocityOnVertices, const std::vector<int>& indexToVertexID, bool initialize_velocity);
-
-extern void velocity_solver_export_l1l2_velocity__(const std::vector<double>& layersRatio, const std::vector<double>& elevationData, const std::vector<double>& regulThk,
-    const std::vector<int>& mpasIndexToVertexID, MPI_Comm reducedComm);
-
-#endif
-
 extern void velocity_solver_set_physical_parameters__(double const& gravity, double const& ice_density, double const& ocean_density, double const& sea_level, double const& flowParamA, 
                         double const& flowLawExponent, double const& dynamic_thickness, bool const& useGLP, double const& clausiusClapeyronCoeff); 
 
@@ -182,31 +160,13 @@ extern void velocity_solver_solve_fo__(int nLayers, int nGlobalVertices,
     const std::vector<double>& smbData,
     const std::vector<double>& stiffnessFactorData,
     const std::vector<double>& effecPressData,
-    const std::vector<double>& temperatureOnTetra,
-    std::vector<double>& dissipationHeatOnTetra,
+    const std::vector<double>& temperatureDataOnPrisms,
+    std::vector<double>& dissipationHeatOnPrisms,
     std::vector<double>& velocityOnVertices,
     int& error,
     const double& deltat = 0.0);
 
-
-#ifdef LIFEV
-extern void velocity_solver_compute_2d_grid__(int nGlobalTriangles,
-     int nGlobalVertices, int nGlobalEdges,
-     const std::vector<int>& indexToVertexID,
-     const std::vector<double>& verticesCoords,
-     const std::vector<bool>& isVertexBoundary,
-     const std::vector<int>& verticesOnTria,
-     const std::vector<bool>& isBoundaryEdge,
-     const std::vector<int>& trianglesOnEdge,
-     const std::vector<int>& trianglesPositionsOnEdge,
-     const std::vector<int>& verticesOnEdge,
-     const std::vector<int>& indexToEdgeID,
-     const std::vector<int>& indexToTriangleID,
-     const std::vector < std::pair<int, int> >& procOnInterfaceEdge);
-
-#else
 extern void velocity_solver_compute_2d_grid__(MPI_Comm);
-#endif
 
 
 extern void velocity_solver_export_2d_data__(MPI_Comm reducedComm,
@@ -215,44 +175,33 @@ extern void velocity_solver_export_2d_data__(MPI_Comm reducedComm,
     const std::vector<double>& betaData,
     const std::vector<int>& indexToVertexID);
 
-extern void velocity_solver_extrude_3d_grid__(int nLayers, int nGlobalTriangles,
-    int nGlobalVertices, int nGlobalEdges, int Ordering, MPI_Comm reducedComm,
+extern void velocity_solver_extrude_3d_grid__(
+    int nLayers, int globalTriangleStride, int globalVertexStride, int globalEdgeStride,
+    int Ordering, MPI_Comm reducedComm,
     const std::vector<int>& indexToVertexID,
-    const std::vector<int>& mpasIndexToVertexID,
+    const std::vector<int>& vertexProcIDs,
     const std::vector<double>& verticesCoords,
-    const std::vector<bool>& isVertexBoundary,
     const std::vector<int>& verticesOnTria,
+    const std::vector<std::vector<int>> procsSharingVertices,
     const std::vector<bool>& isBoundaryEdge,
     const std::vector<int>& trianglesOnEdge,
-    const std::vector<int>& trianglesPositionsOnEdge,
     const std::vector<int>& verticesOnEdge,
     const std::vector<int>& indexToEdgeID,
     const std::vector<int>& indexToTriangleID,
     const std::vector<int>& dirichletNodes,
     const std::vector<int>&floatingEdges);
 
-//extern void velocity_solver_export_l1l2_velocity__();
-
 extern void velocity_solver_export_fo_velocity__(MPI_Comm reducedComm);
 
-
-#ifdef LIFEV
-extern int  velocity_solver_initialize_iceProblem__(bool keep_proc, MPI_Comm reducedComm);
-#endif
-
-//extern void velocity_solver_estimate_SS_SMB__ (const double* u_normal_F, double* sfcMassBal);
-
 exchangeList_Type unpackMpiArray(int const* array);
 
-bool isGhostTriangle(int i, double relTol = 1e-1);
-
 double signedTriangleArea(const double* x, const double* y);
 
 double signedTriangleArea(const double* x, const double* y, const double* z);
 
 void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id);
 
-void import2DFields(std::map<int, int> bdExtensionMap, double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F,
+void importFields(std::map<int, int> bdExtensionMap, double const* bedTopography_F, double const* lowerSurface_F, double const* thickness_F,
     double const* beta_F = 0, double const* stiffnessFactor_F = 0, double const* effecPress_F = 0, double const* temperature_F = 0, double const* smb_F = 0, double eps = 0);
 
 void import2DFieldsObservations(std::map<int, int> bdExtensionMap,
@@ -271,11 +220,6 @@ void write_ascii_mesh_field_int(std::vector<int> fieldData, std::string filename
 
 std::vector<int> extendMaskByOneLayer(int const* verticesMask_F);
 
-void extendMaskByOneLayer(int const* verticesMask_F,
-    std::vector<int>& extendedFVerticesMask);
-
-void importP0Temperature();
-
 void exportDissipationHeat(double * dissipationHeat_F);
 
 void get_prism_velocity_on_FEdges(double* uNormal,
@@ -284,15 +228,9 @@ void get_prism_velocity_on_FEdges(double* uNormal,
 
 int initialize_iceProblem(int nTriangles);
 
-void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F,
-    exchangeList_Type& receiveListReverse_F,
-    const std::vector<int>& fVertexToTriangleID,
-    const std::vector<int>& fCellToVertexID);
-
-void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F,
+void createReverseExchangeLists(exchangeList_Type& sendListReverse_F,
     exchangeList_Type& receiveListReverse_F,
-    const std::vector<int>& fVertexToTriangleID,
-    const std::vector<int>& fEdgeToEdgeID);
+    const std::vector<int>& newProcIds, const int* indexToID_F, exchangeList_Type const * recvList_F);
 
 void mapCellsToVertices(const std::vector<double>& velocityOnCells,
     std::vector<double>& velocityOnVertices, int fieldDim, int numLayers,
@@ -301,9 +239,6 @@ void mapCellsToVertices(const std::vector<double>& velocityOnCells,
 void mapVerticesToCells(const std::vector<double>& velocityOnVertices,
     double* velocityOnCells, int fieldDim, int numLayers, int ordering);
 
-void computeLocalOffset(int nLocalEntities, int& localOffset,
-    int& nGlobalEntities);
-
 void getProcIds(std::vector<int>& field, int const* recvArray);
 
 void getProcIds(std::vector<int>& field, exchangeList_Type const* recvList);
@@ -317,12 +252,6 @@ void allToAll(std::vector<int>& field, exchangeList_Type const* sendList,
 void allToAll(double* field, exchangeList_Type const* sendList,
     exchangeList_Type const* recvList, int fieldDim = 1);
 
-int prismType(long long int const* prismVertexMpasIds, int& minIndex);
-void tetrasFromPrismStructured (long long int const* prismVertexMpasIds, long long int const* prismVertexGIds, long long int tetrasIdsOnPrism[][4]);
-void computeMap();
-
-void setBdFacesOnPrism (const std::vector<std::vector<std::vector<int> > >& prismStruct, const std::vector<int>& prismFaceIds, std::vector<int>& tetraPos, std::vector<int>& facePos);
-void tetrasFromPrismStructured (int const* prismVertexMpasIds, int const* prismVertexGIds, int tetrasIdsOnPrism[][4]);
 void procsSharingVertex(const int vertex, std::vector<int>& procIds);
 
 bool belongToTria(double const* x, double const* t, double bcoords[3], double eps = 1e-3);
diff --git a/src/core_landice/mode_forward/mpas_li_advection.F b/src/core_landice/mode_forward/mpas_li_advection.F
index 4b728d3d29..4c0576ab89 100644
--- a/src/core_landice/mode_forward/mpas_li_advection.F
+++ b/src/core_landice/mode_forward/mpas_li_advection.F
@@ -257,7 +257,7 @@ subroutine li_advection_thickness_tracers(&
       call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
       call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness)
       call mpas_pool_get_array(geometryPool, 'layerThicknessEdge', layerThicknessEdge)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
       call mpas_pool_get_array(geometryPool, 'dynamicThickening', dynamicThickening)
       call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness)
@@ -630,7 +630,7 @@ subroutine li_update_geometry(geometryPool)
 
      call mpas_pool_get_dimension(geometryPool, 'nCells', nCells)
 
-     call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+     call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
      call mpas_pool_get_array(geometryPool, 'thickness', thickness)
      call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
      call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
diff --git a/src/core_landice/mode_forward/mpas_li_calving.F b/src/core_landice/mode_forward/mpas_li_calving.F
index 41e0523d8d..ac7d032ab0 100644
--- a/src/core_landice/mode_forward/mpas_li_calving.F
+++ b/src/core_landice/mode_forward/mpas_li_calving.F
@@ -402,7 +402,7 @@ subroutine li_restore_calving_front(domain, err)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
          call mpas_pool_get_array(geometryPool, 'restoreThickness', restoreThickness)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          ! get required fields from the thermal pool
          call mpas_pool_get_array(thermalPool, 'temperature', temperature)
@@ -870,7 +870,7 @@ subroutine floating_calving(domain, calvingFraction, err)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          calvingThickness = 0.0_RKIND
 
@@ -929,7 +929,7 @@ subroutine remove_small_islands(meshPool, geometryPool)
       call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
       call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
 
       do iCell = 1, nCellsSolve
@@ -1039,7 +1039,7 @@ subroutine topographic_calving(domain, calvingFraction, err)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
 
          calvingThickness = 0.0_RKIND
@@ -1152,7 +1152,7 @@ subroutine eigencalving(domain, err)
          call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
          call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
          call mpas_pool_get_array(geometryPool, 'eigencalvingParameter', eigencalvingParameter)
@@ -1336,7 +1336,7 @@ subroutine distribute_calving_flux(meshPool, geometryPool, scratchPool, calvingF
       ! get fields
       call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
       call mpas_pool_get_array(meshPool, 'deltat', deltat)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
       call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
       call mpas_pool_get_array(geometryPool, 'requiredCalvingVolumeRate', requiredCalvingVolumeRate)
@@ -1566,7 +1566,7 @@ subroutine calculate_calving_front_mask(meshPool, geometryPool, calvingFrontMask
       call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
       call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
       call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
 
       calvingFrontMask = 0 !initialize
@@ -1656,7 +1656,7 @@ subroutine remove_icebergs(domain)
          call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce)
          call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
          call mpas_pool_get_dimension(geometryPool, 'nCellsSolve', nCellsSolve)
@@ -1710,7 +1710,7 @@ subroutine remove_icebergs(domain)
             call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
             call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
             call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool)
-            call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+            call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
             call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce)
             call mpas_pool_get_array(geometryPool, 'contiguousDynamicIceOld', contiguousDynamicIceOld)
             call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
@@ -1765,7 +1765,7 @@ subroutine remove_icebergs(domain)
       block => domain % blocklist
       do while (associated(block))
          call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
          call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce)
diff --git a/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F b/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F
index ab47f0d226..a8f2664137 100644
--- a/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F
+++ b/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F
@@ -258,10 +258,9 @@ subroutine li_calculate_apparent_diffusivity(meshPool, velocityPool, scratchPool
       call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel=1)
       call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness, timeLevel=1)
       call mpas_pool_get_array(geometryPool, 'apparentDiffusivity', apparentDiffusivity)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(velocityPool, 'uReconstructZonal', uReconstructAxis1)
       call mpas_pool_get_array(velocityPool, 'uReconstructMeridional', uReconstructAxis2)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
       call mpas_pool_get_field(scratchPool, 'workCell', slopeReconstructXField)
       call mpas_allocate_scratch_field(slopeReconstructXField, .true.)
       call mpas_pool_get_field(scratchPool, 'workCell2', slopeReconstructYField)
diff --git a/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F b/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F
index 95f6716fa8..0b99e9d843 100644
--- a/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F
+++ b/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F
@@ -247,7 +247,7 @@ subroutine li_basal_melt_floating_ice(domain, err)
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
 
          ! get fields from the scratch pool
@@ -853,7 +853,7 @@ subroutine calc_iceshelf_draft_info(domain, GLdepth, CFdepth, err)
          call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
          call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
          call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
diff --git a/src/core_landice/mode_forward/mpas_li_sia.F b/src/core_landice/mode_forward/mpas_li_sia.F
index b22001c4f2..d400c12e35 100644
--- a/src/core_landice/mode_forward/mpas_li_sia.F
+++ b/src/core_landice/mode_forward/mpas_li_sia.F
@@ -316,7 +316,7 @@ subroutine li_sia_solve(meshPool, geometryPool, thermalPool, velocityPool, err)
       call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'slopeEdge', slopeEdge)
       call mpas_pool_get_array(geometryPool, 'normalSlopeEdge', normalSlopeEdge)
       call mpas_pool_get_array(thermalPool, 'temperature', temperature)
diff --git a/src/core_landice/mode_forward/mpas_li_statistics.F b/src/core_landice/mode_forward/mpas_li_statistics.F
index 47a7be44e8..99de7c4939 100644
--- a/src/core_landice/mode_forward/mpas_li_statistics.F
+++ b/src/core_landice/mode_forward/mpas_li_statistics.F
@@ -235,7 +235,7 @@ subroutine li_compute_statistics(domain, itimestep)
         ! Geometry arrays
         call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
         call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
-        call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+        call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
         call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
         call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel = 1)
         call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness, timeLevel = 1)
diff --git a/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F b/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F
index 9e2e4f79ca..596119d01c 100644
--- a/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F
+++ b/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F
@@ -163,7 +163,7 @@ subroutine li_SGH_init(domain, err)
          waterPressure = max(0.0_RKIND, waterPressure)
          waterPressure = min(waterPressure, rhoi * gravity * thickness)
          ! set pressure correctly under floating ice and open ocean
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
          where ( (li_mask_is_floating_ice(cellMask)) .or. &
                  ( (.not. li_mask_is_ice(cellMask)) .and. (bedTopography < config_sea_level) ) )
@@ -379,7 +379,7 @@ subroutine li_SGH_solve(domain, err)
          call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
          call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput)
          call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          WtillOld = Wtill
          Wtill = Wtill + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - Cd)
@@ -420,7 +420,7 @@ subroutine li_SGH_solve(domain, err)
          call mpas_pool_get_array(hydroPool, 'waterVelocity', waterVelocity)
          call mpas_pool_get_array(hydroPool, 'waterVelocityCellX', waterVelocityCellX)
          call mpas_pool_get_array(hydroPool, 'waterVelocityCellY', waterVelocityCellY)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          allocate(cellJunk(nCells+1))
          call mpas_reconstruct(meshPool, waterVelocity, waterVelocityCellX, waterVelocityCellY, &
               cellJunk, cellJunk, cellJunk)
@@ -476,7 +476,7 @@ subroutine li_SGH_solve(domain, err)
          call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          divergence(:) = 0.0_RKIND  ! zero div before starting
          ! loop over locally owned cells
@@ -551,7 +551,7 @@ subroutine li_SGH_solve(domain, err)
          call mpas_pool_get_array(hydroPool, 'divergence', divergence)
          call mpas_pool_get_array(hydroPool, 'divergenceChannel', divergenceChannel)
          call mpas_pool_get_array(hydroPool, 'channelAreaChangeCell', channelAreaChangeCell)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
          waterThicknessOld = waterThickness
          waterThickness = waterThicknessOld + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - divergence  &
@@ -755,7 +755,7 @@ subroutine calc_edge_quantities(block, err)
       call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeNormal', hydropotentialBaseSlopeNormal)
       call mpas_pool_get_array(hydroPool, 'waterPressureSlopeNormal', waterPressureSlopeNormal)
       call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
       call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell)
       call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeTangent', hydropotentialBaseSlopeTangent)
@@ -1217,7 +1217,7 @@ subroutine calc_pressure(block, err)
       call mpas_pool_get_array(velocityPool, 'basalSpeed', basalSpeed)
       call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
 
       openingRate(:) = bedRough * basalSpeed(:) * (bedRoughMax - waterThickness(:))
@@ -1338,7 +1338,7 @@ subroutine calc_pressure_diag_vars(block, err)
       call mpas_pool_get_array(hydroPool, 'hydropotentialBase', hydropotentialBase)
       call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness)
       call mpas_pool_get_array(hydroPool, 'hydropotential', hydropotential)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
       effectivePressure = rhoi * gravity * thickness - waterPressure
          ! < this should evalute to 0 for floating ice if Pw set correctly there.
diff --git a/src/core_landice/mode_forward/mpas_li_thermal.F b/src/core_landice/mode_forward/mpas_li_thermal.F
index 383ba22dae..6a4807a1dd 100644
--- a/src/core_landice/mode_forward/mpas_li_thermal.F
+++ b/src/core_landice/mode_forward/mpas_li_thermal.F
@@ -831,7 +831,7 @@ subroutine li_thermal_solver(domain, err)
          ! get fields from the geometry pool
          call mpas_pool_get_array(geometryPool, 'thickness', thickness)
          call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'groundedBasalMassBal', groundedBasalMassBal)
          call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
          call mpas_pool_get_array(geometryPool, 'basalWaterThickness', basalWaterThickness)
diff --git a/src/core_landice/mode_forward/mpas_li_velocity.F b/src/core_landice/mode_forward/mpas_li_velocity.F
index 4f5f6a2820..59a7ed5e2e 100644
--- a/src/core_landice/mode_forward/mpas_li_velocity.F
+++ b/src/core_landice/mode_forward/mpas_li_velocity.F
@@ -268,8 +268,9 @@ subroutine li_velocity_solve(domain, solveVelo, err)
       integer, pointer :: nEdges
       integer, pointer :: nVertInterfaces
       integer, pointer :: nCells
-      integer, dimension(:), pointer :: edgeMask, cellMask, vertexMask, vertexMaskOld
+      integer, dimension(:), pointer :: edgeMask, cellMask, cellMaskOld, vertexMask, vertexMaskOld
       integer, dimension(:,:), pointer :: dirichletVelocityMaskOld, dirichletVelocityMaskNew
+      integer, dimension(:), pointer :: floatingEdgesOld, floatingEdgesNew
       real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, normalVelocityInitial
       real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ, &
          uReconstructZonal, uReconstructMeridional
@@ -287,7 +288,9 @@ subroutine li_velocity_solve(domain, solveVelo, err)
       real (kind=RKIND), dimension(:), pointer :: upperSurface
       integer, dimension(:), pointer :: indexToEdgeID
       integer, pointer :: anyDynamicVertexMaskChanged
+      integer, pointer :: cellIceMaskChanged
       integer, pointer :: dirichletMaskChanged
+      integer, pointer :: floatingEdgesChanged
       ! truly local variables
       integer :: cell1, cell2
       integer :: cell3, cell4, thisCell
@@ -300,7 +303,9 @@ subroutine li_velocity_solve(domain, solveVelo, err)
       real (kind=RKIND) :: maxThicknessOnProc, maxThicknessAllProcs
       real (kind=RKIND) :: xVelEdge, yVelEdge
       integer :: blockDynamicVertexMaskChanged, procDynamicVertexMaskChanged
+      integer :: blockCellIceMaskChanged, procCellIceMaskChanged
       integer :: blockDirichletMaskChanged, procDirichletMaskChanged
+      integer :: blockFloatingEdgesChanged, procFloatingEdgesChanged
 
       call mpas_timer_start("velocity solve")
 
@@ -385,7 +390,9 @@ subroutine li_velocity_solve(domain, solveVelo, err)
          ! Don't bother checking this with SIA because it requires an extra global reduce
          maxThicknessOnProc = 0.0_RKIND  ! initialize to
          procDynamicVertexMaskChanged = 0
+         procCellIceMaskChanged = 0
          procDirichletMaskChanged = 0
+         procFloatingEdgesChanged = 0
 
          block => domain % blocklist
          do while (associated(block))
@@ -397,7 +404,7 @@ subroutine li_velocity_solve(domain, solveVelo, err)
             maxThicknessOnProc = max(maxThicknessOnProc, maxval(thickness))
 
             ! The interface expects an array where 1's are floating edges and 0's are non-floating edges.
-            call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges)
+            call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1)
             call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
             floatingEdges = li_mask_is_floating_ice_int(edgeMask)
             call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1)
@@ -407,7 +414,7 @@ subroutine li_velocity_solve(domain, solveVelo, err)
             ! TODO:  There may be some aspects of the mask that are ok change for external dycores,
             !        but for now just check the whole thing.
             call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMaskOld, timeLevel=2)
-            if ( sum(li_mask_is_dynamic_ice_int(vertexMask) - li_mask_is_dynamic_ice_int(vertexMaskOld)) /= 0 ) then
+            if ( sum(abs(li_mask_is_dynamic_ice_int(vertexMask) - li_mask_is_dynamic_ice_int(vertexMaskOld))) /= 0 ) then
                 blockDynamicVertexMaskChanged = 1
             else
                 blockDynamicVertexMaskChanged = 0
@@ -416,11 +423,26 @@ subroutine li_velocity_solve(domain, solveVelo, err)
             ! Determine if any blocks on this processor had a change to the vertex mask
             procDynamicVertexMaskChanged = max(procDynamicVertexMaskChanged, blockDynamicVertexMaskChanged)
             !print *,'procVertexMaskChanged', procVertexMaskChanged
+            
+            ! Determine if the cell mask changed during this time step for this block (needed for external dycores)
+            ! TODO:  There may be some aspects of the mask that are ok change for external dycores,
+            !        but for now just check the whole thing.
+            call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1)
+            call mpas_pool_get_array(geometryPool, 'cellMask', cellMaskOld, timeLevel=2)
+            if (( sum(abs(li_mask_is_dynamic_ice_int(cellMask) - li_mask_is_dynamic_ice_int(cellMaskOld))) /= 0 ) .or. &
+                ( sum(abs(li_mask_is_ice_int(cellMask) - li_mask_is_ice_int(cellMaskOld))) /= 0 )) then
+                blockCellIceMaskChanged = 1
+            else
+                blockCellIceMaskChanged = 0
+            endif
+      
+            ! Determine if any blocks on this processor had a change to the cell mask
+            procCellIceMaskChanged = max(procCellIceMaskChanged, blockCellIceMaskChanged)
 
-            ! Also check to see if the Dirichlet b.c. mask has changed
+            ! Check to see if the Dirichlet b.c. mask has changed
             call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMaskOld, timeLevel=2)
             call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMaskNew, timeLevel=1)
-            if ( sum(dirichletVelocityMaskNew - dirichletVelocityMaskOld) /= 0 ) then
+            if ( sum(abs(dirichletVelocityMaskNew - dirichletVelocityMaskOld)) /= 0 ) then
                 blockDirichletMaskChanged = 1
             else
                 blockDirichletMaskChanged = 0
@@ -428,6 +450,17 @@ subroutine li_velocity_solve(domain, solveVelo, err)
             ! Determine if any blocks on this processor had a change to the vertex mask
             procDirichletMaskChanged = max(procDirichletMaskChanged, blockDirichletMaskChanged)
 
+            ! Check to see if the Floating Edges Ids have changed
+            call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdgesOld, timeLevel=2)
+            call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdgesNew, timeLevel=1)
+            if ( sum(abs(floatingEdgesNew - floatingEdgesOld)) /= 0 ) then
+                blockFloatingEdgesChanged = 1
+            else
+                blockFloatingEdgesChanged = 0
+            endif
+            ! Determine if any blocks on this processor had a change to the vertex mask
+            procFloatingEdgesChanged = max(procFloatingEdgesChanged, blockFloatingEdgesChanged)
+
             ! Set beta for solver to use
             ! (this could potentially be applied to an SIA solver, so it is calculated in this module)
             call calculate_beta(block, err_tmp)
@@ -451,10 +484,16 @@ subroutine li_velocity_solve(domain, solveVelo, err)
          call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged)
          call mpas_dmpar_max_int(domain % dminfo, procDynamicVertexMaskChanged, anyDynamicVertexMaskChanged)
          !print *,'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged
+         ! Do the same for the cell mask
+         call mpas_pool_get_array(velocityPool, 'cellIceMaskChanged', cellIceMaskChanged)
+         call mpas_dmpar_max_int(domain % dminfo, procCellIceMaskChanged, cellIceMaskChanged)
          ! Do the same for the Dirichlet b.c. mask
          call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged)
          call mpas_dmpar_max_int(domain % dminfo, procDirichletMaskChanged, dirichletMaskChanged)
-         !print *,'dirichletMaskChanged', dirichletMaskChanged
+         ! Do the same for the Floating edges
+         call mpas_pool_get_array(velocityPool, 'floatingEdgesChanged', floatingEdgesChanged)
+         call mpas_dmpar_max_int(domain % dminfo, procFloatingEdgesChanged, floatingEdgesChanged)
+         !print *,'floatingEdgesChanged', floatingEdgesChanged
 
       endif
 
@@ -478,7 +517,7 @@ subroutine li_velocity_solve(domain, solveVelo, err)
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex)
          call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge)
-         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
          call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
 
 
@@ -823,7 +862,7 @@ subroutine calculate_beta(block, err)
 
       call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve)
       call mpas_pool_get_array(velocityPool, 'beta', beta)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
 
       betaSolve = beta
 
@@ -999,7 +1038,7 @@ subroutine remove_fast_ice(meshPool, uReconstructX, uReconstructY, normalVelocit
       call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
       call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
       call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'thickness', thickness)
       call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
       call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
diff --git a/src/core_landice/mode_forward/mpas_li_velocity_external.F b/src/core_landice/mode_forward/mpas_li_velocity_external.F
index a214560efe..f1fa8bd529 100644
--- a/src/core_landice/mode_forward/mpas_li_velocity_external.F
+++ b/src/core_landice/mode_forward/mpas_li_velocity_external.F
@@ -327,7 +327,7 @@ subroutine li_velocity_external_block_init(block, err)
       call velocity_solver_set_grid_data(nCells, nEdges, nVertices, nVertInterfaces, &
               nCellsSolve, nEdgesSolve, nVerticesSolve, maxNEdgesOnCell, radius, &
               cellsOnEdge, cellsOnVertex, verticesOnCell, verticesOnEdge, edgesOnCell, &
-              nEdgesOnCell, indexToCellID, &
+              nEdgesOnCell, indexToCellID, indexToEdgeID, indexToVertexID, &
               xCell, yCell, zCell, xVertex, yVertex, zVertex, areaTriangle, &
               sendCellsArray, recvCellsArray, &
               sendEdgesArray, recvEdgesArray, &
@@ -452,7 +452,9 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
       logical, pointer :: config_output_external_velocity_solver_data
       real (kind=RKIND), pointer :: config_ice_density
       integer, pointer :: anyDynamicVertexMaskChanged
+      integer, pointer :: cellIceMaskChanged
       integer, pointer :: dirichletMaskChanged
+      integer, pointer :: floatingEdgesChanged
       integer, pointer :: nEdges
       integer, pointer :: nVertLevels
       integer, pointer :: timestepNumber
@@ -486,7 +488,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
       call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface)
       call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface)
       call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
       call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
 
@@ -501,9 +503,11 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
       call mpas_pool_get_array(velocityPool, 'uReconstructZ', uReconstructZ)
       call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve)
       call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged)
+      call mpas_pool_get_array(velocityPool, 'cellIceMaskChanged', cellIceMaskChanged)
       call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged)
       call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1)
-      call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges)
+      call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1)
+      call mpas_pool_get_array(velocityPool, 'floatingEdgesChanged', floatingEdgesChanged)
       call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)
 
       ! Hydro variables
@@ -521,7 +525,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro
       ! Note these functions will always be called on the first solve because we
       ! initialize vertexMask to garbage which sets anyDynamicVertexMaskChanged to 1.
       if ((anyDynamicVertexMaskChanged == 1) .or. (config_always_compute_fem_grid) .or. &
-          (dirichletMaskChanged == 1) ) then
+          (dirichletMaskChanged == 1) .or. (floatingEdgesChanged == 1) .or. (cellIceMaskChanged == 1)) then
          call mpas_log_write("Generating new external velocity solver FEM grid.", flushNow=.true.)
          call generate_fem_grid(config_velocity_solver, vertexMask, cellMask, dirichletVelocityMask, &
               floatingEdges, layerThicknessFractions, lowerSurface, thickness, err)
@@ -779,10 +783,6 @@ subroutine li_velocity_external_write_albany_mesh(domain)
          call mpas_log_write("config_velocity solver needs to be set to 'FO' for config_write_albany_ascii_mesh to work.", &
             MPAS_LOG_CRIT)
       endif
-      ! check nProcs
-      if (domain % dminfo % nProcs /= 1) then
-         call mpas_log_write("config_write_albany_ascii_mesh currently only works on 1 processor.", MPAS_LOG_CRIT)
-      endif
       ! check nBlocks
       if (domain % dminfo % total_blocks /= 1) then
          call mpas_log_write("config_write_albany_ascii_mesh currently only works on 1 block per processor.", MPAS_LOG_CRIT)
@@ -809,13 +809,13 @@ subroutine li_velocity_external_write_albany_mesh(domain)
       call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
       call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
       call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1)
-      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+      call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1)
       call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
       call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness)
 
       ! Velocity variables
       call mpas_pool_get_array(velocityPool, 'beta', beta)
-      call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges)
+      call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1)
       call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1)
       call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)
 
@@ -1022,7 +1022,7 @@ subroutine generate_fem_grid(config_velocity_solver, vertexMask, cellMask, diric
           case ('FO') ! ===============================================
 #ifdef USE_EXTERNAL_FIRSTORDER
               call mpas_timer_start("velocity_solver_extrude_3d_grid")
-              call velocity_solver_extrude_3d_grid(layerThicknessFractions, lowerSurface, thickness)
+              call velocity_solver_extrude_3d_grid(layerThicknessFractions)
               call mpas_timer_stop("velocity_solver_extrude_3d_grid")
               call mpas_timer_start("velocity_solver_init_FO")
               call velocity_solver_init_FO(layerThicknessFractions)
diff --git a/testing_and_setup/compass/landice/MISMIP+/albany_input.yaml b/testing_and_setup/compass/landice/MISMIP+/albany_input.yaml
index c89dcb14ce..9676e6bbb9 100644
--- a/testing_and_setup/compass/landice/MISMIP+/albany_input.yaml
+++ b/testing_and_setup/compass/landice/MISMIP+/albany_input.yaml
@@ -35,7 +35,7 @@ ANONYMOUS:
     Body Force: 
       Type: FO INTERP SURF GRAD
   Discretization: 
-    Cubature Degree: 1
+    Element Shape: Prism
   Piro: 
     LOCA: 
       Bifurcation: { }
diff --git a/testing_and_setup/compass/landice/Thwaites_variability/albany_input.yaml b/testing_and_setup/compass/landice/Thwaites_variability/albany_input.yaml
index af2a6dd48c..a41f919c8d 100644
--- a/testing_and_setup/compass/landice/Thwaites_variability/albany_input.yaml
+++ b/testing_and_setup/compass/landice/Thwaites_variability/albany_input.yaml
@@ -31,8 +31,8 @@ ANONYMOUS:
         Type: Lateral
         Immersed Ratio: 8.85000000000000009e-01
   Discretization: 
-    Cubature Degree: 11
-    Interleaved Ordering: true
+    Element Shape: Tetrahedron
+    #Exodus Output File Name: albany_output.exo
   Piro: 
     LOCA: 
       Bifurcation: { }
diff --git a/testing_and_setup/compass/landice/circular-shelf/albany_input.yaml b/testing_and_setup/compass/landice/circular-shelf/albany_input.yaml
index 029369ad0c..3ac436fb05 100644
--- a/testing_and_setup/compass/landice/circular-shelf/albany_input.yaml
+++ b/testing_and_setup/compass/landice/circular-shelf/albany_input.yaml
@@ -17,8 +17,6 @@ ANONYMOUS:
       Type: FO INTERP SURF GRAD
   Discretization: 
     Exodus Output File Name: albany_output.exo
-    Cubature Degree: 1
-    Interleaved Ordering: false
   Piro: 
     LOCA: 
       Bifurcation: { }
diff --git a/testing_and_setup/compass/landice/confined-shelf/albany_input.yaml b/testing_and_setup/compass/landice/confined-shelf/albany_input.yaml
index 52809be2bf..1c75c642b0 100644
--- a/testing_and_setup/compass/landice/confined-shelf/albany_input.yaml
+++ b/testing_and_setup/compass/landice/confined-shelf/albany_input.yaml
@@ -17,8 +17,6 @@ ANONYMOUS:
       Type: FO INTERP SURF GRAD
   Discretization: 
     Exodus Output File Name: albany_output.exo
-    Cubature Degree: 1
-    Interleaved Ordering: false
   Piro: 
     LOCA: 
       Bifurcation: { }
diff --git a/testing_and_setup/compass/landice/dome/albany_input.yaml b/testing_and_setup/compass/landice/dome/albany_input.yaml
index 01ec65683a..3357ddccc3 100644
--- a/testing_and_setup/compass/landice/dome/albany_input.yaml
+++ b/testing_and_setup/compass/landice/dome/albany_input.yaml
@@ -11,6 +11,7 @@ ANONYMOUS:
       'Glen''s Law Homotopy Parameter': 1.0
   Discretization: 
     Exodus Output File Name: albany_output.exo
+    Element Shape: Prism
   Piro: 
     LOCA: 
       Bifurcation: { }
diff --git a/testing_and_setup/compass/landice/greenland/albany_input.yaml b/testing_and_setup/compass/landice/greenland/albany_input.yaml
index d2fb7eab2a..67b76efabe 100644
--- a/testing_and_setup/compass/landice/greenland/albany_input.yaml
+++ b/testing_and_setup/compass/landice/greenland/albany_input.yaml
@@ -17,8 +17,7 @@ ANONYMOUS:
       Type: FO INTERP SURF GRAD
   Discretization: 
     Exodus Output File Name: albany_output.exo
-    Cubature Degree: 1
-    Interleaved Ordering: false
+    Element Shape: Tetrahedron
   Piro: 
     LOCA: 
       Bifurcation: { }
diff --git a/testing_and_setup/compass/landice/greenland/albany_schoof_input.yaml b/testing_and_setup/compass/landice/greenland/albany_schoof_input.yaml
index 3e1fa88dc5..045feb0c49 100644
--- a/testing_and_setup/compass/landice/greenland/albany_schoof_input.yaml
+++ b/testing_and_setup/compass/landice/greenland/albany_schoof_input.yaml
@@ -39,8 +39,6 @@ ANONYMOUS:
         Immersed Ratio: 8.8521e-01
   Discretization: 
     Exodus Output File Name: albany_output.exo
-    Cubature Degree: 1
-    Interleaved Ordering: false
   Piro: 
     LOCA: 
       Bifurcation: { }