Velocity remap dimensional consistency

[iangrooms]

The dimensional consistency tests were failing with the old code because of the presence of a literal constant on these lines. I replaced those lines by logic to avoid divide-by-zero, and the new code passes the dimensional consistency test. I also cleaned up the unit scaling to make it more easily legible.

[iangrooms]

PS this will change answers whenever the remap correction is turned on. Dimensional tests are failing, but it's because of #302. If you set FULL_DEPTH_KHTR_MIN=False the dimensional consistency tests pass.

[alperaltuntas]

This is passing our pr_mom suite, but I think we should discuss the CI regression test failures before I merge. @gustavo-marques, @mnlevy1981 .

[mnlevy1981]

My calendar is up to date, and I'm happy to chat any time I'm free today... but in general, I am in favor of answer-changes when we move from "add epsilon to a non-negative denominator to avoid dividing by zero" to "do division when denominator is non-zero / evaluate limit when denominator is zero"

[mnlevy1981]

@alperaltuntas reminded me that we need to use the answer-date construct in MOM6 to preserve answers for GFDL; he's going to update this PR with an option to maintain the previous answers.

I'm still a little confused about when rescale_coef was first introduced, though. It seems like the answer-preserving formulation will fail dimensional scaling tests, so it's hard to believe that computation was ever used in a production run at GFDL. Are we sure this answer-changing commit needs a mechanism to re-create the old answers as well?

[iangrooms]

The original PR #277 was merged May 24, 2024. I don't think this was ever used at GFDL, but I don't know that for sure.

[alperaltuntas]

I am noting here that reverting the changes in lines 383 and 388 fixes the regression errors, which indicate yet another dimensional consistency issue.

[mnlevy1981]

I am noting here that reverting the changes in lines 383 and 388 fixes the regression errors, which indicate yet another dimensional consistency issue.

If the units are really W m-2, then I think the correct unit scaling is US%RZ3_T3_to_W_m2 (rather than the original US%RZ3_T3_to_W_m2 * US%L_to_Z**2 or the new GV%H_to_kg_m2 * US%L_T_to_m_s**2 * US%s_to_T).

[mnlevy1981]

I am noting here that reverting the changes in lines 383 and 388 fixes the regression errors, which indicate yet another dimensional consistency issue.

If the units are really W m-2, then I think the correct unit scaling is US%RZ3_T3_to_W_m2 (rather than the original US%RZ3_T3_to_W_m2 * US%L_to_Z**2 or the new GV%H_to_kg_m2 * US%L_T_to_m_s**2 * US%s_to_T).

@klindsay28 points out that this comment is an over-simplification; my thinking was that all of these computations should be done in the model dimensions (RZ3_T3) and converted to W_m2 in the post_diag call, but it ignores the fact that these intermediate computations might have some terms already in physical units or that need to be converted to density space or a whole host of other issues that could be causing problems

[iangrooms]

The quantity is velocity squared times vertical thickness divided by dt, so the dimensions are H L^2 /T^3. For output it needs to also be multiplied by density. US%RZ3_T3_to_W_m2 is not appropriate; it would assume that the dimensions of velocity are RZ2_T2 where they are in fact L2_T2. After discussing with Bob Hallberg I decided that it makes the most sense to put all the dimensional/unit scaling in one place, namely lines 383 and 388. Previously the code had dimensional constants scattered in more than one place. Also, after talking with Bob I realized that the variable GV%H_to_kg_m2 is a better way to handle the Boussinesq and non-Boussinesq cases simultaneously.

[mnlevy1981]

The quantity is velocity squared times vertical thickness divided by dt, so the dimensions are H L^2 /T^3. For output it needs to also be multiplied by density. US%RZ3_T3_to_W_m2 is not appropriate; it would assume that the dimensions of velocity are RZ2_T2 where they are in fact L2_T2. After discussing with Bob Hallberg I decided that it makes the most sense to put all the dimensional/unit scaling in one place, namely lines 383 and 388. Previously the code had dimensional constants scattered in more than one place. Also, after talking with Bob I realized that the variable GV%H_to_kg_m2 is a better way to handle the Boussinesq and non-Boussinesq cases simultaneously.

I think I understand. With your new scaling factor of GV%H_to_kg_m2 * US%L_T_to_m_s**2 * US%s_to_T, do you need to drop the GV%H_to_RZ terms when computing du2h_tot and dv2h_tot? It seems like you're double-counting the H term as-is

[iangrooms]

Yes, that should be dropped. I thought I did it, but evidently forgot to remove it for this PR. I will make that change asap

[mnlevy1981]

As discussed in the developer's meeting, I think the current test failures (round-off level changes in ocean_model-ale_u2 and ocean_model-ale_v2) are expected due to refactoring of the scaling coefficients. They won't effect GFDL bit-for-bit requirements since these variables are not currently available on dev/gfdl

[alperaltuntas]

pr_mom test is ongoing and should be done in 30 mins.