Removing potential divide by zero error in Orlanski package (attempt #2)

[garrettdreyfus]

This PR is a bug fix for issue #822 . It removes the possibility of a divide by zero error in the obcs/orlanski_[north,east,south,west].F files.

There is an ongoing discussion about whether when the denominator is zero, CL should be set to 0 or CMAX.

This pull request has not yet been tested in a set up with the orlanski condition enabled in each direction.

[mjlosch]

This looks good to me, thanks for providing the bug fix. The changes do not affect the verification experiments, but setting CL=CMAX when denom=0 will.

I have a few cosmetic suggestions (header, indentation, etc) that the not really related to this PR, but will push them once we are happy with the actual code changes.

[mjlosch]

I agree with @garrettdreyfus that for continuity reasons, CL=CMAX for denom==0. makes more sense as the denominator can asymptotically approach zero for slowly varying flow; maybe with an appropriate mask to avoid nonzero values on land points. Comments are welcome.

[mjlosch]

@samarkhatiwala, do you have any opinion on the CL=CMAX vs O question?

[mjlosch]

We decided to leave as it is now (including setting CL=0 for denom=0), to keep results unchanged and pick up the issue of CL to be 0 or CMAX in a separate issue or PR where we have more time to test this. @garrettdreyfus are you OK with this? Also would you be willing to provide some (future?) experience with using 0 or CMAX for this (in a new issue)?

[garrettdreyfus]

Hey @mjlosch ,

Sounds great! I can run some experiments setting CL to 0 or CMAX when denom is 0. Any suggestions on experimental setups that would help tease out differences between the two choices? My suspicion is that denom equaling 0 is a rare enough case that the choice won't make much of a difference regardless of experimental setup.

[jm-c]

@garrettdreyfus and @mjlosch I am going to do a quick check before/towards merging this PR soon.

[mjlosch]

Any suggestions on experimental setups that would help tease out differences between the two choices? My suspicion is that denom equaling 0 is a rare enough case that the choice won't make much of a difference regardless of experimental setup.

Unfortunately, no, but you could try to run your setup where you discovered the bug, and show the differences between CL=0 and CMAX. Maybe there are hardly any, but maybe there's a bit of indication, if the choices affect the solution in a bad way. We test the western and eastern orlanski code in verification/dome and verification/tutorial_plume_on_slope.
Change CL=0 to CL=CMAX changes the results from (output of ./testreport -devel -j 6 -t 'dome tutorial_plume_on_slope' run from verification):

Y Y Y Y>14<16 16 16 16 16 16 16 16 16 16 16 16 16 13 14 16  .  .  .  .  .  .  .  . pass  dome
Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 22 22 22 22  .  .  .  .  .  .  .  . pass  tutorial_plume_on_slope
Y Y Y Y>16<16 16 16 16 16 16 16 16 16 16 16 16 22 22 22 22  .  .  .  .  .  .  .  . pass  tutorial_plume_on_slope.roughBot

to

Y Y Y Y> 4<16 16  8  8 16 11 16 11  7  9  4  5 16  2  3  5  .  .  .  .  .  .  .  . FAIL  dome
Y Y Y Y> 4< 7  8  8  6 16 16 16 16  6  5  1  4 22 22 22 22  .  .  .  .  .  .  .  . FAIL  tutorial_plume_on_slope
Y Y Y Y> 4< 7  8  8  6 16 16 16 16  6  6  1  4 22 22 22 22  .  .  .  .  .  .  .  . FAIL  tutorial_plume_on_slope.roughBot

and it's unclear why this is happening.

[jm-c]

@garrettdreyfus I made minor adjustments (some related to your changes and some minor ones to improve a little bit how this code reads). As expected, the results are unaffected (given it's only minor changes). If you are OK with those we could go and merge this PR soon.

[garrettdreyfus]

@mjlosch Thanks for the lead! I'll check that out when I dig into it.

@jm-c Sounds great thanks for all the help!

[jm-c]

@garrettdreyfus I will merge your PR very soon (later today) and thanks for your contribution.