diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index 0d3da607ab..36ebf9c5b5 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -72,6 +72,8 @@ module MOM_hor_visc
   logical :: use_Leithy      !< If true, use a biharmonic form of 2D Leith
                              !! nonlinear eddy viscosity with harmonic backscatter.
                              !! Ah is the background. Leithy = Leith+E
+  logical :: nonsym_Leithy   !< If true, use sums that do not obey rotational symmetry in the
+                             !! smoothing that occurs when LEITHY = True.
   real    :: c_K             !< Fraction of energy dissipated by the biharmonic term
                              !! that gets backscattered in the Leith+E scheme. [nondim]
   logical :: smooth_Ah       !< If true (default), then Ah and m_leithy are smoothed.
@@ -297,6 +299,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     dudx, dvdy, &    ! components in the horizontal tension [T-1 ~> s-1]
     dudx_smooth, dvdy_smooth, & ! components in the horizontal tension from smoothed velocity [T-1 ~> s-1]
     GME_effic_h, &  ! The filtered efficiency of the GME terms at h points [nondim]
+    m_leithy, &   ! Kh=m_leithy*Ah in Leith+E parameterization [L-2 ~> m-2]
+    Ah_sq, &      ! The square of the biharmonic viscosity [L8 T-2 ~> m8 s-2]
     htot          ! The total thickness of all layers [Z ~> m]
   real :: Del2vort_h ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1]
   real :: grad_vel_mag_bt_h ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2]
@@ -322,8 +326,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     grad_div_mag_q, &  ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1]
     hq, &          ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2]
                    ! This form guarantees that hq/hu < 4.
-    GME_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim]
-    m_leithy       ! Kh=m_leithy*Ah in Leith+E parameterization [L-2 ~> m-2]
+    GME_effic_q    ! The filtered efficiency of the GME terms at q points [nondim]
   real :: grad_vel_mag_bt_q ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2]
   real :: boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim]
 
@@ -402,6 +405,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   logical :: apply_OBC = .false.
   logical :: use_MEKE_Ku
   logical :: use_MEKE_Au
+  integer :: is_vort, ie_vort, js_vort, je_vort
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
   integer :: i, j, k, n
   real :: inv_PI3, inv_PI2, inv_PI6 ! Powers of the inverse of pi [nondim]
@@ -465,6 +469,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
                           "RES_SCALE_MEKE_VISC is True.")
   endif
 
+  ! Set the halo sizes used for the vorticity calculations.
+  if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
+    js_vort = js-3 ; je_vort = Jeq+2 ; is_vort = is-3 ; ie_vort = Ieq+2
+    if ((G%isc-G%isd < 3) .or. (G%isc-G%isd < 3)) call MOM_error(FATAL, &
+          "The minimum halo size is 3 when a Leith viscosity is being used.")
+  else
+    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)
 
@@ -483,7 +496,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     call pass_var(h, G%domain, halo=2)
 
     ! Calculate the barotropic horizontal tension
-    do J=js-2,je+2 ; do I=is-2,ie+2
+    do j=js-2,je+2 ; do i=is-2,ie+2
       dudx_bt(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * ubtav(I,j) - &
                                      G%IdyCu(I-1,j) * ubtav(I-1,j))
       dvdy_bt(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * vbtav(i,J) - &
@@ -502,11 +515,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     enddo ; enddo
 
     if (CS%no_slip) then
-      do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2
+      do J=js-2,je+1 ; do I=is-2,ie+1
         sh_xy_bt(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
       enddo ; enddo
     else
-      do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2
+      do J=js-2,je+1 ; do I=is-2,ie+1
         sh_xy_bt(I,J) = G%mask2dBu(I,J) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
       enddo ; enddo
     endif
@@ -564,7 +577,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       ! One call applies the filter twice
       u_smooth(:,:,k) = u(:,:,k)
       v_smooth(:,:,k) = v(:,:,k)
-      call smooth_x9(G, field_u=u_smooth(:,:,k), field_v=v_smooth(:,:,k), zero_land=.false.)
+      call smooth_x9(G, field_u=u_smooth(:,:,k), field_v=v_smooth(:,:,k), zero_land=.false., &
+                     nonsym_sums=CS%nonsym_Leithy)
     enddo
     call pass_vector(u_smooth, v_smooth, G%Domain)
   endif
@@ -572,7 +586,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   !$OMP parallel do default(none) &
   !$OMP shared( &
   !$OMP   CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, &
-  !$OMP   is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &
+  !$OMP   is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, is_vort, ie_vort, js_vort, je_vort, &
   !$OMP   apply_OBC, rescale_Kh, legacy_bound, find_FrictWork, &
   !$OMP   use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, &
   !$OMP   backscat_subround, GME_effic_h, GME_effic_q, &
@@ -598,7 +612,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   !$OMP   dudx_smooth, dudy_smooth, dvdx_smooth, dvdy_smooth, &
   !$OMP   vort_xy_smooth, vort_xy_dx_smooth, vort_xy_dy_smooth, &
   !$OMP   sh_xx_smooth, sh_xy_smooth, &
-  !$OMP   vert_vort_mag_smooth, m_leithy, AhLthy &
+  !$OMP   vert_vort_mag_smooth, m_leithy, Ah_sq, AhLthy &
   !$OMP )
   do k=1,nz
 
@@ -622,18 +636,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     enddo ; enddo
 
     ! Components for the shearing strain
-    do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2
+    do J=js_vort,je_vort ; do I=is_vort,ie_vort
       dvdx(I,J) = CS%DY_dxBu(I,J)*(v(i+1,J,k)*G%IdyCv(i+1,J) - v(i,J,k)*G%IdyCv(i,J))
       dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k)*G%IdxCu(I,j+1) - u(I,j,k)*G%IdxCu(I,j))
     enddo ; enddo
-    if ((CS%Leith_Ah) .or. (CS%use_Leithy)) then
-      call pass_var(dvdx, G%Domain, halo=2, position=CORNER, complete=.false.)
-      call pass_var(dudy, G%Domain, halo=2, position=CORNER, complete=.true.)
-    endif
 
     if (CS%use_Leithy) then
       ! Calculate horizontal tension from smoothed velocity
