diff --git a/components/omega/doc/devGuide/VertCoord.md b/components/omega/doc/devGuide/VertCoord.md
index d6e2edd2016e..e5a60115a0f0 100644
--- a/components/omega/doc/devGuide/VertCoord.md
+++ b/components/omega/doc/devGuide/VertCoord.md
@@ -4,11 +4,23 @@
 
 ### Initialization
 
-The default `VertCoord` instance is  created by the `init` method, which assumes that `Decomp` has already been initialized.
+The default `VertCoord` instance is created by the `init` method, which assumes that `Decomp` has already been initialized.
 ```
 Decomp::init();
 VertCoord::init();
 ```
+The `init` method accepts two optional arguments with default values:
+```
+VertCoord::init(const bool ReadStream = true, const int NVertLayers = 0)
+```
+These arguments provide flexibility, particularly for unit testing. When `init(false)` is used, the method skips reading
+the `InitialVertCoord` stream. This is useful for unit tests that rely on meshes lacking the fields required by that stream.
+In this case, the min/max layer arrays are initialized with default values based on the number of vertical layers, and
+`bottomDepth` remains uninitialized. Some tests may utilize a specific number of vertical layers that differs from what is
+defined in the mesh file. To handle this, the `VertCoord` can be initialized with `init(false, LocNVertLayers)` to explicitly
+set the number of vertical layers. An argument for `ReadStream` must be provided in order to explicitly set `NVertLayers`.
+For example, `init(42)` is invalid, `init(false, 42)` must be called instead.
+
 The default instance can be retrieved by:
 ```
 auto *DefVertCoord = VertCoord::getDefault();
@@ -16,9 +28,12 @@ auto *DefVertCoord = VertCoord::getDefault();
 
 Additional instances can be created by calling the `create` method.
 ```
-VertCoord *VertCoord::create(const std::string &Name,  ///< [in] Name for vertical coordinate
-                             const Decomp *MeshDecomp, ///< [in] Decomp for mesh
-                             Config *Options)          ///< [in] Confiuration options
+VertCoord *VertCoord::create(const std::string &Name,      ///< [in] Name for vertical coordinate
+                             const Decomp *MeshDecomp,     ///< [in] Decomp for mesh
+                             Config *Options,              ///< [in] Confiuration options
+                             const bool ReadStream = true, ///< [in] optional logical to read stream
+                             const int NVertLayers = 0     ///< [in] optional int to set vertical dim
+)
 ```
 This calls the constructor and places the instance in a map that can be used to retrieve the instance by name:
 ```
