diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index 2091c51429e2..a21d1de90a15 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -326,6 +326,14 @@
 					description="If true, the N2 limiting is applied to the horizontal diffusion term"
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_Redi_use_quasi_monotone_limiter" type="logical" default_value=".false."
+					description="If true, fluxes are reduced to prevent tracers from violating monotonicity. Cross-term fluxes are scaled toward zero to prevent tracers from under/overshooting the min/max values in adjacent cells and layers"
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_Redi_quasi_monotone_safety_factor" type="real" default_value="1.0"
+					description="A safety factor applied to flux scaling when monotonicity is violated. Smaller values scale fluxes toward zero more aggressively."
+					possible_values="A value between 0 and 1"
+		/>
 		<nml_option name="config_Redi_min_layers_diag_terms" type="integer" default_value="6"
 					description="Redi diagonal terms (2 and 3) are turned off from layer 1 through config_Redi_min_layers_diag_terms-1, and on from config_Redi_min_layers_diag_terms to nVertLevels. The Redi diagonal terms are not guaranteed to produce bounded tracer fields, and in practice produce growing temperatures in a few columns with fewer than 5 vertical cells. Redi is meant for isopycnal mixing in the deep ocean, so not applying Redi diagonal terms in very shallow regions is an acceptable solution."
 					possible_values="any integer between 0 (all layers on) and nVertLevels (all layers off)"
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
index 83dd0986a239..ba4a50cf6438 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_hmix_redi.F
@@ -54,9 +54,13 @@ module ocn_tracer_hmix_Redi
    !--------------------------------------------------------------------
 
    real(kind=RKIND) :: term1TaperFactor
+   real(kind=RKIND) :: limiter_safety_factor
+   logical :: use_quasi_monotone
    real(kind=RKIND), dimension(:), allocatable :: rediKappaCell
    real(kind=RKIND), dimension(:, :, :), allocatable :: redi_term2_edge, redi_term3_topOfCell
    real(kind=RKIND), dimension(:, :, :), allocatable :: redi_term1, redi_term2, redi_term3
+   real(kind=RKIND), dimension(:, :, :), allocatable :: redi_flux_scale
+   real(kind=RKIND), dimension(:, :, :), allocatable :: minimumBnd, maximumBnd
 
 !***********************************************************************
 
@@ -127,18 +131,20 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
       !
       !-----------------------------------------------------------------
 
-      integer :: iCell, iEdge, cell1, cell2, iCellSelf
-      integer :: i, k, iTr, nTracers, nCells, nEdges, nCellsP1
+      integer :: iCell, jCell, iEdge, cell1, cell2, iCellSelf
+      integer :: i, k, kUp, kDn, iTr, nTracers, nCells, nEdges, nCellsP1
       logical :: use_Redi_diag_terms
       integer, pointer :: nVertLevels
       integer, dimension(:), pointer :: nCellsArray, nEdgesArray
 
-      integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, nEdgesOnCell, minLevelCell, maxLevelCell
+      integer, dimension(:), pointer :: minLevelEdgeBot, maxLevelEdgeTop, &
+                                        nEdgesOnCell, minLevelCell, maxLevelCell
       integer, dimension(:, :), pointer :: cellsOnEdge, edgesOnCell, edgeSignOnCell, cellsOnCell
 
       real(kind=RKIND) :: tempTracer, invAreaCell
       real(kind=RKIND) :: flux, flux_term2, flux_term3, dTracerDx, coef
       real(kind=RKIND) :: r_tmp, tracer_turb_flux, kappaRediEdgeVal
+      real(kind=RKIND) :: over, diff, scal, scalUp, scalDn
       real(kind=RKIND), dimension(:), pointer :: areaCell, dvEdge, dcEdge
       real(kind=RKIND), dimension(:), allocatable :: minimumVal
       real(kind=RKIND), dimension(:, :), allocatable :: fluxRediZTop
@@ -176,6 +182,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
       allocate (redi_term3(nTracers, nVertLevels, nCells))
       allocate (redi_term2_edge(nTracers, nVertLevels, nEdges))
       allocate (redi_term3_topOfCell(nTracers, nVertLevels + 1, nCells))