-      do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         dudx_smooth(i,j) = CS%DY_dxT(i,j)*(G%IdyCu(I,j) * u_smooth(I,j,k) - &
                                            G%IdyCu(I-1,j) * u_smooth(I-1,j,k))
         dvdy_smooth(i,j) = CS%DX_dyT(i,j)*(G%IdxCv(i,J) * v_smooth(i,J,k) - &
@@ -642,18 +652,16 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       enddo ; enddo
 
       ! Components for the shearing strain from smoothed velocity
-      do J=Jsq,Jeq ; do I=Isq,Ieq
+      do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
         dvdx_smooth(I,J) = CS%DY_dxBu(I,J) * &
                          (v_smooth(i+1,J,k)*G%IdyCv(i+1,J) - v_smooth(i,J,k)*G%IdyCv(i,J))
         dudy_smooth(I,J) = CS%DX_dyBu(I,J) * &
                          (u_smooth(I,j+1,k)*G%IdxCu(I,j+1) - u_smooth(I,j,k)*G%IdxCu(I,j))
       enddo ; enddo
-      call pass_var(dvdx_smooth, G%Domain, halo=2, position=CORNER, complete=.false.)
-      call pass_var(dudy_smooth, G%Domain, halo=2, position=CORNER, complete=.true.)
     endif ! use Leith+E
 
     if (CS%id_normstress > 0) then
-      do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+      do j=js,je ; do i=is,ie
         NoSt(i,j,k) = sh_xx(i,j)
       enddo ; enddo
     endif
@@ -664,17 +672,17 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     ! 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
-      do j=js-2,je+2 ; do I=Isq-1,Ieq+1
+      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
-      do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2
+      do J=js-2,Jeq+1 ; do i=is-2,ie+2
         h_v(i,J) = 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i,j+1)*h(i,j+1,k))
       enddo ; enddo
     else
-      do j=js-2,je+2 ; do I=Isq-1,Ieq+1
+      do j=js-2,je+2 ; do I=is-2,Ieq+1
         h_u(I,j) = 0.5 * (h(i,j,k) + h(i+1,j,k))
       enddo ; enddo
-      do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2
+      do J=js-2,Jeq+1 ; do i=is-2,ie+2
         h_v(i,J) = 0.5 * (h(i,j,k) + h(i,j+1,k))
       enddo ; enddo
     endif
@@ -744,25 +752,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
         ! OBC projections, but they might not be necessary if the accelerations
         ! are always zeroed out at OBC points, in which case the i-loop below
         ! becomes do i=is-1,ie+1. -RWH
-        if ((J >= Jsq-1) .and. (J <= Jeq+1)) then
+        if ((J >= js-2) .and. (J <= Jeq+1)) then
           do i = max(is-2,OBC%segment(n)%HI%isd), min(ie+2,OBC%segment(n)%HI%ied)
             h_v(i,J) = h(i,j,k)
           enddo
         endif
       elseif (OBC%segment(n)%direction == OBC_DIRECTION_S) then
-        if ((J >= Jsq-1) .and. (J <= Jeq+1)) then
+        if ((J >= js-2) .and. (J <= Jeq+1)) then
           do i = max(is-2,OBC%segment(n)%HI%isd), min(ie+2,OBC%segment(n)%HI%ied)
             h_v(i,J) = h(i,j+1,k)
           enddo
         endif
       elseif (OBC%segment(n)%direction == OBC_DIRECTION_E) then
-        if ((I >= Isq-1) .and. (I <= Ieq+1)) then
+        if ((I >= is-2) .and. (I <= Ieq+1)) then
           do j = max(js-2,OBC%segment(n)%HI%jsd), min(je+2,OBC%segment(n)%HI%jed)
             h_u(I,j) = h(i,j,k)
           enddo
         endif
       elseif (OBC%segment(n)%direction == OBC_DIRECTION_W) then
-        if ((I >= Isq-1) .and. (I <= Ieq+1)) then
+        if ((I >= is-2) .and. (I <= Ieq+1)) then
           do j = max(js-2,OBC%segment(n)%HI%jsd), min(je+2,OBC%segment(n)%HI%jed)
             h_u(I,j) = h(i+1,j,k)
           enddo
@@ -774,25 +782,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       J = OBC%segment(n)%HI%JsdB ; I = OBC%segment(n)%HI%IsdB
       if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
         if ((J >= js-2) .and. (J <= je)) then
-          do I = max(Isq-1,OBC%segment(n)%HI%IsdB), min(Ieq+1,OBC%segment(n)%HI%IedB)
+          do I = max(is-2,OBC%segment(n)%HI%IsdB), min(Ieq+1,OBC%segment(n)%HI%IedB)
             h_u(I,j+1) = h_u(I,j)
           enddo
         endif
       elseif (OBC%segment(n)%direction == OBC_DIRECTION_S) then
         if ((J >= js-1) .and. (J <= je+1)) then
-          do I = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+1,OBC%segment(n)%HI%ied)
+          do I = max(is-2,OBC%segment(n)%HI%isd), min(Ieq+1,OBC%segment(n)%HI%ied)
             h_u(I,j) = h_u(I,j+1)
           enddo
         endif
       elseif (OBC%segment(n)%direction == OBC_DIRECTION_E) then
         if ((I >= is-2) .and. (I <= ie)) then
-          do J = max(Jsq-1,OBC%segment(n)%HI%jsd), min(Jeq+1,OBC%segment(n)%HI%jed)
+          do J = max(js-2,OBC%segment(n)%HI%jsd), min(Jeq+1,OBC%segment(n)%HI%jed)
             h_v(i+1,J) = h_v(i,J)
           enddo
         endif
       elseif (OBC%segment(n)%direction == OBC_DIRECTION_W) then
         if ((I >= is-1) .and. (I <= ie+1)) then
-          do J = max(Jsq-1,OBC%segment(n)%HI%jsd), min(Jeq+1,OBC%segment(n)%HI%jed)
+          do J = max(js-2,OBC%segment(n)%HI%jsd), min(Jeq+1,OBC%segment(n)%HI%jed)
             h_v(i,J) = h_v(i+1,J)
           enddo
         endif
@@ -817,11 +825,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       ! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask).
       ! dudy_smooth and dvdx_smooth do not (yet) include modifications at OBCs from above.
       if (CS%no_slip) then
-        do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+        do J=js-1,Jeq ; do I=is-1,Ieq
           sh_xy_smooth(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_smooth(I,J) + dudy_smooth(I,J) )
         enddo ; enddo
       else
-        do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
+        do J=js-1,Jeq ; do I=is-1,Ieq
           sh_xy_smooth(I,J) = G%mask2dBu(I,J) * ( dvdx_smooth(I,J) + dudy_smooth(I,J) )
         enddo ; enddo
       endif
