MOM6:(*)+Reproducing tracer stocks

src/core/MOM_checksum_packages.F90

@@ -253,8 +253,8 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe
   real, dimension(G%isc:G%iec, G%jsc:G%jec) :: &
     tmp_A, &  ! The area per cell [m2] (unscaled to permit reproducing sum).
     tmp_V, &  ! The column-integrated volume [m3] (unscaled to permit reproducing sum)
-    tmp_T, &  ! The column-integrated temperature [degC m3]
-    tmp_S     ! The column-integrated salinity [ppt m3]
+    tmp_T, &  ! The column-integrated temperature [degC m3] (unscaled to permit reproducing sum)
+    tmp_S     ! The column-integrated salinity [ppt m3] (unscaled to permit reproducing sum)
   real :: Vol, dV    ! The total ocean volume and its change [m3] (unscaled to permit reproducing sum).
   real :: Area       ! The total ocean surface area [m2] (unscaled to permit reproducing sum).
   real :: h_minimum  ! The minimum layer thicknesses [H ~> m or kg m-2]
@@ -294,6 +294,8 @@ subroutine MOM_state_stats(mesg, u, v, h, Temp, Salt, G, GV, US, allowChange, pe
         T%average = T%average + dV*Temp(i,j,k)
         S%minimum = min( S%minimum, Salt(i,j,k) ) ; S%maximum = max( S%maximum, Salt(i,j,k) )
         S%average = S%average + dV*Salt(i,j,k)
+        tmp_T(i,j) = tmp_T(i,j) + dV*Temp(i,j,k)
+        tmp_S(i,j) = tmp_S(i,j) + dV*Salt(i,j,k)
       endif
       if (h_minimum > h(i,j,k)) h_minimum = h(i,j,k)
     endif

src/diagnostics/MOM_spatial_means.F90

@@ -19,7 +19,7 @@ module MOM_spatial_means
 public :: global_i_mean, global_j_mean
 public :: global_area_mean, global_area_mean_u, global_area_mean_v, global_layer_mean
 public :: global_area_integral
-public :: global_volume_mean, global_mass_integral
+public :: global_volume_mean, global_mass_integral, global_mass_int_EFP
 public :: adjust_area_mean_to_zero
 
 contains
@@ -234,6 +234,49 @@ function global_mass_integral(h, G, GV, var, on_PE_only, scale)
 
 end function global_mass_integral
 
+!> Find the global mass-weighted order invariant integral of a variable in mks units,
+!! returning the value as an EFP_type. This uses reproducing sums.
+function global_mass_int_EFP(h, G, GV, var, on_PE_only, scale)
+  type(ocean_grid_type),   intent(in)  :: G    !< The ocean's grid structure
+  type(verticalGrid_type), intent(in)  :: GV   !< The ocean's vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(in)  :: h    !< Layer thicknesses [H ~> m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                 optional, intent(in)  :: var  !< The variable being integrated
+  logical,       optional, intent(in)  :: on_PE_only  !< If present and true, the sum is only done
+                                               !! on the local PE, but it is still order invariant.
+  real,          optional, intent(in)  :: scale !< A rescaling factor for the variable
+  type(EFP_type) :: global_mass_int_EFP  !< The mass-weighted integral of var (or 1) in
+                                         !! kg times the units of var
+
+  ! Local variables
+  real, dimension(SZI_(G), SZJ_(G)) :: tmpForSum
+  real :: scalefac  ! An overall scaling factor for the areas and variable.
+  integer :: i, j, k, is, ie, js, je, nz, isr, ier, jsr, jer
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+  isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1)
+
+  scalefac = GV%H_to_kg_m2 * G%US%L_to_m**2
+  if (present(scale)) scalefac = scale * scalefac
+
+  tmpForSum(:,:) = 0.0
+  if (present(var)) then
+    do k=1,nz ; do j=js,je ; do i=is,ie
+      tmpForSum(i,j) = tmpForSum(i,j) + var(i,j,k) * &
+                ((scalefac * h(i,j,k)) * (G%areaT(i,j) * G%mask2dT(i,j)))
+    enddo ; enddo ; enddo
+  else
+    do k=1,nz ; do j=js,je ; do i=is,ie
+      tmpForSum(i,j) = tmpForSum(i,j) + &
+                ((scalefac * h(i,j,k)) * (G%areaT(i,j) * G%mask2dT(i,j)))
+    enddo ; enddo ; enddo
+  endif
+
+  global_mass_int_EFP = reproducing_sum_EFP(tmpForSum, isr, ier, jsr, jer, only_on_PE=on_PE_only)
+
+end function global_mass_int_EFP
+
 
 !> Determine the global mean of a field along rows of constant i, returning it
 !! in a 1-d array using the local indexing. This uses reproducing sums.

