+Refactor the spatial mean calculations

  Refactored the spatial mean calculations in the functions global_area_mean(),
global_area_mean_u(), global_area_mean_v(), global_layer_mean(),
global_volume_mean() and adjust_area_mean_to_zero() to work in rescaled units by
making use of the unscale arguments to the reproducing_sum routines.  Comments
were also added to explain what the code does when both tmp_scale and unscale
arguments are provided, which effectively acts as to allow for changes in the
rescaled units.

  Global_area_integral() and global_mass_integral() were also similarly
refactored, but with the added difference that when a tmp_scale argument is
provided, the areas or masses in the integrals have units of [L2 A ~> m2 a] or
[R Z L2 A ~> kg a] instead of the mixed units of [m2 A ~> m2 a] or [kg A ~> kg
a].  As a result the code surrounding the 4 instances where global_area_integral
or global_mass_integral were being called with tmp_scale arguments (in MOM.F90
and MOM_ice_shelf.F90) also had to be modified in this same commit.

  This commit also includes a rescaling in the units of the areaT_global and
IareaT_global elements of the ocean_grid_type and dyn_horgrid_type to [L2 ~> m2]
and [L-2 ~> m-2], respectively.  Although the dyn_horgrid_type is shared between
MOM6 and SIS2, these elements are not used in SIS2.

  A total of 12 rescaling factors were eliminated or moved into unscale arguments
as a result of these changes.

  All answers are bitwise identical, but there are changes in the rescaled units
of two elements each in two transparent types, and changes to the rescaling
behavior of two publicly visible routines when they are called with tmp_scale
arguments.

src/core/MOM.F90

@@ -4161,7 +4161,7 @@ subroutine get_ocean_stocks(CS, mass, heat, salt, on_PE_only)
   if (present(mass)) &
     mass = global_mass_integral(CS%h, CS%G, CS%GV, on_PE_only=on_PE_only)
   if (present(heat)) &
-    heat = CS%US%Q_to_J_kg*CS%tv%C_p * &
+    heat = CS%US%Q_to_J_kg*CS%US%RZL2_to_kg * CS%tv%C_p * &
            global_mass_integral(CS%h, CS%G, CS%GV, CS%tv%T, on_PE_only=on_PE_only, tmp_scale=CS%US%C_to_degC)
   if (present(salt)) &
     salt = 1.0e-3 * global_mass_integral(CS%h, CS%G, CS%GV, CS%tv%S, on_PE_only=on_PE_only, unscale=CS%US%S_to_ppt)

src/core/MOM_grid.F90

@@ -177,9 +177,9 @@ module MOM_grid
     df_dx, &      !< Derivative d/dx f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
     df_dy         !< Derivative d/dy f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
 
-  ! These variables are global sums that are useful for 1-d diagnostics and should not be rescaled.
-  real :: areaT_global  !< Global sum of h-cell area [m2]
-  real :: IareaT_global !< Global sum of inverse h-cell area (1/areaT_global) [m-2].
+  ! These variables are global sums that are useful for 1-d diagnostics.
+  real :: areaT_global  !< Global sum of h-cell area [L2 ~> m2]
+  real :: IareaT_global !< Global sum of inverse h-cell area (1/areaT_global) [L-2 ~> m-2].
 
   type(unit_scale_type), pointer :: US => NULL() !< A dimensional unit scaling type