+Rotationally symmetric epipycnal diffusion option

Hallberg-NOAA · marshallward · commit a3fd1f3a7466 · 2024-04-05T14:22:08.000-04:00
Added the option to do epipycnal tracer diffusion between a bulk mixed layer
and the interior ocean with expressions that satisfy rotational symmetry.  This
option is enabled by setting the new runtime parameter HOR_DIFF_ANSWER_DATE to
be greater than 20240330.  By default, this parameter is set to use the previous
expressions, although the default may be changed later to follow
DEFAULT_ANSWER_DATE.  Also corrected two bugs with the tracer limits used to
repartition the fluxes between layers within tracer_epipycnal_ML_diff; this
correction is enables by setting the new HOR_DIFF_LIMIT_BUG to .false., but to
retain previous answers by default it is set to true. By default all answers are
bitwise identical, but there are changes to some MOM_parameter_doc files due to
the introduction of two new runtime parameters.
diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90
@@ -65,6 +65,14 @@ module MOM_tracer_hor_diff
                                    !! tracer_hor_diff.
   logical :: recalc_neutral_surf   !< If true, recalculate the neutral surfaces if CFL has been
                                    !! exceeded
+  logical :: limit_bug             !< If true and the answer date is 20240330 or below, use a
+                                   !! rotational symmetry breaking bug when limiting the tracer
+                                   !! properties in tracer_epipycnal_ML_diff.
+  integer :: answer_date           !< The vintage of the order of arithmetic to use for the tracer
+                                   !! diffusion.  Values of 20240330 or below recover the answers
+                                   !! from the original form of this code, while higher values use
+                                   !! mathematically equivalent expressions that recover rotational symmetry
+                                   !! when DIFFUSE_ML_TO_INTERIOR is true.
   type(neutral_diffusion_CS), pointer :: neutral_diffusion_CSp => NULL() !< Control structure for neutral diffusion.
   type(hbd_CS), pointer    :: hor_bnd_diffusion_CSp => NULL() !< Control structure for
                                                               !! horizontal boundary diffusion.
@@ -678,7 +686,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
   real, dimension(SZI_(G),SZJB_(G)),        intent(in)    :: khdt_epi_y !< Meridional epipycnal diffusivity times
                                                            !! a time step and the ratio of the open face width over
                                                            !! the distance between adjacent tracer points [L2 ~> m2]
-  type(unit_scale_type),                    intent(in)    :: US !< A dimensional unit scaling type
+  type(unit_scale_type),                    intent(in)    :: US         !< A dimensional unit scaling type
   type(tracer_hor_diff_CS),                 intent(inout) :: CS         !< module control structure
   type(thermo_var_ptrs),                    intent(in)    :: tv         !< thermodynamic structure
   integer,                                  intent(in)    :: num_itts   !< number of iterations (usually=1)
@@ -706,13 +714,16 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
     k0b_Lv, k0a_Lv, &  ! The original k-indices of the layers that participate
     k0b_Rv, k0a_Rv     ! in each pair of mixing at v-faces.
 
-  !### Accumulating the converge into this array one face at a time may lead to a lack of rotational symmetry.
-  real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: &
-    tr_flux_conv  ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
+    tr_flux_N, &      ! The tracer flux through the northern face [conc H L2 ~> conc m3 or conc kg]
+    tr_flux_S, &      ! The tracer flux through the southern face [conc H L2 ~> conc m3 or conc kg]
+    tr_flux_E, &      ! The tracer flux through the eastern face [conc H L2 ~> conc m3 or conc kg]
+    tr_flux_W, &      ! The tracer flux through the western face [conc H L2 ~> conc m3 or conc kg]
+    tr_flux_conv      ! The flux convergence of tracers [conc H L2 ~> conc m3 or conc kg]
 
   ! The following 3-d arrays were created in 2014 in MOM6 PR#12 to facilitate openMP threading
-  ! on an i-loop, which might have been ill advised.  The k-size extents here might also be problematic.
-  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
+  ! on an i-loop, which might have been ill advised.
+  real, dimension(SZI_(G),SZJB_(G),SZK_(GV)*2) :: &
     Tr_flux_3d, &     ! The tracer flux through pairings at meridional faces [conc H L2 ~> conc m3 or conc kg]
     Tr_adj_vert_L, &  ! Vertical adjustments to which layer the fluxes go into in the southern
                       ! columns at meridional face [conc H L2 ~> conc m3 or conc kg]