src/diagnostics/MOM_sum_output.F90

@@ -733,10 +733,6 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
   enddo ; enddo ; enddo
 
   call sum_across_PEs(CS%ntrunc)
-  !   Sum the various quantities across all the processors.  This sum is NOT
-  ! guaranteed to be bitwise reproducible, even on the same decomposition.
-  !   The sum of Tr_stocks should be reimplemented using the reproducing sums.
-  if (nTr_stocks > 0) call sum_across_PEs(Tr_stocks,nTr_stocks)
 
   call max_across_PEs(max_CFL, 2)
 

src/tracer/MOM_CFC_cap.F90

@@ -4,6 +4,7 @@ module MOM_CFC_cap
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use MOM_coms,            only : EFP_type
 use MOM_diag_mediator,   only : diag_ctrl, register_diag_field, post_data
 use MOM_error_handler,   only : MOM_error, FATAL, WARNING
 use MOM_file_parser,     only : get_param, log_param, log_version, param_file_type
@@ -14,6 +15,7 @@ module MOM_CFC_cap
 use MOM_io,              only : vardesc, var_desc, query_vardesc, stdout
 use MOM_open_boundary,   only : ocean_OBC_type
 use MOM_restart,         only : query_initialized, MOM_restart_CS
+use MOM_spatial_means,   only : global_mass_int_EFP
 use MOM_time_manager,    only : time_type
 use time_interp_external_mod, only : init_external_field, time_interp_external
 use MOM_tracer_registry, only : register_tracer, tracer_registry_type
@@ -341,14 +343,13 @@ end subroutine CFC_cap_column_physics
 !> Calculates the mass-weighted integral of all tracer stocks,
 !! returning the number of stocks it has calculated.  If the stock_index
 !! is present, only the stock corresponding to that coded index is returned.
-function CFC_cap_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
+function CFC_cap_stock(h, stocks, G, GV, CS, names, units, stock_index)
   type(ocean_grid_type),           intent(in)    :: G      !< The ocean's grid structure.
   type(verticalGrid_type),         intent(in)    :: GV     !< The ocean's vertical grid structure.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                    intent(in)    :: h      !< Layer thicknesses [H ~> m or kg m-2].
-  real, dimension(:),              intent(out)   :: stocks !< the mass-weighted integrated amount of each
-                                                           !! tracer, in kg times concentration units [kg conc].
-  type(unit_scale_type),           intent(in)    :: US     !< A dimensional unit scaling type
+  type(EFP_type), dimension(:),    intent(out)   :: stocks !< The mass-weighted integrated amount of each
+                                                           !! tracer, in kg times concentration units [kg conc]
   type(CFC_cap_CS),                pointer       :: CS     !< The control structure returned by a
                                                            !! previous call to register_CFC_cap.
   character(len=*), dimension(:),  intent(out)   :: names  !< The names of the stocks calculated.
@@ -357,11 +358,6 @@ function CFC_cap_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
                                                                 !! stock being sought.
   integer                                        :: CFC_cap_stock !< The number of stocks calculated here.
 
-  ! Local variables
-  real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or 1]
-  real :: mass        ! The cell volume or mass [H L2 ~> m3 or kg]
-  integer :: i, j, k, is, ie, js, je, nz
-  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
   CFC_cap_stock = 0
   if (.not.associated(CS)) return
@@ -377,15 +373,8 @@ function CFC_cap_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   call query_vardesc(CS%CFC12_desc, name=names(2), units=units(2), caller="CFC_cap_stock")
   units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg"
 
-  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
-  stocks(1) = 0.0 ; stocks(2) = 0.0
-  do k=1,nz ; do j=js,je ; do i=is,ie
-    mass = G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k)
-    stocks(1) = stocks(1) + CS%CFC11(i,j,k) * mass
-    stocks(2) = stocks(2) + CS%CFC12(i,j,k) * mass
-  enddo ; enddo ; enddo
-  stocks(1) = stock_scale * stocks(1)
-  stocks(2) = stock_scale * stocks(2)
+  stocks(1) = global_mass_int_EFP(h, G, GV, CS%CFC11, on_PE_only=.true.)
+  stocks(2) = global_mass_int_EFP(h, G, GV, CS%CFC12, on_PE_only=.true.)
 
   CFC_cap_stock = 2
 

src/tracer/MOM_OCMIP2_CFC.F90

