Use dt for vertvisc_remnant() in predictor

Add an option to use `dt` instead of `dt_pred` in the calculation of
vertvisc_remnant() at the end of predictor stage. The change affects
the following `continuity`() call and `btstep`() call in corrector step.

Combined with the changes from the previous commit, the new options
should solve the issue of inconsistent barotropic velocities from
barotropic solver and baroclinic step.

A runtime flag VISC_REM_TIMESTEP_FIX is added to control the behavior of
 this commit. The current default is false (not applying the fix).

src/core/MOM_barotropic.F90

@@ -4614,8 +4614,8 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
 
   call get_param(param_file, mdl, "BAROTROPIC_VERTICAL_WEIGHT_FIX", CS%wt_uv_fix, &
                  "If true, use a normalized weight function for vertical averages of "//&
-                 "baroclinic velocity and forcing.", &
-                 default=.false.)
+                 "baroclinic velocity and forcing. This flag should be used with "//&
+                 "VISC_REM_TIMESTEP_FIX.", default=.false.)
   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, &

src/core/MOM_dynamics_split_RK2.F90

@@ -177,6 +177,7 @@ module MOM_dynamics_split_RK2
   logical :: debug_OBC !< If true, do debugging calls for open boundary conditions.
   logical :: fpmix = .false.                 !< If true, applies profiles of momentum flux magnitude and direction.
   logical :: module_is_initialized = .false. !< Record whether this module has been initialized.
+  logical :: visc_rem_dt_fix = .false. !<If true, use dt rather than dt_pred for vertvisc_rem at the end of predictor.
 
   !>@{ Diagnostic IDs
   integer :: id_uold   = -1, id_vold   = -1
@@ -725,7 +726,11 @@ subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, f
     call start_group_pass(CS%pass_uvp, G%Domain, clock=id_clock_pass)
     call cpu_clock_begin(id_clock_vertvisc)
   endif
-  call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt_pred, G, GV, US, CS%vertvisc_CSp)
+  if (CS%visc_rem_dt_fix) then
+    call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
+  else
+    call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt_pred, G, GV, US, CS%vertvisc_CSp)
+  endif
   call cpu_clock_end(id_clock_vertvisc)
 
   call do_group_pass(CS%pass_visc_rem, G%Domain, clock=id_clock_pass)
@@ -1414,6 +1419,11 @@ 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_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 "//&
+                 "BAROTROPIC_VERTICAL_WEIGHT_FIX,", default=.false.)
 
   allocate(CS%taux_bot(IsdB:IedB,jsd:jed), source=0.0)
   allocate(CS%tauy_bot(isd:ied,JsdB:JedB), source=0.0)