@@ -854,23 +862,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     endif
 
     ! Vorticity
-    if (CS%no_slip) then
-      do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2
-        vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) )
-      enddo ; enddo
-    else
-      do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2
-        vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) )
-      enddo ; enddo
+    if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy) .or. (CS%id_vort_xy_q>0)) then
+      if (CS%no_slip) then
+        do J=js_vort,je_vort ; do I=is_vort,ie_vort
+          vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) )
+        enddo ; enddo
+      else
+        do J=js_vort,je_vort ; do I=is_vort,ie_vort
+          vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) )
+        enddo ; enddo
+      endif
     endif
 
     if (CS%use_Leithy) then
       if (CS%no_slip) then
-        do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1
+        do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
           vort_xy_smooth(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_smooth(I,J) - dudy_smooth(I,J) )
         enddo ; enddo
       else
-        do J=Jsq-2,Jeq+1 ; do I=Isq-2,Ieq+1
+        do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
           vort_xy_smooth(I,J) = G%mask2dBu(I,J) * ( dvdx_smooth(I,J) - dudy_smooth(I,J) )
         enddo ; enddo
       endif
@@ -884,25 +894,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
 
       ! Vorticity gradient
-      do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2
+      do J=js-2,Jeq+1 ; do i=is-2,Ieq+2
         DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
         vort_xy_dx(i,J) = DY_dxBu * (vort_xy(I,J) * G%IdyCu(I,j) - vort_xy(I-1,J) * G%IdyCu(I-1,j))
       enddo ; enddo
 
-      do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+1
+      do j=js-2,Jeq+2 ; do I=is-2,Ieq+1
         DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
         vort_xy_dy(I,j) = DX_dyBu * (vort_xy(I,J) * G%IdxCv(i,J) - vort_xy(I,J-1) * G%IdxCv(i,J-1))
       enddo ; enddo
 
       if (CS%use_Leithy) then
         ! Gradient of smoothed vorticity
-        do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2
+        do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
           DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
           vort_xy_dx_smooth(i,J) = DY_dxBu * &
                       (vort_xy_smooth(I,J) * G%IdyCu(I,j) - vort_xy_smooth(I-1,J) * G%IdyCu(I-1,j))
         enddo ; enddo
 
-        do j=Jsq-1,Jeq+2 ; do I=Isq-1,Ieq+1
+        do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
           DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
           vort_xy_dy_smooth(I,j) = DX_dyBu * &
                       (vort_xy_smooth(I,J) * G%IdxCv(i,J) - vort_xy_smooth(I,J-1) * G%IdxCv(i,J-1))
@@ -910,7 +920,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       endif ! If Leithy
 
       ! Laplacian of vorticity
-      do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1
+      do J=js-2,Jeq+1 ; do I=is-2,Ieq+1
         DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
         DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
 
@@ -929,24 +939,24 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
         enddo ; enddo
 
         ! Magnitude of divergence gradient
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           grad_div_mag_h(i,j) = sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2 + &
                                      (0.5*(div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2)
         enddo ; enddo
-        do j=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1
+        do J=js-1,Jeq ; do I=is-1,Ieq
           grad_div_mag_q(I,J) = sqrt((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2 + &
                                      (0.5*(div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2)
         enddo ; enddo
 
       else
 
-        do j=Jsq-1,Jeq+2 ; do I=is-2,Ieq+1
+        do j=js-2,Jeq+2 ; do I=is-2,Ieq+1
           div_xx_dx(I,j) = 0.0
         enddo ; enddo
-        do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+2
+        do J=js-2,Jeq+1 ; do i=is-2,Ieq+2
           div_xx_dy(i,J) = 0.0
         enddo ; enddo
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           grad_div_mag_h(i,j) = 0.0
         enddo ; enddo
         do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1
@@ -967,7 +977,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
       if (CS%use_QG_Leith_visc) then
 
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           grad_vort_mag_h_2d(i,j) = SQRT((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i,J-1)))**2 + &
                                          (0.5*(vort_xy_dy(I,j) + vort_xy_dy(I-1,j)))**2 )
         enddo ; enddo
@@ -982,7 +992,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
       endif
 
-      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
         grad_vort_mag_h(i,j) = SQRT((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i,J-1)))**2 + &
                                     (0.5*(vort_xy_dy(I,j) + vort_xy_dy(I-1,j)))**2 )
       enddo ; enddo
@@ -992,7 +1002,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       enddo ; enddo
 
       if (CS%use_Leithy) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           vert_vort_mag_smooth(i,j) = SQRT((0.5*(vort_xy_dx_smooth(i,J) + &
                                                  vort_xy_dx_smooth(i,J-1)))**2 + &
                                            (0.5*(vort_xy_dy_smooth(I,j) + &
@@ -1003,7 +1013,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     endif ! CS%Leith_Kh
 
     if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then
-      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
         sh_xx_sq = sh_xx(i,j)**2
         sh_xy_sq = 0.25 * ( (sh_xy(I-1,J-1)**2 + sh_xy(I,J)**2) &
                           + (sh_xy(I-1,J)**2 + sh_xy(I,J-1)**2) )
@@ -1012,13 +1022,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
     endif
 
     if (CS%better_bound_Ah .or. CS%better_bound_Kh) then
-      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
         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))
       enddo ; enddo
 
       if (CS%better_bound_Kh) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           visc_bound_rem(i,j) = 1.0
         enddo ; enddo
       endif
@@ -1031,26 +1041,26 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
       if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
         if (CS%use_QG_Leith_visc) then
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             grad_vort = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)
             grad_vort_qg = 3. * grad_vort_mag_h_2d(i,j)
             vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg)
           enddo ; enddo
         else
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             vert_vort_mag(i,j) = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)
           enddo ; enddo
         endif
       endif
 
       ! Static (pre-computed) background viscosity
-      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
         Kh(i,j) = CS%Kh_bg_xx(i,j)
       enddo ; enddo
 
       ! NOTE: The following do-block can be decomposed and vectorized after the
       !   stack size has been reduced.
-      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
         if (CS%add_LES_viscosity) then
           if (CS%Smagorinsky_Kh) &
             Kh(i,j) = Kh(i,j) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j)