@@ -815,6 +826,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
   do k=2,nkmb ; do j=js-2,je+2 ; do i=is-2,ie+2
     if (Rml_max(i,j) < rho_coord(i,j,k)) Rml_max(i,j) = rho_coord(i,j,k)
   enddo ; enddo ; enddo
+
   !   Use bracketing and bisection to find the k-level that the densest of the
   ! mixed and buffer layer corresponds to, such that:
   !     GV%Rlay(max_kRho-1) < Rml_max <= GV%Rlay(max_kRho)
@@ -1191,12 +1203,7 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
 
   endif ; enddo ; enddo ! i- & j- loops over meridional faces.
 
-! The tracer-specific calculations start here.
-
-  ! Zero out tracer tendencies.
-  do k=1,PEmax_kRho ; do j=js-1,je+1 ; do i=is-1,ie+1
-    tr_flux_conv(i,j,k) = 0.0
-  enddo ; enddo ; enddo
+  ! The tracer-specific calculations start here.
 
   do itt=1,max_itt
 
@@ -1205,12 +1212,19 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
     endif
 
     do m=1,ntr
-!$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPu,m,max_kRho,nz,h,h_exclude, &
-!$OMP                                  k0b_Lu,k0b_Ru,deep_wt_Lu,k0a_Lu,deep_wt_Ru,k0a_Ru,   &
-!$OMP                                  hP_Lu,hP_Ru,I_maxitt,khdt_epi_x,tr_flux_conv,Idt) &
-!$OMP                          private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, &
-!$OMP                                     Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, &
-!$OMP                                     Tr_flux,Tr_adj_vert,wt_a,vol)
+      ! Zero out tracer tendencies.
+      if (CS%answer_date <= 20240330) then
+        tr_flux_conv(:,:,:) = 0.0
+      else
+        tr_flux_N(:,:,:) = 0.0 ; tr_flux_S(:,:,:) = 0.0
+        tr_flux_E(:,:,:) = 0.0 ; tr_flux_W(:,:,:) = 0.0
+      endif
+      tr_flux_3d(:,:,:) = 0.0
+      tr_adj_vert_R(:,:,:) = 0.0 ; tr_adj_vert_L(:,:,:) = 0.0
+
+      !$OMP parallel do default(shared) private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,Tr_La, &
+      !$OMP                                     Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,h_L,h_R, &
+      !$OMP                                     Tr_flux,Tr_adj_vert,wt_a,vol)
       do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then
         ! Determine the fluxes through the zonal faces.
 
@@ -1230,7 +1244,11 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
           kRb = kRa ; if (max_kRho(i+1,j) < nz) kRb = max_kRho(i+1,j)+1
           Tr_La = Tr_min_face ; Tr_Lb = Tr_La ; Tr_Ra = Tr_La ; Tr_Rb = Tr_La
           if (h(i,j,kLa) > h_exclude) Tr_La = Tr(m)%t(i,j,kLa)
-          if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb)
+          if ((CS%answer_date <= 20240330) .and. CS%limit_bug) then
+            if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb)
+          else
+            if (h(i,j,kLb) > h_exclude) Tr_Lb = Tr(m)%t(i,j,kLb)
+          endif
           if (h(i+1,j,kRa) > h_exclude) Tr_Ra = Tr(m)%t(i+1,j,kRa)
           if (h(i+1,j,kRb) > h_exclude) Tr_Rb = Tr(m)%t(i+1,j,kRb)
           Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
@@ -1264,12 +1282,20 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
           endif
 
           h_L = hP_Lu(j)%p(I,k) ; h_R = hP_Ru(j)%p(I,k)
-          Tr_flux = I_maxitt * khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R) * &
-            ((2.0 * h_L * h_R) / (h_L + h_R))
-
+          if (CS%answer_date <= 20240330) then
+            Tr_flux = I_maxitt * khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R) * &
+                      ((2.0 * h_L * h_R) / (h_L + h_R))
+          else
+            Tr_flux = I_maxitt * ((2.0 * h_L * h_R) / (h_L + h_R)) * &
+                      khdt_epi_x(I,j) * (Tr_av_L - Tr_av_R)
+          endif
 
           if (deep_wt_Lu(j)%p(I,k) >= 1.0) then
-            tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux
+            if (CS%answer_date <= 20240330) then
+              tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux
+            else
+              tr_flux_E(i,j,kLb) = tr_flux_E(i,j,kLb) + Tr_flux
+            endif
           else
             Tr_adj_vert = 0.0
             wt_b = deep_wt_Lu(j)%p(I,k) ; wt_a = 1.0 - wt_b
@@ -1299,12 +1325,21 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
               endif
             endif
 
