Fix the bug that Rayleigh_[uv] was not initialized in MOM_barotropic

[herrwang0]

This PR fixes the bug reported in Issue https://github.com/NOAA-GFDL/MOM6/issues/1326.

In the barotropic solver, Rayleigh_[uv] are not initialized, and their values are only assigned at points where CS%lin_drag_[uv] (which are the linear wave drag piston velocity at u/v point, averaged from the input file) are not zero.

Note that this only affects runs with BT_LINEAR_WAVE_DRAG = True

[codecov]

Codecov Report

Merging #1329 (ac88645) into dev/gfdl (d8733fe) will decrease coverage by 0.00%.
The diff coverage is 47.36%.

@@             Coverage Diff              @@
##           dev/gfdl    #1329      +/-   ##
============================================
- Coverage     45.75%   45.75%   -0.01%     
============================================
  Files           234      234              
  Lines         72542    72555      +13     
============================================
+ Hits          33189    33194       +5     
- Misses        39353    39361       +8     

Impacted Files	Coverage Δ
src/core/MOM_barotropic.F90	68.02% <47.36%> (-0.17%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update d8733fe...ac88645. Read the comment docs.

[marshallward]

Thank you for catching this and fixing it. We would not have caught this, as it is currently untested (https://codecov.io/gh/NOAA-GFDL/MOM6/src/dev%2Fgfdl/src/core/MOM_barotropic.F90). (Yet another reminder that our current coverage needs serious attention!)

I would suggest the following alternative, however:

if (CS%linear_wave_drag) then
    do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw
        Rayleigh_u(I,j) = 0.
    enddo ; enddo
    
    do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw
        Rayleigh_v(i,J) = 0.
    enddo ; enddo
endif

This will make two changes:

1. Reduce the stack when Rayleigh_[uv] is unused
   If Rayleigh_[uv] is untouched, then it will never be allocated (in normal circumstances, at least) so this initialization should be conditional.

2. Vectorize the initialization
   Embedding an if-block inside a do-loop can impede pipelining and vectorization speedups, especially when part of a more complex iteration. Overall it is better to structure if-blocks outside of loops.

[herrwang0]

Thank you for catching this and fixing it. We would not have caught this, as it is currently untested (https://codecov.io/gh/NOAA-GFDL/MOM6/src/dev%2Fgfdl/src/core/MOM_barotropic.F90). (Yet another reminder that our current coverage needs serious attention!)

I would suggest the following alternative, however:

if (CS%linear_wave_drag) then
    do j=CS%jsdw,CS%jedw ; do I=CS%isdw-1,CS%iedw
        Rayleigh_u(I,j) = 0.
    enddo ; enddo
    
    do J=CS%jsdw-1,CS%jedw ; do i=CS%isdw,CS%iedw
        Rayleigh_v(i,J) = 0.
    enddo ; enddo
endif

This will make two changes:

1. Reduce the stack when Rayleigh_[uv] is unused
   If Rayleigh_[uv] is untouched, then it will never be allocated (in normal circumstances, at least) so this initialization should be conditional.
2. Vectorize the initialization
   Embedding an if-block inside a do-loop can impede pipelining and vectorization speedups, especially when part of a more complex iteration. Overall it is better to structure if-blocks outside of loops.

Hi Marshall,
Thanks for the comment! The alternative change makes more sense.

I have a few questions related:

1. Rayleigh_[uv] is not an allocatable variable, does a conditional initialization still have an impact on stack?
2. IF we initialize Rayleigh_[uv] regardless, is it better to remove the IF-block here (inside a loop) in https://github.com/NOAA-GFDL/MOM6/blob/45b336670988773c050e3d59e3bd0ebe7b948ecb/src/core/MOM_barotropic.F90#L1980-L1985 and keep Rayleigh_[uv] in the equation? Since it would be zero if CS%linear_wave_drag is False? I guess in that case we are trading an IF-clock with an additional operation. Not sure which is preferred.

It is probably not a big deal in any case. These are more of questions out of curiosity.

[marshallward]

Very reasonable questions! I was trying not to bore with the details.

1. It's not an explicit allocation, you are right, but there is usually an implicit allocation of memory for stack variables (either via malloc or stack pointer adjustment). This is not a hard rule, and depends on the compiler, libc's malloc, and the OS, but it is usually a very good assumption.

   Another way to think of it is that we are giving the compiler and OS memory manager the freedom to do what's best, which in most cases is to do nothing.

2. You are right, that if-block is also problematic, and leaving Rayleigh_[uv] as zero might end up being faster than using the if-block. Even better would be to split the loop and set up two separate equations, one with and without the drag.

   But this would also be a bit more work, and might be better handled in a broader cleanup of the whole function.

I don't expect any of this to make a huge difference (one loop in one function is never going to matter too much) so it's up to you how far to go with it. There will be other opportunities to resolve it.

[herrwang0]

I decided to make the changes we discussed, so that hopefully the code related to the Rayleigh drag term is in an optimized state.

Specifically, the initialization is now conditional and split from other variable.
Adding Rayleigh terms is also in a separate loop

One thing I am not sure about is the about openmp nowait. It looks like it is not used in the first variable updates for [uv]bt, [uv]bt_trans and [uv]_accel_bt (v in odd number steps and u in even number steps). I'm not sure if this is intentional or not, so just followed what was there.

In addition, there seems to be bug related to resetting [uv]bt to zero with CS%vel_underflow. Since this feels like a different issue and I am not familiar with this configuration, I put this in a separate issue https://github.com/NOAA-GFDL/MOM6/issues/1336

[marshallward]

Thank you, this looks good. I'll merge this to dev/gfdl (pending testing) once the main candidate has been merged.

[marshallward]

https://gitlab.gfdl.noaa.gov/ogrp/MOM6/-/pipelines/12118 ✔️