@@ -1067,38 +1077,38 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       ! All viscosity contributions above are subject to resolution scaling
 
       if (rescale_Kh) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           Kh(i,j) = VarMix%Res_fn_h(i,j) * Kh(i,j)
         enddo ; enddo
       endif
 
       if (legacy_bound) then
         ! Older method of bounding for stability
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           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=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
         Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min)
       enddo ; enddo
 
       if (use_MEKE_Ku) then
         ! *Add* the MEKE contribution (which might be negative)
         if (CS%res_scale_MEKE) then
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j)
           enddo ; enddo
         else
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j)
           enddo ; enddo
         endif
       endif
 
       if (CS%anisotropic) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           ! *Add* the tension component of anisotropic viscosity
           Kh(i,j) = Kh(i,j) + CS%Kh_aniso * (1. - CS%n1n2_h(i,j)**2)
         enddo ; enddo
@@ -1106,7 +1116,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
       ! Newer method of bounding for stability
       if (CS%better_bound_Kh) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           if (Kh(i,j) >= hrat_min(i,j) * CS%Kh_Max_xx(i,j)) then
             visc_bound_rem(i,j) = 0.0
             Kh(i,j) = hrat_min(i,j) * CS%Kh_Max_xx(i,j)
@@ -1119,19 +1129,19 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       ! In Leith+E parameterization Kh is computed after Ah in the biharmonic loop.
       ! The harmonic component of str_xx is added in the biharmonic loop.
       if (CS%use_Leithy) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           Kh(i,j) = 0.
         enddo ; enddo
-      end if
+      endif
 
       if (CS%id_Kh_h>0 .or. CS%debug) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           Kh_h(i,j,k) = Kh(i,j)
         enddo ; enddo
       endif
 
       if (CS%id_grid_Re_Kh>0) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js,je ; do i=is,ie
           KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2)
           grid_Kh = max(Kh(i,j), CS%min_grid_Kh)
           grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j))) / grid_Kh
@@ -1139,13 +1149,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       endif
 
       if (CS%id_div_xx_h>0) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js,je ; do i=is,ie
           div_xx_h(i,j,k) = div_xx(i,j)
         enddo ; enddo
       endif
 
       if (CS%id_sh_xx_h>0) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js,je ; do i=is,ie
           sh_xx_h(i,j,k) = sh_xx(i,j)
         enddo ; enddo
       endif
@@ -1172,21 +1182,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       ! Determine the biharmonic viscosity at h points, using the
       ! largest value from several parameterizations. Also get the
       ! biharmonic component of str_xx.
-      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
         Ah(i,j) = CS%Ah_bg_xx(i,j)
       enddo ; enddo
 
       if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
         if (CS%Smagorinsky_Ah) then
           if (CS%bound_Coriolis) then
-            do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+           do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
               AhSm = Shear_mag(i,j) * (CS%Biharm_const_xx(i,j) &
                   + CS%Biharm_const2_xx(i,j) * Shear_mag(i,j) &
               )
               Ah(i,j) = max(Ah(i,j), AhSm)
             enddo ; enddo
           else
-            do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+            do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
               AhSm = CS%Biharm_const_xx(i,j) * Shear_mag(i,j)
               Ah(i,j) = max(Ah(i,j), AhSm)
             enddo ; enddo
@@ -1194,7 +1204,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
         endif
 
         if (CS%Leith_Ah) then
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + &
                                  (Del2vort_q(I-1,J) + Del2vort_q(I,J-1)))
             AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h) * inv_PI6
@@ -1204,7 +1214,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
         if (CS%use_Leithy) then
           ! Get m_leithy
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          if (CS%smooth_Ah) m_leithy(:,:) = 0.0 ! This is here to initialize domain edge halo values.
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + &
                                  (Del2vort_q(I-1,J) + Del2vort_q(I,J-1)))
             AhLth  = CS%Biharm6_const_xx(i,j) * inv_PI6 * abs(Del2vort_h)
@@ -1220,13 +1231,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
           enddo ; enddo
 
           if (CS%smooth_Ah) then
-            ! Smooth m_leithy. A single call smoothes twice.
-            call pass_var(m_leithy, G%Domain, halo=2, position=CORNER)
-            call smooth_x9(G, field_q=m_leithy)
-            call pass_var(m_leithy, G%Domain, position=CORNER)
+            ! Smooth m_leithy.  A single call smoothes twice.
+            call pass_var(m_leithy, G%Domain, halo=2)
+            call smooth_x9(G, field_h=m_leithy, nonsym_sums=CS%nonsym_Leithy)
+            call pass_var(m_leithy, G%Domain)
           endif
           ! Get Ah
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + &
                                  (Del2vort_q(I-1,J) + Del2vort_q(I,J-1)))
             AhLthy = CS%Biharm6_const_xx(i,j) * inv_PI6 * &
@@ -1234,28 +1245,28 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
             Ah(i,j) = max(CS%Ah_bg_xx(i,j), AhLthy)
           enddo ; enddo
           if (CS%smooth_Ah) then
-            ! Smooth Ah before applying upper bound
-            ! square, then smooth, then square root
+            ! Smooth Ah before applying upper bound square, then smooth, then square root.
             ! A single call smoothes twice.
-            do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-              Ah(i,j) = Ah(i,j)**2
+            Ah_sq(:,:) = 0.0 ! This is here to initialize domain edge halo values.
+            do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
+              Ah_sq(i,j) = Ah(i,j)**2
             enddo ; enddo
-            call pass_var(Ah, G%Domain, halo=2, position=CORNER)
-            call smooth_x9(G, field_q=Ah, zero_land=.false.)
-            call pass_var(Ah, G%Domain, position=CORNER)
-            do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-              Ah_h(i,j,k) = max(CS%Ah_bg_xx(i,j), sqrt(max(0., Ah(i,j))))
+            call pass_var(Ah_sq, G%Domain, halo=2)
+            call smooth_x9(G, field_h=Ah_sq, zero_land=.false., nonsym_sums=CS%nonsym_Leithy)
+            call pass_var(Ah_sq, G%Domain)
+            do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
+              Ah_h(i,j,k) = max(CS%Ah_bg_xx(i,j), sqrt(max(0., Ah_sq(i,j))))
               Ah(i,j)     = Ah_h(i,j,k)
             enddo ; enddo
           else
-            do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+            do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
               Ah_h(i,j,k) = Ah(i,j)
             enddo ; enddo
           endif
         endif
 
         if (CS%bound_Ah .and. .not. CS%better_bound_Ah) then
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xx(i,j))
           enddo ; enddo
         endif
