diff --git a/Fltk/GUI_Extras.cpp b/Fltk/GUI_Extras.cpp
index cdc602bbcb21a36557dc63f4b28a603221c7f4b2..896c0625f056b896b599dae9a45611a1d440e161 100644
--- a/Fltk/GUI_Extras.cpp
+++ b/Fltk/GUI_Extras.cpp
@@ -1181,8 +1181,8 @@ struct CGNSWriteDialog
       roundButton1GCatFace->value();
     CTX.mesh.cgns_options.writeBC = checkButtonWriteBC->value();
     CTX.mesh.cgns_options.bocoLocation = roundButton1BCatFace->value();
-    CTX.mesh.cgns_options.writeNormals = checkButtonWriteNormals->value();
-    CTX.mesh.cgns_options.normalSource = roundButton1NormalElem->value();
+    CTX.mesh.cgns_options.normalSource = (checkButtonWriteNormals->value()) ?
+       roundButton1NormalElem->value() + 1 : 0;
     CTX.mesh.cgns_options.vectorDim = choiceVecDim->value() + 2;
     CTX.mesh.cgns_options.writeUserDef = checkButtonUnknownUserDef->value();
   }
@@ -1194,7 +1194,7 @@ struct CGNSWriteDialog
     inputInterfaceName->value(CTX.mesh.cgns_options.interfaceName.c_str());
     inputPatchName->value(CTX.mesh.cgns_options.patchName.c_str());
     checkButtonWriteBC->value(CTX.mesh.cgns_options.writeBC);
-    checkButtonWriteNormals->value(CTX.mesh.cgns_options.writeNormals);
+    checkButtonWriteNormals->value(CTX.mesh.cgns_options.normalSource);
     choiceVecDim->value(CTX.mesh.cgns_options.vectorDim - 2);
     checkButtonUnknownUserDef->value(CTX.mesh.cgns_options.writeUserDef);
 
@@ -1202,14 +1202,15 @@ struct CGNSWriteDialog
     cgnsw_gc_location_cb
       ((CTX.mesh.cgns_options.gridConnectivityLocation) ?
        roundButton1GCatFace : roundButton0GCatVertex, this);
-    cgnsw_write_dummy_bc_cb(checkButtonWriteBC, this);
+    // The order of the next 4 is important
+    cgnsw_normal_source_cb
+      ((CTX.mesh.cgns_options.normalSource == 2) ?
+       roundButton1NormalElem : roundButton0NormalGeo, this);
+    cgnsw_write_normals_cb(checkButtonWriteNormals, this);
     cgnsw_bc_location_cb
       ((CTX.mesh.cgns_options.bocoLocation) ?
        roundButton1BCatFace : roundButton0BCatVertex, this);
-    cgnsw_write_normals_cb(checkButtonWriteNormals, this);
-    cgnsw_normal_source_cb
-      ((CTX.mesh.cgns_options.normalSource) ?
-       roundButton1NormalElem : roundButton0NormalGeo, this);
+    cgnsw_write_dummy_bc_cb(checkButtonWriteBC, this);
   }
 };
 
diff --git a/Geo/CGNSOptions.h b/Geo/CGNSOptions.h
index 549d7bdf256b9cfa5d8764a8cc27011b3701d4aa..bede76d071f163e512d3f6379e2a1d18080bb85c 100644
--- a/Geo/CGNSOptions.h
+++ b/Geo/CGNSOptions.h
@@ -27,16 +27,15 @@ struct CGNSOptions
   int bocoLocation;                     // Location of BC (values
                                         // CGNSLocationType)
   int normalSource;                     // Source for BC normal data
-                                        // 0 - geometry
-                                        // 1 - elements
+                                        // 0 - do not write normals
+                                        // 1 - geometry
+                                        // 2 - elements
   int vectorDim;                        // Number of dimensions in a vector
                                         // (only relevant for a 2D mesh)
   bool writeBC;
-  bool writeNormals;
   bool writeUserDef;                    // T - write user-defined elements for
                                         //     element types unsupported by CGNS
 
-
 //--Constructor
 
   CGNSOptions()