-            tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux + Tr_adj_vert)
-            tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux - Tr_adj_vert)
+            if (CS%answer_date <= 20240330) then
+              tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux + Tr_adj_vert)
+              tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux - Tr_adj_vert)
+            else
+              tr_flux_E(i,j,kLa) = tr_flux_E(i,j,kLa) + (wt_a*Tr_flux + Tr_adj_vert)
+              tr_flux_E(i,j,kLb) = tr_flux_E(i,j,kLb) + (wt_b*Tr_flux - Tr_adj_vert)
+            endif
           endif
 
           if (deep_wt_Ru(j)%p(I,k) >= 1.0) then
-            tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + Tr_flux
+            if (CS%answer_date <= 20240330) then
+              tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + Tr_flux
+            else
+              tr_flux_W(i+1,j,kRb) = tr_flux_W(i+1,j,kRb) + Tr_flux
+            endif
           else
             Tr_adj_vert = 0.0
             wt_b = deep_wt_Ru(j)%p(I,k) ; wt_a = 1.0 - wt_b
@@ -1334,23 +1369,22 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
               endif
             endif
 
-            tr_flux_conv(i+1,j,kRa) = tr_flux_conv(i+1,j,kRa) + &
-                                            (wt_a*Tr_flux - Tr_adj_vert)
-            tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + &
-                                            (wt_b*Tr_flux + Tr_adj_vert)
+            if (CS%answer_date <= 20240330) then
+              tr_flux_conv(i+1,j,kRa) = tr_flux_conv(i+1,j,kRa) + (wt_a*Tr_flux - Tr_adj_vert)
+              tr_flux_conv(i+1,j,kRb) = tr_flux_conv(i+1,j,kRb) + (wt_b*Tr_flux + Tr_adj_vert)
+            else
+              tr_flux_W(i+1,j,kRa) = tr_flux_W(i+1,j,kRa) + (wt_a*Tr_flux - Tr_adj_vert)
+              tr_flux_W(i+1,j,kRb) = tr_flux_W(i+1,j,kRb) + (wt_b*Tr_flux + Tr_adj_vert)
+            endif
           endif
           if (associated(Tr(m)%df2d_x)) &
             Tr(m)%df2d_x(I,j) = Tr(m)%df2d_x(I,j) + Tr_flux * Idt
         enddo ! Loop over pairings at faces.
       endif ; enddo ; enddo ! i- & j- loops over zonal faces.
 
-!$OMP parallel do default(none) shared(is,ie,js,je,G,Tr,nkmb,nPv,m,max_kRho,nz,h,h_exclude, &
-!$OMP                                  k0b_Lv,k0b_Rv,deep_wt_Lv,k0a_Lv,deep_wt_Rv,k0a_Rv,   &
-!$OMP                                  hP_Lv,hP_Rv,I_maxitt,khdt_epi_y,Tr_flux_3d,          &
-!$OMP                                  Tr_adj_vert_L,Tr_adj_vert_R,Idt)                     &
-!$OMP                          private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,             &
-!$OMP                                  Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,        &
-!$OMP                                  h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
+      !$OMP parallel do default(shared) private(Tr_min_face,Tr_max_face,kLa,kLb,kRa,kRb,             &
+      !$OMP                                  Tr_La,Tr_Lb,Tr_Ra,Tr_Rb,Tr_av_L,wt_b,Tr_av_R,        &
+      !$OMP                                  h_L,h_R,Tr_flux,Tr_adj_vert,wt_a,vol)
       do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then
         ! Determine the fluxes through the meridional faces.
 
@@ -1370,7 +1404,11 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
           kRb = kRa ; if (max_kRho(i,j+1) < nz) kRb = max_kRho(i,j+1)+1
           Tr_La = Tr_min_face ; Tr_Lb = Tr_La ; Tr_Ra = Tr_La ; Tr_Rb = Tr_La
           if (h(i,j,kLa) > h_exclude) Tr_La = Tr(m)%t(i,j,kLa)
-          if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb)
+          if ((CS%answer_date <= 20240330) .and. CS%limit_bug) then
+            if (h(i,j,kLb) > h_exclude) Tr_La = Tr(m)%t(i,j,kLb)
+          else
+            if (h(i,j,kLb) > h_exclude) Tr_Lb = Tr(m)%t(i,j,kLb)
+          endif
           if (h(i,j+1,kRa) > h_exclude) Tr_Ra = Tr(m)%t(i,j+1,kRa)
           if (h(i,j+1,kRb) > h_exclude) Tr_Rb = Tr(m)%t(i,j+1,kRb)
           Tr_min_face = min(Tr_min_face, Tr_La, Tr_Lb, Tr_Ra, Tr_Rb)