@@ -1263,13 +1274,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
       if (use_MEKE_Au) then
         ! *Add* the MEKE contribution
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           Ah(i,j) = Ah(i,j) + MEKE%Au(i,j)
         enddo ; enddo
       endif
 
       if (CS%Re_Ah > 0.0) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           KE = 0.125*((u(I,j,k)+u(I-1,j,k))**2 + (v(i,J,k)+v(i,J-1,k))**2)
           Ah(i,j) = sqrt(KE) * CS%Re_Ah_const_xx(i,j)
         enddo ; enddo
@@ -1277,18 +1288,18 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
       if (CS%better_bound_Ah) then
         if (CS%better_bound_Kh) then
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             Ah(i,j) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * CS%Ah_Max_xx(i,j))
           enddo ; enddo
         else
-          do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+          do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
             Ah(i,j) = min(Ah(i,j), hrat_min(i,j) * CS%Ah_Max_xx(i,j))
           enddo ; enddo
         endif
       endif
 
       if ((CS%id_Ah_h>0) .or. CS%debug .or. CS%use_Leithy) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
           Ah_h(i,j,k) = Ah(i,j)
         enddo ; enddo
       endif
@@ -1296,14 +1307,14 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       if (CS%use_Leithy) then
         ! Compute Leith+E Kh after bounds have been applied to Ah
         ! and after it has been smoothed. Kh = -m_leithy * Ah
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-            Kh(i,j) = -m_leithy(i,j) * Ah(i,j)
-            Kh_h(i,j,k) = Kh(i,j)
+        do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
+          Kh(i,j) = -m_leithy(i,j) * Ah(i,j)
+          Kh_h(i,j,k) = Kh(i,j)
         enddo ; enddo
       endif
 
       if (CS%id_grid_Re_Ah>0) then
-        do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+        do j=js,je ; do i=is,ie
           KE = 0.125 * ((u(I,j,k) + u(I-1,j,k))**2 + (v(i,J,k) + v(i,J-1,k))**2)
           grid_Ah = max(Ah(i,j), CS%min_grid_Ah)
           grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j)) / grid_Ah
@@ -1497,7 +1508,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
         ! Leith+E doesn't recompute Kh at q points, it just interpolates it from h to q points
         if (CS%use_Leithy) then
-          Kh(I,J) = Kh_h(i+1,j+1,k)
+          Kh(I,J) = 0.25 * ((Kh_h(i,j,k) + Kh_h(i+1,j+1,k)) + (Kh_h(i,j+1,k) + Kh_h(i+1,j,k)))
         end if
 
         if (CS%id_Kh_q>0 .or. CS%debug) &
@@ -1604,7 +1615,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
       ! Leith+E doesn't recompute Ah at q points, it just interpolates it from h to q points
       if (CS%use_Leithy) then
         do J=js-1,Jeq ; do I=is-1,Ieq
-           Ah(I,J) = Ah_h(i+1,j+1,k)
+          Ah(I,J) = 0.25 * ((Ah_h(i,j,k) + Ah_h(i+1,j+1,k)) + (Ah_h(i,j+1,k) + Ah_h(i+1,j,k)))
         enddo ; enddo
       end if
 
@@ -1668,7 +1679,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
     else ! .not. use_GME
       ! This changes the units of str_xx from [L2 T-2 ~> m2 s-2] to [H L2 T-2 ~> m3 s-2 or kg s-2].
-      do J=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+      do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j))
       enddo ; enddo
 
@@ -2240,7 +2251,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
   if (.not.CS%Laplacian) CS%use_Kh_bg_2d = .false.
   call get_param(param_file, mdl, "KH_BG_2D_BUG", CS%Kh_bg_2d_bug, &
                  "If true, retain an answer-changing horizontal indexing bug in setting "//&
-                 "the corner-point viscosities when USE_KH_BG_2D=True. This is"//&
+                 "the corner-point viscosities when USE_KH_BG_2D=True.  This is "//&
                  "not recommended.", default=.false., do_not_log=.not.CS%use_Kh_bg_2d)
 
   call get_param(param_file, mdl, "USE_GME", CS%use_GME, &
@@ -2250,17 +2261,21 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
                  "Use the split time stepping if true.", default=.true., do_not_log=.true.)
   if (CS%use_Leithy) then
     if (.not.(CS%biharmonic .and. CS%Laplacian)) then
-                   call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init:"//&
+                   call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init: "//&
                    "LAPLACIAN and BIHARMONIC must both be True when USE_LEITHY=True.")
     endif
-    call get_param(param_file, mdl, "LEITHY_CK", CS%c_K, &
-                   "Fraction of biharmonic dissipation that gets backscattered, "//&
-                   "in Leith+E.", units="nondim", default=1.0)
-    call get_param(param_file, mdl, "SMOOTH_AH", CS%smooth_Ah, &
-                   "If true, Ah and m_leithy are smoothed within Leith+E. This requires"//&
-                   "lots of blocking communications, which can be expensive", &
-                   default=.true., do_not_log=.not.CS%use_Leithy)
   endif
+  call get_param(param_file, mdl, "LEITHY_CK", CS%c_K, &
+                 "Fraction of biharmonic dissipation that gets backscattered, "//&
+                 "in Leith+E.", units="nondim", default=1.0, do_not_log=.not.CS%use_Leithy)
+  call get_param(param_file, mdl, "SMOOTH_AH", CS%smooth_Ah, &
+                 "If true, Ah and m_leithy are smoothed within Leith+E.  This requires "//&
+                 "lots of blocking communications, which can be expensive", &
+                 default=.true., do_not_log=.not.CS%use_Leithy)
+  call get_param(param_file, mdl, "LEITHY_NONSYMMETRIC_SUMS", CS%nonsym_Leithy, &
+                 "If true, use sums that do not obey rotational symmetry in the smoothing that "//&
+                 "occurs when LEITHY = True.", default=.true., do_not_log=.not.CS%use_Leithy)
+  !### The default for LEITHY_NONSYMMETRIC_SUMS should be changed to .false.
 
   if (CS%use_GME .and. .not.split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// &
                                            "cannot be used with SPLIT=False.")
@@ -2397,7 +2412,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
     CS%dx2q(I,J) = G%dxBu(I,J)*G%dxBu(I,J) ; CS%dy2q(I,J) = G%dyBu(I,J)*G%dyBu(I,J)
     CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
   enddo ; enddo
-  do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
+  do j=js-2,Jeq+2 ; do i=is-2,Ieq+2
     CS%dx2h(i,j) = G%dxT(i,j)*G%dxT(i,j) ; CS%dy2h(i,j) = G%dyT(i,j)*G%dyT(i,j)
     CS%DX_dyT(i,j) = G%dxT(i,j)*G%IdyT(i,j) ; CS%DY_dxT(i,j) = G%dyT(i,j)*G%IdxT(i,j)
   enddo ; enddo
@@ -2438,7 +2453,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
     ! Calculate and store the background viscosity at h-points
 
     min_grid_sp_h2 = huge(1.)
-    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+    do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
       ! Static factors in the Smagorinsky and Leith schemes
       grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j) + CS%dy2h(i,j))
       CS%grid_sp_h2(i,j) = grid_sp_h2
@@ -2497,11 +2512,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
     enddo ; enddo
   endif
   if (CS%biharmonic) then
-    do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1
+    do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
       CS%Idx2dyCu(I,j) = (G%IdxCu(I,j)*G%IdxCu(I,j)) * G%IdyCu(I,j)
       CS%Idxdy2u(I,j) = G%IdxCu(I,j) * (G%IdyCu(I,j)*G%IdyCu(I,j))
     enddo ; enddo
-    do J=Jsq-1,Jeq+1 ; do i=is-1,Ieq+1
+    do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
       CS%Idx2dyCv(i,J) = (G%IdxCv(i,J)*G%IdxCv(i,J)) * G%IdyCv(i,J)
       CS%Idxdy2v(i,J) = G%IdxCv(i,J) * (G%IdyCv(i,J)*G%IdyCv(i,J))
     enddo ; enddo
@@ -2513,7 +2528,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
       BoundCorConst = 1.0 / (5.0*(bound_Cor_vel*bound_Cor_vel))
 
     min_grid_sp_h4 = huge(1.)
-    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+    do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
       grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j)+CS%dy2h(i,j))
       grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
       CS%grid_sp_h3(i,j) = grid_sp_h3
@@ -2571,7 +2586,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
   endif
   ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1.
   if (CS%Laplacian .and. CS%better_bound_Kh) then
-    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+    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)) * &
           max(G%IdyCu(I,j)*G%IareaCu(I,j), G%IdyCu(I-1,j)*G%IareaCu(I-1,j)) ), &
@@ -2599,7 +2614,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
   ! 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
-    do j=js-1,Jeq+1 ; do I=Isq-1,Ieq+1
+    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)) ) + &
                  CS%Idx2dyCu(I,j)*(CS%dx2q(I,J) * CS%DX_dyBu(I,J) * (G%IdxCu(I,j+1) + G%IdxCu(I,j)) + &
@@ -2609,7 +2624,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
                  CS%Idx2dyCu(I,j)*(CS%dx2q(I,J) * CS%DY_dxBu(I,J) * (G%IdyCv(i+1,J) + G%IdyCv(i,J))   + &
                                    CS%dx2q(I,J-1)*CS%DY_dxBu(I,J-1)*(G%IdyCv(i+1,J-1) + G%IdyCv(i,J-1)) ) )
     enddo ; enddo
-    do J=Jsq-1,Jeq+1 ; do i=is-1,Ieq+1
+    do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
       v0u(i,J) = (CS%Idxdy2v(i,J)*(CS%dy2q(I,J) * CS%DX_dyBu(I,J) * (G%IdxCu(I,j+1) + G%IdxCu(I,j))       + &
                                    CS%dy2q(I-1,J)*CS%DX_dyBu(I-1,J)*(G%IdxCu(I-1,j+1) + G%IdxCu(I-1,j)) ) + &
                  CS%Idx2dyCv(i,J)*(CS%dx2h(i,j+1)*CS%DY_dxT(i,j+1)*(G%IdyCu(I,j+1) + G%IdyCu(I-1,j+1))   + &
@@ -2619,7 +2634,7 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
                  CS%Idx2dyCv(i,J)*(CS%dx2h(i,j+1)*CS%DX_dyT(i,j+1)*(G%IdxCv(i,J+1) + G%IdxCv(i,J))   + &
                                    CS%dx2h(i,j) * CS%DX_dyT(i,j) * (G%IdxCv(i,J) + G%IdxCv(i,J-1)) ) )
     enddo ; enddo
-    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
+    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)*u0u(I,j) + G%IdyCu(I-1,j)*u0u(I-1,j))  + &
@@ -2904,28 +2919,32 @@ end subroutine smooth_GME
 !! horizontal_viscosity subroutine because they are not supposed to be halo-updated.
 !! But you _can_ apply them to Kh_h and Ah_h. Also note that it assumes indices
 !! is-2:ie+2, js-2:je+2 are correct on input.
-subroutine smooth_x9(G, field_h, field_u, field_v, field_q, zero_land)
+subroutine smooth_x9(G, field_h, field_u, field_v, field_q, zero_land, nonsym_sums)
   type(ocean_grid_type),                        intent(in)    :: G         !< Ocean grid
   real, dimension(SZI_(G),SZJ_(G)), optional,   intent(inout) :: field_h   !< field to be smoothed
-                                                              !! at h points
+                                                                           !! at h points [arbitrary]
   real, dimension(SZIB_(G),SZJ_(G)), optional,  intent(inout) :: field_u   !< field to be smoothed
-                                                              !! at u points
+                                                                           !! at u points [arbitrary]
   real, dimension(SZI_(G),SZJB_(G)), optional,  intent(inout) :: field_v   !< field to be smoothed
-                                                              !! at v points
+                                                                           !! at v points [arbitrary]
   real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: field_q   !< field to be smoothed
-                                                              !! at q points
+                                                                           !! at q points [arbitrary]
   logical, optional, intent(in)                               :: zero_land !< An optional argument
