diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90
index 73740e5412..dd1f902651 100644
--- a/src/parameterizations/lateral/MOM_internal_tides.F90
+++ b/src/parameterizations/lateral/MOM_internal_tides.F90
@@ -1505,6 +1505,7 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
           threshold_renorm_N, & ! Maximum allowable error on N profile [H T-1 ~> m s-1 or kg m-2 s-1]
           threshold_verif       ! Maximum allowable error on verification [nondim]
 
+  logical :: apply_Kd_max
   logical :: non_Bous ! fully Non-Boussinesq
   integer :: i, k, is, ie, nz
 
@@ -1512,6 +1513,8 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
 
   non_Bous = .not.(GV%Boussinesq .or. GV%semi_Boussinesq)
 
+  apply_Kd_max = (Kd_max >= 0.0)
+
   h_d = CS%Int_tide_decay_scale
   h_s = CS%Int_tide_decay_scale_slope
   I_h_d = 1 / h_d
@@ -1703,7 +1706,13 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
           Kd_leak_lay(k) = 0.
         endif
         ! add to total Kd in layer
-        if (CS%update_Kd) Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_leak_lay(k), Kd_max)
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_leak_lay(k), Kd_max)
+          else
+            Kd_lay(i,k) = Kd_lay(i,k) + Kd_leak_lay(k)
+          endif
+        endif
       enddo
     endif
 
@@ -1725,7 +1734,13 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
           Kd_Froude_lay(k) = 0.
         endif
         ! add to total Kd in layer
-        if (CS%update_Kd) Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_Froude_lay(k), Kd_max)
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_Froude_lay(k), Kd_max)
+          else
+            Kd_lay(i,k) = Kd_lay(i,k) + Kd_Froude_lay(k)
+          endif
+        endif
       enddo
     endif
 
@@ -1747,7 +1762,13 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
           Kd_itidal_lay(k) = 0.
         endif
         ! add to total Kd in layer
-        if (CS%update_Kd) Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_itidal_lay(k), Kd_max)
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_itidal_lay(k), Kd_max)
+          else
+            Kd_lay(i,k) = Kd_lay(i,k) + Kd_itidal_lay(k)
+          endif
+        endif
       enddo
     endif
 
@@ -1769,7 +1790,13 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
           Kd_slope_lay(k) = 0.
         endif
         ! add to total Kd in layer
-        if (CS%update_Kd) Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_slope_lay(k), Kd_max)
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_slope_lay(k), Kd_max)
+          else
+            Kd_lay(i,k) = Kd_lay(i,k) + Kd_slope_lay(k)
+          endif
+        endif
       enddo
     endif
 
@@ -1791,7 +1818,13 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
           Kd_quad_lay(k) = 0.
         endif
         ! add to total Kd in layer
-        if (CS%update_Kd) Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_quad_lay(k), Kd_max)
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_lay(i,k) = Kd_lay(i,k) + min(Kd_quad_lay(k), Kd_max)
+          else
+            Kd_lay(i,k) = Kd_lay(i,k) + Kd_quad_lay(k)
+          endif
+        endif
       enddo
     endif
 
@@ -1801,7 +1834,13 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
         if (k>1)    Kd_leak(i,K) = 0.5*Kd_leak_lay(k-1)
         if (k<nz+1) Kd_leak(i,K) = Kd_leak(i,K) + 0.5*Kd_leak_lay(k)
         ! add to Kd_int
-        if (CS%update_Kd) Kd_int(i,K) = Kd_int(i,K) + min(Kd_leak(i,K), Kd_max)
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_int(i,K) = Kd_int(i,K) + min(Kd_leak(i,K), Kd_max)
+          else
+            Kd_int(i,K) = Kd_int(i,K)
+          endif
+        endif
       enddo
     endif
 
@@ -1810,8 +1849,14 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
         if (k>1)    Kd_itidal(i,K) = 0.5*Kd_itidal_lay(k-1)
         if (k<nz+1) Kd_itidal(i,K) = Kd_itidal(i,K) + 0.5*Kd_itidal_lay(k)
         ! add to Kd_int
-        if (CS%update_Kd) Kd_int(i,K) = Kd_int(i,K) + min(Kd_itidal(i,K), Kd_max)
-      enddo
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_int(i,K) = Kd_int(i,K) + min(Kd_itidal(i,K), Kd_max)
+          else
+            Kd_int(i,K) = Kd_int(i,K)
+          endif
+        endif
+     enddo
     endif
 
     if (CS%apply_Froude_drag) then
@@ -1819,8 +1864,14 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
         if (k>1)    Kd_Froude(i,K) = 0.5*Kd_Froude_lay(k-1)
         if (k<nz+1) Kd_Froude(i,K) = Kd_Froude(i,K) + 0.5*Kd_Froude_lay(k)
         ! add to Kd_int
-        if (CS%update_Kd) Kd_int(i,K) = Kd_int(i,K) + min(Kd_Froude(i,K), Kd_max)
-      enddo
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_int(i,K) = Kd_int(i,K) + min(Kd_Froude(i,K), Kd_max)
+          else
+            Kd_int(i,K) = Kd_int(i,K)
+          endif
+        endif
+     enddo
     endif
 
     if (CS%apply_residual_drag) then
