From 62d74dfc1ee20edc590c976ec18a076abb026f38 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Dec 2024 14:07:19 -0500 Subject: [PATCH 1/4] *Update defaults for EQN_OF_STATE and 14 parameters Updated the default values of 15 runtime parameters, as agreed upon in a MOM6 consortium wide conversations on July 29, 2024 and December 16, 2024. The most prominent of these is that the default equation of state is now EQN_OF_STATE = "WRIGHT_FULL", in place of the buggy previous default of "WRIGHT". The 8 answer date parameters REGRIDDING_ANSWER_DATE, TIDES_ANSWER_DATE, WAVE_INTERFACE_ANSWER_DATE, MEKE_GM_SRC_ANSWER_DATE, NDIFF_ANSWER_DATE, KPP%ANSWER_DATE, HOR_DIFF_ANSWER_DATE and LOTW_BBL_ANSWER_DATE all now take their default values from DEFAULT_ANSWER_DATE. The bug-retention parameters MEKE_GM_SRC_ALT_SLOPE_BUG, HOR_DIFF_LIMIT_BUG, BACKSCATTER_UNDERBOUND, DETERMINE_TEMP_CONVERGENCE_BUG, LA_MISALIGNMENT_BUG and IDL_HURR_SCM_EDGE_TAPER_BUG are now false by default. The MOM_input files for the test cases in the `.testing/tc[01234]` directories were updated to explicitly set all of these parameters to their previous default values if they are used in the relevant cases, as well as parameters that will be change following discussions from June 2, 2205. Bitwise identical answers are recovered if all 15 of these parameters are set explicitly, but answers can change if any of these parameters take default values. All of these default changes are the consensus decision of consortium-wise MOM6 dev calls. --- .testing/tc0/MOM_input | 2 +- .testing/tc1.a/MOM_tc_variant | 2 +- .testing/tc1.b/MOM_tc_variant | 2 +- .testing/tc1/MOM_input | 14 ++++++++++++++ .testing/tc2/MOM_input | 13 +++++++++++++ .testing/tc3/MOM_input | 5 +++++ .testing/tc4/MOM_input | 5 ++++- src/ALE/MOM_regridding.F90 | 2 +- src/core/MOM_PressureForce_FV.F90 | 9 ++++++--- src/equation_of_state/MOM_EOS.F90 | 2 +- src/parameterizations/lateral/MOM_hor_visc.F90 | 3 +-- .../lateral/MOM_thickness_diffuse.F90 | 8 ++------ src/parameterizations/vertical/MOM_CVMix_KPP.F90 | 2 +- .../vertical/MOM_set_diffusivity.F90 | 3 +-- src/tracer/MOM_neutral_diffusion.F90 | 2 +- src/tracer/MOM_tracer_Z_init.F90 | 6 +----- src/tracer/MOM_tracer_hor_diff.F90 | 9 ++------- src/user/Idealized_Hurricane.F90 | 6 +----- src/user/MOM_wave_interface.F90 | 9 ++------- 19 files changed, 59 insertions(+), 45 deletions(-) diff --git a/.testing/tc0/MOM_input b/.testing/tc0/MOM_input index 7a107486b2..eacf4143de 100644 --- a/.testing/tc0/MOM_input +++ b/.testing/tc0/MOM_input @@ -235,7 +235,7 @@ DIAG_AS_CHKSUM = True DEBUG = True GRID_ROTATION_ANGLE_BUGS = True ! [Boolean] default = True USE_GM_WORK_BUG = True ! [Boolean] default = True -FIX_UNSPLIT_DT_VISC_BUG = False ! [Boolean] default = False USE_LAND_MASK_FOR_HVISC = False ! [Boolean] default = False GUST_CONST = 0.02 ! [Pa] default = 0.02 FIX_USTAR_GUSTLESS_BUG = False ! [Boolean] default = False + diff --git a/.testing/tc1.a/MOM_tc_variant b/.testing/tc1.a/MOM_tc_variant index 26407baf50..ff2dabe065 100644 --- a/.testing/tc1.a/MOM_tc_variant +++ b/.testing/tc1.a/MOM_tc_variant @@ -1,2 +1,2 @@ #override SPLIT=False -#override FIX_UNSPLIT_DT_VISC_BUG = False ! [Boolean] default = False +#override UNSPLIT_DT_VISC_BUG = True ! [Boolean] default = False diff --git a/.testing/tc1.b/MOM_tc_variant b/.testing/tc1.b/MOM_tc_variant index 173196f164..878e582546 100644 --- a/.testing/tc1.b/MOM_tc_variant +++ b/.testing/tc1.b/MOM_tc_variant @@ -1,3 +1,3 @@ #override SPLIT=False #override USE_RK2=True -#override FIX_UNSPLIT_DT_VISC_BUG = False ! [Boolean] default = False +#override UNSPLIT_DT_VISC_BUG = True ! [Boolean] default = False diff --git a/.testing/tc1/MOM_input b/.testing/tc1/MOM_input index 098952ccc2..ea16da62a8 100644 --- a/.testing/tc1/MOM_input +++ b/.testing/tc1/MOM_input @@ -595,3 +595,17 @@ BULKML_CONV_MOMENTUM_BUG = True ! [Boolean] default = True PEN_SW_ABSORB_MINTHICK = 0.001 ! [m] default = 0.001 GUST_CONST = 0.02 ! [Pa] default = 0.02 FIX_USTAR_GUSTLESS_BUG = False ! [Boolean] default = False + + +! Explicitly use the defaults from late 2024 +EQN_OF_STATE = "WRIGHT" ! default = "WRIGHT_FULL" +HOR_DIFF_ANSWER_DATE = 20240101 +HOR_DIFF_LIMIT_BUG = True + +! Explicitly use the defaults from early 2025 +VISC_REM_BUG = True +DRAG_DIFFUSIVITY_ANSWER_DATE = 20250101 +FRICTWORK_BUG = True +HOR_DIFF_ANSWER_DATE = 20240101 +HOR_DIFF_LIMIT_BUG = True +MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP = False diff --git a/.testing/tc2/MOM_input b/.testing/tc2/MOM_input index c8aad58e92..77a2a92678 100644 --- a/.testing/tc2/MOM_input +++ b/.testing/tc2/MOM_input @@ -627,3 +627,16 @@ USE_MLD_ITERATION = False ! [Boolean] default = False PEN_SW_ABSORB_MINTHICK = 0.001 ! [m] default = 0.001 GUST_CONST = 0.02 ! [Pa] default = 0.02 FIX_USTAR_GUSTLESS_BUG = False ! [Boolean] default = False + +! Explicitly use the defaults from late 2024 +EQN_OF_STATE = "WRIGHT" ! default = "WRIGHT_FULL" +TIDES_ANSWER_DATE = 20230630 +NDIFF_ANSWER_DATE = 20240101 +BACKSCATTER_UNDERBOUND = True + +! Explicitly use the defaults from early 2025 +VISC_REM_BUG = True +DRAG_DIFFUSIVITY_ANSWER_DATE = 20250101 +FRICTWORK_BUG = True +NDIFF_ANSWER_DATE = 20240101 +MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP = False diff --git a/.testing/tc3/MOM_input b/.testing/tc3/MOM_input index 6a1238ee96..41ccf286aa 100644 --- a/.testing/tc3/MOM_input +++ b/.testing/tc3/MOM_input @@ -480,3 +480,8 @@ KAPPA_SHEAR_ITER_BUG = True ! [Boolean] default = True KAPPA_SHEAR_ALL_LAYER_TKE_BUG = True ! [Boolean] default = True GUST_CONST = 0.02 ! [Pa] default = 0.02 FIX_USTAR_GUSTLESS_BUG = False ! [Boolean] default = False + +! Explicitly use the defaults from early 2025 +VISC_REM_BUG = True +DRAG_DIFFUSIVITY_ANSWER_DATE = 20250101 +FRICTWORK_BUG = True diff --git a/.testing/tc4/MOM_input b/.testing/tc4/MOM_input index b985b8e082..94ac6a7be8 100644 --- a/.testing/tc4/MOM_input +++ b/.testing/tc4/MOM_input @@ -411,9 +411,12 @@ DEBUG = True INTERPOLATE_RES_FN = True ! [Boolean] default = True GILL_EQUATORIAL_LD = False ! [Boolean] default = False -FIX_UNSPLIT_DT_VISC_BUG = False ! [Boolean] default = False USE_LAND_MASK_FOR_HVISC = False ! [Boolean] default = False KAPPA_SHEAR_ITER_BUG = True ! [Boolean] default = True KAPPA_SHEAR_ALL_LAYER_TKE_BUG = True ! [Boolean] default = True USE_MLD_ITERATION = False ! [Boolean] default = False +! Explicitly use the defaults from early 2025 +VISC_REM_BUG = True +FRICTWORK_BUG = True +MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP = False diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index c29b88286e..926f1f741b 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -313,7 +313,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m "Values below 20190101 result in the use of older, less accurate expressions "//& "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & - default=20181231, do_not_log=.not.GV%Boussinesq) ! ### change to default=default_answer_date) + default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) regrid_answer_date = max(regrid_answer_date, 20230701) call set_regrid_params(CS, regrid_answer_date=regrid_answer_date) endif diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index fffe35104b..1014143965 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -2076,8 +2076,10 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, ADp, SAL default=enable_bugs) call get_param(param_file, mdl, "TIDES", CS%tides, & "If true, apply tidal momentum forcing.", default=.false.) - call get_param(param_file, '', "DEFAULT_ANSWER_DATE", default_answer_date, default=99991231) - if (CS%tides) & + if (CS%tides) then + call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & + "This sets the default value for the various _ANSWER_DATE parameters.", & + default=99991231) call get_param(param_file, mdl, "TIDES_ANSWER_DATE", CS%tides_answer_date, "The vintage of "//& "self-attraction and loading (SAL) and tidal forcing calculations. Setting "//& "dates before 20230701 recovers old answers (Boussinesq and non-Boussinesq "//& @@ -2085,7 +2087,8 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, ADp, SAL "difference is only at bit level and due to a reordered summation. Setting "//& "dates before 20250201 recovers answers (Boussinesq mode) that interface "//& "heights are modified before pressure force integrals are calculated.", & - default=20230630, do_not_log=(.not.CS%tides)) + default=default_answer_date, do_not_log=(.not.CS%tides)) + endif call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, & "If true, calculate self-attraction and loading.", default=CS%tides) if (CS%calculate_SAL) & diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 0227ef78e7..ee49bd282d 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -179,7 +179,7 @@ module MOM_EOS character*(12), parameter :: EOS_ROQUET_RHO_STRING = "ROQUET_RHO" !< A string for specifying the equation of state character*(12), parameter :: EOS_ROQUET_SPV_STRING = "ROQUET_SPV" !< A string for specifying the equation of state character*(12), parameter :: EOS_JACKETT06_STRING = "JACKETT_06" !< A string for specifying the equation of state -character*(12), parameter :: EOS_DEFAULT = EOS_WRIGHT_STRING !< The default equation of state +character*(12), parameter :: EOS_DEFAULT = EOS_WRIGHT_FULL_STRING !< The default equation of state integer, parameter :: TFREEZE_LINEAR = 1 !< A named integer specifying a freezing point expression integer, parameter :: TFREEZE_MILLERO = 2 !< A named integer specifying a freezing point expression diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index c61b4bf54a..61ea69c783 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -2584,8 +2584,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "biharmonic viscosity when no Laplacian viscosity is applied. The default "//& "is true for historical reasons, but this option probably should not be used "//& "because it can contribute to numerical instabilities.", & - default=.true., do_not_log=.not.((CS%better_bound_Kh).and.(CS%better_bound_Ah))) - !### The default for BACKSCATTER_UNDERBOUND should be false. + default=.false., do_not_log=.not.((CS%better_bound_Kh).and.(CS%better_bound_Ah))) call get_param(param_file, mdl, "SMAG_BI_CONST",Smag_bi_const, & "The nondimensional biharmonic Smagorinsky constant, "//& diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index cc21e20e00..1bfa3d340b 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -2193,8 +2193,6 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) ! as the vertical structure of thickness diffusivity. ! Used to determine if FULL_DEPTH_KHTH_MIN should be ! available. - logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to - ! recreate the bugs, or if false bugs are only used if actively selected. logical :: use_meke = .false. ! If true, use the MEKE formulation for the thickness diffusivity. integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. integer :: i, j @@ -2355,14 +2353,12 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) "MEKE_GM_SRC_ALT is true. Values below 20240601 recover the answers from the "//& "original implementation, while higher values use expressions that satisfy "//& "rotational symmetry.", & - default=20240101, do_not_log=.not.CS%GM_src_alt) ! ### Change default to default_answer_date. - call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & - default=.true., do_not_log=.true.) ! This is logged from MOM.F90. + default=default_answer_date, do_not_log=.not.CS%GM_src_alt) call get_param(param_file, mdl, "MEKE_GM_SRC_ALT_SLOPE_BUG", CS%MEKE_src_slope_bug, & "If true, use a bug that limits the positive values, but not the negative values, "//& "of the slopes used when MEKE_GM_SRC_ALT is true. When this is true, it breaks "//& "all of the symmetry rules that MOM6 is supposed to obey.", & - default=enable_bugs, do_not_log=.not.CS%GM_src_alt) + default=.false., do_not_log=.not.CS%GM_src_alt) call get_param(param_file, mdl, "MEKE_GEOMETRIC", CS%MEKE_GEOMETRIC, & "If true, uses the GM coefficient formulation from the GEOMETRIC "//& diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 8e42694b36..bf8de60a27 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -523,7 +523,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive) "The vintage of the order of arithmetic in the CVMix KPP calculations. Values "//& "below 20240501 recover the answers from early in 2024, while higher values "//& "use expressions that have been refactored for rotational symmetry.", & - default=20240101) !### Change to: default=default_answer_date) + default=default_answer_date) call closeParameterBlock(paramFile) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index c87a99ba44..bfee584055 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2468,8 +2468,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ "calculations. Values below 20240630 recover the original answers, while "//& "higher values use more accurate expressions. This only applies when "//& "USE_LOTW_BBL_DIFFUSIVITY is true.", & - default=20190101, do_not_log=.not.CS%use_LOTW_BBL_diffusivity) - !### Set default as default=default_answer_date, or use SET_DIFF_ANSWER_DATE. + default=default_answer_date, do_not_log=.not.CS%use_LOTW_BBL_diffusivity) call get_param(param_file, mdl, "DRAG_DIFFUSIVITY_ANSWER_DATE", CS%drag_diff_answer_date, & "The vintage of the order of arithmetic in the drag diffusivity calculations. "//& "Values above 20250301 use less confusing expressions to set the bottom-drag "//& diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 6e9f9c9f06..32d7ec7291 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -219,7 +219,7 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, "Values of 20240330 or below recover the answers from the original form of the "//& "neutral diffusion code, while higher values use mathematically equivalent "//& "expressions that recover rotational symmetry.", & - default=20240101) !### Change this default later to default_answer_date. + default=default_answer_date) ! Initialize and configure remapping if ( .not.CS%continuous_reconstruction ) then diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 4aacb66409..8b2c12fd9b 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -617,8 +617,6 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, character(len=40) :: mdl = "determine_temperature" ! This subroutine's name. logical :: domore(SZK_(GV)) ! Records which layers need additional iterations logical :: adjust_salt, fit_together, convergence_bug, do_any - logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to - ! recreate the bugs, or if false bugs are only used if actively selected. integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, k, is, ie, js, je, nz, itt @@ -633,14 +631,12 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, "based on the ratio of the thermal and haline coefficients. Otherwise try to "//& "match the density by only adjusting temperatures within a maximum range before "//& "revising estimates of the salinity.", default=.false., do_not_log=just_read) - call get_param(PF, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & - default=.true., do_not_log=.true.) ! This is logged from MOM.F90. call get_param(PF, mdl, "DETERMINE_TEMP_CONVERGENCE_BUG", convergence_bug, & "If true, use layout-dependent tests on the changes in temperature and salinity "//& "to determine when the iterations have converged when DETERMINE_TEMP_ADJUST_T_AND_S "//& "is false. For realistic equations of state and the default values of the "//& "various tolerances, this bug does not impact the solutions.", & - default=enable_bugs, do_not_log=just_read) + default=.false., do_not_log=just_read) call get_param(PF, mdl, "DETERMINE_TEMP_T_MIN", T_min, & "The minimum temperature that can be found by determine_temperature.", & diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index fc6a6f7b1e..9a10826627 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -1641,8 +1641,6 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_tracer_hor_diff" ! This module's name. - logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to - ! recreate the bugs, or if false bugs are only used if actively selected. integer :: default_answer_date if (associated(CS)) then @@ -1718,14 +1716,11 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic "along-isopycnal mixed layer to interior mixing code, while higher values use "//& "mathematically equivalent expressions that recover rotational symmetry "//& "when DIFFUSE_ML_TO_INTERIOR is true.", & - default=20240101, do_not_log=.not.CS%Diffuse_ML_interior) - !### Change the default later to default_answer_date. - call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & - default=.true., do_not_log=.true.) ! This is logged from MOM.F90. + default=default_answer_date, do_not_log=.not.CS%Diffuse_ML_interior) call get_param(param_file, mdl, "HOR_DIFF_LIMIT_BUG", CS%limit_bug, & "If true and the answer date is 20240330 or below, use a rotational symmetry "//& "breaking bug when limiting the tracer properties in tracer_epipycnal_ML_diff.", & - default=enable_bugs, do_not_log=((.not.CS%Diffuse_ML_interior).or.(CS%answer_date>=20240331))) + default=.false., do_not_log=((.not.CS%Diffuse_ML_interior).or.(CS%answer_date>=20240331))) CS%ML_KhTR_scale = 1.0 if (CS%Diffuse_ML_interior) then call get_param(param_file, mdl, "ML_KHTR_SCALE", CS%ML_KhTR_scale, & diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index d8da7553d3..635f94d940 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -132,8 +132,6 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) logical :: continuous_Cd ! If true, use a continuous form for the simple drag coefficient as a ! function of wind speed with the idealized hurricane. When this is false, the ! linear shape for the mid-range wind speeds is specified separately. - logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to - ! recreate the bugs, or if false bugs are only used if actively selected. ! This include declares and sets the variable "version". # include "version_variable.h" @@ -240,13 +238,11 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) call get_param(param_file, mdl, "IDL_HURR_SCM", CS%SCM_mode, & "Single Column mode switch used in the SCM idealized hurricane wind profile.", & default=.false.) - call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & - default=.true., do_not_log=.true.) ! This is logged from MOM.F90. call get_param(param_file, mdl, "IDL_HURR_SCM_EDGE_TAPER_BUG", CS%edge_taper_bug, & "If true and IDL_HURR_SCM is true, use a bug that does all of the tapering and "//& "inflow angle calculations for radii between RAD_EDGE and RAD_AMBIENT as though "//& "they were at RAD_EDGE.", & - default=CS%SCM_mode.and.enable_bugs, do_not_log=.not.CS%SCM_mode) + default=.false., do_not_log=.not.CS%SCM_mode) if (.not.CS%SCM_mode) CS%edge_taper_bug = .false. call get_param(param_file, mdl, "IDL_HURR_SCM_LOCY", CS%dy_from_center, & "Y distance of station used in the SCM idealized hurricane wind profile.", & diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 1d8d00f130..38740dc709 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -293,8 +293,6 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) character*(7), parameter :: COUPLER_STRING = "COUPLER" character*(5), parameter :: INPUT_STRING = "INPUT" integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags - logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to - ! recreate the bugs, or if false bugs are only used if actively selected. logical :: use_waves logical :: StatisticalWaves @@ -337,8 +335,7 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) "\t >= 20230101 - More robust expressions for Update_Stokes_Drift\n"//& "\t >= 20230102 - More robust expressions for get_StokesSL_LiFoxKemper\n"//& "\t >= 20230103 - More robust expressions for ust_2_u10_coare3p5", & - default=20221231, do_not_log=.not.GV%Boussinesq) - !### In due course change the default to default=default_answer_date) + default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701) ! Langmuir number Options @@ -538,12 +535,10 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag) call get_param(param_file, mdl, "LA_MISALIGNMENT", CS%LA_Misalignment, & "Flag (logical) if using misalignment between shear and waves in LA", & default=.false.) - call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & - default=.true., do_not_log=.true.) ! This is logged from MOM.F90. call get_param(param_file, mdl, "LA_MISALIGNMENT_BUG", CS%LA_misalign_bug, & "If true, use a code with a sign error when calculating the misalignment between "//& "the shear and waves when LA_MISALIGNMENT is true.", & - default=CS%LA_Misalignment.and.enable_bugs, do_not_log=.not.CS%LA_Misalignment) + default=.false., do_not_log=.not.CS%LA_Misalignment) call get_param(param_file, mdl, "MIN_LANGMUIR", CS%La_min, & "A minimum value for all Langmuir numbers that is not physical, "//& "but is likely only encountered when the wind is very small and "//& From 05311b5359539a8f600779276f0e41c02ff16720 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 17 Jul 2025 16:19:49 -0400 Subject: [PATCH 2/4] *Update defaults for 6 parameters (VISC_REM_BUG, ) Updated the default values of 6 runtime parameters, as agreed upon in a MOM6 consortium wide conversations on June 2, 2025. VISC_REM_BUG and FRICTWORK_BUG are now false by default. MASS_WEIGHT_IN_PRESSURE_GRADIENT_TOP now takes its default value from the setting for MASS_WEIGHT_IN_PRESSURE_GRADIENT. Similarly, the default for DRAG_DIFFUSIVITY_ANSWER_DATE is now set to follow SET_DIFF_ANSWER_DATE. KELVIN_WAVE_VEL_NUDGING_TIMESCALE is now set to a negative value by default, forcing the user to provide a valid value at runtime. SHELFWAVE_CORRECT_AMPLITUDE now defaults to True. --- src/core/MOM_PressureForce_FV.F90 | 2 +- src/core/MOM_barotropic.F90 | 2 +- src/core/MOM_dynamics_split_RK2.F90 | 2 +- src/core/MOM_dynamics_split_RK2b.F90 | 2 +- src/parameterizations/lateral/MOM_hor_visc.F90 | 2 +- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 3 +-- src/user/Kelvin_initialization.F90 | 3 +-- src/user/shelfwave_initialization.F90 | 2 +- 8 files changed, 8 insertions(+), 10 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 1014143965..b727c30595 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -2131,7 +2131,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, ADp, SAL "If true and MASS_WEIGHT_IN_PRESSURE_GRADIENT is true, use mass weighting when "//& "interpolating T/S for integrals near the top of the water column in FV "//& "pressure gradient calculations. ", & - default=.false.) !### Change Default to MASS_WEIGHT_IN_PRESSURE_GRADIENT? + default=useMassWghtInterp) call get_param(param_file, mdl, "MASS_WEIGHT_IN_PGF_NONBOUS_BUG", MassWghtInterp_NonBous_bug, & "If true, use a masking bug in non-Boussinesq calculations with mass weighting "//& "when interpolating T/S for integrals near the bathymetry in FV pressure "//& diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 88baff2c9b..d611367706 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -5526,7 +5526,7 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & default=.true., do_not_log=.true.) ! This is logged from MOM.F90. - call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, default=enable_bugs, do_not_log=.true.) + call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, default=.false., do_not_log=.true.) call get_param(param_file, mdl, "VISC_REM_BT_WEIGHT_BUG", CS%wt_uv_bug, & "If true, recover a bug in barotropic solver that uses an unnormalized weight "//& "function for vertical averages of baroclinic velocity and forcing. Default "//& diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index d11af637a1..f4c68eb83f 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1488,7 +1488,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p "If true, visc_rem_[uv] in split mode is incorrectly calculated or accounted "//& "for in two places. This parameter controls the defaults of two individual "//& "flags, VISC_REM_TIMESTEP_BUG in MOM_dynamics_split_RK2(b) and "//& - "VISC_REM_BT_WEIGHT_BUG in MOM_barotropic.", default=enable_bugs) + "VISC_REM_BT_WEIGHT_BUG in MOM_barotropic.", default=.false.) call get_param(param_file, mdl, "VISC_REM_TIMESTEP_BUG", CS%visc_rem_dt_bug, & "If true, recover a bug that uses dt_pred rather than dt in "//& "vertvisc_remnant() at the end of predictor stage for the following "//& diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index 9bfbff5191..c77c6d41e2 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -1374,7 +1374,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, "If true, visc_rem_[uv] in split mode is incorrectly calculated or accounted "//& "for in two places. This parameter controls the defaults of two individual "//& "flags, VISC_REM_TIMESTEP_BUG in MOM_dynamics_split_RK2(b) and "//& - "VISC_REM_BT_WEIGHT_BUG in MOM_barotropic.", default=enable_bugs) + "VISC_REM_BT_WEIGHT_BUG in MOM_barotropic.", default=.false.) call get_param(param_file, mdl, "VISC_REM_TIMESTEP_BUG", CS%visc_rem_dt_bug, & "If true, recover a bug that uses dt_pred rather than dt in "//& "vertvisc_remnant() at the end of predictor stage for the following "//& diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 61ea69c783..df7c1e6f86 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -2669,7 +2669,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) call get_param(param_file, mdl, "FRICTWORK_BUG", CS%FrictWork_bug, & "If true, retain an answer-changing bug in calculating the FrictWork, "//& "which cancels the h in thickness flux and the h at velocity point. This is"//& - "not recommended.", default=enable_bugs) + "not recommended.", default=.false.) call get_param(param_file, mdl, "USE_GME", CS%use_GME, & "If true, use the GM+E backscatter scheme in association \n"//& diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index bfee584055..caf0555d28 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2473,8 +2473,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ "The vintage of the order of arithmetic in the drag diffusivity calculations. "//& "Values above 20250301 use less confusing expressions to set the bottom-drag "//& "generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false. ", & - default=20250101, do_not_log=CS%use_LOTW_BBL_diffusivity.or.(CS%BBL_effic<=0.0)) - !### Set default as default=default_answer_date, or use SET_DIFF_ANSWER_DATE. + default=CS%answer_date, do_not_log=CS%use_LOTW_BBL_diffusivity.or.(CS%BBL_effic<=0.0)) CS%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_BBL', diag%axesTi, Time, & 'Bottom Boundary Layer Diffusivity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s) diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 4993651d61..4e035d43b0 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -123,8 +123,7 @@ logical function register_Kelvin_OBC(param_file, CS, US, OBC_Reg) "The timescale with which the inflowing open boundary velocities are nudged toward "//& "their intended values with the Kelvin wave test case, or a negative value to keep "//& "the value that is set when the OBC segments are initialized.", & - units="s", default=1.0/(0.3*86400.), scale=US%s_to_T) - !### Change the default nudging timescale to -1. or another value? + units="s", default=-1.0, scale=US%s_to_T) call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, & default=.true., do_not_log=.true.) ! This is logged from MOM.F90. call get_param(param_file, mdl, "KELVIN_SET_OBC_INDEXING_BUGS", CS%indexing_bugs, & diff --git a/src/user/shelfwave_initialization.F90 b/src/user/shelfwave_initialization.F90 index 1258bd405e..93ed4b2c87 100644 --- a/src/user/shelfwave_initialization.F90 +++ b/src/user/shelfwave_initialization.F90 @@ -81,7 +81,7 @@ function register_shelfwave_OBC(param_file, CS, G, US, OBC_Reg) units="nondim", default=1.) call get_param(param_file, mdl, "SHELFWAVE_CORRECT_AMPLITUDE", CS%shelfwave_correct_amplitude, & "If true, SHELFWAVE_AMPLITUDE gives the actual inflow velocity, rather than giving "//& - "an overall scaling factor for the flow.", default=.false.) !### Make the default .true.? + "an overall scaling factor for the flow.", default=.true.) default_amp = 1.0 ; if (CS%shelfwave_correct_amplitude) default_amp = 0.1 call get_param(param_file, mdl, "SHELFWAVE_AMPLITUDE", CS%my_amp, & "Amplitude of the open boundary current inflows in the shelfwave configuration.", & From 5ac077fadc26cfa7bb12d214cbe3312e0323abba Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 17 Jul 2025 16:28:03 -0400 Subject: [PATCH 3/4] +Obsoleted 7 parameters, including BETTER_BOUND_KH Obsoleted 7 runtime parameters that always take the same value, as agreed upon in a MOM6 consortium wide conversations on June 2, 2025, and added tests to catch these parameters in MOM_obsolete_params. The parameters that were obsoleted include BETTER_BOUND_KH, BETTER_BOUND_AH, USE_DIABATIC_TIME_BUG, FIX_UNSPLIT_DT_VISC_BUG, CFL_BASED_TRUNCATIONS and KD_BACKGROUND_VIA_KDML_BUG. The runtime parameter MAXVEL still exists, but it is logged in a different order than before, changing the contents of the MOM_parameter_doc files. All answers are bitwise identical in cases that run, but cases that use these parameters will experience a fatal error. --- src/core/MOM.F90 | 34 +--- src/core/MOM_dynamics_unsplit.F90 | 30 --- src/core/MOM_dynamics_unsplit_RK2.F90 | 30 --- src/diagnostics/MOM_obsolete_params.F90 | 13 ++ .../lateral/MOM_hor_visc.F90 | 98 ++-------- .../vertical/MOM_bkgnd_mixing.F90 | 29 +-- .../vertical/MOM_vert_friction.F90 | 178 ++++++------------ 7 files changed, 99 insertions(+), 313 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 6e5a4b43d4..dcace89642 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -314,8 +314,6 @@ module MOM logical :: useMEKE !< If true, call the MEKE parameterization. logical :: use_stochastic_EOS !< If true, use the stochastic EOS parameterizations. logical :: useWaves !< If true, update Stokes drift - logical :: use_diabatic_time_bug !< If true, uses the wrong calendar time for diabatic processes, - !! as was done in MOM6 versions prior to February 2018. real :: dtbt_reset_period !< The time interval between dynamic recalculation of the !! barotropic time step [T ~> s]. If this is negative dtbt is never !! calculated, and if it is 0, dtbt is calculated every step. @@ -840,15 +838,9 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS rel_time = 0.0 do n=1,n_max - if (CS%use_diabatic_time_bug) then - ! This wrong form of update was used until Feb 2018, recovered with CS%use_diabatic_time_bug=T. - CS%Time = Time_start + real_to_time(US%T_to_s*int(floor(rel_time+0.5*dt+0.5))) - rel_time = rel_time + dt - else - rel_time = rel_time + dt ! The relative time at the end of the step. - ! Set the universally visible time to the middle of the time step. - CS%Time = Time_start + real_to_time(US%T_to_s*(rel_time - 0.5*dt)) - endif + rel_time = rel_time + dt ! The relative time at the end of the step. + ! Set the universally visible time to the middle of the time step. + CS%Time = Time_start + real_to_time(US%T_to_s*(rel_time - 0.5*dt)) ! Set the local time to the end of the time step. Time_local = Time_start + real_to_time(US%T_to_s*rel_time) @@ -878,16 +870,12 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS endif end_time_thermo = Time_local - if (dtdia > dt .and. .not. CS%use_diabatic_time_bug) then + if (dtdia > dt) then ! If necessary, temporarily reset CS%Time to the center of the period covered ! by the call to step_MOM_thermo, noting that they begin at the same time. - ! This step was missing prior to Feb 2018, and is skipped with CS%use_diabatic_time_bug=T. CS%Time = CS%Time + real_to_time(0.5*US%T_to_s*(dtdia-dt)) - endif - if (dtdia > dt .or. CS%use_diabatic_time_bug) then ! The end-time of the diagnostic interval needs to be set ahead if there ! are multiple dynamic time steps worth of thermodynamics applied here. - ! This line was not conditional prior to Feb 2018, recovered with CS%use_diabatic_time_bug=T. end_time_thermo = Time_local + real_to_time(US%T_to_s*(dtdia-dt)) endif @@ -903,8 +891,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS CS%t_dyn_rel_thermo = -dtdia if (showCallTree) call callTree_waypoint("finished diabatic_first (step_MOM)") - if (dtdia > dt .and. .not. CS%use_diabatic_time_bug) & ! Reset CS%Time to its previous value. - ! This step was missing prior to Feb 2018, recovered with CS%use_diabatic_time_bug=T. + if (dtdia > dt) & ! Reset CS%Time to its previous value. CS%Time = Time_start + real_to_time(US%T_to_s*(rel_time - 0.5*dt)) endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))" @@ -1004,8 +991,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS ! If necessary, temporarily reset CS%Time to the center of the period covered ! by the call to step_MOM_thermo, noting that they end at the same time. - ! This step was missing prior to Feb 2018, and is skipped with CS%use_diabatic_time_bug=T. - if (dtdia > dt .and. .not. CS%use_diabatic_time_bug) & + if (dtdia > dt) & CS%Time = CS%Time - real_to_time(0.5*US%T_to_s*(dtdia-dt)) ! Apply diabatic forcing, do mixing, and regrid. @@ -1024,8 +1010,7 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS endif ! Reset CS%Time to its previous value. - ! This step was missing prior to Feb 2018, and is skipped with CS%use_diabatic_time_bug=T. - if (dtdia > dt .and. .not. CS%use_diabatic_time_bug) & + if (dtdia > dt) & CS%Time = Time_start + real_to_time(US%T_to_s*(rel_time - 0.5*dt)) endif @@ -2733,11 +2718,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & default=default_answer_date, do_not_log=non_Bous) if (non_Bous) CS%answer_date = 99991231 - call get_param(param_file, "MOM", "USE_DIABATIC_TIME_BUG", CS%use_diabatic_time_bug, & - "If true, uses the wrong calendar time for diabatic processes, as was "//& - "done in MOM6 versions prior to February 2018. This is not recommended.", & - default=.false.) - call get_param(param_file, "MOM", "SAVE_INITIAL_CONDS", save_IC, & "If true, write the initial conditions to a file given "//& "by IC_OUTPUT_FILE.", default=.false.) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index d4d3356c3d..658480faba 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -634,9 +634,6 @@ subroutine initialize_dyn_unsplit(u, v, h, tv, Time, G, GV, US, param_file, diag character(len=48) :: flux_units ! This include declares and sets the variable "version". # include "version_variable.h" - logical :: use_correct_dt_visc - logical :: test_value ! This is used to determine whether a logical parameter is being set explicitly. - logical :: explicit_bug, explicit_fix ! These indicate which parameters are set explicitly. integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB @@ -655,33 +652,6 @@ subroutine initialize_dyn_unsplit(u, v, h, tv, Time, G, GV, US, param_file, diag call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "UNSPLIT_DT_VISC_BUG", CS%dt_visc_bug, & - "If false, use the correct timestep in the viscous terms applied in the first "//& - "predictor step with the unsplit time stepping scheme, and in the calculation "//& - "of the turbulent mixed layer properties for viscosity with unsplit or "//& - "unsplit_RK2. If true, an older incorrect value is used.", & - default=.false., do_not_log=.true.) - ! This is used to test whether UNSPLIT_DT_VISC_BUG is being actively set. - call get_param(param_file, mdl, "UNSPLIT_DT_VISC_BUG", test_value, default=.true., do_not_log=.true.) - explicit_bug = CS%dt_visc_bug .eqv. test_value - call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", use_correct_dt_visc, & - "If true, use the correct timestep in the viscous terms applied in the first "//& - "predictor step with the unsplit time stepping scheme, and in the calculation "//& - "of the turbulent mixed layer properties for viscosity with unsplit or "//& - "unsplit_RK2.", default=.true., do_not_log=.true.) - call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", test_value, default=.false., do_not_log=.true.) - explicit_fix = use_correct_dt_visc .eqv. test_value - - if (explicit_bug .and. explicit_fix .and. (use_correct_dt_visc .eqv. CS%dt_visc_bug)) then - ! UNSPLIT_DT_VISC_BUG is being explicitly set, and should not be changed. - call MOM_error(FATAL, "UNSPLIT_DT_VISC_BUG and FIX_UNSPLIT_DT_VISC_BUG are both being set "//& - "with inconsistent values. FIX_UNSPLIT_DT_VISC_BUG is an obsolete "//& - "parameter and should be removed.") - elseif (explicit_fix) then - call MOM_error(WARNING, "FIX_UNSPLIT_DT_VISC_BUG is an obsolete parameter. "//& - "Use UNSPLIT_DT_VISC_BUG instead (noting that it has the opposite sense).") - CS%dt_visc_bug = .not.use_correct_dt_visc - endif - call log_param(param_file, mdl, "UNSPLIT_DT_VISC_BUG", CS%dt_visc_bug, & "If false, use the correct timestep in the viscous terms applied in the first "//& "predictor step with the unsplit time stepping scheme, and in the calculation "//& "of the turbulent mixed layer properties for viscosity with unsplit or "//& diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 7bdae9ab20..e064526a63 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -581,9 +581,6 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, tv, Time, G, GV, US, param_file, character(len=48) :: flux_units ! This include declares and sets the variable "version". # include "version_variable.h" - logical :: use_correct_dt_visc - logical :: test_value ! This is used to determine whether a logical parameter is being set explicitly. - logical :: explicit_bug, explicit_fix ! These indicate which parameters are set explicitly. integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB @@ -618,33 +615,6 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, tv, Time, G, GV, US, param_file, units="nondim", default=0.0) call get_param(param_file, mdl, "UNSPLIT_DT_VISC_BUG", CS%dt_visc_bug, & - "If false, use the correct timestep in the viscous terms applied in the first "//& - "predictor step with the unsplit time stepping scheme, and in the calculation "//& - "of the turbulent mixed layer properties for viscosity with unsplit or "//& - "unsplit_RK2. If true, an older incorrect value is used.", & - default=.false., do_not_log=.true.) - ! This is used to test whether UNSPLIT_DT_VISC_BUG is being explicitly set. - call get_param(param_file, mdl, "UNSPLIT_DT_VISC_BUG", test_value, default=.true., do_not_log=.true.) - explicit_bug = CS%dt_visc_bug .eqv. test_value - call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", use_correct_dt_visc, & - "If true, use the correct timestep in the viscous terms applied in the first "//& - "predictor step with the unsplit time stepping scheme, and in the calculation "//& - "of the turbulent mixed layer properties for viscosity with unsplit or "//& - "unsplit_RK2.", default=.true., do_not_log=.true.) - call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", test_value, default=.false., do_not_log=.true.) - explicit_fix = use_correct_dt_visc .eqv. test_value - - if (explicit_bug .and. explicit_fix .and. (use_correct_dt_visc .eqv. CS%dt_visc_bug)) then - ! UNSPLIT_DT_VISC_BUG is being explicitly set, and should not be changed. - call MOM_error(FATAL, "UNSPLIT_DT_VISC_BUG and FIX_UNSPLIT_DT_VISC_BUG are both being set "//& - "with inconsistent values. FIX_UNSPLIT_DT_VISC_BUG is an obsolete "//& - "parameter and should be removed.") - elseif (explicit_fix) then - call MOM_error(WARNING, "FIX_UNSPLIT_DT_VISC_BUG is an obsolete parameter. "//& - "Use UNSPLIT_DT_VISC_BUG instead (noting that it has the opposite sense).") - CS%dt_visc_bug = .not.use_correct_dt_visc - endif - call log_param(param_file, mdl, "UNSPLIT_DT_VISC_BUG", CS%dt_visc_bug, & "If false, use the correct timestep in the viscous terms applied in the first "//& "predictor step with the unsplit time stepping scheme, and in the calculation "//& "of the turbulent mixed layer properties for viscosity with unsplit or "//& diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index f81f2a7574..bfb621c95c 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -108,10 +108,23 @@ subroutine find_obsolete_params(param_file) call obsolete_real(param_file, "MIN_Z_DIAG_INTERVAL") call obsolete_char(param_file, "Z_OUTPUT_GRID_FILE") + call obsolete_logical(param_file, "CFL_BASED_TRUNCATIONS", .true.) + call obsolete_logical(param_file, "KD_BACKGROUND_VIA_KDML_BUG", .false.) + call obsolete_logical(param_file, "USE_DIABATIC_TIME_BUG", .false.) + call read_param(param_file, "INTERPOLATE_SPONGE_TIME_SPACE", test_logic) call obsolete_logical(param_file, "NEW_SPONGES", warning_val=test_logic, & hint="Use INTERPOLATE_SPONGE_TIME_SPACE instead.") + test_logic = .true. ; call read_param(param_file, "BOUND_KH", test_logic) + call obsolete_logical(param_file, "BETTER_BOUND_KH", warning_val=test_logic, hint="Use BOUND_KH alone.") + test_logic = .true. ; call read_param(param_file, "BOUND_AH", test_logic) + call obsolete_logical(param_file, "BETTER_BOUND_AH", warning_val=test_logic, hint="Use BOUND_AH alone.") + + test_logic = .false. ; call read_param(param_file, "UNSPLIT_DT_VISC_BUG", test_logic) + call obsolete_logical(param_file, "FIX_UNSPLIT_DT_VISC_BUG", warning_val=(.not.test_logic), & + hint="Use UNSPLIT_DT_VISC_BUG instead, but with the reversed meaning.") + call obsolete_logical(param_file, "SMOOTH_RI", hint="Instead use N_SMOOTH_RI.") call obsolete_logical(param_file, "INTERNAL_TIDE_CORNER_ADVECT", .false.) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index df7c1e6f86..965ed12055 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -49,15 +49,11 @@ module MOM_hor_visc !! are not implemented with the biharmonic viscosity. logical :: bound_Kh !< If true, the Laplacian coefficient is locally !! limited to guarantee stability. - logical :: better_bound_Kh !< If true, use a more careful bounding of the - !! Laplacian viscosity to guarantee stability. logical :: EY24_EBT_BS !! If true, use an equivalent barotropic backscatter !! with a stabilizing kill switch in MEKE, !< developed by Yankovsky et al. 2024 logical :: bound_Ah !< If true, the biharmonic coefficient is locally !! limited to guarantee stability. - logical :: better_bound_Ah !< If true, use a more careful bounding of the - !! biharmonic viscosity to guarantee stability. real :: Re_Ah !! If nonzero, the biharmonic coefficient is scaled !< so that the biharmonic Reynolds number is equal to this [nondim]. real :: bound_coef !< The nondimensional coefficient of the ratio of @@ -441,7 +437,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, real :: grid_Kh ! Laplacian viscosity bound by grid [L2 T-1 ~> m2 s-1] real :: grid_Ah ! Biharmonic viscosity bound by grid [L4 T-1 ~> m4 s-1] - logical :: rescale_Kh, legacy_bound + logical :: rescale_Kh logical :: find_FrictWork logical :: apply_OBC = .false. logical :: use_MEKE_Ku @@ -553,9 +549,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, js_vort = js-2 ; je_vort = Jeq+1 ; is_vort = is-2 ; ie_vort = Ieq+1 endif - legacy_bound = (CS%Smagorinsky_Kh .or. CS%Leith_Kh) .and. & - (CS%bound_Kh .and. .not.CS%better_bound_Kh) - if (CS%use_GME) then ! Initialize diagnostic arrays with zeros @@ -672,7 +665,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, & !$OMP is_vort, ie_vort, js_vort, je_vort, & !$OMP is_Kh, ie_Kh, js_Kh, je_Kh, & - !$OMP apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, use_kh_struct, skeb_use_frict, & + !$OMP apply_OBC, rescale_Kh, find_FrictWork, use_kh_struct, skeb_use_frict, & !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, slope_x, slope_y, dz, & !$OMP backscat_subround, GME_effic_h, GME_effic_q, & !$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, & @@ -1121,7 +1114,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - if (CS%better_bound_Ah .or. CS%better_bound_Kh) then + if (CS%bound_Ah .or. CS%bound_Kh) then do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh h_min = min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1)) hrat_min(i,j) = min(1.0, h_min / (h(i,j,k) + h_neglect)) @@ -1176,13 +1169,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - if (legacy_bound) then - ! Older method of bounding for stability - do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh - Kh(i,j) = min(Kh(i,j), CS%Kh_Max_xx(i,j)) - enddo ; enddo - endif - ! Place a floor on the viscosity, if desired. do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min) @@ -1221,7 +1207,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif ! Newer method of bounding for stability - if ((CS%better_bound_Kh) .and. (CS%better_bound_Ah)) then + if ((CS%bound_Kh) .and. (CS%bound_Ah)) then do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh visc_bound_rem(i,j) = 1.0 Kh_max_here = hrat_min(i,j) * CS%Kh_Max_xx(i,j) @@ -1232,7 +1218,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, visc_bound_rem(i,j) = 1.0 - Kh(i,j) / Kh_max_here endif enddo ; enddo - elseif (CS%better_bound_Kh) then + elseif (CS%bound_Kh) then do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh Kh(i,j) = min(Kh(i,j), hrat_min(i,j) * CS%Kh_Max_xx(i,j)) enddo ; enddo @@ -1377,11 +1363,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, endif endif - if (CS%bound_Ah .and. .not. CS%better_bound_Ah) then - do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh - Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xx(i,j)) - enddo ; enddo - endif endif ! Smagorinsky_Ah or Leith_Ah or Leith+E if (use_MEKE_Au) then @@ -1398,8 +1379,8 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - if (CS%better_bound_Ah) then - if (CS%better_bound_Kh) then + if (CS%bound_Ah) then + if (CS%bound_Kh) then do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh Ah(i,j) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * CS%Ah_Max_xx(i,j)) enddo ; enddo @@ -1534,7 +1515,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, / (h_neglect3 + (h2uq + h2vq) * ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J)))) enddo ; enddo - if (CS%better_bound_Ah .or. CS%better_bound_Kh) then + if (CS%bound_Ah .or. CS%bound_Kh) then do J=js-1,Jeq ; do I=is-1,Ieq h_min = min(h_u(I,j), h_u(I,j+1), h_v(i,J), h_v(i+1,J)) hrat_min(I,J) = min(1.0, h_min / (hq(I,J) + h_neglect)) @@ -1627,13 +1608,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - if (legacy_bound) then - ! Older method of bounding for stability - do J=js-1,Jeq ; do I=is-1,Ieq - Kh(I,J) = min(Kh(I,J), CS%Kh_Max_xy(I,J)) - enddo ; enddo - endif - do J=js-1,Jeq ; do I=is-1,Ieq Kh(I,J) = max(Kh(I,J), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired. enddo ; enddo @@ -1671,7 +1645,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, do J=js-1,Jeq ; do I=is-1,Ieq ! Newer method of bounding for stability - if ((CS%better_bound_Kh) .and. (CS%better_bound_Ah)) then + if ((CS%bound_Kh) .and. (CS%bound_Ah)) then visc_bound_rem(I,J) = 1.0 Kh_max_here = hrat_min(I,J) * CS%Kh_Max_xy(I,J) if (Kh(I,J) >= Kh_max_here) then @@ -1680,7 +1654,7 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, elseif ((Kh(I,J) > 0.0) .or. (CS%backscatter_underbound .and. (Kh_max_here > 0.0))) then visc_bound_rem(I,J) = 1.0 - Kh(I,J) / Kh_max_here endif - elseif (CS%better_bound_Kh) then + elseif (CS%bound_Kh) then Kh(I,J) = min(Kh(I,J), hrat_min(I,J) * CS%Kh_Max_xy(I,J)) endif enddo ; enddo @@ -1764,12 +1738,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, Ah(I,J) = max(Ah(I,J), AhLth) enddo ; enddo endif - - if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then - do J=js-1,Jeq ; do I=is-1,Ieq - Ah(I,J) = min(Ah(I,J), CS%Ah_Max_xy(I,J)) - enddo ; enddo - endif endif ! Smagorinsky_Ah or Leith_Ah if (use_MEKE_Au) then @@ -1787,8 +1755,8 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, enddo ; enddo endif - if (CS%better_bound_Ah) then - if (CS%better_bound_Kh) then + if (CS%bound_Ah) then + if (CS%bound_Kh) then do J=js-1,Jeq ; do I=is-1,Ieq Ah(I,J) = min(Ah(I,J), visc_bound_rem(I,J) * hrat_min(I,J) * CS%Ah_Max_xy(I,J)) enddo ; enddo @@ -2484,16 +2452,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) call get_param(param_file, mdl, "BOUND_KH", CS%bound_Kh, & "If true, the Laplacian coefficient is locally limited "//& "to be stable.", default=.true., do_not_log=.not.CS%Laplacian) - call get_param(param_file, mdl, "BETTER_BOUND_KH", CS%better_bound_Kh, & - "If true, the Laplacian coefficient is locally limited "//& - "to be stable with a better bounding than just BOUND_KH.", & - default=CS%bound_Kh, do_not_log=.not.CS%Laplacian) call get_param(param_file, mdl, "EY24_EBT_BS", CS%EY24_EBT_BS, & "If true, use the the backscatter scheme (EBT mode with kill switch)"//& "developed by Yankovsky et al. (2024). ", & default=.false., do_not_log=.not.CS%Laplacian) if (.not.CS%Laplacian) CS%bound_Kh = .false. - if (.not.CS%Laplacian) CS%better_bound_Kh = .false. if (.not.(CS%Laplacian.and.use_MEKE)) CS%EY24_EBT_BS = .false. call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, & "If true, allow anistropic viscosity in the Laplacian "//& @@ -2566,12 +2529,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) call get_param(param_file, mdl, "BOUND_AH", CS%bound_Ah, & "If true, the biharmonic coefficient is locally limited "//& "to be stable.", default=.true., do_not_log=.not.CS%biharmonic) - call get_param(param_file, mdl, "BETTER_BOUND_AH", CS%better_bound_Ah, & - "If true, the biharmonic coefficient is locally limited "//& - "to be stable with a better bounding than just BOUND_AH.", & - default=CS%bound_Ah, do_not_log=.not.CS%biharmonic) if (.not.CS%biharmonic) CS%bound_Ah = .false. - if (.not.CS%biharmonic) CS%better_bound_Ah = .false. call get_param(param_file, mdl, "RE_AH", CS%Re_Ah, & "If nonzero, the biharmonic coefficient is scaled "//& "so that the biharmonic Reynolds number is equal to this.", & @@ -2584,7 +2542,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "biharmonic viscosity when no Laplacian viscosity is applied. The default "//& "is true for historical reasons, but this option probably should not be used "//& "because it can contribute to numerical instabilities.", & - default=.false., do_not_log=.not.((CS%better_bound_Kh).and.(CS%better_bound_Ah))) + default=.false., do_not_log=.not.((CS%bound_Kh).and.(CS%bound_Ah))) call get_param(param_file, mdl, "SMAG_BI_CONST",Smag_bi_const, & "The nondimensional biharmonic Smagorinsky constant, "//& @@ -2643,7 +2601,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) "The nondimensional coefficient of the ratio of the "//& "viscosity bounds to the theoretical maximum for "//& "stability without considering other terms.", units="nondim", & - default=0.8, do_not_log=.not.(CS%better_bound_Ah .or. CS%better_bound_Kh)) + default=0.8, do_not_log=.not.(CS%bound_Ah .or. CS%bound_Kh)) call get_param(param_file, mdl, "KILL_SWITCH_COEF", CS%KS_coef, & "A nondimensional coefficient on the biharmonic viscosity that "// & "sets the kill switch for backscatter. Default is 1.0.", units="nondim", & @@ -2743,7 +2701,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%grid_sp_h2(isd:ied,jsd:jed)) ; CS%grid_sp_h2(:,:) = 0.0 ALLOC_(CS%Kh_bg_xx(isd:ied,jsd:jed)) ; CS%Kh_bg_xx(:,:) = 0.0 ALLOC_(CS%Kh_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_bg_xy(:,:) = 0.0 - if (CS%bound_Kh .or. CS%better_bound_Kh .or. CS%EY24_EBT_BS) then + if (CS%bound_Kh .or. CS%EY24_EBT_BS) then allocate(CS%Kh_Max_xx(Isd:Ied,Jsd:Jed), source=0.0) allocate(CS%Kh_Max_xy(IsdB:IedB,JsdB:JedB), source=0.0) endif @@ -2801,7 +2759,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) ALLOC_(CS%Ah_bg_xx(isd:ied,jsd:jed)) ; CS%Ah_bg_xx(:,:) = 0.0 ALLOC_(CS%Ah_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_bg_xy(:,:) = 0.0 ALLOC_(CS%grid_sp_h3(isd:ied,jsd:jed)) ; CS%grid_sp_h3(:,:) = 0.0 - if (CS%bound_Ah .or. CS%better_bound_Ah) then + if (CS%bound_Ah) then allocate(CS%Ah_Max_xx(isd:ied,jsd:jed), source=0.0) allocate(CS%Ah_Max_xy(IsdB:IedB,JsdB:JedB), source=0.0) endif @@ -2908,11 +2866,6 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) slat_fn = abs( sin( deg2rad * G%geoLatT(i,j) ) ) ** Kh_pwr_of_sine CS%Kh_bg_xx(i,j) = MAX(Kh_sin_lat * slat_fn, CS%Kh_bg_xx(i,j)) endif - if (CS%bound_Kh .and. .not.CS%better_bound_Kh) then - ! Limit the background viscosity to be numerically stable - CS%Kh_Max_xx(i,j) = Kh_Limit * grid_sp_h2 - CS%Kh_bg_xx(i,j) = MIN(CS%Kh_bg_xx(i,j), CS%Kh_Max_xx(i,j)) - endif min_grid_sp_h2 = min(grid_sp_h2, min_grid_sp_h2) enddo ; enddo call min_across_PEs(min_grid_sp_h2) @@ -2943,11 +2896,6 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) slat_fn = abs( sin( deg2rad * G%geoLatBu(I,J) ) ) ** Kh_pwr_of_sine CS%Kh_bg_xy(I,J) = MAX(Kh_sin_lat * slat_fn, CS%Kh_bg_xy(I,J)) endif - if (CS%bound_Kh .and. .not.CS%better_bound_Kh) then - ! Limit the background viscosity to be numerically stable - CS%Kh_Max_xy(I,J) = Kh_Limit * grid_sp_q2 - CS%Kh_bg_xy(I,J) = MIN(CS%Kh_bg_xy(I,J), CS%Kh_Max_xy(I,J)) - endif enddo ; enddo endif if (CS%biharmonic) then @@ -2962,7 +2910,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) CS%Ah_bg_xy(:,:) = 0.0 ! The 0.3 below was 0.4 in HIM 1.10. The change in hq requires ! this to be less than 1/3, rather than 1/2 as before. - if (CS%better_bound_Ah .or. CS%bound_Ah) Ah_Limit = 0.3 / (dt*64.0) + if (CS%bound_Ah) Ah_Limit = 0.3 / (dt*64.0) if (CS%Smagorinsky_Ah .and. CS%bound_Coriolis) & BoundCorConst = 1.0 / (5.0*(bound_Cor_vel*bound_Cor_vel)) @@ -2992,10 +2940,6 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xx(i,j) = grid_sp_h3 / CS%Re_Ah if (Ah_time_scale > 0.) CS%Ah_bg_xx(i,j) = & MAX(CS%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / Ah_time_scale) - if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then - CS%Ah_Max_xx(i,j) = Ah_Limit * (grid_sp_h2 * grid_sp_h2) - CS%Ah_bg_xx(i,j) = MIN(CS%Ah_bg_xx(i,j), CS%Ah_Max_xx(i,j)) - endif min_grid_sp_h4 = min(grid_sp_h2**2, min_grid_sp_h4) enddo ; enddo call min_across_PEs(min_grid_sp_h4) @@ -3017,14 +2961,10 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xy(i,j) = grid_sp_q3 / CS%Re_Ah if (Ah_time_scale > 0.) CS%Ah_bg_xy(i,j) = & MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale) - if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then - CS%Ah_Max_xy(I,J) = Ah_Limit * (grid_sp_q2 * grid_sp_q2) - CS%Ah_bg_xy(I,J) = MIN(CS%Ah_bg_xy(I,J), CS%Ah_Max_xy(I,J)) - endif enddo ; enddo endif ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1. - if (CS%Laplacian .and. CS%better_bound_Kh) then + if (CS%Laplacian .and. CS%bound_Kh) then do j=js-1,Jeq+1 ; do i=is-1,Ieq+1 denom = max( & (CS%dy2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j)) * & @@ -3052,7 +2992,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp) endif ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but ! empirically work for CS%bound_coef <~ 1.0 - if (CS%biharmonic .and. CS%better_bound_Ah) then + if (CS%biharmonic .and. CS%bound_Ah) then do j=js-1,Jeq+1 ; do I=is-2,Ieq+1 u0u(I,j) = ((CS%Idxdy2u(I,j)*((CS%dy2h(i+1,j)*CS%DY_dxT(i+1,j)*(G%IdyCu(I+1,j) + G%IdyCu(I,j))) + & (CS%dy2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j))) )) + & diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index b42bd3a8ad..d6c51201a6 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -85,9 +85,6 @@ module MOM_bkgnd_mixing !! here is to assume that the in-situ stratification is the same as the reference stratificaiton. logical :: physical_OBL_scheme !< If true, a physically-based scheme is used to determine mixing in the !! ocean's surface boundary layer, such as ePBL, KPP, or a refined bulk mixed layer scheme. - logical :: Kd_via_Kdml_bug !< If true and KDML /= KD and a number of other higher precedence - !! options are not used, the background diffusivity is set incorrectly using a - !! bug that was introduced in March, 2018. logical :: debug !< If true, turn on debugging in this module ! Diagnostic handles and pointers type(diag_ctrl), pointer :: diag => NULL() !< A structure that regulates diagnostic output @@ -305,16 +302,6 @@ subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS, physical_OBL if (CS%Henyey_IGW_background .and. CS%Kd_tanh_lat_fn) call MOM_error(FATAL, & "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.") - CS%Kd_via_Kdml_bug = .false. - if ((CS%Kd /= CS%Kd_tot_ml) .and. .not.(CS%Kd_tanh_lat_fn .or. CS%physical_OBL_scheme .or. & - CS%Henyey_IGW_background .or. & - CS%horiz_varying_background .or. CS%Bryan_Lewis_diffusivity)) then - call get_param(param_file, mdl, "KD_BACKGROUND_VIA_KDML_BUG", CS%Kd_via_Kdml_bug, & - "If true and KDML /= KD and several other conditions apply, the background "//& - "diffusivity is set incorrectly using a bug that was introduced in March, 2018.", & - default=.false.) ! This parameter should be obsoleted. - endif - ! call closeParameterBlock(param_file) end subroutine bkgnd_mixing_init @@ -479,20 +466,10 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G, do i=is,ie ; depth(i) = 0.0 ; enddo do k=1,nz ; do i=is,ie depth_c = depth(i) + 0.5*h(i,j,k) - if (CS%Kd_via_Kdml_bug) then - ! These two lines should update Kd_lay, not Kd_int. They were correctly working on the - ! same variables until MOM6 commit 7a818716 (PR#750), which was added on March 26, 2018. - if (depth_c <= CS%Hmix) then ; Kd_int(i,K) = CS%Kd_tot_ml - elseif (depth_c >= 2.0*CS%Hmix) then ; Kd_int(i,K) = Kd_sfc(i) - else - Kd_lay(i,k) = ((Kd_sfc(i) - CS%Kd_tot_ml) * I_Hmix) * depth_c + (2.0*CS%Kd_tot_ml - Kd_sfc(i)) - endif + if (depth_c <= CS%Hmix) then ; Kd_lay(i,k) = CS%Kd_tot_ml + elseif (depth_c >= 2.0*CS%Hmix) then ; Kd_lay(i,k) = Kd_sfc(i) else - if (depth_c <= CS%Hmix) then ; Kd_lay(i,k) = CS%Kd_tot_ml - elseif (depth_c >= 2.0*CS%Hmix) then ; Kd_lay(i,k) = Kd_sfc(i) - else - Kd_lay(i,k) = ((Kd_sfc(i) - CS%Kd_tot_ml) * I_Hmix) * depth_c + (2.0*CS%Kd_tot_ml - Kd_sfc(i)) - endif + Kd_lay(i,k) = ((Kd_sfc(i) - CS%Kd_tot_ml) * I_Hmix) * depth_c + (2.0*CS%Kd_tot_ml - Kd_sfc(i)) endif depth(i) = depth(i) + h(i,j,k) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index c5a085298e..c14b2f0052 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -77,7 +77,6 @@ module MOM_vert_friction !! viscosity via Kv_gl90 = alpha_gl90 * f^2. Note that the implied !! Kv_gl90 corresponds to a kappa_gl90 that scales as N^2 with depth. !! [H Z T ~> m2 s or kg s m-1] - real :: maxvel !< Velocity components greater than maxvel are truncated [L T-1 ~> m s-1]. real :: vel_underflow !< Velocity components smaller than vel_underflow !! are set to 0 [L T-1 ~> m s-1]. logical :: CFL_based_trunc !< If true, base truncations on CFL numbers, not @@ -2944,9 +2943,6 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure ! Local variables - - real :: maxvel ! Velocities components greater than maxvel are truncated [L T-1 ~> m s-1] - real :: truncvel ! The speed to which velocity components greater than maxvel are set [L T-1 ~> m s-1] real :: CFL ! The local CFL number [nondim] real :: H_report ! A thickness below which not to report truncations [H ~> m or kg m-2] real :: vel_report(SZIB_(G),SZJB_(G)) ! The velocity to report [L T-1 ~> m s-1] @@ -2957,8 +2953,6 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - maxvel = CS%maxvel - truncvel = 0.9*maxvel H_report = 6.0 * GV%Angstrom_H if (len_trim(CS%u_trunc_file) > 0) then @@ -2966,36 +2960,26 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS do j=js,je trunc_any = .false. do I=Isq,Ieq ; dowrite(I,j) = .false. ; enddo - if (CS%CFL_based_trunc) then - do I=Isq,Ieq ; vel_report(i,j) = 3.0e8*US%m_s_to_L_T ; enddo ! Speed of light default. - do k=1,nz ; do I=Isq,Ieq - if (abs(u(I,j,k)) < CS%vel_underflow) u(I,j,k) = 0.0 - if (u(I,j,k) < 0.0) then - CFL = (-u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i+1,j)) - else - CFL = (u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j)) - endif - if (CFL > CS%CFL_trunc) trunc_any = .true. - if (CFL > CS%CFL_report) then - dowrite(I,j) = .true. - vel_report(I,j) = MIN(vel_report(I,j), abs(u(I,j,k))) - endif - enddo ; enddo - else - do I=Isq,Ieq; vel_report(I,j) = maxvel; enddo - do k=1,nz ; do I=Isq,Ieq - if (abs(u(I,j,k)) < CS%vel_underflow) then ; u(I,j,k) = 0.0 - elseif (abs(u(I,j,k)) > maxvel) then - dowrite(I,j) = .true. ; trunc_any = .true. - endif - enddo ; enddo - endif + do I=Isq,Ieq ; vel_report(i,j) = 3.0e8*US%m_s_to_L_T ; enddo ! Speed of light default. + do k=1,nz ; do I=Isq,Ieq + if (abs(u(I,j,k)) < CS%vel_underflow) u(I,j,k) = 0.0 + if (u(I,j,k) < 0.0) then + CFL = (-u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i+1,j)) + else + CFL = (u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j)) + endif + if (CFL > CS%CFL_trunc) trunc_any = .true. + if (CFL > CS%CFL_report) then + dowrite(I,j) = .true. + vel_report(I,j) = MIN(vel_report(I,j), abs(u(I,j,k))) + endif + enddo ; enddo do I=Isq,Ieq ; if (dowrite(I,j)) then u_old(I,j,:) = u(I,j,:) endif ; enddo - if (trunc_any) then ; if (CS%CFL_based_trunc) then + if (trunc_any) then do k=1,nz ; do I=Isq,Ieq if ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i+1,j) < -CS%CFL_trunc) then u(I,j,k) = (-0.9*CS%CFL_trunc) * (G%areaT(i+1,j) / (dt * G%dy_Cu(I,j))) @@ -3005,36 +2989,20 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS if (h(i,j,k) + h(i+1,j,k) > H_report) CS%ntrunc = CS%ntrunc + 1 endif enddo ; enddo - else - do k=1,nz ; do I=Isq,Ieq ; if (abs(u(I,j,k)) > maxvel) then - u(I,j,k) = SIGN(truncvel,u(I,j,k)) - if (h(i,j,k) + h(i+1,j,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - endif ; enddo ; enddo - endif ; endif + endif enddo ! j-loop else ! Do not report accelerations leading to large velocities. - if (CS%CFL_based_trunc) then - !$OMP parallel do default(shared) - do k=1,nz ; do j=js,je ; do I=Isq,Ieq - if (abs(u(I,j,k)) < CS%vel_underflow) then ; u(I,j,k) = 0.0 - elseif ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i+1,j) < -CS%CFL_trunc) then - u(I,j,k) = (-0.9*CS%CFL_trunc) * (G%areaT(i+1,j) / (dt * G%dy_Cu(I,j))) - if (h(i,j,k) + h(i+1,j,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - elseif ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i,j) > CS%CFL_trunc) then - u(I,j,k) = (0.9*CS%CFL_trunc) * (G%areaT(i,j) / (dt * G%dy_Cu(I,j))) - if (h(i,j,k) + h(i+1,j,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - endif - enddo ; enddo ; enddo - else - !$OMP parallel do default(shared) - do k=1,nz ; do j=js,je ; do I=Isq,Ieq - if (abs(u(I,j,k)) < CS%vel_underflow) then ; u(I,j,k) = 0.0 - elseif (abs(u(I,j,k)) > maxvel) then - u(I,j,k) = SIGN(truncvel, u(I,j,k)) - if (h(i,j,k) + h(i+1,j,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - endif - enddo ; enddo ; enddo - endif + !$OMP parallel do default(shared) + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + if (abs(u(I,j,k)) < CS%vel_underflow) then ; u(I,j,k) = 0.0 + elseif ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i+1,j) < -CS%CFL_trunc) then + u(I,j,k) = (-0.9*CS%CFL_trunc) * (G%areaT(i+1,j) / (dt * G%dy_Cu(I,j))) + if (h(i,j,k) + h(i+1,j,k) > H_report) CS%ntrunc = CS%ntrunc + 1 + elseif ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i,j) > CS%CFL_trunc) then + u(I,j,k) = (0.9*CS%CFL_trunc) * (G%areaT(i,j) / (dt * G%dy_Cu(I,j))) + if (h(i,j,k) + h(i+1,j,k) > H_report) CS%ntrunc = CS%ntrunc + 1 + endif + enddo ; enddo ; enddo endif if (len_trim(CS%u_trunc_file) > 0) then @@ -3050,36 +3018,26 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS do J=Jsq,Jeq trunc_any = .false. do i=is,ie ; dowrite(i,J) = .false. ; enddo - if (CS%CFL_based_trunc) then - do i=is,ie ; vel_report(i,J) = 3.0e8*US%m_s_to_L_T ; enddo ! Speed of light default. - do k=1,nz ; do i=is,ie - if (abs(v(i,J,k)) < CS%vel_underflow) v(i,J,k) = 0.0 - if (v(i,J,k) < 0.0) then - CFL = (-v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1)) - else - CFL = (v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j)) - endif - if (CFL > CS%CFL_trunc) trunc_any = .true. - if (CFL > CS%CFL_report) then - dowrite(i,J) = .true. - vel_report(i,J) = MIN(vel_report(i,J), abs(v(i,J,k))) - endif - enddo ; enddo - else - do i=is,ie ; vel_report(i,J) = maxvel ; enddo - do k=1,nz ; do i=is,ie - if (abs(v(i,J,k)) < CS%vel_underflow) then ; v(i,J,k) = 0.0 - elseif (abs(v(i,J,k)) > maxvel) then - dowrite(i,J) = .true. ; trunc_any = .true. - endif - enddo ; enddo - endif + do i=is,ie ; vel_report(i,J) = 3.0e8*US%m_s_to_L_T ; enddo ! Speed of light default. + do k=1,nz ; do i=is,ie + if (abs(v(i,J,k)) < CS%vel_underflow) v(i,J,k) = 0.0 + if (v(i,J,k) < 0.0) then + CFL = (-v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1)) + else + CFL = (v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j)) + endif + if (CFL > CS%CFL_trunc) trunc_any = .true. + if (CFL > CS%CFL_report) then + dowrite(i,J) = .true. + vel_report(i,J) = MIN(vel_report(i,J), abs(v(i,J,k))) + endif + enddo ; enddo do i=is,ie ; if (dowrite(i,J)) then v_old(i,J,:) = v(i,J,:) endif ; enddo - if (trunc_any) then ; if (CS%CFL_based_trunc) then + if (trunc_any) then do k=1,nz ; do i=is,ie if ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j+1) < -CS%CFL_trunc) then v(i,J,k) = (-0.9*CS%CFL_trunc) * (G%areaT(i,j+1) / (dt * G%dx_Cv(i,J))) @@ -3089,36 +3047,20 @@ subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS if (h(i,j,k) + h(i,j+1,k) > H_report) CS%ntrunc = CS%ntrunc + 1 endif enddo ; enddo - else - do k=1,nz ; do i=is,ie ; if (abs(v(i,J,k)) > maxvel) then - v(i,J,k) = SIGN(truncvel,v(i,J,k)) - if (h(i,j,k) + h(i,j+1,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - endif ; enddo ; enddo - endif ; endif + endif enddo ! J-loop else ! Do not report accelerations leading to large velocities. - if (CS%CFL_based_trunc) then - !$OMP parallel do default(shared) - do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - if (abs(v(i,J,k)) < CS%vel_underflow) then ; v(i,J,k) = 0.0 - elseif ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j+1) < -CS%CFL_trunc) then - v(i,J,k) = (-0.9*CS%CFL_trunc) * (G%areaT(i,j+1) / (dt * G%dx_Cv(i,J))) - if (h(i,j,k) + h(i,j+1,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - elseif ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j) > CS%CFL_trunc) then - v(i,J,k) = (0.9*CS%CFL_trunc) * (G%areaT(i,j) / (dt * G%dx_Cv(i,J))) - if (h(i,j,k) + h(i,j+1,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - endif - enddo ; enddo ; enddo - else - !$OMP parallel do default(shared) - do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - if (abs(v(i,J,k)) < CS%vel_underflow) then ; v(i,J,k) = 0.0 - elseif (abs(v(i,J,k)) > maxvel) then - v(i,J,k) = SIGN(truncvel, v(i,J,k)) - if (h(i,j,k) + h(i,j+1,k) > H_report) CS%ntrunc = CS%ntrunc + 1 - endif - enddo ; enddo ; enddo - endif + !$OMP parallel do default(shared) + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + if (abs(v(i,J,k)) < CS%vel_underflow) then ; v(i,J,k) = 0.0 + elseif ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j+1) < -CS%CFL_trunc) then + v(i,J,k) = (-0.9*CS%CFL_trunc) * (G%areaT(i,j+1) / (dt * G%dx_Cv(i,J))) + if (h(i,j,k) + h(i,j+1,k) > H_report) CS%ntrunc = CS%ntrunc + 1 + elseif ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j) > CS%CFL_trunc) then + v(i,J,k) = (0.9*CS%CFL_trunc) * (G%areaT(i,j) / (dt * G%dx_Cv(i,J))) + if (h(i,j,k) + h(i,j+1,k) > H_report) CS%ntrunc = CS%ntrunc + 1 + endif + enddo ; enddo ; enddo endif if (len_trim(CS%v_trunc_file) > 0) then @@ -3386,12 +3328,6 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & "near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//& "defined but LINEAR_DRAG is not.", & units="m", fail_if_missing=.true., scale=US%m_to_Z) - call get_param(param_file, mdl, "MAXVEL", CS%maxvel, & - "The maximum velocity allowed before the velocity components are truncated.", & - units="m s-1", default=3.0e8, scale=US%m_s_to_L_T) - call get_param(param_file, mdl, "CFL_BASED_TRUNCATIONS", CS%CFL_based_trunc, & - "If true, base truncations on the CFL number, and not an absolute speed.", & - default=.true.) call get_param(param_file, mdl, "CFL_TRUNCATE", CS%CFL_trunc, & "The value of the CFL number that will cause velocity "//& "components to be truncated; instability can occur past 0.5.", & @@ -3710,9 +3646,9 @@ end subroutine vertvisc_end !! side. Both of these thickness estimates are second order !! accurate. Above this the arithmetic mean thickness is used. !! -!! In addition, vertvisc truncates any velocity component that -!! exceeds maxvel to truncvel. This basically keeps instabilities -!! spatially localized. The number of times the velocity is +!! In addition, vertvisc truncates any velocity component that exceeds a +!! maximum CFL number to a fraction of this value. This basically keeps +!! instabilities spatially localized. The number of times the velocity is !! truncated is reported each time the energies are saved, and if !! exceeds CS%Maxtrunc the model will stop itself and change the time !! to a large value. This has proven very useful in (1) diagnosing From c5bb2bcf5f52fc48b07433a49902609d40d1c661 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 Aug 2025 17:00:46 -0400 Subject: [PATCH 4/4] Set intent(out) stub MOM_generic_tracer arguments Added code to set all intent(out) arguments in the stub version of the MOM_generic_tracer routines MOM_generic_tracer_stock(), MOM_generic_tracer_min_max() and MOM_generic_tracer_get() to avoid compile-time warnings. These stub routines are never supposed to be used, and there would be a fatal error during the registration phase if they are used. The values that are set in these stub routines are extremely large or otherwise nonsensical so that it would make it easier to debug them if they were used accidentally. All answers are bitwise identical in any cases that run successfully. --- .../GFDL_ocean_BGC/MOM_generic_tracer.F90 | 26 ++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90 b/config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90 index 80d209fa6e..808cace1e2 100644 --- a/config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90 +++ b/config_src/external/GFDL_ocean_BGC/MOM_generic_tracer.F90 @@ -16,7 +16,7 @@ module MOM_generic_tracer ! ### These imports should not reach into FMS directly ### use MOM_ALE_sponge, only : ALE_sponge_CS -use MOM_coms, only : EFP_type +use MOM_coms, only : EFP_type, real_to_EFP use MOM_diag_mediator, only : diag_ctrl use MOM_error_handler, only : MOM_error, FATAL use MOM_file_parser, only : param_file_type @@ -195,8 +195,14 @@ function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_inde integer :: MOM_generic_tracer_stock !< Return value, the !! number of stocks calculated here. + integer :: m MOM_generic_tracer_stock = 0 + ! These should never be used, but they are set to avoid compile-time warnings + do m=1,size(names) ; names(m) = "" ; enddo + do m=1,size(units) ; units(m) = "" ; enddo + do m=1,size(stocks) ; stocks(m) = real_to_EFP(0.0) ; enddo + end function MOM_generic_tracer_stock !> This subroutine finds the global min and max of either of all available @@ -227,8 +233,24 @@ function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, G, CS, na integer :: MOM_generic_tracer_min_max !< Return value, the !! number of tracers done here. + integer :: m + MOM_generic_tracer_min_max = 0 + ! These should never be used, but they are set to avoid compile-time warnings. Note that the minimum values + ! are delibarately set to be larger than the maximum values. + got_minmax(:) = .false. + gmax(:) = -huge(gmax) + gmin(:) = huge(gmin) + do m=1,size(names) ; names(m) = "" ; enddo + do m=1,size(units) ; units(m) = "" ; enddo + if (present(xgmin)) xgmin(:) = 0.0 + if (present(ygmin)) ygmin(:) = 0.0 + if (present(zgmin)) zgmin(:) = 0.0 + if (present(xgmax)) xgmax(:) = 0.0 + if (present(ygmax)) ygmax(:) = 0.0 + if (present(zgmax)) zgmax(:) = 0.0 + end function MOM_generic_tracer_min_max !> This subroutine calculates the surface state and sets coupler values for @@ -271,6 +293,8 @@ subroutine MOM_generic_tracer_get(name,member,array, CS) ! arbitrary units [A] character(len=128), parameter :: sub_name = 'MOM_generic_tracer_get' + array(:,:,:) = huge(array) + end subroutine MOM_generic_tracer_get !> This subroutine deallocates the memory owned by this module.