Which layer thicknesses to use for thickness-weighted momentum diagnostics and vorticity budgets

[hmkhatri]

@StephenGriffies and I have been using depth-integrated 2D momentum diagnostics for vorticity budget analysis in global MOM6 simulations and we observed that some math expressions in vorticity budgets do not agree when compared against these diagnostics. This is probably due to the numerical treatment of individual terms in the model.

An important issue is - which layer thicknesses to use for thickness-weighted momentum diagnostics such that the discrepancy between the online diagnostics and theoretical expressions is minimal. This would make it much easier to analyze various budget terms from the model output and compare them against theories.

I added 2D momentum diagnostics last year and, in the current implementation of 2D momentum diagnostics, the layer thicknesses at u/v points (diag_hu / diag_hv) are returned in the btstep call as part of the ADp control structure.

https://github.com/NOAA-GFDL/MOM6/blob/2912058b17b5e8636c8687ed236ac68aba36f609/src/core/MOM_dynamics_split_RK2.F90#L761-L764

https://github.com/NOAA-GFDL/MOM6/blob/2912058b17b5e8636c8687ed236ac68aba36f609/src/core/MOM_barotropic.F90#L2678-L2687

With this choice of h, some terms look unrealistically large. For example, consider the following equation

Using the new 2D momentum diagnotics, one can compute the LHS of the equation above as
d/dx [ intz_CAv_2d - intz_rvxu_2d - intz_gKEv_2d ] - d/dy [intz_CAu_2d - intz_rvxv_2d - intz_gKEu_2d ] and this should roughly be equal to - beta x V (others terms are very small), which one can compute independently using vmo_2d / vmo diagnostic.

As seen in the figure, fields from the two calculations differ by more than an order of magnitude. Thus, it is critical to use thicknesses consistent with numerics so that thickness-weighted diagnostics result in correct magnitudes of terms in higher-order budgets.

I was thinking maybe we could use h_av, which is used in the CorAdcalc call. Since MOM6 uses [ f / h_av] x (hu) for the Coriolis term, multiplying h_av to CAu/CAv should result in (f + zeta) x thickness flux. If we sum f x thickness flux over all layers, then the curl of the resultant field should be close to - beta x V.

https://github.com/NOAA-GFDL/MOM6/blob/2912058b17b5e8636c8687ed236ac68aba36f609/src/core/MOM_dynamics_split_RK2.F90#L719-L720

@Hallberg-NOAA @adcroft Would appreciate your thoughts on this issue and which layer thicknesses to use for these diagnostics. If you agree that h_av is a good place to start with, could you please guide me on where I should add code for calculating thicknesses hu / hv for 2D diagnostics? I see that h_av is computed multiple times in MOM_dynamics_split_RK2.F90 around predictor and corrector steps. I am not sure which h_av to use.

@NoraLoose @gustavo-marques are looking into thickness weighted momentum diagnostics in the CPT project. I hope discussing this issue will be helpful for their analysis too.