diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml
index 94aa1a3b3daf..efc04e521057 100644
--- a/components/omega/configs/Default.yml
+++ b/components/omega/configs/Default.yml
@@ -66,6 +66,7 @@ Omega:
     VelocityVertAdvTendencyEnable: true
     TracerVertAdvTendencyEnable: true
     PressureGradTendencyEnable: true
+    FrazilTendencyEnable: true
     VelVertMixTendencyEnable: true
     TracerVertMixTendencyEnable: true
   ManufacturedSolution:
@@ -81,6 +82,12 @@ Omega:
       DRhoDT: -0.2
       DRhoDS: 0.8
       RhoT0S0: 1000.0
+  Frazil:
+    FrazilType: basic
+    MassLimit: 0.1
+    Phi: 0.75
+    DepthLimit: -1.0
+    ConservationCheck: false
   VertMix:
     Background:
       Diffusivity: 1.0e-5
diff --git a/components/omega/doc/devGuide/Frazil.md b/components/omega/doc/devGuide/Frazil.md
new file mode 100644
index 000000000000..ee13c3dd99ce
--- /dev/null
+++ b/components/omega/doc/devGuide/Frazil.md
@@ -0,0 +1,101 @@
+(omega-dev-frazil)=
+
+# Frazil
+
+This page describes frazil design and implementation details in Omega,
+including both the `basic` and `teos` pathways.
+
+## Purpose and coupling points
+
+Frazil computes phase-change-related tendencies that modify:
+
+- pseudo-thickness tendency
+- temperature tracer tendency
+- salinity tracer tendency
+
+The tendency hook-up is implemented through `FrazilOnCell` in the tracer
+tendency phase, where frazil contributions are added to the accumulated
+`PseudoThicknessTend` and `TracerTend` arrays.
+
+## Data flow and call sequence
+
+1. `Tendencies::computeTracerTendenciesOnly` checks
+   `Tendencies.FrazilTendencyEnable`.
+2. If enabled, `FrazilOnCell::operator()` retrieves the default `Frazil`
+   object and zeros frazil tendency/accumulator arrays.
+3. `FrazilOnCell` extracts `Temperature` and `Salinity` tracer subviews,
+   then calls `Frazil::computeFrazil(CT, SA, PressureMid, PseudoThickness)`.
+4. `Frazil::computeFrazil` dispatches to `computeFrazilBasicImpl` or
+   `computeFrazilTeosImpl` based on `FrazilType`.
+5. Returned frazil tendencies are added into `PseudoThicknessTend` and
+   `TracerTend` for all active cell layers.
+
+## Configuration coupling
+
+Frazil behavior is configured with:
+
+- `Omega.Tendencies.FrazilTendencyEnable`
+  - global switch for applying frazil tendency terms
+- `Omega.Frazil.FrazilType`
+  - implementation choice (`basic` or `teos`)
+- `Omega.Frazil.MassLimit`
+  - per-layer mass/thickness limiter used by frazil formation/melt pathways
+- `Omega.Frazil.Phi`
+  - teos pathway liquid-fraction parameter for new frazil partitioning
+- `Omega.Frazil.DepthLimit`
+  - optional depth cutoff for frazil activity; negative means no cutoff
+- `Omega.Frazil.ConservationCheck`
+  - optional post-compute column conservation diagnostic/logging
+
+## Physics/algorithm summary
+
+### Common behavior
+
+- Freezing-point checks are based on conservative temperature and absolute
+  salinity with pressure-dependent freezing temperature.
+- Vertical accumulation order is bottom-to-top within each active column.
+- Frazil tendencies are not time-step scaled inside `Frazil`; they are
+  accumulated as tendency contributions.
+
+### Basic pathway (`FrazilType: basic`)
+
+- Uses simplified energetics: the energy of the super-cooled water sets the
+amount of solid ice formed (used constant latent heat of fusion of fresh ice).
+Salt is added based on a constant bulk salinity `IceRefSal` (default) or a
+manual toggle (for now) using the local salinity and the frazil porosity.
+Melting of existing frazil is set by the amount of pure ice that can be melted
+by the warm layer, and the enthalpy of melt water at the local freezing point
+is added to the temperature tendency.
+  temperature comparisons and fractional thickness limits.
+- Computes local layer tendencies (`HTend`, `TTend`, `STend`) and updates
+  accumulated frazil stores.
+- Applies surface salt redistribution adjustment and converts accumulators to
+  coupler units at the end of the column loop.
+  -Warnings: 1) basic frazil formation/melt is does not conserve energy.
+  2) Using porosity to set the salt content includes a redistribution of excess
+  salt in the surface layer, which can be very significant.
+
+### Teos pathway (`FrazilType: teos`)
+
+- Uses teos-10 Gibbs SeaWater routines for frazil formation and melt state
+  transitions.
+- Formation uses `Phi` and `MassLimit` to partition and limit newly formed
+  frazil contributions.
+- Melt computes fraction melted subject to available thermodynamic energy and
+  mass-limit constraints.
+
+## Existing ctest coverage
+
+The existing frazil test driver is in
+`components/omega/test/ocn/FrazilTest.cpp` and covers:
+
+- teos frazil formation in cold and warm single-layer states
+- basic frazil formation in cold and warm single-layer states
+- mixed warm/cold column behavior with sign checks for branch switching
+- depth-limit behavior ensuring excluded deep layers have zero frazil tendency
+
+
+## Related pages
+
+- User-facing options: [User Frazil Guide](../userGuide/Frazil.md)
+- Tendency container/hook-up: [Tendencies](Tendencies.md)
diff --git a/components/omega/doc/devGuide/Tendencies.md b/components/omega/doc/devGuide/Tendencies.md
index 443e8e6ddcd5..6f2e8d2e594c 100644
--- a/components/omega/doc/devGuide/Tendencies.md
+++ b/components/omega/doc/devGuide/Tendencies.md
@@ -64,6 +64,23 @@ Tendencies.computeVelocityTendencies(State, AuxState, ThickTimeLevel, VelTimeLev
 To call only the tracer tendency terms:
 Tendencies.computeTracerTendencies(State, AuxState, TracerArray, ThickTimeLevel, VelTimeLevel);
 
+## Frazil tendencies
+
+Frazil tendencies are calculated by the Frazil class, not the TendencyTerms.
+The tendencies are integrated in the tracer tendency compute path and read in from
+the frazil computation. When
+`Omega.Tendencies.FrazilTendencyEnable` is true, `computeTracerTendenciesOnly`
+invokes `FrazilOnCell`, which calls the default `Frazil` object and adds frazil
+contributions to:
+
+- `PseudoThicknessTend`
+- `TracerTend` for `Temperature`
+- `TracerTend` for `Salinity`
+
+The frazil implementation choice and parameters are read from the
+`Omega.Frazil` block. Detailed algorithm/physics notes are in
+[Frazil](Frazil.md).
+
 ## Removal of tendencies
 To erase a specific named tendencies instance use `erase`
 ```c++
diff --git a/components/omega/doc/index.md b/components/omega/doc/index.md
index 5cb09351f8cf..8b93d2cb593c 100644
--- a/components/omega/doc/index.md
+++ b/components/omega/doc/index.md
@@ -54,6 +54,7 @@ userGuide/Timing
 userGuide/VerticalMixingCoeff
 userGuide/VertAdv
 userGuide/Forcing
+userGuide/Frazil
 userGuide/SfcCoupling
 ```
 
@@ -101,6 +102,7 @@ devGuide/Timing
 devGuide/VerticalMixingCoeff
 devGuide/VertAdv
 devGuide/Forcing
+devGuide/Frazil
 devGuide/SfcCoupling
 ```
 
diff --git a/components/omega/doc/userGuide/Frazil.md b/components/omega/doc/userGuide/Frazil.md
new file mode 100644
index 000000000000..1109cefe31dc
--- /dev/null
+++ b/components/omega/doc/userGuide/Frazil.md
@@ -0,0 +1,67 @@
+(omega-user-frazil)=
+
+# Frazil
+
+This page describes user-facing configuration for frazil tendencies in Omega.
+Frazil physics is used to represent the formation and melt of frazil ice within
+the ocean water column. It impacts the local layer pseudo-thickness (i.e. mass),
+temperature, and salinity tendencies. The vertical sum of the frazil energy,
+mass of water and mass of salt are passed to the coupler (if coupled) or
+discarded (in ocean standalone mode).
+
+## Configuration overview
+
+Frazil behavior is controlled by one enable switch in `Tendencies` and one
+`Frazil` configuration block:
+
+```yaml
+Omega:
+  Tendencies:
+    FrazilTendencyEnable: true
+
+  Frazil:
+    FrazilType: basic
+    MassLimit: 0.1
+    Phi: 0.75
+    DepthLimit: -1.0
+    ConservationCheck: false
+```
+
+- `Tendencies.FrazilTendencyEnable`
+  - enables/disables application of frazil tendency contributions
+- `Frazil.FrazilType`
+  - selects frazil option
+  - supported options in current code: `basic` and `teos`
+  - `simple` exists as a placeholder name but is not supported
+- `Frazil.MassLimit`
+  - limits per-layer frazil mass/thickness tendency magnitude (applied to formation and melt)
+- `Frazil.Phi`
+  - liquid-mass fraction of frazil (used by the teos frazil formation, or by the basic frazil when using porosity rather than constant salinity).
+- `Frazil.DepthLimit`
+  - limits the depth range over which frazil is computed
+  - negative values mean no depth limit
+- `Frazil.ConservationCheck`
+  - enables a column-level diagnostic conservation check with logging
+
+## Available frazil options
+
+Omega currently includes two active frazil pathways:
+
+- `basic`
+  - freezing is based on the formation of fresh solid ice, to which salt is added (similar the mpas-ocean implementation).
+- `teos`
+  - teos-10-based option using Gibbs SeaWater thermodynamic routines.
+
+Both pathways contribute to:
+
+- pseudo-thickness tendency
+- temperature tracer tendency
+- salinity tracer tendency
+
+## Notes
+
+- Frazil tendencies are applied through the `Tendencies` tracer-step workflow.
+- The frazil tendency hook assumes tracer names include `Temperature` and
+  `Salinity`.
+- For implementation and algorithm details, see
+  [Developer Frazil Guide](../devGuide/Frazil.md).
diff --git a/components/omega/doc/userGuide/Tendencies.md b/components/omega/doc/userGuide/Tendencies.md
index 6fc1803b176d..11cca7c4ca25 100644
--- a/components/omega/doc/userGuide/Tendencies.md
+++ b/components/omega/doc/userGuide/Tendencies.md
@@ -6,3 +6,9 @@ The `Tendencies` class provides a container for the [tendency terms](#omega-user
 Upon creation of an `Tendencies` instance, these functors are initialized and arrays for the
 accumulated tendencies are allocated.
 There are no user-configurable options beyond those for the tendency term functors.
+
+Frazil is one of the configurable tendency contributions. It is controlled by
+`Omega.Tendencies.FrazilTendencyEnable` together with the `Omega.Frazil`
+configuration block (`FrazilType`, `MassLimit`, `Phi`, `DepthLimit`, and
+`ConservationCheck`). For operational guidance and options, see
+[Frazil](Frazil.md).
diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h
index a530549d493d..effaefec3ada 100644
--- a/components/omega/src/ocn/Eos.h
+++ b/components/omega/src/ocn/Eos.h
@@ -354,8 +354,8 @@ class Teos10Eos {
 
    /// Calculates freezing Conservative Temperature using TEOS-10 polynomial
    /// (polynomial error in [-5e-4, 6e-4] K, from GSW package)
-   KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P,
-                                       const Real SaturationFract) const {
+   KOKKOS_FUNCTION static Real calcCtFreezing(const Real Sa, const Real P,
+                                              const Real SaturationFract) {
       constexpr Real Sso = 35.16504;
       constexpr Real C0  = 0.017947064327968736;
       constexpr Real C1  = -6.076099099929818;
@@ -756,6 +756,20 @@ class Eos {
    /// Convert potential temperature to Conservative Temperature
    Real calcCtFromPt(const Real &Sa, const Real &Pt) const;
 
+   /// Calculate Conservative Temperature at freezing point.
+   /// P is expected in dbar to match TEOS polynomial convention.
+   KOKKOS_FUNCTION static Real calcCtFreezing(const Real Sa, const Real P,
+                                              const Real SaturationFract,
+                                              const EosType Choice) {
+      if (Choice == EosType::Teos10Eos) {
+         return Teos10Eos::calcCtFreezing(Sa, P, SaturationFract);
+      }
+
+      Kokkos::abort("Eos::calcCtFreezing: CtFreezing not implemented for "
+                    "non-TEOS-10 EOS");
+      return 0.0_Real;
+   }
+
    /// Initialize EOS from config and mesh
    static void init();
 
diff --git a/components/omega/src/ocn/Frazil.cpp b/components/omega/src/ocn/Frazil.cpp
new file mode 100644
index 000000000000..9f062a7db42c
--- /dev/null
+++ b/components/omega/src/ocn/Frazil.cpp
@@ -0,0 +1,535 @@
+//===-- ocn/Frazil.cpp - Frazil Ice Formation -----------------*- C++ -*-===//
+//
+// The Frazil class manages frazil tendencies and accumulators.
+// This initial implementation only has a teos-10 configuration.
+//
+//===----------------------------------------------------------------------===//
+
+#include "Frazil.h"
+#include "Eos.h"
+#include "Error.h"
+#include "Logging.h"
+
+namespace OMEGA {
+
+namespace {
+
+KOKKOS_INLINE_FUNCTION bool isApprox(const Real X, const Real Y,
+                                     const Real RTol, const Real ATol = 0) {
+   if (Kokkos::isnan(X) || Kokkos::isnan(Y) || Kokkos::isinf(X) ||
+       Kokkos::isinf(Y)) {
+      return false;
+   }
+
+   return Kokkos::abs(X - Y) <=
+          Kokkos::max(ATol, RTol * Kokkos::max(Kokkos::abs(X), Kokkos::abs(Y)));
+}
+
+} // namespace
+
+Frazil *Frazil::DefaultFrazil = nullptr;
+std::map<std::string, std::unique_ptr<Frazil>> Frazil::AllFrazil;
+
+/// Constructor for FrazilFormation
+FrazilFormation::FrazilFormation() {}
+
+/// Constructor for FrazilMelt
+FrazilMelt::FrazilMelt() {}
+
+/// Constructor for BasicFrazilFormation
+BasicFrazilFormation::BasicFrazilFormation() {}
+
+/// Constructor for BasicFrazilMelt
+BasicFrazilMelt::BasicFrazilMelt() {}
+
+void Frazil::init() {
+
+   if (!HorzMesh::getDefault() or !VertCoord::getDefault()) {
+      ABORT_ERROR("Frazil::init: HorzMesh and VertCoord must be initialized");
+   }
+
+   // Frazil freezing-temperature calculations depend on EOS configuration.
+   Eos::init();
+
+   if (!DefaultFrazil) {
+      Error Err;
+      bool FrazilTendencyEnable = false;
+      Config *OmegaConfig       = Config::getOmegaConfig();
+      Config TendConfig("Tendencies");
+
+      Err += OmegaConfig->get(TendConfig);
+      CHECK_ERROR_ABORT(Err,
+                        "Frazil::init: Tendencies group not found in Config");
+
+      Err += TendConfig.get("FrazilTendencyEnable", FrazilTendencyEnable);
+      CHECK_ERROR_ABORT(
+          Err, "Frazil::init: FrazilTendencyEnable not found in Tendencies");
+
+      if (!FrazilTendencyEnable) {
+         LOG_INFO("Frazil::init: Frazil tendency disabled; skipping default "
+                  "frazil object creation");
+         LOG_INFO("All frazil is off - frazil parameters will be ignored");
+         return;
+      }
+
+      DefaultFrazil = create("Default");
+   }
+}
+
+Frazil::Frazil(const HorzMesh *Mesh, const VertCoord *VCoord)
+    : frazilChoice(FrazilType::TeosFrazil), computeFrazilFormation(),
+      computeFrazilMelt(), NCellsAll(Mesh->NCellsAll),
+      NChunks((VCoord->NVertLayers + VecLength - 1) / VecLength), MeshPtr(Mesh),
+      VCoordPtr(VCoord) {
+
+   FrazilTTend =
+       Array2DReal("FrazilTTend", Mesh->NCellsSize, VCoord->NVertLayers);
+   FrazilSTend =
+       Array2DReal("FrazilSTend", Mesh->NCellsSize, VCoord->NVertLayers);
+   FrazilHTend =
+       Array2DReal("FrazilHTend", Mesh->NCellsSize, VCoord->NVertLayers);
+
+   AccMIce  = Array1DReal("AccMIce", Mesh->NCellsSize);
+   AccEIce  = Array1DReal("AccEIce", Mesh->NCellsSize);
+   AccMLiq  = Array1DReal("AccMLiq", Mesh->NCellsSize);
+   AccELiq  = Array1DReal("AccELiq", Mesh->NCellsSize);
+   AccMSalt = Array1DReal("AccMSalt", Mesh->NCellsSize);
+
+   deepCopy(FrazilTTend, 0.0_Real);
+   deepCopy(FrazilSTend, 0.0_Real);
+   deepCopy(FrazilHTend, 0.0_Real);
+   deepCopy(AccMIce, 0.0_Real);
+   deepCopy(AccEIce, 0.0_Real);
+   deepCopy(AccMLiq, 0.0_Real);
+   deepCopy(AccELiq, 0.0_Real);
+   deepCopy(AccMSalt, 0.0_Real);
+}
+
+Frazil::~Frazil() {}
+
+Frazil *Frazil::create(const std::string &Name) {
+   if (AllFrazil.find(Name) != AllFrazil.end()) {
+      LOG_ERROR("Attempted to create Frazil {} but it already exists", Name);
+      return nullptr;
+   }
+
+   auto *NewFrazil =
+       new Frazil(HorzMesh::getDefault(), VertCoord::getDefault());
+   AllFrazil.emplace(Name, NewFrazil);
+
+   Error Err;
+   Config *OmegaConfig = Config::getOmegaConfig();
+   Config FrazilConfig("Frazil");
+   Err += OmegaConfig->get(FrazilConfig);
+   CHECK_ERROR_ABORT(Err, "Frazil::create: Frazil group not found in Config");
+
+   std::string FrazilTypeStr;
+   Err += FrazilConfig.get("FrazilType", FrazilTypeStr);
+   CHECK_ERROR_ABORT(Err,
+                     "Frazil::create: FrazilType not found in Frazil config");
+
+   if ((FrazilTypeStr == "Basic") or (FrazilTypeStr == "basic") or
+       (FrazilTypeStr == "BasicFrazil")) {
+      NewFrazil->frazilChoice = FrazilType::BasicFrazil;
+   } else if ((FrazilTypeStr == "Simple") or (FrazilTypeStr == "simple") or
+              (FrazilTypeStr == "SimpleFrazil")) {
+      NewFrazil->frazilChoice = FrazilType::SimpleFrazil;
+      ABORT_ERROR("Frazil::create: SimpleFrazil not supported yet");
+   } else if ((FrazilTypeStr == "teos") or (FrazilTypeStr == "Teos") or
+              (FrazilTypeStr == "TEOS") or (FrazilTypeStr == "Teos10") or
+              (FrazilTypeStr == "teos10") or (FrazilTypeStr == "TEOS10")) {
+      NewFrazil->frazilChoice = FrazilType::TeosFrazil;
+   } else {
+      ABORT_ERROR("Frazil::create: Unknown FrazilType requested");
+   }
+
+   Err += FrazilConfig.get("MassLimit",
+                           NewFrazil->computeFrazilFormation.massLimit);
+   Err += FrazilConfig.get("MassLimit", NewFrazil->computeFrazilMelt.massLimit);
+   CHECK_ERROR_ABORT(Err,
+                     "Frazil::create: MassLimit not found in Frazil config");
+
+   NewFrazil->computeBasicFrazilFormation.massLimit =
+       NewFrazil->computeFrazilFormation.massLimit;
+   NewFrazil->computeBasicFrazilMelt.massLimit =
+       NewFrazil->computeFrazilFormation.massLimit;
+
+   Err += FrazilConfig.get("Phi", NewFrazil->computeFrazilFormation.phi);
+   CHECK_ERROR_ABORT(Err, "Frazil::create: Phi not found in Frazil config");
+
+   Err += FrazilConfig.get("ConservationCheck", NewFrazil->conservationCheck);
+   CHECK_ERROR_ABORT(
+       Err, "Frazil::create: ConservationCheck not found in Frazil config");
+
+   Err += FrazilConfig.get("DepthLimit", NewFrazil->depthLimit);
+   CHECK_ERROR_ABORT(Err,
+                     "Frazil::create: DepthLimit not found in Frazil config");
+
+   if (Name == "Default") {
+      DefaultFrazil = NewFrazil;
+   }
+
+   return NewFrazil;
+}
+
+Frazil *Frazil::getDefault() { return DefaultFrazil; }
+
+Frazil *Frazil::get(const std::string &Name) {
+   auto it = AllFrazil.find(Name);
+   if (it != AllFrazil.end()) {
+      return it->second.get();
+   }
+
+   LOG_ERROR("Frazil::get: Attempted to retrieve non-existent Frazil {}", Name);
+   return nullptr;
+}
+
+void Frazil::erase(std::string InName) {
+   auto *ToErase = get(InName);
+   AllFrazil.erase(InName);
+   if (ToErase == DefaultFrazil) {
+      DefaultFrazil = nullptr;
+   }
+}
+
+void Frazil::clear() {
+   AllFrazil.clear();
+   DefaultFrazil = nullptr;
+}
+
+void Frazil::checkColumnConservation() const {
+   auto MinLayerCellH = createHostMirrorCopy(VCoordPtr->MinLayerCell);
+   auto MaxLayerCellH = createHostMirrorCopy(VCoordPtr->MaxLayerCell);
+   auto FrazilHTendH  = createHostMirrorCopy(FrazilHTend);
+   auto FrazilTTendH  = createHostMirrorCopy(FrazilTTend);
+   auto FrazilSTendH  = createHostMirrorCopy(FrazilSTend);
+   auto AccMIceH      = createHostMirrorCopy(AccMIce);
+   auto AccMLiqH      = createHostMirrorCopy(AccMLiq);
+   auto AccMSaltH     = createHostMirrorCopy(AccMSalt);
+   auto AccELiqH      = createHostMirrorCopy(AccELiq);
+   auto AccEIceH      = createHostMirrorCopy(AccEIce);
+
+   constexpr Real RTol = 1.0e-10_Real;
+
+   for (I4 ICell = 0; ICell < NCellsAll; ++ICell) {
+      const I4 KMin = MinLayerCellH(ICell);
+      const I4 KMax = MaxLayerCellH(ICell);
+
+      Real MassTend   = 0.0_Real;
+      Real EnergyTend = 0.0_Real;
+      Real SaltTend   = 0.0_Real;
+
+      for (I4 K = KMin; K <= KMax; ++K) {
+         MassTend += FrazilHTendH(ICell, K);
+         EnergyTend += FrazilTTendH(ICell, K);
+         SaltTend += FrazilSTendH(ICell, K);
+      }
+
+      const Real MassTotal   = AccMIceH(ICell) + AccMLiqH(ICell);
+      const Real EnergyTotal = AccELiqH(ICell) + AccEIceH(ICell);
+      const Real SaltTotal   = AccMSaltH(ICell);
+      if (ICell == 0) {
+         LOG_INFO("Frazil column conservation check: cell {} MassTend={} "
+                  "MassTotal={} "
+                  "EnergyTend={} EnergyTotal={} SaltTend={} SaltTotal={}",
+                  ICell, MassTend * RhoSw, MassTotal,
+                  EnergyTend * Cp0Sw * RhoSw, EnergyTotal,
+                  SaltTend * RhoSw * PPt2Salt, SaltTotal);
+         LOG_INFO("Frazil column conservation check: cell {} EpsMass={} "
+                  "EpsE={} EpsS={} ",
+                  ICell, MassTend * RhoSw + MassTotal,
+                  EnergyTend * Cp0Sw * RhoSw + EnergyTotal,
+                  SaltTend * RhoSw * PPt2Salt + SaltTotal);
+      }
+
+      if (!isApprox(-MassTend * RhoSw, MassTotal, RTol)) {
+         LOG_INFO(
+             "Frazil column mass check failed: cell {} tendency={} total={}",
+             ICell, -MassTend * RhoSw, MassTotal);
+      }
+      if (!isApprox(-EnergyTend * Cp0Sw * RhoSw, EnergyTotal, RTol)) {
+         LOG_INFO(
+             "Frazil column energy check failed: cell {} tendency={} total={}",
+             ICell, -EnergyTend * Cp0Sw * RhoSw, EnergyTotal);
+      }
+      if (!isApprox(-SaltTend * RhoSw * PPt2Salt, SaltTotal, RTol)) {
+         LOG_INFO(
+             "Frazil column salt check failed: cell {} tendency={} total={}",
+             ICell, -SaltTend * RhoSw * PPt2Salt, SaltTotal);
+      }
+   }
+}
+
+void Frazil::computeFrazilBasicImpl(const Array2DReal &CT,
+                                    const Array2DReal &SA, const Array2DReal &P,
+                                    const Array2DReal &LayerH) {
+   const EosType LocEosChoice = Eos::getInstance()->EosChoice;
+   const Real LocDepthLimit   = depthLimit;
+
+   OMEGA_SCOPE(MinLayerCell, VCoordPtr->MinLayerCell);
+   OMEGA_SCOPE(MaxLayerCell, VCoordPtr->MaxLayerCell);
+   OMEGA_SCOPE(LocGeomZMid, VCoordPtr->GeomZMid);
+
+   OMEGA_SCOPE(LocComputeBasicFrazilFormation, computeBasicFrazilFormation);
+   OMEGA_SCOPE(LocComputeBasicFrazilMelt, computeBasicFrazilMelt);
+   OMEGA_SCOPE(LocFrazilTTend, FrazilTTend);
+   OMEGA_SCOPE(LocFrazilSTend, FrazilSTend);
+   OMEGA_SCOPE(LocFrazilHTend, FrazilHTend);
+   OMEGA_SCOPE(LocAccMIce, AccMIce);
+   OMEGA_SCOPE(LocAccEIce, AccEIce);
+   OMEGA_SCOPE(LocAccMLiq, AccMLiq);
+   OMEGA_SCOPE(LocAccELiq, AccELiq);
+   OMEGA_SCOPE(LocAccMSalt, AccMSalt);
+   OMEGA_SCOPE(LocIceRefSal, IceRefSal);
+   OMEGA_SCOPE(LocLatIce, LatIce);
+
+   parallelFor(
+       {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) {
+          const I4 KMin = MinLayerCell(ICell);
+          const I4 KMax = MaxLayerCell(ICell);
+
+          I4 Klim          = KMax;
+          bool HasKlim     = true;
+          const bool Limit = (LocDepthLimit >= 0.0_Real);
+
+          if (Limit) {
+             HasKlim = false;
+             for (I4 K = KMax; K >= KMin; --K) {
+                if (Kokkos::abs(LocGeomZMid(ICell, K)) <= LocDepthLimit) {
+                   Klim    = K;
+                   HasKlim = true;
+                   break;
+                }
+             }
+          }
+
+          // Explicit accumulation order: bottom layer to top layer.
+          for (I4 K = KMax; K >= KMin; --K) {
+             if (!HasKlim || K > Klim) {
+                LocFrazilHTend(ICell, K) = 0.0_Real;
+                LocFrazilTTend(ICell, K) = 0.0_Real;
+                LocFrazilSTend(ICell, K) = 0.0_Real;
+                continue;
+             }
+
+             const Real SAIn = SA(ICell, K);
+             const Real CTIn = CT(ICell, K);
+             const Real PIn  = P(ICell, K);
+             const Real PDb  = PIn * Pa2Db;
+             const Real H    = LayerH(ICell, K);
+
+             const Real Tfrz =
+                 Eos::calcCtFreezing(SAIn, PDb, 0.0_Real, LocEosChoice);
+
+             Real HTend = 0.0_Real;
+             Real TTend = 0.0_Real;
+             Real STend = 0.0_Real;
+
+             if (CTIn < Tfrz) {
+                LocComputeBasicFrazilFormation(
+                    SAIn, CTIn, PDb, H, LocAccMIce(ICell), LocAccMSalt(ICell),
+                    LocAccEIce(ICell), HTend, TTend, STend, Tfrz);
+             } else if (LocAccMIce(ICell) > 0.0_Real) {
+                LocComputeBasicFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell),
+                                          LocAccMSalt(ICell), LocAccEIce(ICell),
+                                          HTend, TTend, STend, Tfrz);
+             }
+
+             // temporary log -- TBRemoved
+             // if  (ICell == 0) {
+             // LOG_INFO("computeFrazil cell = {}, SAIn = {}, CTIn = {}, PIn = "
+             //          "{}, H = {}, Tfrz = {}",
+             //          ICell, SAIn, CTIn, PIn, H, Tfrz);
+             // LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} "
+             //          "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}",
+             //          ICell, K, (CTIn < Tfrz), LocAccMIce(ICell),
+             //          LocAccMLiq(ICell), LocAccMSalt(ICell),
+             //          LocAccELiq(ICell), LocAccEIce(ICell));
+             // LOG_INFO("                                     HTend= {} "
+             //          "TTend= {} STend= {}",
+             //          HTend, TTend, STend);
+             //}
+
+             LocFrazilHTend(ICell, K) = HTend; // not scaled by dt
+             LocFrazilTTend(ICell, K) = TTend;
+             LocFrazilSTend(ICell, K) = STend;
+          } // end of vertical loop
+
+          //  // temporary log -- TBRemoved
+          //  if (ICell == 0) {
+          //     LOG_INFO(
+          //         "Frazil::computeFrazil: (basic={}), cell={} AccMIce={} "
+          //         "AccMSaltCalc={} AccMSaltCpl={} "
+          //         "AccEIceCalc={} AccEIceCpl={}",
+          //         true, ICell, LocAccMIce(ICell),
+          //         LocAccMSalt(ICell), LocAccMIce(ICell) * LocIceRefSal,
+          //         LocAccEIce(ICell), -LocAccMIce(ICell) * LocLatIce);
+          //  }
+
+          // Redistribute excess salt at the surface. No treatment of low
+          // salinity frazil for now.
+          LocFrazilSTend(ICell, KMin) += Kokkos::max(
+              0.0_Real, LocAccMIce(ICell) * LocIceRefSal - LocAccMSalt(ICell));
+          // hijack total terms before the coupling
+          LocAccMSalt(ICell) = LocAccMIce(ICell) * LocIceRefSal;
+          LocAccEIce(ICell)  = -LocAccMIce(ICell) * LocLatIce;
+
+          // Convert to coupler units
+          LocAccMIce(ICell)  = LocAccMIce(ICell) * RhoSw;
+          LocAccMLiq(ICell)  = LocAccMLiq(ICell) * RhoSw;
+          LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt;
+          LocAccELiq(ICell)  = LocAccELiq(ICell) * RhoSw;
+          LocAccEIce(ICell)  = LocAccEIce(ICell) * RhoSw;
+       }); // end of NCells loop
+}
+
+void Frazil::computeFrazilTeosImpl(const Array2DReal &CT, const Array2DReal &SA,
+                                   const Array2DReal &P,
+                                   const Array2DReal &LayerH) {
+   const EosType LocEosChoice = Eos::getInstance()->EosChoice;
+   const Real LocDepthLimit   = depthLimit;
+
+   OMEGA_SCOPE(MinLayerCell, VCoordPtr->MinLayerCell);
+   OMEGA_SCOPE(MaxLayerCell, VCoordPtr->MaxLayerCell);
+   OMEGA_SCOPE(LocGeomZMid, VCoordPtr->GeomZMid);
+
+   OMEGA_SCOPE(LocComputeFrazilFormation, computeFrazilFormation);
+   OMEGA_SCOPE(LocComputeFrazilMelt, computeFrazilMelt);
+   OMEGA_SCOPE(LocFrazilTTend, FrazilTTend);
+   OMEGA_SCOPE(LocFrazilSTend, FrazilSTend);
+   OMEGA_SCOPE(LocFrazilHTend, FrazilHTend);
+   OMEGA_SCOPE(LocAccMIce, AccMIce);
+   OMEGA_SCOPE(LocAccEIce, AccEIce);
+   OMEGA_SCOPE(LocAccMLiq, AccMLiq);
+   OMEGA_SCOPE(LocAccELiq, AccELiq);
+   OMEGA_SCOPE(LocAccMSalt, AccMSalt);
+
+   parallelFor(
+       {NCellsAll}, KOKKOS_LAMBDA(I4 ICell) {
+          const I4 KMin = MinLayerCell(ICell);
+          const I4 KMax = MaxLayerCell(ICell);
+
+          I4 Klim          = KMax;
+          bool HasKlim     = true;
+          const bool Limit = (LocDepthLimit >= 0.0_Real);
+
+          // calculates the depth limit based on geometric height
+          if (Limit) {
+             HasKlim = false;
+             for (I4 K = KMax; K >= KMin; --K) {
+                if (Kokkos::abs(LocGeomZMid(ICell, K)) <= LocDepthLimit) {
+                   Klim    = K;
+                   HasKlim = true;
+                   break;
+                }
+             }
+          }
+
+          // Explicit accumulation order: bottom layer to top layer.
+          for (I4 K = KMax; K >= KMin; --K) {
+             if (!HasKlim || K > Klim) {
+                LocFrazilHTend(ICell, K) = 0.0_Real;
+                LocFrazilTTend(ICell, K) = 0.0_Real;
+                LocFrazilSTend(ICell, K) = 0.0_Real;
+                continue;
+             }
+
+             const Real SAIn = SA(ICell, K);
+             const Real CTIn = CT(ICell, K);
+             const Real PIn  = P(ICell, K);
+             const Real PDb  = PIn * Pa2Db;
+             const Real H    = LayerH(ICell, K);
+
+             const Real Tfrz =
+                 Eos::calcCtFreezing(SAIn, PDb, 0.0_Real, LocEosChoice);
+
+             Real HTend = 0.0_Real;
+             Real TTend = 0.0_Real;
+             Real STend = 0.0_Real;
+
+             if (CTIn < Tfrz) {
+                LocComputeFrazilFormation(SAIn, CTIn, PDb, H, LocAccMIce(ICell),
+                                          LocAccMLiq(ICell), LocAccMSalt(ICell),
+                                          LocAccELiq(ICell), LocAccEIce(ICell),
+                                          HTend, TTend, STend);
+             } else if (LocAccMIce(ICell) > 0.0_Real) {
+                LocComputeFrazilMelt(SAIn, CTIn, PDb, H, LocAccMIce(ICell),
+                                     LocAccMLiq(ICell), LocAccMSalt(ICell),
+                                     LocAccELiq(ICell), LocAccEIce(ICell),
+                                     HTend, TTend, STend);
+             }
+
+             //  // temporary kernel logging -- TBRemoved
+             //  // used for debugging - only on CPUs.
+             //  // can be removed once we have single-column testing
+             //  if (ICell == 0) {
+             //     const Real Hf  = H + HTend;
+             //     const Real SAf = (H * SAIn + STend) / Hf;
+             //     const Real Tf  = (H * CTIn + TTend) / Hf;
+             //     // LOG_INFO("computeFrazil cell = {}, SAIn = {}, CTIn = {},
+             //     PIn
+             //     // = "
+             //     //          "{}, H = {}, Tfrz = {}",
+             //     //          ICell, SAIn, CTIn, PIn, H, Tfrz);
+             //     LOG_INFO("computeFrazil cell={} K={} (cold={}) AccMIce={} "
+             //              "AccMLiq={} AccMSalt={} AccELiq={} AccEIce={}",
+             //              ICell, K, (CTIn < Tfrz), LocAccMIce(ICell),
+             //              LocAccMLiq(ICell), LocAccMSalt(ICell),
+             //              LocAccELiq(ICell), LocAccEIce(ICell));
+             //     LOG_INFO("                                     HTend= {} "
+             //              "TTend= {} STend= {}",
+             //              HTend, TTend, STend);
+             //     LOG_INFO("                                    "
+             //              "H= {}-->{}, T= {}-->{}, S= {}-->{}",
+             //              H, Hf, CTIn, Tf, SAIn, SAf);
+             //  }
+
+             LocFrazilHTend(ICell, K) = HTend; // not scaled by dt
+             LocFrazilTTend(ICell, K) = TTend;
+             LocFrazilSTend(ICell, K) = STend;
+          } // end of vertical loop
+
+          // Convert to coupler units
+          LocAccMIce(ICell)  = LocAccMIce(ICell) * RhoSw;
+          LocAccMLiq(ICell)  = LocAccMLiq(ICell) * RhoSw;
+          LocAccMSalt(ICell) = LocAccMSalt(ICell) * RhoSw * PPt2Salt;
+          LocAccELiq(ICell)  = LocAccELiq(ICell) * RhoSw;
+          LocAccEIce(ICell)  = LocAccEIce(ICell) * RhoSw;
+       }); // end of NCells loop
+}
+
+void Frazil::computeFrazil(const Array2DReal &CT, const Array2DReal &SA,
+                           const Array2DReal &P, const Array2DReal &LayerH) {
+   Eos *DefEos = Eos::getInstance();
+   if (!DefEos) {
+      ABORT_ERROR("Frazil::computeFrazil: Eos must be initialized before "
+                  "computeFrazil");
+   }
+
+   const EosType LocEosChoice = DefEos->EosChoice;
+   if (LocEosChoice != EosType::Teos10Eos) {
+      ABORT_ERROR("Frazil::computeFrazil: CtFreezing not implemented for "
+                  "non-TEOS-10 EOS");
+   }
+
+   switch (frazilChoice) {
+   case FrazilType::BasicFrazil:
+      computeFrazilBasicImpl(CT, SA, P, LayerH);
+      break;
+   case FrazilType::TeosFrazil:
+      computeFrazilTeosImpl(CT, SA, P, LayerH);
+      break;
+   case FrazilType::SimpleFrazil:
+      ABORT_ERROR("Frazil::computeFrazil: SimpleFrazil not supported yet");
+      break;
+   default:
+      ABORT_ERROR("Frazil::computeFrazil: Unknown frazilChoice");
+      break;
+   }
+
+   if (conservationCheck) {
+      checkColumnConservation();
+   }
+} // end of computeFrazil
+
+} // namespace OMEGA
diff --git a/components/omega/src/ocn/Frazil.h b/components/omega/src/ocn/Frazil.h
new file mode 100644
index 000000000000..8eb6f4a9102a
--- /dev/null
+++ b/components/omega/src/ocn/Frazil.h
@@ -0,0 +1,378 @@
+#ifndef OMEGA_FRAZIL_H
+#define OMEGA_FRAZIL_H
+//===-- ocn/Frazil.h - Frazil Ice Formation -------------------*- C++ -*-===//
+//
+// The Frazil class manages frazil tendencies and accumulators.
+// This initial implementation only has a teos-10 configuration.
+// but carries scaffolding for other implementations.
+//
+//===----------------------------------------------------------------------===//
+
+#include "Config.h"
+// #include "DataTypes.h"
+#include "GlobalConstants.h"
+#include "HorzMesh.h"
+#include "OmegaKokkos.h"
+#include "VertCoord.h"
+
+#include <map>
+#include <memory>
+#include <string>
+
+#include <gswteos-10.h>
+
+namespace OMEGA {
+
+enum class FrazilType {
+   BasicFrazil,  ///< MPAS-O style basic frazil option
+   SimpleFrazil, ///< Placeholder simple frazil option
+   TeosFrazil    ///< Placeholder TEOS frazil option
+};
+
+class BasicFrazilFormation {
+ public:
+   BasicFrazilFormation();
+
+   Real massLimit         = 0.1_Real;  // to do:  remove default
+   Real FrazilIceSalinity = IceRefSal; // Global constant
+   Real LatFrazil         = LatIce;    // Global constant
+
+   Real FrazilPorosity =
+       1.0_Real; // Internal for now; can move to config later.
+
+   KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real PDb,
+                                   const Real H, Real &SumIceThickness,
+                                   Real &SumSalt, Real &SumEnergy, Real &HTend,
+                                   Real &TTend, Real &STend,
+                                   const Real Tfrz) const {
+
+      const Real potential      = H * Cp0Sw * RhoSw * (CT - Tfrz);
+      const Real freezingEnergy = Kokkos::max(0.0_Real, -potential);
+
+      HTend = 0.0_Real;
+      TTend = 0.0_Real;
+      STend = 0.0_Real;
+
+      Real newFrzThickness =
+          freezingEnergy /
+          (LatFrazil * RhoSw); // frazil (ice) mass in pseudo-thickness terms
+
+      newFrzThickness   = Kokkos::min(newFrzThickness, H * massLimit);
+      Real newFrzEnergy = -newFrzThickness * LatFrazil; // (<0; enthalpy of ice)
+
+      // MANUAL TOGGLE: uncomment line below to use porosity
+      // const Real FrazilIceSalinity = FrazilPorosity * SA;
+
+      const Real frazilSalinity = Kokkos::min(FrazilIceSalinity, SA);
+      const Real newSaltContent =
+          newFrzThickness * frazilSalinity; // in m.(g/kg)
+
+      HTend = -newFrzThickness - newSaltContent;
+      // TTend below should include the enthalpy associated with the mass flux
+      // this is intentionally not added here to match the mpas-o implementation
+      TTend =
+          -(newFrzEnergy) / (Cp0Sw); // (E< 0 so TTend>0) // scaled to h.CT tend
+      STend = -newSaltContent;
+
+      SumIceThickness += newFrzThickness;
+      // non-conservation between Sum and HTend by construction in the original
+      SumSalt += newSaltContent;
+      SumEnergy += newFrzEnergy;
+   }
+};
+
+// TO-DO: is there a new for salt reconciliation at the surface?
+// frazilSalinityTendency(minLevelCell(iCell),iCell) =
+// frazilSalinityTendency(minLevelCell(iCell),iCell) + &
+//     max(0.0_RKIND,(sumNewThicknessWeightedSaltContent -
+//     newThicknessWeightedSaltContent) ) / dt
+// accumulatedFrazilIceSalinityNew(iCell) =
+// accumulatedFrazilIceSalinityOld(iCell) + newThicknessWeightedSaltContent
+
+class BasicFrazilMelt {
+ public:
+   BasicFrazilMelt();
+
+   Real massLimit         = 0.1_Real;  // to do:  remove default
+   Real FrazilIceSalinity = IceRefSal; // Global constant
+   Real LatFrazil         = LatIce;    // Global constant
+
+   KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real PDb,
+                                   const Real H, Real &SumIceThickness,
+                                   Real &SumSalt, Real &SumEnergy, Real &HTend,
+                                   Real &TTend, Real &STend,
+                                   const Real Tfrz) const {
+      constexpr Real Eps = 1.0e-12_Real;
+
+      if (SumIceThickness <= Eps) { // potential leak if we dont redistribute
+         SumIceThickness = 0.0_Real;
+         SumSalt         = 0.0_Real;
+         HTend           = 0.0_Real;
+         TTend           = 0.0_Real;
+         STend           = 0.0_Real;
+         return;
+      }
+      const Real potential       = H * Cp0Sw * RhoSw * (CT - Tfrz);
+      const Real availableEnergy = Kokkos::max(0.0_Real, potential);
+
+      HTend = 0.0_Real;
+      TTend = 0.0_Real;
+      STend = 0.0_Real;
+
+      Real meltThickness =
+          availableEnergy /
+          (LatFrazil * RhoSw); // mass in pseudo-thickness units
+      meltThickness = Kokkos::min(meltThickness, SumIceThickness);
+      meltThickness = Kokkos::min(meltThickness,
+                                  H * massLimit); // also 0.1h lim on added mass
+      const Real meltAverageSalinity = SumSalt / SumIceThickness;
+      const Real meltingEnergy =
+          meltThickness * LatFrazil; // (energy needed to melt)
+
+      HTend = +meltThickness * (1 + meltAverageSalinity); // (>0 so HTend>0)
+      TTend =
+          -meltingEnergy /
+              (Cp0Sw) + // (TTend < 0 when melting for phase change)
+          meltThickness *
+              Tfrz; // (added the enthalpy of the melt water at freezing temp)
+      STend = +meltAverageSalinity * meltThickness; // (STend > 0 when melting)
+
+      SumIceThickness -= meltThickness;
+      SumSalt -= meltAverageSalinity * meltThickness;
+      SumEnergy -=
+          meltingEnergy; // necessarily non-conservative by construction
+   }
+};
+
+class FrazilMelt {
+ public:
+   /// constructor declaration
+   FrazilMelt();
+
+   // masslimit parameter (set in config)
+   Real massLimit;
+
+   //   The functor for FrazilMelt takes as inputs:
+   //   the local ocean layer state (SA, CT, P, h),
+   //   the accumulated frazil solid and liquid mass, energy, and salt
+   //   and outputs the frazil tendencies (HTend, TTend, STend) and updated
+   //   accumulators.
+   KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real P,
+                                   const Real h, Real &AccMIce, Real &AccMLiq,
+                                   Real &AccMSalt, Real &AccELiq, Real &AccEIce,
+                                   Real &HTend, Real &TTend,
+                                   Real &STend) const {
+
+      constexpr Real Eps = 1.0e-12_Real;
+
+      // this check on AccMIce and E should be done in the calling function, but
+      // is here for safety we can do a better implementation of the checks
+      if (AccMIce <= Eps ||
+          AccEIce >= Eps) { // potential leak if we dont redistribute
+         AccMIce = 0.0_Real;
+         AccEIce = 0.0_Real;
+         HTend   = 0.0_Real;
+         TTend   = 0.0_Real;
+         STend   = 0.0_Real;
+         return;
+      }
+
+      if (AccMLiq <= Eps ||
+          AccELiq >= Eps) { // potential leak if we dont redistribute
+         AccMLiq = 0.0_Real;
+         AccELiq = 0.0_Real;
+         HTend   = 0.0_Real;
+         TTend   = 0.0_Real;
+         STend   = 0.0_Real;
+         return;
+      }
+
+      // 1. we start by adding the solid ice to the ocean layer (no brine yet)
+      const Real potEnthalpyIce = AccEIce / AccMIce;
+      const Real newLayerMass =
+          Kokkos::max(h + AccMIce, Eps); // max unnecessary but for safety
+      const Real newLayerIceFraction = AccMIce / newLayerMass;
+
+      // 2. we calculate the (mass- and energy-conserving) ocean layer evolution
+      // ... how much (pure) ice does this layer melt?
+
+      // typecasting for now but will be simplified once gsw functions are
+      // ported
+      const double SA_d     = static_cast<double>(SA);
+      const double CT_d     = static_cast<double>(CT);
+      const double P_d      = static_cast<double>(P);
+      const double wIhIn_d  = static_cast<double>(newLayerIceFraction);
+      const double pt0Ice_d = gsw_pt_from_pot_enthalpy_ice_poly(
+          static_cast<double>(potEnthalpyIce));
+      const double tIce_d = gsw_t_from_pt0_ice(pt0Ice_d, P_d);
+      double SAnew_d      = SA_d;
+      double CTnew_d      = CT_d;
+      double wIhOut_d     = wIhIn_d;
+
+      gsw_melting_ice_into_seawater(SA_d, CT_d, P_d, wIhIn_d, tIce_d, &SAnew_d,
+                                    &CTnew_d, &wIhOut_d);
+
+      const Real wIhOut = static_cast<Real>(
+          wIhOut_d); // by def 0<= wIhOut <= 1 ;  To-do: check function behavior
+
+      // 3. we calculate the mass fraction of the frazil (pure) ice that was
+      // melted - limited by a total mass limit of 0.1h
+      const Real solidMassMelted = Kokkos::max(
+          0.0_Real,
+          AccMIce - wIhOut * newLayerMass); // original - left-over solid ice,
+      const Real frazilFractionMelted = Kokkos::min(
+          solidMassMelted / AccMIce,
+          h * massLimit / (AccMIce + AccMLiq)); // total added mass < 0.1h
+
+      // the frazil fraction based on the solid ice also sets the (proportional)
+      // contributions from the frazil brine
+      HTend = +(frazilFractionMelted * (AccMLiq + AccMIce));
+      TTend = +(frazilFractionMelted * (AccELiq + AccEIce)) / Cp0Sw;
+      STend = +(frazilFractionMelted * AccMSalt);
+
+      const Real frazilFractionLeft =
+          Kokkos::max(0.0_Real, 1.0_Real - frazilFractionMelted);
+
+      AccMIce  = frazilFractionLeft * AccMIce;
+      AccMLiq  = frazilFractionLeft * AccMLiq;
+      AccMSalt = frazilFractionLeft * AccMSalt;
+      AccELiq  = frazilFractionLeft * AccELiq;
+      AccEIce  = frazilFractionLeft * AccEIce;
+   }
+};
+
+class FrazilFormation {
+ public:
+   /// constructor declaration
+   FrazilFormation();
+
+   /// Parameters for FrazilFormation (set by yaml file)
+   Real phi;       ///< liquid mass fraction of new frazil (0 < Phi < 1)
+   Real massLimit; ///< layer mass fraction limit for thickness tendency
+
+   //   The functor for FrazilFormation takes as inputs:
+   //   the local ocean layer state (SA, CT, P, h),
+   //   the accumulated frazil solid and liquid mass, energy, and salt
+   //   and outputs the frazil tendencies (HTend, TTend, STend) and updated
+   //   accumulators.
+   KOKKOS_FUNCTION void operator()(const Real SA, const Real CT, const Real P,
+                                   const Real h, Real &AccMIce, Real &AccMLiq,
+                                   Real &AccMSalt, Real &AccELiq, Real &AccEIce,
+                                   Real &HTend, Real &TTend,
+                                   Real &STend) const {
+
+      Real CTnew;
+      Real SAnew;
+      Real wIh            = 0.0_Real;
+      Real solidMass      = 0.0_Real;
+      Real liquidMass     = 0.0_Real;
+      Real solidEnthalpy  = 0.0_Real;
+      Real liquidEnthalpy = 0.0_Real;
+
+      // type casting for now because gsw expects double
+      // and Real can be either float or double depending on build configuration
+      double SAnew_d = 0.0;
+      double CTnew_d = 0.0;
+      double wIh_d   = 0.0;
+      // double PTnew_d = 0.0;
+
+      gsw_frazil_properties_potential_poly(
+          static_cast<double>(SA), static_cast<double>(Cp0Sw * CT),
+          static_cast<double>(P), &SAnew_d, &CTnew_d, &wIh_d);
+      double PTnew_d =
+          gsw_pt_from_ct(SAnew_d, CTnew_d); // convert to potential temperature
+
+      SAnew = static_cast<Real>(SAnew_d);
+      CTnew = static_cast<Real>(CTnew_d);
+      wIh   = static_cast<Real>(wIh_d);
+
+      const Real oneMinusPhi = Kokkos::max(1.0e-12_Real, 1.0_Real - phi);
+      // anything called mass below is in pseudo-thickness units (m) and needs
+      // to be scaled by RhoSw for coupling
+      solidMass      = h * Kokkos::min(wIh, oneMinusPhi * massLimit);
+      liquidMass     = (phi / oneMinusPhi) * solidMass;
+      solidEnthalpy  = solidMass * gsw_pot_enthalpy_from_pt_ice_poly(PTnew_d);
+      liquidEnthalpy = liquidMass * Cp0Sw * CTnew;
+      // per timestep (not scaled by dt here)
+      HTend = -(solidMass + liquidMass);
+      TTend = -(liquidEnthalpy + solidEnthalpy) / Cp0Sw;
+      STend = -(liquidMass * SAnew);
+
+      // Local unit of mass is pseudo thickness (m)
+      // these all need a RhoSw factor before coupling
+      AccMIce += solidMass;           // m
+      AccMLiq += liquidMass;          // m
+      AccMSalt += liquidMass * SAnew; // (m)(g/kg)
+      AccELiq += liquidEnthalpy;      // (m)(J/kg)
+      AccEIce += solidEnthalpy;       // (m)(J/kg)
+   }
+};
+
+class Frazil {
+ public:
+   static void init();
+   /// Creates a new frazil object and stores it in the AllFrazil map.
+   static Frazil *create(const std::string &Name);
+
+   /// Retrieve frazil object by name.
+   static Frazil *get(const std::string &Name);
+
+   /// Retrieve default frazil object.
+   static Frazil *getDefault();
+
+   /// Destructor
+   ~Frazil();
+
+   /// Deallocates arrays
+   static void clear();
+
+   /// Remove frazil object by name.
+   static void erase(std::string InName); ///< [in] name to remove
+
+   Array2DReal FrazilTTend;
+   Array2DReal FrazilSTend;
+   Array2DReal FrazilHTend;
+   Array1DReal AccMIce;
+   Array1DReal AccEIce;
+   Array1DReal AccMLiq;
+   Array1DReal AccELiq;
+   Array1DReal AccMSalt;
+
+   void computeFrazil(const Array2DReal &CT, const Array2DReal &SA,
+                      const Array2DReal &P, const Array2DReal &H);
+   void computeFrazilBasicImpl(const Array2DReal &CT, const Array2DReal &SA,
+                               const Array2DReal &P, const Array2DReal &LayerH);
+   void computeFrazilTeosImpl(const Array2DReal &CT, const Array2DReal &SA,
+                              const Array2DReal &P, const Array2DReal &LayerH);
+   bool conservationCheck = false;
+   Real depthLimit        = -1.0_Real;
+
+ private:
+   static Frazil *DefaultFrazil;
+   static std::map<std::string, std::unique_ptr<Frazil>> AllFrazil;
+
+   Frazil(const HorzMesh *Mesh, const VertCoord *VCoord);
+
+   // Forbid copy and move construction/assignment.
+   Frazil(const Frazil &)            = delete;
+   Frazil &operator=(const Frazil &) = delete;
+   Frazil(Frazil &&)                 = delete;
+   Frazil &operator=(Frazil &&)      = delete;
+
+   FrazilType frazilChoice;
+   BasicFrazilFormation computeBasicFrazilFormation;
+   BasicFrazilMelt computeBasicFrazilMelt;
+   FrazilFormation computeFrazilFormation;
+   FrazilMelt computeFrazilMelt;
+   I4 NCellsAll;
+   I4 NChunks;
+
+   const HorzMesh *MeshPtr;
+   const VertCoord *VCoordPtr;
+
+   void checkColumnConservation() const;
+};
+
+} // namespace OMEGA
+
+#endif
diff --git a/components/omega/src/ocn/OceanFinal.cpp b/components/omega/src/ocn/OceanFinal.cpp
index 1ae7061bd79f..edd131dfce5c 100644
--- a/components/omega/src/ocn/OceanFinal.cpp
+++ b/components/omega/src/ocn/OceanFinal.cpp
@@ -10,6 +10,7 @@
 #include "Eos.h"
 #include "Field.h"
 #include "Forcing.h"
+#include "Frazil.h"
 #include "Halo.h"
 #include "HorzMesh.h"
 #include "IO.h"
@@ -48,6 +49,7 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time
    AuxiliaryState::clear();
    Forcing::clear();
    OceanState::clear();
+   Frazil::clear();
    VertAdv::clear();
    VertCoord::clear();
    Dimension::clear();
@@ -61,7 +63,6 @@ int ocnFinalize(const TimeInstant &CurrTime ///< [in] current sim time
    // destroy singletons that use raw new/delete
    Eos::destroyInstance();
    VertMix::destroyInstance();
-
    return RetVal;
 } // end ocnFinalize
 
diff --git a/components/omega/src/ocn/OceanInit.cpp b/components/omega/src/ocn/OceanInit.cpp
index d1c7d5d227e9..633cf9726c28 100644
--- a/components/omega/src/ocn/OceanInit.cpp
+++ b/components/omega/src/ocn/OceanInit.cpp
@@ -15,6 +15,7 @@
 #include "Error.h"
 #include "Field.h"
 #include "Forcing.h"
+#include "Frazil.h"
 #include "Halo.h"
 #include "HorzMesh.h"
 #include "IO.h"
@@ -206,6 +207,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) {
    Forcing::init();
    AuxiliaryState::init();
    Eos::init();
+   Frazil::init();
    PressureGrad::init();
    VertMix::init();
    Tendencies::init();
diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp
index a85e24d37dc6..447fa24aaf5d 100644
--- a/components/omega/src/ocn/Tendencies.cpp
+++ b/components/omega/src/ocn/Tendencies.cpp
@@ -266,6 +266,10 @@ void Tendencies::readConfig(Config *OmegaConfig ///< [in] Omega config
    CHECK_ERROR_ABORT(
        Err, "Tendencies: PressureGradTendencyEnable not found in TendConfig");
 
+   Err += TendConfig.get("FrazilTendencyEnable", this->FrazilTerm.Enabled);
+   CHECK_ERROR_ABORT(
+       Err, "Tendencies: FrazilTendencyEnable not found in TendConfig");
+
    Err += TendConfig.get("SurfaceTracerRestoringEnable",
                          this->SurfaceTracerRestoring.Enabled);
    CHECK_ERROR_ABORT(
@@ -420,7 +424,7 @@ Tendencies::Tendencies(const std::string &Name_, ///< [in] Name for tendencies
       SfcStressForcing(Mesh, VCoord), BottomDrag(Mesh, VCoord),
       TracerDiffusion(Mesh, VCoord), TracerHyperDiff(Mesh, VCoord),
       TracerHorzAdv(Mesh, VCoord), SurfaceTracerRestoring(Mesh),
-      CustomThicknessTend(InCustomThicknessTend),
+      FrazilTerm(Mesh, VCoord), CustomThicknessTend(InCustomThicknessTend),
       CustomVelocityTend(InCustomVelocityTend), EqState(EqState), PGrad(PGrad),
       VMix(VMix) {
 
@@ -850,6 +854,15 @@ void Tendencies::computeTracerTendenciesOnly(
       Pacer::stop("Tend:surfaceTracerRestoring", 2);
    }
 
+   if (FrazilTerm.Enabled) {
+      Pacer::start("Tend:frazil", 2);
+      const auto &PressureMid     = VCoord->PressureMid;
+      Array2DReal PseudoThickness = State->getPseudoThickness(ThickTimeLevel);
+      FrazilTerm(PseudoThicknessTend, TracerTend, TracerArray, PressureMid,
+                 PseudoThickness);
+      Pacer::stop("Tend:frazil", 2);
+   }
+
    Pacer::stop("Tend:computeTracerTendenciesOnly", 1);
 } // end tracer tendency compute
 
diff --git a/components/omega/src/ocn/Tendencies.h b/components/omega/src/ocn/Tendencies.h
index ac826f915a85..20915d143058 100644
--- a/components/omega/src/ocn/Tendencies.h
+++ b/components/omega/src/ocn/Tendencies.h
@@ -73,6 +73,7 @@ class Tendencies {
    TracerDiffOnCell TracerDiffusion;
    TracerHyperDiffOnCell TracerHyperDiff;
    SurfaceTracerRestoringOnCell SurfaceTracerRestoring;
+   FrazilOnCell FrazilTerm;
 
    std::string Name;
 
diff --git a/components/omega/src/ocn/TendencyTerms.cpp b/components/omega/src/ocn/TendencyTerms.cpp
index 5b142888142c..a79e1c5523ef 100644
--- a/components/omega/src/ocn/TendencyTerms.cpp
+++ b/components/omega/src/ocn/TendencyTerms.cpp
@@ -11,6 +11,7 @@
 #include "TendencyTerms.h"
 #include "AuxiliaryState.h"
 #include "DataTypes.h"
+#include "Error.h"
 #include "HorzMesh.h"
 #include "HorzOperators.h"
 #include "OceanState.h"
@@ -113,6 +114,73 @@ TracerHyperDiffOnCell::TracerHyperDiffOnCell(const HorzMesh *Mesh,
 SurfaceTracerRestoringOnCell::SurfaceTracerRestoringOnCell(
     const HorzMesh *Mesh) {}
 
+FrazilOnCell::FrazilOnCell(const HorzMesh *Mesh, const VertCoord *VCoord)
+    : NCellsAll(Mesh->NCellsAll), TempTracerIndex(-1), SaltTracerIndex(-1),
+      MinLayerCell(VCoord->MinLayerCell), MaxLayerCell(VCoord->MaxLayerCell) {
+   Tracers::getIndex(TempTracerIndex, "Temperature");
+   Tracers::getIndex(SaltTracerIndex, "Salinity");
+
+   OMEGA_REQUIRE(TempTracerIndex >= 0,
+                 "FrazilOnCell: Temperature tracer index is undefined");
+   OMEGA_REQUIRE(SaltTracerIndex >= 0,
+                 "FrazilOnCell: Salinity tracer index is undefined");
+}
+
+void FrazilOnCell::operator()(const Array2DReal &PseudoThicknessTend,
+                              const Array3DReal &TracerTend,
+                              const Array3DReal &TracerArray,
+                              const Array2DReal &PressureMid,
+                              const Array2DReal &PseudoThickness) const {
+   auto *Frazil = Frazil::getDefault();
+   if (!Enabled || !Frazil) {
+      return;
+   }
+
+   deepCopy(Frazil->FrazilTTend, 0.0_Real);
+   deepCopy(Frazil->FrazilSTend, 0.0_Real);
+   deepCopy(Frazil->FrazilHTend, 0.0_Real);
+   deepCopy(Frazil->AccMIce, 0.0_Real);
+   deepCopy(Frazil->AccEIce, 0.0_Real);
+   deepCopy(Frazil->AccMLiq, 0.0_Real);
+   deepCopy(Frazil->AccELiq, 0.0_Real);
+   deepCopy(Frazil->AccMSalt, 0.0_Real);
+
+   const auto ConservTemp =
+       Kokkos::subview(TracerArray, TempTracerIndex, Kokkos::ALL, Kokkos::ALL);
+   const auto AbsSalinity =
+       Kokkos::subview(TracerArray, SaltTracerIndex, Kokkos::ALL, Kokkos::ALL);
+
+   Frazil->computeFrazil(ConservTemp, AbsSalinity, PressureMid,
+                         PseudoThickness);
+
+   const auto FrazilHTend = Frazil->FrazilHTend;
+   const auto FrazilTTend = Frazil->FrazilTTend;
+   const auto FrazilSTend = Frazil->FrazilSTend;
+   const I4 TempIndex     = TempTracerIndex;
+   const I4 SaltIndex     = SaltTracerIndex;
+
+   OMEGA_SCOPE(LocPseudoThicknessTend, PseudoThicknessTend);
+   OMEGA_SCOPE(LocTracerTend, TracerTend);
+   OMEGA_SCOPE(LocFrazilHTend, FrazilHTend);
+   OMEGA_SCOPE(LocFrazilTTend, FrazilTTend);
+   OMEGA_SCOPE(LocFrazilSTend, FrazilSTend);
+   OMEGA_SCOPE(LocMinLayerCell, MinLayerCell);
+   OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell);
+
+   parallelForOuter(
+       {NCellsAll}, KOKKOS_LAMBDA(int ICell, const TeamMember &Team) {
+          const int KMin = LocMinLayerCell(ICell);
+          const int KMax = LocMaxLayerCell(ICell);
+
+          parallelForInner(
+              Team, Range{KMin, KMax}, INNER_LAMBDA(int K) {
+                 LocPseudoThicknessTend(ICell, K) += LocFrazilHTend(ICell, K);
+                 LocTracerTend(TempIndex, ICell, K) += LocFrazilTTend(ICell, K);
+                 LocTracerTend(SaltIndex, ICell, K) += LocFrazilSTend(ICell, K);
+              });
+       });
+}
+
 void TracerHorzAdvOnCell::init() {
    const HorzMesh *Mesh    = this->HorzontalMesh;
    const VertCoord *VCoord = this->VerticalCoord;
diff --git a/components/omega/src/ocn/TendencyTerms.h b/components/omega/src/ocn/TendencyTerms.h
index 389d6bb1aafb..0d46101ec398 100644
--- a/components/omega/src/ocn/TendencyTerms.h
+++ b/components/omega/src/ocn/TendencyTerms.h
@@ -11,6 +11,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "AuxiliaryState.h"
+#include "Frazil.h"
 #include "GlobalConstants.h"
 #include "HorzMesh.h"
 #include "MachEnv.h"
@@ -599,5 +600,26 @@ class SurfaceTracerRestoringOnCell {
    }
 };
 
+/// Frazil tendency hook-up term
+class FrazilOnCell {
+ public:
+   bool Enabled = false;
+
+   FrazilOnCell(const HorzMesh *Mesh, const VertCoord *VCoord);
+
+   void operator()(const Array2DReal &PseudoThicknessTend,
+                   const Array3DReal &TracerTend,
+                   const Array3DReal &TracerArray,
+                   const Array2DReal &PressureMid,
+                   const Array2DReal &PseudoThickness) const;
+
+ private:
+   I4 NCellsAll;
+   I4 TempTracerIndex;
+   I4 SaltTracerIndex;
+   Array1DI4 MinLayerCell;
+   Array1DI4 MaxLayerCell;
+};
+
 } // namespace OMEGA
 #endif
diff --git a/components/omega/test/CMakeLists.txt b/components/omega/test/CMakeLists.txt
index 76df1e94a435..341897e02156 100644
--- a/components/omega/test/CMakeLists.txt
+++ b/components/omega/test/CMakeLists.txt
@@ -512,6 +512,13 @@ add_omega_test(
     "-n;2"
 )
 
+add_omega_test(
+  FRAZIL_TEST
+  testFrazil.exe
+  ocn/FrazilTest.cpp
+  "-n;2"
+)
+
 ################
 # VertCoord test
 ################
diff --git a/components/omega/test/infra/IOStreamTest.cpp b/components/omega/test/infra/IOStreamTest.cpp
index 7d8b37a50a38..475d242ff00c 100644
--- a/components/omega/test/infra/IOStreamTest.cpp
+++ b/components/omega/test/infra/IOStreamTest.cpp
@@ -116,12 +116,13 @@ void initIOStreamTest(Clock *&ModelClock // Model clock
    // Initialize VertMix
    VertMix::init();
 
+   // Initialize Tracers before Tendencies so tracer indices are available
+   // during FrazilOnCell construction inside Tendencies::init().
+   Tracers::init();
+
    // Intialize Tendencies
    Tendencies::init();
 
-   // Initialize Tracers
-   Tracers::init();
-
    // IOStream::validateAll() depends on Forcing::init() so Forcing fields
    // are registered before stream validation.
    Forcing::init();
diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp
new file mode 100644
index 000000000000..9120e681d002
--- /dev/null
+++ b/components/omega/test/ocn/FrazilTest.cpp
@@ -0,0 +1,642 @@
+//===-- Test driver for OMEGA FrazilFormation ---------------------------*- C++
+//-*-===//
+//
+/// \file
+/// \brief Minimal test driver for OMEGA frazil formation functor
+//
+//===-----------------------------------------------------------------------===/
+
+#include "Frazil.h"
+#include "Config.h"
+#include "DataTypes.h"
+#include "Decomp.h"
+#include "Dimension.h"
+#include "Field.h"
+#include "Halo.h"
+#include "HorzMesh.h"
+#include "IO.h"
+#include "IOStream.h"
+#include "Logging.h"
+#include "MachEnv.h"
+#include "OceanTestCommon.h"
+#include "OmegaKokkos.h"
+#include "Pacer.h"
+#include "TimeMgr.h"
+#include "VertCoord.h"
+#include "mpi.h"
+
+using namespace OMEGA;
+
+constexpr int NVertLayers = 60;
+
+void initFrazilTest(const std::string &mesh) {
+   MachEnv::init(MPI_COMM_WORLD);
+   MachEnv *DefEnv  = MachEnv::getDefault();
+   MPI_Comm DefComm = DefEnv->getComm();
+
+   initLogging(DefEnv);
+   LOG_INFO("------ Frazil Unit Tests ------");
+
+   Config("Omega");
+   Config::readAll("omega.yml");
+
+   Calendar::init("No Leap");
+   TimeInstant StartTime(0, 1, 1, 0, 0, 0.0);
+   TimeInterval TimeStep(1, TimeUnits::Hours);
+   Clock ModelClockTmp(StartTime, TimeStep);
+   Clock *ModelClock = &ModelClockTmp;
+
+   IO::init(DefComm);
+   Decomp::init(mesh);
+   Field::init(ModelClock);
+   IOStream::init(ModelClock);
+   Halo::init();
+   HorzMesh::init(ModelClock);
+   VertCoord::init(false);
+   Frazil::init();
+}
+
+void finalizeFrazilTest() {
+   Frazil::clear();
+   VertCoord::clear();
+   HorzMesh::clear();
+   Halo::clear();
+   Decomp::clear();
+   Field::clear();
+   Dimension::clear();
+   IOStream::finalize();
+   MachEnv::removeAll();
+}
+
+// this test only excercises the frazil formation functor (no melt)
+// in a cold case: the frazil terms should be
+// - strictly positive for ice, liquid, and salt mass
+//- strictly negative for ice and liquid energy
+// - positive for T tendency and negative for S, H tendencies
+void testFrazilFormationCold() {
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
+
+   VCoord->NVertLayers = NVertLayers;
+
+   const Real SAIn = 35.0_Real;
+   const Real CTIn = -2.0_Real;
+   const Real PIn  = 100.0_Real;
+   const Real h    = 10.0_Real;
+   const Real RTol = 1e-10_Real;
+
+   (void)Mesh;
+
+   FrazilFormation ComputeFrazilFormation;
+   ComputeFrazilFormation.phi       = 0.75_Real;
+   ComputeFrazilFormation.massLimit = 0.1_Real;
+
+   Real AccMIce  = 0.0_Real;
+   Real AccMLiq  = 0.0_Real;
+   Real AccMSalt = 0.0_Real;
+   Real AccELiq  = 0.0_Real;
+   Real AccEIce  = 0.0_Real;
+
+   Real HTend = 0.0_Real;
+   Real TTend = 0.0_Real;
+   Real STend = 0.0_Real;
+
+   ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMLiq, AccMSalt,
+                          AccELiq, AccEIce, HTend, TTend, STend);
+
+   if (AccMIce <= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilFormationTestCold: accumulated ice mass is non-positive: {}",
+          AccMIce);
+   }
+   if (AccMLiq <= 0.0_Real) {
+      ABORT_ERROR("FrazilFormationTestCold: accumulated liquid mass is "
+                  "non-positive: {}",
+                  AccMLiq);
+   }
+   if (AccMSalt <= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilFormationTestCold: accumulated salt mass is non-positive: {}",
+          AccMSalt);
+   }
+
+   if (AccELiq >= 0.0_Real) {
+      ABORT_ERROR("FrazilFormationTestCold: accumulated liquid energy is "
+                  "positive (exp. negative): {}",
+                  AccELiq);
+   }
+   if (AccEIce >= 0.0_Real) {
+      ABORT_ERROR("FrazilFormationTestCold: accumulated ice energy is positive "
+                  "(exp. negative): {}",
+                  AccEIce);
+   }
+   if (HTend >= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilFormationTestCold: HTend is positive (exp. negative): {}",
+          HTend);
+   }
+   if (TTend <= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilFormationTestCold: TTend is negative (exp. positive): {}",
+          TTend);
+   }
+   if (STend >= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilFormationTestCold: STend is positive (exp. negative): {}",
+          STend);
+   }
+   LOG_INFO(
+       "FrazilFormationTestCold: AccMIce = {}, AccMLiq = {}, AccMSalt = {}, "
+       "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}",
+       AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend);
+}
+
+// this test only excercises the frazil formation functor (no melt)
+// in a warm case: the frazil FORMATION terms should all be zero
+void testFrazilFormationWarm() {
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
+
+   VCoord->NVertLayers = NVertLayers;
+
+   const Real SAIn = 35.0_Real;
+   const Real CTIn = 10.0_Real;
+   const Real PIn  = 100.0_Real;
+   const Real h    = 10.0_Real;
+   const Real RTol = 1e-10_Real;
+
+   (void)Mesh;
+
+   FrazilFormation ComputeFrazilFormation;
+   ComputeFrazilFormation.phi       = 0.75_Real;
+   ComputeFrazilFormation.massLimit = 0.1_Real;
+
+   Real AccMIce  = 0.0_Real;
+   Real AccMLiq  = 0.0_Real;
+   Real AccMSalt = 0.0_Real;
+   Real AccELiq  = 0.0_Real;
+   Real AccEIce  = 0.0_Real;
+
+   Real HTend = 0.0_Real;
+   Real TTend = 0.0_Real;
+   Real STend = 0.0_Real;
+
+   ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMLiq, AccMSalt,
+                          AccELiq, AccEIce, HTend, TTend, STend);
+
+   if (!isApprox(AccMIce, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilFormationTest warm: expected zero AccMIce, got {}",
+                  AccMIce);
+   }
+
+   if (!isApprox(AccMSalt, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilFormationTest warm: expected zero AccMSalt, got {}",
+                  AccMSalt);
+   }
+   if (!isApprox(AccMLiq, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilFormationTest warm: expected zero AccMLiq, got {}",
+                  AccMLiq);
+   }
+
+   if (!isApprox(AccELiq, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilFormationTest warm: expected zero AccELiq, got {}",
+                  AccELiq);
+   }
+   if (!isApprox(AccEIce, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilFormationTest warm: expected zero AccEIce, got {}",
+                  AccEIce);
+   }
+   if (!isApprox(HTend, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilFormationTest warm: expected zero HTend, got {}",
+                  HTend);
+   }
+
+   if (!isApprox(TTend, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilFormationTest warm: expected zero TTend, got {}",
+                  TTend);
+   }
+
+   if (!isApprox(STend, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilFormationTest warm: expected zero STend, got {}",
+                  STend);
+   }
+   LOG_INFO(
+       "FrazilFormationTestWarm: AccMIce = {}, AccMLiq = {}, AccMSalt = {}, "
+       "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}",
+       AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend);
+}
+
+// this test only excercises the frazil formation functor (no melt)
+// in a warm case: the frazil FORMATION terms should all be zero
+void testBasicFrazilFormationWarm() {
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
+
+   VCoord->NVertLayers = NVertLayers;
+
+   const Real SAIn = 35.0_Real;
+   const Real CTIn = 10.0_Real;
+   const Real PIn  = 100.0_Real;
+   const Real h    = 10.0_Real;
+   const Real RTol = 1e-10_Real;
+
+   (void)Mesh;
+
+   BasicFrazilFormation ComputeFrazilFormation;
+
+   Real AccMIce  = 0.0_Real;
+   Real AccMLiq  = 0.0_Real;
+   Real AccMSalt = 0.0_Real;
+   Real AccELiq  = 0.0_Real;
+   Real AccEIce  = 0.0_Real;
+
+   Real HTend = 0.0_Real;
+   Real TTend = 0.0_Real;
+   Real STend = 0.0_Real;
+
+   Real CTfrz = gsw_ct_freezing_poly(SAIn, PIn, 0.0_Real);
+
+   ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMSalt, AccEIce, HTend,
+                          TTend, STend, CTfrz);
+
+   if (!isApprox(AccMIce, 0.0_Real, RTol)) {
+      ABORT_ERROR(
+          "FrazilBasicFormationTest warm: expected zero AccMIce, got {}",
+          AccMIce);
+   }
+
+   if (!isApprox(AccMSalt, 0.0_Real, RTol)) {
+      ABORT_ERROR(
+          "FrazilBasicFormationTest warm: expected zero AccMSalt, got {}",
+          AccMSalt);
+   }
+
+   if (!isApprox(AccEIce, 0.0_Real, RTol)) {
+      ABORT_ERROR(
+          "FrazilBasicFormationTest warm: expected zero AccEIce, got {}",
+          AccEIce);
+   }
+   if (!isApprox(HTend, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilBasicFormationTest warm: expected zero HTend, got {}",
+                  HTend);
+   }
+
+   if (!isApprox(TTend, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilBasicFormationTest warm: expected zero TTend, got {}",
+                  TTend);
+   }
+
+   if (!isApprox(STend, 0.0_Real, RTol)) {
+      ABORT_ERROR("FrazilBasicFormationTest warm: expected zero STend, got {}",
+                  STend);
+   }
+   LOG_INFO("FrazilBasicFormationTestWarm: AccMIce = {}, AccMLiq = {}, "
+            "AccMSalt = {}, "
+            "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}",
+            AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend);
+}
+
+// this test only excercises the frazil formation functor (no melt)
+// in a cold case: the frazil terms should be
+// - strictly positive for ice, liquid, and salt mass
+//- strictly negative for ice and liquid energy
+// - positive for T tendency and negative for S, H tendencies
+void testBasicFrazilFormationCold() {
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
+
+   VCoord->NVertLayers = NVertLayers;
+
+   const Real SAIn = 35.0_Real;
+   const Real CTIn = -2.0_Real;
+   const Real PIn  = 100.0_Real;
+   const Real h    = 10.0_Real;
+   const Real RTol = 1e-10_Real;
+
+   (void)Mesh;
+
+   BasicFrazilFormation ComputeFrazilFormation;
+
+   Real AccMIce  = 0.0_Real;
+   Real AccMLiq  = 0.0_Real;
+   Real AccMSalt = 0.0_Real;
+   Real AccELiq  = 0.0_Real;
+   Real AccEIce  = 0.0_Real;
+
+   Real HTend = 0.0_Real;
+   Real TTend = 0.0_Real;
+   Real STend = 0.0_Real;
+
+   Real CTfrz = gsw_ct_freezing_poly(SAIn, PIn, 0.0_Real);
+
+   ComputeFrazilFormation(SAIn, CTIn, PIn, h, AccMIce, AccMSalt, AccEIce, HTend,
+                          TTend, STend, CTfrz);
+
+   if (AccMIce <= 0.0_Real) {
+      ABORT_ERROR("FrazilBasicFormationTestCold: accumulated ice mass is "
+                  "non-positive: {}",
+                  AccMIce);
+   }
+
+   if (AccMSalt <= 0.0_Real) {
+      ABORT_ERROR("FrazilBasicFormationTestCold: accumulated salt mass is "
+                  "non-positive: {}",
+                  AccMSalt);
+   }
+
+   if (AccEIce >= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilBasicFormationTestCold: accumulated ice energy is positive "
+          "(exp. negative): {}",
+          AccEIce);
+   }
+   if (HTend >= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilBasicFormationTestCold: HTend is positive (exp. negative): {}",
+          HTend);
+   }
+   if (TTend <= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilBasicFormationTestCold: TTend is negative (exp. positive): {}",
+          TTend);
+   }
+   if (STend >= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilBasicFormationTestCold: STend is positive (exp. negative): {}",
+          STend);
+   }
+   LOG_INFO("FrazilBasicFormationTestCold: AccMIce = {}, AccMLiq = {}, "
+            "AccMSalt = {}, "
+            "AccELiq = {}, AccEIce = {}, HTend = {}, TTend = {}, STend = {}",
+            AccMIce, AccMLiq, AccMSalt, AccELiq, AccEIce, HTend, TTend, STend);
+}
+
+// this test exercises the frazil formation and melt functors
+// in a column of water with both cold and warm layers.
+// It turns to frazil column conservation check.
+// In dev, there is extra verbose logging in the frazil code (TBRemoved).
+void testComputeFrazilColumn() {
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
+   auto *TestFrazil  = Frazil::getDefault();
+
+   if (!TestFrazil) {
+      ABORT_ERROR("FrazilTestColumn: default frazil object is null");
+   }
+
+   const Real RTol    = 1e-10_Real;
+   const Real SACold  = 35.0_Real;
+   const Real PRef    = 100.0_Real;
+   const Real HRef    = 10.0_Real;
+   const Real CTCold  = -2.0_Real;
+   const Real CTWarm  = 0.0_Real;
+   const Real CTWarm2 = -1.9_Real;
+
+   Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers);
+   Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers);
+   Array2DReal P("P", Mesh->NCellsSize, NVertLayers);
+   Array2DReal H("H", Mesh->NCellsSize, NVertLayers);
+
+   deepCopy(SA, SACold);
+   deepCopy(CT, CTWarm);
+   deepCopy(P, PRef);
+   deepCopy(H, HRef);
+
+   deepCopy(TestFrazil->AccMIce, 0.0_Real);
+   deepCopy(TestFrazil->AccMLiq, 0.0_Real);
+   deepCopy(TestFrazil->AccMSalt, 0.0_Real);
+   deepCopy(TestFrazil->AccELiq, 0.0_Real);
+   deepCopy(TestFrazil->AccEIce, 0.0_Real);
+   deepCopy(TestFrazil->FrazilHTend, 0.0_Real);
+   deepCopy(TestFrazil->FrazilTTend, 0.0_Real);
+   deepCopy(TestFrazil->FrazilSTend, 0.0_Real);
+
+   auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell);
+   auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell);
+
+   const I4 ICell = 0;
+   const I4 KMin  = MinLayerCellH(ICell);
+   const I4 KMax  = MaxLayerCellH(ICell);
+   if ((KMax - KMin + 1) < 4) {
+      ABORT_ERROR("FrazilTestColumn: cell {} has fewer than 4 active layers",
+                  ICell);
+   }
+
+   const I4 KBottom0  = KMax;
+   const I4 KBottom1  = KMax - 1;
+   const I4 KWarm     = KMax - 2;
+   const I4 KTopCold  = KMax - 3;
+   const I4 KCold2    = KMin + 3;
+   const I4 KCold3    = KMin + 2;
+   const I4 KWarm2    = KMin + 1;
+   const I4 KTopCold2 = KMin;
+
+   auto CTH               = createHostMirrorCopy(CT);
+   CTH(ICell, KBottom0)   = CTCold;
+   CTH(ICell, KBottom1)   = CTCold;
+   CTH(ICell, KWarm)      = CTWarm;
+   CTH(ICell, KTopCold)   = CTCold;
+   CTH(ICell, KCold2 + 2) = CTCold;
+   CTH(ICell, KCold2 + 1) = CTCold - .5_Real;
+   CTH(ICell, KCold2)     = CTCold;
+   CTH(ICell, KCold3)     = CTCold;
+   CTH(ICell, KWarm2)     = CTWarm2;
+   CTH(ICell, KTopCold2)  = CTCold;
+   deepCopy(CT, CTH);
+
+   const bool SavedConservationCheck = TestFrazil->conservationCheck;
+   TestFrazil->conservationCheck     = true;
+   // TestFrazil->conservationCheck = false; // manual toggle for BasicFrazil
+   TestFrazil->computeFrazil(CT, SA, P, H);
+   TestFrazil->conservationCheck = SavedConservationCheck;
+
+   auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend);
+   auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend);
+   auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend);
+
+   if (HTendH(ICell, KBottom0) >= 0.0_Real ||
+       TTendH(ICell, KBottom0) <= 0.0_Real ||
+       STendH(ICell, KBottom0) >= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilTestColumn: bottom cold layer sign check failed (HTend<0, "
+          "TTend>0, STend<0 expected)");
+   }
+
+   if (HTendH(ICell, KBottom1) >= 0.0_Real ||
+       TTendH(ICell, KBottom1) <= 0.0_Real ||
+       STendH(ICell, KBottom1) >= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilTestColumn: second cold layer sign check failed (HTend<0, "
+          "TTend>0, STend<0 expected)");
+   }
+
+   if (HTendH(ICell, KWarm) < 0.0_Real || TTendH(ICell, KWarm) > 0.0_Real ||
+       STendH(ICell, KWarm) < 0.0_Real) {
+      ABORT_ERROR("FrazilTestColumn: warm layer sign check failed (HTend>=0, "
+                  "TTend<=0, STend>=0 expected)");
+   }
+
+   if (HTendH(ICell, KTopCold) >= 0.0_Real ||
+       TTendH(ICell, KTopCold) <= 0.0_Real ||
+       STendH(ICell, KTopCold) >= 0.0_Real) {
+      ABORT_ERROR(
+          "FrazilTestColumn: top cold layer sign check failed (HTend<0, "
+          "TTend>0, STend<0 expected)");
+   }
+
+   LOG_INFO(
+       "FrazilTestColumn: XTend branch-switching checks passed for ICell={}",
+       ICell);
+}
+
+// this test exercises the frazil formation and melt functors
+// with a depth limit set. Layers deeper than the depth limit
+// should have zero frazil tendencies.
+void testComputeFrazilDepthLimit() {
+   const auto Mesh   = HorzMesh::getDefault();
+   const auto VCoord = VertCoord::getDefault();
+   auto *TestFrazil  = Frazil::getDefault();
+
+   if (!TestFrazil) {
+      ABORT_ERROR("FrazilTestColumn: default frazil object is null");
+   }
+
+   const Real RTol    = 1e-12_Real;
+   const Real SACold  = 35.0_Real;
+   const Real PRef    = 100.0_Real;
+   const Real HRef    = 10.0_Real;
+   const Real CTCold  = -2.0_Real;
+   const Real CTWarm  = 0.0_Real;
+   const Real CTWarm2 = -1.9_Real;
+
+   Array2DReal SA("SA", Mesh->NCellsSize, NVertLayers);
+   Array2DReal CT("CT", Mesh->NCellsSize, NVertLayers);
+   Array2DReal P("P", Mesh->NCellsSize, NVertLayers);
+   Array2DReal H("H", Mesh->NCellsSize, NVertLayers);
+
+   deepCopy(SA, SACold);
+   deepCopy(CT, CTWarm);
+   deepCopy(P, PRef);
+   deepCopy(H, HRef);
+
+   deepCopy(TestFrazil->AccMIce, 0.0_Real);
+   deepCopy(TestFrazil->AccMLiq, 0.0_Real);
+   deepCopy(TestFrazil->AccMSalt, 0.0_Real);
+   deepCopy(TestFrazil->AccELiq, 0.0_Real);
+   deepCopy(TestFrazil->AccEIce, 0.0_Real);
+   deepCopy(TestFrazil->FrazilHTend, 0.0_Real);
+   deepCopy(TestFrazil->FrazilTTend, 0.0_Real);
+   deepCopy(TestFrazil->FrazilSTend, 0.0_Real);
+
+   auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell);
+   auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell);
+
+   const I4 ICell = 0;
+   const I4 KMin  = MinLayerCellH(ICell);
+   const I4 KMax  = MaxLayerCellH(ICell);
+   if ((KMax - KMin + 1) < 10) {
+      ABORT_ERROR("FrazilTestColumn: cell {} has fewer than 10 active layers",
+                  ICell);
+   }
+
+   const I4 KBottom0  = KMax;
+   const I4 KBottom1  = KMax - 1;
+   const I4 KWarm     = KMax - 2;
+   const I4 KTopCold  = KMax - 3;
+   const I4 KCold2    = KMin + 3;
+   const I4 KCold3    = KMin + 2;
+   const I4 KWarm2    = KMin + 1;
+   const I4 KTopCold2 = KMin;
+
+   auto CTH               = createHostMirrorCopy(CT);
+   CTH(ICell, KBottom0)   = CTCold;
+   CTH(ICell, KBottom1)   = CTCold;
+   CTH(ICell, KWarm)      = CTWarm;
+   CTH(ICell, KTopCold)   = CTCold;
+   CTH(ICell, KCold2 + 2) = CTCold;
+   CTH(ICell, KCold2 + 1) = CTCold - .5_Real;
+   CTH(ICell, KCold2)     = CTCold;
+   CTH(ICell, KCold3)     = CTCold;
+   CTH(ICell, KWarm2)     = CTWarm2;
+   CTH(ICell, KTopCold2)  = CTCold;
+   deepCopy(CT, CTH);
+
+   const bool SavedConservationCheck = TestFrazil->conservationCheck;
+   const bool SavedDepthLimit        = TestFrazil->depthLimit;
+   const Real TestDepthLimit         = 35.0_Real; // this needs to be positive
+   // if TestDepthLimit is negative, test will fail:
+   // - the code assume depthlimit < 0 mean no limit (i.e. full depth frazil)
+   // - the test below will exclude all layers and fail because Tend !=0.
+
+   // Populate GeomZMid explicitly for the test column.
+   auto GeomZMidH = createHostMirrorCopy(VCoord->GeomZMid);
+   for (I4 K = KMin; K <= KMax; ++K) {
+      GeomZMidH(ICell, K) = -10.0_Real * (K - KMin + 1);
+   }
+   deepCopy(VCoord->GeomZMid, GeomZMidH);
+
+   TestFrazil->conservationCheck = true;
+   TestFrazil->depthLimit        = TestDepthLimit;
+   TestFrazil->computeFrazil(CT, SA, P, H);
+   TestFrazil->conservationCheck = SavedConservationCheck;
+   TestFrazil->depthLimit        = SavedDepthLimit;
+
+   auto HTendH = createHostMirrorCopy(TestFrazil->FrazilHTend);
+   auto TTendH = createHostMirrorCopy(TestFrazil->FrazilTTend);
+   auto STendH = createHostMirrorCopy(TestFrazil->FrazilSTend);
+
+   bool FoundExcludedLayer = false;
+   for (I4 K = KMin; K <= KMax; ++K) {
+      const Real Depth    = GeomZMidH(ICell, K);
+      const Real AbsDepth = Depth < 0.0_Real ? -Depth : Depth;
+
+      if (AbsDepth > TestDepthLimit) {
+         FoundExcludedLayer = true;
+         if (!isApprox(HTendH(ICell, K), 0.0_Real, RTol) ||
+             !isApprox(TTendH(ICell, K), 0.0_Real, RTol) ||
+             !isApprox(STendH(ICell, K), 0.0_Real, RTol)) {
+            ABORT_ERROR("FrazilDepthLimitTest: excluded layer K={} has "
+                        "non-zero tendencies (HTend={}, TTend={}, STend={})",
+                        K, HTendH(ICell, K), TTendH(ICell, K),
+                        STendH(ICell, K));
+         }
+      }
+   }
+   if (!FoundExcludedLayer) {
+      ABORT_ERROR("FrazilDepthLimitTest: no layers were excluded for ICell={} "
+                  "with depthLimit={}",
+                  ICell, TestDepthLimit);
+   }
+
+   LOG_INFO("FrazilDepthLimitTest: depthLimit={} exclusion check passed for "
+            "ICell={}",
+            TestDepthLimit, ICell);
+}
+
+void frazilTest(const std::string &MeshFile = "OmegaMesh.nc") {
+   initFrazilTest(MeshFile);
+   testFrazilFormationCold();
+   testFrazilFormationWarm();
+   testBasicFrazilFormationCold();
+   testBasicFrazilFormationWarm();
+   testComputeFrazilColumn();
+   testComputeFrazilDepthLimit();
+   finalizeFrazilTest();
+}
+
+int main(int argc, char *argv[]) {
+   MPI_Init(&argc, &argv);
+   Kokkos::initialize(argc, argv);
+   Pacer::initialize(MPI_COMM_WORLD);
+   Pacer::setPrefix("Omega:");
+
+   frazilTest();
+
+   LOG_INFO("------ Frazil Unit Tests Successful ------");
+
+   Pacer::finalize();
+   Kokkos::finalize();
+   MPI_Finalize();
+
+   return 0;
+}
diff --git a/components/omega/test/ocn/TendenciesTest.cpp b/components/omega/test/ocn/TendenciesTest.cpp
index 5268d0a436ea..0099c17aec8d 100644
--- a/components/omega/test/ocn/TendenciesTest.cpp
+++ b/components/omega/test/ocn/TendenciesTest.cpp
@@ -139,6 +139,7 @@ int initTendenciesTest(const std::string &mesh) {
    VertCoord::init();
 
    Tracers::init();
+   Frazil::init();
    VertAdv::init();
    PressureGrad::init();
    Eos::init();
@@ -244,6 +245,7 @@ int testTendencies() {
        "BaselineNormalVelocityTend", Mesh->NEdgesSize, VCoord->NVertLayers);
 
    DefTendencies->SfcStressForcing.Enabled = false;
+   DefTendencies->FrazilTerm.Enabled = false; // frazil needs to be off for now
    DefTendencies->computeAllTendencies(State, AuxState, TracerArray,
                                        ThickTimeLevel, VelTimeLevel,
                                        TracerTimeLevel, Time, Interval);
@@ -305,6 +307,145 @@ int testTendencies() {
 
    DefTendencies->SfcStressForcing.Enabled = OrigSfcStressEnabled;
 
+   // Verify frazil tendencies are plumbed into pseudo-thickness and tracer
+   // tendencies by comparing runs with FrazilTerm disabled/enabled.
+   I4 TempTracerIndex = -1;
+   I4 SaltTracerIndex = -1;
+   Tracers::getIndex(TempTracerIndex, "Temperature");
+   Tracers::getIndex(SaltTracerIndex, "Salinity");
+
+   if (TempTracerIndex < 0 || SaltTracerIndex < 0) {
+      Err++;
+      LOG_ERROR("TendenciesTest: missing Temperature/Salinity tracer indices");
+   } else {
+      auto MinLayerCellH = createHostMirrorCopy(VCoord->MinLayerCell);
+      auto MaxLayerCellH = createHostMirrorCopy(VCoord->MaxLayerCell);
+      auto TracerArrayH  = createHostMirrorCopy(TracerArray);
+
+      // Set Salinity in active layers to [31, 33], top-to-bottom ramp.
+      for (I4 ICell = 0; ICell < Mesh->NCellsAll; ++ICell) {
+         const I4 KMin = MinLayerCellH(ICell);
+         const I4 KMax = MaxLayerCellH(ICell);
+
+         for (I4 K = KMin; K <= KMax; ++K) {
+            Real frac = 0.0_Real;
+            if (KMax > KMin) {
+               frac =
+                   static_cast<Real>(K - KMin) / static_cast<Real>(KMax - KMin);
+            }
+            TracerArrayH(SaltTracerIndex, ICell, K) =
+                31.0_Real + 2.0_Real * frac;
+         }
+      }
+
+      // Optional: keep Temperature warm everywhere first, then add one cold
+      // point.
+      for (I4 ICell = 0; ICell < Mesh->NCellsAll; ++ICell) {
+         const I4 KMin = MinLayerCellH(ICell);
+         const I4 KMax = MaxLayerCellH(ICell);
+         for (I4 K = KMin; K <= KMax; ++K) {
+            TracerArrayH(TempTracerIndex, ICell, K) = 8.0_Real;
+         }
+      }
+      const I4 FrazilTestCell  = 0;
+      const I4 FrazilTestLayer = MinLayerCellH(FrazilTestCell);
+      TracerArrayH(TempTracerIndex, FrazilTestCell, FrazilTestLayer) =
+          -2.0_Real;
+      deepCopy(TracerArray, TracerArrayH);
+
+      LOG_INFO("TendenciesTest: minLayerCell = {}", FrazilTestLayer);
+
+      Array2DReal BaselinePseudoThicknessTend(
+          "BaselinePseudoThicknessTend", Mesh->NCellsSize, VCoord->NVertLayers);
+      Array3DReal BaselineTracerTend("BaselineTracerTend",
+                                     Tracers::getNumTracers(), Mesh->NCellsSize,
+                                     VCoord->NVertLayers);
+
+      const bool OrigFrazilEnabled = DefTendencies->FrazilTerm.Enabled;
+
+      DefTendencies->FrazilTerm.Enabled = false;
+      DefTendencies->computeAllTendencies(State, AuxState, TracerArray,
+                                          ThickTimeLevel, VelTimeLevel,
+                                          TracerTimeLevel, Time, Interval);
+      deepCopy(BaselinePseudoThicknessTend, DefTendencies->PseudoThicknessTend);
+      deepCopy(BaselineTracerTend, DefTendencies->TracerTend);
+
+      deepCopy(DefTendencies->PseudoThicknessTend, 0.0_Real);
+      deepCopy(DefTendencies->NormalVelocityTend, 0.0_Real);
+      deepCopy(DefTendencies->TracerTend, 0.0_Real);
+
+      DefTendencies->FrazilTerm.Enabled = true;
+      DefTendencies->computeAllTendencies(State, AuxState, TracerArray,
+                                          ThickTimeLevel, VelTimeLevel,
+                                          TracerTimeLevel, Time, Interval);
+
+      Array2DReal PseudoThicknessTendDiff(
+          "PseudoThicknessTendDiff", Mesh->NCellsSize, VCoord->NVertLayers);
+      Array2DReal TempTracerTendDiff("TempTracerTendDiff", Mesh->NCellsSize,
+                                     VCoord->NVertLayers);
+      Array2DReal SaltTracerTendDiff("SaltTracerTendDiff", Mesh->NCellsSize,
+                                     VCoord->NVertLayers);
+
+      OMEGA_SCOPE(LocPseudoThicknessTendDiff, PseudoThicknessTendDiff);
+      OMEGA_SCOPE(LocTempTracerTendDiff, TempTracerTendDiff);
+      OMEGA_SCOPE(LocSaltTracerTendDiff, SaltTracerTendDiff);
+      OMEGA_SCOPE(LocPseudoThicknessTend, DefTendencies->PseudoThicknessTend);
+      OMEGA_SCOPE(LocBaselinePseudoThicknessTend, BaselinePseudoThicknessTend);
+      OMEGA_SCOPE(LocTracerTend, DefTendencies->TracerTend);
+      OMEGA_SCOPE(LocBaselineTracerTend, BaselineTracerTend);
+      OMEGA_SCOPE(LocMinLayerCell, VCoord->MinLayerCell);
+      OMEGA_SCOPE(LocMaxLayerCell, VCoord->MaxLayerCell);
+
+      parallelForOuter(
+          "TendenciesTest:FrazilTendDiff", {Mesh->NCellsAll},
+          KOKKOS_LAMBDA(int ICell, const TeamMember &Team) {
+             const int KMin = LocMinLayerCell(ICell);
+             const int KMax = LocMaxLayerCell(ICell);
+
+             parallelForInner(
+                 Team, Range{KMin, KMax}, INNER_LAMBDA(int K) {
+                    LocPseudoThicknessTendDiff(ICell, K) =
+                        Kokkos::abs(LocPseudoThicknessTend(ICell, K) -
+                                    LocBaselinePseudoThicknessTend(ICell, K));
+                    LocTempTracerTendDiff(ICell, K) = Kokkos::abs(
+                        LocTracerTend(TempTracerIndex, ICell, K) -
+                        LocBaselineTracerTend(TempTracerIndex, ICell, K));
+                    LocSaltTracerTendDiff(ICell, K) = Kokkos::abs(
+                        LocTracerTend(SaltTracerIndex, ICell, K) -
+                        LocBaselineTracerTend(SaltTracerIndex, ICell, K));
+                 });
+          });
+
+      const Real PseudoDelta = sum(PseudoThicknessTendDiff, Mesh->NCellsOwned,
+                                   VCoord->MinLayerCell, VCoord->MaxLayerCell);
+      const Real TempDelta   = sum(TempTracerTendDiff, Mesh->NCellsOwned,
+                                   VCoord->MinLayerCell, VCoord->MaxLayerCell);
+      const Real SaltDelta   = sum(SaltTracerTendDiff, Mesh->NCellsOwned,
+                                   VCoord->MinLayerCell, VCoord->MaxLayerCell);
+
+      constexpr Real FrazilDeltaATol = 1e-12_Real;
+      if (!Kokkos::isfinite(PseudoDelta) ||
+          isApprox(PseudoDelta, 0._Real, 0._Real, FrazilDeltaATol)) {
+         Err++;
+         LOG_ERROR("TendenciesTest: Frazil did not change "
+                   "PseudoThicknessTend");
+      }
+      if (!Kokkos::isfinite(TempDelta) ||
+          isApprox(TempDelta, 0._Real, 0._Real, FrazilDeltaATol)) {
+         Err++;
+         LOG_ERROR("TendenciesTest: Frazil did not change Temperature tracer "
+                   "tendency");
+      }
+      if (!Kokkos::isfinite(SaltDelta) ||
+          isApprox(SaltDelta, 0._Real, 0._Real, FrazilDeltaATol)) {
+         Err++;
+         LOG_ERROR("TendenciesTest: Frazil did not change Salinity tracer "
+                   "tendency");
+      }
+      // add a log info if no errors?
+      DefTendencies->FrazilTerm.Enabled = OrigFrazilEnabled;
+   }
+
    // check that everything got computed correctly
    int NCellsOwned = Mesh->NCellsOwned;
    int NEdgesOwned = Mesh->NEdgesOwned;
@@ -326,12 +467,12 @@ int testTendencies() {
       LOG_ERROR("TendenciesTest: NormVelTendSum FAIL");
    }
 
-   const Real TraceTendSum =
+   const Real TracerTendSum =
        sum(DefTendencies->TracerTend, NTracers, NCellsOwned,
            VCoord->MinLayerCell, VCoord->MaxLayerCell);
-   if (!Kokkos::isfinite(TraceTendSum) || TraceTendSum == 0) {
+   if (!Kokkos::isfinite(TracerTendSum) || TracerTendSum == 0) {
       Err++;
-      LOG_ERROR("TendenciesTest: TraceTendSum FAIL");
+      LOG_ERROR("TendenciesTest: TracerTendSum FAIL");
    }
 
    Tendencies::clear();
@@ -342,6 +483,7 @@ int testTendencies() {
 void finalizeTendenciesTest() {
    Forcing::clear();
    Tracers::clear();
+   Frazil::clear();
    PressureGrad::clear();
    VertMix::destroyInstance();
    Eos::destroyInstance();