-                                                              !! indicating whether to set values
-                                                              !! on land to zero (not present or
-                                                              !! .true.) or whether to ignore land
-                                                              !!  values (.false.)
-  ! local variables. It would be good to make the _original variables allocatable.
-  real, dimension(SZI_(G),SZJ_(G))   :: field_h_original
-  real, dimension(SZIB_(G),SZJ_(G))  :: field_u_original
-  real, dimension(SZI_(G),SZJB_(G))  :: field_v_original
-  real, dimension(SZIB_(G),SZJB_(G)) :: field_q_original
-  real, dimension(3,3) :: weights, local_weights ! averaging weights for smoothing, nondimensional
+                                                                           !! indicating whether to set values
+                                                                           !! on land to zero (not present or
+                                                                           !! .true.) or whether to ignore land
+                                                                           !! values (.false.)
+  logical, optional, intent(in)                               :: nonsym_sums !< If present and true, use sums
+                                                                           !! that do not respect rotational symmetry.
+  ! Local variables.  It would be good to make the _prev variables allocatable.
+  real, dimension(SZI_(G),SZJ_(G))   :: fh_prev  ! The value of the h-point field at the previous iteration [arbitrary]
+  real, dimension(SZIB_(G),SZJ_(G))  :: fu_prev  ! The value of the u-point field at the previous iteration [arbitrary]
+  real, dimension(SZI_(G),SZJB_(G))  :: fv_prev  ! The value of the v-point field at the previous iteration [arbitrary]
+  real, dimension(SZIB_(G),SZJB_(G)) :: fq_prev  ! The value of the q-point field at the previous iteration [arbitrary]
+  real, dimension(3,3) :: weights, local_weights ! averaging weights for smoothing [nondim]
+  real :: Iwts ! The inverse of the sum of the weights [nondim]
   logical :: zero_land_val ! actual value of zero_land optional argument
+  logical :: broken_symmetry ! If true, use expressions that do not obey rotational symmetry.
   integer :: i, j, s
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq
 
@@ -2934,6 +2953,7 @@ subroutine smooth_x9(G, field_h, field_u, field_v, field_q, zero_land)
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
 
   weights = reshape([1., 2., 1., 2., 4., 2., 1., 2., 1.],shape(weights))/16.
+  broken_symmetry = .false. ; if (present(nonsym_sums)) broken_symmetry = nonsym_sums
 
   if (present(zero_land)) then
     zero_land_val = zero_land
@@ -2942,59 +2962,127 @@ subroutine smooth_x9(G, field_h, field_u, field_v, field_q, zero_land)
   endif
 
   if (present(field_h)) then
-    do s=1,0,-1
-      field_h_original(:,:) = field_h(:,:)
-      ! apply smoothing on field_h
+    if (broken_symmetry) then ; do s=1,0,-1
+      fh_prev(:,:) = field_h(:,:)
+      ! apply smoothing on field_h using the original non-rotationally symmetric expressions.
       do j=js-s,je+s ; do i=is-s,ie+s
         ! skip land points
         if (G%mask2dT(i,j)==0.) cycle
         ! compute local weights
         local_weights = weights*G%mask2dT(i-1:i+1,j-1:j+1)
         if (.not. zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16)
-        field_h(i,j) =  sum(local_weights*field_h_original(i-1:i+1,j-1:j+1))
+        field_h(i,j) =  sum(local_weights*fh_prev(i-1:i+1,j-1:j+1))
       enddo ; enddo
-    enddo
+    enddo ; else ; do s=1,0,-1
+      fh_prev(:,:) = field_h(:,:)
+      ! apply smoothing on field_h using rotationally symmetric expressions.
+      do j=js-s,je+s ; do i=is-s,ie+s ; if (G%mask2dT(i,j) > 0.0) then
+        Iwts = 0.0625
+        if (.not. zero_land_val) &
+          Iwts = 1.0 / ( (4.0*G%mask2dT(i,j) + &
+                          ( 2.0*((G%mask2dT(i-1,j) + G%mask2dT(i+1,j)) + &
+                                 (G%mask2dT(i,j-1) + G%mask2dT(i,j+1))) + &
+                           ((G%mask2dT(i-1,j-1) + G%mask2dT(i+1,j+1)) + &
+                            (G%mask2dT(i-1,j+1) + G%mask2dT(i+1,j-1))) ) ) + 1.0e-16 )
+        field_h(i,j) = Iwts * ( 4.0*G%mask2dT(i,j) * fh_prev(i,j) &
+                              + (2.0*((G%mask2dT(i-1,j) * fh_prev(i-1,j) + G%mask2dT(i+1,j) * fh_prev(i+1,j)) + &
+                                      (G%mask2dT(i,j-1) * fh_prev(i,j-1) + G%mask2dT(i,j+1) * fh_prev(i,j+1))) &
+                                + ((G%mask2dT(i-1,j-1) * fh_prev(i-1,j-1) + G%mask2dT(i+1,j+1) * fh_prev(i+1,j+1)) + &
+                                   (G%mask2dT(i-1,j+1) * fh_prev(i-1,j+1) + G%mask2dT(i+1,j-1) * fh_prev(i-1,j-1))) ))
+      endif ; enddo ; enddo
+    enddo ; endif
   endif
 
-  if (present(field_u)) then
-    do s=1,0,-1
-      field_u_original(:,:) = field_u(:,:)
-      ! apply smoothing on field_u
+  if (present(field_u) .and. present(field_v)) then
+    if (broken_symmetry) then ; do s=1,0,-1
+      fu_prev(:,:) = field_u(:,:)
+      ! apply smoothing on field_u using the original non-rotationally symmetric expressions.
       do j=js-s,je+s ; do I=Isq-s,Ieq+s
         ! skip land points
         if (G%mask2dCu(I,j)==0.) cycle
         ! compute local weights
         local_weights = weights*G%mask2dCu(I-1:I+1,j-1:j+1)
         if (.not. zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16)
-        field_u(I,j) =  sum(local_weights*field_u_original(I-1:I+1,j-1:j+1))
+        field_u(I,j) =  sum(local_weights*fu_prev(I-1:I+1,j-1:j+1))
       enddo ; enddo
 
-      field_v_original(:,:) = field_v(:,:)
-      ! apply smoothing on field_v
+      fv_prev(:,:) = field_v(:,:)
+      ! apply smoothing on field_v using the original non-rotationally symmetric expressions.
       do J=Jsq-s,Jeq+s ; do i=is-s,ie+s
         ! skip land points
         if (G%mask2dCv(i,J)==0.) cycle
         ! compute local weights
         local_weights = weights*G%mask2dCv(i-1:i+1,J-1:J+1)
         if (.not. zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16)
-        field_v(i,J) =  sum(local_weights*field_v_original(i-1:i+1,J-1:J+1))
+        field_v(i,J) =  sum(local_weights*fv_prev(i-1:i+1,J-1:J+1))
       enddo ; enddo
