From 1bf82205a981fa0ae7390d0e513cdef8315bc23c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 16 Dec 2021 08:04:13 -0500 Subject: [PATCH] (*)+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 | 4 - src/tracer/MOM_tracer_flow_control.F90 | 100 +++++++++++++++---------- src/tracer/advection_test_tracer.F90 | 23 +++--- src/tracer/ideal_age_example.F90 | 19 ++--- 4 files changed, 77 insertions(+), 69 deletions(-) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 668c297658..a7cae98620 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/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) diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 2ae72a3270..1941096832 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/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 ! 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. diff --git a/src/tracer/advection_test_tracer.F90 b/src/tracer/advection_test_tracer.F90 index 8fdb525b4a..b37822823a 100644 --- a/src/tracer/advection_test_tracer.F90 +++ b/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 diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90 index d5c813b3d0..5913251b14 100644 --- a/src/tracer/ideal_age_example.F90 +++ b/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