@@ -1464,42 +1502,69 @@ subroutine tracer_epipycnal_ML_diff(h, dt, Tr, ntr, khdt_epi_x, khdt_epi_y, G, &
             Tr(m)%df2d_y(i,J) = Tr(m)%df2d_y(i,J) + Tr_flux * Idt
         enddo ! Loop over pairings at faces.
       endif ; enddo ; enddo ! i- & j- loops over meridional faces.
-!$OMP parallel do default(none) shared(is,ie,js,je,G,nPv,k0b_Lv,k0b_Rv,deep_wt_Lv,  &
-!$OMP                                  tr_flux_conv,Tr_flux_3d,k0a_Lv,Tr_adj_vert_L,&
-!$OMP                                  deep_wt_Rv,k0a_Rv,Tr_adj_vert_R) &
-!$OMP                          private(kLa,kLb,kRa,kRb,wt_b,wt_a)
-      do i=is,ie ; do J=js-1,je ; if (G%mask2dCv(i,J) > 0.0) then
+
+      !$OMP parallel do default(shared) private(kLa,kLb,kRa,kRb,wt_b,wt_a)
+      do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then
         ! The non-stride-1 loop order here is to facilitate openMP threading. However, it might be
         ! suboptimal when openMP threading is not used, at which point it might be better to fuse
-        ! these loope with those that precede it and thereby eliminate the need for three 3-d arrays.
-        do k=1,nPv(i,J)
-          kLb = k0b_Lv(J)%p(i,k); kRb = k0b_Rv(J)%p(i,k)
-          if (deep_wt_Lv(J)%p(i,k) >= 1.0) then
-            tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux_3d(i,J,k)
-          else
-            kLa = k0a_Lv(J)%p(i,k)
-            wt_b = deep_wt_Lv(J)%p(i,k) ; wt_a = 1.0 - wt_b
-            tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k))
-            tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k))
-          endif
-          if (deep_wt_Rv(J)%p(i,k) >= 1.0) then
-            tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + Tr_flux_3d(i,J,k)
-          else
-            kRa = k0a_Rv(J)%p(i,k)
-            wt_b = deep_wt_Rv(J)%p(i,k) ; wt_a = 1.0 - wt_b
-            tr_flux_conv(i,j+1,kRa) = tr_flux_conv(i,j+1,kRa) + &
-                                            (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k))
-            tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + &
-                                            (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k))
-          endif
-        enddo
+        ! this loop with those that precede it and thereby eliminate the need for three 3-d arrays.
+        if (CS%answer_date <= 20240330) then
+          do k=1,nPv(i,J)
+            kLb = k0b_Lv(J)%p(i,k); kRb = k0b_Rv(J)%p(i,k)
+            if (deep_wt_Lv(J)%p(i,k) >= 1.0) then
+              tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - Tr_flux_3d(i,J,k)
+            else
+              kLa = k0a_Lv(J)%p(i,k)
+              wt_b = deep_wt_Lv(J)%p(i,k) ; wt_a = 1.0 - wt_b
+              tr_flux_conv(i,j,kLa) = tr_flux_conv(i,j,kLa) - (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k))
+              tr_flux_conv(i,j,kLb) = tr_flux_conv(i,j,kLb) - (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k))
+            endif
+            if (deep_wt_Rv(J)%p(i,k) >= 1.0) then
+              tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + Tr_flux_3d(i,J,k)
+            else
+              kRa = k0a_Rv(J)%p(i,k)
+              wt_b = deep_wt_Rv(J)%p(i,k) ; wt_a = 1.0 - wt_b
+              tr_flux_conv(i,j+1,kRa) = tr_flux_conv(i,j+1,kRa) + &
+                                              (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k))
+              tr_flux_conv(i,j+1,kRb) = tr_flux_conv(i,j+1,kRb) + &
+                                              (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k))
+            endif
+          enddo
+        else
+          do k=1,nPv(i,J)
+            kLb = k0b_Lv(J)%p(i,k); kRb = k0b_Rv(J)%p(i,k)
+            if (deep_wt_Lv(J)%p(i,k) >= 1.0) then
+              tr_flux_N(i,j,kLb) = tr_flux_N(i,j,kLb) + Tr_flux_3d(i,J,k)
+            else
+              kLa = k0a_Lv(J)%p(i,k)
+              wt_b = deep_wt_Lv(J)%p(i,k) ; wt_a = 1.0 - wt_b
+              tr_flux_N(i,j,kLa) = tr_flux_N(i,j,kLa) + (wt_a*Tr_flux_3d(i,J,k) + Tr_adj_vert_L(i,J,k))
+              tr_flux_N(i,j,kLb) = tr_flux_N(i,j,kLb) + (wt_b*Tr_flux_3d(i,J,k) - Tr_adj_vert_L(i,J,k))
+            endif
+            if (deep_wt_Rv(J)%p(i,k) >= 1.0) then
+              tr_flux_S(i,j+1,kRb) = tr_flux_S(i,j+1,kRb) + Tr_flux_3d(i,J,k)
+            else
+              kRa = k0a_Rv(J)%p(i,k)
+              wt_b = deep_wt_Rv(J)%p(i,k) ; wt_a = 1.0 - wt_b
+              tr_flux_S(i,j+1,kRa) = tr_flux_S(i,j+1,kRa) + (wt_a*Tr_flux_3d(i,J,k) - Tr_adj_vert_R(i,J,k))
+              tr_flux_S(i,j+1,kRb) = tr_flux_S(i,j+1,kRb) + (wt_b*Tr_flux_3d(i,J,k) + Tr_adj_vert_R(i,J,k))
+            endif
+          enddo
+        endif
       endif ; enddo ; enddo
