+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.

  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

@@ -4147,7 +4147,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
 

src/diagnostics/MOM_spatial_means.F90

@@ -55,25 +55,23 @@ function global_area_mean(var, G, scale, tmp_scale, unscale)
   ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
   ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
   real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! An unscaled cell integral [a m2]
-  real :: scalefac   ! An overall scaling factor for the areas and variable [a m2 A-1 L-2 ~> 1]
+  real :: scalefac   ! A scaling factor for the variable that is not reversed [a A-1 ~> 1]
   real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
   integer :: i, j, is, ie, js, je
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
   temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
-  scalefac = G%US%L_to_m**2*temp_scale
-  if (present(unscale)) then ; scalefac = scalefac * unscale
-  elseif (present(scale)) then ; scalefac = scalefac * scale ; endif
+
+  scalefac = 1.0
+  if (present(unscale)) then ; scalefac = unscale
+  elseif (present(scale)) then ; scalefac = scale ; endif
 
   tmpForSumming(:,:) = 0.
   do j=js,je ; do i=is,ie
     tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j))
   enddo ; enddo
 
-  global_area_mean = reproducing_sum(tmpForSumming) * G%IareaT_global
-
-  if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
-    global_area_mean = global_area_mean / temp_scale
+  global_area_mean = reproducing_sum(tmpForSumming, unscale=temp_scale*G%US%L_to_m**2) * G%IareaT_global
 
 end function global_area_mean
 
@@ -94,26 +92,22 @@ function global_area_mean_v(var, G, tmp_scale)
   ! Local variables
   ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
   ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
-  real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! An unscaled cell integral [a m2]
-  real :: scalefac   ! An overall scaling factor for the areas and variable [a m2 A-1 L-2 ~> 1]
+  real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! An unscaled cell integral [A L2 ~> a m2]
   real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
   integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB
 
   temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
-  scalefac = G%US%L_to_m**2*temp_scale
 
   tmpForSumming(:,:) = 0.
   do j=js,je ; do i=is,ie
-    tmpForSumming(i,j) = G%areaT(i,j) * scalefac * &
+    tmpForSumming(i,j) = G%areaT(i,j) * &
              (var(i,J) * G%mask2dCv(i,J) + var(i,J-1) * G%mask2dCv(i,J-1)) / &
              max(1.e-20, G%mask2dCv(i,J)+G%mask2dCv(i,J-1))
   enddo ; enddo
-  global_area_mean_v = reproducing_sum(tmpForSumming) * G%IareaT_global
-  if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
-    global_area_mean_v = global_area_mean_v / temp_scale
+  global_area_mean_v = reproducing_sum(tmpForSumming, unscale=G%US%L_to_m**2*temp_scale) * G%IareaT_global
 
 end function global_area_mean_v
 
@@ -133,31 +127,29 @@ function global_area_mean_u(var, G, tmp_scale)
   ! Local variables
   ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
   ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
-  real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! An unscaled cell integral [a m2]
-  real :: scalefac   ! An overall scaling factor for the areas and variable [a m2 A-1 L-2 ~> 1]
+  real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming ! An unscaled cell integral [A L2 ~> a m2]
   real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
   integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB
 
   temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
-  scalefac = G%US%L_to_m**2*temp_scale
 
   tmpForSumming(:,:) = 0.
   do j=js,je ; do i=is,ie
-    tmpForSumming(i,j) = G%areaT(i,j) * scalefac * &
+    tmpForSumming(i,j) = G%areaT(i,j) * &
              (var(I,j) * G%mask2dCu(I,j) + var(I-1,j) * G%mask2dCu(I-1,j)) / &
              max(1.e-20, G%mask2dCu(I,j)+G%mask2dCu(I-1,j))
   enddo ; enddo
-  global_area_mean_u = reproducing_sum(tmpForSumming) * G%IareaT_global
-  if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
-    global_area_mean_u = global_area_mean_u / temp_scale
+  global_area_mean_u = reproducing_sum(tmpForSumming, unscale=G%US%L_to_m**2*temp_scale) * G%IareaT_global
 
 end function global_area_mean_u
 
 !> Return the global area integral of a variable, by default using the masked area from the
-!! grid, but an alternate could be used instead.  This uses reproducing sums.
+!! grid, but an alternate could be used instead.  The presence of the optional tmp_scale
+!! argument determines whether the returned value is in scaled (if it is present) or unscaled units
+!! for both the variable itself and for the area in the integral.  This function uses reproducing sums.
 function global_area_integral(var, G, scale, area, tmp_scale, unscale)
   type(ocean_grid_type),            intent(in)  :: G     !< The ocean's grid structure
   real, dimension(SZI_(G),SZJ_(G)), intent(in)  :: var   !< The variable to integrate in
