diff --git a/m3dc1_scorec/api/m3dc1_scorec.cc b/m3dc1_scorec/api/m3dc1_scorec.cc
index be0289a22..02ab5dcfd 100644
--- a/m3dc1_scorec/api/m3dc1_scorec.cc
+++ b/m3dc1_scorec/api/m3dc1_scorec.cc
@@ -7,36 +7,694 @@
   BSD license as described in the LICENSE file in the top-level directory.
  
 *******************************************************************************/
-#include "m3dc1_scorec.h"
-#include "m3dc1_matrix.h"
-#include "m3dc1_model.h"
-#include "m3dc1_mesh.h"
-#include "m3dc1_field.h"
 #include <mpi.h>
-#include <PCU.h>
-#include "gmi_null.h" // FIXME: should be deleted later on since it's added temporarily for null model
-#include <gmi_analytic.h>
 #include <map>
 #include <cstring>
 #include <iomanip> // setprecision
 #include <fstream> // file input
-#include "apfMDS.h"
+// #include <stdio.h>
+// #include <stdlib.h>
+
+// headers from SCOREC/Core
+#include <PCU.h>
+#include <pcu_util.h>
+#include <lionPrint.h>
+#include <gmi_null.h> // FIXME: should be deleted later on since it's added temporarily for null model
+#include <gmi_analytic.h>
+#include <apfMDS.h>
+#include <apf.h>
+#include <apfShape.h>
+#include <pumi.h>
+#include <spr.h>
+
+// local headers
+#include "m3dc1_scorec.h"
+#include "m3dc1_matrix.h"
+#include "m3dc1_model.h"
+#include "m3dc1_mesh.h"
+#include "m3dc1_field.h"
 #include "Expression.h"
 #include "m3dc1_slnTransfer.h"
 #include "m3dc1_sizeField.h"
 #include "ReducedQuinticImplicit.h"
-#include "pumi.h"
-// #include <stdio.h>
-// #include <stdlib.h>
+
 #ifdef M3DC1_TRILINOS
 #include "m3dc1_ls.h"
 #endif
 #include <alloca.h>
 
+const int dofNode = C1TRIDOFNODE;
+
 #ifdef DEBUG
-int begin_numVert;
+static const char* get_field_name_from_id(const FieldID id)
+{
+   m3dc1_field* mf = (*m3dc1_mesh::instance()->field_container)[id];
+   return getName(mf->get_field());
+}
 #endif
 
+
+// static functions used for spr-adapt
+static apf::Field* get_field_at_index(apf::Mesh2* m, apf::Field* inField, int index, int numDofs)
+{
+  int numComps = apf::countComponents(inField);
+  int numFields = numComps/numDofs;
+  PCU_ALWAYS_ASSERT(index <= numFields);
+  PCU_ALWAYS_ASSERT(index > 0);
+
+  apf::Field* targetField = apf::createPackedField(m, "target_field", numDofs);
+
+  apf::NewArray<double> allDofs(numComps);
+  apf::MeshEntity* v;
+  apf::MeshIterator* it = m->begin(0);
+  while ( (v = m->iterate(it)) )
+  {
+    apf::getComponents(inField, v, 0, &(allDofs[0]));
+    apf::setComponents(targetField, v, 0, &(allDofs[(index-1)*numDofs]));
+  }
+  m->end(it);
+  return targetField;
+}
+
+// Get a scaler field (single or scaler combination of any of six DOFs)
+static apf::Field* get_scalerComponent_of_field(apf::Mesh2* m, apf::Field* in, int comp)
+{
+  int comps = countComponents(in);
+  PCU_ALWAYS_ASSERT(comp < comps);
+  apf::Field* out = apf::createFieldOn(m, "in_field_comp", apf::SCALAR);
+
+  apf::MeshEntity* e;
+  apf::MeshIterator* it = m->begin(0);
+
+  double dofs[FIXSIZEBUFF];
+  while ( (e = m->iterate(it)) )
+  {
+    getComponents(in, e, 0, &dofs[0]);
+    double dof_comp = dofs[comp];
+    setComponents(out, e, 0, &dof_comp);
+  }
+  m->end(it);
+  return out;
+}
+
+// Get a vector of desired field 
+static apf::Field* get_vectorComponent_of_field(apf::Mesh2* m, apf::Field* in)
+{
+  apf::Field* out = apf::createFieldOn(m, "in_field_comp", apf::VECTOR);
+  apf::MeshEntity* e;
+  apf::MeshIterator* it = m->begin(0);
+
+  double dofs[FIXSIZEBUFF];
+  while ( (e = m->iterate(it)) )
+  {
+    getComponents(in, e, 0, &dofs[0]);
+    double comp_dR = dofs[1];
+    double comp_dZ = dofs[2];
+    double comp_dR2 = dofs[3];
+    double comp_dZ2 = dofs[5];
+    double vField[3] = {comp_dR, comp_dZ, 0.0};
+    setComponents(out, e, 0, &vField[0]);
+  }
+  m->end(it);
+  return out;
+}
+
+static apf::Field* get_ip_field(apf::Mesh2* m, apf::Field* in)
+{
+  ReducedQuinticImplicit shape;
+  int numComps = apf::countComponents(in);
+  assert(numComps == dofNode);
+  int dim = m->getDimension();
+  assert(dim == 2);
+  int order = 2;
+  apf::Field* ip = apf::createIPField(m, "ip_field", apf::VECTOR, order);
+
+  apf::MeshEntity* e;
+  apf::MeshIterator* it = m->begin(dim);
+
+  while ( (e = m->iterate(it)) )
+  {
+    // setup the ReducedQuintic Related Info
+    apf::MeshEntity* dvs[3];
+    int nd = m->getDownward(e, 0, dvs);
+    double coords[3][2];
+    for (int i = 0; i < 3; i++) {
+      apf::Vector3 p;
+      m->getPoint(dvs[i], 0, p);
+      coords[i][0] = p[0];
+      coords[i][1] = p[1];
+    }
+    shape.setCoord(coords);
+
+    apf::NewArray<double> values(3*dofNode);
+    apf::NewArray<double> dofAtXi(dofNode);
+    for (int i = 0; i < 3; i++)
+      apf::getComponents(in, dvs[i], 0, &(values[dofNode*i]));
+
+    shape.setDofs(&(values[0]));
+
+
+    apf::MeshElement* me = apf::createMeshElement(m, e);
+    for (int i = 0; i < apf::countIntPoints(me, order); i++) {
+      apf::Vector3 xi; // parametric coords of the point in e at which we are evaluating the field
+      apf::getIntPoint(me, order, i, xi);
+      apf::Vector3 p; // physical coords of the point in e at which we are evaluating the field
+      apf::mapLocalToGlobal(me, xi, p);
+      double pArray[3];
+      p.toArray(pArray);
+      shape.eval_g(pArray, &(dofAtXi[0]));
+      apf::Vector3 grad(dofAtXi[1], dofAtXi[2], 0.0);
+      apf::setVector(ip, e, i, grad);
+    }
+    apf::destroyMeshElement(me);
+  }
+  m->end(it);
+  return ip;
+}
+
+static void process_size_field(apf::Mesh2* m, apf::Field* in_size, int ts,
+    double max_size, int refine_level, int coarsen_level)
+{
+  //compute both average and min current size at each vertex
+  apf::Field* sum_field = apf::createFieldOn(m, "sum_field", apf::SCALAR);
+  apf::Field* min_field = apf::createFieldOn(m, "min_field", apf::SCALAR);
+  apf::Field* cnt_field = apf::createFieldOn(m, "cnt_field", apf::SCALAR);
+
+  apf::MeshEntity* v;
+  apf::MeshIterator* it = m->begin(0);
+  while ( (v = m->iterate(it)) )
+  {
+    double current_min = 1.e32;
+    double current_sum = 0.;
+    double current_cnt = 0.;
+    for (int i = 0; i < m->countUpward(v); i++) {
+      double edge_length = apf::measure(m, m->getUpward(v, i));
+      if (edge_length < current_min)
+      	current_min = edge_length;
+      current_sum += edge_length;
+      current_cnt += 1.;
+    }
+    apf::setScalar(sum_field, v, 0, current_sum);
+    apf::setScalar(min_field, v, 0, current_min);
+    apf::setScalar(cnt_field, v, 0, current_cnt);
+  }
+  m->end(it);
+  // accumulate sum and cnt fields and compute the averages
+  apf::accumulate(sum_field);
+  apf::accumulate(cnt_field);
+  it = m->begin(0);
+  while ( (v = m->iterate(it)) )
+  {
+    double total_sum = apf::getScalar(sum_field, v, 0);
+    double total_cnt = apf::getScalar(cnt_field, v, 0);
+    apf::setScalar(sum_field, v, 0, total_sum/total_cnt);
+  }
+  m->end(it);
+
+  // update the min field using share reduction with min op
+  apf::sharedReduction(min_field, 0, false, apf::ReductionMin<double>());
+
+  // compute min/max of avg_size over the whole mesh
+  double mesh_min = 1.e16;
+  double mesh_max = -1.e16;
+  it = m->begin(0);
+  while ( (v = m->iterate(it)) )
+  {
+    double s = apf::getScalar(sum_field, v, 0);
+    if (s > mesh_max)
+      mesh_max = s;
+    if (s < mesh_min)
+      mesh_min = s;
+  }
+  m->end(it);
+
+  PCU_Min_Doubles(&mesh_min, 1);
+  PCU_Max_Doubles(&mesh_max, 1);
+
+  if (!PCU_Comm_Self()) {
+    printf("min/max of current avg size at time step %d: %f/%f\n", ts, mesh_min, mesh_max);
+    printf("user requested max_size at time step %d: %f\n", ts, max_size);
+  }
+
+  int bdim = m->getDimension() - 1;
+
+  double refine_factor = 1.;
+  for (int i = 0; i < refine_level; i++)
+    refine_factor *= 2.;	
+
+  double coarsen_factor = 1.;
+  for (int i = 0; i < coarsen_level; i++)
+    coarsen_factor *= 2.;
+
+  it = m->begin(0);
+  while ( (v = m->iterate(it)) )
+  {
+    double asked_size = apf::getScalar(in_size, v, 0);
+    double avg_size = apf::getScalar(sum_field, v, 0);
+    double min_size = apf::getScalar(min_field, v, 0);
+    int mtype = m->getModelType(m->toModel(v));
+    if (mtype == bdim) {
+      asked_size = avg_size;
+    }
+    else {
+      if (asked_size < min_size / refine_factor) // cap refinement by refine_factor (:=2^refine_level)
+      	asked_size = min_size / refine_factor;
+      if (asked_size > min_size * coarsen_factor) // cap coarsening by coarsen_factor (:=2^coarsen_level)
+      	asked_size = min_size * coarsen_factor;
+    }
+    // cap the biggest size to user specified max_size,
+    // thus never allowing the mesh to get coarser that max_size
+    if (asked_size > max_size)
+      asked_size = max_size;
+    apf::setScalar(in_size, v, 0, asked_size);
+  }
+  m->end(it);
+  apf::synchronize(in_size);
+
+  // compute min/max of avg_size over the whole mesh
+  double asked_min = 1.e16;
+  double asked_max = -1.e16;
+  it = m->begin(0);
+  while ( (v = m->iterate(it)) )
+  {
+    double s = apf::getScalar(in_size, v, 0);
+    if (s > asked_max)
+      asked_max = s;
+    if (s < asked_min)
+      asked_min = s;
+  }
+  m->end(it);
+
+  PCU_Min_Doubles(&asked_min, 1);
+  PCU_Max_Doubles(&asked_max, 1);
+
+  if (!PCU_Comm_Self())
+    printf("min/max of asked size at time step %d: %f/%f\n", ts, asked_min, asked_max);
+
+  // clean up
+  m->removeField(sum_field);
+  m->removeField(min_field);
+  m->removeField(cnt_field);
+  apf::destroyField(sum_field);
+  apf::destroyField(min_field);
+  apf::destroyField(cnt_field);
+}
+
+// for 3D
+// typedefs
+typedef std::vector<apf::Field*> MultiField; // fields on different planes
+
+// the naming convention of multi-plane fields is "name_pnnn" where
+//   "name" is the original name of the field
+//   "nnn" is the plane number, zero-padded to have length equal to zero_pad_length
+const int zero_pad_length = 3; // this would allow number of planes to be between 0 and 999
+const int max_num_plane = 999; // this should be equal (10^zero_pad_length - 1)
+
+static const char* get_plane_name_format(int pad = zero_pad_length)
+{
+  PCU_ALWAYS_ASSERT(pad >= 3 && pad < 7);
+  static const char* format[7] =
+  {
+    "", // 0
+    "", // 1
+    "", // 2
+    "%s_p%03d", // 3
+    "%s_p%04d", // 3
+    "%s_p%05d", // 3
+    "%s_p%06d", // 3
+  };
+  return format[pad];
+}
+
+static void get_original_field_name(const char* plane_name, char* original_name)
+{
+  std::string plane_name_str(plane_name);
+  std::size_t s = plane_name_str.size();
+  s -= (zero_pad_length+2);
+  plane_name_str.resize(s);
+
+  std::strcpy(original_name, plane_name_str.c_str());
+}
+
+// helper functions
+static bool is_in_plane(apf::MeshEntity* e)
+{
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+  PCU_ALWAYS_ASSERT(m->getType(e) == apf::Mesh::EDGE);
+
+  double tol = 2. * M3DC1_PI / m3dc1_model::instance()->num_plane / 1.e6;
+  apf::Vector3 p[2];
+  apf::MeshEntity* v[2];
+  m->getDownward(e, 0, v);
+  for (int i = 0; i < 2; i++) {
+    m->getPoint(v[i], 0, p[i]);
+  }
+  return std::fabs(p[0][2] - p[1][2]) < tol;
+}
+
+static void destroyElement(apf::Mesh2* m, apf::MeshEntity* e)
+{
+  int dim = apf::getDimension(m,e);
+  if (dim < m->getDimension())
+  { //destruction is a no-op if this entity still supports
+    //higher-order ones
+    if (m->hasUp(e))
+      return;
+  }
+  if (dim < m->getDimension())
+  {
+    int etype = m->getType(e);
+    if (etype == apf::Mesh::TRIANGLE) return;
+    if (etype == apf::Mesh::EDGE)
+    {
+      if (is_in_plane(e)) return;
+    }
+    /* if (etype == apf::Mesh::VERTEX) */
+  }
+  apf::Downward down;
+  int nd = 0;
+  if (dim > 0)
+    nd = m->getDownward(e,dim-1,down);
+  m->destroy(e);
+  /* destruction applies recursively to the closure of the entity */
+  if (dim > 0)
+    for (int i=0; i < nd; ++i)
+      destroyElement(m,down[i]);
+}
+
+static void remove_all_wedges()
+{
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+
+  apf::MeshEntity* e;
+  apf::MeshIterator* it;
+
+  it = m->begin(3);
+  while ( (e = m->iterate(it)) )
+  {
+    destroyElement(m, e);
+  }
+  m->acceptChanges();
+  changeMdsDimension(m, 2);
+}
+
+
+
+static void transfer_field_data_on_plane(apf::Field* in, apf::Field* out, int p)
+{
+  int np = m3dc1_model::instance()->num_plane;
+  int local_planeid = m3dc1_model::instance()->local_planeid;
+  PCU_ALWAYS_ASSERT_VERBOSE(p >= 0, "p must be strictly positive!");
+  PCU_ALWAYS_ASSERT_VERBOSE(p < np, "p must be strictly less than number of planes!");
+  int nc = apf::countComponents(in);
+  apf::NewArray<double> dofs(nc);
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+  apf::MeshEntity* e;
+  apf::MeshIterator* it = m->begin(0);
+  while ( (e = m->iterate(it)) )
+  {
+    // if not on the plane continue
+    if (p != local_planeid) continue;
+    // if not owned continue
+    /* if (!m->isOwned(e)) continue; */
+    apf::getComponents(in, e, 0, &dofs[0]);
+    apf::setComponents(out, e, 0, &dofs[0]);
+  }
+  m->end(it);
+  synchronize_field(out);
+}
+
+static void move_field_data_down(apf::Field* f, int startp)
+{
+  int np = m3dc1_model::instance()->num_plane;
+  int local_planeid = m3dc1_model::instance()->local_planeid;
+  PCU_ALWAYS_ASSERT_VERBOSE(startp > 0, "startp must be strictly positive!");
+  PCU_ALWAYS_ASSERT_VERBOSE(startp < np, "startp must be strictly less than number of planes!");
+  int nc = apf::countComponents(f);
+  apf::NewArray<double> zeros(nc);
+  for (int i = 0; i < nc; i++)
+    zeros[i] = 0.;
+
+
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+  apf::MeshEntity* e;
+  apf::MeshIterator* it = m->begin(1);
+  while ( (e = m->iterate(it)) )
+  {
+    // continue if the edge is in any of the poloidal planes
+    if (is_in_plane(e)) continue;
+    // continue if the edge is not on local plane
+    if (startp-1 != local_planeid) continue;
+    apf::MeshEntity* vs[2];
+    m->getDownward(e, 0, vs);
+    apf::Vector3 ps[2];
+    for (int i = 0; i < 2; i++)
+      m->getPoint(vs[i], 0, ps[i]);
+    // make sure the second vert in vs is on lower plane (i.e. has lower z component)
+    if (ps[1][2] > ps[0][2])
+    {
+      std::swap(ps[0], ps[1]);
+      std::swap(vs[0], vs[1]);
+    }
+    apf::NewArray<double> dofs(nc);
+    apf::getComponents(f, vs[0], 0, &dofs[0]);
+    apf::setComponents(f, vs[1], 0, &dofs[0]);
+    apf::setComponents(f, vs[0], 0, &zeros[0]);
+  }
+  m->end(it);
+  synchronize_field(f);
+  /* PCU_Barrier(); */
+}
+
+static void transfer_field_to_main(apf::Field* f, MultiField& pfields)
+{
+  apf::MeshEntity* e;
+  apf::MeshIterator* it;
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+  int np = m3dc1_model::instance()->num_plane;
+  int nc = apf::countComponents(f);
+
+  PCU_ALWAYS_ASSERT(np <= max_num_plane);
+
+  for (int i = 0; i < np; i++) {
+    char pfieldname[128];
+    sprintf(pfieldname, get_plane_name_format(), apf::getName(f), i);
+    apf::Field* tmp = createPackedField(m, pfieldname, nc, apf::getShape(f));
+    apf::zeroField(tmp);
+    pfields.push_back(tmp);
+  }
+
+
+  for (int i = 0; i < (int)pfields.size(); i++)
+    transfer_field_data_on_plane(f, pfields[i], i);
+
+
+  for (int i = 1; i < (int)pfields.size(); i++)
+    for (int j = 0; j < i; j++)
+      move_field_data_down(pfields[i], i-j);
+}
+
+void transfer_field_from_main(MultiField& pfields)
+{
+  apf::MeshEntity* e;
+  apf::MeshIterator* it;
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+  int np = m3dc1_model::instance()->num_plane;
+  PCU_ALWAYS_ASSERT(np == (int)pfields.size());
+  int lpid = m3dc1_model::instance()->local_planeid;
+  int nc = apf::countComponents(pfields[0]);
+
+  // get the original field corresponding to multiplane pfields
+  char original_name[128];
+  get_original_field_name(apf::getName(pfields[0]), original_name);
+  apf::Field* f = m->findField(original_name);
+  PCU_ALWAYS_ASSERT_VERBOSE(f, "was not able to find the original field!");
+
+
+  /* // first transfer everything to the last plane */
+  it = m->begin(1);
+  while ( (e = m->iterate(it)) )
+  {
+    if (lpid != np-1) continue;
+    if (is_in_plane(e)) continue;
+    apf::MeshEntity* v[2];
+    apf::Vector3 p[2];
+    m->getDownward(e, 0, v);
+    for (int i = 0; i < 2; i++)
+      m->getPoint(v[i], 0, p[i]);
+    if (p[0][2] > p[1][2])
+    {
+      std::swap(v[0], v[1]);
+      std::swap(p[0], p[1]);
+    }
+    apf::NewArray<double> dofs(nc);
+    for (int j = 0; j < (int)pfields.size(); j++) {
+      apf::getComponents(pfields[j], v[0], 0, &dofs[0]);
+      apf::setComponents(pfields[j], v[1], 0, &dofs[0]);
+    }
+  }
+  m->end(it);
+
+  for (int i = 0; i < (int)pfields.size(); i++)
+    synchronize_field(pfields[i]);
+
+  for (int i = 1; i < np-1 ; i++)
+    for (int j = np-1; j > i; j--)
+      move_field_data_down(pfields[i], j);
+
+  for (int i = 0; i < (int)pfields.size(); i++)
+    transfer_field_data_on_plane(pfields[i], f, i);
+}
+
+static apf::Field* compute_multiplane_size_field(const MultiField& pfields, double ar)
+{
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+  int np = m3dc1_model::instance()->num_plane;
+  PCU_ALWAYS_ASSERT(np == (int)pfields.size());
+
+  int nc = apf::countComponents(pfields[0]);
+
+  apf::MeshEntity* e;
+  apf::MeshIterator* it;
+
+  MultiField sizefields; // size fields computed for each plane
+
+  for (int i = 0; i < np; i++) {
+    apf::Field* ip = get_ip_field(m, pfields[i]);
+    apf::Field* size = spr::getSPRSizeField(ip, ar);
+    char szname[128];
+    sprintf(szname, "sz_%s", apf::getName(pfields[i]));
+    apf::Field* sz = apf::createField(m, szname, apf::SCALAR, apf::getShape(size));
+    apf::copyData(sz, size);
+    sizefields.push_back(sz);
+    m->removeField(size);
+    apf::destroyField(size);
+    m->removeField(ip);
+    apf::destroyField(ip);
+    /* if (i != np-1) { */
+    /*   m->removeField(ip); */
+    /*   apf::destroyField(ip); */
+    /* } */
+  }
+
+  apf::Field* mastersize = apf::createField(m, "size", apf::SCALAR, apf::getShape(sizefields[0]));
+
+  it = m->begin(0);
+  while ( (e = m->iterate(it)) )
+  {
+    double minsize = 1.e16;
+    for (int i = 0; i < np; i++) {
+      double s = apf::getScalar(sizefields[i], e, 0);
+      if (s < minsize)
+	minsize = s;
+    }
+    apf::setScalar(mastersize, e, 0, minsize);
+  }
+  m->end(it);
+  /* apf::writeVtkFiles("03_mesh_with_all_sizes", m); */
+  for (int i = 0; i < np; i++) {
+    m->removeField(sizefields[i]);
+    apf::destroyField(sizefields[i]);
+  }
+  synchronize_field(mastersize);
+  /* apf::writeVtkFiles("04_mesh_with_master_size", m); */
+  return mastersize;
+}
+
+static void zero_fields_on_non_master(const MultiField& pfields)
+{
+  int np = m3dc1_model::instance()->num_plane;
+  PCU_ALWAYS_ASSERT(np == (int)pfields.size());
+
+  int nc = apf::countComponents(pfields[0]);
+  int local_planeid = m3dc1_model::instance()->local_planeid;
+
+  double tol = 2. * M3DC1_PI / np / 1.e6;
+
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+
+  // get the original field corresponding to multiplane pfields
+  char original_name[128];
+  get_original_field_name(apf::getName(pfields[0]), original_name);
+  apf::Field* f = m->findField(original_name);
+  PCU_ALWAYS_ASSERT_VERBOSE(f, "was not able to find the original field!");
+
+
+  apf::NewArray<double> zeros(nc);
+  for (int i = 0; i < nc; i++)
+    zeros[i] = 0.;
+
+  apf::MeshEntity* e;
+  apf::MeshIterator* it = m->begin(0);
+  while ( (e = m->iterate(it)) )
+  {
+    apf::setComponents(f, e, 0, &zeros[0]);
+
+    if (local_planeid == 0) continue;
+
+    for (int j = 0; j < (int)pfields.size(); j++)
+      apf::setComponents(pfields[j], e, 0, &zeros[0]);
+  }
+  m->end(it);
+
+  it = m->begin(1);
+  while ( (e = m->iterate(it)) )
+  {
+    if (local_planeid != 0) continue;
+    if (is_in_plane(e)) continue;
+    apf::MeshEntity* v[2];
+    apf::Vector3 p[2];
+    m->getDownward(e, 0, v);
+    for (int i = 0; i < 2; i++) {
+      m->getPoint(v[i], 0, p[i]);
+    }
+
+    if (std::fabs(p[0][2]) > tol)
+    {
+      for (int j = 0; j < (int)pfields.size(); j++)
+	apf::setComponents(pfields[j], v[0], 0, &zeros[0]);
+    }
+    else
+    {
+      PCU_ALWAYS_ASSERT(std::fabs(p[1][2]) > tol);
+      for (int j = 0; j < (int)pfields.size(); j++)
+	apf::setComponents(pfields[j], v[1], 0, &zeros[0]);
+    }
+  }
+  m->end(it);
+  for (int j = 0; j < np; j++)
+    apf::synchronize(pfields[j]);
+}
+
+// the following should be called only on master plane 0
+static void zero_fields_on_master(apf::Field* f)
+{
+  int np = m3dc1_model::instance()->num_plane;
+  int nc = apf::countComponents(f);
+  int local_planeid = m3dc1_model::instance()->local_planeid;
+  PCU_ALWAYS_ASSERT_VERBOSE(local_planeid==0, "function can only be called on master palne!");
+
+  double tol = 2. * M3DC1_PI / np / 1.e6;
+
+  apf::Mesh2* m = m3dc1_mesh::instance()->mesh;
+
+  apf::NewArray<double> zeros(nc);
+  for (int i = 0; i < nc; i++)
+    zeros[i] = 0.;
+
+  apf::MeshEntity* e;
+  apf::MeshIterator* it = m->begin(0);
+  while ( (e = m->iterate(it)) )
+    apf::setComponents(f, e, 0, &zeros[0]);
+  m->end(it);
+  for (int j = 0; j < np; j++)
+    apf::synchronize(f);
+}
+// end of static functions used for spr-adapt
+
+int begin_numVert;
+//#endif
+
 double begin_mem, begin_time;
 // helper routines
 void group_complex_dof (apf::Field* field, int option);
