Skip to content

Commit 79f2f40

Browse files
alicebarthelalicebarthel
committed
added p,T,S limits (with CtFreezing calc) to eos Poly75t and a check on
CtFreezing calculation (against GSW-C)
1 parent e73ce82 commit 79f2f40

File tree

2 files changed

+129
-7
lines changed

2 files changed

+129
-7
lines changed

‎components/omega/src/ocn/Eos.h‎

Lines changed: 100 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,11 @@ enum class EosType {
2525
Teos10Eos /// Roquet et al. 2015 75 term expansion
2626
};
2727

28+
struct Range {
29+
Real lo;
30+
Real hi;
31+
};
32+
2833
/// TEOS10 75-term Polynomial Equation of State
2934
class Teos10Eos {
3035
public:
@@ -51,7 +56,7 @@ class Teos10Eos {
5156
/// Calculate the local specific volume polynomial pressure
5257
/// coefficients
5358
calcPCoeffs(LocSpecVolPCoeffs, KVec, ConservTemp(ICell, K),
54-
AbsSalinity(ICell, K));
59+
AbsSalinity(ICell, K), Pressure(ICell, K));
5560

5661
/// Calculate the specific volume at the given pressure
5762
/// If KDisp is 0, we use the current pressure, otherwise we use the
@@ -77,12 +82,15 @@ class Teos10Eos {
7782
/// TEOS-10 helpers
7883
/// Calculate pressure polynomial coefficients for TEOS-10
7984
KOKKOS_FUNCTION void calcPCoeffs(Array2DReal SpecVolPCoeffs, const I4 K,
80-
const Real Ct, const Real Sa) const {
85+
const Real Ct, const Real Sa, const Real P) const {
8186
constexpr Real SAu = 40.0 * 35.16504 / 35.0;
8287
constexpr Real CTu = 40.0;
8388
constexpr Real DeltaS = 24.0;
84-
Real Ss = Kokkos::sqrt((Sa + DeltaS) / SAu);
85-
Real Tt = Ct / CTu;
89+
Range Srange = calcSLimits(P);
90+
Range Trange = calcTLimits(Sa, P);
91+
Real SaInFunnel = Kokkos::clamp(Sa, Srange.lo, Srange.hi); // Salt limited to Poly75t valid range
92+
Real Ss = Kokkos::sqrt((SaInFunnel + DeltaS) / SAu);
93+
Real Tt = Kokkos::clamp(Ct, Trange.lo, Trange.hi) / CTu;
8694

8795
/// Coefficients for the polynomial expansion
8896
constexpr Real V000 = 1.0769995862e-03;
@@ -205,7 +213,8 @@ class Teos10Eos {
205213
const Real P) const {
206214

207215
constexpr Real Pu = 1e4;
208-
Real Pp = P / Pu;
216+
constexpr Real Pmax = 8000.0;
217+
Real Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range
209218

210219
Real Delta = ((((SpecVolPCoeffs(5, K) * Pp + SpecVolPCoeffs(4, K)) * Pp +
211220
SpecVolPCoeffs(3, K)) *
@@ -219,22 +228,106 @@ class Teos10Eos {
219228
}
220229

221230
/// Calculate reference profile for TEOS-10
222-
KOKKOS_FUNCTION Real calcRefProfile(Real P) const {
231+
KOKKOS_FUNCTION Real calcRefProfile(const Real P) const {
223232
constexpr Real Pu = 1e4;
224233
constexpr Real V00 = -4.4015007269e-05;
225234
constexpr Real V01 = 6.9232335784e-06;
226235
constexpr Real V02 = -7.5004675975e-07;
227236
constexpr Real V03 = 1.7009109288e-08;
228237
constexpr Real V04 = -1.6884162004e-08;
229238
constexpr Real V05 = 1.9613503930e-09;
230-
Real Pp = P / Pu;
239+
constexpr Real Pmax = 8000.0;
240+
Real Pp = Kokkos::min(P, Pmax) / Pu; // P limited to Poly75t valid range
231241

232242
Real V0 =
233243
(((((V05 * Pp + V04) * Pp + V03) * Pp + V02) * Pp + V01) * Pp + V00) *
234244
Pp;
235245
return V0;
236246
}
237247

248+
/// Calculate S limits of validity given pressure p
249+
KOKKOS_FUNCTION Range calcSLimits(const Real P) const {
250+
Real lo = 0.0;
251+
Real hi = 42.0;
252+
Real lo2 = P * 5e-3 - 2.5;
253+
Real lo3 = 30.0;
254+
255+
if (P >= 500.0 && P < 6500.0) {
256+
lo = Kokkos::max(lo, lo2);
257+
}
258+
else if (P >= 6500.0) {
259+
lo = Kokkos::max(lo, lo3);
260+
}
261+
return {lo, hi};
262+
}
263+
264+
/// Calculate T limits of validity given p,Sa
265+
KOKKOS_FUNCTION Range calcTLimits(const Real Sa, const Real P) const {
266+
Real lo = - 15.0;
267+
Real hi = 95.0;
268+
269+
if (P < 500.0) {
270+
lo = calcCtFreezing(Sa, P, Real{0.0});
271+
} else if (P < 6500.0) {
272+
lo = calcCtFreezing(Sa, Real{500.0}, Real{0.0});
273+
hi = 31.66666666666667 - P * 3.333333333333334e-3;
274+
} else {
275+
lo = calcCtFreezing(Sa, Real{500.0}, Real{0.0});
276+
hi = 10.0;
277+
}
278+
return {lo, hi};
279+
}
280+
281+
/// Calculates freezing CTemperature (polynomial error in [-5e-4,6e-4] K)
282+
KOKKOS_FUNCTION Real calcCtFreezing(const Real Sa, const Real P, const Real SaturationFract) const {
283+
constexpr Real Sso = 35.16504;
284+
constexpr Real c0 = 0.017947064327968736;
285+
constexpr Real c1 = -6.076099099929818;
286+
constexpr Real c2 = 4.883198653547851;
287+
constexpr Real c3 = -11.88081601230542;
288+
constexpr Real c4 = 13.34658511480257;
289+
constexpr Real c5 = -8.722761043208607;
290+
constexpr Real c6 = 2.082038908808201;
291+
constexpr Real c7 = -7.389420998107497;
292+
constexpr Real c8 = -2.110913185058476;
293+
constexpr Real c9 = 0.2295491578006229;
294+
constexpr Real c10 = -0.9891538123307282;
295+
constexpr Real c11 = -0.08987150128406496;
296+
constexpr Real c12 = 0.3831132432071728;
297+
constexpr Real c13 = 1.054318231187074;
298+
constexpr Real c14 = 1.065556599652796;
299+
constexpr Real c15 = -0.7997496801694032;
300+
constexpr Real c16 = 0.3850133554097069;
301+
constexpr Real c17 = -2.078616693017569;
302+
constexpr Real c18 = 0.8756340772729538;
303+
constexpr Real c19 = -2.079022768390933;
304+
constexpr Real c20 = 1.596435439942262;
305+
constexpr Real c21 = 0.1338002171109174;
306+
constexpr Real c22 = 1.242891021876471;
307+
308+
// Note: a = 0.502500117621 / Sso
309+
constexpr Real a = 0.014289763856964;
310+
constexpr Real b = 0.057000649899720;
311+
312+
Real CtFreez;
313+
314+
Real Sar = Sa*1.0e-2;
315+
Real X = Kokkos::sqrt(Sar);
316+
Real Pr = P*1.0e-4;
317+
318+
CtFreez = c0
319+
+ Sar*(c1 + X*(c2 + X*(c3 + X*(c4 + X*(c5 + c6*X)))))
320+
+ Pr*(c7 + Pr*(c8 + c9*Pr)) + Sar*Pr*(c10 + Pr*(c12
321+
+ Pr*(c15 + c21*Sar)) + Sar*(c13 + c17*Pr + c19*Sar)
322+
+ X*(c11 + Pr*(c14 + c18*Pr) + Sar*(c16 + c20*Pr + c22*Sar)));
323+
324+
/* Adjust for the effects of dissolved air */
325+
CtFreez = CtFreez - SaturationFract*
326+
(1e-3)*(2.4 - a*Sa)*(1.0 + b*(1.0 - Sa/Sso));
327+
328+
return CtFreez;
329+
}
330+
238331
private:
239332
const int NVertLevels;
240333
};

‎components/omega/test/ocn/EosTest.cpp‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -315,6 +315,34 @@ int checkValueGswcSpecVol() {
315315
return Err;
316316
}
317317

318+
/// Test that the external GSW-C library returns the expected specific volume
319+
int checkValueCtFreezing() {
320+
int Err = 0;
321+
const Real RTol = 1e-10;
322+
Teos10Eos TestEos(1);
323+
// TestEos->EosChoice = EosType::Teos10Eos;
324+
constexpr Real SaturationFrac = 0.0;
325+
constexpr Real P = 500.0;
326+
constexpr Real Sa = 32.0;
327+
328+
///
329+
/// Get specific volume from GSW-C library
330+
double CtFreezGswc = gsw_ct_freezing_poly(Sa, P, SaturationFrac);
331+
double CtFreez = TestEos.calcCtFreezing(Sa, P, SaturationFrac);
332+
/// Check the value against the expected TEOS-10 value
333+
bool Check = isApprox(CtFreezGswc, CtFreez, RTol);
334+
if (!Check) {
335+
Err++;
336+
LOG_ERROR(
337+
"checkValueCtFreezing: CtFreez isApprox FAIL, expected {}, got {}",
338+
CtFreezGswc, CtFreez);
339+
}
340+
if (Err == 0) {
341+
LOG_INFO("checkValueCtFreezing: PASS");
342+
}
343+
return Err;
344+
}
345+
318346
// the main test (all in one to have the same log)
319347
// Single value test:
320348
// --> test calls the external GSW-C library
@@ -333,6 +361,7 @@ int eosTest(const std::string &MeshFile = "OmegaMesh.nc") {
333361

334362
LOG_INFO("Single value checks:");
335363
Err += checkValueGswcSpecVol();
364+
Err += checkValueCtFreezing();
336365

337366
LOG_INFO("Full array checks:");
338367
Err += testEosLinear();

0 commit comments

Comments
 (0)