@@ -175,7 +167,7 @@ function global_area_integral(var, G, scale, area, tmp_scale, unscale)
                                                          !! Here scale and unscale are synonymous, but unscale
                                                          !! is preferred and takes precedence if both are present.
   real :: global_area_integral !< The returned area integral, usually in the units of var times an area,
-                               !! [a m2] or [A m2 ~> a m2] depending on which optional arguments are provided
+                               !! [a m2] or [A L2 ~> a m2] depending on which optional arguments are provided
 
   ! Local variables
   ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
@@ -186,8 +178,13 @@ function global_area_integral(var, G, scale, area, tmp_scale, unscale)
   integer :: i, j, is, ie, js, je
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
-  temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
-  scalefac = G%US%L_to_m**2*temp_scale
+  if (present(tmp_scale)) then
+    temp_scale = G%US%L_to_m**2 * tmp_scale
+    scalefac = 1.0
+  else
+    temp_scale = 1.0
+    scalefac = G%US%L_to_m**2
+  endif
   if (present(unscale)) then ; scalefac = scalefac * unscale
   elseif (present(scale)) then ; scalefac = scalefac * scale ; endif
 
@@ -202,10 +199,7 @@ function global_area_integral(var, G, scale, area, tmp_scale, unscale)
     enddo ; enddo
   endif
 
-  global_area_integral = reproducing_sum(tmpForSumming)
-
-  if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
-    global_area_integral = global_area_integral / temp_scale
+  global_area_integral = reproducing_sum(tmpForSumming, unscale=temp_scale)
 
 end function global_area_integral
 
@@ -252,12 +246,14 @@ function global_layer_mean(var, h, G, GV, scale, tmp_scale, unscale)
   tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0.
 
   do k=1,nz ; do j=js,je ; do i=is,ie
-    weight(i,j,k)  =  (GV%H_to_MKS * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
+    weight(i,j,k)  =  (GV%H_to_MKS * h(i,j,k)) * (G%areaT(i,j) * G%mask2dT(i,j))
     tmpForSumming(i,j,k) =  scalefac * var(i,j,k) * weight(i,j,k)
   enddo ; enddo ; enddo
 
-  global_temp_scalar = reproducing_sum(tmpForSumming, EFP_lay_sums=laysums(1:nz), only_on_PE=.true.)
-  global_weight_scalar = reproducing_sum(weight, EFP_lay_sums=laysums(nz+1:2*nz), only_on_PE=.true.)
+  global_temp_scalar = reproducing_sum(tmpForSumming, EFP_lay_sums=laysums(1:nz), only_on_PE=.true., &
+                                       unscale=G%US%L_to_m**2)
+  global_weight_scalar = reproducing_sum(weight, EFP_lay_sums=laysums(nz+1:2*nz), only_on_PE=.true., &
+                                       unscale=G%US%L_to_m**2)
   call EFP_sum_across_PEs(laysums, 2*nz)
 
   do k=1,nz
@@ -300,23 +296,26 @@ function global_volume_mean(var, h, G, GV, scale, tmp_scale, unscale)
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
   temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
-  scalefac = temp_scale
+
+  scalefac = 1.0
   if (present(unscale)) then ; scalefac = temp_scale * unscale
   elseif (present(scale)) then ; scalefac = temp_scale * scale ; endif
   tmpForSumming(:,:) = 0. ; sum_weight(:,:) = 0.
 
   do k=1,nz ; do j=js,je ; do i=is,ie
-    weight_here  =  (GV%H_to_MKS * h(i,j,k)) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j))
+    weight_here  =  (GV%H_to_MKS * h(i,j,k)) * (G%areaT(i,j) * G%mask2dT(i,j))
     tmpForSumming(i,j) = tmpForSumming(i,j) + scalefac * var(i,j,k) * weight_here
     sum_weight(i,j) = sum_weight(i,j) + weight_here
   enddo ; enddo ; enddo
-  global_volume_mean = (reproducing_sum(tmpForSumming)) / &
-                       (temp_scale * reproducing_sum(sum_weight))
+  global_volume_mean = (reproducing_sum(tmpForSumming, unscale=temp_scale*G%US%L_to_m**2)) / &
+                       (reproducing_sum(sum_weight, unscale=G%US%L_to_m**2))
 
 end function global_volume_mean
 
 
-!> Find the global mass-weighted integral of a variable. This uses reproducing sums.
+!> Find the global mass-weighted integral of a variable.  The presence of the optional tmp_scale
+!! argument determines whether the returned value is in scaled (if it is present) or unscaled units
+!! for both the variable itself and for the mass in the integral.  This function uses reproducing sums.
 function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale, unscale)
   type(ocean_grid_type),   intent(in)  :: G    !< The ocean's grid structure
   type(verticalGrid_type), intent(in)  :: GV   !< The ocean's vertical grid structure