@@ -3,25 +3,28 @@ module MOM_OCMIP2_CFC
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
-use MOM_coupler_types, only : extract_coupler_type_data, set_coupler_type_data
-use MOM_coupler_types, only : atmos_ocn_coupler_flux
-use MOM_diag_mediator, only : diag_ctrl
-use MOM_error_handler, only : MOM_error, FATAL, WARNING
-use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
-use MOM_forcing_type, only : forcing
-use MOM_hor_index, only : hor_index_type
-use MOM_grid, only : ocean_grid_type
-use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc
-use MOM_open_boundary, only : ocean_OBC_type
-use MOM_restart, only : query_initialized, MOM_restart_CS
-use MOM_sponge, only : set_up_sponge_field, sponge_CS
-use MOM_time_manager, only : time_type
+use MOM_coms,            only : EFP_type
+use MOM_coupler_types,   only : extract_coupler_type_data, set_coupler_type_data
+use MOM_coupler_types,   only : atmos_ocn_coupler_flux
+use MOM_diag_mediator,   only : diag_ctrl
+use MOM_error_handler,   only : MOM_error, FATAL, WARNING
+use MOM_file_parser,     only : get_param, log_param, log_version, param_file_type
+use MOM_forcing_type,    only : forcing
+use MOM_hor_index,       only : hor_index_type
+use MOM_grid,            only : ocean_grid_type
+use MOM_io,              only : file_exists, MOM_read_data, slasher
+use MOM_io,              only : vardesc, var_desc, query_vardesc
+use MOM_open_boundary,   only : ocean_OBC_type
+use MOM_restart,         only : query_initialized, MOM_restart_CS
+use MOM_spatial_means,   only : global_mass_int_EFP
+use MOM_sponge,          only : set_up_sponge_field, sponge_CS
+use MOM_time_manager,    only : time_type
 use MOM_tracer_registry, only : register_tracer, tracer_registry_type
 use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
-use MOM_tracer_Z_init, only : tracer_Z_init
-use MOM_unit_scaling, only : unit_scale_type
-use MOM_variables, only : surface
-use MOM_verticalGrid, only : verticalGrid_type
+use MOM_tracer_Z_init,   only : tracer_Z_init
+use MOM_unit_scaling,    only : unit_scale_type
+use MOM_variables,       only : surface
+use MOM_verticalGrid,    only : verticalGrid_type
 
 implicit none ; private
 
@@ -478,14 +481,13 @@ end subroutine OCMIP2_CFC_column_physics
 !> This function calculates the mass-weighted integral of all tracer stocks,
 !! returning the number of stocks it has calculated.  If the stock_index
 !! is present, only the stock corresponding to that coded index is returned.
-function OCMIP2_CFC_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
+function OCMIP2_CFC_stock(h, stocks, G, GV, CS, names, units, stock_index)
   type(ocean_grid_type),           intent(in)    :: G      !< The ocean's grid structure.
   type(verticalGrid_type),         intent(in)    :: GV     !< The ocean's vertical grid structure.
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                                    intent(in)    :: h      !< Layer thicknesses [H ~> m or kg m-2].
-  real, dimension(:),              intent(out)   :: stocks !< the mass-weighted integrated amount of each
-                                                           !! tracer, in kg times concentration units [kg conc].
-  type(unit_scale_type),           intent(in)    :: US     !< A dimensional unit scaling type
+  type(EFP_type), dimension(:),    intent(out)   :: stocks !< The mass-weighted integrated amount of each
+                                                           !! tracer, in kg times concentration units [kg conc]
   type(OCMIP2_CFC_CS),             pointer       :: CS     !< The control structure returned by a
                                                            !! previous call to register_OCMIP2_CFC.
   character(len=*), dimension(:),  intent(out)   :: names  !< The names of the stocks calculated.
@@ -494,11 +496,6 @@ function OCMIP2_CFC_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
                                                                 !! stock being sought.
   integer                                        :: OCMIP2_CFC_stock !< The number of stocks calculated here.
 
-  ! Local variables
-  real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or 1]
-  real :: mass        ! The cell volume or mass [H L2 ~> m3 or kg]
-  integer :: i, j, k, is, ie, js, je, nz
-  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
   OCMIP2_CFC_stock = 0
   if (.not.associated(CS)) return
@@ -514,15 +511,8 @@ function OCMIP2_CFC_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   call query_vardesc(CS%CFC12_desc, name=names(2), units=units(2), caller="OCMIP2_CFC_stock")
   units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg"
 