diff --git a/components/omega/src/infra/IOStream.cpp b/components/omega/src/infra/IOStream.cpp
index 6d3c3dfb0c8a..802e82f6e928 100644
--- a/components/omega/src/infra/IOStream.cpp
+++ b/components/omega/src/infra/IOStream.cpp
@@ -1578,11 +1578,19 @@ Error IOStream::readFieldData(
    // lower case
    std::string OldFieldName = FieldName;
    OldFieldName[0]          = std::tolower(OldFieldName[0]);
-   bool OnHost              = FieldPtr->isOnHost();
-   bool IsDistributed       = FieldPtr->isDistributed();
-   bool IsTimeDependent     = FieldPtr->isTimeDependent();
-   ArrayDataType MyType     = FieldPtr->getType();
-   int NDims                = FieldPtr->getNumDims();
+   // Check for "Layer" in field name, backwards compatibility requires
+   // replacing "Layer" with "Level"
+   std::string OmegaSubStr = "Layer";
+   std::string MPASSubStr  = "Level";
+   size_t pos              = OldFieldName.find(OmegaSubStr);
+   if (pos != std::string::npos) {
+      OldFieldName.replace(pos, OmegaSubStr.length(), MPASSubStr);
+   }
+   bool OnHost          = FieldPtr->isOnHost();
+   bool IsDistributed   = FieldPtr->isDistributed();
+   bool IsTimeDependent = FieldPtr->isTimeDependent();
+   ArrayDataType MyType = FieldPtr->getType();
+   int NDims            = FieldPtr->getNumDims();
    if (NDims < 0)
       ABORT_ERROR("IOStream readFieldData: "
                   "Invalid number of dimensions for Field {}",
diff --git a/components/omega/src/ocn/AuxiliaryState.cpp b/components/omega/src/ocn/AuxiliaryState.cpp
index 64cdad749f25..2164fad2d6a7 100644
--- a/components/omega/src/ocn/AuxiliaryState.cpp
+++ b/components/omega/src/ocn/AuxiliaryState.cpp
@@ -18,14 +18,15 @@ static std::string stripDefault(const std::string &Name) {
 // Constructor. Constructs the member auxiliary variables and registers their
 // fields with IOStreams
 AuxiliaryState::AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
-                               Halo *MeshHalo, int NVertLayers, int NTracers)
-    : Mesh(Mesh), MeshHalo(MeshHalo), Name(stripDefault(Name)),
+                               const VertCoord *VCoord, Halo *MeshHalo,
+                               int NVertLayers, int NTracers)
+    : Mesh(Mesh), VCoord(VCoord), MeshHalo(MeshHalo), Name(stripDefault(Name)),
       KineticAux(stripDefault(Name), Mesh, NVertLayers),
       LayerThicknessAux(stripDefault(Name), Mesh, NVertLayers),
       VorticityAux(stripDefault(Name), Mesh, NVertLayers),
-      VelocityDel2Aux(stripDefault(Name), Mesh, NVertLayers),
+      VelocityDel2Aux(stripDefault(Name), Mesh, VCoord, NVertLayers),
       WindForcingAux(stripDefault(Name), Mesh),
-      TracerAux(stripDefault(Name), Mesh, NVertLayers, NTracers) {
+      TracerAux(stripDefault(Name), Mesh, VCoord, NVertLayers, NTracers) {
 
    GroupName = "AuxiliaryState";
    if (Name != "Default") {
@@ -192,7 +193,8 @@ void AuxiliaryState::computeAll(const OceanState *State,
 
 // Create a non-default auxiliary state
 AuxiliaryState *AuxiliaryState::create(const std::string &Name,
-                                       const HorzMesh *Mesh, Halo *MeshHalo,
+                                       const HorzMesh *Mesh,
+                                       const VertCoord *VCoord, Halo *MeshHalo,
                                        int NVertLayers, const int NTracers) {
    if (AllAuxStates.find(Name) != AllAuxStates.end()) {
       LOG_ERROR("Attempted to create a new AuxiliaryState with name {} but it "
@@ -202,7 +204,7 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name,
    }
 
    auto *NewAuxState =
-       new AuxiliaryState(Name, Mesh, MeshHalo, NVertLayers, NTracers);
+       new AuxiliaryState(Name, Mesh, VCoord, MeshHalo, NVertLayers, NTracers);
    AllAuxStates.emplace(Name, NewAuxState);
 
    return NewAuxState;
@@ -211,14 +213,15 @@ AuxiliaryState *AuxiliaryState::create(const std::string &Name,
 // Create the default auxiliary state. Assumes that HorzMesh, VertCoord and
 // Halo have been initialized.
 void AuxiliaryState::init() {
-   const HorzMesh *DefMesh = HorzMesh::getDefault();
-   Halo *DefHalo           = Halo::getDefault();
+   const HorzMesh *DefMesh    = HorzMesh::getDefault();
+   const VertCoord *DefVCoord = VertCoord::getDefault();
+   Halo *DefHalo              = Halo::getDefault();
 
    int NVertLayers = VertCoord::getDefault()->NVertLayers;
    int NTracers    = Tracers::getNumTracers();
 
    AuxiliaryState::DefaultAuxState = AuxiliaryState::create(
-       "Default", DefMesh, DefHalo, NVertLayers, NTracers);
+       "Default", DefMesh, DefVCoord, DefHalo, NVertLayers, NTracers);
 
    Config *OmegaConfig = Config::getOmegaConfig();
    DefaultAuxState->readConfigOptions(OmegaConfig);
diff --git a/components/omega/src/ocn/AuxiliaryState.h b/components/omega/src/ocn/AuxiliaryState.h
index 92eee116e6ca..509aee6dcf22 100644
--- a/components/omega/src/ocn/AuxiliaryState.h
+++ b/components/omega/src/ocn/AuxiliaryState.h
@@ -50,7 +50,8 @@ class AuxiliaryState {
 
    // Create a non-default auxiliary state
    static AuxiliaryState *create(const std::string &Name, const HorzMesh *Mesh,
-                                 Halo *MeshHalo, int NVertLayers, int NTracers);
+                                 const VertCoord *VCoord, Halo *MeshHalo,
+                                 int NVertLayers, int NTracers);
 
    /// Get the default auxiliary state
    static AuxiliaryState *getDefault();
@@ -82,13 +83,15 @@ class AuxiliaryState {
                    int TimeLevel) const;
 
  private:
-   AuxiliaryState(const std::string &Name, const HorzMesh *Mesh, Halo *MeshHalo,
-                  int NVertLayers, int NTracers);
+   AuxiliaryState(const std::string &Name, const HorzMesh *Mesh,
+                  const VertCoord *VCoord, Halo *MeshHalo, int NVertLayers,
+                  int NTracers);
 
    AuxiliaryState(const AuxiliaryState &) = delete;
    AuxiliaryState(AuxiliaryState &&)      = delete;
 
    const HorzMesh *Mesh;
+   const VertCoord *VCoord;
    Halo *MeshHalo;
    static AuxiliaryState *DefaultAuxState;
    static std::map<std::string, std::unique_ptr<AuxiliaryState>> AllAuxStates;
diff --git a/components/omega/src/ocn/HorzMesh.cpp b/components/omega/src/ocn/HorzMesh.cpp
index f08b18458e6d..0c626b49905e 100644
--- a/components/omega/src/ocn/HorzMesh.cpp
+++ b/components/omega/src/ocn/HorzMesh.cpp
@@ -37,18 +37,15 @@ void HorzMesh::init() {
    // Retrieve the default decomposition
    Decomp *DefDecomp = Decomp::getDefault();
 
-   I4 NVertLayers = VertCoord::getDefault()->NVertLayers;
-
    // Create the default mesh and set pointer to it
-   HorzMesh::DefaultHorzMesh = create("Default", DefDecomp, NVertLayers);
+   HorzMesh::DefaultHorzMesh = create("Default", DefDecomp);
 }
 
 //------------------------------------------------------------------------------
 // Construct a new local mesh given a decomposition
 
 HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
-                   Decomp *MeshDecomp,      //< [in] Decomp for the new mesh
-                   I4 InNVertLayers         //< [in} num vertical layers
+                   Decomp *MeshDecomp       //< [in] Decomp for the new mesh
 ) {
 
    MeshName = Name;
@@ -56,9 +53,6 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
    // Retrieve mesh files name from Decomp
    MeshFileName = MeshDecomp->MeshFileName;
 
-   // Set NVertLayers
-   NVertLayers = InNVertLayers;
-
    // Retrieve mesh cell/edge/vertex totals from Decomp
    NCellsHalo  = MeshDecomp->NCellsHalo;
    NCellsHaloH = MeshDecomp->NCellsHaloH;
@@ -141,9 +135,6 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
    // Compute EdgeSignOnCells and EdgeSignOnVertex
    computeEdgeSign();
 
-   // TODO: implement setMasks during Mesh constructor
-   setMasks(NVertLayers);
-
    // set mesh scaling coefficients
    setMeshScaling();
 
@@ -152,8 +143,7 @@ HorzMesh::HorzMesh(const std::string &Name, //< [in] Name for new mesh
 /// Creates a new mesh by calling the constructor and puts it in the
 /// AllHorzMeshes map
 HorzMesh *HorzMesh::create(const std::string &Name, //< [in] Name for new mesh
-                           Decomp *MeshDecomp, //< [in] Decomp for the new mesh
-                           I4 InNVertLayers    //< [in] num vertical layers
+                           Decomp *MeshDecomp //< [in] Decomp for the new mesh
 ) {
    // Check to see if a mesh of the same name already exists and
    // if so, exit with an error
@@ -166,7 +156,7 @@ HorzMesh *HorzMesh::create(const std::string &Name, //< [in] Name for new mesh
 
    // create a new mesh on the heap and put it in a map of
    // unique_ptrs, which will manage its lifetime
-   auto *NewHorzMesh = new HorzMesh(Name, MeshDecomp, InNVertLayers);
+   auto *NewHorzMesh = new HorzMesh(Name, MeshDecomp);
    AllHorzMeshes.emplace(Name, NewHorzMesh);
 
    return NewHorzMesh;
@@ -574,33 +564,6 @@ void HorzMesh::computeEdgeSign() {
    EdgeSignOnVertexH = createHostMirrorCopy(EdgeSignOnVertex);
 } // end computeEdgeSign
 
-//------------------------------------------------------------------------------
-// set computational masks for mesh elements
-// TODO: this is just a placeholder, implement actual masks for edges, cells,
-// and vertices
-void HorzMesh::setMasks(int NVertLayers) {
-
-   EdgeMask = Array2DReal("EdgeMask", NEdgesSize, NVertLayers);
-
-   OMEGA_SCOPE(LocEdgeMask, EdgeMask);
-   OMEGA_SCOPE(LocCellsOnEdge, CellsOnEdge);
-   OMEGA_SCOPE(LocNCellsAll, NCellsAll);
-
-   deepCopy(EdgeMask, 1.0);
-   parallelFor(
-       {NEdgesAll, NVertLayers}, KOKKOS_LAMBDA(int Edge, int K) {
-          const I4 Cell1 = LocCellsOnEdge(Edge, 0);
-          const I4 Cell2 = LocCellsOnEdge(Edge, 1);
-          if (!(Cell1 >= 0 and Cell1 < LocNCellsAll) or
-              !(Cell2 >= 0 and Cell2 < LocNCellsAll)) {
-             LocEdgeMask(Edge, K) = 0.0;
-          }
-       });
-
-   EdgeMaskH = createHostMirrorCopy(EdgeMask);
-
-} // end setMasks
-
 //------------------------------------------------------------------------------
 // Set mesh scaling coefficients for mixing terms in momentum and tracer
 // equations so viscosity and diffusion scale with mesh.
diff --git a/components/omega/src/ocn/HorzMesh.h b/components/omega/src/ocn/HorzMesh.h
index a1a73a0db83a..e6c3d375079f 100644
--- a/components/omega/src/ocn/HorzMesh.h
+++ b/components/omega/src/ocn/HorzMesh.h
@@ -14,7 +14,6 @@
 #include "Decomp.h"
 #include "MachEnv.h"
 #include "OmegaKokkos.h"
-#include "VertCoord.h"
 
 #include <memory>
 #include <string>
@@ -68,8 +67,7 @@ class HorzMesh {
 
    /// Construct a new local mesh for a given decomposition
    HorzMesh(const std::string &Name, ///< [in] Name for mesh
-            Decomp *Decomp,          ///< [in] Decomposition for mesh
-            I4 InNVertLayers         ///< [in] num vertical layers
+            Decomp *Decomp           ///< [in] Decomposition for mesh
    );
 
    // Forbid copy and move construction
@@ -81,8 +79,6 @@ class HorzMesh {
    // private function.
    void computeEdgeSign();
 
-   void setMasks(int NVertLayers);
-
    void setMeshScaling();
 
    // Variables
@@ -97,8 +93,6 @@ class HorzMesh {
    // Note that all sizes are actual counts (1-based) so that loop extents
    // should always use the 0:NCellsXX-1 form.
 
-   I4 NVertLayers; ///< number of vertical layers
-
    Array1DI4 NCellsHalo;      ///< num cells owned+halo for halo layer
    HostArray1DI4 NCellsHaloH; ///< num cells owned+halo for halo layer
    I4 NCellsOwned;            ///< Number of cells owned by this task
@@ -272,8 +266,7 @@ class HorzMesh {
    /// Creates a new mesh by calling the constructor and puts it in the
    /// AllHorzMeshes map
    static HorzMesh *create(const std::string &Name, ///< [in] Name for mesh
-                           Decomp *Decomp,  ///< [in] Decomposition for mesh
-                           I4 InNVertLayers ///< [in] num vertival layers
+                           Decomp *Decomp ///< [in] Decomposition for mesh
    );
 
    /// Destructor - deallocates all memory and deletes a HorzMesh
diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp
index c07a3d4aeba1..2f4c725d3aeb 100644
--- a/components/omega/src/ocn/OceanInit.cpp
+++ b/components/omega/src/ocn/OceanInit.cpp
@@ -128,9 +128,8 @@ int initOmegaModules(MPI_Comm Comm) {
       ABORT_ERROR("ocnInit: Error initializing default halo");
    }
 
-   VertCoord::init1();
    HorzMesh::init();
-   VertCoord::init2();
+   VertCoord::init();
    Tracers::init();
    AuxiliaryState::init();
    Tendencies::init();
diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp
index 7d3efab243bd..3a256402d007 100644
--- a/components/omega/src/ocn/Tendencies.cpp
+++ b/components/omega/src/ocn/Tendencies.cpp
@@ -221,11 +221,12 @@ Tendencies::Tendencies(const std::string &Name, ///< [in] Name for tendencies
                        Config *Options,         ///< [in] Configuration options
                        CustomTendencyType InCustomThicknessTend,
                        CustomTendencyType InCustomVelocityTend)
-    : ThicknessFluxDiv(Mesh), PotientialVortHAdv(Mesh), KEGrad(Mesh),
-      SSHGrad(Mesh), VelocityDiffusion(Mesh), VelocityHyperDiff(Mesh),
-      WindForcing(Mesh), BottomDrag(Mesh, VCoord), TracerHorzAdv(Mesh),
-      TracerDiffusion(Mesh), TracerHyperDiff(Mesh),
-      CustomThicknessTend(InCustomThicknessTend),
+    : ThicknessFluxDiv(Mesh), PotientialVortHAdv(Mesh, VCoord),
+      KEGrad(Mesh, VCoord), SSHGrad(Mesh, VCoord),
+      VelocityDiffusion(Mesh, VCoord), VelocityHyperDiff(Mesh, VCoord),
+      WindForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord),
+      TracerHorzAdv(Mesh, VCoord), TracerDiffusion(Mesh, VCoord),
+      TracerHyperDiff(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend),
       CustomVelocityTend(InCustomVelocityTend) {
 
    // Tendency arrays
diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp
index aa4570d1f297..53c79c1c063a 100644
--- a/components/omega/src/ocn/TendencyTerms.cpp
+++ b/components/omega/src/ocn/TendencyTerms.cpp
@@ -23,53 +23,60 @@ ThicknessFluxDivOnCell::ThicknessFluxDivOnCell(const HorzMesh *Mesh)
       DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell),
       EdgeSignOnCell(Mesh->EdgeSignOnCell) {}
 
-PotentialVortHAdvOnEdge::PotentialVortHAdvOnEdge(const HorzMesh *Mesh)
+PotentialVortHAdvOnEdge::PotentialVortHAdvOnEdge(const HorzMesh *Mesh,
+                                                 const VertCoord *VCoord)
     : NEdgesOnEdge(Mesh->NEdgesOnEdge), EdgesOnEdge(Mesh->EdgesOnEdge),
-      WeightsOnEdge(Mesh->WeightsOnEdge), EdgeMask(Mesh->EdgeMask) {}
+      WeightsOnEdge(Mesh->WeightsOnEdge), EdgeMask(VCoord->EdgeMask) {}
 
-KEGradOnEdge::KEGradOnEdge(const HorzMesh *Mesh)
+KEGradOnEdge::KEGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord)
     : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge),
-      EdgeMask(Mesh->EdgeMask) {}
+      EdgeMask(VCoord->EdgeMask) {}
 
-SSHGradOnEdge::SSHGradOnEdge(const HorzMesh *Mesh)
+SSHGradOnEdge::SSHGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord)
     : CellsOnEdge(Mesh->CellsOnEdge), DcEdge(Mesh->DcEdge),
-      EdgeMask(Mesh->EdgeMask) {}
+      EdgeMask(VCoord->EdgeMask) {}
 
-VelocityDiffusionOnEdge::VelocityDiffusionOnEdge(const HorzMesh *Mesh)
+VelocityDiffusionOnEdge::VelocityDiffusionOnEdge(const HorzMesh *Mesh,
+                                                 const VertCoord *VCoord)
     : CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge),
       DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge),
-      MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(Mesh->EdgeMask) {}
+      MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(VCoord->EdgeMask) {}
 
-VelocityHyperDiffOnEdge::VelocityHyperDiffOnEdge(const HorzMesh *Mesh)
+VelocityHyperDiffOnEdge::VelocityHyperDiffOnEdge(const HorzMesh *Mesh,
+                                                 const VertCoord *VCoord)
     : CellsOnEdge(Mesh->CellsOnEdge), VerticesOnEdge(Mesh->VerticesOnEdge),
       DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge),
-      MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {}
+      MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(VCoord->EdgeMask) {}
 
-WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh)
-    : Enabled(false), EdgeMask(Mesh->EdgeMask) {}
+WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh,
+                                     const VertCoord *VCoord)
+    : Enabled(false), EdgeMask(VCoord->EdgeMask) {}
 
 BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh,
                                    const VertCoord *VCoord)
     : Enabled(false), Coeff(0), CellsOnEdge(Mesh->CellsOnEdge),
-      NVertLayers(VCoord->NVertLayers), EdgeMask(Mesh->EdgeMask) {}
+      NVertLayers(VCoord->NVertLayers), EdgeMask(VCoord->EdgeMask) {}
 
-TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh)
+TracerHorzAdvOnCell::TracerHorzAdvOnCell(const HorzMesh *Mesh,
+                                         const VertCoord *VCoord)
     : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
-      DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell), EdgeMask(Mesh->EdgeMask) {
-}
+      DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell),
+      EdgeMask(VCoord->EdgeMask) {}
 
-TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh)
+TracerDiffOnCell::TracerDiffOnCell(const HorzMesh *Mesh,
+                                   const VertCoord *VCoord)
     : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
       DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AreaCell(Mesh->AreaCell),
-      MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(Mesh->EdgeMask) {}
+      MeshScalingDel2(Mesh->MeshScalingDel2), EdgeMask(VCoord->EdgeMask) {}
 
-TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh)
+TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh,
+                                             const VertCoord *VCoord)
     : NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
       DvEdge(Mesh->DvEdge), DcEdge(Mesh->DcEdge), AreaCell(Mesh->AreaCell),