@@ -76,6 +734,14 @@ int m3dc1_scorec_init()
   return M3DC1_SUCCESS; 
 }
 
+//*******************************************************
+int m3dc1_scorec_verbosity(int* l)
+//*******************************************************
+{ 
+  lion_set_verbosity(*l);
+  return M3DC1_SUCCESS; 
+}
+
 //*******************************************************
 int m3dc1_scorec_finalize()
 //*******************************************************
@@ -435,6 +1101,8 @@ int m3dc1_mesh_load(char* mesh_file)
         apf::removeTagFromDimension(mesh, tags[i], idim);
       mesh->destroyTag(tags[i]);
     }
+
+    // TODO check and make sure all fields/tags are gone
   } 
   else // non-master plane
     m3dc1_mesh::instance()->mesh = pumi_mesh_create(pumi::instance()->model, 2, false);
@@ -532,6 +1200,321 @@ void m3dc1_dir_import(double* dir, int ts)
   fclose(fp);
 }
 
+
+int m3dc1_spr_then_adapt (FieldID* field_id, int* index, int* ts,
+    double* ar, double* max_size, int* refine_level, int* coarsen_level, bool* update)
+{
+  char filename[256];
+#ifdef DEBUG
+  if (!PCU_Comm_Self())
+    std::cout<<"[M3D-C1 INFO] "<<__func__<<" field id "<<*field_id<<" , name "<<get_field_name_from_id(*field_id)<<"\n";
+#endif
+
+  // delete all the matrix
+#ifdef M3DC1_TRILINOS
+  while (m3dc1_ls::instance()-> matrix_container->size())
+  {
+    std::map<int, m3dc1_epetra*>::iterator mat_it = m3dc1_ls::instance()->matrix_container->begin();
+    mat_it->second->destroy();
+    delete mat_it->second;
+    m3dc1_ls::instance()->matrix_container->erase(mat_it);
+  }
+#endif
+#ifdef M3DC1_PETSC
+  while (m3dc1_solver::instance()-> matrix_container->size())
+  {
+    std::map<int, m3dc1_matrix*> :: iterator mat_it = m3dc1_solver::instance()-> matrix_container->begin();
+    delete mat_it->second;
+    m3dc1_solver::instance()->matrix_container->erase(mat_it);
+  }
+#endif
+
+
+  apf::Mesh2* mesh = m3dc1_mesh::instance()->mesh;
+  int np = m3dc1_model::instance()->num_plane;
+
+  // in_filed will hold all the dofs of all the fields (num being the total number of fields)
+  // at each vertex. e.g.
+  // f1_1, f1_2, f1_3, f1_4, f1_5, f1_6, ! dofs of 1st field
+  // f2_1, f2_2, f2_3, f2_4, f2_5, f2_6, ! dofs of 2nd field
+  // ...
+  //
+  // findex_1, findex_2, findex_3, findex_4, findex_5, findex_6, ! dofs of index'th field
+  // ...
+  // fnum_1, fnum_2, fnum_3, fnum_5, fnum_5, fnum_6 ! dofs of num'th (last) field
+  apf::Field* inField = (*m3dc1_mesh::instance()->field_container)[*field_id]->get_field();
+  PCU_ALWAYS_ASSERT_VERBOSE(inField, "pointer is empty!");
+  if (!PCU_Comm_Self())
+    std::cout << "received field with name " << apf::getName(inField) << "to run spr on" << std::endl;
+
+  // the following call will extract the ones at index
+  apf::Field* targetField = get_field_at_index(mesh, inField, *index, dofNode);
+
+  // transfer targetField (for spr_computations) and all the fields
+  // that need to be transfered for the next solve step  onto the master-plane
+  // the vector pFields is only used in 3D
+  // Preprocessing 3D mesh. This involves the following steps
+  // 1- copy the fields that are needed for solution transfer onto the master plane
+  // 2- convert the mesh to 2D
+
+
+  // the vector pFields and zFields are only used for 3D
+  // pFields holds the multi-plane fields that need to be transfered during adapt
+  // zFields holds the other fields so they can be zero-ed out after adapt+3D mesh reconstruction
+  std::vector<MultiField> pFields;
+  pFields.clear();
+
+  std::vector<apf::Field*> zFields;
+  zFields.clear();
+
+  if (np > 1) // 3D
+  {
+    MultiField targetMultiField;
+    transfer_field_to_main(targetField, targetMultiField);
+    pFields.push_back(targetMultiField);
+    mesh->removeField(targetField);
+    apf::destroyField(targetField);
+
+    std::map<FieldID, m3dc1_field*>::iterator it = m3dc1_mesh::instance()->field_container->begin();
+    while(it!=m3dc1_mesh::instance()->field_container->end())
+    {
+      apf::Field* field = it->second->get_field();
+      int complexType = it->second->get_value_type();
+      /* assert(valueType==complexType); */
+      if (complexType) group_complex_dof(field, 1);
+      if (isFrozen(field)) unfreeze(field);
+      if (it->second->should_transfer())
+      {
+      	MultiField mf;
+      	transfer_field_to_main(field, mf);
+	pFields.push_back(mf);
+      }
+      else
+      	zFields.push_back(field);
+      it++;
+    }
+    // also add the targetMultiField fields to be zero-ed out after adapt
+    for (int i = 0; i < (int)targetMultiField.size() ; i++)
+      zFields.push_back(targetMultiField[i]);
+
+    m3dc1_mesh::instance()->remove3D();
+    m3dc1_mesh::instance()->rebuildPointersOnNonMasterPlane(pFields, zFields);
+    // Important: remvoe3D+rebuildPointersOnNonMasterPlane modifie the mesh pointer
+    // in m3dc1_mesh::instance(). Therefore the local variable pointing to mesh pointer
+    // has to be updated to reflect that.
+    mesh = m3dc1_mesh::instance()->mesh;
+  }
+
+  int valueType = (*(m3dc1_mesh::instance()->field_container))[*field_id]->get_value_type();
+  vector<apf::Field*> fields;
+  std::map<FieldID, m3dc1_field*> :: iterator it=m3dc1_mesh::instance()->field_container->begin();
+
+  if (m3dc1_model::instance()->num_plane == 1) // 2D
+  {
+    while(it!=m3dc1_mesh::instance()->field_container->end())
+    {
+      apf::Field* field = it->second->get_field();
+      int complexType = it->second->get_value_type();
+      assert(valueType==complexType);
+      if (complexType) group_complex_dof(field, 1);
+      if (isFrozen(field)) unfreeze(field);
+      if (it->second->should_transfer())
+      {
+	if (!PCU_Comm_Self()) std::cout<<"[M3D-C1 INFO] "<<__func__<<": field with name "<<getName(field)<< " is added to solution transfer\n";
+	fields.push_back(field);
+      }
+      it++;
+    }
+  }
+  else // 3D
+  {
+    // ignoring the first one since that is only used for size calculations
+    // hence the reason for the following loop start from 1
+    for (int i = 1; i < (int)pFields.size(); i++) {
+      PCU_ALWAYS_ASSERT((int)pFields[i].size() == np);
+      for (int j = 0; j < np; j++)
+      	fields.push_back(pFields[i][j]);
+    }
+  }
+
+  if (!PCU_Comm_Self()) std::cout<<"[M3D-C1 INFO] "<<__func__<<": "<<fields.size()<<" fields have been added to solution transfer\n";
+  for (int i = 0; i < (int)fields.size(); i++) {
+    if (!PCU_Comm_Self()) std::cout<<"[M3D-C1 INFO] "<<__func__<< " field with name " << apf::getName(fields[i]) <<  "fields have been added to solution transfer\n";
+  }
+
+  while(mesh->countNumberings())
+  {
+    apf::Numbering* n = mesh->getNumbering(0);
+    if (!PCU_Comm_Self()) std::cout<<"[M3D-C1 INFO] "<<__func__<<": numbering "<<getName(n)<<" deleted\n";
+    apf::destroyNumbering(n);
+  }
+
+  // compute the size field here and remove the ip and target fields afterwards
+  apf::Field* size_field = 0;
+  ma::Input* in;
+
+  if (np == 1) // 2D
+  {
+    //apf::Field* ip = get_ip_field(mesh, targetField);
+    apf::Field* targetField0 = get_vectorComponent_of_field(mesh, targetField);
+    apf::Field* ip = spr::getGradIPField(targetField0, "ip", 2);
+
+    size_field = spr::getSPRSizeField(ip, *ar);
+    process_size_field(mesh, size_field, *ts, *max_size, *refine_level, *coarsen_level);
+    fields.push_back(size_field);
+
+    mesh->removeField(ip);
+    mesh->removeField(targetField);
+    mesh->removeField(targetField0);
+    destroyField(ip);
+    destroyField(targetField);
+    destroyField(targetField0);
+
+    ReducedQuinticImplicit shape;
+    ReducedQuinticTransfer slnTrans(mesh,fields, &shape);
+    in = ma::makeAdvanced(ma::configure(mesh, size_field, &slnTrans));
+
+    in->shouldSnap=false;
+    in->shouldTransferParametric=false;
+    in->shouldRunPostZoltan = true;
+    in->goodQuality = 0.5;
+    in->maximumIterations = (*refine_level) + 1;
+
+    // turn off coarsening if coarsen_level is negative
+    if (coarsen_level < 0)
+      in->shouldCoarsen=false;
+
+    ma::adapt(in);
+    mesh->removeField(size_field);
+    apf::destroyField(size_field);
+    apf::reorderMdsMesh(mesh);
+
+    m3dc1_mesh::instance()->initialize();
+    // Note: These are only needed for 2D. For 3D these are called
+    // at the end of restore3D
+    compute_globalid(m3dc1_mesh::instance()->mesh, 0);
+    compute_globalid(m3dc1_mesh::instance()->mesh, m3dc1_mesh::instance()->mesh->getDimension());
+  }
+  else // 3D
+  {
+    // change the comm
+    MPI_Comm groupComm;
+    int groupSize = PCU_Comm_Peers()/np;
+    int lpid = PCU_Comm_Self()/groupSize;
+    int grnk = PCU_Comm_Self()%groupSize;
+    MPI_Comm_split(m3dc1_model::instance()->oldComm, lpid, grnk, &groupComm);
+    PCU_Switch_Comm(groupComm);
+    // size filed computation and adapt applied only on the master plane
+    if (m3dc1_model::instance()->local_planeid == 0)
+    {
+      // update the 2D part.smb if "update" is true
+      if (*update)
+      {
+      	int s = 1;
+      	int ns = -1;
+      	m3dc1_mesh_write("part", &s, &ns);
+      }
+      // NOTE: pFields[0] holds the target fields we need to run spr on
+      size_field = compute_multiplane_size_field(pFields[0], *ar);
+      /* return 0; */
+      process_size_field(mesh, size_field, *ts, *max_size, *refine_level, *coarsen_level);
+      fields.push_back(size_field);
+
+      ReducedQuinticImplicit shape;
+      ReducedQuinticTransfer slnTrans(mesh,fields, &shape);
+      // the commented lines are for debugging
+      /* ma::Input* in = ma::makeAdvanced(ma::configureIdentity(mesh, 0, &slnTrans)); */
+      /* ma::Input* in = ma::makeAdvanced(ma::configureUniformRefine(mesh, 1, &slnTrans)); */
+      ma::Input* in = ma::makeAdvanced(ma::configure(mesh, size_field, &slnTrans));
+
+      in->shouldSnap=false;
+      in->shouldFixShape = true;
+      in->shouldTransferParametric=false;
+      in->shouldRunPostZoltan = true;
+      in->goodQuality = 0.5;
+      in->maximumIterations = (*refine_level);
+
+      // turn off coarsening if coarsen_level is negative
+      if (coarsen_level < 0)
+	in->shouldCoarsen=false;
+
+      ma::adapt(in);
+
+      for (int i = 0; i < (int)zFields.size(); i++)
+	zero_fields_on_master(zFields[i]);
+
+      mesh->removeField(size_field);
+      apf::destroyField(size_field);
+
+      // remove numberings
+      while (mesh->countNumberings())
+      {
+	apf::Numbering* n = mesh->getNumbering(0);
+	mesh->removeNumbering(n);
+	apf::destroyNumbering(n);
+      }
+
+      for (int i = 1; i < (int)pFields.size(); i++)
+	for (int j = 0; j < np; j++)
+	  synchronize_field(pFields[i][j]);
+      apf::reorderMdsMesh(mesh);
+
+    }
+
+    // switch comm back to original
+    PCU_Switch_Comm(m3dc1_model::instance()->oldComm);
+    MPI_Comm_free(&groupComm);
+
+    m3dc1_mesh::instance()->restore3D();
+
+    for (int i = 1; i < (int)pFields.size(); i++)
+      zero_fields_on_non_master(pFields[i]);
+
+    for (int i = 1; i < (int)pFields.size(); i++)
+      transfer_field_from_main(pFields[i]);
+
+    // clean up multi-plane fields
+    for (int i = 1; i < (int)pFields.size(); i++)
+      for (int j = 0; j < (int)pFields[i].size(); j++) {
+      	mesh->removeField(pFields[i][j]);
+	apf::destroyField(pFields[i][j]);
+      }
+
+    // zero out all the fields that are not transfered during adapt
+    for (int i = 0; i < (int)zFields.size(); i++)
+      apf::zeroField(zFields[i]);
+
+    // delete pFields[0] fields here
+    for (int i = 0; i < (int)pFields[0].size(); i++) {
+      mesh->removeField(pFields[0][i]);
+      apf::destroyField(pFields[0][i]);
+    }
+  }
+
+  it=m3dc1_mesh::instance()->field_container->begin();
+  while(it!=m3dc1_mesh::instance()->field_container->end())
+  {
+    apf::Field* field = it->second->get_field();
+    int complexType = it->second->get_value_type();
+    if (complexType) group_complex_dof(field, 0);
+    if (!isFrozen(field)) freeze(field);
+#ifdef DEBUG
+    int isnan;
+    int fieldId= it->first;
+    m3dc1_field_isnan(&fieldId, &isnan);
+    assert(isnan==0);
+#endif
+    synchronize_field(field);
+
+#ifdef DEBUG
+    m3dc1_field_isnan(&fieldId, &isnan);
+    assert(isnan==0);
+#endif
+    it++;
+  }
+  return M3DC1_SUCCESS;
+}
+
 /* new mesh adaptation */
 /* Input Parameters
  * field_id_h1, field_id_h2: removed before adaptation so it won't be available after adaptation
@@ -557,6 +1540,7 @@ void m3dc1_mesh_adapt(int* field_id_h1, int* field_id_h2, double* dir)
 #endif
 }
 
+
 /* ghosting functions */
 //*******************************************************
 int m3dc1_ghost_create (int* num_layer )
