diff --git a/components/omega/configs/Default.yml b/components/omega/configs/Default.yml
index 9db42e178466..86b032a5f77c 100644
--- a/components/omega/configs/Default.yml
+++ b/components/omega/configs/Default.yml
@@ -63,6 +63,7 @@ Omega:
     VelocityVertAdvTendencyEnable: true
     TracerVertAdvTendencyEnable: true
     PressureGradTendencyEnable: true
+    FrazilTendencyEnable: true
   ManufacturedSolution:
     WavelengthX: 5.0e6
     WavelengthY: 4.33013e6
@@ -76,6 +77,12 @@ Omega:
       DRhoDT: -0.2
       DRhoDS: 0.8
       RhoT0S0: 1000.0
+  Frazil:
+    FrazilType: teos
+    MassLimit: 0.1
+    Phi: 0.75
+    DepthLimit: -1.0
+    ConservationCheck: false
   VertMix:
     Background:
       Diffusivity: 1.0e-5
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..9a323f51b76f
--- /dev/null
+++ b/components/omega/src/ocn/Frazil.cpp
@@ -0,0 +1,385 @@
+//===-- 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() {}
+
+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;
+      ABORT_ERROR("Frazil::create: BasicFrazil not supported yet");
+   } 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");
+
+   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::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");
+   }
+
+   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 = (depthLimit >= 0.0_Real);
+
+          if (Limit) {
+             HasKlim = false;
+             for (I4 K = KMax; K >= KMin; --K) {
+                if (Kokkos::abs(LocGeomZMid(ICell, K)) <= depthLimit) {
+                   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);
+                }
+                //    else {
+                //       temporary kernel logging - TBRemoved
+                //        LOG_INFO("warm layer but no ice to melt");
+                //    }
+             }
+
+             //  // temporary kernel logging -- TBRemoved
+             //  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;
+          }
+
+          // 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 parallelFor
+
+   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..bf99f3caedd5
--- /dev/null
+++ b/components/omega/src/ocn/Frazil.h
@@ -0,0 +1,257 @@
+#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 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);
+   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;
+   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 ceabc0ede033..64f975570b10 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"
@@ -199,6 +200,7 @@ static int initOmegaModulesImpl(MPI_Comm Comm) {
    Forcing::init();
    AuxiliaryState::init();
    Eos::init();
+   Frazil::init();
    PressureGrad::init();
    Tendencies::init();
 
diff --git a/components/omega/src/ocn/Tendencies.cpp b/components/omega/src/ocn/Tendencies.cpp
index ba0df5dd692a..2dd13399a464 100644
--- a/components/omega/src/ocn/Tendencies.cpp
+++ b/components/omega/src/ocn/Tendencies.cpp
@@ -263,6 +263,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(
@@ -392,7 +396,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) {
 
    // Tendency arrays
@@ -823,6 +827,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 5d26ca041458..866f1794b982 100644
--- a/components/omega/src/ocn/Tendencies.h
+++ b/components/omega/src/ocn/Tendencies.h
@@ -72,6 +72,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 4d7ce499941f..f62f14430dc3 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"
@@ -115,6 +116,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 50c42bf0a2f0..a7c949da1f19 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 1678a08470f6..51620505568c 100644
--- a/components/omega/test/infra/IOStreamTest.cpp
+++ b/components/omega/test/infra/IOStreamTest.cpp
@@ -108,12 +108,13 @@ void initIOStreamTest(Clock *&ModelClock // Model clock
 
    PressureGrad::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();
-
    // Initialize Aux State
    AuxiliaryState::init();
 
diff --git a/components/omega/test/ocn/FrazilTest.cpp b/components/omega/test/ocn/FrazilTest.cpp
new file mode 100644
index 000000000000..b118b0f12d02
--- /dev/null
+++ b/components/omega/test/ocn/FrazilTest.cpp
@@ -0,0 +1,494 @@
+//===-- 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 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->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();
+   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 5c8bb8d3d3f0..e0b151cf4b88 100644
--- a/components/omega/test/ocn/TendenciesTest.cpp
+++ b/components/omega/test/ocn/TendenciesTest.cpp
@@ -138,6 +138,7 @@ int initTendenciesTest(const std::string &mesh) {
    VertCoord::init();
 
    Tracers::init();
+   Frazil::init();
    VertAdv::init();
    PressureGrad::init();
    Eos::init();
@@ -241,6 +242,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);
@@ -302,6 +304,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;
@@ -323,12 +464,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();
@@ -339,6 +480,7 @@ int testTendencies() {
 void finalizeTendenciesTest() {
    Forcing::clear();
    Tracers::clear();
+   Frazil::clear();
    PressureGrad::clear();
    Eos::destroyInstance();
    AuxiliaryState::clear();