+      allocate (minimumBnd(nTracers, nVertLevels, nCells))
+      allocate (maximumBnd(nTracers, nVertLevels, nCells))
+      allocate (redi_flux_scale(nTracers, nVertLevels, nCells))
 
       ! RediKappa changes every time step if either:
       !   config_Redi_horizontal_taper == 'RossbyRadius'  
@@ -183,6 +192,18 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
       ! in that case, RediKappa is updated every time step in mpas_ocn_gm.F
       ! For static RediKappa, it is set on init and remains unchanged.
 
+      !$omp parallel
+      !$omp do schedule(runtime) private(k,iTr)
+      do iCell=1,nCells
+         do k=1,nVertLevels
+            do iTr = 1,nTracers
+               redi_flux_scale(iTr,k,iCell) = 1.0_RKIND
+            end do
+         end do
+      end do
+      !$omp end do
+      !$omp end parallel
+
       nCells = nCellsArray(2)
       !$omp parallel
       !$omp do schedule(runtime) private(i, iEdge)
@@ -200,13 +221,15 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
       nCells = nCellsArray(2)
       nCellsP1 = nCellsArray(size(nCellsArray)) + 1
 
-      ! Term 1: this is the "standard" horizontal del2 term, but with RediKappa coefficient.
-      ! \kappa_2 \nabla \phi on edge
+      ! Term 1: the "standard" horizontal del^2 term, but with RediKappa coefficient.
+      ! Term 2: a cross-flux: kappa div( h S d/dz psi )
+      ! Term 3: a cross-flux: kappa d/dz div( S grad psi )
+      ! Term 4: the vertical flux (done in vmix) d/dz kappa S^2 d/dz psi
       !$omp parallel
       !$omp do schedule(runtime) &
       !$omp private(fluxRediZTop, invAreaCell, i, iEdge, cell1, cell2, iCellSelf, r_tmp, coef, k, &
       !$omp         kappaRediEdgeVal, iTr, tracer_turb_flux, flux, flux_term2, flux_term3, &
-      !$omp         dTracerDx, use_Redi_diag_terms)
+      !$omp         dTracerDx, use_Redi_diag_terms, kUp, kDn)
       do iCell = 1, nCells
          if (maxLevelCell(iCell) .ge. config_Redi_min_layers_diag_terms) then
             use_Redi_diag_terms = .true.
@@ -214,6 +237,24 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
             use_Redi_diag_terms = .false.
          endif
 
+         minimumBnd(:, :, iCell) = 0.0_RKIND
+         maximumBnd(:, :, iCell) = 0.0_RKIND
+         if (use_quasi_monotone) then
+            ! init. min/max bnds in column
+            do k = minLevelCell(iCell), maxLevelCell(iCell)
+               kUp = max(minLevelCell(iCell), k - 1)
+               kDn = min(maxLevelCell(iCell), k + 1)
+               do iTr = 1, nTracers
+                  minimumBnd(iTr, k, iCell) = min(tracers(iTr, kUp, iCell), &
+                                                  tracers(iTr, k, iCell), &
+                                                  tracers(iTr, kDn, iCell))
+                  maximumBnd(iTr, k, iCell) = max(tracers(iTr, kUp, iCell), &
+                                                  tracers(iTr, k, iCell), &
+                                                  tracers(iTr, kDn, iCell))
+               end do
+            end do
+         end if
+
          redi_term1(:, :, iCell) = 0.0_RKIND
          redi_term2(:, :, iCell) = 0.0_RKIND
          redi_term3(:, :, iCell) = 0.0_RKIND
@@ -224,6 +265,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
             ! Check if neighboring cell exists
             if (cellsOnCell(i, iCell) .eq. nCellsP1) cycle
             iEdge = edgesOnCell(i, iCell)
+            jCell = cellsOnCell(i, iCell)
             cell1 = cellsOnEdge(1, iEdge)
             cell2 = cellsOnEdge(2, iEdge)
             if (cell1 == iCell) then
@@ -232,6 +274,24 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
                iCellSelf = 2
             endif
 