@@ -948,9 +1932,22 @@ int m3dc1_ent_getgeomclass (int* /* in */ ent_dim, int* /* in */ ent_id,
 { 
   apf::MeshEntity* ent = getMdsEntity(m3dc1_mesh::instance()->mesh, *ent_dim, *ent_id);
   assert(ent);
+  /* apf::Vector3 p; */
+  /* apf::Mesh2* m = m3dc1_mesh::instance()->mesh; */
+  /* if (m->getType(ent) == apf::Mesh::VERTEX) */
+  /* { */
+  /*   m->getPoint(ent, 0, p); */
+  /*   printf("~~~~~ id%d vid%d is owned_%d, coords (%f,%f,%f) mtag ", */
+	/* PCU_Comm_Self(), *ent_id, m->isOwned(ent), p[0], p[1], p[2]); */
+  /*   printf("%d\n", m->getModelTag(m->toModel(ent))); */
+  /* } */
   gmi_ent* gent= (gmi_ent*)(m3dc1_mesh::instance()->mesh->toModel(ent));
   *geom_class_dim = gmi_dim(m3dc1_model::instance()->model,gent);
   *geom_class_id = gmi_tag(m3dc1_model::instance()->model,gent);
+
+  /* if (!PCU_Comm_Self()) */
+  /*   std::cout << "~~~~~ in before 3d " << __func__ << " " */
+  /*             << "geomclassdim/id " << *geom_class_dim << "/" << *geom_class_id << std::endl; */
   // if 3D mesh, need to return the classification on the original plane
   if ( m3dc1_mesh::instance()->mesh->getDimension() ==3 )
   {
@@ -967,6 +1964,9 @@ int m3dc1_ent_getgeomclass (int* /* in */ ent_dim, int* /* in */ ent_id,
     }
     *geom_class_id+=1;
   }
+  /* if (!PCU_Comm_Self()) */
+  /*   std::cout << "~~~~~ in after 3d " << __func__ << " " */
+  /*             << "geomclassdim/id " << *geom_class_dim << "/" << *geom_class_id << std::endl; */
   return M3DC1_SUCCESS;
 }
 
@@ -1512,6 +2512,26 @@ int* /*in*/ scalar_type, int* /*in*/ num_dofs_per_value)
   return M3DC1_SUCCESS;
 }
 
+//*******************************************************
+int m3dc1_mark_for_solutiontransfer (FieldID* /*in*/ field_id)
+//*******************************************************
+{
+  if (!m3dc1_mesh::instance()->field_container)
+    return M3DC1_FAILURE;
+  if (!m3dc1_mesh::instance()->field_container->count(*field_id))
+    return M3DC1_FAILURE;
+
+  apf::Field* f = (*m3dc1_mesh::instance()->field_container)[*field_id]->get_field();
+#ifdef DEBUG
+  if (!PCU_Comm_Self()) 
+    std::cout<<"[M3D-C1 INFO] "<<__func__<<": field "<<*field_id<<", name "<<getName(f)<<" is being marked for solution transfer\n";
+#endif
+
+  m3dc1_field* mf = (*m3dc1_mesh::instance()->field_container)[*field_id];
+  mf->mark_for_solutiontransfer();
+  return M3DC1_SUCCESS;
+}
+
 //*******************************************************
 int m3dc1_field_delete (FieldID* /*in*/ field_id)
 //*******************************************************
