Add the Eos function to convert CT into PotT

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:

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() {
 

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();
 

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

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();