-  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
-  stocks(1) = 0.0 ; stocks(2) = 0.0
-  do k=1,nz ; do j=js,je ; do i=is,ie
-    mass = G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k)
-    stocks(1) = stocks(1) + CS%CFC11(i,j,k) * mass
-    stocks(2) = stocks(2) + CS%CFC12(i,j,k) * mass
-  enddo ; enddo ; enddo
-  stocks(1) = stock_scale * stocks(1)
-  stocks(2) = stock_scale * stocks(2)
+  stocks(1) = global_mass_int_EFP(h, G, GV, CS%CFC11, on_PE_only=.true.)
+  stocks(2) = global_mass_int_EFP(h, G, GV, CS%CFC12, on_PE_only=.true.)
 
   OCMIP2_CFC_stock = 2
 

src/tracer/MOM_generic_tracer.F90

@@ -29,7 +29,7 @@ module MOM_generic_tracer
   use g_tracer_utils,   only: g_tracer_get_pointer,g_tracer_get_alias,g_tracer_set_csdiag
 
   use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS
-  use MOM_coms, only : max_across_PEs, min_across_PEs, PE_here
+  use MOM_coms, only : EFP_type, max_across_PEs, min_across_PEs, PE_here
   use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
   use MOM_diag_mediator, only : diag_ctrl, get_diag_time_end
   use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe
@@ -40,7 +40,7 @@ module MOM_generic_tracer
   use MOM_io, only : file_exists, MOM_read_data, slasher
   use MOM_open_boundary, only : ocean_OBC_type
   use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS
-  use MOM_spatial_means, only : global_area_mean
+  use MOM_spatial_means, only : global_area_mean, global_mass_int_EFP
   use MOM_sponge, only : set_up_sponge_field, sponge_CS
   use MOM_time_manager, only : time_type, set_time
   use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
@@ -568,13 +568,12 @@ end subroutine MOM_generic_tracer_column_physics
   !! being requested specifically, returning the number of stocks it has
   !! calculated. If the stock_index is present, only the stock corresponding
   !! to that coded index is returned.
-  function MOM_generic_tracer_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
+  function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index)
     type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure
     type(verticalGrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure
     real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
-    real, dimension(:),                 intent(out)   :: stocks !< The mass-weighted integrated amount of each
-                                                                !! tracer, in kg times concentration units [kg conc].
-    type(unit_scale_type),              intent(in)    :: US     !< A dimensional unit scaling type
+    type(EFP_type), dimension(:),       intent(out)   :: stocks !< The mass-weighted integrated amount of each
+                                                                !! tracer, in kg times concentration units [kg conc]
     type(MOM_generic_tracer_CS),        pointer       :: CS     !< Pointer to the control structure for this module.
     character(len=*), dimension(:),     intent(out)   :: names  !< The names of the stocks calculated.
     character(len=*), dimension(:),     intent(out)   :: units  !< The units of the stocks calculated.
@@ -584,14 +583,12 @@ function MOM_generic_tracer_stock(h, stocks, G, GV, US, CS, names, units, stock_
                                                                      !! number of stocks calculated here.
 
     ! Local variables
-    real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or 1]
     type(g_tracer_type), pointer  :: g_tracer, g_tracer_next
     real, dimension(:,:,:,:), pointer   :: tr_field
     real, dimension(:,:,:), pointer     :: tr_ptr
     character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock'
 
-    integer :: i, j, k, is, ie, js, je, nz, m
-    is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+    integer :: m
 
     MOM_generic_tracer_stock = 0
     if (.not.associated(CS)) return
@@ -605,20 +602,15 @@ function MOM_generic_tracer_stock(h, stocks, G, GV, US, CS, names, units, stock_
 
     if (.NOT. associated(CS%g_tracer_list)) return ! No stocks.
 
-    stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
     m=1 ; g_tracer=>CS%g_tracer_list
     do
       call g_tracer_get_alias(g_tracer,names(m))
       call g_tracer_get_values(g_tracer,names(m),'units',units(m))
       units(m) = trim(units(m))//" kg"
       call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field)
 
-      stocks(m) = 0.0
       tr_ptr => tr_field(:,:,:,1)
-      do k=1,nz ; do j=js,je ; do i=is,ie
-        stocks(m) = stocks(m) + tr_ptr(i,j,k) * (G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k))
-      enddo ; enddo ; enddo
-      stocks(m) = stock_scale * stocks(m)
+      stocks(m) = global_mass_int_EFP(h, G, GV, tr_ptr, on_PE_only=.true.)
 
       !traverse the linked list till hit NULL
       call g_tracer_get_next(g_tracer, g_tracer_next)