diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index fb343978fc..2d463b909c 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1062,7 +1062,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, Iwt_u_tot(I,j) = Iwt_u_tot(I,j) + wt_u(I,j,k) enddo ; enddo ; enddo do j=js,je ; do I=is-1,ie - Iwt_u_tot(I,j) = G%mask2dCu(I,j) / (Iwt_u_tot(I,j) + h_neglect) + if (abs(Iwt_u_tot(I,j)) > 0.0 ) Iwt_u_tot(I,j) = G%mask2dCu(I,j) / Iwt_u_tot(I,j) enddo ; enddo do k=1,nz ; do j=js,je ; do I=is-1,ie wt_u(I,j,k) = wt_u(I,j,k) * Iwt_u_tot(I,j) @@ -1073,7 +1073,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, Iwt_v_tot(i,J) = Iwt_v_tot(i,J) + wt_v(i,J,k) enddo ; enddo ; enddo do J=js-1,je ; do i=is,ie - Iwt_v_tot(i,J) = G%mask2dCv(i,J) / (Iwt_v_tot(i,J) + h_neglect) + if (abs(Iwt_v_tot(i,J)) > 0.0 ) Iwt_v_tot(i,J) = G%mask2dCv(i,J) / Iwt_v_tot(i,J) enddo ; enddo do k=1,nz ; do J=js-1,je ; do i=is,ie wt_v(i,J,k) = wt_v(i,J,k) * Iwt_v_tot(i,J) @@ -4464,6 +4464,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: use_BT_cont_type logical :: use_tides + logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG. character(len=48) :: thickness_units, flux_units character*(40) :: hvel_str integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -4612,10 +4613,12 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701) + call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, default=.true., do_not_log=.true.) call get_param(param_file, mdl, "VISC_REM_BT_WEIGHT_FIX", CS%wt_uv_fix, & "If true, use a normalized weight function for vertical averages of "//& - "baroclinic velocity and forcing. This flag should be used with "//& - "VISC_REM_TIMESTEP_FIX.", default=.false.) + "baroclinic velocity and forcing. Default of this flag is set by "//& + "VISC_REM_BUG. This flag should be used with VISC_REM_TIMESTEP_FIX.", & + default=.not.visc_rem_bug) call get_param(param_file, mdl, "TIDES", use_tides, & "If true, apply tidal momentum forcing.", default=.false.) call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, & diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 31464c22f0..cedcdc573b 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -2713,6 +2713,7 @@ subroutine continuity_PPM_init(Time, G, GV, US, param_file, diag, CS) !> This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_continuity_PPM" ! This module's name. + logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG. CS%initialized = .true. @@ -2773,9 +2774,11 @@ subroutine continuity_PPM_init(Time, G, GV, US, param_file, diag, CS) "If true, use the marginal face areas from the continuity "//& "solver for use as the weights in the barotropic solver. "//& "Otherwise use the transport averaged areas.", default=.true.) + call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, default=.true., do_not_log=.true.) call get_param(param_file, mdl, "VISC_REM_CONT_HVEL_FIX", CS%visc_rem_hvel_fix, & "If true, velocity cell thickness h_[uv] from the continuity solver "//& - "is not multiplied by visc_rem_[uv].", default=.false.) + "is not multiplied by visc_rem_[uv]. Default of this flag is set by "//& + "VISC_REM_BUG.", default=.not.visc_rem_bug) CS%diag => diag id_clock_reconstruct = cpu_clock_id('(Ocean continuity reconstruction)', grain=CLOCK_ROUTINE) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 57e118b47f..daabc8d135 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1353,6 +1353,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p type(group_pass_type) :: pass_av_h_uvh logical :: debug_truncations logical :: read_uv, read_h2 + logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB @@ -1420,11 +1421,18 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p call get_param(param_file, mdl, "DEBUG_OBC", CS%debug_OBC, default=.false.) call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, & default=.false.) + call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, & + "If true, visc_rem_[uv] in split mode is incorrectly calculated or accounted "//& + "for in three places. This parameter controls the defaults of three individual "//& + "flags, VISC_REM_TIMESTEP_FIX in MOM_dynamics_split_RK2(b), "//& + "VISC_REM_BT_WEIGHT_FIX in MOM_barotropic, and VISC_REM_CONT_HVEL_FIX in "//& + "MOM_continuity_PPM. Eventually, the three individual flags should be removed "//& + "after tests and the default of VISC_REM_BUG should be to False.", default=.true.) call get_param(param_file, mdl, "VISC_REM_TIMESTEP_FIX", CS%visc_rem_dt_fix, & "If true, use dt rather than dt_pred in vertvisc_remnant() at the end of "//& "predictor stage for the following continuity() call and btstep() call "//& "in the corrector step. This flag should be used with "//& - "VISC_REM_BT_WEIGHT_FIX.", default=.false.) + "VISC_REM_BT_WEIGHT_FIX.", default=.not.visc_rem_bug) allocate(CS%taux_bot(IsdB:IedB,jsd:jed), source=0.0) allocate(CS%tauy_bot(isd:ied,JsdB:JedB), source=0.0) diff --git a/src/core/MOM_dynamics_split_RK2b.F90 b/src/core/MOM_dynamics_split_RK2b.F90 index cadd976455..e16c10cde6 100644 --- a/src/core/MOM_dynamics_split_RK2b.F90 +++ b/src/core/MOM_dynamics_split_RK2b.F90 @@ -1277,6 +1277,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, character(len=48) :: thickness_units, flux_units, eta_rest_name logical :: debug_truncations logical :: read_uv, read_h2 + logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB @@ -1335,11 +1336,18 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, call get_param(param_file, mdl, "DEBUG_OBC", CS%debug_OBC, default=.false.) call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, & default=.false.) + call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, & + "If true, visc_rem_[uv] in split mode is incorrectly calculated or accounted "//& + "for in three places. This parameter controls the defaults of three individual "//& + "flags, VISC_REM_TIMESTEP_FIX in MOM_dynamics_split_RK2(b), "//& + "VISC_REM_BT_WEIGHT_FIX in MOM_barotropic, and VISC_REM_CONT_HVEL_FIX in "//& + "MOM_continuity_PPM. Eventually, the three individual flags should be removed "//& + "after tests and the default of VISC_REM_BUG should be to False.", default=.true.) call get_param(param_file, mdl, "VISC_REM_TIMESTEP_FIX", CS%visc_rem_dt_fix, & "If true, use dt rather than dt_pred in vertvisc_remnant() at the end of "//& "predictor stage for the following continuity() call and btstep() call "//& - "in the corrector step. This flag should be used with "//& - "VISC_REM_BT_WEIGHT_FIX.", default=.false.) + "in the corrector step. Default of this flag is set by VISC_REM_BUG. "//& + "This flag should be used with VISC_REM_BT_WEIGHT_FIX.", default=.not.visc_rem_bug) allocate(CS%taux_bot(IsdB:IedB,jsd:jed), source=0.0) allocate(CS%tauy_bot(isd:ied,JsdB:JedB), source=0.0)