@@ -49,15 +48,14 @@ struct CGNSOptions
   void setDefaults()
   {
     baseName = "Base_1";
-    zoneName = "Zone_&I&";
-    interfaceName = "Interface_&I&";
-    patchName = "Patch_&I&";
+    zoneName = "Zone_&I%4&";
+    interfaceName = "Interface_&I%4&";
+    patchName = "Patch_&I%4&";
     gridConnectivityLocation = 0;
     bocoLocation = 0;
-    normalSource = 0;
+    normalSource = 1;
     vectorDim = 2;
     writeBC = true;
-    writeNormals = true;
     writeUserDef = false;
   }
 };
diff --git a/Geo/CustomContainer.h b/Geo/CustomContainer.h
index e018f6c73b40ce0db06b0ebe46459559edcad32e..ef9cf2c5ffe4e850fea458964c2e3e1208123a2c 100644
--- a/Geo/CustomContainer.h
+++ b/Geo/CustomContainer.h
@@ -11,6 +11,8 @@
 #include <cstdlib>
 #include <cstring>
 
+#include "Message.h"
+
 /*******************************************************************************
  *
  * This lightweight pool class attempts to reduce memory usage by having the
@@ -104,6 +106,8 @@ class Pool
   void free_memory()
   {
     if(numUsedElement == 0) delete_all_blocks();
+    else Msg::Debug("Request to delete pool with used elements in "
+                    "CustomContainer.h");
   }
 
  private:
@@ -135,6 +139,7 @@ class Pool
       tailBlock = block->prev;
       delete block;
     }
+    tailElement = 0;
   }
 
   // Copy and assignment are not permitted
diff --git a/Geo/GModelIO_CGNS.cpp b/Geo/GModelIO_CGNS.cpp
index 550acb944d53669ac7a68336854518bce205a3b0..d620d0a00b372446ab4b14c5ac02ce210b3e49ee 100644
--- a/Geo/GModelIO_CGNS.cpp
+++ b/Geo/GModelIO_CGNS.cpp
@@ -39,7 +39,7 @@
 
 int cgnsErr(const int cgIndexFile = -1)
 {
-  Message::Error("This is a CGNS error\n");
+  Message::Error("Error detected by CGNS library\n");
   Message::Error(cg_get_error());
   if(cgIndexFile != -1)
     if(cg_close(cgIndexFile)) Message::Error("Unable to close CGNS file");
@@ -989,56 +989,61 @@ int write_CGNS_zones(GModel &model, const int zoneDefinition, const int numZone,
  * Write the remaining unconnected vertices as boundary conditions
  *--------------------------------------------------------------------*/
 
