Better error messages in interpolate_for_nondim_position

[Hallberg-NOAA]

Cleaned up the error messages when there are fatal problems in interpolate_for_nondim_position() in the MOM_neutral_diffusion module to explicitly give an indications of all the problems. These had previously triggered fatal errors (or should have), but had less explicit or incorrect error messages in some cases, so all solutions and behavior are identical in any cases that worked previously.

[marshallward]

I believe that the conditionals can be simplified due to the boundedness of interpolate_for_nondim_position, and the error message can be consolidated into a single block.

While this does not need to be addressed now, we should explore if this function might be written without branching, and without MOM_error calls. Both coule be making vectorization and GPU migration more difficult.

[marshallward]

This could be written with any(...).

However, I think the interpolate_for_nondim_position checks for > 1. and < 0. are too conservative. I don't see any way that these conditions could ever be true.

• Ppos < Pneg is an error, Ppos == Pneg is 0.5, Ppos > Pneg is primarily a flag to inspect dRhoPos and dRhoNeg.
• dRhoPos < dRhoNeg is an error, dRhoPos > dRhoNeg is a mathematical function that must be between 0 and 1 (inclusve). dRhoPos == dRhoNeg is explicitly either 0, 1, or 0.5 (the final case being double-expressed for Ppos > Pneg and Ppos == Pneg.)

So I don't believe this needs to be checked.

An alternative form is something like

if (Ppos - Pneg > 0. .and. dRhoPos - dRhoNeg >= 0.) then
  if (dRhoPos - dRhoNeg > 0.) then
    r = min(1., max(0., ...)
  else ! dRhoPos == dRhoNeg
    if (dRhoNeg > 0.) then
      r = 0.
    elseif (dRhoNeg < 0.) then
      r = 1.
    else ! dRhoPos == dRhoNeg == 0
      r = 0.5
    endif
  endif
elseif (Ppos - Pneg == 0.) then
  r = 0.5
else ! Ppos < Pneg .or. dRpos < dRhoNeg
  call MOM_error()
endif

BTW I am not sure if Ppos - Pneg == 0 (for example) was deliberately chosen rather than Ppos == Pneg. (I find the latter more readable, but the results could be different for Pneg / Ppos < 10**(-16)). (Nevermind, I don't think this is an issue.) If the exisiting forms were deliberate, maybe this should be noted in a comment. Either way, it should not change the bounds of interpolate_for_nondim_position.

[Hallberg-NOAA]

These are all good questions. Probably the right person to answer them would be @ashao, who wrote the original version of this code.

I agree with your analysis of how the logic could be simplified. Everything in this commit is retaining the original logic, but just making the messages more explicit. I would be happy to modify this as suggested, but before doing so, I would prefer to hear what @ashao thought about this.

[Hallberg-NOAA]

This PR has been revised as suggested.

[marshallward]

Gaea regression: https://gitlab.gfdl.noaa.gov/ogrp/mom6ci/MOM6/-/pipelines/24069 ✔️