KE-conserving correction to velocity remap

[iangrooms]

This PR does two things:

• Adds an option to correct the remapped velocity so that the KE is conserved. This can be used in combination with any remapping scheme. It is not a new remapping scheme, it is a correction to the remapped velocity.
• Adds two diagnostics, ale_u2 and ale_v2, that measure the rate of depth-integrated KE loss (or gain) associated with remapping velocity. The rates are computed separately for the u and v components, and the output units are Watts per square meter.

The correction is described in more detail here ke_conserving_remap.pdf

[Hallberg-NOAA]

I agree that du2h_tot and dv2h_tot should have units of [R Z L2 T-3 ~> W m-2], but as currently written they actually have units of R H L2 T-3 ~> W m-2 or kg2 m-5 s-3], as is reflected by the consistent conversion factors in the register_diag_field calls for these diagnostics on lines 365 and 371 above.

This is easily corrected, though, if instead of using GV%Rho0 where du2h_tot and dv2h_tot are calculated at lines 1167 and 1235, we were to use GV%H_to_RZ. In Boussinesq mode, GV%H_to_RZ = GV%Rho0*GV%H_to_m*US%m_to_Z, but in non-Boussinesq mode the thicknesses already include the masses and GV%H_to_RZ is just the (constant) rescaling factor to convert them to units of [kg m-2].

This suggested change has the added benefit of not introducing the Boussinesq reference density (GV%Rho0) into non-Boussinesq calculations. In non-Boussinesq mode, I have been working to avoid the use of GV%H_to_m or GV%H_to_Z whereever possible because they imply division by the arbitrary Boussinesq reference density. In this case the various factors would cancel out mathematically, but multiplying and then diving by 1035 kg m-3 changes answers at roundoff, so they can not be cancelled out after people start using the code without changing answers.

[Hallberg-NOAA]

Please revise this expression from du2h_tot(I,j) = GV%Rho0 * u2h_tot / dt to du2h_tot(I,j) = GV%H_to_RZ * u2h_tot * I_dt with the line if (present(dt)) I_dt = 1.0 / dt added somewhere outside the i- and j- loops for this routine.

Changing out GV%Rho0 for this rescaling factor will place du2h_tot into the documented units of [R Z L2 T-3 ~> W m-2], and it avoids using the Boussinesq reference density in non-Boussinesq cases.

An optimizing compiler will replace the divisions by dt with a multiplication by its reciprocal but without optimizations this won't happen. If we explicitly do this replacement ourselves, we help to ensure that the optimized and non-optimized (debugging) versions of the code are closer to giving the same answers.

[Hallberg-NOAA]

Please revise this expression from dv2h_tot(i,J) = GV%Rho0 * v2h_tot / dt to dv2h_tot(i,J) = GV%H_to_RZ * v2h_tot * I_dt, for the reasons discussed above with the corresponding suggested changes to du2h_tot.

[Hallberg-NOAA]

If the changes suggested below are made, these conversion factors would become conversion=US%R_to_kg_m3*US%Z_to_m*US%L_T_to_m_s**2*US%s_to_T or equivalently conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2. The same comment also applies to the register_diag_field() call for CS%id_remap_delta_integ_v2.

Also, it would be more convenient for verifying the consistency between the conversion and units arguments could be kept entirely on the same line and not split the conversion factor with line continuation. In this case, this could easily be accommodated within the 100 character recommended line length by moving 'applied .` up to the preceding line. Also note that the spaces before and after the
continuation between lines 369 and 370 will lead to a double space that we might not want.

[iangrooms]

Thanks @Hallberg-NOAA for your review! I am still working on mastering the unit scalings and Boussinesq vs non aspects of MOM6. I think my latest commit addresses all your comments, but please take a look.

[Hallberg-NOAA]

With these changes that you just introduced, it is all looking good to me!

[gustavo-marques]

Thanks, @iangrooms, for this contribution, and thanks, @Hallberg-NOAA, for reviewing this PR.

The pr_mom tests are passing. I only have the minor comment raised above regarding a new logical parameter (allow_preserve_variance) in subroutine ALE_remap_velocities.

New MOM_input parameter:

REMAP_VEL_CONSERVE_KE = False !   [Boolean] default = False
                                 ! If true, a correction is applied to the baroclinic component of velocity after
                                 ! remapping so that total KE is conserved. KE may not be conserved when
                                 ! (CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)

New available diagnostics:

 "ale_u2"  [Unused]
     ! modules: {ocean_model,ocean_model_d2}
     ! long_name: Rate of change in half rho0 times depth integral of squared zonal velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then this measures the change before the KE-conserving correction is applied.
     ! units: W m-2
     ! cell_methods: xq:point yh:mean

 "ale_v2"  [Unused]
     ! modules: {ocean_model,ocean_model_d2}
     ! long_name: Rate of change in half rho0 times depth integral of squared meridional velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then this measures the change before the KE-conserving correction is applied.
     ! units: W m-2
     ! cell_methods: xh:mean yq:point

[iangrooms]

When I was first implementing this I noticed that the subroutine ALE_remap_velocities gets called in more than one context; e.g. this subroutine is applied to remap diffusive tendencies. If we just control the action with a MOM_input parameter then the correction either happens every time the subroutine is called or not at all, and we don't necessarily want to apply this KE-conserving correction every time the subroutine is called - only when remapping the velocity during the time stepping loop.

At the time I discussed this with Alistair, who suggested that I add the logical parameter ALE_remap_velocities. I then hard-code the value to 'true' but only in the places where I want the correction to be allowed. In other places it is 'false' by default. (Alistair also helped me significantly by providing code snippets for various parts of the algorithm.)

In short, we have 2 logical parameters controlling this feature. Only one parameter can be set by the user, and it turns the feature on or off. The other parameter is only internal, and serves to ensure that the correction only happens in the places we want it to happen.

[gustavo-marques]

Thanks for the clarification!

[gustavo-marques]

Is this logical parameter necessary?
Looks like CS%conserve_ke is enough to determine if the correction should be applied.