Adds option to homogenize forces and fluxes fields

[breichl]

• Adds functions to do global averages on U and V grids in MOM_spatial_means
• Adds functionality to average all forcing and fluxes fields in MOM_forcing_types
• Adds flag to average all forcing and fluxes in MOM.F90
• This new feature is primarily for running in single column like configurations with the coupler, which requires perfectly equal forcing across all cells.

[breichl]

This code is meant to replace #1480 @Hallberg-NOAA let me know what you think of these changes. I still have testing and some functionality to complete, but I wanted to get your opinion on these changes compared to the original PR.

One thing I'm unsure of is the best way to average ustar in the fluxes routine. In forces it is fine, since we have taux and tauy. The options I can think of so far to get the correct ustar in fluxes are (1) to recompute the fluxes ustar from forces by passing forces to the fluxes homogenization routine, (2) to add the stress information (perhaps at h-points) to fluxes structure. I also need to either add functionality to read the gustiness or add the gustiness factor into the forces control structure.

[codecov]

Codecov Report

Merging #1541 (533dccc) into dev/gfdl (1324656) will decrease coverage by 0.06%.
The diff coverage is 1.81%.

❗ Current head 533dccc differs from pull request most recent head 7e3f66a. Consider uploading reports for the commit 7e3f66a to get more accurate results

@@             Coverage Diff              @@
##           dev/gfdl    #1541      +/-   ##
============================================
- Coverage     29.15%   29.09%   -0.07%     
============================================
  Files           239      239              
  Lines         71533    71698     +165     
============================================
+ Hits          20855    20858       +3     
- Misses        50678    50840     +162     

Impacted Files	Coverage Δ
src/core/MOM_forcing_type.F90	43.04% <0.00%> (-4.18%)	⬇️
src/diagnostics/MOM_spatial_means.F90	11.87% <0.00%> (-1.60%)	⬇️
src/core/MOM.F90	58.89% <37.50%> (-0.13%)	⬇️

---

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 1324656...7e3f66a. Read the comment docs.

[Hallberg-NOAA]

I think that this PR is a great improvement over its predecessor, and I like the general direction it is going. There are a few comments or issues that I see with this.

1. The new routines global_area_mean_u() and global_area_mean_v() will not reproduce across PE layouts when used in symmetric memory mode, and they will not reproduce across rotation by 180 degrees when used in non-symmetric memory mode. The basic problem is that they will double-count the velocity points at the edges of the PEs, or only count the ones on the one edge. In addition, it is using the integrated area of the unmasked tracer points, not the velocity points. For these routines, I think that the solution may be to do an appropriate interpolation of the velocities or stresses to the tracer points (perhaps with land masking and area-weighted averaging), and then do the global average on the tracer points.

2. What is wrong with taking the area-weighted average of ustar at the tracer points directly? Yes, it will be larger than the ustar that would be based on the area-averaged stresses, but in the same way that the sub-gridscale gustiness is something that we regularly include when calculating ustar in the first place, this difference is the legitimate result of the averaging operator not commuting with the nonlinear conversion of the velocities into ustar. As long as we are explicit about how the averaged ustar is calculated, why do we care if it can not be directly derived from the averaged stresses? If we do care about ustar being exactly the average of taux and tauy, that would probably be handled best in a separate subroutine to copy the information from a mech_forcing type to a forcing type, with the call controlled by another run-time flag to determine whether to calculate ustar that way.

3. There is a lot of repetition of nearly identical blocks of about 6 lines to take the averages in the new homogenize_forcing() and homogenize_mech_forcing() routines. (There are 53 such blocks for tracer points by my count.) The code might be a lot shorter (by perhaps 300 lines) and less error prone if a small local subroutine were introduced to perform this task of replacing a field with its average in a single-line call.

[breichl]

I caught a couple formatting issues here and wanted to raise two other issues here.

[breichl]

I noticed this is a repeat of code above. I think it is unintentional.

[breichl]

I don't know if this is worth addressing before taking this PR?

[breichl]

Fix formatting.

[Hallberg-NOAA]

In this expression, I think that we would prefer to have division at the end, or else put parentheses around the area divided by the masks, so that someone reading the code does not need to remember the order of precedence in Fortran between multiplication and division (they are the same), or whether the Fortran order of operations goes left to right or right to left. The same comment would apply in global_area_mean_u().

[breichl]

The most recent code push should address this.

[breichl]

I think that this PR is a great improvement over its predecessor, and I like the general direction it is going. There are a few comments or issues that I see with this.

1. The new routines `global_area_mean_u()` and `global_area_mean_v()` will not reproduce across PE layouts when used in  symmetric memory mode, and they will not reproduce across rotation by 180 degrees when used in non-symmetric memory mode.  The basic problem is that they will double-count the velocity points at the edges of the PEs, or only count the ones on the one edge.  In addition, it is using the integrated area of the unmasked tracer points, not the velocity points. For these routines, I think that the solution may be to do an appropriate interpolation of the velocities or stresses to the tracer points (perhaps with land masking and area-weighted averaging), and then do the global average on the tracer points.

2. What is wrong with taking the area-weighted average of ustar at the tracer points directly?  Yes, it will be larger than the ustar that would be based on the area-averaged stresses, but in the same way that the sub-gridscale gustiness is something that we regularly include when calculating ustar in the first place, this difference is the legitimate result of the averaging operator not commuting with the nonlinear conversion of the velocities into ustar.  As long as we are explicit about how the averaged ustar is calculated, why do we care if it can not be directly derived from the averaged stresses?  If we do care about ustar being exactly the average of taux and tauy, that would probably be handled best in a separate subroutine to copy the information from a `mech_forcing` type to a `forcing` type, with the call controlled by another run-time flag to determine whether to calculate ustar that way.

3. There is a lot of repetition of nearly identical blocks of about 6 lines to take the averages in the new `homogenize_forcing()` and `homogenize_mech_forcing()` routines.  (There are 53 such blocks for tracer points by my count.)  The code might be a lot shorter (by perhaps 300 lines) and less error prone if a small local subroutine were introduced to perform this task of replacing a field with its average in a single-line call.

Thanks for the review and suggestions, @Hallberg-NOAA.

I think I've fixed (1) based on our discussion offline. It now interpolates u/v point variables to tracer points and does the averaging there. I don't have a test for this yet.

I've added options for (2) to either compute ustar from the mean tau or as the mean of ustar. We cannot include gustiness at this point; it could be added but would require either a new get_param or passing it to the forces module via each cap. Since it is not clear how to do it and not clear that we need it I've put that decision off for now.

For (3) I've added three new subroutines to do the averaging now, which indeed reduces the file lengths.

[breichl]

@nikizadehgfdl can you try this at some point in your workflow? It is a little different than the previous attempt, but I believe will have the same functionality. The replacement parameter to activate homogenized forcing with this code will be "HOMOGENIZE_FORCINGS = True"