From 1862ec88e5165bba6be69d6731a7f7650b04436c Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 9 Mar 2026 14:35:26 -0700 Subject: [PATCH 1/5] added Ct/Pt conversions; added global constants - not hooked yet - compile on cpu --- components/omega/src/ocn/Eos.h | 119 +++++++++++++++++++++ components/omega/src/ocn/GlobalConstants.h | 15 ++- 2 files changed, 130 insertions(+), 4 deletions(-) diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 06c9695f7267..738dd9302f21 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -238,6 +238,125 @@ class Teos10Eos { return V0; } + /// Calculate 2nd derivative of Gibbs wrt pot temp at ref P for TEOS-10 + KOKKOS_FUNCTION Real calcGibbsDerivPt0Pt0(Real Sa, Real P) const { + constexpr Real gsw_sfac = 0.0248826675584615; // replace + Real x2 = gsw_sfac * Sa; + Real x = Kokkos::sqrt(x2); + Real y = P * 0.025; + + Real g03 = + -24715.571866078 + + y * (4420.4472249096725 + + y * (-1778.231237203896 + + y * (1160.5182516851419 + + y * (-569.531539542516 + y * 128.13429152494615)))); + + Real g08 = + x2 * + (1760.062705994408 + + x * (-86.1329351956084 + + x * (-137.1145018408982 + + y * (296.20061691375236 + + y * (-205.67709290374563 + 49.9394019139016 * y))) + + y * (-60.136422517125 + y * 10.50720794170734)) + + y * (-1351.605895580406 + + y * (1097.1125373015109 + + y * (-433.20648175062206 + 63.905091254154904 * y)))); + + return ((g03 + g08) * 0.000625); + } + + /// Calculate Pot Temmperature from Conservative Temp + KOKKOS_FUNCTION Real calcPtFromCt(Real Sa, Real Ct) const { + constexpr Real gsw_t0 = 273.15; // replace with TkFrz + constexpr Real gsw_cp0 = 3991.86795711963; // replace + constexpr Real gsw_ups = 35.16504 / 35.0; // replace + constexpr Real a0 = -1.446013646344788e-2; + constexpr Real a1 = -3.305308995852924e-3; + constexpr Real a2 = 1.062415929128982e-4; + constexpr Real a3 = 9.477566673794488e-1; + constexpr Real a4 = 2.166591947736613e-3; + constexpr Real a5 = 3.828842955039902e-3; + constexpr Real b0 = 1.000000000000000e0; + constexpr Real b1 = 6.506097115635800e-4; + constexpr Real b2 = 3.830289486850898e-3; + constexpr Real b3 = 1.247811760368034e-6; + Real a5ct, b3ct, ct_factor, pt_num, pt_recden, ct_diff; + Real pt, pt_old, ptm, dpt_dct, s1; + + s1 = Sa / gsw_ups; + + a5ct = a5 * Ct; + b3ct = b3 * Ct; + + ct_factor = (a3 + a4 * s1 + a5ct); + pt_num = a0 + s1 * (a1 + a2 * s1) + Ct * ct_factor; + pt_recden = 1.0 / (b0 + b1 * s1 + Ct * (b2 + b3ct)); + pt = pt_num * pt_recden; + + dpt_dct = pt_recden * (ct_factor + a5ct - (b2 + b3ct + b3ct) * pt); + + ct_diff = calcCtFromPt(Sa, pt) - Ct; + pt_old = pt; + pt = pt_old - ct_diff * dpt_dct; + ptm = 0.5 * (pt + pt_old); + + dpt_dct = -gsw_cp0 / ((ptm + gsw_t0) * calcGibbsDerivPt0Pt0(Sa, ptm)); + + pt = pt_old - ct_diff * dpt_dct; + ct_diff = calcCtFromPt(Sa, pt) - Ct; + pt_old = pt; + return (pt_old - ct_diff * dpt_dct); + } + + /// Calculate Conservative Temmperature from Potential Temp + KOKKOS_FUNCTION Real calcCtFromPt(Real Sa, Real Pt) const { + constexpr Real gsw_cp0 = 3991.86795711963; // replace + constexpr Real gsw_sfac = 0.0248826675584615; // replace + Real x2, x, y, pot_enthalpy; + + x2 = gsw_sfac * Sa; + x = Kokkos::sqrt(x2); + y = Pt * 0.025e0; /*! normalize for F03 and F08 */ + pot_enthalpy = + 61.01362420681071e0 + + y * (168776.46138048015e0 + + y * (-2735.2785605119625e0 + + y * (2574.2164453821433e0 + + y * (-1536.6644434977543e0 + + y * (545.7340497931629e0 + + (-50.91091728474331e0 - + 18.30489878927802e0 * y) * + y))))) + + x2 * + (268.5520265845071e0 + + y * (-12019.028203559312e0 + + y * (3734.858026725145e0 + + y * (-2046.7671145057618e0 + + y * (465.28655623826234e0 + + (-0.6370820302376359e0 - + 10.650848542359153e0 * y) * + y)))) + + x * (937.2099110620707e0 + + y * (588.1802812170108e0 + y * (248.39476522971285e0 + + (-3.871557904936333e0 - + 2.6268019854268356e0 * y) * + y)) + + x * (-1687.914374187449e0 + + x * (246.9598888781377e0 + + x * (123.59576582457964e0 - + 48.5891069025409e0 * x)) + + y * (936.3206544460336e0 + + y * (-942.7827304544439e0 + + y * (369.4389437509002e0 + + (-33.83664947895248e0 - + 9.987880382780322e0 * y) * + y)))))); + + return (pot_enthalpy / gsw_cp0); + } + private: Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; diff --git a/components/omega/src/ocn/GlobalConstants.h b/components/omega/src/ocn/GlobalConstants.h index b13a4e804acb..bef552130b39 100644 --- a/components/omega/src/ocn/GlobalConstants.h +++ b/components/omega/src/ocn/GlobalConstants.h @@ -69,6 +69,8 @@ constexpr Real CpFw = pcd::pure_water_specific_heat_capacity_reference; constexpr Real CpSw = pcd::seawater_specific_heat_capacity_reference; // Specific heat capacity of seawater ~ J/(kg*K) (from Physical Constants // Dictionary) +constexpr Real Cp0Sw = + 3991.86795711963; // Specific heat capacity of seawater for use with constexpr Real CpIce = pcd::sea_ice_specific_heat_capacity_reference; // Specific heat capacity of ice ~ J/(kg*K) (from Physical Constants Dictionary) constexpr Real LatIce = pcd::latent_heat_of_fusion_reference; @@ -98,7 +100,11 @@ constexpr Real Rad2Deg = pcd::radian; // Radians to degrees (from Physical Constants Dictionary) constexpr Real Deg2Rad = pcd::degree; // Degrees to radians (from Physical Constants Dictionary) -constexpr Real Salt2PPt = 1000.0; // Salinity (kg/kg) to parts per thousand +constexpr Real Salt2PPt = 1000.0; // Salinity (kg/kg) to parts per thousand +constexpr Real SS0 = 35.16504; // Standard Ocean Reference Salinity (g/kg) +constexpr Real Psu2Gpkg = + SS0 / 35.0; // unit conversion factor for salinities (psu --> g/kg) +constexpr Real Sfac = 0.0248826675584615; // sfac = 1/(40*Psu2Gpkg) constexpr Real PPt2Salt = 1.0e-3; // Parts per thousand to salinity (kg/kg) constexpr Real Mass2Sv = 1.0e-12; // Mass flux (kg/s) to Sverdrup constexpr Real Heat2Pw = 4.186e-15; // Heat flux (W) to Petawatt @@ -109,11 +115,12 @@ constexpr Real Pa2Db = 1.0e-4; // Pascal to Decibar constexpr Real Cm2M = 1.0e-2; // Centimeters to meters constexpr Real M2Cm = 1.0e2; // Meters to centimeters constexpr Real HFluxFac = - 1.0 / (RhoSw * CpSw); // Heat flux (W/m^2) to temp flux (C*m/s) + 1.0 / (RhoSw * Cp0Sw); // Heat flux (W/m^2) to Conserv Temp flux (C*m/s) constexpr Real FwFluxFac = 1.e-6; // Fw flux (kg/m^2/s) to salt((msu/psu)*m/s) constexpr Real SaltFac = - -OcnRefSal * FwFluxFac; // Fw flux (kg/m^2/s) to salt flux (msu*m/s) -constexpr Real SFluxFac = 1.0; // Salt flux (kg/m^2/s) to salt flux (msu*m/s) + -OcnRefSal * FwFluxFac; // Fw flux (kg/m^2/s) to salt flux (msu*m/s) +constexpr Real SFluxFac = + 1.e3 / RhoSw; // Salt flux (kg/m^2/s) to salinity flux (m*(g/kg)/s) } // namespace OMEGA #endif From 0835de0df7a03ca91f6b4187074ef2b76088beaa Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 9 Mar 2026 14:57:31 -0700 Subject: [PATCH 2/5] switched to use global constants --- components/omega/src/ocn/Eos.h | 40 +++++++++++++++------------------- 1 file changed, 17 insertions(+), 23 deletions(-) diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index 738dd9302f21..d7f8f8ae2db3 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -240,10 +240,9 @@ class Teos10Eos { /// Calculate 2nd derivative of Gibbs wrt pot temp at ref P for TEOS-10 KOKKOS_FUNCTION Real calcGibbsDerivPt0Pt0(Real Sa, Real P) const { - constexpr Real gsw_sfac = 0.0248826675584615; // replace - Real x2 = gsw_sfac * Sa; - Real x = Kokkos::sqrt(x2); - Real y = P * 0.025; + Real x2 = Sfac * Sa; + Real x = Kokkos::sqrt(x2); + Real y = P * 0.025; Real g03 = -24715.571866078 + @@ -269,23 +268,20 @@ class Teos10Eos { /// Calculate Pot Temmperature from Conservative Temp KOKKOS_FUNCTION Real calcPtFromCt(Real Sa, Real Ct) const { - constexpr Real gsw_t0 = 273.15; // replace with TkFrz - constexpr Real gsw_cp0 = 3991.86795711963; // replace - constexpr Real gsw_ups = 35.16504 / 35.0; // replace - constexpr Real a0 = -1.446013646344788e-2; - constexpr Real a1 = -3.305308995852924e-3; - constexpr Real a2 = 1.062415929128982e-4; - constexpr Real a3 = 9.477566673794488e-1; - constexpr Real a4 = 2.166591947736613e-3; - constexpr Real a5 = 3.828842955039902e-3; - constexpr Real b0 = 1.000000000000000e0; - constexpr Real b1 = 6.506097115635800e-4; - constexpr Real b2 = 3.830289486850898e-3; - constexpr Real b3 = 1.247811760368034e-6; + constexpr Real a0 = -1.446013646344788e-2; + constexpr Real a1 = -3.305308995852924e-3; + constexpr Real a2 = 1.062415929128982e-4; + constexpr Real a3 = 9.477566673794488e-1; + constexpr Real a4 = 2.166591947736613e-3; + constexpr Real a5 = 3.828842955039902e-3; + constexpr Real b0 = 1.000000000000000e0; + constexpr Real b1 = 6.506097115635800e-4; + constexpr Real b2 = 3.830289486850898e-3; + constexpr Real b3 = 1.247811760368034e-6; Real a5ct, b3ct, ct_factor, pt_num, pt_recden, ct_diff; Real pt, pt_old, ptm, dpt_dct, s1; - s1 = Sa / gsw_ups; + s1 = Sa / Psu2Gpkg; a5ct = a5 * Ct; b3ct = b3 * Ct; @@ -302,7 +298,7 @@ class Teos10Eos { pt = pt_old - ct_diff * dpt_dct; ptm = 0.5 * (pt + pt_old); - dpt_dct = -gsw_cp0 / ((ptm + gsw_t0) * calcGibbsDerivPt0Pt0(Sa, ptm)); + dpt_dct = -Cp0Sw / ((ptm + TkFrz) * calcGibbsDerivPt0Pt0(Sa, ptm)); pt = pt_old - ct_diff * dpt_dct; ct_diff = calcCtFromPt(Sa, pt) - Ct; @@ -312,11 +308,9 @@ class Teos10Eos { /// Calculate Conservative Temmperature from Potential Temp KOKKOS_FUNCTION Real calcCtFromPt(Real Sa, Real Pt) const { - constexpr Real gsw_cp0 = 3991.86795711963; // replace - constexpr Real gsw_sfac = 0.0248826675584615; // replace Real x2, x, y, pot_enthalpy; - x2 = gsw_sfac * Sa; + x2 = Sfac * Sa; x = Kokkos::sqrt(x2); y = Pt * 0.025e0; /*! normalize for F03 and F08 */ pot_enthalpy = @@ -354,7 +348,7 @@ class Teos10Eos { 9.987880382780322e0 * y) * y)))))); - return (pot_enthalpy / gsw_cp0); + return (pot_enthalpy / Cp0Sw); } private: From 314718987942f86485cce54c77cd0b96ae5fbaf9 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Mon, 9 Mar 2026 15:54:50 -0700 Subject: [PATCH 3/5] drafted tests for conversion functions - DEBUG --- components/omega/test/ocn/EosTest.cpp | 49 +++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index d5655033d82a..46956364dbaa 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -690,6 +690,53 @@ void checkValueGswcN2() { return; } +/// +void checkValueGswcCtFromPt() { + /// Get mesh and coordinate info + // const auto Mesh = HorzMesh::getDefault(); + // const auto VCoord = VertCoord::getDefault(); + // VCoord->NVertLayers = NVertLayers; + // I4 NCellsAll = Mesh->NCellsAll; + /// Get Eos instance to test + Eos *TestEos = Eos::getInstance(); + TestEos->EosChoice = EosType::Teos10Eos; + + /// Get Ct reference from GSW-C library + Real CtExpValue = gsw_ct_from_pt(Sa, Ct); + /// Get Ct from our TEOS + Real CtTeos = TestEos->calcCtFromPt(Sa, Ct); + /// Check the value against the expected TEOS-10 value + bool Check = isApprox(CtTeos, CtExpValue, RTol); + if (!Check) { + ABORT_ERROR("checkValueGswcCtfromPt: Ct calc FAIL, expected {}, got {}", + CtExpValue, CtTeos); + } + return; +} + +void checkValueGswcPtFromCt() { + /// Get mesh and coordinate info + // const auto Mesh = HorzMesh::getDefault(); + // const auto VCoord = VertCoord::getDefault(); + // VCoord->NVertLayers = NVertLayers; + // I4 NCellsAll = Mesh->NCellsAll; + /// Get Eos instance to test + Eos *TestEos = Eos::getInstance(); + TestEos->EosChoice = EosType::Teos10Eos; + + /// Get Ct reference from GSW-C library + Real PtExpValue = gsw_pt_from_ct(Sa, Ct); + /// Get Ct from our TEOS + Real PtTeos = TestEos->calcPtFromCt(Sa, Ct); + /// Check the value against the expected TEOS-10 value + bool Check = isApprox(PtTeos, PtExpValue, RTol); + if (!Check) { + ABORT_ERROR("checkValueGswcPtfromCt: Pt calc FAIL, expected {}, got {}", + PtExpValue, PtTeos); + } + return; +} + // the main tests (all in one to have the same log): // Single value test: // --> test calls the external GSW-C library @@ -709,6 +756,8 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { checkValueGswcSpecVol(); checkValueGswcN2(); + checkValueGswcCtFromPt(); + checkValueGswcPtFromCt(); testEosLinear(); testEosLinearDisplaced(); From 66ce849b7bbd64c70f67088e87f7b8f68ea85053 Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Tue, 10 Mar 2026 12:44:42 -0700 Subject: [PATCH 4/5] added the tests for single cell value and added .cpp declarations --- components/omega/src/ocn/Eos.cpp | 15 +++++++++++++++ components/omega/src/ocn/Eos.h | 7 +++++++ components/omega/test/ocn/EosTest.cpp | 18 ++++++------------ 3 files changed, 28 insertions(+), 12 deletions(-) diff --git a/components/omega/src/ocn/Eos.cpp b/components/omega/src/ocn/Eos.cpp index 1728595939cb..b92d2c6be73e 100644 --- a/components/omega/src/ocn/Eos.cpp +++ b/components/omega/src/ocn/Eos.cpp @@ -270,6 +270,21 @@ void Eos::computeBruntVaisalaFreqSq(const Array2DReal &ConservTemp, } } +Real Eos::calcPtFromCt(const Real &Sa, const Real &Ct) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcPtFromCt(Sa, Ct); + } else if (EosChoice == EosType::LinearEos) { + return Ct; + } +} + +Real Eos::calcCtFromPt(const Real &Sa, const Real &Pt) const { + if (EosChoice == EosType::Teos10Eos) { + return ComputeSpecVolTeos10.calcCtFromPt(Sa, Pt); + } else if (EosChoice == EosType::LinearEos) { + return Pt; + } +} /// Define IO fields and metadata for output void Eos::defineFields() { diff --git a/components/omega/src/ocn/Eos.h b/components/omega/src/ocn/Eos.h index d7f8f8ae2db3..51f39eb7e4aa 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -680,6 +680,13 @@ class Eos { const Array2DReal &AbsSalinity, const Array2DReal &Pressure, const Array2DReal &SpecVol); + + /// Compute Pot Temp from Consev Temp for ONE cell + Real calcPtFromCt(const Real &AbsSalinity, const Real &ConservTemp) const; + + /// Compute Consev Temp from Pot Temp for ONE cell + Real calcCtFromPt(const Real &AbsSalinity, const Real &PotTemp) const; + /// Initialize EOS from config and mesh static void init(); diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index 46956364dbaa..e6459224d009 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -692,11 +692,6 @@ void checkValueGswcN2() { /// void checkValueGswcCtFromPt() { - /// Get mesh and coordinate info - // const auto Mesh = HorzMesh::getDefault(); - // const auto VCoord = VertCoord::getDefault(); - // VCoord->NVertLayers = NVertLayers; - // I4 NCellsAll = Mesh->NCellsAll; /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::Teos10Eos; @@ -705,21 +700,18 @@ void checkValueGswcCtFromPt() { Real CtExpValue = gsw_ct_from_pt(Sa, Ct); /// Get Ct from our TEOS Real CtTeos = TestEos->calcCtFromPt(Sa, Ct); - /// Check the value against the expected TEOS-10 value + /// Check the produced value against the GSW-C value bool Check = isApprox(CtTeos, CtExpValue, RTol); if (!Check) { ABORT_ERROR("checkValueGswcCtfromPt: Ct calc FAIL, expected {}, got {}", CtExpValue, CtTeos); } + /// LOG_INFO("checkValueGswcCtfromPt: Ct calc PASS, expected {}, got {}", + /// CtExpValue, CtTeos); return; } void checkValueGswcPtFromCt() { - /// Get mesh and coordinate info - // const auto Mesh = HorzMesh::getDefault(); - // const auto VCoord = VertCoord::getDefault(); - // VCoord->NVertLayers = NVertLayers; - // I4 NCellsAll = Mesh->NCellsAll; /// Get Eos instance to test Eos *TestEos = Eos::getInstance(); TestEos->EosChoice = EosType::Teos10Eos; @@ -728,12 +720,14 @@ void checkValueGswcPtFromCt() { Real PtExpValue = gsw_pt_from_ct(Sa, Ct); /// Get Ct from our TEOS Real PtTeos = TestEos->calcPtFromCt(Sa, Ct); - /// Check the value against the expected TEOS-10 value + /// Check the produced value against the GSW-C value bool Check = isApprox(PtTeos, PtExpValue, RTol); if (!Check) { ABORT_ERROR("checkValueGswcPtfromCt: Pt calc FAIL, expected {}, got {}", PtExpValue, PtTeos); } + /// LOG_INFO("checkValueGswcPtfromCt: Pt calc PASS, expected {}, got {}", + /// PtExpValue, PtTeos); return; } From 8aba3e3e07ce714db6345f7ce26d1f44fe61e06a Mon Sep 17 00:00:00 2001 From: Alice Barthel Date: Tue, 17 Mar 2026 14:59:12 -0700 Subject: [PATCH 5/5] added a note about conversion function in the devGuide doc --- components/omega/doc/devGuide/EOS.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/components/omega/doc/devGuide/EOS.md b/components/omega/doc/devGuide/EOS.md index 0762727a25b5..727b819230ef 100644 --- a/components/omega/doc/devGuide/EOS.md +++ b/components/omega/doc/devGuide/EOS.md @@ -66,6 +66,14 @@ volume arrays, do Eos.computeBruntVaisalaFreqSq(ConservTemp, AbsSalinity, Pressure, SpecVol); ``` +## Helper functions for conversion + +To provide the surface in-situ temperature that the coupler requires, use the helper function that converts conservative temperature (the state variable when using Teos-10) into potential temperature, given the absolute salinity and pressure arrays: + +```c++ +Eos.calcPtFromCt(ConsrvTemp, AbsSalinity, Pressure); +``` + ## Removal of Eos To clear the Eos instance do: