Use Adcroft_reciprocal for Iwt_[uv]_tot

* Correct a bug using dimensional h_neglect in calculating
non-dimensional Iwt_[uv]_tot, by replacing the division with
an Adcroft_reciprocal style division.

* Add a parameter VISC_REM_BUG to control the defaults of all three
viscous remnant bug related parameters. The individual parameters should
 be removed in the future.

Author: herrwang0

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, &

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)

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)

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)