@@ -337,8 +336,8 @@ function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale, unsca
                                                !! units to enable the use of the reproducing sums.
                                                !! Here scale and unscale are synonymous, but unscale
                                                !! is preferred and takes precedence if both are present.
-  real :: global_mass_integral  !< The mass-weighted integral of var (or 1) in
-                                !! kg times the arbitrary units of var [kg a] or [kg A ~> kg a]
+  real :: global_mass_integral  !< The mass-weighted integral of var (or 1) in kg times the arbitrary
+                                !! units of var [kg a] or in [R Z L2 A ~> kg a] if tmp_scale is present
 
   ! Local variables
   ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
@@ -350,36 +349,38 @@ function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale, unsca
   integer :: i, j, k, is, ie, js, je, nz
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
-  temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
-  scalefac = G%US%L_to_m**2*temp_scale
+  if (present(tmp_scale)) then
+    temp_scale = G%US%RZL2_to_kg * tmp_scale
+    scalefac = 1.0
+  else
+    temp_scale = 1.0
+    scalefac = G%US%RZL2_to_kg
+  endif
   if (present(unscale)) then ; scalefac = scalefac * unscale
   elseif (present(scale)) then ; scalefac = scalefac * scale ; endif
-  tmpForSumming(:,:) = 0.0
 
+  tmpForSumming(:,:) = 0.0
   if (present(var)) then
     do k=1,nz ; do j=js,je ; do i=is,ie
       tmpForSumming(i,j) = tmpForSumming(i,j) + var(i,j,k) * &
-                ((GV%H_to_kg_m2 * h(i,j,k)) * (scalefac*G%areaT(i,j) * G%mask2dT(i,j)))
+                ((GV%H_to_RZ * h(i,j,k)) * (scalefac*G%areaT(i,j) * G%mask2dT(i,j)))
     enddo ; enddo ; enddo
   else
     do k=1,nz ; do j=js,je ; do i=is,ie
       tmpForSumming(i,j) = tmpForSumming(i,j) + &
-                ((GV%H_to_kg_m2 * h(i,j,k)) * (scalefac*G%areaT(i,j) * G%mask2dT(i,j)))
+                ((GV%H_to_RZ * h(i,j,k)) * (scalefac*G%areaT(i,j) * G%mask2dT(i,j)))
     enddo ; enddo ; enddo
   endif
   global_sum = .true. ; if (present(on_PE_only)) global_sum = .not.on_PE_only
   if (global_sum) then
-    global_mass_integral = reproducing_sum(tmpForSumming)
+    global_mass_integral = reproducing_sum(tmpForSumming, unscale=temp_scale)
   else
     global_mass_integral = 0.0
     do j=js,je ; do i=is,ie
       global_mass_integral = global_mass_integral + tmpForSumming(i,j)
     enddo ; enddo
   endif
 
-  if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
-    global_mass_integral = global_mass_integral / temp_scale
-
 end function global_mass_integral
 
 !> Find the global mass-weighted order invariant integral of a variable in mks units,
@@ -654,11 +655,12 @@ subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale, unscale)
   ! Local variables
   ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units of the
   ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
-  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: posVals, negVals ! The positive or negative values in a cell or 0 [a]
-  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: areaXposVals, areaXnegVals ! The cell area integral of the values [m2 a]
+  real :: posVals(G%isc:G%iec, G%jsc:G%jec) ! The positive values in a cell or 0 [A ~> a]
+  real :: negVals(G%isc:G%iec, G%jsc:G%jec) ! The negative values in a cell or 0 [A ~> a]
+  real :: areaXposVals(G%isc:G%iec, G%jsc:G%jec) ! The cell area integral of the positive values [L2 A ~> m2 a]
+  real :: areaXnegVals(G%isc:G%iec, G%jsc:G%jec) ! The cell area integral of the negative values [L2 A ~> m2 a]
   type(EFP_type), dimension(2) :: areaInt_EFP ! An EFP version integral of the values on the current PE [m2 a]
   real :: scalefac  ! A scaling factor for the variable [a A-1 ~> 1]
-  real :: I_scalefac ! The Adcroft reciprocal of scalefac [A a-1 ~> 1]
   real :: areaIntPosVals, areaIntNegVals ! The global area integral of the positive and negative values [m2 a]
   real :: posScale, negScale ! The scaling factor to apply to positive or negative values [nondim]
   integer :: i,j
@@ -667,21 +669,19 @@ subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale, unscale)
   if (present(unscale)) then ; scalefac = unscale
   elseif (present(unit_scale)) then ; scalefac = unit_scale ; endif
 
