+Unscales area before taking global sum

  Undoes the dimensional scaling of the cell areas before taking their global
sum, so that the reproducing sum does not overflow when there is dimensional
rescaling.  All answers are bitwise identical when there is no rescaling, but
this eliminates a source of inadvertent overflows or underflows in the global
sums, and there is a new optional argument to compute_global_grid_integrals.

src/framework/MOM_spatial_means.F90

@@ -43,7 +43,7 @@ function global_area_mean(var, G, scale)
   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%US%m_to_L**2 * G%IareaT_global)
+  global_area_mean = reproducing_sum(tmpForSumming) * G%IareaT_global
 
 end function global_area_mean
 

src/initialization/MOM_fixed_initialization.F90

@@ -159,7 +159,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
   call initialize_grid_rotation_angle(G, PF)
 
 ! Compute global integrals of grid values for later use in scalar diagnostics !
-  call compute_global_grid_integrals(G)
+  call compute_global_grid_integrals(G, US=US)
 
 ! Write out all of the grid data used by this run.
   if (write_geom) call write_ocean_geometry_file(G, PF, output_dir, US=US)

src/initialization/MOM_shared_initialization.F90

@@ -1145,17 +1145,21 @@ end subroutine set_velocity_depth_min
 ! -----------------------------------------------------------------------------
 !> Pre-compute global integrals of grid quantities (like masked ocean area) for
 !! later use in reporting diagnostics
-subroutine compute_global_grid_integrals(G)
-  type(dyn_horgrid_type), intent(inout) :: G  !< The dynamic horizontal grid
+subroutine compute_global_grid_integrals(G, US)
+  type(dyn_horgrid_type),          intent(inout) :: G  !< The dynamic horizontal grid
+  type(unit_scale_type), optional, intent(in)    :: US !< A dimensional unit scaling type
 
   ! Local variables
   real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
+  real :: area_scale  ! A scaling factor for area into MKS units
   integer :: i,j
 
+  area_scale = 1.0 ; if (present(US)) area_scale = US%L_to_m**2
+
   tmpForSumming(:,:) = 0.
   G%areaT_global = 0.0 ; G%IareaT_global = 0.0
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    tmpForSumming(i,j) = G%areaT(i,j) * G%mask2dT(i,j)
+    tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j)
   enddo ; enddo
   G%areaT_global = reproducing_sum(tmpForSumming)
 

src/initialization/MOM_state_initialization.F90

@@ -1889,16 +1889,19 @@ end subroutine set_velocity_depth_max
 
 !> Subroutine to pre-compute global integrals of grid quantities for
 !! later use in reporting diagnostics
-subroutine compute_global_grid_integrals(G)
+subroutine compute_global_grid_integrals(G, US)
   type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
+  type(unit_scale_type), intent(in)    :: US !< A dimensional unit scaling type
   ! Local variables
   real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
+  real :: area_scale
   integer :: i,j
 
+  area_scale = US%L_to_m**2
   tmpForSumming(:,:) = 0.
   G%areaT_global = 0.0 ; G%IareaT_global = 0.0
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
-    tmpForSumming(i,j) = G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j)
+    tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j)
   enddo ; enddo
   G%areaT_global = reproducing_sum(tmpForSumming)
   G%IareaT_global = 1. / (G%areaT_global)