+
+      if (CS%answer_date >= 20240331) then
+        !$OMP parallel do default(shared)
+        do k=1,PEmax_kRho ; do j=js,je ; do i=is,ie
+          tr_flux_conv(i,j,k) = ((tr_flux_W(i,j,k) - tr_flux_E(i,j,k)) + &
+                                 (tr_flux_S(i,j,k) - tr_flux_N(i,j,k)))
+        enddo ; enddo ; enddo
+      endif
+
       !$OMP parallel do default(shared)
       do k=1,PEmax_kRho ; do j=js,je ; do i=is,ie
         if ((G%mask2dT(i,j) > 0.0) .and. (h(i,j,k) > 0.0)) then
-          Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / &
-                                            (h(i,j,k)*G%areaT(i,j))
-          tr_flux_conv(i,j,k) = 0.0
+          Tr(m)%t(i,j,k) = Tr(m)%t(i,j,k) + tr_flux_conv(i,j,k) / (h(i,j,k)*G%areaT(i,j))
         endif
       enddo ; enddo ; enddo
 
@@ -1546,6 +1611,7 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic
   ! This include declares and sets the variable "version".
 # include "version_variable.h"
   character(len=40)  :: mdl = "MOM_tracer_hor_diff" ! This module's name.
+  integer :: default_answer_date
 
   if (associated(CS)) then
     call MOM_error(WARNING, "tracer_hor_diff_init called with associated control structure.")
@@ -1604,6 +1670,21 @@ subroutine tracer_hor_diff_init(Time, G, GV, US, param_file, diag, EOS, diabatic
                  "If true, then recalculate the neutral surfaces if the \n"//&
                  "diffusive CFL is exceeded. If false, assume that the  \n"//&
                  "positions of the surfaces do not change \n", default=.false.)
+  call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
+                 "This sets the default value for the various _ANSWER_DATE parameters.", &
+                 default=99991231, do_not_log=.true.)
+  call get_param(param_file, mdl, "HOR_DIFF_ANSWER_DATE", CS%answer_date, &
+                 "The vintage of the order of arithmetic to use for the tracer diffusion.  "//&
+                 "Values of 20240330 or below recover the answers from the original form of the "//&
+                 "along-isopycnal mixed layer to interior mixing code, while higher values use "//&
+                 "mathematically equivalent expressions that recover rotational symmetry "//&
+                 "when DIFFUSE_ML_TO_INTERIOR is true.", &
+                 default=20240101, do_not_log=.not.CS%Diffuse_ML_interior)
+                 !### Change the default later to default_answer_date.
+  call get_param(param_file, mdl, "HOR_DIFF_LIMIT_BUG", CS%limit_bug, &
+                 "If true and the answer date is 20240330 or below, use a rotational symmetry "//&
+                 "breaking bug when limiting the tracer properties in tracer_epipycnal_ML_diff.", &
+                 default=.true., do_not_log=((.not.CS%Diffuse_ML_interior).or.(CS%answer_date>=20240331)))
   CS%ML_KhTR_scale = 1.0
   if (CS%Diffuse_ML_interior) then
     call get_param(param_file, mdl, "ML_KHTR_SCALE", CS%ML_KhTR_scale, &
