Add temporary global constants header file

components/omega/configs/Default.yml

@@ -28,6 +28,7 @@ Omega:
   WindStress:
     InterpType: Isotropic
   VertCoord:
+    Density0: 1026.0
     MovementWeightType: Uniform
   Tendencies:
     ThicknessFluxTendencyEnable: true
@@ -40,7 +41,6 @@ Omega:
     ViscDel4: 1.2e11
     DivFactor: 1.0
     WindForcingTendencyEnable: false
-    Density0: 1026.0
     BottomDragTendencyEnable: false
     BottomDragCoeff: 0.0
     TracerHorzAdvTendencyEnable: true

components/omega/src/CMakeLists.txt

@@ -12,6 +12,7 @@ target_include_directories(
     ${OMEGA_SOURCE_DIR}/src/timeStepping
     ${E3SM_ROOT}/share/timing
     ${E3SM_ROOT}/share/pacer
+    ${E3SM_ROOT}/share/util_cxx
     ${PIOC_SOURCE_DIR}
     ${PIOC_BINARY_DIR}
     ${Parmetis_INCLUDE_DIRS}

components/omega/src/ocn/CustomTendencyTerms.cpp

@@ -8,6 +8,7 @@
 
 #include "CustomTendencyTerms.h"
 #include "Config.h"
+#include "GlobalConstants.h"
 #include "TimeStepper.h"
 
 namespace OMEGA {
@@ -79,11 +80,9 @@ void ManufacturedSolution::init() {
    R8 H0                 = DefHorzMesh->BottomDepthH(0);
 
    // Define and compute common constants
-   R8 Grav    = 9.80665_Real;                          // Gravity acceleration
-   R8 Pii     = 3.141592653589793_Real;                // Pi
-   R8 Kx      = 2.0_Real * Pii / WavelengthX;          // Wave in X-dir
-   R8 Ky      = 2.0_Real * Pii / WavelengthY;          // Wave in Y-dir
-   R8 AngFreq = sqrt(H0 * Grav * (Kx * Kx + Ky * Ky)); // Angular frequency
+   R8 Kx      = TwoPi / WavelengthX;                      // Wave in X-dir
+   R8 Ky      = TwoPi / WavelengthY;                      // Wave in Y-dir
+   R8 AngFreq = sqrt(H0 * Gravity * (Kx * Kx + Ky * Ky)); // Angular frequency
 
    // Assign constants for thickness tendency function
    ManufacturedThickTend.H0      = H0;
@@ -93,7 +92,7 @@ void ManufacturedSolution::init() {
    ManufacturedThickTend.AngFreq = AngFreq;
 
    // Assign constants for velocity tendency function
-   ManufacturedVelTend.Grav    = Grav;
+   ManufacturedVelTend.Grav    = Gravity;
    ManufacturedVelTend.Eta0    = Amplitude;
    ManufacturedVelTend.Kx      = Kx;
    ManufacturedVelTend.Ky      = Ky;
@@ -160,7 +159,7 @@ void ManufacturedSolution::ManufacturedVelocityTendency::operator()(
    Array1DReal YEdge     = Mesh->YEdge;
    Array1DReal AngleEdge = Mesh->AngleEdge;
 
-   OMEGA_SCOPE(LocGrav, Grav);
+   OMEGA_SCOPE(LocGrav, Gravity);
    OMEGA_SCOPE(LocEta0, Eta0);
    OMEGA_SCOPE(LocKx, Kx);
    OMEGA_SCOPE(LocKy, Ky);

components/omega/src/ocn/Eos.h

@@ -12,6 +12,7 @@
 
 #include "AuxiliaryState.h"
 #include "Config.h"
+#include "GlobalConstants.h"
 #include "HorzMesh.h"
 #include "MachEnv.h"
 #include "OmegaKokkos.h"
@@ -205,8 +206,7 @@ class Teos10Eos {
    KOKKOS_FUNCTION Real calcDelta(const Array2DReal &SpecVolPCoeffs, const I4 K,
                                   const Real P) const {
 
-      constexpr Real Pu = 1e4;
-      Real Pp           = P / Pu;
+      Real Pp = P * Pa2Db;
 
       Real Delta = ((((SpecVolPCoeffs(5, K) * Pp + SpecVolPCoeffs(4, K)) * Pp +
                       SpecVolPCoeffs(3, K)) *
@@ -221,14 +221,13 @@ class Teos10Eos {
 
    /// Calculate reference profile for TEOS-10
    KOKKOS_FUNCTION Real calcRefProfile(Real P) const {
-      constexpr Real Pu  = 1e4;
       constexpr Real V00 = -4.4015007269e-05;
       constexpr Real V01 = 6.9232335784e-06;
       constexpr Real V02 = -7.5004675975e-07;
       constexpr Real V03 = 1.7009109288e-08;
       constexpr Real V04 = -1.6884162004e-08;
       constexpr Real V05 = 1.9613503930e-09;
-      Real Pp            = P / Pu;
+      Real Pp            = P * Pa2Db;
 
       Real V0 =
           (((((V05 * Pp + V04) * Pp + V03) * Pp + V02) * Pp + V01) * Pp + V00) *
@@ -244,9 +243,9 @@ class Teos10Eos {
 class LinearEos {
  public:
    /// Coefficients for LinearEos (overwritten by config file if set there)
-   Real DRhodT  = -0.2;   ///< Thermal expansion coefficient (kg m^-3 degC^-1)
-   Real DRhodS  = 0.8;    ///< Haline contraction coefficient (kg m^-3)
-   Real RhoT0S0 = 1000.0; ///< Reference density (kg m^-3) at (T,S)=(0,0)
+   Real DRhodT  = -0.2;  ///< Thermal expansion coefficient (kg m^-3 degC^-1)
+   Real DRhodS  = 0.8;   ///< Haline contraction coefficient (kg m^-3)
+   Real RhoT0S0 = RhoFw; ///< Reference density (kg m^-3) at (T,S)=(0,0)
 
    /// constructor declaration
    LinearEos();

components/omega/src/ocn/GlobalConstants.h

@@ -0,0 +1,100 @@
+#ifndef OMEGA_GLOBALCONSTANTS_H
+#define OMEGA_GLOBALCONSTANTS_H
+//===-- ocn/GlobalConstants.h - Global Constants --------------------*- C++
+//-*-===//
+//
+/// \file
+/// \brief Contains global constants for use across Omega
+///
+/// This header defines constants to be used across Omega (to be replaced by a
+/// shared E3SM constants file when that becomes available)
+//
+//===----------------------------------------------------------------------===//
+
+#include "DataTypes.h"
+#include "pcd_const.h"
+#include <cmath>
+
+namespace OMEGA {
+
+// Mathematical constants
+constexpr Real Pi = pcd::pi; // Pi (from Physical Constants Dictionary)
+constexpr Real TwoPi =
+    2.0 * pcd::pi; // 2*Pi (from Physical Constants Dictionary)
+constexpr Real Radian =
+    pcd::radian; // Radians in one degree (from Physical Constants Dictionary)
+constexpr Real Degree =
+    pcd::degree; // Degrees in one radian (from Physical Constants Dictionary)
+constexpr Real SqrtTwo =
+    pcd::square_root_of_2; // Square root of 2 (from Physical Constants
+                           // Dictionary)
+constexpr Real SqrtThree =
+    pcd::square_root_of_3; // Square root of 3 (from Physical Constants
+                           // Dictionary)
+
+// Earth constants
+constexpr Real Gravity =
+    pcd::standard_acceleration_of_gravity; // Acceleration due to gravity ~
+                                           // m/s^2 (from Physical Constants
+                                           // Dictionary)
+constexpr Real CDay = 86400.0;             // Seconds in a calendar day ~ sec
+constexpr Real SDay = 86164.0;             // Seconds in a sidereal day ~ sec
+constexpr Real Omega =
+    2.0 * Pi / SDay; // Angular velocity of the Earth ~ rad/sec
+constexpr Real REarth = pcd::mean_radius; // Mean radius of the Earth ~ m (from
+                                          // Physical Constants Dictionary)
+
+/// Physical constants
+constexpr Real TkTrip =
+    pcd::water_triple_point_temperature; // Triple point of fresh water ~ K
+                                         // (from Physical Constants Dictionary)
+constexpr Real TkFrz   = 273.15;         // Freezing point of fresh water ~ K
+constexpr Real TkFrzSw = TkFrz - 1.8;    // Freezing point of seawater ~ K
+constexpr Real RhoAir  = 1.2;            // Density of air ~ kg/m^3
+constexpr Real RhoFw   = 1.000e3;        // Density of fresh water ~ kg/m^3
+constexpr Real RhoSw   = 1.026e3;        // Density of seawater ~ kg/m^3
+constexpr Real RhoIce  = 0.917e3;        // Density of ice ~ kg/m^3
+constexpr Real CpAir   = 1.005e3; // Specific heat capacity of air ~ J/(kg*K)
+constexpr Real CpFw =
+    4.188e3; // Specific heat capacity of fresh water ~ J/(kg*K)
+constexpr Real CpSw  = 3.996e3; // Specific heat capacity of seawater ~ J/(kg*K)
+constexpr Real CpIce = 2.108e3; // Specific heat capacity of ice ~ J/(kg*K)
+constexpr Real LatIce    = 3.337e5; // Latent heat of fusion ~ J/kg
+constexpr Real LatVap    = 2.501e6; // Latent heat of vaporization ~ J/kg
+constexpr Real LatSub    = LatIce + LatVap; // Latent heat of sublimation ~ J/kg
+constexpr Real CondIce   = 2.1;   // Universal gas constant ~ J/(mol*K)
+constexpr Real OcnRefSal = 34.7;  // Reference ocean salinity ~ psu
+constexpr Real IceRefSal = 4.0;   // Reference ice salinity ~ psu
+constexpr Real Sound     = 1.5e2; // Speed of sound ~ m/s
+constexpr Real VonKar    = 0.4;   // Von Karman constant ~ dimensionless
+constexpr Real Emiss     = 1.0;   // Emissivity ~ dimensionless
+constexpr Real AtmRefP =
+    pcd::standard_atmosphere; // Reference atmospheric pressure ~ Pa (from
+                              // Physical Constants Dictionary)
+
+// Conversion factors
+constexpr Real Sec2Day = 1.0 / 86400.0; // Seconds to days
+constexpr Real Day2Sec = 86400.0;       // Days to seconds
+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 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
+constexpr Real Salt2SvPpt = 1.0e-9;    // Salt flux (kg/s) to Sv*ppt
+constexpr Real Salt2MmDay = 3.1536e+5; // Salt flux to water flux (mm/day)
+constexpr Real Db2Pa      = 1.0e4;     // Decibar to Pascal
+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)
+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)
+
+} // namespace OMEGA
+#endif

components/omega/src/ocn/Tendencies.cpp

@@ -179,9 +179,6 @@ void Tendencies::readTendConfig(
    CHECK_ERROR_ABORT(
        Err, "Tendencies: WindForcingTendencyEnable not found in TendConfig");
 
-   Err += TendConfig->get("Density0", this->WindForcing.SaltWaterDensity);
-   CHECK_ERROR_ABORT(Err, "Tendencies: Density0 not found in TendConfig");
-
    Err += TendConfig->get("BottomDragTendencyEnable", this->BottomDrag.Enabled);
    CHECK_ERROR_ABORT(
        Err, "Tendencies: BottomDragTendencyEnable not found in TendConfig");

components/omega/src/ocn/TendencyTerms.cpp

@@ -46,7 +46,7 @@ VelocityHyperDiffOnEdge::VelocityHyperDiffOnEdge(const HorzMesh *Mesh)
       MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {}
 
 WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh)
-    : Enabled(false), EdgeMask(Mesh->EdgeMask) {}
+    : Enabled(false), LocRhoSw(RhoSw), EdgeMask(Mesh->EdgeMask) {}
 
 BottomDragOnEdge::BottomDragOnEdge(const HorzMesh *Mesh,
                                    const VertCoord *VCoord)

components/omega/src/ocn/TendencyTerms.h

@@ -11,6 +11,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "AuxiliaryState.h"
+#include "GlobalConstants.h"
 #include "HorzMesh.h"
 #include "MachEnv.h"
 #include "OceanState.h"
@@ -166,14 +167,13 @@ class SSHGradOnEdge {
 
       for (int KVec = 0; KVec < VecLength; ++KVec) {
          const I4 K = KStart + KVec;
-         Tend(IEdge, K) -= EdgeMask(IEdge, K) * Grav *
+         Tend(IEdge, K) -= EdgeMask(IEdge, K) * Gravity *
                            (SshCell(ICell1, K) - SshCell(ICell0, K)) *
                            InvDcEdge;
       }
    }
 
  private:
-   Real Grav = 9.80665_Real;
    Array2DI4 CellsOnEdge;
    Array1DReal DcEdge;
    Array2DReal EdgeMask;
@@ -281,7 +281,7 @@ class VelocityHyperDiffOnEdge {
 class WindForcingOnEdge {
  public:
    bool Enabled;
-   Real SaltWaterDensity;
+   Real LocRhoSw;
 
    /// constructor declaration
    WindForcingOnEdge(const HorzMesh *Mesh);
@@ -296,7 +296,7 @@ class WindForcingOnEdge {
 
          const Real InvThickEdge = 1._Real / LayerThickEdge(IEdge, K);
          Tend(IEdge, K) += EdgeMask(IEdge, K) * InvThickEdge *
-                           NormalStressEdge(IEdge) / SaltWaterDensity;
+                           NormalStressEdge(IEdge) / LocRhoSw;
       }
    }
 

components/omega/src/ocn/VertCoord.cpp

@@ -10,6 +10,7 @@
 #include "VertCoord.h"
 #include "Dimension.h"
 #include "Field.h"
+#include "GlobalConstants.h"
 #include "IO.h"
 #include "IOStream.h"
 #include "OmegaKokkos.h"
@@ -216,14 +217,14 @@ void VertCoord::completeSetup(Config *Options //< [in] configuration options
    MinLayerCellH = createHostMirrorCopy(MinLayerCell);
    BottomDepthH  = createHostMirrorCopy(BottomDepth);
 
-   // Fetch reference desnity from Config
-   Config TendConfig("Tendencies");
+   // Fetch reference density from Config
+   Config VCoordConfig("VertCoord");
    Err.reset();
-   Err += Options->get(TendConfig);
-   CHECK_ERROR_ABORT(Err, "VertCoord: Tendencies group not found in Config");
+   Err += Options->get(VCoordConfig);
+   CHECK_ERROR_ABORT(Err, "VertCoord: VertCoord group not found in Config");
 
-   Err += TendConfig.get("Density0", Rho0);
-   CHECK_ERROR_ABORT(Err, "VertCoord: Density0 not found in TendConfig");
+   Err += VCoordConfig.get("Density0", Rho0);
+   CHECK_ERROR_ABORT(Err, "VertCoord: Density0 not found in VertCoord");
 
 } // end completeSetup
 
@@ -656,8 +657,6 @@ void VertCoord::computePressure(
     const Array1DReal &SurfacePressure // [in] surface pressure
 ) {
 
-   Real Gravity = 9.80616_Real; // gravitationl acceleration
-
    OMEGA_SCOPE(LocRho0, Rho0);
    OMEGA_SCOPE(LocMinLayerCell, MinLayerCell);
    OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell);
@@ -745,8 +744,6 @@ void VertCoord::computeGeopotential(
     const Array1DReal &SelfAttractionLoading // [in] self attraction and loading
 ) {
 
-   Real Gravity = 9.80616_Real; // gravitationl acceleration
-
    OMEGA_SCOPE(LocMinLayerCell, MinLayerCell);
    OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell);
    OMEGA_SCOPE(LocGeopotMid, GeopotentialMid);
@@ -784,8 +781,6 @@ void VertCoord::computeGeopotential(
 // reductions and a parallel_for over the active layers within a column.
 void VertCoord::computeTargetThickness() {
 
-   Real Gravity = 9.80616_Real; // gravitationl acceleration
-
    OMEGA_SCOPE(LocRho0, Rho0);
    OMEGA_SCOPE(LocMinLayerCell, MinLayerCell);
    OMEGA_SCOPE(LocMaxLayerCell, MaxLayerCell);

components/omega/test/ocn/AuxiliaryStateTest.cpp

@@ -4,6 +4,7 @@
 #include "Decomp.h"
 #include "Dimension.h"
 #include "Field.h"
+#include "GlobalConstants.h"
 #include "Halo.h"
 #include "HorzMesh.h"
 #include "IO.h"
@@ -23,7 +24,7 @@
 using namespace OMEGA;
 
 struct TestSetup {
-   Real Radius = 6371220;
+   Real Radius = REarth;
 
    KOKKOS_FUNCTION Real layerThickness(Real Lon, Real Lat) const {
       return (2 + std::cos(Lon) * std::pow(std::cos(Lat), 4));