-  I_scalefac = 0.0 ; if (scalefac /= 0.0) I_scalefac = 1.0 / scalefac
-
   ! areaXposVals(:,:) = 0.  ! This zeros out halo points.
   ! areaXnegVals(:,:) = 0.  ! This zeros out halo points.
 
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    posVals(i,j) = max(0., scalefac*array(i,j))
-    areaXposVals(i,j) = G%US%L_to_m**2*G%areaT(i,j) * posVals(i,j)
-    negVals(i,j) = min(0., scalefac*array(i,j))
-    areaXnegVals(i,j) = G%US%L_to_m**2*G%areaT(i,j) * negVals(i,j)
+    posVals(i,j) = max(0., array(i,j))
+    areaXposVals(i,j) = G%areaT(i,j) * posVals(i,j)
+    negVals(i,j) = min(0., array(i,j))
+    areaXnegVals(i,j) = G%areaT(i,j) * negVals(i,j)
   enddo ; enddo
 
   ! Combining the sums like this avoids separate blocking global sums.
-  areaInt_EFP(1) = reproducing_sum_EFP( areaXposVals, only_on_PE=.true. )
-  areaInt_EFP(2) = reproducing_sum_EFP( areaXnegVals, only_on_PE=.true. )
+  areaInt_EFP(1) = reproducing_sum_EFP( areaXposVals, only_on_PE=.true., unscale=scalefac*G%US%L_to_m**2 )
+  areaInt_EFP(2) = reproducing_sum_EFP( areaXnegVals, only_on_PE=.true., unscale=scalefac*G%US%L_to_m**2 )
   call EFP_sum_across_PEs(areaInt_EFP, 2)
   areaIntPosVals = EFP_to_real( areaInt_EFP(1) )
   areaIntNegVals = EFP_to_real( areaInt_EFP(2) )
@@ -691,12 +691,12 @@ subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale, unscale)
     if (areaIntPosVals>-areaIntNegVals) then ! Scale down positive values
       posScale = - areaIntNegVals / areaIntPosVals
       do j=G%jsc,G%jec ; do i=G%isc,G%iec
-        array(i,j) = ((posScale * posVals(i,j)) + negVals(i,j)) * I_scalefac
+        array(i,j) = (posScale * posVals(i,j)) + negVals(i,j)
       enddo ; enddo
     elseif (areaIntPosVals<-areaIntNegVals) then ! Scale down negative values
       negScale = - areaIntPosVals / areaIntNegVals
       do j=G%jsc,G%jec ; do i=G%isc,G%iec
-        array(i,j) = (posVals(i,j) + (negScale * negVals(i,j))) * I_scalefac
+        array(i,j) = posVals(i,j) + (negScale * negVals(i,j))
       enddo ; enddo
     endif
   endif

src/framework/MOM_dyn_horgrid.F90

@@ -175,9 +175,9 @@ module MOM_dyn_horgrid
     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]
 
   ! These parameters are run-time parameters that are used during some
   ! initialization routines (but not all)

src/ice_shelf/MOM_ice_shelf.F90

@@ -1139,9 +1139,9 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step)
   ! local variables
   real :: frac_shelf       !< The fractional area covered by the ice shelf [nondim].
   real :: frac_open        !< The fractional area of the ocean that is not covered by the ice shelf [nondim].
-  real :: delta_mass_shelf !< Change in ice shelf mass over one time step [R Z m2 T-1 ~> kg s-1]
+  real :: delta_mass_shelf !< Change in ice shelf mass over one time step [R Z L2 T-1 ~> kg s-1]
   real :: balancing_flux   !< The fresh water flux that balances the integrated melt flux [R Z T-1 ~> kg m-2 s-1]
-  real :: balancing_area   !< total area where the balancing flux is applied [m2]
+  real :: balancing_area   !< total area where the balancing flux is applied [L2 ~> m2]
   type(time_type) :: dTime !< The time step as a time_type
   type(time_type) :: Time0 !< The previous time (Time-dt)
   real, dimension(SZDI_(G),SZDJ_(G)) :: bal_frac  !< Fraction of the cell where the mass flux
@@ -1315,7 +1315,7 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step)
       endif
     enddo ; enddo
 
-    balancing_area = global_area_integral(bal_frac, G, area=G%areaT)
+    balancing_area = global_area_integral(bal_frac, G, area=G%areaT, tmp_scale=1.0)
     if (balancing_area > 0.0) then
       balancing_flux = ( global_area_integral(ISS%water_flux, G, tmp_scale=US%RZ_T_to_kg_m2s, &
                                               area=ISS%area_shelf_h) + &