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.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_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index fffe35104b..b727c30595 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) & @@ -2128,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/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/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..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,8 +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=.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%bound_Kh).and.(CS%bound_Ah))) call get_param(param_file, mdl, "SMAG_BI_CONST",Smag_bi_const, & "The nondimensional biharmonic Smagorinsky constant, "//& @@ -2644,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", & @@ -2670,7 +2627,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"//& @@ -2744,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 @@ -2802,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 @@ -2909,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) @@ -2944,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 @@ -2963,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)) @@ -2993,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) @@ -3018,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)) * & @@ -3053,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/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_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_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index c87a99ba44..caf0555d28 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -2468,14 +2468,12 @@ 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 "//& "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/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 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/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/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 "//& 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.", &