Omega/add eos clamp

[alicebarthel]

This PR constrains the P, Sa, T ocean state fields to the valid range as defined by the ''oceanographic funnel'' of the GSW function gsw_infunnel (originally based on McDougall et al. 2003).
The input values are pre-processed against the valid P,Sa,T bounds and the specific volume calculation is performed on the modified values (if necessary). The original state fields are not modified.
The functions to calculate the valid bounds include the polynomial approximation of the freezing temperature. It is faster but less accurate than the full GSW/Teos-10 freezing temperature calculation.

I included a unit test to check the freezing temperature against the equivalent function in GSW-C.
The code compiled and passed tests on pm-cpu (gnu) and pm-gpu.

Checklist

• Documentation:
  • User's Guide has been updated
  • Developer's Guide has been updated
  • Documentation has been built locally and changes look as expected
• Building
  • CMake build does not produce any new warnings from changes in this PR
• Testing
  • A comment in the PR documents testing used to verify the changes including any tests that are added/modified/impacted.
  • Unit tests have passed. Please provide a relevant CDash build entry for verification.

[alicebarthel]

It seems like a forgot to run the pre-commit checks. I'll update.

[philipwjones]

This looks fine and also passes Ctests on Frontier cpu/gpu. I had one minor change but will go ahead and approve this.

[alicebarthel]

Adding a comment for discussion here. @katsmith133 raised the good point of whether the T,S bounds should apply to the clamped values or the original values. I wrote everything assuming that we calculate, say the freezing temperature bounds, on the actual (P,S) fields.
The conundrum is: how strongly do we want to constrain our values to the original "funnel"? Its original definition limits it to P < 8,000 dbar, for example, though the Poly75t was evaluated down to 10,000dbar:

For context, the 75-t polynomial was tested against the full TEOS-10 in the funnel (which was used for the fit), but also in WOA and the "cube" range (SA=0-42; CT= -2-40 C; P=0-10,000 dbar) (Table3 of Roquet et al 2015)

[alicebarthel]

As far as I understand, the TEOS-10 toolbox documents the range of validity (and offers a check) but does not restrict the range to the funnel.
e.g.

def specvol(SA, CT, p):
    """
    Calculates specific volume from Absolute Salinity, Conservative
    Temperature and pressure, using the computationally-efficient 75-term
    polynomial expression for specific volume (Roquet et al., 2015).

    Parameters
    ----------
    SA : array-like
        Absolute Salinity, g/kg
    CT : array-like
        Conservative Temperature (ITS-90), degrees C
    p : array-like
        Sea pressure (absolute pressure minus 10.1325 dbar), dbar

    Returns
    -------
    specvol : array-like, m^3/kg
        specific volume


    Notes
    -----
    Note that the 75-term equation has been fitted in a restricted range of
    parameter space, and is most accurate inside the "oceanographic funnel"
    described in McDougall et al. (2003).  The GSW library function
    "gsw_infunnel(SA,CT,p)" is available to be used if one wants to test if
    some of one's data lies outside this "funnel".

[vanroekel]

@alicebarthel I'm not sure I'm tracking what the discussion point is. bounds and clamping are confusing me. Let me see if I understand and please correct me where I'm off. Right now you compute bounds based on the state, but the question is if we should compute bounds on the valid range? If this is right it seems cyclical to me. Is this a question of allowing say P > 8000 in the calculation, which would loosen the T/S range some?

If I'm tracking, how different do we expect these options to be? The ranges you show are already quite broad it seems.

[alicebarthel]

Thanks @vanroekel. It's a detail of the implementation.

The easiest example might be:
if S = 44 (above the 42 g/kg funnel bound), should the T bounds be based on gsw_ct_freezing(sa, p, 0) or gsw_ct_freezing(42, p, 0)?

In reality, I don't expect it to apply to many (or any?) locations at all, but perhaps some edge cases. We can have locations with S > 42 g/kg (e.g. the Red Sea and sometimes the Mediterranean), but they tend to be warm.
As I said, I see no indication in the gsw toolbox that they limited the T,S,P range when using the eos, except for documenting the validity range.

[vanroekel]

Thanks @alicebarthel I feel like I wasn't too far off in my understanding. My follow up question would be how different would those two results be in your example? Just trying to see how important of a consideration this distinction is.

[alicebarthel]

@vanroekel For my example above, the difference in density is about O(0.01).

Similar order of magnitude than calculating rho based on the S bound (42) rather than S itself:

[sbrus89]

This looks good to me, @alicebarthel and @katsmith133. I just have one suggestion.

[katsmith133]

Approving based upon visual inspection and my own testing. Thanks, @alicebarthel!

[vanroekel]

approving based on visual inspection. Discussion of clamp implication will be moved to a separate discussion.

Thanks @alicebarthel

[alicebarthel]

@vanroekel I added a short description of the underlying functions for teos-10. It's not meant to be comprehensive but a pointer for developers to understand the building blocks and their relative cost (hint: reuse the coefficients if you can). I hope that this additional info is helpful.

[alicebarthel]

The code compiled and passed tests on pm-cpu (gnu) and pm-gpu.

[philipwjones]

Unfortunately, this will not build on Frontier because the LOG_WARN macros for clamping are a host function inside a device function. And I don't think we want I/O or Logging from inside these loops anyway for performance reasons. The easiest fix might be to have a separate function (with the loops inside) that performs clamping in a reduction loop that counts the number of times things get clamped and has the log warning outside the loop. This would be called before the loop where the calcPCoeffs, calcDelta and calcRefProfile functions are called. There are other possible ways around this if that option in unpalatable.

While I was looking at LOG statements, I also found a couple of other minor log-related changes.

[philipwjones]

Change to ABORT_ERROR

[alicebarthel]

Thanks. Out of curiosity, is the use case for LOG_ERROR when it's an error but we don't want to exit? Can you think of a use case for it?
I understand that ABORT_ERROR is the most appropriate to exit when encountering a critical error?

[philipwjones]

Yes, once we created the Error handler, some of the LOG macros are superseded (though the Error handler uses the LOG macros underneath). You can still use LOG_ERROR to write a message with an error prefix to the log and not abort/return. A possible use case might be when you have a lot of information to write and multiple messages might be cleaner than a single large message with multiple lines. And yes, ABORT_ERROR is preferred for critical errors. I've been going through and replacing LOG macros with ERROR macros when I can but haven't made it up to the ocn directory yet and just happened to see these.

[philipwjones]

While this might be useful info, this will get printed every time the routine is called and wil generate too much output.

[alicebarthel]

Makes sense, this was mostly useful in the dev stage, probably unhelpful for production.

[alicebarthel]

Thanks @philipwjones, that's very helpful to learn more about the LOG_WARN macro. If I understand correctly, the proposed structure would be a pre-processing step where we i) check the full arrays for out-of-bounds value, sum the number of locations, and write out the warning, ii) with clamping enabled, this step makes a copy of the state arrays, i.e. PClamped, TClamped, SClamped, which we then run the eos on.

[philipwjones]

Yes, that's what I was proposing. Another approach would be to have all of those routines pass a count variable and make sure the calling loop is a reduction loop. Then do a check of that count outside the loop for the warning. And a third approach might be to have an internal flag that gets tripped if any variable is clamped - it would need to be an atomic update since it's within the loop (Kokkos docs have a section on atomic operations).