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: 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 06c9695f7267..51f39eb7e4aa 100644 --- a/components/omega/src/ocn/Eos.h +++ b/components/omega/src/ocn/Eos.h @@ -238,6 +238,119 @@ 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 { + Real x2 = 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 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 / Psu2Gpkg; + + 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 = -Cp0Sw / ((ptm + TkFrz) * 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 { + Real x2, x, y, pot_enthalpy; + + x2 = 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 / Cp0Sw); + } + private: Array1DI4 MinLayerCell; Array1DI4 MaxLayerCell; @@ -567,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/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 diff --git a/components/omega/test/ocn/EosTest.cpp b/components/omega/test/ocn/EosTest.cpp index d5655033d82a..e6459224d009 100644 --- a/components/omega/test/ocn/EosTest.cpp +++ b/components/omega/test/ocn/EosTest.cpp @@ -690,6 +690,47 @@ void checkValueGswcN2() { return; } +/// +void checkValueGswcCtFromPt() { + /// 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 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 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 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; +} + // the main tests (all in one to have the same log): // Single value test: // --> test calls the external GSW-C library @@ -709,6 +750,8 @@ void eosTest(const std::string &MeshFile = "OmegaMesh.nc") { checkValueGswcSpecVol(); checkValueGswcN2(); + checkValueGswcCtFromPt(); + checkValueGswcPtFromCt(); testEosLinear(); testEosLinearDisplaced();