diff --git a/components/omega/src/ocn/KPPConstants.h b/components/omega/src/ocn/KPPConstants.h index 41d2ee9f7de9..5b3f802c7cfc 100644 --- a/components/omega/src/ocn/KPPConstants.h +++ b/components/omega/src/ocn/KPPConstants.h @@ -49,11 +49,11 @@ constexpr Real AMoS = -28.8617393793_Real; constexpr Real CMoS = 98.9545350148_Real; /// Surface value of the momentum shape function (dimensionless) -constexpr Real HuOn = 0.03; +constexpr Real HuOn = 0.03_Real; /// Floor on friction velocity (m/s). Keeps the turbulent velocity scales and /// the Monin-Obukhov length finite in near-calm conditions. -constexpr Real MinUStar = 1.0e-4; +constexpr Real MinUStar = 1.0e-4_Real; // =========================================================================== // Numerical Safeguards @@ -81,16 +81,16 @@ constexpr Real MinStokesDrift = 1.0e-8_Real; /// /// Critical bulk Richardson number defining the OSBL base /// REFERENCES: Large et al. (1994) Eq. (21) -constexpr Real CriticalRi = 0.25; +constexpr Real CriticalRi = 0.25_Real; /// Minimum OSBL depth under sea ice (m), applied above IceSuppressThresh -constexpr Real MinOSBLUnderIce = 5.0; +constexpr Real MinOSBLUnderIce = 5.0_Real; /// Ice fraction above which Langmuir enhancement is disabled -constexpr Real IceFracThresh = 0.05; +constexpr Real IceFracThresh = 0.05_Real; /// Ice fraction above which the minimum OSBL depth is enforced -constexpr Real IceSuppressThresh = 0.15; +constexpr Real IceSuppressThresh = 0.15_Real; // ========================================================================== // Surface Layer Parameters @@ -99,11 +99,11 @@ constexpr Real IceSuppressThresh = 0.15; /// Surface layer extent as a fraction of the trial OSBL depth (epsilon in /// Large et al. 1994). Reference values entering the bulk Richardson number /// are averaged over the top SurfaceLayerExtent * d of the column. -constexpr Real SurfaceLayerExtent = 0.1; +constexpr Real SurfaceLayerExtent = 0.1_Real; /// Empirical convective velocity coefficient in the turbulent velocity scale /// w_s = (u*^3 + ConvectiveVelFac * max(-B_0,0) * h)^(1/3) -constexpr Real ConvectiveVelFac = 0.35; +constexpr Real ConvectiveVelFac = 0.35_Real; // ========================================================================== // KPP Shape and Stability Functions @@ -125,10 +125,10 @@ constexpr Real ConvectiveVelFac = 0.35; /// @return Dimensionless shape value KOKKOS_INLINE_FUNCTION Real kppShapeMomentum(Real Sigma) { - Sigma = Kokkos::fmax(-1.0, Kokkos::fmin(0.0, Sigma)); + Sigma = Kokkos::fmax(-1.0_Real, Kokkos::fmin(0.0_Real, Sigma)); const Real SigmaMu = -Sigma; - return SigmaMu * (1.0 - SigmaMu) * (1.0 - SigmaMu); + return SigmaMu * (1.0_Real - SigmaMu) * (1.0_Real - SigmaMu); } /// @brief Matched KPP gradient shape function. @@ -139,11 +139,11 @@ Real kppShapeMomentum(Real Sigma) { /// preserving SimpleShapes behavior when ShapeAtBase is zero. KOKKOS_INLINE_FUNCTION Real kppShapeMatched(Real Sigma, Real ShapeAtBase) { - Sigma = Kokkos::fmax(-1.0, Kokkos::fmin(0.0, Sigma)); + Sigma = Kokkos::fmax(-1.0_Real, Kokkos::fmin(0.0_Real, Sigma)); const Real SigmaMu = -Sigma; - const Real Simple = SigmaMu * (1.0 - SigmaMu) * (1.0 - SigmaMu); - const Real Smooth = SigmaMu * SigmaMu * (3.0 - 2.0 * SigmaMu); + const Real Simple = SigmaMu * (1.0_Real - SigmaMu) * (1.0_Real - SigmaMu); + const Real Smooth = SigmaMu * SigmaMu * (3.0_Real - 2.0_Real * SigmaMu); return Simple + ShapeAtBase * Smooth; } @@ -178,10 +178,10 @@ Real kppPhiInvMomentum(Real Zeta) { /// @return Dimensionless shape value KOKKOS_INLINE_FUNCTION Real kppShapeScalar(Real Sigma) { - Sigma = Kokkos::fmax(-1.0, Kokkos::fmin(0.0, Sigma)); + Sigma = Kokkos::fmax(-1.0_Real, Kokkos::fmin(0.0_Real, Sigma)); const Real SigmaMu = -Sigma; - return SigmaMu * (1.0 - SigmaMu) * (1.0 - SigmaMu); + return SigmaMu * (1.0_Real - SigmaMu) * (1.0_Real - SigmaMu); } /// @brief phi_s^{-1}(zeta) - Inverse scalar Monin-Obukhov stability function @@ -216,7 +216,7 @@ Real kppPhiInvScalar(Real Zeta) { KOKKOS_INLINE_FUNCTION Real kppSurfaceMomentumScale(Real Sigma) { // Linear decay from HuOn at the surface to zero at the OSBL base - return HuOn * (1.0 + Sigma); + return HuOn * (1.0_Real + Sigma); } // ========================================================================== @@ -245,12 +245,12 @@ Real kppSurfaceMomentumScale(Real Sigma) { /// @return Stokes drift velocity scale (m/s) KOKKOS_INLINE_FUNCTION Real estimateStokesDriftSL(Real Wind10m, Real HBL) { - Wind10m = Kokkos::fmax(0.0, Wind10m); - HBL = Kokkos::fmax(1.0, HBL); + Wind10m = Kokkos::fmax(0.0_Real, Wind10m); + HBL = Kokkos::fmax(1.0_Real, HBL); - const Real UStokes = 0.016 * Wind10m; + const Real UStokes = 0.016_Real * Wind10m; - return Kokkos::fmax(0.0, UStokes); + return Kokkos::fmax(0.0_Real, UStokes); } /// @brief Turbulent Langmuir number from friction velocity and Stokes drift @@ -277,18 +277,18 @@ Real computeLangmuirNumber(Real UStar, Real UStokes) { KOKKOS_INLINE_FUNCTION Real computeLangmuirEnhancement(Real Wind10m, Real UStar, Real HBL) { UStar = Kokkos::fmax(MinUStar, UStar); - Wind10m = Kokkos::fmax(0.0, Wind10m); - HBL = Kokkos::fmax(1.0, HBL); + Wind10m = Kokkos::fmax(0.0_Real, Wind10m); + HBL = Kokkos::fmax(1.0_Real, HBL); const Real UStokes = estimateStokesDriftSL(Wind10m, HBL); const Real La = computeLangmuirNumber(UStar, UStokes); // R_L = sqrt(1 + 0.5/La^2); La is floored at 0.5 so that the weak-wave // limit returns the unenhanced scales rather than diverging. - const Real LaInv = 1.0 / Kokkos::fmax(0.5, La); - const Real RL = Kokkos::sqrt(1.0 + 0.5 * LaInv * LaInv); + const Real LaInv = 1.0_Real / Kokkos::fmax(0.5_Real, La); + const Real RL = Kokkos::sqrt(1.0_Real + 0.5_Real * LaInv * LaInv); - return Kokkos::fmin(2.0, Kokkos::fmax(1.0, RL)); + return Kokkos::fmin(2.0_Real, Kokkos::fmax(1.0_Real, RL)); } // ========================================================================== @@ -325,7 +325,7 @@ Real constrainOSBLDepth(Real HOSBL, Real LayerThickness, Real WaterDepth, } // Upper bound: cannot exceed water depth - HOSBL = Kokkos::fmin(HOSBL, WaterDepth * 0.95); + HOSBL = Kokkos::fmin(HOSBL, WaterDepth * 0.95_Real); return HOSBL; } diff --git a/components/omega/src/ocn/VertMix.cpp b/components/omega/src/ocn/VertMix.cpp index cb0a563e1fc8..5a09695ceff4 100644 --- a/components/omega/src/ocn/VertMix.cpp +++ b/components/omega/src/ocn/VertMix.cpp @@ -235,17 +235,12 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, OMEGA_SCOPE(LocBackVisc, BackVisc); OMEGA_SCOPE(MinLayerCell, VCoord->MinLayerCell); OMEGA_SCOPE(MaxLayerCell, VCoord->MaxLayerCell); - const Real LocConvDiff = LocComputeVertMixConv.ConvDiff; - const Real LocConvTriggerBVF = LocComputeVertMixConv.ConvTriggerBVF; - Array1DI4 KPPBoundaryLayerIndex("VertMix-KPPBoundaryLayerIndex", - Mesh->NCellsSize); - deepCopy(KPPBoundaryLayerIndex, -1); + Array1DI4 KPPBoundaryLayerIndex; KPPMix *KPPInstance = KPPMix::getInstance(); const bool LocKPPEnabled = (KPPInstance && KPPInstance->Enabled); if (LocKPPEnabled) { - deepCopy(KPPBoundaryLayerIndex, KPPInstance->OSBLDepthIndex); + KPPBoundaryLayerIndex = KPPInstance->OSBLDepthIndex; } - OMEGA_SCOPE(LocKPPBoundaryLayerIndex, KPPBoundaryLayerIndex); /// First, initialize VertDiff and VertVisc to background values parallelForOuter( @@ -358,7 +353,7 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, parallelFor( "VertMix-KPP", {Mesh->NCellsAll, NVertLayers + 1}, KOKKOS_LAMBDA(I4 ICell, I4 K) { - if (K <= LocKPPBoundaryLayerIndex(ICell) + 1) { + if (K <= KPPBoundaryLayerIndex(ICell) + 1) { if (KPPMergeMode == 1) { LocVertDiff(ICell, K) = LocKPPVertDiff(ICell, K); LocVertVisc(ICell, K) = LocKPPVertVisc(ICell, K); @@ -381,19 +376,9 @@ void VertMix::computeVertMix(const Array2DReal &NormalVelocity, parallelForInner( Team, KRange, INNER_LAMBDA(int KChunk) { - const I4 KStart = chunkStart(KChunk, KMin); - const I4 KLen = chunkLength(KChunk, KStart, KMax); - for (int KVec = 0; KVec < KLen; ++KVec) { - const I4 K = KStart + KVec; - const bool ApplyConv = - (!LocKPPEnabled) || - (K > LocKPPBoundaryLayerIndex(ICell) + 1); - if (ApplyConv && - BruntVaisalaFreqSq(ICell, K) < LocConvTriggerBVF) { - LocVertDiff(ICell, K) += LocConvDiff; - LocVertVisc(ICell, K) += LocConvDiff; - } - } + LocComputeVertMixConv(LocVertDiff, LocVertVisc, ICell, + KChunk, BruntVaisalaFreqSq, + KPPBoundaryLayerIndex); }); }); } diff --git a/components/omega/src/ocn/VertMix.h b/components/omega/src/ocn/VertMix.h index 03c97e852132..34d2f8d7896a 100644 --- a/components/omega/src/ocn/VertMix.h +++ b/components/omega/src/ocn/VertMix.h @@ -36,7 +36,8 @@ class ConvectiveMix { KOKKOS_FUNCTION void operator()(Array2DReal VertDiff, Array2DReal VertVisc, I4 ICell, I4 KChunk, - const Array2DReal &BruntVaisalaFreqSq) const { + const Array2DReal &BruntVaisalaFreqSq, + const Array1DI4 &KPPBoundaryLayerIndex) const { const I4 KStart = chunkStart(KChunk, MinLayerCell(ICell) + 1); const I4 KLen = chunkLength(KChunk, KStart, MaxLayerCell(ICell)); @@ -44,7 +45,12 @@ class ConvectiveMix { for (int KVec = 0; KVec < KLen; ++KVec) { const I4 K = KStart + KVec; - if (BruntVaisalaFreqSq(ICell, K) < ConvTriggerBVF) { + // KPP is enabled when the index array is not empty + const bool KPPEnabled = KPPBoundaryLayerIndex.data(); + const bool ApplyConv = + (!KPPEnabled) || (K > KPPBoundaryLayerIndex(ICell) + 1); + + if (ApplyConv && BruntVaisalaFreqSq(ICell, K) < ConvTriggerBVF) { VertDiff(ICell, K) += ConvDiff; VertVisc(ICell, K) += ConvDiff; } diff --git a/components/omega/test/ocn/VertMixTest.cpp b/components/omega/test/ocn/VertMixTest.cpp index 7ed3601c068e..e64a70019642 100644 --- a/components/omega/test/ocn/VertMixTest.cpp +++ b/components/omega/test/ocn/VertMixTest.cpp @@ -553,11 +553,12 @@ void testConvVertMix() { /// Compute only convective vertical viscosity and diffusivity OMEGA_SCOPE(ComputeVertMixConv, TestVertMix->ComputeVertMixConv); + Array1DI4 KPPBoundaryLayerIndex; // Empty, so KPP is disabled for this test parallelFor( "ApplyVertMixConv", {Mesh->NCellsAll, NChunks}, KOKKOS_LAMBDA(I4 ICell, I4 KChunk) { ComputeVertMixConv(VertDiffOut, VertViscOut, ICell, KChunk, - BruntVaisalaFreqSqIn); + BruntVaisalaFreqSqIn, KPPBoundaryLayerIndex); }); const auto &MinLayerCell = VCoord->MinLayerCell;