-      MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {}
+      MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(VCoord->EdgeMask) {}
 
 } // end namespace OMEGA
 
diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h
index a6a07574a931..fefebbd4f995 100644
--- a/components/omega/src/ocn/TendencyTerms.h
+++ b/components/omega/src/ocn/TendencyTerms.h
@@ -72,7 +72,7 @@ class PotentialVortHAdvOnEdge {
    bool Enabled;
 
    /// constructor declaration
-   PotentialVortHAdvOnEdge(const HorzMesh *Mesh);
+   PotentialVortHAdvOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    /// The functor takes edge index, vertical chunk index, and arrays for
    /// normalized relative vorticity, normalized planetary vorticity, layer
@@ -120,7 +120,7 @@ class KEGradOnEdge {
    bool Enabled;
 
    /// constructor declaration
-   KEGradOnEdge(const HorzMesh *Mesh);
+   KEGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    /// The functor takes edge index, vertical chunk index, and kinetic energy
    /// array as inputs, outputs the tendency array
@@ -152,7 +152,7 @@ class SSHGradOnEdge {
    bool Enabled;
 
    /// constructor declaration
-   SSHGradOnEdge(const HorzMesh *Mesh);
+   SSHGradOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    /// The functor takes edge index, vertical chunk index, and array of
    /// layer thickness/SSH, outputs tendency array
@@ -187,7 +187,7 @@ class VelocityDiffusionOnEdge {
    Real ViscDel2;
 
    /// constructor declaration
-   VelocityDiffusionOnEdge(const HorzMesh *Mesh);
+   VelocityDiffusionOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    /// The functor takes edge index, vertical chunk index, and arrays for
    /// divergence of horizontal velocity (defined at cell centers) and relative
@@ -236,7 +236,7 @@ class VelocityHyperDiffOnEdge {
    Real DivFactor;
 
    /// Constructor declaration
-   VelocityHyperDiffOnEdge(const HorzMesh *Mesh);
+   VelocityHyperDiffOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    /// The functor takes the edge index, vertical chunk index, and arrays for
    /// the laplacian of divergence of horizontal velocity and the laplacian of
@@ -284,7 +284,7 @@ class WindForcingOnEdge {
    Real SaltWaterDensity;
 
    /// constructor declaration
-   WindForcingOnEdge(const HorzMesh *Mesh);
+   WindForcingOnEdge(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    /// The functor takes the edge index, vertical chunk index, and arrays for
    /// normal wind stress and edge layer thickness, outputs tendency array
@@ -344,7 +344,7 @@ class TracerHorzAdvOnCell {
  public:
    bool Enabled;
 
-   TracerHorzAdvOnCell(const HorzMesh *Mesh);
+   TracerHorzAdvOnCell(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell,
                                    I4 KChunk, const Array2DReal &NormVelEdge,
@@ -389,7 +389,7 @@ class TracerDiffOnCell {
 
    Real EddyDiff2;
 
-   TracerDiffOnCell(const HorzMesh *Mesh);
+   TracerDiffOnCell(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    KOKKOS_FUNCTION void
    operator()(const Array3DReal &Tend, I4 L, I4 ICell, I4 KChunk,
@@ -444,7 +444,7 @@ class TracerHyperDiffOnCell {
 
    Real EddyDiff4;
 
-   TracerHyperDiffOnCell(const HorzMesh *Mesh);
+   TracerHyperDiffOnCell(const HorzMesh *Mesh, const VertCoord *VCoord);
 
    KOKKOS_FUNCTION void operator()(const Array3DReal &Tend, I4 L, I4 ICell,
                                    I4 KChunk,
diff --git a/components/omega/src/ocn/VertCoord.cpp b/components/omega/src/ocn/VertCoord.cpp
index 7f3133aa22fd..9e3f8f17d70c 100644
--- a/components/omega/src/ocn/VertCoord.cpp
+++ b/components/omega/src/ocn/VertCoord.cpp
@@ -23,33 +23,64 @@ VertCoord *VertCoord::DefaultVertCoord = nullptr;
 std::map<std::string, std::unique_ptr<VertCoord>> VertCoord::AllVertCoords;
 
 //------------------------------------------------------------------------------
-// Begin initialization of default vertical coordinate, requires prior
-// initialization of Decomp.
-void VertCoord::init1() {
+// Initialize the default vertical coordinate, requires prior initialization
+// of Decomp. The optional arguments simplify the use of the VertCoord in some
+// unit tests. If ReadStream is false, the InitialVertCoord stream will not be
+// read during construction. If a value for NVertLayers is passed as argument,
+// the dimension will not be read from the mesh file.
+void VertCoord::init(
+    const bool
+        ReadStream, //< [in] optional argument to read stream, true by default
+    const int
+        InNVertLayers //< [in] optional argument to set NVertLayers explicitly
+                      //< instead of reading dimension from mesh file
+) {
 
    Decomp *DefDecomp = Decomp::getDefault();
 
-   VertCoord::DefaultVertCoord = create("Default", DefDecomp);
-
-} // end init1
-
-//------------------------------------------------------------------------------
-// Complete initialization of default vertical coordinate, requires prior
-// initialization of HorzMesh.
-void VertCoord::init2() {
-
    Config *OmegaConfig = Config::getOmegaConfig();
 
-   DefaultVertCoord->completeSetup(OmegaConfig);
+   VertCoord::DefaultVertCoord =
+       create("Default", DefDecomp, OmegaConfig, ReadStream, InNVertLayers);
 
-} // end init2
+} // end init
 
 //------------------------------------------------------------------------------
-// Construct a new VertCoord instance given a Decomp. New object is incomplete
-// and completeSetup must be called afterwards
+// Construct a new VertCoord instance given a Decomp.
 VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
-                     const Decomp *Decomp      //< [in] associated Decomp
+                     const Decomp *Decomp,     //< [in] associated Decomp
+                     Config *Options,          //< [in] configuration options
+                     const bool ReadStream,    //< [in] logical to read stream
+                     const int InNVertLayers   //< [in] int to set vertical dim
 ) {
+   Error Err; // Error code
+
+   // Read Config for movement weight type, store in enum
+   Config VCoordConfig("VertCoord");
+   Err += Options->get(VCoordConfig);
+   CHECK_ERROR_ABORT(Err, "VertCoord: VertCoord group not found in Config");
+
+   std::string MovementWeightStr;
+   Err += VCoordConfig.get("MovementWeightType", MovementWeightStr);
+   CHECK_ERROR_ABORT(Err,
+                     "VertCoord: MovementWeightType not found in VertCoord");
+
+   if (MovementWeightStr == "Fixed") {
+      MvmtWgtChoice = MovementWeightType::Fixed;
+   } else if (MovementWeightStr == "Uniform") {
+      MvmtWgtChoice = MovementWeightType::Uniform;
+   } else {
+      ABORT_ERROR("VertCoord: Unknown MovementWeightType requested");
+   }
+
+   // Fetch reference desnity from Config
+   Config TendConfig("Tendencies");
+   Err.reset();
+   Err += Options->get(TendConfig);
+   CHECK_ERROR_ABORT(Err, "VertCoord: Tendencies group not found in Config");
+
+   Err += TendConfig.get("Density0", Rho0);
+   CHECK_ERROR_ABORT(Err, "VertCoord: Density0 not found in TendConfig");
 
    // Store name suffix
    Name = Name_;
@@ -57,19 +88,36 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
    // Retrieve mesh filename from Decomp
    MeshFileName = Decomp->MeshFileName;
 
-   // Open the mesh file for reading (assume IO has already been initialized)
-   IO::openFile(MeshFileID, MeshFileName, IO::ModeRead);
-
-   // Set NVertLayers and NVertLayersP1 and create the vertical dimensions
-   Error Err; // Error code
-   I4 NVertLayersID;
-   Err = IO::getDimFromFile(MeshFileID, "nVertLevels", NVertLayersID,
-                            NVertLayers);
-   if (!Err.isSuccess()) {
-      LOG_WARN("VertCoord: error reading nVertLevels from mesh file, "
-               "using NVertLayers = 1");
-      NVertLayers = 1;
+   // If InNVertlayers is 0 (default), attempt to read the dimension from the
+   // mesh file. If the mesh file does not define a vertical dimension, use
+   // NVertLayers = 1. If a value for InNVertLayers is explicitly provided,
+   // use that value is instead.
+   if (InNVertLayers == 0) {
+
+      std::string OmegaDimName = "NVertLayers";
+      std::string MPASDimName  = "nVertLevels";
+
+      // Open the mesh file for reading (assume IO has already been initialized)
+      IO::openFile(MeshFileID, MeshFileName, IO::ModeRead);
+
+      // Set NVertLayers
+      I4 NVertLayersID;
+      Err = IO::getDimFromFile(MeshFileID, OmegaDimName, NVertLayersID,
+                               NVertLayers);
+      if (Err.isFail()) { // dim not found, try again with older MPAS name
+         Err = IO::getDimFromFile(MeshFileID, MPASDimName, NVertLayersID,
+                                  NVertLayers);
+         if (Err.isFail()) {
+            LOG_INFO("VertCoord: vertical dimension not found in mesh file, "
+                     "using NVertLayers = 1");
+            NVertLayers = 1;
+         }
+      }
+   } else {
+      NVertLayers = InNVertLayers;
    }
+
+   // Create the dimensions
    NVertLayersP1 = NVertLayers + 1;
 
    std::string DimName   = "NVertLayers";
@@ -96,7 +144,7 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
    NVerticesSize  = Decomp->NVerticesSize;
    VertexDegree   = Decomp->VertexDegree;
 
-   // Retrieve connectivity arrays from HorzMesh
+   // Retrieve connectivity arrays from Decomp
    CellsOnEdge   = Decomp->CellsOnEdge;
    CellsOnVertex = Decomp->CellsOnVertex;
 
@@ -124,114 +172,32 @@ VertCoord::VertCoord(const std::string &Name_, //< [in] Name for new VertCoord
    LayerThicknessTargetH = createHostMirrorCopy(LayerThicknessTarget);
    RefLayerThicknessH    = createHostMirrorCopy(RefLayerThickness);
 
-} // end constructor
-
-//------------------------------------------------------------------------------
-// Complete construction of new VertCoord instance
-void VertCoord::completeSetup(Config *Options //< [in] configuration options
-) {
-
    // Define field metadata
    defineFields();
 
-   I4 FillValueI4     = -1;
-   Real FillValueReal = -999._Real;
-
-   deepCopy(MinLayerCell, FillValueI4);
-   deepCopy(MaxLayerCell, FillValueI4);
-   deepCopy(BottomDepth, FillValueReal);
-
-   OMEGA_SCOPE(LocMinLayerCell, MinLayerCell);
-   OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell);
-   OMEGA_SCOPE(LocBottomDepth, BottomDepth);
-
-   // Fetch input stream and validate
-   std::string StreamName = "InitialVertCoord";
-   if (Name != "Default") {
-      StreamName.append(Name);
-   }
-
-   auto VCoordStream = IOStream::get(StreamName);
-
-   bool IsValidated = VCoordStream->validate();
-
-   // Read InitialVertCoord stream
-   Error Err; // Error code
-   if (IsValidated) {
-      Err = IOStream::read(StreamName);
-      if (!Err.isSuccess()) {
-         LOG_WARN("VertCoord: Error reading {} stream", StreamName);
-         I4 Sum1 = 0;
-         parallelReduce(
-             {MinLayerCell.extent_int(0)},
-             KOKKOS_LAMBDA(int I, int &Accum) { Accum += LocMinLayerCell(I); },
-             Sum1);
-         if (Sum1 < 0) {
-            LOG_WARN("VertCoord: Error reading minLevelCell from {}, "
-                     "using MinLayerCell = 0",
-                     StreamName);
-            deepCopy(MinLayerCell, 1);
-         }
-         I4 Sum2 = 0;
-         parallelReduce(
-             {MaxLayerCell.extent_int(0)},
-             KOKKOS_LAMBDA(int I, int &Accum) { Accum += LocMaxLayerCell(I); },
-             Sum2);
-         if (Sum2 < 0) {
-            LOG_WARN("VertCoord: Error reading maxLevelCell from {}, "
-                     "using MaxLayerCell = NVertLayers - 1",
-                     StreamName);
-            deepCopy(MaxLayerCell, NVertLayers);
-         }
-         Real Sum3 = 0.;
-         parallelReduce(
-             {BottomDepth.extent_int(0)},
-             KOKKOS_LAMBDA(int I, Real &Accum) { Accum += LocBottomDepth(I); },
-             Sum3);
-         if (Sum3 < 0.) {
-            ABORT_ERROR("VertCoord: Error reading bottomDepth from {}",
-                        StreamName);
-         }
-      }
-   } else {
-      ABORT_ERROR("Error validating IO stream {}", StreamName);
-   }
-
-   // Subtract 1 to convert to zero-based indexing
-   parallelFor(
-       {NCellsAll}, KOKKOS_LAMBDA(int ICell) {
-          LocMinLayerCell(ICell) -= 1;
-          LocMaxLayerCell(ICell) -= 1;
-       });
+   // Set MinLayerCell, MaxLayerCell, and BottomDepth arrays
+   setStreamArrays(ReadStream);
 
    // Compute Edge and Vertex vertical ranges
    minMaxLayerEdge();
    minMaxLayerVertex();
 
-   // Initialize movement weights
-   initMovementWeights(Options);
-
-   // Make host copies for device arrays read from mesh file
-   MaxLayerCellH = createHostMirrorCopy(MaxLayerCell);
-   MinLayerCellH = createHostMirrorCopy(MinLayerCell);
-   BottomDepthH  = createHostMirrorCopy(BottomDepth);
-
-   // Fetch reference desnity from Config
-   Config TendConfig("Tendencies");
-   Err.reset();
-   Err += Options->get(TendConfig);
-   CHECK_ERROR_ABORT(Err, "VertCoord: Tendencies group not found in Config");
+   // Set computational masks
+   setMasks();
 
-   Err += TendConfig.get("Density0", Rho0);
-   CHECK_ERROR_ABORT(Err, "VertCoord: Density0 not found in TendConfig");
+   // Initialize movement weights
+   initMovementWeights();
 
-} // end completeSetup
+} // end constructor
 
 //------------------------------------------------------------------------------
 // Calls the VertCoord constructor and places it in the AllVertCoords map
-VertCoord *
-VertCoord::create(const std::string &Name, // [in] name for new VertCoord
-                  const Decomp *Decomp     // [in] associated Decomp
+VertCoord *VertCoord::create(
+    const std::string &Name, // [in] name for new VertCoord
+    const Decomp *Decomp,    // [in] associated Decomp
+    Config *Options,         // [in] configuration options
+    const bool ReadStream,   // [in] optional logical to read stream
+    const int InNVertLayers  // [in] optional int to set vertical dim
 ) {
    // Check to see if a VertCoord of the same name already exists and, if so,
    // exit with an error
@@ -244,7 +210,8 @@ VertCoord::create(const std::string &Name, // [in] name for new VertCoord
 
    // create a new VertCoord on the heap and put it in a map of unique_ptrs,
    // which will manage its lifetime
-   auto *NewVertCoord = new VertCoord(Name, Decomp);
+   auto *NewVertCoord =
+       new VertCoord(Name, Decomp, Options, ReadStream, InNVertLayers);
    AllVertCoords.emplace(Name, NewVertCoord);
 
    return NewVertCoord;
@@ -255,8 +222,8 @@ VertCoord::create(const std::string &Name, // [in] name for new VertCoord
 void VertCoord::defineFields() {
 
    // Set field names (append Name if not default)
-   MinLayerCellFldName   = "MinLevelCell";
-   MaxLayerCellFldName   = "MaxLevelCell";
+   MinLayerCellFldName   = "MinLayerCell";
+   MaxLayerCellFldName   = "MaxLayerCell";
    BottomDepthFldName    = "BottomDepth";
    PressInterfFldName    = "PressureInterface";
    PressMidFldName       = "PressureMid";
@@ -478,6 +445,110 @@ void VertCoord::clear() {
 
 } // end clear
 
+//------------------------------------------------------------------------------
+// If ReadStream = true, read MinLayerCell, MaxLayerCell, and BottomDepth from
+// the initial stream. If ReadStream = false, default values are used.
+void VertCoord::setStreamArrays(const bool ReadStream) {
+
+   Error Err; // Error code
+
+   OMEGA_SCOPE(LocMinLayerCell, MinLayerCell);
+   OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell);
+   OMEGA_SCOPE(LocBottomDepth, BottomDepth);
+
+   // If ReadStream is true (default) attempt to read values for MinLayerCell,
+   // MaxLayerCell, and BottomDepth from the InitialVertCoord stream. Otherwise,
+   // MinLayerCell and MaxLayerCell will be set to the first and last indices of
+   // the vertical range, BottomDepth will remain uninitialized and will need
+   // to be initialized explicitly if needed.
+   if (ReadStream) {
+
+      I4 FillValueI4     = -1;
+      Real FillValueReal = -999._Real;
+
+      deepCopy(MinLayerCell, FillValueI4);
+      deepCopy(MaxLayerCell, FillValueI4);
+      deepCopy(BottomDepth, FillValueReal);
+
+      // Fetch input stream and validate
+      std::string StreamName = "InitialVertCoord";
+      if (Name != "Default") {
+         StreamName.append(Name);
+      }
+
+      // Validate InitalVertCoord stream
+      auto VCoordStream = IOStream::get(StreamName);
+      bool IsValidated  = VCoordStream->validate();
+
+      // Read InitialVertCoord stream
+      if (IsValidated) {
+         // Attempt to read stream, an error will be raised if any field fails
+         // to be read. Determine which fields may not have been read properly.
+         // If MinLayerCell or MaxLayerCell were not read properly default
+         // values will be used. If BottomDepth was not read properly, abort
+         // with error.
+         Err = IOStream::read(StreamName);
+         if (Err.isFail()) {
+            LOG_INFO("VertCoord: Error while reading {} stream", StreamName);
+            I4 Sum1 = 0;
+            parallelReduce(
+                {MinLayerCell.extent_int(0)},
+                KOKKOS_LAMBDA(int I, int &Accum) {
+                   Accum += LocMinLayerCell(I);
+                },
+                Sum1);
+            if (Sum1 < 0) {
+               LOG_INFO("VertCoord: Error reading MinLayerCell from {}, "
+                        "using MinLayerCell = 0",
+                        StreamName);
+               deepCopy(MinLayerCell, 1);
+            }
+            I4 Sum2 = 0;
+            parallelReduce(
+                {MaxLayerCell.extent_int(0)},
+                KOKKOS_LAMBDA(int I, int &Accum) {
+                   Accum += LocMaxLayerCell(I);
+                },
+                Sum2);
+            if (Sum2 < 0) {
+               LOG_INFO("VertCoord: Error reading MaxLayerCell from {}, "
+                        "using MaxLayerCell = NVertLayers - 1",
+                        StreamName);
+               deepCopy(MaxLayerCell, NVertLayers);
+            }
+            Real Sum3 = 0.;
+            parallelReduce(
+                {BottomDepth.extent_int(0)},
+                KOKKOS_LAMBDA(int I, Real &Accum) {
+                   Accum += LocBottomDepth(I);
+                },
+                Sum3);
+            if (Sum3 < 0.) {
+               ABORT_ERROR("VertCoord: Error reading bottomDepth from {}",
+                           StreamName);
+            }
+         }
+      } else {
+         ABORT_ERROR("Error validating IO stream {}", StreamName);
+      }
+   } else {
+      deepCopy(MinLayerCell, 1);
+      deepCopy(MaxLayerCell, NVertLayers);
+   }
+
+   // Subtract 1 to convert to zero-based indexing
+   parallelFor(
+       {NCellsAll}, KOKKOS_LAMBDA(int ICell) {
+          LocMinLayerCell(ICell) -= 1;
+          LocMaxLayerCell(ICell) -= 1;
+       });
+
+   // Make host copies for device arrays read from mesh file
+   MaxLayerCellH = createHostMirrorCopy(MaxLayerCell);
+   MinLayerCellH = createHostMirrorCopy(MinLayerCell);
+   BottomDepthH  = createHostMirrorCopy(BottomDepth);
+}
+
 //------------------------------------------------------------------------------
 // Compute min and max layer indices for edges based on MinLayerCell and
 // MaxLayerCell
@@ -610,36 +681,99 @@ void VertCoord::minMaxLayerVertex() {
 } // end MinMaxLayerVertex
 
 //------------------------------------------------------------------------------
-// Store VertCoordMovementWeights based on config choice
-void VertCoord::initMovementWeights(
-    Config *Options // [in] configuration options
-) {
+// set computational masks for mesh elements
+void VertCoord::setMasks() {
 
-   Error Err; // default successful error code
+   EdgeMask = Array2DReal("EdgeMask", NEdgesSize, NVertLayers);
 
-   Config VCoordConfig("VertCoord");
-   Err += Options->get(VCoordConfig);
-   CHECK_ERROR_ABORT(Err, "VertCoord: VertCoord group not found in Config");
+   OMEGA_SCOPE(LocEdgeMask, EdgeMask);
+   OMEGA_SCOPE(LocMinLyrEdgeBot, MinLayerEdgeBot);
+   OMEGA_SCOPE(LocMaxLyrEdgeTop, MaxLayerEdgeTop);
 
-   std::string MovementWeightType;
-   Err += VCoordConfig.get("MovementWeightType", MovementWeightType);
-   CHECK_ERROR_ABORT(Err,
-                     "VertCoord: MovementWeightType not found in VertCoord");
+   // EdgeMask = 1 if active layers on both sides, 0 otherwise.
+   deepCopy(EdgeMask, 0.);
+   parallelForOuter(
+       {NEdgesAll}, KOKKOS_LAMBDA(int IEdge, const TeamMember &Team) {
+          const I4 KMin = LocMinLyrEdgeBot(IEdge);
+          const I4 KMax = LocMaxLyrEdgeTop(IEdge);
+
+          parallelForInner(
+              Team, KMax - KMin + 1, INNER_LAMBDA(int K) {
+                 I4 KLyr = KMin + K;
+
+                 LocEdgeMask(IEdge, KLyr) = 1._Real;
+              });
+       });
+
+   EdgeMaskH = createHostMirrorCopy(EdgeMask);
+
+   CellMask = Array2DReal("CellMask", NCellsSize, NVertLayers);
+
+   OMEGA_SCOPE(LocCellMask, CellMask);
+   OMEGA_SCOPE(LocMinLyrCell, MinLayerCell);
+   OMEGA_SCOPE(LocMaxLyrCell, MaxLayerCell);
+
+   // CellMask = 1 in active layers, 0 otherwise.
+   deepCopy(CellMask, 0.);
+   parallelForOuter(
+       {NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) {
+          const I4 KMin = LocMinLyrCell(ICell);
+          const I4 KMax = LocMaxLyrCell(ICell);
+
+          parallelForInner(
+              Team, KMax - KMin + 1, INNER_LAMBDA(int K) {
+                 I4 KLyr = KMin + K;
+
+                 LocCellMask(ICell, KLyr) = 1._Real;
+              });
+       });
+
+   CellMaskH = createHostMirrorCopy(CellMask);
+
+   VertexMask = Array2DReal("VertexMask", NVerticesSize, NVertLayers);
+
+   OMEGA_SCOPE(LocVrtxMask, VertexMask);
+   OMEGA_SCOPE(LocMinLyrVrtxTop, MinLayerVertexTop);
+   OMEGA_SCOPE(LocMaxLyrVrtxBot, MaxLayerVertexBot);
+
+   // VertexMask = 1 if at least 1 surrounding cell layer is active,
+   // 0 otherwise.
+   deepCopy(VertexMask, 0.);
+   parallelForOuter(
+       {NVerticesAll}, KOKKOS_LAMBDA(int IVertex, const TeamMember &Team) {
+          const I4 KMin = LocMinLyrVrtxTop(IVertex);
+          const I4 KMax = LocMaxLyrVrtxBot(IVertex);
+
+          parallelForInner(
+              Team, KMax - KMin + 1, INNER_LAMBDA(int K) {
+                 I4 KLyr = KMin + K;
+
+                 LocVrtxMask(IVertex, KLyr) = 1._Real;
+              });
+       });
+
+   VertexMaskH = createHostMirrorCopy(VertexMask);
+
+} // end setMasks()
+
+//------------------------------------------------------------------------------
+// Store VertCoordMovementWeights based on config choice
+void VertCoord::initMovementWeights() {
+
+   Error Err; // default successful error code
 
    VertCoordMovementWeights =
        Array1DReal("VertCoordMovementWeights", NVertLayers);
 
    OMEGA_SCOPE(LocVertCoordMovementWeights, VertCoordMovementWeights);
-   if (MovementWeightType == "Fixed") {
+   if (MvmtWgtChoice == MovementWeightType::Fixed) {
       deepCopy(VertCoordMovementWeights, 0._Real);
       parallelFor(
           {1}, KOKKOS_LAMBDA(const int &) {
              LocVertCoordMovementWeights(0) = 1._Real;
           });
-   } else if (MovementWeightType == "Uniform") {
+   } else if (MvmtWgtChoice == MovementWeightType::Uniform) {
       deepCopy(VertCoordMovementWeights, 1._Real);
-   } else {
-      ABORT_ERROR("VertCoord: Unknown MovementWeightType requested");
    }
 
    VertCoordMovementWeightsH = createHostMirrorCopy(VertCoordMovementWeights);
diff --git a/components/omega/src/ocn/VertCoord.h b/components/omega/src/ocn/VertCoord.h
index 6423c5fa9cf1..04576d45104a 100644
--- a/components/omega/src/ocn/VertCoord.h
+++ b/components/omega/src/ocn/VertCoord.h
@@ -25,6 +25,11 @@
 
 namespace OMEGA {
 
+enum class MovementWeightType {
+   Fixed,  /// Distribute perturbations in top level
+   Uniform /// Uniform stretching
+};
+
 class VertCoord {
 
  private:
@@ -42,6 +47,9 @@ class VertCoord {
    Array2DI4 CellsOnEdge;
    Array2DI4 CellsOnVertex;
 
+   // Choice of VertCoorMovementWeight type
+   MovementWeightType MvmtWgtChoice;
+
    std::string MeshFileName;
    int MeshFileID;
 
@@ -51,8 +59,11 @@ class VertCoord {
    // methods
 
    /// construct a new vertical coordinate object
-   VertCoord(const std::string &Name, ///< [in] Name for new VertCoord
-             const Decomp *MeshDecomp ///< [in] associated Decomp
+   VertCoord(const std::string &Name,  ///< [in] Name for new VertCoord
+             const Decomp *MeshDecomp, ///< [in] associated Decomp
+             Config *Options,          ///< [in] configuration options
+             const bool ReadStream,    ///< [in] logical to read stream
+             const int NVertLayers     ///< [in] int to set vertical dim
    );
 
    /// define field metadata
@@ -112,6 +123,20 @@ class VertCoord {
    HostArray1DReal VertCoordMovementWeightsH;
    HostArray2DReal RefLayerThicknessH;
 
+   // Masks
+   Array2DReal EdgeMask;        ///< Mask to determine if computations should be
+                                ///  done on edge
+   HostArray2DReal EdgeMaskH;   ///< Mask to determine if computations should be
+                                ///  done on edge
+   Array2DReal CellMask;        ///< Mask to determine if computations should be
+                                ///  done on cell
+   HostArray2DReal CellMaskH;   ///< Mask to determine if computations should be
+                                ///  done on cell
+   Array2DReal VertexMask;      ///< Mask to determine if computations should be
+                                ///  done on vertex
+   HostArray2DReal VertexMaskH; ///< Mask to determine if computations should be
+                                ///  done on vertex
+
    // BottomDepth read from mesh file
    Array1DReal BottomDepth;
 
@@ -137,22 +162,17 @@ class VertCoord {
 
    // methods
 
-   /// 1st phase of initialization for default vertical coordinate
-   static void init1();
-
-   /// 2nd phase of initialization for default vertical coordinate
-   static void init2();
+   /// Initialization of default vertical coordinate
+   static void init(const bool ReadStream = true, const int NVertLayers = 0);
 
    /// Creates a new vertical coordinate object by calling the constructor and
-   /// puts it in the AllVertCoords map. This object is mostly empty and must
-   /// completed by completeSetup.
+   /// puts it in the AllVertCoords map.
    static VertCoord *
-   create(const std::string &Name, /// [in] name for new VertCoord
-          const Decomp *MeshDecomp /// [in] associated Decomp
-   );
-
-   /// Read InitialVertCoord stream and complete initialization
-   void completeSetup(Config *Options /// [in] configuration options
+   create(const std::string &Name,      /// [in] name for new VertCoord
+          const Decomp *MeshDecomp,     /// [in] associated Decomp
+          Config *Options,              /// [in] configuration options
+          const bool ReadStream = true, /// [in] optional logical to read stream
+          const int NVertLayers = 0 /// [in] optional int to set vertical dim
    );
 
    /// Copy member arrays from device to host
@@ -176,11 +196,14 @@ class VertCoord {
    /// Retreive a VertCoord by name
    static VertCoord *get(std::string name);
 
-   // Variable initialization methods
+   // Array initialization methods
+   void setStreamArrays(const bool ReadStream);
    void minMaxLayerEdge();
    void minMaxLayerVertex();
-   void initMovementWeights(Config *Options ///< [in] configuration options
-   );
+   void initMovementWeights();
+
+   /// Initialize computational masks
+   void setMasks();
 
    /// Sums the mass thickness times g from the top layer down, starting with
    /// the surface pressure
diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp
index 1ef7eefa6394..5a6dbeea706c 100644
--- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp
+++ b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.cpp
@@ -6,8 +6,8 @@
 namespace OMEGA {
 
 TracerAuxVars::TracerAuxVars(const std::string &AuxStateSuffix,
-                             const HorzMesh *Mesh, const I4 NVertLayers,
-                             const I4 NTracers)
+                             const HorzMesh *Mesh, const VertCoord *VCoord,
+                             const I4 NVertLayers, const I4 NTracers)
     : HTracersEdge("ThickTracersEdge" + AuxStateSuffix, NTracers,
                    Mesh->NEdgesSize, NVertLayers),
       Del2TracersCell("Del2TracerCell" + AuxStateSuffix, NTracers,
@@ -15,7 +15,7 @@ TracerAuxVars::TracerAuxVars(const std::string &AuxStateSuffix,
       NEdgesOnCell(Mesh->NEdgesOnCell), EdgesOnCell(Mesh->EdgesOnCell),
       CellsOnEdge(Mesh->CellsOnEdge), EdgeSignOnCell(Mesh->EdgeSignOnCell),
       DcEdge(Mesh->DcEdge), DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell),
-      EdgeMask(Mesh->EdgeMask) {}
+      EdgeMask(VCoord->EdgeMask) {}
 
 void TracerAuxVars::registerFields(const std::string &AuxGroupName,
                                    const std::string &MeshName) const {
diff --git a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h
index 39e6f5bf232a..8eb8eeea7219 100644
--- a/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h
+++ b/components/omega/src/ocn/auxiliaryVars/TracerAuxVars.h
@@ -5,6 +5,7 @@
 #include "Field.h"
 #include "HorzMesh.h"
 #include "OmegaKokkos.h"
+#include "VertCoord.h"
 
 #include <string>
 
@@ -20,7 +21,8 @@ class TracerAuxVars {
    FluxTracerEdgeOption TracersOnEdgeChoice;
 
    TracerAuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh,
-                 const I4 NVertLayers, const I4 NTracers);
+                 const VertCoord *VCoord, const I4 NVertLayers,
+                 const I4 NTracers);
 
    KOKKOS_FUNCTION void computeVarsOnEdge(int L, int IEdge, int KChunk,
                                           const Array2DReal &NormalVelEdge,
diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp
index 8140723b83f9..cc7a77170fc0 100644
--- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp
+++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.cpp
@@ -7,7 +7,9 @@
 namespace OMEGA {
 
 VelocityDel2AuxVars::VelocityDel2AuxVars(const std::string &AuxStateSuffix,
-                                         const HorzMesh *Mesh, int NVertLayers)
+                                         const HorzMesh *Mesh,
+                                         const VertCoord *VCoord,
+                                         int NVertLayers)
     : Del2Edge("VelDel2Edge" + AuxStateSuffix, Mesh->NEdgesSize, NVertLayers),
       Del2DivCell("VelDel2DivCell" + AuxStateSuffix, Mesh->NCellsSize,
                   NVertLayers),
@@ -17,7 +19,7 @@ VelocityDel2AuxVars::VelocityDel2AuxVars(const std::string &AuxStateSuffix,
       EdgeSignOnCell(Mesh->EdgeSignOnCell), DcEdge(Mesh->DcEdge),
       DvEdge(Mesh->DvEdge), AreaCell(Mesh->AreaCell),
       EdgesOnVertex(Mesh->EdgesOnVertex), CellsOnEdge(Mesh->CellsOnEdge),
-      VerticesOnEdge(Mesh->VerticesOnEdge), EdgeMask(Mesh->EdgeMask),
+      VerticesOnEdge(Mesh->VerticesOnEdge), EdgeMask(VCoord->EdgeMask),
       EdgeSignOnVertex(Mesh->EdgeSignOnVertex),
       AreaTriangle(Mesh->AreaTriangle), VertexDegree(Mesh->VertexDegree) {}
 
diff --git a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h
index 468d5ea9279e..4cd76115a5d4 100644
--- a/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h
+++ b/components/omega/src/ocn/auxiliaryVars/VelocityDel2AuxVars.h
@@ -4,6 +4,7 @@
 #include "DataTypes.h"
 #include "HorzMesh.h"
 #include "OmegaKokkos.h"
+#include "VertCoord.h"
 
 #include <string>
 
@@ -16,7 +17,7 @@ class VelocityDel2AuxVars {
    Array2DReal Del2RelVortVertex;
 
    VelocityDel2AuxVars(const std::string &AuxStateSuffix, const HorzMesh *Mesh,
-                       int NVertLayers);
+                       const VertCoord *VCoord, int NVertLayers);
 
    KOKKOS_FUNCTION void
    computeVarsOnEdge(int IEdge, int KChunk, const Array2DReal &VelocityDivCell,
diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp
index dfe9710ea1e1..a193670185f7 100644
--- a/components/omega/test/infra/IOStreamTest.cpp
+++ b/components/omega/test/infra/IOStreamTest.cpp
@@ -91,21 +91,12 @@ void initIOStreamTest(Clock *&ModelClock // Model clock
    // Initialize IOStreams
    IOStream::init(ModelClock);
 
-   // Initialize the vertical coordinate (phase 1)
-   VertCoord::init1();
-
-   // Reset vertical layers dimension
-   Dimension::destroy("NVertLayers");
-   I4 NVertLayers = 60;
-   std::shared_ptr<Dimension> VertDim =
-       Dimension::create("NVertLayers", NVertLayers);
-
    // Initialize HorzMesh - this should read Mesh stream
    HorzMesh::init();
    HorzMesh *DefMesh = HorzMesh::getDefault();
 
-   // Initialize the vertical coordinate (phase 2)
-   VertCoord::init2();
+   // Initialize the vertical coordinate
+   VertCoord::init();
 
    // Initialize State
    TmpErr = OceanState::init();
diff --git a/components/omega/test/ocn/AuxiliaryStateTest.cpp b/components/omega/test/ocn/AuxiliaryStateTest.cpp
index 217c21de46df..28675ddde0b0 100644
--- a/components/omega/test/ocn/AuxiliaryStateTest.cpp
+++ b/components/omega/test/ocn/AuxiliaryStateTest.cpp
@@ -108,9 +108,10 @@ int initAuxStateTest(const std::string &mesh) {
       LOG_ERROR("AuxStateTest: error initializing default halo");
    }
 
-   VertCoord::init1();
    HorzMesh::init();
 
+   VertCoord::init(false);
+
    Tracers::init();
    int StateErr = OceanState::init();
    if (StateErr != 0) {
@@ -138,10 +139,11 @@ int testAuxState() {
       return -1;
    }
 
-   const auto *Mesh = HorzMesh::getDefault();
-   auto *MeshHalo   = Halo::getDefault();
+   const auto *Mesh   = HorzMesh::getDefault();
+   const auto *VCoord = VertCoord::getDefault();
+   auto *MeshHalo     = Halo::getDefault();
    // test creation of another auxiliary state
-   AuxiliaryState::create("AnotherAuxState", Mesh, MeshHalo, 12, 3);
+   AuxiliaryState::create("AnotherAuxState", Mesh, VCoord, MeshHalo, 12, 3);
 
    // test retrievel of another
    if (AuxiliaryState::get("AnotherAuxState")) {
@@ -279,11 +281,10 @@ void finalizeAuxStateTest() {
    Tracers::clear();
    OceanState::clear();
    VertCoord::clear();
+   HorzMesh::clear();
    Field::clear();
    Dimension::clear();
    TimeStepper::clear();
-   HorzMesh::clear();
-   VertCoord::clear();
    Halo::clear();
    Decomp::clear();
    MachEnv::removeAll();
diff --git a/components/omega/test/ocn/AuxiliaryVarsTest.cpp b/components/omega/test/ocn/AuxiliaryVarsTest.cpp
index 00505498439c..d87b4c0e73aa 100644
--- a/components/omega/test/ocn/AuxiliaryVarsTest.cpp
+++ b/components/omega/test/ocn/AuxiliaryVarsTest.cpp
@@ -615,7 +615,8 @@ int testVelocityDel2AuxVars(Real RTol) {
 
    const auto Decomp = Decomp::getDefault();
    const auto Mesh   = HorzMesh::getDefault();
-   VelocityDel2AuxVars VelocityDel2Aux("", Mesh, NVertLayers);
+   const auto VCoord = VertCoord::getDefault();
+   VelocityDel2AuxVars VelocityDel2Aux("", Mesh, VCoord, NVertLayers);
 
    // Use analytical expressions to compute inputs
 
@@ -721,9 +722,10 @@ int testTracerAuxVars(const Array2DReal &LayerThickCell,
    TestSetup Setup;
    int Err = 0;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
-   TracerAuxVars TracerAux("", Mesh, NVertLayers, NTracers);
+   TracerAuxVars TracerAux("", Mesh, VCoord, NVertLayers, NTracers);
    TracerAux.TracersOnEdgeChoice = FluxTracerEdgeOption::Upwind;
 
    // Set input arrays
@@ -823,25 +825,20 @@ int initAuxVarsTest(const std::string &mesh) {
       LOG_ERROR("AuxVarsTest: error initializing default halo");
    }
 
-   VertCoord::init1();
-   // Reset NVertLayers to the test value
-   auto *DefVertCoord        = VertCoord::getDefault();
-   DefVertCoord->NVertLayers = NVertLayers;
-   Dimension::destroy("NVertLayers");
-   std::shared_ptr<Dimension> VertDim =
-       Dimension::create("NVertLayers", NVertLayers);
-
    HorzMesh::init();
 
+   // initialize vertical coordinate, do not read stream and use local
+   // NVertLayers value
+   VertCoord::init(false, NVertLayers);
+
    return Err;
 }
 
 void finalizeAuxVarsTest() {
    VertCoord::clear();
+   HorzMesh::clear();
    Field::clear();
    Dimension::clear();
-   HorzMesh::clear();
-   VertCoord::clear();
    Halo::clear();
    Decomp::clear();
    MachEnv::removeAll();
diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp
index 028d136e9c16..db7aa7ea421c 100644
--- a/components/omega/test/ocn/EosTest.cpp
+++ b/components/omega/test/ocn/EosTest.cpp
@@ -74,12 +74,12 @@ I4 initEosTest(const std::string &mesh) {
    /// Initialize decomposition
    Decomp::init(mesh);
 
-   /// Initialize vertical coordinate (phase 1)
-   VertCoord::init1();
-
    /// Initialize mesh
    HorzMesh::init();
 
+   /// Initialize vertical coordinate
+   VertCoord::init(false);
+
    /// Initialize Eos
    Eos::init();
 
@@ -286,8 +286,8 @@ int testEosTeos10Displaced() {
 
 /// Finalize and clean up all test infrastructure
 void finalizeEosTest() {
-   HorzMesh::clear();
    VertCoord::clear();
+   HorzMesh::clear();
    Decomp::clear();
    Field::clear();
    Dimension::clear();
diff --git a/components/omega/test/ocn/HorzMeshTest.cpp b/components/omega/test/ocn/HorzMeshTest.cpp
index 5713bb4c77a0..0c218325b100 100644
--- a/components/omega/test/ocn/HorzMeshTest.cpp
+++ b/components/omega/test/ocn/HorzMeshTest.cpp
@@ -20,7 +20,6 @@
 #include "MachEnv.h"
 #include "OmegaKokkos.h"
 #include "Pacer.h"
-#include "VertCoord.h"
 #include "mpi.h"
 
 #include <iostream>
@@ -60,9 +59,6 @@ int initHorzMeshTest() {
    if (Err != 0)
       LOG_ERROR("HorzMeshTest: error initializing default halo");
 
-   // Initialize the vertical coordinate (phase 1)
-   VertCoord::init1();
-
    // Initialize the default mesh
    HorzMesh::init();
 
@@ -773,7 +769,6 @@ int main(int argc, char *argv[]) {
       }
       // Finalize Omega objects
       HorzMesh::clear();
-      VertCoord::clear();
       Dimension::clear();
       Halo::clear();
       Decomp::clear();
diff --git a/components/omega/test/ocn/HorzOperatorsTest.cpp b/components/omega/test/ocn/HorzOperatorsTest.cpp
index ca57e8e3ed8e..86fc7c30fce4 100644
--- a/components/omega/test/ocn/HorzOperatorsTest.cpp
+++ b/components/omega/test/ocn/HorzOperatorsTest.cpp
@@ -11,7 +11,6 @@
 #include "OceanTestCommon.h"
 #include "OmegaKokkos.h"
 #include "Pacer.h"
-#include "VertCoord.h"
 #include "mpi.h"
 
 #include <cmath>
@@ -450,8 +449,6 @@ int initOperatorsTest(const std::string &MeshFile) {
       LOG_ERROR("OperatorsTest: error initializing default halo");
    }
 
-   VertCoord::init1();
-
    HorzMesh::init();
 
    return Err;
@@ -459,7 +456,6 @@ int initOperatorsTest(const std::string &MeshFile) {
 
 void finalizeOperatorsTest() {
    HorzMesh::clear();
-   VertCoord::clear();
    Dimension::clear();
    Halo::clear();
    Decomp::clear();
diff --git a/components/omega/test/ocn/StateTest.cpp b/components/omega/test/ocn/StateTest.cpp
index 2ae92f1641f9..b83333785455 100644
--- a/components/omega/test/ocn/StateTest.cpp
+++ b/components/omega/test/ocn/StateTest.cpp
@@ -78,14 +78,11 @@ void initStateTest() {
    if (Err != 0)
       ABORT_ERROR("State: error initializing default halo");
 
-   // Initialize the vertical coordinate (phase 1)
-   VertCoord::init1();
-
    // Initialize the default mesh
    HorzMesh::init();
 
-   // Initialize the vertical coordinate (phase 2)
-   VertCoord::init2();
+   // Initialize the vertical coordinate
+   VertCoord::init();
 
    // Initialize tracers
    Tracers::init();
diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp
index 2785876987ae..5239a243e538 100644
--- a/components/omega/test/ocn/TendenciesTest.cpp
+++ b/components/omega/test/ocn/TendenciesTest.cpp
@@ -110,8 +110,8 @@ int initTendenciesTest(const std::string &mesh) {
       LOG_ERROR("TendenciesTest: error initializing default halo");
    }
 
-   VertCoord::init1();
    HorzMesh::init();
+   VertCoord::init(false);
    Tracers::init();
 
    int StateErr = OceanState::init();
@@ -221,11 +221,10 @@ void finalizeTendenciesTest() {
    AuxiliaryState::clear();
    OceanState::clear();
    VertCoord::clear();
+   HorzMesh::clear();
    Field::clear();
    Dimension::clear();
    TimeStepper::clear();
-   HorzMesh::clear();
-   VertCoord::clear();
    Halo::clear();
    Decomp::clear();
    MachEnv::removeAll();
diff --git a/components/omega/test/ocn/TendencyTermsTest.cpp b/components/omega/test/ocn/TendencyTermsTest.cpp
index e537e9e6341a..468feeb9c3b0 100644
--- a/components/omega/test/ocn/TendencyTermsTest.cpp
+++ b/components/omega/test/ocn/TendencyTermsTest.cpp
@@ -393,7 +393,8 @@ int testPotVortHAdv(int NVertLayers, Real RTol) {
    int Err = 0;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    // Compute exact result
    Array2DReal ExactPotVortHAdv("ExactPotVortHAdv", Mesh->NEdgesOwned,
@@ -441,7 +442,7 @@ int testPotVortHAdv(int NVertLayers, Real RTol) {
    // Compute numerical result
    Array2DReal NumPotVortHAdv("NumPotVortHAdv", Mesh->NEdgesOwned, NVertLayers);
 
-   PotentialVortHAdvOnEdge PotVortHAdvOnE(Mesh);
+   PotentialVortHAdvOnEdge PotVortHAdvOnE(Mesh, VCoord);
    parallelFor(
        {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int KLayer) {
           PotVortHAdvOnE(NumPotVortHAdv, IEdge, KLayer, NormRelVortEdge,
@@ -469,7 +470,8 @@ int testKEGrad(int NVertLayers, Real RTol) {
    int Err = 0;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    // Compute exact result
    Array2DReal ExactKEGrad("ExactKEGrad", Mesh->NEdgesOwned, NVertLayers);
@@ -491,7 +493,7 @@ int testKEGrad(int NVertLayers, Real RTol) {
    // Compute numerical result
    Array2DReal NumKEGrad("NumKEGrad", Mesh->NEdgesOwned, NVertLayers);
 
-   KEGradOnEdge KEGradOnE(Mesh);
+   KEGradOnEdge KEGradOnE(Mesh, VCoord);
    parallelFor(
        {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int KLayer) {
           KEGradOnE(NumKEGrad, IEdge, KLayer, KECell);
@@ -517,7 +519,8 @@ int testSSHGrad(int NVertLayers, Real RTol) {
    int Err = 0;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    // Compute exact result
    Array2DReal ExactSSHGrad("ExactSSHGrad", Mesh->NEdgesOwned, NVertLayers);
@@ -539,7 +542,7 @@ int testSSHGrad(int NVertLayers, Real RTol) {
    // Compute numerical result
    Array2DReal NumSSHGrad("NumSSHGrad", Mesh->NEdgesOwned, NVertLayers);
 
-   SSHGradOnEdge SSHGradOnE(Mesh);
+   SSHGradOnEdge SSHGradOnE(Mesh, VCoord);
    parallelFor(
        {Mesh->NEdgesOwned, NVertLayers}, KOKKOS_LAMBDA(int IEdge, int KLayer) {
           SSHGradOnE(NumSSHGrad, IEdge, KLayer, SSHCell);
@@ -562,14 +565,15 @@ int testVelDiff(int NVertLayers, Real RTol) {
    Error Err1;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    Config *OmegaConfig = Config::getOmegaConfig();
    Config TendConfig("Tendencies");
    Err1 += OmegaConfig->get(TendConfig);
    CHECK_ERROR_ABORT(Err1, "Tendencies: Tendencies group not found in Config");
 
-   VelocityDiffusionOnEdge VelDiffOnE(Mesh);
+   VelocityDiffusionOnEdge VelDiffOnE(Mesh, VCoord);
    Err1 += TendConfig.get("ViscDel2", VelDiffOnE.ViscDel2);
    CHECK_ERROR_ABORT(Err1, "Tendencies: ViscDel2 not found in TendConfig");
 
@@ -627,14 +631,15 @@ int testVelHyperDiff(int NVertLayers, Real RTol) {
    Error Err1;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    Config *OmegaConfig = Config::getOmegaConfig();
    Config TendConfig("Tendencies");
    Err1 += OmegaConfig->get(TendConfig);
    CHECK_ERROR_ABORT(Err1, "Tendencies: Tendencies group not found in Config");
 
-   VelocityHyperDiffOnEdge VelHyperDiffOnE(Mesh);
+   VelocityHyperDiffOnEdge VelHyperDiffOnE(Mesh, VCoord);
    Err1 += TendConfig.get("ViscDel4", VelHyperDiffOnE.ViscDel4);
    CHECK_ERROR_ABORT(Err1, "Tendencies: ViscDel4 not found in TendConfig");
 
@@ -697,7 +702,8 @@ int testWindForcing(int NVertLayers) {
    int Err = 0;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    const Real SaltWaterDensity = 0.987654321;
 
@@ -737,7 +743,7 @@ int testWindForcing(int NVertLayers) {
    // Compute numerical result
    Array2DReal NumWindForcing("NumWindForcing", Mesh->NEdgesOwned, NVertLayers);
 
-   WindForcingOnEdge WindForcingOnE(Mesh);
+   WindForcingOnEdge WindForcingOnE(Mesh, VCoord);
    WindForcingOnE.SaltWaterDensity = SaltWaterDensity;
 
    parallelFor(
@@ -840,7 +846,8 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) {
    I4 Err = 0;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    // Compute exact result
    Array3DReal ExactTrFluxDiv("ExactTrFluxDiv", NTracers, Mesh->NCellsOwned,
@@ -869,7 +876,7 @@ int testTracerHorzAdvOnCell(int NVertLayers, int NTracers, Real RTol) {
    // Compute numerical result
    Array3DReal NumTrFluxDiv("NumTrFluxDiv", NTracers, Mesh->NCellsOwned,
                             NVertLayers);
-   TracerHorzAdvOnCell TrHorzAdvOnC(Mesh);
+   TracerHorzAdvOnCell TrHorzAdvOnC(Mesh, VCoord);
    parallelFor(
        {NTracers, Mesh->NCellsOwned, NVertLayers},
        KOKKOS_LAMBDA(int L, int ICell, int KLayer) {
@@ -896,7 +903,8 @@ int testTracerDiffOnCell(int NVertLayers, int NTracers, Real RTol) {
    I4 Err = 0;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    // Compute exact result
    Array3DReal ExactTracerDiff("ExactTracerDiff", NTracers, Mesh->NCellsOwned,
@@ -923,7 +931,7 @@ int testTracerDiffOnCell(int NVertLayers, int NTracers, Real RTol) {
    // Compute numerical result
    Array3DReal NumTracerDiff("NumTracerDiff", NTracers, Mesh->NCellsOwned,
                              NVertLayers);
-   TracerDiffOnCell TrDiffOnC(Mesh);
+   TracerDiffOnCell TrDiffOnC(Mesh, VCoord);
    TrDiffOnC.EddyDiff2 = 1._Real;
 
    parallelFor(
@@ -952,7 +960,8 @@ int testTracerHyperDiffOnCell(int NVertLayers, int NTracers, Real RTol) {
    I4 Err = 0;
    TestSetup Setup;
 
-   const auto Mesh = HorzMesh::getDefault();
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
 
    // Compute exact result
    Array3DReal ExactTracerHyperDiff("ExactTracerHyperDiff", NTracers,
@@ -973,7 +982,7 @@ int testTracerHyperDiffOnCell(int NVertLayers, int NTracers, Real RTol) {
    // Compute numerical result
    Array3DReal NumTracerHyperDiff("NumTracerHyperDiff", NTracers,
                                   Mesh->NCellsOwned, NVertLayers);
-   TracerHyperDiffOnCell TrHypDiffOnC(Mesh);
+   TracerHyperDiffOnCell TrHypDiffOnC(Mesh, VCoord);
    TrHypDiffOnC.EddyDiff4 = 1._Real;
    parallelFor(
        {NTracers, Mesh->NCellsOwned, NVertLayers},
@@ -1019,17 +1028,12 @@ void initTendTest(const std::string &MeshFile, int NVertLayers) {
       ABORT_ERROR("TendencyTermsTest: error initializing default halo");
    }
 
-   VertCoord::init1();
-
-   // Reset NVertLayers to the test value
-   auto *DefVertCoord        = VertCoord::getDefault();
-   DefVertCoord->NVertLayers = NVertLayers;
-   Dimension::destroy("NVertLayers");
-   std::shared_ptr<Dimension> VertDim =
-       Dimension::create("NVertLayers", NVertLayers);
-
    HorzMesh::init();
 
+   // initialize vertical coordinate, do not read stream and use local
+   // NVertLayers value
+   VertCoord::init(false, NVertLayers);
+
 } // end initTendTest
 
 void finalizeTendTest() {
diff --git a/components/omega/test/ocn/TracersTest.cpp b/components/omega/test/ocn/TracersTest.cpp
index 7e1f541c1bb3..daf85e6f3536 100644
--- a/components/omega/test/ocn/TracersTest.cpp
+++ b/components/omega/test/ocn/TracersTest.cpp
@@ -69,12 +69,12 @@ I4 initTracersTest() {
       return Err;
    }
 
-   // Initialize the vertical coordinate (phase 1)
-   VertCoord::init1();
-
    // Initialize the default mesh
    HorzMesh::init();
 
+   // Initialize the vertical coordinate
+   VertCoord::init(false);
+
    return 0;
 }
 
diff --git a/components/omega/test/ocn/VertCoordTest.cpp b/components/omega/test/ocn/VertCoordTest.cpp
index 827c288f0e4a..3417d75a75ef 100644
--- a/components/omega/test/ocn/VertCoordTest.cpp
+++ b/components/omega/test/ocn/VertCoordTest.cpp
@@ -25,6 +25,7 @@
 #include "mpi.h"
 
 #include <algorithm>
+#include <random>
 
 using namespace OMEGA;
 
@@ -60,14 +61,11 @@ int initVertCoordTest() {
    if (Err != 0)
       LOG_ERROR("VertCoordTest: error initializing default halo");
 
-   // Begin initialization of the default vertical coordinate
-   VertCoord::init1();
-
    // Initialize the default mesh
    HorzMesh::init();
 
-   // Complete initialization of the default vertical coordinate
-   VertCoord::init2();
+   // Initialize the default vertical coordinate
+   VertCoord::init();
 
    return Err;
 } // end initVertCoordTest
@@ -586,6 +584,95 @@ int main(int argc, char *argv[]) {
              Error(ErrorCode::Fail, "VertCoordTest: minMaxLayerVertex FAIL");
       }
 
+      // Tests for Masks
+
+      /// Setup random values for MinLayerCell and MaxLayerCell near the top
+      /// and bottom of each column
+      I4 Seed = 12345;
+      std::mt19937 Generator(Seed);
+
+      I4 MinMin = 0;
+      I4 MaxMin = 10;
+      I4 MinMax = NVertLayers - 11;
+      I4 MaxMax = NVertLayers - 1;
+
+      std::uniform_int_distribution<int> MinDist(MinMin, MaxMin);
+      std::uniform_int_distribution<int> MaxDist(MinMax, MaxMax);
+
+      for (int ICell = 0; ICell < NCellsAll; ICell++) {
+         DefVertCoord->MinLayerCellH(ICell) = MinDist(Generator);
+         DefVertCoord->MaxLayerCellH(ICell) = MaxDist(Generator);
+      }
+
+      deepCopy(DefVertCoord->MinLayerCell, DefVertCoord->MinLayerCellH);
+      deepCopy(DefVertCoord->MaxLayerCell, DefVertCoord->MaxLayerCellH);
+      Kokkos::fence();
+
+      /// Reset min/max layer arrays and set masks
+      DefVertCoord->minMaxLayerEdge();
+      DefVertCoord->minMaxLayerVertex();
+      DefVertCoord->setMasks();
+
+      Err = 0;
+      for (int ICell = 0; ICell < NCellsAll; ++ICell) {
+         /// The sum of the masks in a column is the number of active
+         /// layers in that column
+         Real Expected = DefVertCoord->MaxLayerCellH(ICell) -
+                         DefVertCoord->MinLayerCellH(ICell) + 1._Real;
+         Real Sum = 0.;
+         for (int K = 0; K < NVertLayers; ++K) {
+            Sum += DefVertCoord->CellMaskH(ICell, K);
+         }
+
+         Real Diff = std::abs(Sum - Expected);
+         if (Diff > 1e-10) {
+            Err += 1;
+         }
+      }
+
+      for (int IEdge = 0; IEdge < NEdgesAll; ++IEdge) {
+         Real Expected;
+
+         /// The sum of the masks along an edge is the number of layers with
+         /// active cells on both sides.
+         if ((DefMesh->CellsOnEdgeH(IEdge, 1) == NCellsAll) ||
+             (DefMesh->CellsOnEdgeH(IEdge, 0) == NCellsAll)) {
+            Expected = 0.;
+         } else {
+            Expected = DefVertCoord->MaxLayerEdgeTopH(IEdge) -
+                       DefVertCoord->MinLayerEdgeBotH(IEdge) + 1._Real;
+         }
+         Real Sum = 0.;
+         for (int K = 0; K < NVertLayers; ++K) {
+            Sum += DefVertCoord->EdgeMaskH(IEdge, K);
+         }
+         Real Diff = std::abs(Sum - Expected);
+         if (Diff > 1e-10) {
+            Err += 1;
+         }
+      }
+
+      for (int IVertex = 0; IVertex < NVerticesAll; ++IVertex) {
+         Real Expected = DefVertCoord->MaxLayerVertexBotH(IVertex) -
+                         DefVertCoord->MinLayerVertexTopH(IVertex) + 1._Real;
+
+         Real Sum = 0.;
+         for (int K = 0; K < NVertLayers; ++K) {
+            Sum += DefVertCoord->VertexMaskH(IVertex, K);
+         }
+         Real Diff = std::abs(Sum - Expected);
+         if (Diff > 1e-10) {
+            Err += 1;
+         }
+      }
+
+      /// Determine test pass/fail
+      if (Err == 0) {
+         LOG_INFO("VertCoordTest: setMasks PASS");
+      } else {
+         ErrAll += Error(ErrorCode::Fail, "VertCoordTest: setMasks FAIL");
+      }
+
       // Finalize Omega objects
       IOStream::finalize();
       TimeStepper::clear();
diff --git a/components/omega/test/timeStepping/TimeStepperTest.cpp b/components/omega/test/timeStepping/TimeStepperTest.cpp
index 7a280d0bcf1e..53b3dc1cb96c 100644
--- a/components/omega/test/timeStepping/TimeStepperTest.cpp
+++ b/components/omega/test/timeStepping/TimeStepperTest.cpp
@@ -37,7 +37,6 @@
 
 #include <cmath>
 #include <iomanip>
-#include <iostream>
 
 using namespace OMEGA;
 
@@ -170,15 +169,13 @@ int initTimeStepperTest(const std::string &mesh) {
       LOG_ERROR("TimeStepperTest: error initializing default halo");
    }
 
-   // Initialize the vertical coordinate and reset NVertLayers to 1
-   VertCoord::init1();
-   auto *DefVertCoord        = VertCoord::getDefault();
-   DefVertCoord->NVertLayers = 1;
-   Dimension::destroy("NVertLayers");
-   std::shared_ptr<Dimension> VertDim =
-       Dimension::create("NVertLayers", NVertLayers);
-
    HorzMesh::init();
+
+   // initialize vertical coordinate, do not read stream and use local
+   // NVertLayers value
+   VertCoord::init(false, NVertLayers);
+   auto *DefVertCoord = VertCoord::getDefault();
+
    Tracers::init();
    AuxiliaryState::init();
    Tendencies::init();
@@ -207,8 +204,8 @@ int initTimeStepperTest(const std::string &mesh) {
       LOG_ERROR("TimeStepperTest: error creating test state");
    }
 
-   auto *TestAuxState = AuxiliaryState::create("TestAuxState", DefMesh, DefHalo,
-                                               NVertLayers, NTracers);
+   auto *TestAuxState = AuxiliaryState::create(
+       "TestAuxState", DefMesh, DefVertCoord, DefHalo, NVertLayers, NTracers);
 
    Config *OmegaConfig = Config::getOmegaConfig();
    TestAuxState->readConfigOptions(OmegaConfig);
@@ -282,10 +279,9 @@ void finalizeTimeStepperTest() {
    AuxiliaryState::clear();
    OceanState::clear();
    VertCoord::clear();
+   HorzMesh::clear();
    Dimension::clear();
    Field::clear();
-   HorzMesh::clear();
-   VertCoord::clear();
    Halo::clear();
    Decomp::clear();
    MachEnv::removeAll();