Fix USE_CONT_THICKNESS bug

Optionally correct a bug due to a missing halo exchange.  This likely
isn't needed when compiled for nonsymmetric memory, but the added halo
exchange does no harm in such cases.  The pass_vector call could probably
be replaced with fill_symmetric_edges, except there is no subroutine
fill_vector_symmetric_edges_3d.

USE_CONT_THICKNESS is not yet widely used, so rather than preserving the old
(incorrect) solutions by default, this bug is corrected by default.  However,
the previous answers can be recovered by setting the new runtime parameter
USE_CONT_THICKNESS_BUG to true.  This parameter is only used (and logged)
when USE_CONT_THICKNESS set to true.

By default, this commit does change answers in symmetric memory cases with
USE_CONT_THICKNESS = True, and there is a new runtime parameter in such cases.

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -123,6 +123,7 @@ module MOM_hor_visc
   real    :: min_grid_Ah     !< Minimun horizontal biharmonic viscosity used to
                              !! limit grid Reynolds number [L4 T-1 ~> m4 s-1]
   logical :: use_cont_thick  !< If true, thickness at velocity points adopts h[uv] in BT_cont from continuity solver.
+  logical :: use_cont_thick_bug  !< If true, retain an answer-changing bug for thickness at velocity points.
   type(ZB2020_CS) :: ZB2020  !< Zanna-Bolton 2020 control structure.
   logical :: use_ZB2020      !< If true, use Zanna-Bolton 2020 parameterization.
 
@@ -282,9 +283,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
   type(thickness_diffuse_CS), optional, intent(in) :: TD !< Thickness diffusion control structure
   type(accel_diag_ptrs), optional, intent(in)  :: ADp  !< Acceleration diagnostics
   real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
-                        optional, intent(in) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2].
+                    optional, intent(inout) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2].
   real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
-                        optional, intent(in) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2].
+                    optional, intent(inout) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2].
 
   ! Local variables
   real, dimension(SZIB_(G),SZJ_(G)) :: &
@@ -477,6 +478,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
   use_MEKE_Au = allocated(MEKE%Au)
 
   use_cont_huv = CS%use_cont_thick .and. present(hu_cont) .and. present(hv_cont)
+  if (use_cont_huv .and. .not.CS%use_cont_thick_bug) then
+    call pass_vector(hu_cont, hv_cont, G%domain, To_All+Scalar_Pair, halo=2)
+  endif
 
   rescale_Kh = .false.
   if (VarMix%use_variable_mixing) then
@@ -709,7 +713,14 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
     ! in OBCs, which are not ordinarily be necessary, and might not be necessary
     ! even with OBCs if the accelerations are zeroed at OBC points, in which
     ! case the j-loop for h_u could collapse to j=js=1,je+1. -RWH
-    if (CS%use_land_mask) then
+    if (use_cont_huv) then
+      do j=js-2,je+2 ; do I=Isq-1,Ieq+1
+        h_u(I,j) = hu_cont(I,j,k)
+      enddo ; enddo
+      do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2
+        h_v(i,J) = hv_cont(i,J,k)
+      enddo ; enddo
+    elseif (CS%use_land_mask) then
       do j=js-2,je+2 ; do I=is-2,Ieq+1
         h_u(I,j) = 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i+1,j)*h(i+1,j,k))
       enddo ; enddo
@@ -725,16 +736,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
       enddo ; enddo
     endif
 
-    ! The following should obviously be combined with the previous block if adopted.
-    if (use_cont_huv) then
-      do j=js-2,je+2 ; do I=Isq-1,Ieq+1
-        h_u(I,j) = hu_cont(I,j,k)
-      enddo ; enddo
-      do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2
-        h_v(i,J) = hv_cont(i,J,k)
-      enddo ; enddo
-    endif
-
     ! Adjust contributions to shearing strain and interpolated values of
     ! thicknesses on open boundaries.
     if (apply_OBC) then ; do n=1,OBC%number_of_segments
@@ -2140,6 +2141,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
   call get_param(param_file, mdl, "USE_CONT_THICKNESS", CS%use_cont_thick, &
                  "If true, use thickness at velocity points from continuity solver. This option "//&
                  "currently only works with split mode.", default=.false.)
+  call get_param(param_file, mdl, "USE_CONT_THICKNESS_BUG", CS%use_cont_thick_bug, &
+                   "If true, retain an answer-changing halo update bug when "//&
+                   "USE_CONT_THICKNESS=True.  This is not recommended.", &
+                   default=.false., do_not_log=.not.CS%use_cont_thick)
+
   call get_param(param_file, mdl, "LAPLACIAN", CS%Laplacian, &
                  "If true, use a Laplacian horizontal viscosity.", &
                  default=.false.)