-      ZoneBoVec zoneBoVec;              // from 'MZoneBoundary.h'
-      if(mZoneBoundary.exteriorBoundaryVertices(zoneBoVec) == 0) {  // Have BC
-
-        // Sort by zone index and then entity index
-        const int numBoVert = zoneBoVec.size();
-        std::vector<int> iZBV(numBoVert);
-        for(int i = 0; i != numBoVert; ++i) iZBV[i] = i;
-        std::sort<int*, ZoneBoVecSort>(&iZBV[0], &iZBV[numBoVert],
-                                       ZoneBoVecSort(zoneBoVec));
-
-        dBuffer.reserve(1024);
-        iBuffer1.reserve(1024);
-
-        int iVert = 0;
-        while(iVert != numBoVert) {
-          dBuffer.clear();
-          iBuffer1.clear();
-          const int zoneIndex = zoneBoVec[iZBV[iVert]].zoneIndex;
-          const int patchIndex = zoneBoVec[iZBV[iVert]].bcPatchIndex;
-          const int iVertStart = iVert;
-          while(iVert != numBoVert &&
-                zoneBoVec[iZBV[iVert]].zoneIndex == zoneIndex &&
-                zoneBoVec[iZBV[iVert]].bcPatchIndex == patchIndex) {
-            VertexBoundary &vertBo = zoneBoVec[iZBV[iVert]];
-            dBuffer.push_back(vertBo.normal[0]);
-            dBuffer.push_back(vertBo.normal[1]);
-            if(vectorDim == 3) dBuffer.push_back(vertBo.normal[2]);
-            iBuffer1.push_back(vertBo.vertexIndex);
-            ++iVert;
-          }
-          const int numBCVert = iVert - iVertStart;
-
-          int cgIndexBoco;
-          std::string patchName = options.patchName;
-          expand_name(patchName, patchIndex-1, "Patch");
-          if(patchName.length() == 0) {
-            patchName = "Patch_";
-            patchName += CGNSNameStr(patchIndex).c_str();
-          }
-          if(cg_boco_write(cgIndexFile, cgIndexBase,
-                           zoneInfo[zoneIndex].cgIndex, patchName.c_str(),
-                           BCTypeNull, PointList, numBCVert, &iBuffer1[0],
-                           &cgIndexBoco))
-            return cgnsErr();
+      if(options.writeBC) {
+        ZoneBoVec zoneBoVec;            // from 'MZoneBoundary.h'
+        if(mZoneBoundary.exteriorBoundaryVertices
+           (options.normalSource, zoneBoVec) == 0) {
+
+          // Sort by zone index and then entity index
+          const int numBoVert = zoneBoVec.size();
+          std::vector<int> iZBV(numBoVert);
+          for(int i = 0; i != numBoVert; ++i) iZBV[i] = i;
+          std::sort<int*, ZoneBoVecSort>(&iZBV[0], &iZBV[numBoVert],
+                                         ZoneBoVecSort(zoneBoVec));
+
+          dBuffer.reserve(1024);
+          iBuffer1.reserve(1024);
+
+          int iVert = 0;
+          while(iVert != numBoVert) {
+            dBuffer.clear();
+            iBuffer1.clear();
+            const int zoneIndex = zoneBoVec[iZBV[iVert]].zoneIndex;
+            const int patchIndex = zoneBoVec[iZBV[iVert]].bcPatchIndex;
+            const int iVertStart = iVert;
+            while(iVert != numBoVert &&
+                  zoneBoVec[iZBV[iVert]].zoneIndex == zoneIndex &&
+                  zoneBoVec[iZBV[iVert]].bcPatchIndex == patchIndex) {
+              VertexBoundary &vertBo = zoneBoVec[iZBV[iVert]];
+              dBuffer.push_back(vertBo.normal[0]);
+              dBuffer.push_back(vertBo.normal[1]);
+              if(vectorDim == 3) dBuffer.push_back(vertBo.normal[2]);
+              iBuffer1.push_back(vertBo.vertexIndex);
+              ++iVert;
+            }
+            const int numBCVert = iVert - iVertStart;
+
+            int cgIndexBoco;
+            std::string patchName = options.patchName;
+            expand_name(patchName, patchIndex, "Patch");
+            if(patchName.length() == 0) {
+              patchName = "Patch_";
+              patchName += CGNSNameStr(patchIndex+1).c_str();
+            }
+            if(cg_boco_write(cgIndexFile, cgIndexBase,
+                             zoneInfo[zoneIndex].cgIndex, patchName.c_str(),
+                             BCTypeNull, PointList, numBCVert, &iBuffer1[0],
+                             &cgIndexBoco))
+              return cgnsErr();
 
-          int normalIndex;
-          if(cg_boco_normal_write(cgIndexFile, cgIndexBase,
-                                  zoneInfo[zoneIndex].cgIndex, cgIndexBoco,
-                                  &normalIndex, 1, RealDouble, &dBuffer[0]))
-            return cgnsErr();
+            if(options.normalSource) {
+              int normalIndex;
+              if(cg_boco_normal_write(cgIndexFile, cgIndexBase,
+                                      zoneInfo[zoneIndex].cgIndex, cgIndexBoco,
+                                      &normalIndex, 1, RealDouble, &dBuffer[0]))
+                return cgnsErr();
+            }
+          }
         }
       }
 
diff --git a/Geo/MZoneBoundary.cpp b/Geo/MZoneBoundary.cpp
index 07223536945934f4152bcf4c71ffaf7481dcaf87..4b4561f41fa3f5876a27cc35d23e665cd2ef7c43 100644
--- a/Geo/MZoneBoundary.cpp
+++ b/Geo/MZoneBoundary.cpp
@@ -13,6 +13,14 @@
 #include "MZoneBoundary.h"
 #include "Message.h"
 
+//--Types at local scope
+
+enum {
+  NormalSourceGeometry = 1,
+  NormalSourceElement = 2
+};
+
+
 /*******************************************************************************
  *
  * class BCPatchIndex
@@ -94,12 +102,12 @@ class BCPatchIndex
   // Once all entity tags have been added, generate patch indices
   void generatePatchIndices()
   {
-    if(sharedPatch) {  // Don't renumber if not shared to preserve entity
-                       // numbers
-      int c = 1;
+//     if(sharedPatch) {  // Don't renumber if not shared to preserve entity
+//                        // numbers.  Mostly useful for debugging.
+      int c = 0;
       for(PatchDataListIt pDIt = patchData.begin(); pDIt != patchData.end();
           ++pDIt) pDIt->index = c++;
-    }
+//     }
   }
 
   // Get the patch index for an entity (generate patch indices first)
@@ -135,6 +143,7 @@ class BCPatchIndex
  * I/O
  * ===
  *
+ *   normalSource       - (I) source for normals
  *   vertex             - (I) vertex to add to zoneBoVec and to find the normal
  *                            at
  *   zoneIndex          - (I)
@@ -159,7 +168,8 @@ class BCPatchIndex
 
 template <unsigned DIM, typename FaceT>
 void updateBoVec
-(const MVertex *const vertex, const int zoneIndex, const int vertIndex,
+(const int normalSource, const MVertex *const vertex, const int zoneIndex,
+ const int vertIndex,
  const CCon::FaceVector<typename MZoneBoundary<DIM>::template
  GlobalVertexData<FaceT>::FaceDataB> &faces, ZoneBoVec &zoneBoVec,
  BCPatchIndex &patch, bool &warnNormFromElem);
@@ -248,12 +258,15 @@ int edge_normal
 
 template <>
 void updateBoVec<2, MEdge>
-(const MVertex *const vertex, const int zoneIndex, const int vertIndex,
+(const int normalSource, const MVertex *const vertex, const int zoneIndex,
+ const int vertIndex,
  const CCon::FaceVector<MZoneBoundary<2>::GlobalVertexData<MEdge>::FaceDataB>
  &faces, ZoneBoVec &zoneBoVec, BCPatchIndex &patch, bool &warnNormFromElem)
 {
 
-  GEntity *const ent = vertex->onWhat();
+  GEntity *ent;
+  if(normalSource == NormalSourceElement) goto getNormalFromElements;
+  ent = vertex->onWhat();
   if(ent == 0) {
     goto getNormalFromElements;
     // No entity: try to find a normal from the faces
@@ -406,7 +419,7 @@ void updateBoVec<2, MEdge>
  *--------------------------------------------------------------------*/
 
   {
-    if(warnNormFromElem) {
+    if(warnNormFromElem && normalSource == 1) {
       Msg::Warning("Some or all boundary normals were determined from mesh "
                    "elements instead of from the geometry");
       warnNormFromElem = false;
@@ -465,18 +478,22 @@ void updateBoVec<2, MEdge>
  *     splitting of the patch.  Merging of entities into a single patch has not
  *     yet been implemented, but the best way to accomplish this is probably to
  *     split all entities into separate patches and then record which pairs of
- *     entities cannot be merged.  Later, all others and be merged and
+ *     entities cannot be merged.  Later, all others can be merged and
  *     'zoneBoVec' updated accordingly.
  *
  ******************************************************************************/
 
 template <>
 void updateBoVec<3, MFace>
-(const MVertex *const vertex, const int zoneIndex, const int vertIndex,
+(const int normalSource, const MVertex *const vertex, const int zoneIndex,
+ const int vertIndex,
  const CCon::FaceVector<MZoneBoundary<3>::GlobalVertexData<MFace>::FaceDataB>
  &faces, ZoneBoVec &zoneBoVec, BCPatchIndex &patch, bool &warnNormFromElem)
 {
-  GEntity *const ent = vertex->onWhat();
+
+  GEntity *ent;
+  if(normalSource == NormalSourceElement) goto getNormalFromElements;
+  ent = vertex->onWhat();
   if(ent == 0) {
     goto getNormalFromElements;
     // No entity: try to find a normal from the faces
@@ -684,7 +701,7 @@ void updateBoVec<3, MFace>
  *--------------------------------------------------------------------*/
 
   {
-    if(warnNormFromElem) {
+    if(warnNormFromElem && normalSource == 1) {
       Msg::Warning("Some or all boundary normals were determined from mesh "
                    "elements instead of from the geometry");
       warnNormFromElem = false;
@@ -890,19 +907,27 @@ int MZoneBoundary<DIM>::interiorBoundaryVertices
  * I/O
  * ===
  *
- *   zoneBoMap          - (O) Boundary vertices for all zones.
+ *   normalSource       - (I) Source for normals if the are requested.
+ *                            0 - do not obtain normals
+ *                            1 - geometry
+ *                            2 - elements
+ *                            Note that if normals cannot be found from the
+ *                            geometry, the elements will be used.
+ *   zoneBoMap          - (O) boundary vertices for all zones.
  *
  * Notes
  * =====
  *
  *   - Should only be called after all zones have been added via
  *     'interior_boundary_vertices'
+ *   - Unless normals are obtained from the geometry, all boundary vertices will
+ *     be written to one patch
  *
  ******************************************************************************/
 
 template <unsigned DIM>
 int MZoneBoundary<DIM>::exteriorBoundaryVertices
-(ZoneBoVec &zoneBoVec)
+(const int normalSource, ZoneBoVec &zoneBoVec)
 {
 
   if(globalBoVertMap.size() == 0) return 1;
@@ -928,19 +953,31 @@ int MZoneBoundary<DIM>::exteriorBoundaryVertices
 
 //--Try to find an outwards normal for this vertex
 
-      updateBoVec<DIM, FaceT>(vMapIt->first, zoneData.zoneIndex,
-                              zoneData.vertexIndex, vMapIt->second.faces,
-                              zoneBoVec, patch, warnNormFromElem);
+      if(normalSource) {
+        updateBoVec<DIM, FaceT>(normalSource, vMapIt->first, zoneData.zoneIndex,
+                                zoneData.vertexIndex, vMapIt->second.faces,
+                                zoneBoVec, patch, warnNormFromElem);
+      }
+      else {
+        // Keys to 'globalBoVertMap' will not change so const_cast is okay.
+        zoneBoVec.push_back(VertexBoundary(zoneData.zoneIndex, 0, SVector3(0.),
+                                           const_cast<MVertex*>(vMapIt->first),
+                                           zoneData.vertexIndex));
+      }
     }
   }
 
-  // zoneBoVec has entities as a patch index.  Update with the actual patch
-  // index in 'patch'.
-  patch.generatePatchIndices();
-  const int nBoVert = zoneBoVec.size();
-  for(int iBoVert = 0; iBoVert != nBoVert; ++iBoVert) {
-    zoneBoVec[iBoVert].bcPatchIndex =
-      patch.getIndex(zoneBoVec[iBoVert].bcPatchIndex);
+  // If normals were written from the geometry, zoneBoVec currently stores
+  // entities in bcPatchIndex.  Update with the actual patch index.  Why? -
+  // because two entities may have been merged if the interface between them is
+  // continuous.
+  if(normalSource == NormalSourceGeometry) {
+    patch.generatePatchIndices();
+    const int nBoVert = zoneBoVec.size();
+    for(int iBoVert = 0; iBoVert != nBoVert; ++iBoVert) {
+      zoneBoVec[iBoVert].bcPatchIndex =
+        patch.getIndex(zoneBoVec[iBoVert].bcPatchIndex);
+    }
   }
 
   return 0;
diff --git a/Geo/MZoneBoundary.h b/Geo/MZoneBoundary.h
index b0e8a8e4d0edb75f1cdf18b5d00251055680d3e1..8689e3cc1382355cd377201ae2c9195ef570eaf5 100644
--- a/Geo/MZoneBoundary.h
+++ b/Geo/MZoneBoundary.h
@@ -311,7 +311,7 @@ class MZoneBoundary
 //--Return exterior boundary vertices (unconnected vertices at the extent of the
 //--domain
 
-  int exteriorBoundaryVertices(ZoneBoVec &zoneBoVec);
+  int exteriorBoundaryVertices(const int normalSource, ZoneBoVec &zoneBoVec);
 
 //--Memory management