Adding the gradient model into MOM6

[sinakhani]

This is a new PR that addresses comments by @marshallward and @Hallberg-NOAA to the previous version. Thank you for the comments.

[Hallberg-NOAA]

This PR must be revised so that it will not break existing configurations before it could be accepted.

Please be considerate of the Hippocratic principle of community ocean modeling, "Do no harm, and do not change other people's answers without their prior consent." In my view this is the single most important consideration for the long-term viability and cohesiveness of the MOM6 development community.

[sinakhani]

Thank you for the comments, @Hallberg-NOAA . I have prepared a revision based on your comments. Please let me know what you think. I also include @marshallward here. Thank you for your time and consideration.

[Hallberg-NOAA]

I suspect that there are conversion arguments missing for this and the following 3 diagnostics, perhaps conversion=US%s_to_T or conversion=GV%H_to_mks*US%m_to_L*US%s_to_T) for the first two and conversion=US%L_to_m**2 for the latter two.

You can verify all of this (and the dimensional consistency of these new code options in general) by repeating a short run with these new capabilities with various integers set for the 6 debugging runtime parameters [ZLTRQCSH]_RESCALE_POWER and verifying that your new diagnostics and solutions are identical in each case. (See your MOM_parameter_doc.debugging file for a description of these runtime parameters.)

[Hallberg-NOAA]

This comment still applies and has not been addressed by the latest revised version of this PR. Please add the appropriate conversion factors and test that they are correct. Also the one-point indent in the continuation lines of these register_diag_field() calls is not in compliance with the 2/4 point indexing convention articulated in the MOM6 style guide.

[Hallberg-NOAA]

Please remove this added white space and the other white space changes in calc_slope_functions_using_just_e(), at lines 1009 of the original code and lines 1085, 1109, 1113, 1132 and 1137 of the revised code. This PR is no longer making any substantive changes in calc_slope_functions_using_just_e(), so there is no need for these white-space changes either.

(This is very minor, and I would not suggest this were there not already the need for other changes to fix likely problems with this PR.)

[sinakhani]

The suggested changes have been made.

[Hallberg-NOAA]

Not all of the suggested reversions to the original white spacing have been made. Given that this PR is already very extensive, please do not mix in more changes than is absolutely necessary. Please see github.com//pull/937/changes and revert the highlighted changes to calc_slope_functions_using_just_e() that are not related to the substantive changes that are central to this PR.

[Hallberg-NOAA]

This expression will not pass rotational symmetry testing. The first part of the numerator in this expression is staggered at the same point as h(i+1,j), while (G%IareaT(I+1,j) + G%IareaT(I,j)) and the normalized thickness gradient is staggered at the same point as u(I,j), and G%dyT(i,j) is staggered at the same point as h(i,j).

You can verify the rotational consistency of these new code by doing a small single-PE run that exercises it with the run-time settings ROTATE_INDEX = True and INDEX_TURNS set to 0, 1, 2, or 3. If the expressions are rotationally consistent, you will get identical answers in each case. If the answers change (as I suspect they will here) you will get different answers.

[sinakhani]

The gradient slopes Ux_Hx(I,j) and Vx_Hx(I,j) are defined on u-points, and Uy_Hy(i,J) and Vy_Hy(i,J) on v-points. For this purpose, please consider the following definition for Ux_Hx = du/dx * dh/dx:

  Ux_Hx(I,j) = ((G%IdxCu(I+1,j)*G%IdyCu(I+1,j)*uh(I+1,j,K)) - (G%IdxCu(I,j)*G%IdyCu(I,j)*uh(I,j,k))) * (G%IareaT(I+1,j) + G%IareaT(I,j)) * G%dyT(I,j) * ((h(i+1,j,k) - h(i,j,k))/( h(i+1,j,k) + h(i,j,k) + h_neglect))

In which du/dx = ((G%IdxCu(I+1,j)*G%IdyCu(I+1,j)*uh(I+1,j,K)) - (G%IdxCu(I,j)*G%IdyCu(I,j)uh(I,j,k))) / (1/2)(h(i+1,j,k) + h(i,j,k) + h_neglect))