@@ -2717,7 +3737,7 @@ int m3dc1_matrix_create(int* matrix_id, int* matrix_type, int* scalar_type, Fiel
 
 #ifdef DEBUG
   if (!PCU_Comm_Self())
-    std::cout<<"[M3D-C1 INFO] "<<__func__<<": matrix "<<*matrix_id<<", field "<<*field_id<<"\n";
+    std::cout<<"[M3D-C1 INFO] "<<__func__<<": matrix "<<*matrix_id<<", field id "<<*field_id<<" , name "<<get_field_name_from_id(*field_id)<<"\n";
 #endif 
 
   if (*matrix_type==M3DC1_MULTIPLY) // matrix for multiplication
@@ -2951,9 +3971,10 @@ int m3dc1_matrix_setlaplacebc(int * matrix_id, int *row,
 int m3dc1_matrix_solve(int* matrix_id, FieldID* rhs_sol) //solveSysEqu_
 {  
   m3dc1_matrix* mat = m3dc1_solver::instance()->get_matrix(*matrix_id);
-#ifdef DEBUG
+/* #ifdef DEBUG */
+  const char* name = apf::getName((*m3dc1_mesh::instance()->field_container)[*rhs_sol]->get_field());
   if (!PCU_Comm_Self())
-     std::cout <<"[M3D-C1 INFO] "<<__func__<<": matrix "<<* matrix_id<<", field "<<*rhs_sol<<"\n";
+     std::cout <<"[M3D-C1 INFO] "<<__func__<<": matrix "<<* matrix_id<<", field "<<*rhs_sol<<", name "<<name<<"\n";
   if (!mat) 
   {  
     if (!PCU_Comm_Self())
@@ -2967,7 +3988,7 @@ int m3dc1_matrix_solve(int* matrix_id, FieldID* rhs_sol) //solveSysEqu_
       std::cout <<"[M3D-C1 ERROR] "<<__func__<<" not supported with matrix for multiplication (id"<<*matrix_id<<")\n";
     return M3DC1_FAILURE;
   }
-#endif
+/* #endif */
 
   (dynamic_cast<matrix_solve*>(mat))->solve(*rhs_sol);
   addMatHit(*matrix_id);
@@ -2982,7 +4003,7 @@ int m3dc1_matrix_multiply(int* matrix_id, FieldID* inputvecid,
   m3dc1_matrix* mat = m3dc1_solver::instance()->get_matrix(*matrix_id);
 #ifdef DEBUG
   if (!PCU_Comm_Self())
-     std::cout <<"[M3D-C1 INFO] "<<__func__<<": matrix "<<* matrix_id<<", in-field "<<*inputvecid<<", out-field "<<*outputvecid<<"\n";
+     std::cout <<"[M3D-C1 INFO] "<<__func__<<": matrix "<<* matrix_id<<", in-field "<<*inputvecid<<", in_name "<<get_field_name_from_id(*inputvecid)<<", out-field "<<*outputvecid<<", out_name "<<get_field_name_from_id(*outputvecid)<<"\n";
   if (!mat)
   {
     if (!PCU_Comm_Self())
@@ -3351,7 +4372,11 @@ int adapt_by_field (int * fieldId, double* psi0, double * psil)
     assert(valueType==complexType);
     if (complexType) group_complex_dof(field, 1);
     if (isFrozen(field)) unfreeze(field);
-    fields.push_back(field);
+    if (it->second->should_transfer())
+    {
+      if (!PCU_Comm_Self()) std::cout<<"[M3D-C1 INFO] "<<__func__<<": field with name "<<getName(field)<< " is added to solution transfer\n";
+      fields.push_back(field);
+    }
     it++;
   }
   while(mesh->countNumberings())
@@ -3361,11 +4386,7 @@ int adapt_by_field (int * fieldId, double* psi0, double * psil)
     apf::destroyNumbering(n);
   }
   ReducedQuinticTransfer slnTrans(mesh,fields, &shape);
-#ifdef OLDMA
-  ma::Input* in = ma::configure(mesh,&sf,&slnTrans);
-#else
   ma::Input* in = ma::makeAdvanced(ma::configure(mesh,&sf,&slnTrans));
-#endif
   in->maximumIterations = 9;
 
   in->shouldSnap=false;
@@ -3609,8 +4630,11 @@ int adapt_by_error_field (double * errorData, double * errorAimed, int * max_ada
     int complexType = it->second->get_value_type();
     if (complexType) group_complex_dof(field, 1);
     if (isFrozen(field)) unfreeze(field);
-    //if (!PCU_Comm_Self()) std::cout<<"Solution transfer: add field "<<apf::getName(field)<<std::endl;
-    fields.push_back(field);
+    if (it->second->should_transfer())
+    {
+      if (!PCU_Comm_Self()) std::cout<<"[M3D-C1 INFO] "<<__func__<<": field with name "<<getName(field)<< " is added to solution transfer\n";
+      fields.push_back(field);
+    }
     it++;
   }
   while(mesh->countNumberings())
@@ -3624,11 +4648,7 @@ int adapt_by_error_field (double * errorData, double * errorAimed, int * max_ada
   //apf::writeVtkFiles(filename,mesh);
 
   ReducedQuinticTransfer slnTrans(mesh,fields, &shape);
-#ifdef OLDMA
-  ma::Input* in = ma::configure(mesh,&sf,&slnTrans);
-#else
   ma::Input* in = ma::makeAdvanced(ma::configure(mesh,&sf,&slnTrans));
-#endif
   in->maximumIterations = 5;
   in->shouldSnap=false;
   in->shouldTransferParametric=false;
@@ -3717,7 +4737,10 @@ int m3dc1_mesh_write(char* filename, int *option, int* timestep)
   // vtk
   if (*option==0 ||*option==3)
   {
-    sprintf(filename_buff, "ts%d-%s",*timestep,filename);
+    if (*timestep >= 0)
+      sprintf(filename_buff, "ts%04d-%s",*timestep,filename);
+    else
+      sprintf(filename_buff, "%s",filename);
 
     apf::Mesh2* mesh = m3dc1_mesh::instance()->mesh;
     apf::MeshEntity* e;
@@ -3744,7 +4767,10 @@ int m3dc1_mesh_write(char* filename, int *option, int* timestep)
   }
   else // smb
   {
-    sprintf(filename_buff, "ts%d-%s.smb",*timestep,filename);
+    if (*timestep >= 0)
+      sprintf(filename_buff, "ts%d-%s.smb",*timestep,filename);
+    else
+      sprintf(filename_buff, "%s.smb",filename);
 
     int fieldID=12;
     double dofBuff[1024];
@@ -3762,7 +4788,22 @@ int m3dc1_mesh_write(char* filename, int *option, int* timestep)
       mesh->setDoubleTag(e,tag, dofBuff);
     }
     mesh->end(it);
+
+    std::vector<apf::Field*> allFields;
+    allFields.clear();
+    while (mesh->countFields())
+    {
+      apf::Field* f = mesh->getField(0);
+      apf::freeze(f);
+      allFields.push_back(f);
+      mesh->removeField(f);
+    }
     m3dc1_mesh::instance()->mesh->writeNative(filename_buff);
+    for (std::size_t i = 0; i < allFields.size(); i++) {
+      apf::unfreeze(allFields[i]);
+      mesh->addField(allFields[i]);
+    }
+
     apf::removeTagFromDimension(mesh, tag, dim);
     mesh->destroyTag(tag);
     if (!PCU_Comm_Self()) 
diff --git a/m3dc1_scorec/api/m3dc1_scorec.h b/m3dc1_scorec/api/m3dc1_scorec.h
index 36efd119a..677a21d03 100644
--- a/m3dc1_scorec/api/m3dc1_scorec.h
+++ b/m3dc1_scorec/api/m3dc1_scorec.h
@@ -45,6 +45,7 @@ enum m3dc1_matrix_status { /*0*/ M3DC1_NOT_FIXED=0,
 bool m3dc1_double_isequal(double A, double B);
 
 int m3dc1_scorec_init();
+int m3dc1_scorec_verbosity(int*);
 int m3dc1_scorec_finalize();
 
 /** plane functions */
@@ -133,7 +134,8 @@ int m3dc1_field_getnewid (FieldID* /*out*/field_id);
 // is num_dofs input or output?
 // *value_type is either M3DC1_REAL or M3DC1_COMPLEX
 int m3dc1_field_create (FieldID* /*in*/ field_id, const char* /* in */ field_name, int* num_values, int* value_type, int* num_dofs_per_value);
-int m3dc1_field_delete (FieldID* /*in*/ field_id); 
+int m3dc1_mark_for_solutiontransfer (FieldID* /*in*/ field_id);
+int m3dc1_field_delete (FieldID* /*in*/ field_id);
 
 int m3dc1_field_getinfo(FieldID* /*in*/ field_id, char* /* out*/ field_name, int* num_values, int* value_type, int* total_num_dof);
 
@@ -215,6 +217,8 @@ int m3dc1_matrix_print(int* matrix_id);
 #endif // #ifdef M3DC1_PETSC
 
 // adaptation
+int m3dc1_spr_then_adapt (int * fieldId, int * index, int * ts,
+    double * ar, double * max_size, int * refine_level, int * coarsen_level, bool* update);
 int adapt_by_field (int * fieldId, double* psi0, double * psil);
 int set_adapt_p (double * pp);
 int adapt_by_error_field (double * errorField, double * errorAimed, int* max_node, int* option); // option 0: local error control; 1 global
diff --git a/m3dc1_scorec/api/name_convert.h b/m3dc1_scorec/api/name_convert.h
index 4ef061907..ccc97f2c9 100644
--- a/m3dc1_scorec/api/name_convert.h
+++ b/m3dc1_scorec/api/name_convert.h
@@ -1,4 +1,5 @@
 #define m3dc1_scorec_init m3dc1_domain_init_
+#define m3dc1_scorec_verbosity m3dc1_domain_verbosity_
 #define m3dc1_scorec_finalize m3dc1_domain_finalize_
 #define m3dc1_plane_setnum m3dc1_plane_setnum_
 #define m3dc1_plane_getnum m3dc1_plane_getnum_
@@ -16,6 +17,7 @@
 #define m3dc1_mesh_load m3dc1_mesh_load_
 #define m3dc1_mesh_load_3d m3dc1_mesh_load_3d_
 #define m3dc1_mesh_build3d m3dc1_mesh_build3d_
+#define m3dc1_spr_then_adapt m3dc1_spr_then_adapt_
 #define m3dc1_mesh_adapt m3dc1_mesh_adapt_
 #define m3dc1_mesh_write m3dc1_mesh_write_
 #define m3dc1_mesh_getnument m3dc1_mesh_getnument_
@@ -59,6 +61,7 @@
 #define m3dc1_region_getoriginalface m3dc1_region_getoriginalface_
 #define m3dc1_field_getnewid m3dc1_field_genid_
 #define m3dc1_field_create m3dc1_field_create_
+#define m3dc1_mark_for_solutiontransfer m3dc1_mark_for_solutiontransfer_
 #define m3dc1_field_delete m3dc1_field_delete_
 #define m3dc1_field_exist m3dc1_field_exist_
 #define m3dc1_field_sync m3dc1_field_sync_
diff --git a/m3dc1_scorec/include/m3dc1_field.h b/m3dc1_scorec/include/m3dc1_field.h
index 21d140fb1..358773edb 100644
--- a/m3dc1_scorec/include/m3dc1_field.h
+++ b/m3dc1_scorec/include/m3dc1_field.h
@@ -17,19 +17,23 @@ void group_complex_dof (apf::Field* field, int option);
 class m3dc1_field
 {
 public:
-  m3dc1_field (int i, apf::Field* f, int n, int t, int ndof): id(i), field(f), num_value(n), value_type(t),dof_per_value(ndof) {}
+  m3dc1_field (int i, apf::Field* f, int n, int t, int ndof): id(i), field(f), num_value(n), value_type(t),dof_per_value(ndof) {transfer = false;}
   ~m3dc1_field() {}
   apf::Field* get_field() { return field; }
+  void set_field(apf::Field* f) { field = f; }
   int get_id() { return id; }
   int get_num_value() { return num_value; }
   int get_value_type() { return value_type; }
   int get_dof_per_value() {return dof_per_value;}
+  bool should_transfer() {return transfer;}
+  void mark_for_solutiontransfer() {transfer = true;}
 private:
   int id;
   apf::Field* field; // name and #dofs are available from apf::Field
   int num_value;
   int value_type;
   int dof_per_value;
+  bool transfer;
 };
 
 void synchronize_field(apf::Field* f);
diff --git a/m3dc1_scorec/include/m3dc1_mesh.h b/m3dc1_scorec/include/m3dc1_mesh.h
index 4dd9eb465..7e0b9d292 100644
--- a/m3dc1_scorec/include/m3dc1_mesh.h
+++ b/m3dc1_scorec/include/m3dc1_mesh.h
@@ -58,7 +58,9 @@ class m3dc1_mesh
   void restore3D();
 
   void build3d(int num_field, int* field_id, int* num_dofs_per_value);
-  void initialize(); 
+  void initialize();
+  void rebuildPointersOnNonMasterPlane(std::vector<std::vector<apf::Field*>>& pFields,
+              std::vector<apf::Field*>& zFields);
   void set_mcount(); // fill in # local, own, global mesh entity count
   void update_partbdry(apf::MeshEntity** remote_vertices, apf::MeshEntity** remote_edges, 
               apf::MeshEntity** remote_faces, std::vector<apf::MeshEntity*>& btw_plane_edges, 
diff --git a/m3dc1_scorec/src/m3dc1_adapt.cc b/m3dc1_scorec/src/m3dc1_adapt.cc
index 65cbcd357..d669dffcf 100644
--- a/m3dc1_scorec/src/m3dc1_adapt.cc
+++ b/m3dc1_scorec/src/m3dc1_adapt.cc
@@ -14,7 +14,9 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <PCU.h>
+#include <pcu_util.h>
 #include "apfMDS.h"
+#include "apfField.h"
 #include "ReducedQuinticImplicit.h"
 #include "m3dc1_slnTransfer.h"
 #include "apfShape.h" // getLagrange
@@ -115,7 +117,7 @@ void m3dc1_mesh::remove3D()
         mesh->destroy(e);
       mesh->end(ent_it);
     }
-  
+
   mesh->acceptChanges();
   changeMdsDimension(mesh, 2);
 
@@ -125,6 +127,187 @@ void m3dc1_mesh::remove3D()
     std::cout<<"\n*** Wedges and non-master 2D planes removed ***\n";
 }
 
+// this static function only used in m3dc1_mesh::rebuildPointersOnNonMasterPlane
+static int get_id_in_container(
+    std::map<FieldID, m3dc1_field*>& fcontainer,
+    apf::Field* f)
+{
+  typedef std::map<FieldID, m3dc1_field*> fct;
+
+  int id = -1;
+  for (fct::iterator it = fcontainer.begin(); it != fcontainer.end(); ++it)
+  {
+    m3dc1_field* mf = it->second;
+    if (f == mf->get_field())
+    {
+      id = it->first;
+      break;
+    }
+  }
+  return id;
+}
+// this will rebuild the mesh data-structure on non-master plane after remove3D
+//
+// Notes
+// (A) this is done to ensure "restore3D" behaves the same as "build3d".
+// More specifically, the order of entity creation and iteration remains the same
+// on all non-master planes.
+//
+// (B) the field pointers passed in the arrays will be updated for non-master planes
+void m3dc1_mesh::rebuildPointersOnNonMasterPlane(
+    std::vector<std::vector<apf::Field*>>& pFields, // multi-plane fields
+    std::vector<apf::Field*>& zFields)             // fields that need to be zero-d
+{
+  if (!(m3dc1_model::instance()->num_plane)) // if 2D do nothing
+    return;
+
+  // store the names so they can be updated later on
+  std::vector<std::vector<std::string>> pFieldNames; pFieldNames.clear();
+  std::vector<std::string>              zFieldNames; zFieldNames.clear();
+  std::vector<std::string>              tempNames;
+
+  for (int i = 0; i < (int)pFields.size(); i++)
+  {
+    tempNames.clear();
+    for (int j = 0; j < pFields[i].size(); j++)
+    {
+      tempNames.push_back(std::string(apf::getName(pFields[i][j])));
+    }
+    PCU_ALWAYS_ASSERT(pFields[i].size() == tempNames.size());
+    pFieldNames.push_back(tempNames);
+  }
+  PCU_ALWAYS_ASSERT(pFieldNames.size() == pFields.size());
+
+  for (int i = 0; i < (int)zFields.size(); i++)
+    zFieldNames.push_back(std::string(apf::getName(zFields[i])));
+  PCU_ALWAYS_ASSERT(zFieldNames.size() == zFields.size());
+
+
+  apf::Mesh2* mesh = m3dc1_mesh::instance()->mesh;
+  std::vector<int> fncs; fncs.clear();
+  std::vector<int> ids; ids.clear();  // if the field is in the m3dc1_mesh::instance()->field_container hold the id here
+  std::vector<std::string> fnames; fnames.clear();
+  std::vector<apf::FieldShape*> fshapes; fshapes.clear();
+
+  for (int i = 0; i < mesh->countFields(); i++) {
+    apf::Field* f = mesh->getField(i);
+    fncs.push_back(f->countComponents());
+    fnames.push_back(std::string(apf::getName(f)));
+    fshapes.push_back(apf::getShape(f));
+    int id = get_id_in_container(*m3dc1_mesh::instance()->field_container, f);
+    ids.push_back(id);
+  }
+
+  PCU_Barrier();
+
+  // now remove everything on non-master planes, and recreate them
+  // Removal Phase
+  // =============
+  // a-fields
+  // b-numberings
+  // c-tags
+  // d-internal mesh data-structure by calling destroyNative
+  // =============
+  // Recreate Phase
+  // a-create empty mds meshes on non-master planes
+  // b-set local_entid_tag, own_partid_tag, num_global_adj_node_tag, num_own_adj_node_tag to NULL
+  if (m3dc1_model::instance()->local_planeid != 0)
+  {
+    // manually delete all the fields/numberings and associated tags
+    while ( mesh->countFields() )
+    {
+      apf::Field* f = mesh->getField(0);
+      mesh->removeField(f);
+      apf::destroyField(f);
+    }
+
+    while ( mesh->countNumberings() )
+    {
+      apf::Numbering* n = mesh->getNumbering(0);
+      mesh->removeNumbering(n);
+      apf::destroyNumbering(n);
+    }
+
+    apf::DynamicArray<apf::MeshTag*> tags;
+    mesh->getTags(tags);
+    for (int i=0; i<tags.getSize(); i++)
+    {
+      if (mesh->findTag("norm_curv")==tags[i]) continue;
+      for (int idim=0; idim<4; idim++)
+        apf::removeTagFromDimension(mesh, tags[i], idim);
+      mesh->destroyTag(tags[i]);
+    }
+
+    // destroy the native mesh
+    apf::disownMdsModel(mesh);
+    mesh->destroyNative();
+    apf::destroyMesh(mesh);
+    m3dc1_mesh::instance()->mesh = pumi_mesh_create(pumi::instance()->model, 2, false);
+    mesh = m3dc1_mesh::instance()->mesh;
+
+    local_entid_tag = NULL;
+    own_partid_tag = NULL;
+    num_global_adj_node_tag = NULL;
+    num_own_adj_node_tag = NULL;
+  }
+
+  PCU_Barrier();
+
+
+  // update the field pointers in m3dc1_mesh::field_container
+  if (m3dc1_model::instance()->local_planeid != 0)
+  {
+    for (int i = 0; i < (int)fncs.size(); i++)
+      {
+	apf::Field* newf = apf::createPackedField(mesh, fnames[i].c_str(), fncs[i], fshapes[i]);
+	apf::zeroField(newf);
+	if (ids[i] > -1)
+	{
+	  m3dc1_field* mf = (*m3dc1_mesh::instance()->field_container)[ids[i]];
+	  mf->set_field(newf);
+	}
+      }
+  }
+
+  // update the field pointers in pFields
+  if (m3dc1_model::instance()->local_planeid != 0)
+  {
+    for (int i = 0; i < (int)pFieldNames.size(); i++)
+    {
+      for (int j = 0; j < pFieldNames[i].size(); j++)
+      {
+	apf::Field* f = mesh->findField(pFieldNames[i][j].c_str());
+	PCU_ALWAYS_ASSERT(f);
+	pFields[i][j] = f;
+      }
+    }
+  }
+
+
+  // update the field pointers in zFields
+  if (m3dc1_model::instance()->local_planeid != 0)
+  {
+    for (int i = 0; i < (int)zFieldNames.size(); i++)
+    {
+      apf::Field* f = mesh->findField(zFieldNames[i].c_str());
+      PCU_ALWAYS_ASSERT(f);
+      zFields[i] = f;
+    }
+  }
+
+
+  // since we have removed the Linear numbering for non-master planes, we add it beck here
+  apf::Numbering* linnumbering = mesh->findNumbering("Linear");
+  if (!linnumbering)
+    apf::createNumbering(mesh, "Linear", mesh->getShape(), 1);
+
+  mesh->acceptChanges();
+  set_mcount();
+
+  if (!PCU_Comm_Self())
+    std::cout<<"\n*** Data structures have been re-initialized on non-master planes ***\n";
+}
+
 void compute_size_and_frame_fields(apf::Mesh2* m, double* size_1, double* size_2, 
      double* angle, apf::Field* sizefield, apf::Field* framefield)
 {
@@ -356,12 +539,14 @@ void adapt_mesh (int field_id_h1, int field_id_h2, double* dir)
 	
   ReducedQuinticImplicit shape;
   ReducedQuinticTransfer slnTransfer(mesh,fields, &shape);
+  ma::Input* in = ma::makeAdvanced(ma::configure(mesh, size_field, frame_field, &slnTransfer));
+/*
 #ifdef OLDMA
   ma::Input* in = ma::configure(mesh, size_field, frame_field, &slnTransfer);
 #else
   ma::Input* in = ma::makeAdvanced(ma::configure(mesh, size_field, frame_field, &slnTransfer));
 #endif
-	
+*/	
   in->shouldSnap = 0;
   in->shouldTransferParametric = 0;
   in->shouldRunMidZoltan = 1;
@@ -1222,24 +1407,24 @@ void create_localid(apf::Mesh2* mesh, int dim)
 void m3dc1_mesh::restore3D()
 // *********************************************************
 {
-    // compute the fields to copy from master plane
+    /* // compute the fields to copy from master plane */
     int num_field = 0;
     int* field_id = NULL;
     int* num_dofs_per_value = NULL;
 
-    if (!m3dc1_model::instance()->local_planeid && m3dc1_mesh::instance()->field_container)
-    {
-      num_field = m3dc1_mesh::instance()->field_container->size();
-      cout<<__func__<<": #fields to copy from master to non-master plane "<<num_field<<"\n";
-      field_id = new int [num_field];
-      num_dofs_per_value = new int [num_field];
-      for (int i=0; i<num_field; ++i)
-      {
-        m3dc1_field * mf = (*(m3dc1_mesh::instance()->field_container))[*field_id];
-        field_id[i] = i;
-        num_dofs_per_value[i] = mf->get_dof_per_value();
-      }
-    }
+    /* if (!m3dc1_model::instance()->local_planeid && m3dc1_mesh::instance()->field_container) */
+    /* { */
+    /*   num_field = m3dc1_mesh::instance()->field_container->size(); */
+    /*   cout<<__func__<<": #fields to copy from master to non-master plane "<<num_field<<"\n"; */
+    /*   field_id = new int [num_field]; */
+    /*   num_dofs_per_value = new int [num_field]; */
+    /*   for (int i=0; i<num_field; ++i) */
+    /*   { */
+    /*     m3dc1_field * mf = (*(m3dc1_mesh::instance()->field_container))[*field_id]; */
+    /*     field_id[i] = i; */
+    /*     num_dofs_per_value[i] = mf->get_dof_per_value(); */
+    /*   } */
+    /* } */
 
   int local_partid=PCU_Comm_Self();
 
@@ -1487,9 +1672,9 @@ void m3dc1_mesh::restore3D()
   apf::Numbering* local_n = mesh->findNumbering(mesh->getShape()->getName());
   if (local_n) destroyNumbering(local_n);
   
-  // FIXME: re-create the field and copy field data on master process group to non-master
-  for (int i=0; i<num_field; ++i)
-    update_field(field_id[i], num_dofs_per_value[i], num_local_vtx, remote_vertices);
+  /* // FIXME: re-create the field and copy field data on master process group to non-master */
+  /* for (int i=0; i<num_field; ++i) */
+  /*   update_field(field_id[i], num_dofs_per_value[i], num_local_vtx, remote_vertices); */
 
   // clear temp memory
   delete [] remote_vertices;
diff --git a/m3dc1_scorec/src/m3dc1_slnTransfer.cc b/m3dc1_scorec/src/m3dc1_slnTransfer.cc
index 5c04400d0..ac912a5be 100644
--- a/m3dc1_scorec/src/m3dc1_slnTransfer.cc
+++ b/m3dc1_scorec/src/m3dc1_slnTransfer.cc
@@ -7,16 +7,23 @@
   BSD license as described in the LICENSE file in the top-level directory.
  
 *******************************************************************************/
+#include <assert.h>
+
+// headers from SCOREC/Core
+#include <PCU.h>
+#include <pcu_util.h>
+#include <apfMesh.h>
+#include <apfMesh2.h>
+#include <apf.h>
+#include <apfMDS.h>
+#include <apfNumbering.h>
+
+// local headers
 #include "m3dc1_slnTransfer.h"
 #include "m3dc1_scorec.h"
 #include "m3dc1_mesh.h"
 #include "m3dc1_field.h"
-#include "apfMesh.h"
-#include "apfMesh2.h"
-#include "apf.h"
-#include "apfMDS.h"
-#include "apfNumbering.h"
-#include <assert.h>
+
 
 int ReducedQuinticTransfer::dofNode = C1TRIDOFNODE;
 void  ReducedQuinticTransfer::onVertex(apf::MeshElement* parent, ma::Vector const& xi, ma::Entity* vert)
@@ -37,9 +44,9 @@ void  ReducedQuinticTransfer::onVertex(apf::MeshElement* parent, ma::Vector cons
   assert(num_face<=2);
 #endif
 
-  for(int i=0; i<fields.size(); i++)
+  for(int fi=0; fi<fields.size(); fi++)
   {
-    apf::Field* field = fields.at(i);
+    apf::Field* field = fields.at(fi);
     int numComp = apf::countComponents(field);
     vector<vector<double> > dofsVertex(2);
     apf::MeshEntity*  vertices[2];
@@ -47,12 +54,13 @@ void  ReducedQuinticTransfer::onVertex(apf::MeshElement* parent, ma::Vector cons
     apf::Vector3 xyz;
     m3dc1_mesh::instance()->mesh->getPoint(vert, 0, xyz);
     apf::Vector3 xyz2[2];
+    for (int i = 0; i < 2; i++) {
+      m3dc1_mesh::instance()->mesh->getPoint(vertices[i], 0, xyz2[i]);
+    }
     for( int i=0; i<2; i++)
     {
       dofsVertex.at(i).resize(numComp);
-      apf::Element* vertex = apf::createElement(field,vertices[i]);
-      apf::getComponents(vertex,xi,&(dofsVertex[i][0]));
-      m3dc1_mesh::instance()->mesh->getPoint(vertices[i], 0, xyz2[i]);
+      apf::getComponents(field, vertices[i], 0, &(dofsVertex[i][0]));
     }
     double len1= sqrt((xyz2[0][0]-xyz2[1][0])*(xyz2[0][0]-xyz2[1][0])+(xyz2[0][1]-xyz2[1][1])*(xyz2[0][1]-xyz2[1][1]));
     double len2= sqrt((xyz2[0][0]-xyz[0])*(xyz2[0][0]-xyz[0])+(xyz2[0][1]-xyz[1])*(xyz2[0][1]-xyz[1]));
@@ -93,8 +101,7 @@ void  ReducedQuinticTransfer::onVertex(apf::MeshElement* parent, ma::Vector cons
           miss_flag=1;
           break;
         }
-        apf::Element* vertex = apf::createElement(field,vertices[i]);
-        apf::getComponents(vertex,xi,&(value[numComp*i]));
+        apf::getComponents(field, vertices[i], 0, &(value[numComp*i]));
       }
       assert(!miss_flag);
       if(miss_flag)
@@ -148,8 +155,8 @@ void  ReducedQuinticTransfer::onVertex(apf::MeshElement* parent, ma::Vector cons
         if(hasEntity(field,vertices[i]))
         {
           getComponents(field, vertices[i], 0, &dofsBuff[0]);
-          for(int i=0; i<numComp; i++)
-            newdofs[i]+=dofsBuff[i];
+          for(int j=0; j<numComp; j++)
+            newdofs[j]+=dofsBuff[j];
           numVtx++;
         }
       }
diff --git a/unstructured/adapt.f90 b/unstructured/adapt.f90
index 6d372a4b6..24c043ec9 100644
--- a/unstructured/adapt.f90
+++ b/unstructured/adapt.f90
@@ -2,6 +2,12 @@ module adapt
   use vector_mod
   use scorec_adapt
   implicit none
+  integer :: ispradapt
+  integer :: isprrefinelevel
+  integer :: isprcoarsenlevel
+  integer :: isprntime
+  real :: isprmaxsize
+  real :: isprweight
   real :: adapt_ke
   integer :: iadapt_ntime
   real :: adapt_target_error
@@ -64,6 +70,7 @@ subroutine adapt_by_psi
     real :: psib
 
     call create_field(temporary_field)
+    call mark_field_for_solutiontransfer(temporary_field)
     temporary_field = 0.
 
     if(adapt_pellet_delta.gt.0) then
@@ -287,6 +294,98 @@ subroutine adapt_by_psi
   end subroutine adapt_by_psi
 
 
+! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+  subroutine adapt_by_spr(fid,idx,t,ar,maxsize,refinelevel,coarsenlevel,update)
+! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+    use diagnostics
+    use basic
+    use error_estimate
+    use scorec_mesh_mod
+    use basic
+    use mesh_mod
+    use arrays
+    use newvar_mod
+    use sparse
+    use time_step
+    use auxiliary_fields
+    use scorec_mesh_mod
+    use transport_coefficients
+
+    character(len=32) :: mesh_file_name
+    integer, intent(in) :: fid
+    integer, intent(in) :: idx
+    real, intent(in) :: ar
+    real, intent(in) :: maxsize
+    integer, intent(in) :: t
+    integer, intent(in) :: refinelevel
+    integer, intent(in) :: coarsenlevel
+    logical, intent(in) :: update
+
+
+  !  if (myrank .eq. 0) print *, " error exceeds tolerance, start adapting mesh"
+    call straighten_fields()
+    call m3dc1_spr_then_adapt(fid,idx,t,ar,maxsize,refinelevel,coarsenlevel,update)
+    call space(0)
+    call update_nodes_owned()
+    call reset_itris()
+    call tridef
+    call unstraighten_fields()
+
+    call create_newvar_matrices
+    if(irestart .ne. 0 .or. linear .eq. 0) return
+    if (myrank .eq. 0) print *, "reset simulation after adapt .."
+    field_vec = 0.
+    field0_vec = 0.
+    jphi_field = 0.
+    !vor_field = 0.
+    !com_field = 0.
+    !resistivity_field = 0.
+    !kappa_field = 0.
+    !visc_field = 0.
+    !visc_c_field = 0.
+    if(ipforce.gt.0) pforce_field = 0.
+    if(ipforce.gt.0) pmach_field = 0.
+    if(momentum_source) Fphi_field = 0.
+    if(heat_source) Q_field = 0.
+    if(rad_source) then
+      Totrad_field = 0.
+      Linerad_field = 0.
+      Bremrad_field = 0.
+      Ionrad_field = 0.
+      Reckrad_field = 0.
+      Recprad_field = 0.
+    endif
+    if(icd_source.gt.0) cd_field = 0.
+    bf_field(0) = 0.
+    bf_field(1) = 0.
+    bfp_field(0) = 0.
+    bfp_field(1) = 0.
+    if(ibootstrap.gt.0) visc_e_field = 0.
+    psi_coil_field = 0.
+    !call destroy_auxiliary_fields
+    !call create_auxiliary_fields
+
+    call initial_conditions
+    ! combine the equilibrium and perturbed fields of linear=0
+    ! unless eqsubtract = 1
+    if(eqsubtract.eq.0) then
+      call add(field_vec, field0_vec)
+      field0_vec = 0.
+    endif
+    i_control%err_i = 0.
+    i_control%err_p_old = 0.
+    n_control%err_i = 0.
+    n_control%err_p_old = 0.
+    call reset_scalars
+    if(myrank.eq.0 .and. iprint.ge.2) print *, "  transport coefficients"
+    call define_transport_coefficients
+    if(eqsubtract.eq.1) then
+      call derived_quantities(0)
+    end if
+    call derived_quantities(1)
+    meshAdapted =1
+  end subroutine adapt_by_spr
+
 ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   subroutine adapt_by_error
 ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -327,6 +426,7 @@ subroutine adapt_by_error
     call m3dc1_mesh_getnumglobalent (0, num_node_total)
 
     !if (myrank .eq. 0) print*, "time", ntime, "current max", max_val(1), min_val(1), "mesh size before adapt", num_node_total
+
     call m3dc1_field_max(field_vec%id, max_val, min_val)
     maxPhi = max(abs(max_val(1+(u_g-1)*dofs_per_node)),abs(min_val(1+(u_g-1)*dofs_per_node)))
     maxPs =  max(abs(max_val((psi_g-1)*dofs_per_node)+1),abs(min_val((psi_g-1)*dofs_per_node)+1))
diff --git a/unstructured/auxiliary_fields.f90 b/unstructured/auxiliary_fields.f90
index 930f732b0..0961cd25c 100644
--- a/unstructured/auxiliary_fields.f90
+++ b/unstructured/auxiliary_fields.f90
@@ -36,44 +36,45 @@ subroutine create_auxiliary_fields
   use kprad_m3dc1
   implicit none
 
-  call create_field(bdotgradp)
-  call create_field(bdotgradt)
-  call create_field(torque_density_em)
-  call create_field(torque_density_ntv)
-  call create_field(chord_mask)
-  call create_field(mag_reg)
-  call create_field(ef_r)
-  call create_field(ef_phi)
-  call create_field(ef_z)
-  call create_field(ef_par)
-  call create_field(eta_j)
-  call create_field(mesh_zone)
-  call create_field(z_effective)
-  if(ikprad.eq.1) call create_field(kprad_totden)
+  call create_field(bdotgradp, "bdotgradp")
+  call create_field(bdotgradt, "bdotgradt")
+  call create_field(torque_density_em, "torque_density_em")
+  call create_field(torque_density_ntv, "torque_density_ntv")
+  call create_field(chord_mask, "chord_mask")
+  call create_field(mag_reg, "mag_reg")
+  call create_field(ef_r, "ef_r")
+  call create_field(ef_phi, "ef_phi")
+  call create_field(ef_z, "ef_z")
+  call create_field(ef_par, "ef_par")
+  call create_field(eta_j, "eta_j")
+  call create_field(mesh_zone, "mesh_zone")
+  call create_field(z_effective, "z_effective")
+
+  if(ikprad.eq.1) call create_field(kprad_totden, "kprad_totden")
   if(jadv.eq.0) then
-     call create_field(psidot)
-     call create_field(veldif)
-     call create_field(eta_jdb)
-     call create_field(bdgp)
-     call create_field(vlbdgp)
+     call create_field(psidot, "psidot")
+     call create_field(veldif, "veldif")
+     call create_field(eta_jdb, "eta_jdb")
+     call create_field(bdgp, "bdgp")
+     call create_field(vlbdgp, "vlbdgp")
   endif
   if(itemp_plot.eq.1) then
-     call create_field(vdotgradt)
-     call create_field(adv1)
-     call create_field(adv2)
-     call create_field(adv3)
-     call create_field(deldotq_perp)
-     call create_field(deldotq_par)
-     call create_field(eta_jsq)
-  call create_field(pot2_field)
-     call create_field(vpar_field)
-     call create_field(f1vplot)
-     call create_field(f1eplot)
-     call create_field(f2vplot)
-     call create_field(f2eplot)
-     call create_field(f3vplot)
-     call create_field(f3eplot)
-     call create_field(jdbobs)
+     call create_field(vdotgradt, "vdotgradt")
+     call create_field(adv1, "adv1")
+     call create_field(adv2, "adv2")
+     call create_field(adv3, "adv3")
+     call create_field(deldotq_perp, "deldotq_perp")
+     call create_field(deldotq_par, "deldotq_par")
+     call create_field(eta_jsq, "eta_jsq")
+     call create_field(pot2_field, "pot2_field")
+     call create_field(vpar_field, "vpar_field")
+     call create_field(f1vplot, "f1vplot")
+     call create_field(f1eplot, "f1eplot")
+     call create_field(f2vplot, "f2vplot")
+     call create_field(f2eplot, "f2eplot")
+     call create_field(f3vplot, "f3vplot")
+     call create_field(f3eplot, "f3eplot")
+     call create_field(jdbobs, "jdbobs")
   endif
   initialized = .true.
 end subroutine create_auxiliary_fields
@@ -150,8 +151,8 @@ subroutine calculate_pressures(ilin, pe, p, ne, nion, te, ti, ieqsub)
 
   if(myrank.eq.0 .and. iprint.ge.1) print *, ' Calculating pressures'
 
-  call create_field(pe_f)
-  call create_field(p_f)
+  call create_field(pe_f, " pe_f")
+  call create_field(p_f, "p_f")
 
   pe_f = 0.
   p_f = 0.
@@ -268,8 +269,8 @@ subroutine calculate_temperatures(ilin, te, ti, pe, p, ne, nion, ieqsub)
 
   if(myrank.eq.0 .and. iprint.ge.1) print *, ' Calculating Temperatures'
 
-  call create_field(te_f)
-  call create_field(ti_f)
+  call create_field(te_f, "te_f")
+  call create_field(ti_f, "ti_f")
 
   te_f = 0.
   ti_f = 0.
@@ -411,7 +412,7 @@ subroutine calculate_ne(ilin, nion, ne, ieqsub)
      return
   end if
   
-  call create_field(ne_f)
+  call create_field(ne_f, "ne_f")
   ne_f = 0.
 
   numelms = local_elements()
diff --git a/unstructured/field.f90 b/unstructured/field.f90
index ae5f354a8..6ed22f27d 100644
--- a/unstructured/field.f90
+++ b/unstructured/field.f90
@@ -78,16 +78,35 @@ module field
   ! ~~~~~~~~~~~~
   ! creates a field and allocates its associated size-1 vector 
   !====================================================================== 
-   subroutine create_field(f)
+   subroutine create_field(f,prefix)
     implicit none
 
     type(field_type), intent(inout) :: f
+    character(len=*), intent(in), optional :: prefix
 
     allocate(f%vec)
-    call create_vector(f%vec, 1)
+    if (present(prefix)) then
+      call create_vector(f%vec, 1, prefix)
+    else
+      call create_vector(f%vec, 1)
+    end if
     f%index = 1
   end subroutine create_field
 
+  !======================================================================
+  ! mark_field_for_solutiontransfer
+  ! ~~~~~~~~~~~~
+  ! creates a field and allocates its associated size-1 vector 
+  !====================================================================== 
+   subroutine mark_field_for_solutiontransfer(f)
+    implicit none
+
+    type(field_type), intent(inout) :: f
+    call mark_vector_for_solutiontransfer(f%vec)
+  end subroutine mark_field_for_solutiontransfer
+
+
+
   !======================================================================
   ! destroy_field
   ! ~~~~~~~~~~~~~
diff --git a/unstructured/gradshafranov.f90 b/unstructured/gradshafranov.f90
index 4dcde9b59..ffe028284 100644
--- a/unstructured/gradshafranov.f90
+++ b/unstructured/gradshafranov.f90
@@ -1045,13 +1045,13 @@ subroutine gradshafranov_solve
   numelms = local_elements()
 
   ! allocate memory for arrays
-  call create_field(b1vecini_vec)
-  call create_field(b2vecini_vec)
-  call create_field(psi_vec)
-  call create_field(fun1_vec)
-  call create_field(fun2_vec)
-  call create_field(fun3_vec)
-  call create_field(fun4_vec)
+  call create_field(b1vecini_vec, "b1vecini_vec")
+  call create_field(b2vecini_vec, "b2vecini_vec")
+  call create_field(psi_vec, "psi_vec")
+  call create_field(fun1_vec, "fun1_vec")
+  call create_field(fun2_vec, "fun2_vec")
+  call create_field(fun3_vec, "fun3_vec")
+  call create_field(fun4_vec, "fun4_vec")
 
   psi_vec = psi_field(0)
   b1vecini_vec = jphi_field
@@ -1321,8 +1321,8 @@ subroutine gradshafranov_solve
   b1vecini_vec = 0.
   b2vecini_vec = 0.
 
-  call create_field(b3vecini_vec)
-  if(irot.ne.0) call create_field(b4vecini_vec)
+  call create_field(b3vecini_vec, "b3vecini_vec")
+  if(irot.ne.0) call create_field(b4vecini_vec, "b4vecini_vec")
 
   !recalculate omega
 !  if(irot.eq.-1) call calc_omega_profile
diff --git a/unstructured/input.f90 b/unstructured/input.f90
index a485cf1d6..433721f0b 100644
--- a/unstructured/input.f90
+++ b/unstructured/input.f90
@@ -1080,6 +1080,14 @@ subroutine set_defaults
   ! Mesh adaptation
   call add_var_int("iadapt", iadapt, 0, "", adapt_grp)
 
+  ! SPR based adapt parameters
+  call add_var_int("ispradapt", ispradapt, 0, "", adapt_grp);
+  call add_var_int("isprntime", isprntime, 10, "", adapt_grp);
+  call add_var_double("isprweight", isprweight, 0.1,"", adapt_grp)
+  call add_var_double("isprmaxsize", isprmaxsize, 0.05,"", adapt_grp)
+  call add_var_int("isprrefinelevel", isprrefinelevel, 1, "", adapt_grp);
+  call add_var_int("isprcoarsenlevel", isprcoarsenlevel, -1, "", adapt_grp);
+
 
   !Micellaneous parameters or mesh adaptation
   call add_var_int("iadapt_writevtk", iadapt_writevtk, 0, "", adapt_grp)
diff --git a/unstructured/kprad_m3dc1.f90 b/unstructured/kprad_m3dc1.f90
index 61bd0ca58..053ec00d5 100644
--- a/unstructured/kprad_m3dc1.f90
+++ b/unstructured/kprad_m3dc1.f90
@@ -9,7 +9,8 @@ module kprad_m3dc1
   implicit none
 
   type(field_type), allocatable :: kprad_n(:)
-  type(field_type), allocatable, private :: kprad_temp(:)
+  ! type(field_type), allocatable, private :: kprad_temp(:)
+  type(field_type), allocatable :: kprad_temp(:)
   type(field_type) :: kprad_rad      ! power lost to line radiation
   type(field_type) :: kprad_brem     ! power lost to bremsstrahlung
   type(field_type) :: kprad_ion      ! power lost to ionization
@@ -51,6 +52,7 @@ subroutine kprad_init(ierr)
     integer, intent(out) :: ierr
 
     integer :: i
+    character(len=32) :: fname
 
     ierr = 0
     if(ikprad.eq.0) return
@@ -63,18 +65,21 @@ subroutine kprad_init(ierr)
     allocate(kprad_particle_source(0:kprad_z))
     allocate(lp_source_rate(0:kprad_z))
     do i=0, kprad_z
-       call create_field(kprad_n(i))
-       call create_field(kprad_temp(i))
-       call create_field(kprad_particle_source(i))
+       write(fname,"(A5,I2.2,A)")  "kprn", i, 0
+       call create_field(kprad_n(i), trim(fname))
+       write(fname,"(A5,I2.2,A)")  "kprt", i, 0
+       call create_field(kprad_temp(i), trim(fname))
+       write(fname,"(A5,I2.2,A)")  "kprp", i, 0
+       call create_field(kprad_particle_source(i), trim(fname))
        kprad_particle_source(i) = 0.
     end do
-    call create_field(kprad_rad)
-    call create_field(kprad_brem)
-    call create_field(kprad_ion)
-    call create_field(kprad_reck)
-    call create_field(kprad_recp)
-    call create_field(kprad_sigma_e)
-    call create_field(kprad_sigma_i)
+    call create_field(kprad_rad, "kprad_rad")
+    call create_field(kprad_brem, "kprad_brem")
+    call create_field(kprad_ion, "kprad_ion")
+    call create_field(kprad_reck, "kprad_reck")
+    call create_field(kprad_recp, "kprad_recp")
+    call create_field(kprad_sigma_e,"kprad_sigma_e")
+    call create_field(kprad_sigma_i, "kprad_sigma_i")
 
     if(ikprad_min_option.eq.2 .or. ikprad_min_option.eq.3) then
        kprad_nemin = kprad_nemin*n0_norm
diff --git a/unstructured/newpar.f90 b/unstructured/newpar.f90
index 3ea57039d..bf515ada9 100644
--- a/unstructured/newpar.f90
+++ b/unstructured/newpar.f90
@@ -49,6 +49,7 @@ Program Reducedquintic
   character*256 :: arg, solveroption_filename
   integer :: ip
   character(len=32) :: mesh_file_name
+  logical :: update_mesh
 #ifdef _OPENACC
   integer :: num_devices
 #endif
@@ -160,6 +161,14 @@ Program Reducedquintic
   if(myrank.eq.0) print *, ' Reading input'
   call input
 
+!if using SCOREC set adapt verbosity output if iprint.ge.1
+#ifdef USESCOREC
+  if (iprint.ge.1) then
+    call m3dc1_domain_verbosity(1) ! 0 for non-verbose outputs
+  end if
+#endif
+
+
   ! load mesh
   if(myrank.eq.0 .and. iprint.ge.1) print *, ' Loading mesh nplane='
   if(myrank==0 .and. nplanes.gt.1) call parse_solver_options(nplanes, trim(solveroption_filename)//PETSC_NULL_CHARACTER)
@@ -170,7 +179,6 @@ Program Reducedquintic
 #endif
 
   call load_mesh
-  
 !  call print_node_data
 !  call safestop(1)
 
@@ -208,7 +216,7 @@ Program Reducedquintic
 
   call calc_wall_dist
 
-  call init_hyperv_mat
+  !call init_hyperv_mat                ! Dont call this here in MeshAdapt 
   
   
   ! Set initial conditions either from restart file
@@ -341,6 +349,7 @@ Program Reducedquintic
   endif
 
   ! output initial conditions
+  call marker ! mark the fields necessary for solution transfer
   call output
 
   ! if there are no timesteps to calculate, then skip time loop
@@ -361,7 +370,7 @@ Program Reducedquintic
   do ntime=ntime+1,ntimemax
 
      if(myrank.eq.0) print *, 'TIME STEP: ', ntime
-
+     call init_hyperv_mat
      ! check for error
      if(ekin.ne.ekin .or. emag.ne.emag) then
         print *, "Error: energy is NaN"
@@ -432,14 +441,31 @@ Program Reducedquintic
      if(myrank.eq.0 .and. iprint.ge.1) print *, " Writing output."
      call output
 
-      if (iadapt .gt. 1) then
+
+    ! for now call spr adapt every 10 time steps
+    if (ispradapt .eq. 1) then
+      if (mod(ntime, isprntime) .eq. 0) then
+        write(mesh_file_name,"(A11,A)") 'beforeadapt', 0
+        call m3dc1_mesh_write (mesh_file_name,0,ntime)
+        ! if update_mesh is true
+        ! the (2D) part.smb will be updated (overwritten) inside adapt_by_spr
+        ! the frequency of update is the same as frequency of output time slices
+        update_mesh = .false.
+        if(mod(ntime-ntime0,ntimepr).eq.0) then
+          update_mesh = .true.
+        end if
+        call adapt_by_spr(field_vec%id, psi_g, ntime, isprweight, isprmaxsize, isprrefinelevel, isprcoarsenlevel, update_mesh)
+      endif
+    endif
+
+    if (iadapt .gt. 1 .and. ispradapt .eq.0) then
       ! adapt_flag=1 if
       !(1) iadapt_ntime(N)>0 -- run adapt_by_error at the end of every N time steps
       !(2) non-linear & iadapt_ntime=0 -- run adapt_by_error at the end of every time step
       !(3) linear, adapt_ke>0 & ekin>adapt_ke -- run adapt_by_error in this time step  
         call diagnose_adapt(adapt_flag)
        if(adapt_flag .eq. 1) call adapt_by_error
-     endif
+    endif
   enddo ! ntime
 
   if(myrank.eq.0 .and. iprint.ge.1) print *, "Done time loop."
@@ -700,7 +726,7 @@ subroutine derived_quantities(ilin)
      if(linear.eq.1) then 
         if(ntime.eq.ntime0) call lcfs(psi_field(0))
      else
-        call create_field(psi_temp)
+        call create_field(psi_temp, "psi_temp")
         psi_temp = psi_field(0)
         call add_field_to_field(psi_temp, psi_field(1))
         call lcfs(psi_temp)
@@ -724,7 +750,7 @@ subroutine derived_quantities(ilin)
           endif
         endif
      else
-        call create_field(te_temp)
+        call create_field(te_temp, "te_temp")
         te_temp = te_field(0)
         call add_field_to_field(te_temp, te_field(1))
         if(ifixed_temax .eq. 0) then
@@ -1307,7 +1333,7 @@ subroutine space(ifirstcall)
 #ifdef USESCOREC
   if(ifirstcall .eq. 1) then
      do i=1, num_fields
-       write(field_name,"(I2,A)")  i,0
+       write(field_name,"(A3,I0,A)")  "mat", i, 0
 #ifdef USECOMPLEX
        call m3dc1_field_create (i, trim(field_name), i, 1, dofs_per_node)
 #else
@@ -1316,7 +1342,6 @@ subroutine space(ifirstcall)
      end do
   endif ! on firstcall
 #endif
-  
   numelms = local_elements()
 
 ! arrays defined at all vertices
@@ -1325,48 +1350,54 @@ subroutine space(ifirstcall)
      if(myrank.eq.0 .and. iprint.ge.1) print *, 'Allocating...'
 
      ! Physical Variables
-     call create_vector(field_vec , num_fields)
-     call create_vector(field0_vec, num_fields)
+     call create_vector(field_vec , num_fields, "field_vec")
+     call create_vector(field0_vec, num_fields, "field_vec0")
      !if(iadapt .ne. 0) then
-        call create_vector(field_vec_pre, 2)
+        call create_vector(field_vec_pre, 2, "field_vec_pre")
      !end if
 
+     call mark_vector_for_solutiontransfer(field_vec)
+     call mark_vector_for_solutiontransfer(field0_vec)
+     call mark_vector_for_solutiontransfer(field_vec_pre)
+
      ! Auxiliary Variables
-     call create_field(jphi_field)
-     call create_field(resistivity_field)
-     call create_field(kappa_field)
-     call create_field(kappar_field)
-     call create_field(denm_field)
-     call create_field(visc_field)
-     call create_field(visc_c_field)
-     if(ipforce.gt.0) call create_field(pforce_field)
-     if(ipforce.gt.0) call create_field(pmach_field)
-     if(density_source) call create_field(sigma_field)
-     if(momentum_source) call create_field(Fphi_field)
-     if(heat_source) call create_field(Q_field)
-     if(icd_source.gt.0) call create_field(cd_field)
+     call create_field(jphi_field, "jphi")
+     !call create_field(vor_field, "vor")
+     !call create_field(com_field, "com")
+     call create_field(resistivity_field, "resistivity")
+     call create_field(kappa_field, "kappa")
+     call create_field(kappar_field, "kappar")
+     call create_field(denm_field, "denm")
+     call create_field(visc_field, "visc")
+     call create_field(visc_c_field, "visc_c")
+     if(ipforce.gt.0) call create_field(pforce_field, "pforce")
+     if(ipforce.gt.0) call create_field(pmach_field, "pmach")
+     if(density_source) call create_field(sigma_field, "sigma")
+     if(momentum_source) call create_field(Fphi_field, "Fphi")
+     if(heat_source) call create_field(Q_field, "Q")
+     if(icd_source.gt.0) call create_field(cd_field, "cd")
      if(rad_source) then
-        call create_field(Totrad_field)
-        call create_field(Linerad_field)
-        call create_field(Bremrad_field)
-        call create_field(Ionrad_field)
-        call create_field(Reckrad_field)
-        call create_field(Recprad_field)
+        call create_field(Totrad_field, "Torad")
+        call create_field(Linerad_field, "Linerad")
+        call create_field(Bremrad_field, "Bremrad")
+        call create_field(Ionrad_field, "Ionrad")
+        call create_field(Reckrad_field, "Reckrad")
+        call create_field(Recprad_field, "Recprad")
      end if
-     call create_field(bf_field(0))
-     call create_field(bf_field(1))
-     call create_field(bfp_field(0))
-     call create_field(bfp_field(1))
-     if(ibootstrap.gt.0) call create_field(visc_e_field)
+     call create_field(bf_field(0), "bf0")
+     call create_field(bf_field(1), "bf1")
+     call create_field(bfp_field(0), "bfp0")
+     call create_field(bfp_field(1), "bfp1")
+     if(ibootstrap.gt.0) call create_field(visc_e_field, "visc_e")
 
-     call create_field(psi_coil_field)
+     call create_field(psi_coil_field, "psi_coil")
 
      ! create external fields
      if(extsubtract.eq.1) then
-        call create_field(psi_ext)
-        call create_field(bz_ext)
-        call create_field(bf_ext)
-        call create_field(bfp_ext)
+        call create_field(psi_ext, "pis_ext")
+        call create_field(bz_ext, "bz_ext")
+        call create_field(bf_ext, "bf_ext")
+        call create_field(bfp_ext, "bfp_ext")
         use_external_fields = .true.
      end if
 
diff --git a/unstructured/output.f90 b/unstructured/output.f90
index b5af40767..d75da496f 100644
--- a/unstructured/output.f90
+++ b/unstructured/output.f90
@@ -42,6 +42,27 @@ subroutine finalize_output
     if(ier.ne.0) print *, 'Error finalizing HDF5:',ier
   end subroutine finalize_output
 
+  ! ======================================================================
+  ! marker
+  ! ~~~~~~
+  !
+  !
+  ! ======================================================================
+  subroutine marker
+    use basic
+    use hdf5_output
+    use diagnostics
+    use auxiliary_fields
+
+    implicit none
+
+    include 'mpif.h'
+
+    integer :: ier,i
+    call mark_fields(0);
+  end subroutine marker
+
+
   ! ======================================================================
   ! output
   ! ~~~~~~
@@ -98,6 +119,7 @@ subroutine output
        endif
     endif
 
+
     if(itimer.eq.1) then
        if(myrank.eq.0) call second(tstart)
        if(myrank.eq.0 .and. iprint.ge.2) print *, "  writing timings"
@@ -1236,6 +1258,310 @@ subroutine output_fields(time_group_id, equilibrium, error)
 end subroutine output_fields
 
 
+! mark_fields (for solution transfer during adapt)
+! =============
+subroutine mark_fields(equilibrium)
+  use hdf5
+  use hdf5_output
+  use basic
+  use arrays
+  use time_step
+  use auxiliary_fields
+  use transport_coefficients
+  use kprad_m3dc1
+
+  implicit none
+
+  integer :: i, nelms, ilin
+  integer, intent(in) :: equilibrium
+  ilin = 1 - equilibrium
+
+  ! psi_plasma
+  if(icsubtract.eq.1 .or. &
+       (extsubtract.eq.1 .and. (ilin.eq.1 .or. eqsubtract.eq.0))) then
+     !psi_field(ilin)
+  end if
+
+  ! psi
+  if(icsubtract.eq.1 .and. (ilin.eq.0 .or. eqsubtract.eq.0)) then
+     !psi_coil_field
+     call mark_field_for_solutiontransfer(psi_coil_field)
+  endif
+  if(extsubtract.eq.1 .and. (ilin.eq.1 .or. eqsubtract.eq.0)) then
+     !psi_ext
+     call mark_field_for_solutiontransfer(psi_ext)
+  end if
+
+  ! u
+  ! u_field(ilin)
+
+#if defined(USE3D) || defined(USECOMPLEX)
+  ! electrostatic potential
+  if(jadv.eq.0) then
+     ! e_field(ilin)
+  endif
+#endif
+
+  ! I
+  ! bz_field(ilin)
+  if(extsubtract.eq.1 .and. (ilin.eq.1 .or. eqsubtract.eq.0)) then 
+     ! bz_ext
+     call mark_field_for_solutiontransfer(bz_ext)
+  end if
+
+  ! BF
+  if(ifout.eq.1) then
+     ! bf_field(ilin)
+     if(extsubtract.eq.1 .and. (ilin.eq.1 .or. eqsubtract.eq.0)) then 
+        ! bf_ext
+        call mark_field_for_solutiontransfer(bf_ext)
+     end if
+  endif
+
+  ! BFP
+  if(ifout.eq.1) then
+     ! bfp_field(ilin)
+     call mark_field_for_solutiontransfer(bfp_field(ilin))
+     if(extsubtract.eq.1 .and. (ilin.eq.1 .or. eqsubtract.eq.0)) then 
+        ! bfp_ext
+        call mark_field_for_solutiontransfer(bfp_ext)
+     end if
+  endif
+
+  !  vz_field(ilin)
+  !  pe_field(ilin)
+  !   p_field(ilin)
+  ! chi_field(ilin)
+  ! den_field(ilin)
+  !  te_field(ilin)
+  !  ti_field(ilin)
+  !  ne_field(ilin)
+
+  if(icsubtract.eq.1) then
+     ! psi_coil_field
+     call mark_field_for_solutiontransfer(psi_coil_field)
+  end if
+
+#ifdef USEPARTICLES
+  if (kinetic.eq.1) then
+     if (associated(p_i_perp%vec)) then
+        !Perpendicular component of hot ion pressure tensor
+        ! p_i_perp
+        call mark_vector_for_solutiontransfer(p_i_perp%vec)
+     endif
+
+     if (associated(p_i_par%vec)) then
+        !Parallel component of hot ion pressure tensor
+        ! p_i_par
+        call mark_vector_for_solutiontransfer(p_i_par%vec)
+     endif
+  endif
+#endif
+
+  if(use_external_fields) then
+     ! psi_ext
+     ! bz_ext
+     ! bf_ext
+     ! bfp_ext
+     call mark_field_for_solutiontransfer(psi_ext)
+     call mark_field_for_solutiontransfer(bz_ext)
+     call mark_field_for_solutiontransfer(bf_ext)
+     call mark_field_for_solutiontransfer(bfp_ext)
+  endif
+
+  if(ikprad.eq.1) then
+     do i=0, kprad_z
+        ! kprad_n(i)
+        ! kprad_particle_source(i)
+        call mark_field_for_solutiontransfer(kprad_n(i))
+        ! call mark_field_for_solutiontransfer(kprad_temp(i))
+        call mark_field_for_solutiontransfer(kprad_particle_source(i))
+     end do
+     ! kprad_sigma_e,
+     ! kprad_sigma_i,
+     ! kprad_rad,
+     ! kprad_brem,
+     ! kprad_ion,
+     ! kprad_reck,
+     ! kprad_recp,
+     call mark_field_for_solutiontransfer(kprad_sigma_e)
+     call mark_field_for_solutiontransfer(kprad_sigma_i)
+     call mark_field_for_solutiontransfer(kprad_rad)
+     call mark_field_for_solutiontransfer(kprad_brem)
+     call mark_field_for_solutiontransfer(kprad_ion)
+     call mark_field_for_solutiontransfer(kprad_reck)
+     call mark_field_for_solutiontransfer(kprad_recp)
+  end if
+
+  ! transport coefficients do not change with time in linear calculations
+  ! so don't output linear perturbations
+  if(iwrite_transport_coeffs.eq.1 .and. &
+     (linear.eq.0 .or. equilibrium.eq.1)) then
+     !  resistivity_field
+     !  visc_field
+     !  visc_c_field
+     !  kappa_field
+     !  denm_field
+     call mark_field_for_solutiontransfer(resistivity_field)
+     call mark_field_for_solutiontransfer(visc_field)
+     call mark_field_for_solutiontransfer(visc_c_field)
+     call mark_field_for_solutiontransfer(kappa_field)
+     call mark_field_for_solutiontransfer(denm_field)
+
+     ! poloidal force and mach number
+     if(ipforce.gt.0) then
+        ! pforce_field
+        ! pmach_field
+        call mark_field_for_solutiontransfer(pforce_field)
+        call mark_field_for_solutiontransfer(pmach_field)
+     endif
+
+     if(ibootstrap.gt.0) then
+        ! visc_e_field
+        call mark_field_for_solutiontransfer(visc_e_field)
+     endif
+  end if !(iwrite_transport_coeffs.eq.1)
+
+  if(irunaway.ne.0) then
+     ! nre_field
+     call mark_field_for_solutiontransfer(nre_field(ilin))
+  end if
+
+  if(iwrite_aux_vars.eq.1) then
+    ! wall_dist
+    ! jphi_field
+    ! vor_field
+    ! com_field
+    ! torque_density_em
+    ! torque_density_ntv
+    ! bdotgradp, nelms
+    ! bdotgradt, nelms
+    ! z_effective
+    call mark_field_for_solutiontransfer(wall_dist)
+    call mark_field_for_solutiontransfer(jphi_field)
+    !call mark_field_for_solutiontransfer(vor_field)
+    !call mark_field_for_solutiontransfer(com_field)
+    call mark_field_for_solutiontransfer(torque_density_ntv)
+    call mark_field_for_solutiontransfer(bdotgradp)
+    call mark_field_for_solutiontransfer(bdotgradt)
+    call mark_field_for_solutiontransfer(z_effective)
+    if(ikprad.eq.1) then
+       ! kprad_totden
+       call mark_field_for_solutiontransfer(kprad_totden)
+    end if
+    if(itemp_plot .eq. 1) then
+       ! vdotgradt
+       ! adv1
+       ! adv2
+       ! adv3
+       ! deldotq_perp
+       ! deldotq_par
+       ! eta_jsq
+       ! vpar_field
+       ! f1vplot
+       ! f1eplot
+       ! f2vplot
+       ! f2eplot
+       ! f3vplot
+       ! f3eplot
+       ! jdbobs
+       ! pot2_field
+       call mark_field_for_solutiontransfer(vdotgradt)
+       call mark_field_for_solutiontransfer(adv1)
+       call mark_field_for_solutiontransfer(adv2)
+       call mark_field_for_solutiontransfer(adv3)
+       call mark_field_for_solutiontransfer(deldotq_perp)
+       call mark_field_for_solutiontransfer(deldotq_par)
+       call mark_field_for_solutiontransfer(eta_jsq)
+       call mark_field_for_solutiontransfer(vpar_field)
+       call mark_field_for_solutiontransfer(f1vplot)
+       call mark_field_for_solutiontransfer(f1eplot)
+       call mark_field_for_solutiontransfer(f2eplot)
+       call mark_field_for_solutiontransfer(f2vplot)
+       call mark_field_for_solutiontransfer(f3eplot)
+       call mark_field_for_solutiontransfer(f3vplot)
+       call mark_field_for_solutiontransfer(jdbobs)
+       call mark_field_for_solutiontransfer(pot2_field)
+
+       if(jadv.eq.0) then
+          ! psidot
+          ! veldif
+          ! eta_jdb
+          ! bdgp
+          ! vlbdgp
+          call mark_field_for_solutiontransfer(psidot)
+          call mark_field_for_solutiontransfer(veldif)
+          call mark_field_for_solutiontransfer(eta_jdb)
+          call mark_field_for_solutiontransfer(bdgp)
+          call mark_field_for_solutiontransfer(vlbdgp)
+       endif
+
+    endif
+
+    ! sigma
+    if(density_source) then
+       ! sigma_field
+       call mark_field_for_solutiontransfer(sigma_field)
+    endif
+
+    ! momentum source
+    if(momentum_source) then
+       ! Fphi_field
+       call mark_field_for_solutiontransfer(Fphi_field)
+    endif
+
+    ! heat source
+    if(heat_source) then
+       ! Q_field
+       call mark_field_for_solutiontransfer(Q_field)
+    endif
+
+    ! radiation source
+    if(rad_source) then
+       ! Totrad_field
+       ! Linerad_field
+       ! Bremrad_field
+       ! Ionrad_field
+       ! Reckrad_field
+       ! Recprad_field
+       call mark_field_for_solutiontransfer(Totrad_field)
+       call mark_field_for_solutiontransfer(Linerad_field)
+       call mark_field_for_solutiontransfer(Bremrad_field)
+       call mark_field_for_solutiontransfer(Ionrad_field)
+       call mark_field_for_solutiontransfer(Reckrad_field)
+       call mark_field_for_solutiontransfer(Recprad_field)
+    endif
+
+    ! current drive source
+    if(icd_source.gt.0) then
+       ! cd_field
+       call mark_field_for_solutiontransfer(cd_field)
+    endif
+
+    if(xray_detector_enabled.eq.1) then
+       ! chord_mask
+       call mark_field_for_solutiontransfer(chord_mask)
+    end if
+
+    ! magnetic region
+    ! mag_reg
+    call mark_field_for_solutiontransfer(mag_reg)
+
+    ! mesh zone
+    ! mesh_zone
+    call mark_field_for_solutiontransfer(mesh_zone)
+
+    ! electric_field
+    call mark_field_for_solutiontransfer(ef_r)
+    call mark_field_for_solutiontransfer(ef_phi)
+    call mark_field_for_solutiontransfer(ef_z)
+    call mark_field_for_solutiontransfer(ef_par)
+    call mark_field_for_solutiontransfer(eta_j)
+
+ endif !(iwrite_aux_vars.eq.1)
+
+end subroutine mark_fields
+
 ! hdf5_write_keharmonics
 ! ==================
 subroutine hdf5_write_keharmonics(error)
diff --git a/unstructured/romulus.mk b/unstructured/romulus.mk
new file mode 100644
index 000000000..4e49a9466
--- /dev/null
+++ b/unstructured/romulus.mk
@@ -0,0 +1,84 @@
+FOPTS = -c -fdefault-real-8 -fdefault-double-8 -finit-local-zero -Wall -cpp -DPETSC_VERSION=313 -DUSEBLAS $(OPTS) 
+CCOPTS  = -c -O -DPETSC_VERSION=313 -DDEBUG
+
+ifeq ($(OPT), 1)
+  FOPTS  := $(FOPTS) -O2 
+  CCOPTS := $(CCOPTS) -O
+else
+  FOPTS := $(FOPTS) -g noarg_temp_created 
+endif
+
+ifeq ($(PAR), 1)
+  FOPTS := $(FOPTS) -DUSEPARTICLES
+endif
+
+CC = $(MPICH_ROOT)/bin/mpicc
+CPP = $(MPICH_ROOT)/bin/mpicxx
+F90 = $(MPICH_ROOT)/bin/mpif90
+F77 = $(MPICH_ROOT)/bin/mpif90
+LOADER = $(MPICH_ROOT)/bin/mpif90
+FOPTS := $(FOPTS)
+
+F90OPTS = $(F90FLAGS) $(FOPTS) 
+F77OPTS = $(F77FLAGS) $(FOPTS)
+
+PETSC_DIR=/lore/seol/petsc-3.13.5
+ifeq ($(COM), 1)
+  PETSC_ARCH=cplx-gcc4.8.5-v5m6xwi-mpich3.2.1-geowaxe
+else
+  PETSC_ARCH=real-gcc4.8.5-v5m6xwi-mpich3.2.1-geowaxe
+endif
+
+SCOREC_BASE_DIR=/lore/seol/trash
+PUMI_DIR=$(SCOREC_BASE_DIR)
+PUMI_LIB = -lpumi -lapf -lapf_zoltan -lcrv -lsam -lspr -lmth -lgmi -lma -lmds -lparma -lpcu -lph -llion
+
+ifdef SCORECVER
+  SCOREC_DIR=$(SCOREC_BASE_DIR)/$(SCORECVER)
+else
+  SCOREC_DIR=$(SCOREC_BASE_DIR)
+endif
+
+ifeq ($(COM), 1)
+  M3DC1_SCOREC_LIB=-lm3dc1_scorec_complex
+else
+  M3DC1_SCOREC_LIB=-lm3dc1_scorec
+endif
+
+SCOREC_LIB = -L$(SCOREC_DIR)/lib $(M3DC1_SCOREC_LIB) \
+            -Wl,--start-group,-rpath,$(PUMI_DIR)/lib -L$(PUMI_DIR)/lib \
+           $(PUMI_LIB) -Wl,--end-group
+
+ifeq ($(ENABLE_ZOLTAN), 0)
+  ZOLTAN_LIB=
+  BIN_POSTFIX := $(BIN_POSTFIX)-nozoltan
+  CCOPTS := $(CCOPTS) -DDISABLE_ZOLTAN
+else
+  ZOLTAN_LIB=-lzoltan
+endif 
+
+PETSC_WITH_EXTERNAL_LIB = -L${PETSC_DIR}/${PETSC_ARCH}/lib -Wl,-rpath,$(PETSC_DIR)/$(PETSC_ARCH)/lib -L$(PETSC_DIR)/$(PETSC_ARCH)/lib -L/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-4.8.5/mpich-3.2.1-geowaxeufjfz7vgpszmfrm4l6rp4rnxo/lib -Wl,-rpath,/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-rhel7_4.8.5/gcc-4.8.5-v5m6xwipewvkufflu3zvxs2x6jv4libr/lib:/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-rhel7_4.8.5/gcc-4.8.5-v5m6xwipewvkufflu3zvxs2x6jv4libr/lib64 -L/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-rhel7_4.8.5/gcc-4.8.5-v5m6xwipewvkufflu3zvxs2x6jv4libr/lib64 -L/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-rhel7_4.8.5/gcc-4.8.5-v5m6xwipewvkufflu3zvxs2x6jv4libr/lib/gcc/x86_64-unknown-linux-gnu/4.8.5 -L/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-4.8.5/hdf5-1.10.4-wousvnks7b7aplenali7jg5bsh243hyd/lib -L/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-4.8.5/zlib-1.2.11-vhzh5cfaki5lx5sjuth5iuojq5azdkbd/lib -L/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-rhel7_4.8.5/gcc-4.8.5-v5m6xwipewvkufflu3zvxs2x6jv4libr/lib -Wl,-rpath,/opt/scorec/spack/install/linux-rhel7-x86_64/gcc-4.8.5/mpich-3.2.1-geowaxeufjfz7vgpszmfrm4l6rp4rnxo/lib -lpetsc -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lpthread -lscalapack -lsuperlu -lsuperlu_dist -lfftw3_mpi -lfftw3 -lflapack -lfblas $(ZOLTAN_LIB) -lpthread -lX11 -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5 -lparmetis -lmetis -lm -lz -lstdc++ -ldl -lmpifort -lmpi -lgfortran -lm -lgfortran -lm -lgcc_s -lquadmath -lstdc++ -ldl
+
+LIBS = 	\
+	$(SCOREC_LIB) \
+        $(PETSC_WITH_EXTERNAL_LIB) \
+	-L$(GSL_ROOT)/lib -lgsl -lgslcblas 
+
+INCLUDE = -I$(PETSC_DIR)/include \
+        -I$(PETSC_DIR)/$(PETSC_ARCH)/include \
+        -I$(GSL_ROOT)/include
+
+%.o : %.c
+	$(CC)  $(CCOPTS) $(INCLUDE) $< -o $@
+
+%.o : %.cpp
+	$(CPP) $(CCOPTS) $(INCLUDE) $< -o $@
+
+%.o: %.f
+	$(F77) $(F77OPTS) $(INCLUDE) $< -o $@
+
+%.o: %.F
+	$(F77) $(F77OPTS) $(INCLUDE) $< -o $@
+
+%.o: %.f90
+	$(F90) $(F90OPTS) $(INCLUDE) -fpic $< -o $@
diff --git a/unstructured/scorec_vector.f90 b/unstructured/scorec_vector.f90
index c1babb633..75821616f 100644
--- a/unstructured/scorec_vector.f90
+++ b/unstructured/scorec_vector.f90
@@ -33,6 +33,11 @@ module scorec_vector_mod
      module procedure scorec_vector_create
   end interface
 
+  interface mark_vector_for_solutiontransfer
+     module procedure scorec_vector_mark_for_solutiontransfer
+  end interface
+
+
   interface destroy_vector
      module procedure scorec_vector_destroy
   end interface
@@ -482,17 +487,24 @@ end subroutine scorec_vector_insert_complex
   !======================================================================
   ! create
   ! ~~~~~~
-  ! creates a vector of size n
+  ! creates a vector of size n with name prefix (optional)
   !======================================================================
-  subroutine scorec_vector_create(f,n)
+  subroutine scorec_vector_create(f,n,prefix)
     implicit none
 
     type(scorec_vector), intent(inout) :: f
     integer, intent(in) :: n
+    character(len=*), intent(in), optional :: prefix
     integer :: dataType
     character(len=32) :: field_name
     call m3dc1_field_genid (f%id)
-    write(field_name,"(I0,A)") f%id, 0
+
+    if(present(prefix)) then
+      write(field_name,"(A16,A)") prefix, 0
+      field_name = adjustl(field_name)
+    else
+      write(field_name,"(I0,A)") f%id, 0
+    end if
     dataType = 0
 #ifdef USECOMPLEX
     dataType = 1
@@ -501,6 +513,14 @@ subroutine scorec_vector_create(f,n)
     f%isize = n
   end subroutine scorec_vector_create
 
+  subroutine scorec_vector_mark_for_solutiontransfer(f)
+    implicit none
+
+    type(scorec_vector), intent(inout) :: f
+    call m3dc1_mark_for_solutiontransfer(f%id)
+  end subroutine scorec_vector_mark_for_solutiontransfer
+
+
   !======================================================================
   ! destroy
   ! ~~~~~~~
diff --git a/unstructured/stellar.mk b/unstructured/stellar.mk
index 102287ab6..d2431b6cc 100644
--- a/unstructured/stellar.mk
+++ b/unstructured/stellar.mk
@@ -47,21 +47,26 @@ PETSC_WITH_EXTERNAL_LIB = -L${PETSC_DIR}/${PETSC_ARCH}/lib \
 	-lifport -lifcoremt_pic -limf -lsvml -lm -lipgo -lirc -lgcc_s -lirc_s -lquadmath \
 	-lstdc++ -ldl
 
-SCOREC_BASE_DIR=/projects/M3DC1/scorec/$(MPIVER)/$(PETSCVER)/202209
+SCOREC_BASE_DIR=/scratch/gpfs/ur8212/install/core/core_230210
 SCOREC_UTIL_DIR=$(SCOREC_BASE_DIR)/bin
 ifdef SCORECVER
-  SCOREC_DIR=$(SCOREC_BASE_DIR)/$(SCORECVER)
+  SCOREC_DIR=/scratch/gpfs/ur8212/install/M3DC1/M3DC1_230509/m3dc1_scorec
 else
-  SCOREC_DIR=$(SCOREC_BASE_DIR)
+  SCOREC_DIR=/scratch/gpfs/ur8212/install/M3DC1/M3DC1_230509/m3dc1_scorec
 endif
 
 ZOLTAN_LIB=-L$(PETSC_DIR)/$(PETSC_ARCH)/lib -lzoltan
 
-SCOREC_LIBS= -L$(SCOREC_DIR)/lib $(M3DC1_SCOREC_LIB) \
+#SCOREC_LIBS= -L$(SCOREC_DIR)/lib $(M3DC1_SCOREC_LIB) \
              -Wl,--start-group,-rpath,$(SCOREC_DIR)/lib -L$(SCOREC_DIR)/lib \
              -lpumi -lapf -lapf_zoltan -lgmi -llion -lma -lmds -lmth -lparma \
              -lpcu -lph -lsam -lspr -lcrv -Wl,--end-group
 
+SCOREC_LIBS= -L$(SCOREC_DIR)/lib $(M3DC1_SCOREC_LIB) \
+             -Wl,--start-group,-rpath,$(SCOREC_DIR)/lib -L$(SCOREC_BASE_DIR)/lib \
+             -lpumi -lapf -lapf_zoltan -lgmi -llion -lma -lmds -lmth -lparma \
+             -lpcu -lph -lsam -lspr -lcrv -Wl,--end-group
+
 LIBS = 	-L$(I_MPI_ROOT)/lib -lmpicxx\
 	$(SCOREC_LIBS) \
         $(ZOLTAN_LIB) \
diff --git a/unstructured/time_step_split.f90 b/unstructured/time_step_split.f90
index f4d8d827c..21d458c62 100644
--- a/unstructured/time_step_split.f90
+++ b/unstructured/time_step_split.f90
@@ -88,38 +88,59 @@ subroutine initialize_timestep_split
   
 
     ! Vectors
-    call create_vector(phi_vec,      vecsize_phi)
-    call create_vector(phip_vec,     vecsize_phi)
-    call create_vector(q4_vec,       vecsize_phi)
+    call create_vector(phi_vec,      vecsize_phi, "phi_vec")
+    call create_vector(phip_vec,     vecsize_phi, "phip_vec")
+    call create_vector(q4_vec,       vecsize_phi, "q4_vec")
+    call mark_vector_for_solutiontransfer(phi_vec)
+    call mark_vector_for_solutiontransfer(phip_vec)
+    call mark_vector_for_solutiontransfer(q4_vec)
      
-    call create_vector(b1_phi, vecsize_phi)
-    call create_vector(b2_phi, vecsize_phi)
-
-    call create_vector(vel_vec,      vecsize_vel)
-    call create_vector(veln_vec,     vecsize_vel)
-    call create_vector(r4_vec,       vecsize_vel)
+    call create_vector(b1_phi, vecsize_phi, "b1_phi")
+    call create_vector(b2_phi, vecsize_phi, "b2_phi")
+    call mark_vector_for_solutiontransfer(b1_phi)
+    call mark_vector_for_solutiontransfer(b2_phi)
+
+    call create_vector(vel_vec,      vecsize_vel, "vel_vec")
+    call create_vector(veln_vec,     vecsize_vel, "veln_vec")
+    call create_vector(r4_vec,       vecsize_vel, "r4_vec")
+    call mark_vector_for_solutiontransfer(vel_vec)
+    call mark_vector_for_solutiontransfer(veln_vec)
+    call mark_vector_for_solutiontransfer(r4_vec)
 
     if(ipres.eq.1 .or. ipressplit.eq.1) then
-       call create_vector(pres_vec,    vecsize_p)
-       call create_vector(qp4_vec,     vecsize_p)
+       call create_vector(pres_vec,    vecsize_p, "pres_vec")
+       call create_vector(qp4_vec,     vecsize_p, "qp4_vec")
+       call mark_vector_for_solutiontransfer(pres_vec)
+       call mark_vector_for_solutiontransfer(qp4_vec)
     endif
 
-    call create_vector(den_vec,    vecsize_n)
-    call create_vector(denold_vec, vecsize_n)
-    call create_vector(ne_vec,     1)
-    call create_vector(neold_vec,  1)
-    call create_vector(qn4_vec,    vecsize_n)
+    call create_vector(den_vec,    vecsize_n, "den_vec")
+    call create_vector(denold_vec, vecsize_n, "denold_vec")
+    call create_vector(ne_vec,     1, "ne_vec")
+    call create_vector(neold_vec,  1, "neold_vec")
+    call create_vector(qn4_vec,    vecsize_n, "qn4_vec")
+    call mark_vector_for_solutiontransfer(den_vec)
+    call mark_vector_for_solutiontransfer(denold_vec)
+    call mark_vector_for_solutiontransfer(ne_vec)
+    call mark_vector_for_solutiontransfer(neold_vec)
+    call mark_vector_for_solutiontransfer(qn4_vec)
 
     if(irunaway .gt. 0) then
-      call create_vector(nre_vec,    vecsize_n)
-      call create_vector(nreold_vec, vecsize_n)
-      call create_vector(qn5_vec,    vecsize_n)
+      call create_vector(nre_vec,    vecsize_n, "nre_vec")
+      call create_vector(nreold_vec, vecsize_n, "nreold_vec")
+      call create_vector(qn5_vec,    vecsize_n, "qn5_vec")
+      call mark_vector_for_solutiontransfer(nre_vec)
+      call mark_vector_for_solutiontransfer(nreold_vec)
+      call mark_vector_for_solutiontransfer(qn5_vec)
     endif
 
-    call create_vector(b1_vel, vecsize_vel)
-    call create_vector(b2_vel, vecsize_vel)
+    call create_vector(b1_vel, vecsize_vel, "b1_vel")
+    call create_vector(b2_vel, vecsize_vel, "b2_vel")
+    call mark_vector_for_solutiontransfer(b1_vel)
+    call mark_vector_for_solutiontransfer(b2_vel)
 
-    call create_vector(pret_vec,    vecsize_t)
+    call create_vector(pret_vec,    vecsize_t, "pret_vec")
+    call mark_vector_for_solutiontransfer(pret_vec)
 
 call PetscLogStagePush(stageA,jer)
     ! Matrices associated with velocity advance
diff --git a/unstructured/wall.f90 b/unstructured/wall.f90
index 56881d4cc..a9f968148 100644
--- a/unstructured/wall.f90
+++ b/unstructured/wall.f90
@@ -26,7 +26,7 @@ subroutine calc_wall_dist
     call set_matrix_index(wall_matrix, wall_mat_index)
     call create_mat(wall_matrix, 1, 1, icomplex, 1)
 
-    call create_field(wall_dist)
+    call create_field(wall_dist, "wall_dist")
     wall_dist = 0.
 
     ibound = ior(BOUNDARY_DIRICHLET,BOUNDARY_NEUMANN)