-    enddo
+    enddo ; else ; do s=1,0,-1
+      fu_prev(:,:) = field_u(:,:)
+      ! apply smoothing on field_u using the original non-rotationally symmetric expressions.
+      do j=js-s,je+s ; do I=Isq-s,Ieq+s ; if (G%mask2dCu(I,j) > 0.0) then
+        Iwts = 0.0625
+        if (.not. zero_land_val) &
+          Iwts = 1.0 / ( (4.0*G%mask2dCu(I,j) + &
+                          ( 2.0*((G%mask2dCu(I-1,j) + G%mask2dCu(I+1,j)) + &
+                                 (G%mask2dCu(I,j-1) + G%mask2dCu(I,j+1))) + &
+                           ((G%mask2dCu(I-1,j-1) + G%mask2dCu(I+1,j+1)) + &
+                            (G%mask2dCu(I-1,j+1) + G%mask2dCu(I+1,j-1))) ) ) + 1.0e-16 )
+        field_u(I,j) = Iwts * ( 4.0*G%mask2dCu(I,j) * fu_prev(I,j) &
+                              + (2.0*((G%mask2dCu(I-1,j) * fu_prev(I-1,j) + G%mask2dCu(I+1,j) * fu_prev(I+1,j)) + &
+                                      (G%mask2dCu(I,j-1) * fu_prev(I,j-1) + G%mask2dCu(I,j+1) * fu_prev(I,j+1))) &
+                                + ((G%mask2dCu(I-1,j-1) * fu_prev(I-1,j-1) + G%mask2dCu(I+1,j+1) * fu_prev(I+1,j+1)) + &
+                                   (G%mask2dCu(I-1,j+1) * fu_prev(I-1,j+1) + G%mask2dCu(I+1,j-1) * fu_prev(I-1,j-1))) ))
+      endif ; enddo ; enddo
+
+      fv_prev(:,:) = field_v(:,:)
+      ! apply smoothing on field_v using the original non-rotationally symmetric expressions.
+      do J=Jsq-s,Jeq+s ; do i=is-s,ie+s ; if (G%mask2dCv(i,J) > 0.0) then
+        Iwts = 0.0625
+        if (.not. zero_land_val) &
+          Iwts = 1.0 / ( (4.0*G%mask2dCv(i,J) + &
+                          ( 2.0*((G%mask2dCv(i-1,J) + G%mask2dCv(i+1,J)) + &
+                                 (G%mask2dCv(i,J-1) + G%mask2dCv(i,J+1))) + &
+                           ((G%mask2dCv(i-1,J-1) + G%mask2dCv(i+1,J+1)) + &
+                            (G%mask2dCv(i-1,J+1) + G%mask2dCv(i+1,J-1))) ) ) + 1.0e-16 )
+        field_v(i,J) = Iwts * ( 4.0*G%mask2dCv(i,J) * fv_prev(i,J) &
+                              + (2.0*((G%mask2dCv(i-1,J) * fv_prev(i-1,J) + G%mask2dCv(i+1,J) * fv_prev(i+1,J)) + &
+                                      (G%mask2dCv(i,J-1) * fv_prev(i,J-1) + G%mask2dCv(i,J+1) * fv_prev(i,J+1))) &
+                                + ((G%mask2dCv(i-1,J-1) * fv_prev(i-1,J-1) + G%mask2dCv(i+1,J+1) * fv_prev(i+1,J+1)) + &
+                                   (G%mask2dCv(i-1,J+1) * fv_prev(i-1,J+1) + G%mask2dCv(i+1,J-1) * fv_prev(i-1,J-1))) ))
+      endif ; enddo ; enddo
+    enddo ; endif
   endif
 
   if (present(field_q)) then
-    do s=1,0,-1
-      field_q_original(:,:) = field_q(:,:)
-      ! apply smoothing on field_q
+    if (broken_symmetry) then ; do s=1,0,-1
+      fq_prev(:,:) = field_q(:,:)
+      ! apply smoothing on field_q using the original non-rotationally symmetric expressions.
       do J=Jsq-s,Jeq+s ; do I=Isq-s,Ieq+s
         ! skip land points
         if (G%mask2dBu(I,J)==0.) cycle
         ! compute local weights
         local_weights = weights*G%mask2dBu(I-1:I+1,J-1:J+1)
         if (.not. zero_land_val) local_weights = local_weights/(sum(local_weights) + 1.E-16)
-        field_q(I,J) =  sum(local_weights*field_q_original(I-1:I+1,J-1:J+1))
+        field_q(I,J) =  sum(local_weights*fq_prev(I-1:I+1,J-1:J+1))
       enddo ; enddo
-    enddo
+    enddo ; else ; do s=1,0,-1
+      fq_prev(:,:) = field_q(:,:)
+      ! apply smoothing on field_q using rotationally symmetric expressions.
+      do J=Jsq-s,Jeq+s ; do I=Isq-s,Ieq+s ; if (G%mask2dBu(I,J) > 0.0) then
+        Iwts = 0.0625
+        if (.not. zero_land_val) &
+          Iwts = 1.0 / ( (4.0*G%mask2dBu(I,J) + &
+                          ( 2.0*((G%mask2dBu(I-1,J) + G%mask2dBu(I+1,J)) + &
+                                 (G%mask2dBu(I,J-1) + G%mask2dBu(I,J+1))) + &
+                           ((G%mask2dBu(I-1,J-1) + G%mask2dBu(I+1,J+1)) + &
+                            (G%mask2dBu(I-1,J+1) + G%mask2dBu(I+1,J-1))) ) ) + 1.0e-16 )
+        field_q(I,J) = Iwts * ( 4.0*G%mask2dBu(I,J) * fq_prev(I,J) &
+                              + (2.0*((G%mask2dBu(I-1,J) * fq_prev(I-1,J) + G%mask2dBu(I+1,J) * fq_prev(I+1,J)) + &
+                                      (G%mask2dBu(I,J-1) * fq_prev(I,J-1) + G%mask2dBu(I,J+1) * fq_prev(I,J+1))) &
+                                + ((G%mask2dBu(I-1,J-1) * fq_prev(I-1,J-1) + G%mask2dBu(I+1,J+1) * fq_prev(I+1,J+1)) + &
+                                   (G%mask2dBu(I-1,J+1) * fq_prev(I-1,J+1) + G%mask2dBu(I+1,J-1) * fq_prev(I-1,J-1))) ))
+      endif ; enddo ; enddo
+    enddo ; endif
   endif
 
 end subroutine smooth_x9