@@ -1828,8 +1879,14 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
         if (k>1)    Kd_slope(i,K) = 0.5*Kd_slope_lay(k-1)
         if (k<nz+1) Kd_slope(i,K) = Kd_slope(i,K) + 0.5*Kd_slope_lay(k)
         ! add to Kd_int
-        if (CS%update_Kd) Kd_int(i,K) = Kd_int(i,K) + min(Kd_slope(i,K), Kd_max)
-      enddo
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_int(i,K) = Kd_int(i,K) + min(Kd_slope(i,K), Kd_max)
+          else
+            Kd_int(i,K) = Kd_int(i,K)
+          endif
+        endif
+     enddo
     endif
 
     if (CS%apply_bottom_drag) then
@@ -1837,8 +1894,14 @@ subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2
         if (k>1)    Kd_quad(i,K) = 0.5*Kd_quad_lay(k-1)
         if (k<nz+1) Kd_quad(i,K) = Kd_quad(i,K) + 0.5*Kd_quad_lay(k)
         ! add to Kd_int
-        if (CS%update_Kd) Kd_int(i,K) = Kd_int(i,K) + min(Kd_quad(i,K), Kd_max)
-      enddo
+        if (CS%update_Kd) then
+          if (apply_Kd_max) then
+            Kd_int(i,K) = Kd_int(i,K) + min(Kd_quad(i,K), Kd_max)
+          else
+            Kd_int(i,K) = Kd_int(i,K)
+          endif
+        endif
+     enddo
     endif
   enddo ! i-loop
 
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index c6124e0c9a..39b1084ac6 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -634,21 +634,40 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
       if (CS%id_Kd_slope > 0) then ; do K=1,nz+1 ; do i=is,ie
         dd%Kd_slope(i,j,K) = Kd_slope_2d(i,K)
       enddo ; enddo ; endif
-      if (associated (VBF%Kd_leak)) then ; do K=1,nz+1 ; do i=is,ie
-        VBF%Kd_leak(i,j,K) = min(Kd_leak_2d(i,K), CS%Kd_max)
-      enddo ; enddo ; endif
-      if (associated (VBF%Kd_quad)) then ; do K=1,nz+1 ; do i=is,ie
-        VBF%Kd_quad(i,j,K) = min(Kd_quad_2d(i,K), CS%Kd_max)
-      enddo ; enddo ; endif
-      if (associated (VBF%Kd_itidal)) then ; do K=1,nz+1 ; do i=is,ie
-        VBF%Kd_itidal(i,j,K) = min(Kd_itidal_2d(i,K), CS%Kd_max)
-      enddo ; enddo ; endif
-      if (associated (VBF%Kd_Froude)) then ; do K=1,nz+1 ; do i=is,ie
-        VBF%Kd_Froude(i,j,K) = min(Kd_Froude_2d(i,K), CS%Kd_max)
-      enddo ; enddo ; endif
-      if (associated (VBF%Kd_slope)) then ; do K=1,nz+1 ; do i=is,ie
-        VBF%Kd_slope(i,j,K) = min(Kd_slope_2d(i,K), CS%Kd_max)
-      enddo ; enddo ; endif
+      ! Only apply diffusivity cap if Kd_max is positive
+      if (CS%Kd_max < 0.0) then
+        if (associated (VBF%Kd_leak)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_leak(i,j,K) = Kd_leak_2d(i,K)
+        enddo ; enddo ; endif
+        if (associated (VBF%Kd_quad)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_quad(i,j,K) = Kd_quad_2d(i,K)
+        enddo ; enddo ; endif
+        if (associated (VBF%Kd_itidal)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_itidal(i,j,K) = Kd_itidal_2d(i,K)
+        enddo ; enddo ; endif
+        if (associated (VBF%Kd_Froude)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_Froude(i,j,K) = Kd_Froude_2d(i,K)
+        enddo ; enddo ; endif
+        if (associated (VBF%Kd_slope)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_slope(i,j,K) = Kd_slope_2d(i,K)
+        enddo ; enddo ; endif
+      else
+        if (associated (VBF%Kd_leak)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_leak(i,j,K) = min(Kd_leak_2d(i,K), CS%Kd_max)
+        enddo ; enddo ; endif
+        if (associated (VBF%Kd_quad)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_quad(i,j,K) = min(Kd_quad_2d(i,K), CS%Kd_max)
+        enddo ; enddo ; endif
+        if (associated (VBF%Kd_itidal)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_itidal(i,j,K) = min(Kd_itidal_2d(i,K), CS%Kd_max)
+        enddo ; enddo ; endif
+        if (associated (VBF%Kd_Froude)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_Froude(i,j,K) = min(Kd_Froude_2d(i,K), CS%Kd_max)
+        enddo ; enddo ; endif
+        if (associated (VBF%Kd_slope)) then ; do K=1,nz+1 ; do i=is,ie
+          VBF%Kd_slope(i,j,K) = min(Kd_slope_2d(i,K), CS%Kd_max)
+        enddo ; enddo ; endif
+      endif
 
       if (CS%id_prof_leak > 0) then ; do k=1,nz; do i=is,ie
         dd%prof_leak(i,j,k) = prof_leak_2d(i,k)