and dh/dx = (1/2)*(G%IareaT(I+1,j) + G%IareaT(I,j)) * G%dyT(I,j) * ((h(i+1,j,k) - h(i,j,k))

where coefficients (1/2) are cancelled from numerator and denominator. I think the Bold and Italic texts {(h(i+1,j,k) + h(i,j,k) + h_neglect)) and ((h(i+1,j,k) - h(i,j,k))} shows that h(i,j) is defined in u-points in agreement with your staggering points.

Similar definitions are hold for other slops for gradient terms Vx_Hx, Uy_Hy and Vy_Hy.

[Hallberg-NOAA]

From my reading of the code, I still think that there are significant inconsistencies in where the various terms in this expression are discretized, exactly as described in the original comment. This can be verified empirically by rotational consistency testing, which I am confident that this PR as currently written will not show identical answers under grid rotation. @sinakhani, let's please take this offline and set up a virtual conversation to work through these issues.

[sinakhani]

Thank you for the comments, @Hallberg-NOAA. The results for a rotation test are shown under the discussion below (please see them here). I will reach out via email.

[Hallberg-NOAA]

As the code is currently written, this expression is dimensionally inconsistent, especially in non-Boussinesq mode. The left hand side has units of [L2 T-1 ~> m2 s-1], whereas the right hand side has units of [L2 ~> m2] * [H L-1 T-1 ~> s-1 or kg m-3 s-1] which combines into [H L T-1 ~> m2 s-1 or kg m-1 s-1], although I suspect that this will be sorted out when the dimensional consistency issues inside of calc_slope_functions_with_gradient_model() are sorted out.

[sinakhani]

So, the dimension for h is [H ~> m or kg m-2], and the dimension for e is [Z ~> m]. Also, the equation for layer thickness based on interface position is h_i = e_{i-1} - e_{i}. In this case, this equation is dimensionally inconsistent. How this inconsistency is sorted in non-Boussinesq mode?

[Hallberg-NOAA]

In both Boussinesq and non-Boussinesq mode we can convert an array of layer thicknesses (in [H ~> m or kg m-2]) into vertical distances across layers (in [Z ~> m]) by calling thickness_to_dz(). Internally, this routine either does a simple multiplication by the scaling factor GV%H_to_Z in Boussinesq mode or in non-Boussinesq mode it multiplies the thickness by the pre-calculated layer-averaged specific volume times the scaling factor GV%H_to_RZ, along with some tests to ensure that this pre-calculated specific volume is valid for the given halo size.

If you want to find the geometric positions of the interfaces (in [Z ~> m]) as opposed to the layer extents, you can call find_eta(), which also takes the layer thicknesses and a thermo_vars_ptrs type as an argument. Both thickness_to_dz() and find_eta() are in MOM6/src/core/MOM_interface_heights.F90.

[sinakhani]

I do not convert from layer thickness to vertical distance because the gradient model uses the layer thickness. I have updated units of variables to make sure that there no inconsistency on dimensions (I think dimensions are now specified correctly in calc_slope_functions_with_gradient_model). The eddy diffusivity has units of [kg m-1 s-1] and/or [m2 s-1] for non-Boussinesq and Boussinesq modes. Please let me know if this is fine with you.

Also, I have attached the results of testing for rescale-power and rotation. Please see the attached plot.

Please let me know if I need to do any further revisions before we can merge these updates into the main branch. Thank you.

[Hallberg-NOAA]

This revised version is getting much closer, but I think that there are some dimensional inconsistencies and horizontal staggering bugs (as would be manifest in a breaking of rotational symmetry) that should be sorted out. Please verify that you get the same answers and diagnostic output when you use different values of L_RESCALE_POWER and H_RESCALE_POWER, and also when you set ROTATE_INDEX = True and use different values of INDEX_TURNS.

[sinakhani]

Thank you for the second round of the review. I have prepared a revision based on the comments.

Also, I have tested the [LHZTRQCS]_RESCALE_POWER being 0, 1, 2 and ROTATE_INDEX = True with INDEX_TURNS = 0, 1, 2, 3. Please see the attached PDF file for a comparison on u field.
Gradient_Model_Rescale_Rotation_Tests.pdf.

Moreover, I have a question for the comments on the dimension of KH_u and KH_v terms (please see the comment above). I think the dimension inconsistency in these equations are similar to that in the definition of thickness layer from interface position (i.e. h_{i} = e_{i-1} - e_{i} ). h has a dimension of [H ~> m or kg m-2] while the dimension of e is [Z ~> m]. I think I can use similar arguments to resolve the dimension inconsistency in KH_u and KH_v equations to that was used to sort out the inconsistency in h definition based on e. Am I correct? Any comments and thoughts are very welcome.

Thank you for your consideration. @Hallberg-NOAA @marshallward