(*)+Reproducing tracer stocks

  Use reproducing sums for tabulating tracer stocks, and move the global sum for
the tracer stocks form write_energy into call_tracer_stocks.  This involves
changes to the type of an argument (from real to EFP_type) for two arguments to
the internal routine store_stocks.  Existing tracer stock packages will still
work, but to benefit from the reproducing sums, they will also have to change
their reported values from real to EFP_type.  This is demonstrated for two
packages (advection_test_tracer and ideal_age_example), where the stocks are now
found with calls to global_mass_int_EFP(), replacing the previous explicit
sums.  With this change, the reported stock values from these packages are
identical for different PE layouts and can be much more accurate than before,
but they are different from the previously reported values at roundoff (for
positive-definite tracers), but it could be larger for tracers with a near-zero
mean value.  All solutions are bitwise identical, but output changes.

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_tracer_flow_control.F90

@@ -3,21 +3,22 @@ module MOM_tracer_flow_control
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use MOM_coms,          only : EFP_type, assignment(=), EFP_to_real, real_to_EFP, EFP_sum_across_PEs
 use MOM_diag_mediator, only : time_type, diag_ctrl
 use MOM_error_handler, only : MOM_error, FATAL, WARNING
-use MOM_file_parser, only : get_param, log_version, param_file_type, close_param_file
-use MOM_forcing_type, only : forcing, optics_type
-use MOM_get_input, only : Get_MOM_input
-use MOM_grid, only : ocean_grid_type
-use MOM_hor_index, only : hor_index_type
+use MOM_file_parser,   only : get_param, log_version, param_file_type, close_param_file
+use MOM_forcing_type,  only : forcing, optics_type
+use MOM_get_input,     only : Get_MOM_input
+use MOM_grid,          only : ocean_grid_type
+use MOM_hor_index,     only : hor_index_type
 use MOM_open_boundary, only : ocean_OBC_type
-use MOM_restart, only : MOM_restart_CS
-use MOM_sponge, only : sponge_CS
-use MOM_ALE_sponge, only : ALE_sponge_CS
+use MOM_restart,       only : MOM_restart_CS
+use MOM_sponge,        only : sponge_CS
+use MOM_ALE_sponge,    only : ALE_sponge_CS
 use MOM_tracer_registry, only : tracer_registry_type
-use MOM_unit_scaling, only : unit_scale_type
-use MOM_variables, only : surface, thermo_var_ptrs
-use MOM_verticalGrid, only : verticalGrid_type
+use MOM_unit_scaling,  only : unit_scale_type
+use MOM_variables,     only : surface, thermo_var_ptrs
+use MOM_verticalGrid,  only : verticalGrid_type
 #include <MOM_memory.h>
 
 ! Add references to other user-provide tracer modules here.
@@ -582,8 +583,8 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, US, CS, stock_names, stock
   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) :: stock_values !< The integrated amounts of a tracer
-                             !! on the current PE, usually in kg x concentration [kg conc].
+  real, dimension(:),             intent(out) :: stock_values !< The globally mass-integrated
+                                                             !! amount of a tracer [kg conc].
   type(unit_scale_type),          intent(in)  :: US          !< A dimensional unit scaling type
   type(tracer_flow_control_CS),   pointer     :: CS          !< The control structure returned by a
                                                              !! previous call to
@@ -612,7 +613,9 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, US, CS, stock_names, stock
   character(len=200), dimension(MAX_FIELDS_) :: names, units
   character(len=200) :: set_pkg_name
   real, dimension(MAX_FIELDS_) :: values
-  integer :: max_ns, ns_tot, ns, index, pkg, max_pkgs, nn
+  type(EFP_type), dimension(MAX_FIELDS_) :: values_EFP
+  type(EFP_type), dimension(MAX_FIELDS_) :: stock_val_EFP
+  integer :: max_ns, ns_tot, ns, index, pkg, max_pkgs, nn, n
 
   if (.not. associated(CS)) call MOM_error(FATAL, "call_tracer_stocks: "// &
        "Module must be initialized via call_tracer_register before it is used.")
@@ -627,57 +630,66 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, US, CS, stock_names, stock
   if (CS%use_USER_tracer_example) then
     ns = USER_tracer_stock(h, values, G, GV, US, CS%USER_tracer_example_CSp, &
                            names, units, stock_index)