+            if (use_quasi_monotone) then
+               ! expand min/max bnds on edge neighbours
+               do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
+                  kUp = max(minLevelEdgeBot(iEdge), k - 1)
+                  kDn = min(maxLevelEdgeTop(iEdge), k + 1)
+                  do iTr = 1, nTracers
+                     minimumBnd(iTr, k, iCell) = &
+                        min(minimumBnd(iTr, k, iCell), tracers(iTr, kUp, jCell), &
+                                                       tracers(iTr, k, jCell), &
+                                                       tracers(iTr, kDn, jCell))
+                     maximumBnd(iTr, k, iCell) = &
+                        max(maximumBnd(iTr, k, iCell), tracers(iTr, kUp, jCell), &
+                                                       tracers(iTr, k, jCell), &
+                                                       tracers(iTr, kDn, jCell))
+                  end do
+               end do
+            end if
+
             r_tmp = dvEdge(iEdge)/dcEdge(iEdge)
             coef = dvEdge(iEdge)
 
@@ -241,6 +301,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
                                            RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))
 
             do iTr = 1, nTracers
+               ! term 1: div( h grad psi )
                tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1)
                flux = layerThickEdgeMean(k, iEdge)*tracer_turb_flux*r_tmp*RediKappa(iEdge)
                redi_term1(iTr, k, iCell) = redi_term1(iTr, k, iCell) - (1.0_RKIND + term1TaperFactor* &
@@ -248,7 +309,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
 
                if (.not.use_Redi_diag_terms) cycle
 
-               ! Term 2: div( h S dphi/dz)
+               ! term 2: div( h S d/dz psi )
                flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* &
                             0.25_RKIND* &
                             (slopeTriadDown(k, 1, iEdge)*(tracers(iTr, k, cell1) - tracers(iTr, k + 1, cell1)) &
@@ -268,6 +329,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
 
                redi_term2(iTr, k, iCell) = redi_term2(iTr, k, iCell) - edgeSignOnCell(i, iCell)*flux_term2*invAreaCell
                redi_term2_edge(iTr, k, iEdge) = -flux_term2
+
+               ! term 3: zero at surface
             end do
 
             do k = minLevelEdgeBot(iEdge) + 1, maxLevelEdgeTop(iEdge) - 1
@@ -276,17 +339,15 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
                                               RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))
 
                do iTr = 1, nTracers
-                  ! \kappa_2 \nabla \phi on edge
+                  ! term 1: div( h grad psi )
                   tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1)
-
-                  ! div(h \kappa_2 \nabla \phi) at cell center
                   flux = layerThickEdgeMean(k, iEdge)*tracer_turb_flux*r_tmp*RediKappa(iEdge)
-
                   redi_term1(iTr, k, iCell) = redi_term1(iTr, k, iCell) - (1.0_RKIND + term1TaperFactor* &
                                               (kappaRediEdgeVal - 1.0_RKIND))*edgeSignOnCell(i, iCell)*flux*invAreaCell
 
                   if (.not.use_Redi_diag_terms) cycle
 
+                  ! term 2: div( h S d/dz psi )
                   flux_term2 = coef*RediKappa(iEdge)*kappaRediEdgeVal*layerThickEdgeMean(k, iEdge)* &
                                0.25_RKIND* &
                                (slopeTriadUp(k, 1, iEdge)*(tracers(iTr, k - 1, cell1) - tracers(iTr, k, cell1)) &
@@ -301,10 +362,11 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
                                 invAreaCell
                   redi_term2_edge(iTr, k, iEdge) = -flux_term2
 
-                  dTracerDx = (tracers(iTr, k, cell2) - tracers(iTr, k, cell1))*dvEdge(iEdge)*0.25_RKIND
-                  fluxRediZTop(iTr, k) = fluxRediZTop(iTr, k) + &
-                                         0.5_RKIND*(slopeTriadUp(k, iCellSelf, iEdge) + slopeTriadDown(k - 1, iCellSelf, iEdge)) &
-                                         *dTracerDx
+                  ! term 3 : d/dz div( S grad psi )
+                  ! includes weights for edge-to-cell remap
+                  fluxRediZTop(iTr, k) = fluxRediZTop(iTr, k) + 0.125_RKIND*dvEdge(iEdge)*( &
+                     slopeTriadDown(k - 1, iCellSelf, iEdge)*(tracers(iTr, k - 1, cell2) - tracers(iTr, k - 1, cell1)) + &
+                     slopeTriadUp(k, iCellSelf, iEdge)*(tracers(iTr, k, cell2) - tracers(iTr, k, cell1)))
                end do
             end do
 
@@ -313,6 +375,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
                                            RediKappaScaling(k + 1, cell1) + RediKappaScaling(k + 1, cell2))
 
             do iTr = 1, nTracers
+               ! term 1: div( h grad psi )
                tracer_turb_flux = tracers(iTr, k, cell2) - tracers(iTr, k, cell1)
                flux = layerThickEdgeMean(k, iEdge)*tracer_turb_flux*r_tmp*RediKappa(iEdge)
                redi_term1(iTr, k, iCell) = redi_term1(iTr, k, iCell) - (1.0_RKIND + term1TaperFactor* &
@@ -320,6 +383,7 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
 
                if (.not.use_Redi_diag_terms) cycle
 
+               ! term 2: div( h S d/dz psi )
                ! For bottom layer, only use triads pointing up:
                flux_term2 = coef*kappaRediEdgeVal*RediKappa(iEdge)*layerThickEdgeMean(k, iEdge)* &
                             0.25_RKIND* &
@@ -342,32 +406,32 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
                                 *invAreaCell
                redi_term2_edge(iTr, k, iEdge) = -flux_term2
 
-               dTracerDx = (tracers(iTr, k, cell2) - tracers(iTr, k, cell1))*dvEdge(iEdge)*0.25_RKIND
-               fluxRediZTop(iTr, k) = fluxRediZTop(iTr, k) + &
-                                      0.5_RKIND*(slopeTriadUp(k, iCellSelf, iEdge) + slopeTriadDown(k - 1, iCellSelf, iEdge)) &
-                                      *dTracerDx
+               ! term 3: d/dz div( S grad psi )
+               ! includes weights for edge-to-cell remap
+               fluxRediZTop(iTr, k) = fluxRediZTop(iTr, k) + 0.125_RKIND*dvEdge(iEdge)*( &
+                  slopeTriadDown(k - 1, iCellSelf, iEdge)*(tracers(iTr, k - 1, cell2) - tracers(iTr, k - 1, cell1)) + &
+                  slopeTriadUp(k, iCellSelf, iEdge)*(tracers(iTr, k, cell2) - tracers(iTr, k, cell1)))
             end do
          end do ! nEdgesOnCell(iCell)
 
          if ( use_Redi_diag_terms) then
             ! impose no-flux boundary conditions at top and bottom of column
             fluxRediZTop(:, 1:minLevelCell(iCell)) = 0.0_RKIND
-            fluxRediZTop(:, maxLevelCell(iCell) + 1) = 0.0_RKIND
             redi_term3_topOfCell(:, 1:minLevelCell(iCell), iCell) = 0.0_RKIND
             do k = minLevelCell(iCell) + 1, maxLevelCell(iCell)
                redi_term3_topOfCell(:, k, iCell) = fluxRediZTop(:, k) * invAreaCell
             end do
+            fluxRediZTop(:, maxLevelCell(iCell) + 1) = 0.0_RKIND
             redi_term3_topOfCell(:, maxLevelCell(iCell) + 1, iCell) = 0.0_RKIND
 
+            ! take d/dz and remap div( S grad psi ) from edge to cells
             do k = minLevelCell(iCell), maxLevelCell(iCell)
                do iTr = 1, nTracers
-                  ! Add tendency for Term 3: d/dz ( h S grad phi) = ( S grad phi) fluxes
-                  ! 2.0 in next line is because a dot product on a C-grid
-                  ! requires a factor of 1/2 to average to the cell center.
-                  flux_term3 = rediKappaCell(iCell)*2.0_RKIND* &
-                               (RediKappaScaling(k, iCell)*fluxRediZTop(iTr, k)  &
-                               * invAreaCell - RediKappaScaling(k + 1, iCell) &
-                               *fluxRediZTop(iTr, k + 1) * invAreaCell)
+                  flux_term3 = rediKappaCell(iCell) * ( &
+                     RediKappaScaling(k, iCell)* &
+                     redi_term3_topOfCell(iTr, k, iCell) - &
+                     RediKappaScaling(k + 1, iCell)* &
+                     redi_term3_topOfCell(iTr, k + 1, iCell))
 
                   redi_term3(iTr, k, iCell) = redi_term3(iTr, k, iCell) + flux_term3
                end do
@@ -394,26 +458,39 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
       endif
 
       !$omp parallel
-      !$omp do schedule(runtime) private(k, iTr, tempTracer, iEdge)
+      !$omp do schedule(runtime) &
+      !$omp private(k, iTr, tempTracer, iEdge, over, scal, diff)
       do iCell = 1, nCells
          do k = minLevelCell(iCell), maxLevelCell(iCell)
             do iTr = 1, ntracers
-               tempTracer = tracers(iTr, k, iCell) + dt*(redi_term1(iTr, k, iCell) + redi_term2(iTr, k, iCell) + &
-                                                         redi_term3(iTr, k, iCell))/layerThickness(k, iCell)
+               tempTracer = tracers(iTr, k, iCell) + dt*( &
+                  redi_term1(iTr, k, iCell) + redi_term2(iTr, k, iCell) + &
+                  redi_term3(iTr, k, iCell))/layerThickness(k, iCell)
 
-               if (tempTracer < minimumVal(iTr)) then
+               if (tempTracer < minimumVal(iTr) .or. &
+                  (use_quasi_monotone .and. tempTracer < minimumBnd(iTr, k, iCell)) .or. &
+                  (use_quasi_monotone .and. tempTracer > maximumBnd(iTr, k, iCell))) then
 
-                  do i = 1, nEdgesOnCell(iCell)
-                     iEdge = edgesOnCell(i, iCell)
-                     redi_term2_edge(iTr, k, iEdge) = 0.0_RKIND
-                  end do
+                  ! scale cross-fluxes down relative to size of bnds violation
+                  over = max(tempTracer - maximumBnd(iTr, k, iCell), &
+                             minimumBnd(iTr, k, iCell) - tempTracer)
+                  
+                  diff = abs(tempTracer - tracers(iTr, k, iCell))
+
+                  scal = min(1.0_RKIND, max(0.0_RKIND, &
+                     limiter_safety_factor * ( &
+                     1.0_RKIND - over / (diff + 1.0E-16_RKIND))))
+
+                  ! scale to zero if not using monotone limiter
+                  if (.not. use_quasi_monotone) scal = 0.0_RKIND
+
+                  ! cell flux scalings, to be interp. onto edge/levs next pass
+                  redi_flux_scale(iTr, k, iCell) = scal
 
-                  redi_term3_topOfCell(iTr, k, iCell) = 0.0_RKIND
-                  redi_term3_topOfCell(iTr, k + 1, iCell) = 0.0_RKIND
                   if (isActiveTracer) then
                      rediLimiterCount(k, iCell) = rediLimiterCount(k, iCell) + 1
-                  endif
-               endif
+                  end if
+               end if
             end do
          end do
       end do ! iCell
@@ -422,7 +499,8 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
 
       !now go back and reapply all tendencies
       !$omp parallel
-      !$omp do schedule(runtime) private(invAreaCell, k, iTr, i, iEdge)
+      !$omp do schedule(runtime) &
+      !$omp private(invAreaCell, k, kUp, kDn, iTr, i, iEdge, jCell, scal, scalUp, scalDn)
       do iCell = 1, nCells
          invAreaCell = 1.0_RKIND/areaCell(iCell)
          do k = minLevelCell(iCell), maxLevelCell(iCell)
@@ -432,21 +510,48 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
          end do
 
          if (maxLevelCell(iCell) .ge. config_Redi_min_layers_diag_terms) then
+            ! term 3: d/dz div( S grad psi )
             do k = minLevelCell(iCell), maxLevelCell(iCell)
+               kUp = max(minLevelCell(iCell), k - 1)
+               kDn = min(maxLevelCell(iCell), k + 1)
                do iTr = 1, ntracers
-                  tend(iTr, k, iCell) = tend(iTr, k, iCell) + rediKappaCell(iCell)*2.0_RKIND* &
-                                        (RediKappaScaling(k, iCell)* &
-                                         redi_term3_topOfCell(iTr, k, iCell) - &
-                                         RediKappaScaling(k + 1, iCell)*redi_term3_topOfCell(iTr, k + 1, iCell))
+                  ! harmonic-mean; not strictly monotone, but smooth
+                  scalUp = 2.0_RKIND * &
+                     redi_flux_scale(iTr, kUp, iCell) * &
+                     redi_flux_scale(iTr, k, iCell) / ( &
+                     redi_flux_scale(iTr, kUp, iCell) + &
+                     redi_flux_scale(iTr, k, iCell) + 1.0E-16_RKIND)
+
+                  scalDn = 2.0_RKIND * &
+                     redi_flux_scale(iTr, kDn, iCell) * &
+                     redi_flux_scale(iTr, k, iCell) / ( &
+                     redi_flux_scale(iTr, kDn, iCell) + &
+                     redi_flux_scale(iTr, k, iCell) + 1.0E-16_RKIND)
+
+                  tend(iTr, k, iCell) = tend(iTr, k, iCell) + &
+                     rediKappaCell(iCell) * ( &
+                        RediKappaScaling(k, iCell)* &
+                        scalUp*redi_term3_topOfCell(iTr, k, iCell) - &
+                        RediKappaScaling(k + 1, iCell)* &
+                        scalDn*redi_term3_topOfCell(iTr, k + 1, iCell))
                end do
             end do
 
+            ! term 2: div( h S d/dz psi )
             do i = 1, nEdgesOnCell(iCell)
                if (cellsOnCell(i, iCell) .eq. nCellsP1) cycle
                iEdge = edgesOnCell(i, iCell)
+               jCell = cellsOnCell(i, iCell)
                do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
                   do iTr = 1, ntracers
-                     tend(iTr, k, iCell) = tend(iTr, k, iCell) + edgeSignOnCell(i, iCell)*invAreaCell*redi_term2_edge(iTr, k, iEdge)
+                     scal = 2.0_RKIND * &
+                     redi_flux_scale(iTr, k, iCell) * &
+                     redi_flux_scale(iTr, k, jCell) / ( &
+                     redi_flux_scale(iTr, k, iCell) + &
+                     redi_flux_scale(iTr, k, jCell) + 1.0E-16_RKIND)
+
+                     tend(iTr, k, iCell) = tend(iTr, k, iCell) + &
+                        scal*edgeSignOnCell(i, iCell)*invAreaCell*redi_term2_edge(iTr, k, iEdge)
                   end do
                end do
             end do
@@ -463,6 +568,9 @@ subroutine ocn_tracer_hmix_Redi_tend(meshPool, layerThickness, layerThickEdgeMea
       deallocate (redi_term3)
       deallocate (redi_term2_edge)
       deallocate (redi_term3_topOfCell)
+      deallocate (minimumBnd)
+      deallocate (maximumBnd)
+      deallocate (redi_flux_scale)
 
       call mpas_timer_stop("tracer redi")
 
@@ -545,6 +653,12 @@ subroutine ocn_tracer_hmix_Redi_init(domain, err)!{{{
             call mpas_dmpar_finalize(domain%dminfo)
          end if
 
+         ! quasi bounds-preserving limiter
+         use_quasi_monotone = config_Redi_use_quasi_monotone_limiter
+      
+         ! safety factor for limiter: 0 <= fac <= 1
+         limiter_safety_factor = config_Redi_quasi_monotone_safety_factor
+
          ! Initialize horizontal taper
          if (config_Redi_horizontal_taper == 'none' .or. &
              config_Redi_horizontal_taper == 'RossbyRadius') then