-    call store_stocks("tracer_example", ns, names, units, values, index, stock_values, &
+    do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("tracer_example", ns, names, units, values_EFP, index, stock_val_EFP, &
                        set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 ! if (CS%use_DOME_tracer) then
 !   ns = DOME_tracer_stock(h, values, G, GV, CS%DOME_tracer_CSp, &
 !                          names, units, stock_index)
-!   call store_stocks("DOME_tracer", ns, names, units, values, index, stock_values, &
+!   do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+!   call store_stocks("DOME_tracer", ns, names, units, values_EFP, index, stock_val_EFP, &
 !                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
 ! endif
   if (CS%use_ideal_age) then
-    ns = ideal_age_stock(h, values, G, GV, US, CS%ideal_age_tracer_CSp, &
+    ns = ideal_age_stock(h, values_EFP, G, GV, CS%ideal_age_tracer_CSp, &
                          names, units, stock_index)
-    call store_stocks("ideal_age_example", ns, names, units, values, index, &
-           stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    ! do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("ideal_age_example", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
   if (CS%use_regional_dyes) then
     ns = dye_stock(h, values, G, GV, US, CS%dye_tracer_CSp, &
                          names, units, stock_index)
-    call store_stocks("regional_dyes", ns, names, units, values, index, &
-           stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("regional_dyes", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
   if (CS%use_oil) then
     ns = oil_stock(h, values, G, GV, US, CS%oil_tracer_CSp, &
                          names, units, stock_index)
-    call store_stocks("oil_tracer", ns, names, units, values, index, &
-           stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("oil_tracer", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
   if (CS%use_OCMIP2_CFC) then
     ns = OCMIP2_CFC_stock(h, values, G, GV, US, CS%OCMIP2_CFC_CSp, names, units, stock_index)
-    call store_stocks("MOM_OCMIP2_CFC", ns, names, units, values, index, stock_values, &
-                       set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("MOM_OCMIP2_CFC", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
   if (CS%use_CFC_cap) then
     ns = CFC_cap_stock(h, values, G, GV, US, CS%CFC_cap_CSp, names, units, stock_index)
-    call store_stocks("MOM_CFC_cap", ns, names, units, values, index, stock_values, &
-                       set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("MOM_CFC_cap", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
   if (CS%use_advection_test_tracer) then
-    ns = advection_test_stock( h, values, G, GV, US, CS%advection_test_tracer_CSp, &
+    ns = advection_test_stock( h, values_EFP, G, GV, CS%advection_test_tracer_CSp, &
                          names, units, stock_index )
-    call store_stocks("advection_test_tracer", ns, names, units, values, index, &
-           stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    ! do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("advection_test_tracer", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
   if (CS%use_MOM_generic_tracer) then
     ns = MOM_generic_tracer_stock(h, values, G, GV, US, CS%MOM_generic_tracer_CSp, &
                                    names, units, stock_index)
-    call store_stocks("MOM_generic_tracer", ns, names, units, values, index, stock_values, &
-                       set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("MOM_generic_tracer", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
     nn=ns_tot-ns+1
     nn=MOM_generic_tracer_min_max(nn, got_min_max, global_min,  global_max, &
                                   xgmin, ygmin, zgmin, xgmax, ygmax, zgmax ,&
@@ -687,18 +699,26 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, US, CS, stock_names, stock
   if (CS%use_pseudo_salt_tracer) then
     ns = pseudo_salt_stock(h, values, G, GV, US, CS%pseudo_salt_tracer_CSp, &
                          names, units, stock_index)
-    call store_stocks("pseudo_salt_tracer", ns, names, units, values, index, &
-           stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("pseudo_salt_tracer", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
   if (CS%use_boundary_impulse_tracer) then
     ns = boundary_impulse_stock(h, values, G, GV, US, CS%boundary_impulse_tracer_CSp, &
                          names, units, stock_index)
-    call store_stocks("boundary_impulse_tracer", ns, names, units, values, index, &
-           stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+    do n=1,ns ; values_EFP(n) = real_to_EFP(values(n)) ; enddo
+    call store_stocks("boundary_impulse_tracer", ns, names, units, values_EFP, index, stock_val_EFP, &
+                      set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
   endif
 
-  if (ns_tot == 0) stock_values(1) = 0.0
+  !   Sum the various quantities across all the processors.
+  if (ns_tot > 0) then
+    call EFP_sum_across_PEs(stock_val_EFP, ns_tot)
+    do n=1,ns_tot ; stock_values(n) = EFP_to_real(stock_val_EFP(n)) ; enddo
+  else
+    stock_values(1) = 0.0
+  endif
 
   if (present(num_stocks)) num_stocks = ns_tot
 
@@ -713,11 +733,13 @@ subroutine store_stocks(pkg_name, ns, names, units, values, index, stock_values,
                       intent(in)    :: names   !< Diagnostic names to use for each stock.
   character(len=*), dimension(:), &
                       intent(in)    :: units   !< Units to use in the metadata for each stock.
-  real, dimension(:), intent(in)    :: values  !< The values of the tracer stocks
+  type(EFP_type), dimension(:), &
+                      intent(in)    :: values  !< The values of the tracer stocks
   integer,            intent(in)    :: index   !< The integer stock index from
                              !! stocks_constants_mod of the stock to be returned.  If this is
                              !! present and greater than 0, only a single stock can be returned.
-  real, dimension(:), intent(inout) :: stock_values !< The master list of stock values
+  type(EFP_type), dimension(:), &
+                      intent(inout) :: stock_values !< The master list of stock values
   character(len=*),   intent(inout) :: set_pkg_name !< The name of the last tracer package whose
                                                !! stocks were stored for a specific index.  This is
                                                !! used to trigger an error if there are redundant stocks.

src/tracer/advection_test_tracer.F90

@@ -3,16 +3,18 @@ module advection_test_tracer
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use MOM_coms,          only : EFP_type
 use MOM_coupler_types, only : set_coupler_type_data, 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_file_parser,   only : get_param, log_param, log_version, param_file_type
+use MOM_forcing_type,  only : forcing
 use MOM_grid, only : ocean_grid_type
 use MOM_hor_index, only : hor_index_type
 use MOM_io, only : 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_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
@@ -75,8 +77,8 @@ function register_advection_test_tracer(HI, GV, param_file, CS, tr_Reg, restart_
 
   ! Local variables
   character(len=80)  :: name, longname
-! This include declares and sets the variable "version".
-#include "version_variable.h"
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   character(len=40)  :: mdl = "advection_test_tracer" ! This module's name.
   character(len=200) :: inputdir
   character(len=48)  :: flux_units ! The units for tracer fluxes, usually
@@ -344,13 +346,12 @@ end subroutine advection_test_tracer_surface_state
 
 !> Calculate 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 advection_test_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
+function advection_test_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
+  type(EFP_type), 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(advection_test_tracer_CS),     pointer       :: CS     !< The control structure returned by a previous
                                                               !! call to register_advection_test_tracer.
   character(len=*), dimension(:),     intent(out)   :: names  !< the names of the stocks calculated.
@@ -359,7 +360,6 @@ function advection_test_stock(h, stocks, G, GV, US, CS, names, units, stock_inde
   integer                                           :: advection_test_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]
   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
 
@@ -374,14 +374,9 @@ function advection_test_stock(h, stocks, G, GV, US, CS, names, units, stock_inde
     return
   endif ; endif
 
-  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   do m=1,CS%ntr
     call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="advection_test_stock")
-    stocks(m) = 0.0
-    do k=1,nz ; do j=js,je ; do i=is,ie
-      stocks(m) = stocks(m) + CS%tr(i,j,k,m) * (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, CS%tr(:,:,:,m), on_PE_only=.true.)
   enddo
   advection_test_stock = CS%ntr
 

src/tracer/ideal_age_example.F90

@@ -3,6 +3,7 @@ module ideal_age_example
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use MOM_coms,          only : EFP_type
 use MOM_coupler_types, only : set_coupler_type_data, atmos_ocn_coupler_flux
 use MOM_diag_mediator, only : diag_ctrl
 use MOM_error_handler, only : MOM_error, FATAL, WARNING
@@ -13,6 +14,7 @@ module ideal_age_example
 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_spatial_means, only : global_mass_int_EFP
 use MOM_sponge, only : set_up_sponge_field, sponge_CS
 use MOM_time_manager, only : time_type, time_type_to_real
 use MOM_tracer_registry, only : register_tracer, tracer_registry_type
@@ -78,8 +80,8 @@ function register_ideal_age_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
                                                   !! diffusion module
   type(MOM_restart_CS), target, intent(inout) :: restart_CS !< MOM restart control struct
 
-! This include declares and sets the variable "version".
-#include "version_variable.h"
+  ! This include declares and sets the variable "version".
+# include "version_variable.h"
   character(len=40)  :: mdl = "ideal_age_example" ! This module's name.
   character(len=200) :: inputdir ! The directory where the input files are.
   character(len=48)  :: var_name ! The variable's name.
@@ -369,14 +371,13 @@ end subroutine ideal_age_tracer_column_physics
 
 !> Calculates the mass-weighted integral of all tracer stocks, returning the number of stocks it
 !! has calculated.  If stock_index is present, only the stock corresponding to that coded index is found.
-function ideal_age_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
+function ideal_age_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
+  type(EFP_type), 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(ideal_age_tracer_CS),          pointer       :: CS   !< The control structure returned by a previous
                                                             !! call to register_ideal_age_tracer.
   character(len=*), dimension(:),     intent(out)   :: names  !< the names of the stocks calculated.
@@ -386,7 +387,6 @@ function ideal_age_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
   integer                                           :: ideal_age_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]
   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
 
@@ -401,15 +401,10 @@ function ideal_age_stock(h, stocks, G, GV, US, CS, names, units, stock_index)
     return
   endif ; endif
 
-  stock_scale = US%L_to_m**2 * GV%H_to_kg_m2
   do m=1,CS%ntr
     call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="ideal_age_stock")
     units(m) = trim(units(m))//" kg"
-    stocks(m) = 0.0
-    do k=1,nz ; do j=js,je ; do i=is,ie
-      stocks(m) = stocks(m) + CS%tr(i,j,k,m) * (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, CS%tr(:,:,:,m), on_PE_only=.true.)
   enddo
   ideal_age_stock = CS%ntr