diff --git a/src/core_ocean/adc_mixing/Registry_adc_mixing_fields_opts.xml b/src/core_ocean/adc_mixing/Registry_adc_mixing_fields_opts.xml
index fb2f152c5b..3482034ed6 100644
--- a/src/core_ocean/adc_mixing/Registry_adc_mixing_fields_opts.xml
+++ b/src/core_ocean/adc_mixing/Registry_adc_mixing_fields_opts.xml
@@ -11,6 +11,10 @@
 	      description="flag to use buoyancy length scale from previous closure, temporary"
 	      possible_values=".true. or .false."
 	      />
+  <nml_option name="config_use_grad_diff_horiz_mom_flux" type="logical" default_value=".false." units="unitless"
+	      description="flag to use down gradient diffusion closure for horizontal momentum fluxes"
+	      possible_values=".true. or .false."
+	      />
   <nml_option name="config_adc_truncate_tend" type="logical" default_value=".false." units="unitless"
 	      description="flag to truncate tendency terms"
 	      possible_values=".true. or .false."
@@ -66,15 +70,19 @@
   <nml_option name="config_adc_alpha_tracer1" type="real" default_value="0.6" units="unitless"
 	      description="alpha_scalar1 first rapid pressure coefficient in tracer flux budgets (Mironov 2001)"
 	      possible_values="suggested values: 0.6 (Canuto 2001) or equal to alpha_scalar2, >0, <1 (Harcourt 2013)"
-  />
+              />
   <nml_option name="config_adc_alpha_tracer2" type="real" default_value="1.0" units="unitless"
 	      description="alpha_scalar2 rapid pressure coefficient in tracer flux budgets (Mironov 2001)"
 	      possible_values="suggested values: 1.0 (Canuto 2001) or equal to alpha_scalar2, >0, <1 (Harcourt 2013)"
-  />
+              />
   <nml_option name="config_adc_c11" type="real" default_value="0.1" units="unitless"
 	      description="c_11 is the buoyancy contribution in the pressure correlation in w3 eqn"
 	      possible_values="values near 0.1"
 	      />
+  <nml_option name="config_adc_Cmom_flux" type="real" default_value="0.1" units="unitless"
+	      description="horizontal momentum flux down gradient coefficient"
+	      possible_values="values near 0.1"
+	      />
   <nml_option name="config_adc_Cmom" type="real" default_value="0.0" units="unitless"
 	      description="third order down gradient coefficient, from lappen and randall 2001b"
 	      possible_values="values near 0.5"
diff --git a/src/core_ocean/shared/mpas_ocn_adcReconstruct.F b/src/core_ocean/shared/mpas_ocn_adcReconstruct.F
index 3da2158cd8..5903c84907 100644
--- a/src/core_ocean/shared/mpas_ocn_adcReconstruct.F
+++ b/src/core_ocean/shared/mpas_ocn_adcReconstruct.F
@@ -488,11 +488,13 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
                  lenav = 0.5_RKIND*(length(k,iCell) + length(k+1,iCell))
                  diff = C_mom * sqrt(0.5_RKIND*(KE + KEp1)) * lenav
                  dz = ze(k,iCell) - ze(k+1,iCell)
-                 uw2(k,iCell) = -diff*(uw(i1,k,iCell) - uw(i1,k+1,iCell)) / dz
-                 vw2(k,iCell) = -diff*(vw(i1,k,iCell) - vw(i1,k+1,iCell)) / dz
                  u2w(k,iCell) = -diff*(u2(i1,k,iCell) - u2(i1,k+1,iCell)) / dz
                  v2w(k,iCell) = -diff*(v2(i1,k,iCell) - v2(i1,k+1,iCell)) / dz
-                 uvw(k,iCell) = -diff*(uv(i1,k,iCell) - uv(i1,k+1,iCell)) / dz
+                 if (.not. config_use_grad_diff_horiz_mom_flux) then
+                    uw2(k,iCell) = -diff*(uw(i1,k,iCell) - uw(i1,k+1,iCell)) / dz
+                    vw2(k,iCell) = -diff*(vw(i1,k,iCell) - vw(i1,k+1,iCell)) / dz
+                    uvw(k,iCell) = -diff*(uv(i1,k,iCell) - uv(i1,k+1,iCell)) / dz
+                 endif
 
                  diff = C_therm*sqrt(0.5_RKIND*(KE + KEp1)) * lenav
                  uwt(k,iCell) = -diff*(ut(i1,k,iCell) - ut(i1,k+1,iCell)) / dz
@@ -675,6 +677,7 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
                  endif
 
 
+                 if (.not. config_use_grad_diff_horiz_mom_flux) then
 !      CALCULATE UW TENDENCY TERMS
 !      The uw tendency is the sum of the following terms (which are saved individually)
 !                   (1) turbulent transport of uw [ -d(uww)/dz ]
@@ -683,30 +686,30 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
 !                   (4) buoyancy production of uw [ u*buoyancy]
 !                   (5) dissipation of uw [ -kappa_FL*[d(uw)/dz]^2 ]
 !      The tendency term is truncated if config_adc_truncate_tend flag it switched on
-                 uwtend1(k,iCell) = -(uw2(k-1,iCell) - uw2(k,iCell)) / dzmid
+                    uwtend1(k,iCell) = -(uw2(k-1,iCell) - uw2(k,iCell)) / dzmid
 
-                 uwtend2(k,iCell) = -tauVel(k,iCell)*uw(i1,k,iCell) + &
-                   0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND +  &
-                   (alpha_1 - alpha_2)*u2(i1,k,iCell) + (alpha_1 +  &
-                   alpha_2)*w2(i1,k,iCell))*Uz + 0.5_RKIND*(alpha_1 - alpha_2)*    &
-                   uv(i1,k,iCell)*Vz - C_b*grav*(alphaT(k,iCell)*   &
-                   ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell))
+                    uwtend2(k,iCell) = -tauVel(k,iCell)*uw(i1,k,iCell) &
+                      + 0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND +  &
+                      (alpha_1 - alpha_2)*u2(i1,k,iCell) + (alpha_1 +  &
+                      alpha_2)*w2(i1,k,iCell))*Uz + 0.5_RKIND*(alpha_1 - alpha_2)*    &
+                      uv(i1,k,iCell)*Vz - C_b*grav*(alphaT(k,iCell)*   &
+                      ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell))
 
-                 uwtend3(k,iCell) = - w2(i1,k,iCell)*Uz
+                      uwtend3(k,iCell) = - w2(i1,k,iCell)*Uz
 
-                 uwtend4(k,iCell) = grav*(alphaT(k,iCell)*   &
-                   ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell))
+                      uwtend4(k,iCell) = grav*(alphaT(k,iCell)*   &
+                        ut(i1,k,iCell) - betaS(k,iCell)*us(i1,k,iCell))
 
-                 uwtend5(k,iCell) = kappa_FL*(uw(i1,k-1,iCell) - uw(i1,k+1,iCell)) / &
-                   (ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND
+                      uwtend5(k,iCell) = kappa_FL*(uw(i1,k-1,iCell) - uw(i1,k+1,iCell)) / &
+                        (ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND
 !      Sum the above tendency terms to find the total uw tendency.
-                 uwtend(i3_f,k,iCell) = uwtend1(k,iCell) + uwtend2(k,iCell) + &
-                   uwtend3(k,iCell) + uwtend4(k,iCell) + uwtend5(k,iCell)
+                      uwtend(i3_f,k,iCell) = uwtend1(k,iCell) + uwtend2(k,iCell) + &
+                        uwtend3(k,iCell) + uwtend4(k,iCell) + uwtend5(k,iCell)
 
-                 if (config_adc_truncate_tend) then
-                    uwtend(i3_f,k,iCell) = FLOAT (INT(uwtend(i3_f,k,iCell) * adcRound &
-                      + 0.5_RKIND)) / adcRound
-                 endif
+                      if (config_adc_truncate_tend) then
+                         uwtend(i3_f,k,iCell) = FLOAT (INT(uwtend(i3_f,k,iCell) * adcRound &
+                           + 0.5_RKIND)) / adcRound
+                      endif
 
 
 !      CALCULATE VW TENDENCY TERMS
@@ -717,30 +720,30 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
 !                   (4) buoyancy production of vw [ v*buoyancy]
 !                   (5) dissipation of vw [ -kappa_FL*[d(vw)/dz]^2 ]
 !      The tendency term is truncated if config_adc_truncate_tend flag it switched on
-                 vwtend1(k,iCell) = -(vw2(k-1,iCell) - vw2(k,iCell)) / dzmid
+                      vwtend1(k,iCell) = -(vw2(k-1,iCell) - vw2(k,iCell)) / dzmid
 
-                 vwtend2(k,iCell) = -tauVel(k,iCell)*vw(i1,k,iCell) + &
-                   0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND +  &
-                   (alpha_1 - alpha_2)*v2(i1,k,iCell) + (alpha_1 +  &
-                   alpha_2)*w2(i1,k,iCell))*Vz + 0.5_RKIND*(alpha_1 - alpha_2)*    &
-                   uv(i1,k,iCell)*Uz - C_b*grav*(alphaT(k,iCell)*   &
-                   vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell))
+                      vwtend2(k,iCell) = -tauVel(k,iCell)*vw(i1,k,iCell)  &
+                        + 0.5_RKIND*((alpha_0 - 4.0_RKIND*alpha_1/3.0_RKIND)*KE**2.0_RKIND +  &
+                        (alpha_1 - alpha_2)*v2(i1,k,iCell) + (alpha_1 +  &
+                        alpha_2)*w2(i1,k,iCell))*Vz + 0.5_RKIND*(alpha_1 - alpha_2)*    &
+                        uv(i1,k,iCell)*Uz - C_b*grav*(alphaT(k,iCell)*   &
+                        vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell))
 
-                 vwtend3(k,iCell) = - w2(i1,k,iCell)*Vz
+                      vwtend3(k,iCell) = - w2(i1,k,iCell)*Vz
 
-                 vwtend4(k,iCell) = grav*(alphaT(k,iCell)*   &
-                   vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell))
+                      vwtend4(k,iCell) = grav*(alphaT(k,iCell)*   &
+                        vt(i1,k,iCell) - betaS(k,iCell)*vs(i1,k,iCell))
 
-                 vwtend5(k,iCell) = kappa_FL*(vw(i1,k-1,iCell) - vw(i1,k+1,iCell)) / &
-                   (ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND
+                      vwtend5(k,iCell) = kappa_FL*(vw(i1,k-1,iCell) - vw(i1,k+1,iCell)) / &
+                        (ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND
 !      Sum the above tendency terms to find the total vw tendency.
-                 vwtend(i3_f,k,iCell) = vwtend1(k,iCell) + vwtend2(k,iCell) + &
-                   vwtend3(k,iCell) + vwtend4(k,iCell) + vwtend5(k,iCell)
+                      vwtend(i3_f,k,iCell) = vwtend1(k,iCell) + vwtend2(k,iCell) + &
+                        vwtend3(k,iCell) + vwtend4(k,iCell) + vwtend5(k,iCell)
 
-                 if (config_adc_truncate_tend) then
-                    vwtend(i3_f,k,iCell) = FLOAT (INT(vwtend(i3_f,k,iCell) * adcRound &
-                      + 0.5_RKIND)) / adcRound
-                 endif
+                      if (config_adc_truncate_tend) then
+                         vwtend(i3_f,k,iCell) = FLOAT (INT(vwtend(i3_f,k,iCell) * adcRound &
+                           + 0.5_RKIND)) / adcRound
+                      endif
 
 
 !      CALCULATE UV TENDENCY TERMS (these are not saved individually)
@@ -750,20 +753,20 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
 !                   (3) shear production of uv [uw*dV/dz + vw*dU/dz]
 !                   (4) dissipation of uv [ -kappa_FL*[d(uv)/dz]^2 ]
 !      The tendency term is truncated if config_adc_truncate_tend flag it switched on
-                uvtend(i3_f,k,iCell) = -(uvw(k-1,iCell) - uvw(k,iCell)) / dz & ! transport
-                   - tauVel(k,iCell)*uv(i1,k,iCell)  & ! Rotta
-                   + 0.5_RKIND*(alpha_1+alpha_2)* &
-                   (uw(i1,k,iCell)*Vz + vw(i1,k,iCell)*Uz) & ! rapid
-                   - (uw(i1,k,iCell)*Vz + vw(i1,k,iCell)*Uz) & ! shear production
-                   + kappa_VAR*(uv(i1,k-1,iCell) - uv(i1,k+1,iCell)) / &
-                   (ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND ! dissipation
+                      uvtend(i3_f,k,iCell) = -(uvw(k-1,iCell) - uvw(k,iCell)) / dz & ! transport
+                        - tauVel(k,iCell)*uv(i1,k,iCell)  & ! Rotta
+                        + 0.5_RKIND*(alpha_1+alpha_2)* &
+                        (uw(i1,k,iCell)*Vz + vw(i1,k,iCell)*Uz) & ! rapid
+                        - (uw(i1,k,iCell)*Vz + vw(i1,k,iCell)*Uz) & ! shear production
+                        + kappa_VAR*(uv(i1,k-1,iCell) - uv(i1,k+1,iCell)) / &
+                        (ze(k-1,iCell) - ze(k+1,iCell))**2.0_RKIND ! dissipation
 
 !      Truncate the uv tendency term if the truncation flag is switched on
-                 if (config_adc_truncate_tend) then
-                    uvtend(i3_f,k,iCell) = FLOAT (INT(uvtend(i3_f,k,iCell) * adcRound &
-                      + 0.5_RKIND)) / adcRound
-                 endif
-
+                      if (config_adc_truncate_tend) then
+                         uvtend(i3_f,k,iCell) = FLOAT (INT(uvtend(i3_f,k,iCell) * adcRound &
+                           + 0.5_RKIND)) / adcRound
+                      endif
+                   endif
 
 !      CALCULATE U2 TENDENCY TERMS
 !      The u2 tendency is the sum of the following terms (which are saved individually)
@@ -1064,19 +1067,26 @@ subroutine compute_ADC_tends(nCells,nVertLevels, nTracers, dt,activeTracers, uve
                     v2(i2,k,iCell) = epsilon
                  endif
 
-                 uw(i2,k,iCell) = uw(i1,k,iCell) + dt_small*(Cw3 * uwtend(i3_f,k,iCell) + &
-                   Cw2 * uwtend(i2_f,k,iCell) + Cw1 * uwtend(i1_f,k,iCell))
-                 vw(i2,k,iCell) = vw(i1,k,iCell) + dt_small*(Cw3 * vwtend(i3_f,k,iCell) + &
-                   Cw2 * vwtend(i2_f,k,iCell) + Cw1 * vwtend(i1_f,k,iCell))
-                 uv(i2,k,iCell) = uv(i1,k,iCell) + dt_small*(Cw3 * uvtend(i3_f,k,iCell) + &
-                   Cw2 * uvtend(i2_f,k,iCell) + Cw1 * uvtend(i1_f,k,iCell))
+                 if (config_use_grad_diff_horiz_mom_flux) then
+                    Uz = (uvel(k-1,iCell) - uvel(k,iCell)) / dzmid
+                    Vz = (vvel(k-1,iCell) - vvel(k,iCell)) / dzmid
+                    uw(i2,k,iCell) = -C_mom_flux*sqrt(0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell) + w2(i1,k,iCell)))*length(k,iCell)*Uz
+                    vw(i2,k,iCell) = -C_mom_flux*sqrt(0.5_RKIND*(u2(i1,k,iCell) + v2(i1,k,iCell) + w2(i1,k,iCell)))*length(k,iCell)*Vz
+                 else
+                    uw(i2,k,iCell) = uw(i1,k,iCell) + dt_small*(Cw3 * uwtend(i3_f,k,iCell) + &
+                      Cw2 * uwtend(i2_f,k,iCell) + Cw1 * uwtend(i1_f,k,iCell))
+                    vw(i2,k,iCell) = vw(i1,k,iCell) + dt_small*(Cw3 * vwtend(i3_f,k,iCell) + &
+                      Cw2 * vwtend(i2_f,k,iCell) + Cw1 * vwtend(i1_f,k,iCell))
+                    uv(i2,k,iCell) = uv(i1,k,iCell) + dt_small*(Cw3 * uvtend(i3_f,k,iCell) + &
+                      Cw2 * uvtend(i2_f,k,iCell) + Cw1 * uvtend(i1_f,k,iCell))
+                 end if
                  ut(i2,k,iCell) = ut(i1,k,iCell) + dt_small*(Cw3 * uttend(i3_f,k,iCell) + &
                    Cw2 * uttend(i2_f,k,iCell) + Cw1 * uttend(i1_f,k,iCell))
+                 vt(i2,k,iCell) = vt(i1,k,iCell) + dt_small*(Cw3 * vttend(i3_f,k,iCell) + &
+                   Cw2 * vttend(i2_f,k,iCell) + Cw1 * vttend(i1_f,k,iCell))
                  wt(i2,k,iCell) = (wt(i1,k,iCell) + dt_small*(Cw3 * wttend(i3_f,k,iCell) + &
                    Cw2 * wttend(i2_f,k,iCell) + Cw1 * wttend(i1_f,k,iCell))) /              &
                    (1.0_RKIND + dt_small*tau_tracer(k,iCell))
-                 vt(i2,k,iCell) = vt(i1,k,iCell) + dt_small*(Cw3 * vttend(i3_f,k,iCell) + &
-                   Cw2 * vttend(i2_f,k,iCell) + Cw1 * vttend(i1_f,k,iCell))
                  us(i2,k,iCell) = us(i1,k,iCell) + dt_small*(Cw3 * ustend(i3_f,k,iCell) + &
                    Cw2 * ustend(i2_f,k,iCell) + Cw1 * ustend(i1_f,k,iCell))
                  vs(i2,k,iCell) = vs(i1,k,iCell) + dt_small*(Cw3 * vstend(i3_f,k,iCell) + &
diff --git a/src/core_ocean/shared/mpas_ocn_turbulence.F b/src/core_ocean/shared/mpas_ocn_turbulence.F
index 88d89e9a1a..d457e41b0b 100644
--- a/src/core_ocean/shared/mpas_ocn_turbulence.F
+++ b/src/core_ocean/shared/mpas_ocn_turbulence.F
@@ -45,7 +45,7 @@ module ocn_turbulence
 
    real(kind=RKIND),public :: epsilon, &! small value below which w2, v2, u2 cannot go
                        sigmat, Ko, c_b, c_b_tracer, alpha_tracer1, alpha_tracer2, c11, &
-                       alpha_0, alpha_1, alpha_2, B1, Kt, grav, c_mom, &
+                       alpha_0, alpha_1, alpha_2, B1, Kt, grav, c_mom, c_mom_flux, &
                        c_therm, c_mom_w3, c_epsilon, c_slow, c_slow_tracer, slow_w_factor, &
                        kappa_FL, kappa_w3, kappa_VAR, Cww_D, Cww_E, adcRound
 
@@ -162,6 +162,7 @@ subroutine ocn_turbulenceCreate(domain)
       Kt = 0.4_RKIND
       c_b = config_adc_c_b
       c11 = config_adc_c11
+      c_mom_flux = config_adc_Cmom_flux
       c_mom = config_adc_Cmom
       c_therm = config_adc_Ctherm
       c_mom_w3 = config_adc_Cmom_w3
@@ -560,6 +561,7 @@ subroutine ocn_turbulenceCreate(domain)
       !$acc                   Kt,           &
       !$acc                   c_1,          &
       !$acc                   c_2,          &
+      !$acc                   c_mom_flux,   &
       !$acc                   c_mom,        &
       !$acc                   c_therm,      &
       !$acc                   c_mom_w3,     &
@@ -766,6 +768,7 @@ subroutine ocn_turbulenceUpdateFields(domain)
       !$acc                   Kt,           &
       !$acc                   c_1,          &
       !$acc                   c_2,          &
+      !$acc                   c_mom_flux,   &
       !$acc                   c_mom,        &
       !$acc                   c_therm,      &
       !$acc                   c_mom_w3,     &
@@ -1186,6 +1189,7 @@ subroutine ocn_turbulenceDestroy(ierr)
       alpha_tracer2 = 0.0_RKIND
       B1 = 0.0_RKIND
       Kt = 0.0_RKIND
+      c_mom_flux = 0.0_RKIND
       c_mom = 0.0_RKIND
       c_therm = 0.0_RKIND
       c_mom_w3 = 0.0_RKIND
@@ -1220,6 +1224,7 @@ subroutine ocn_turbulenceDestroy(ierr)
       !$acc                   Kt,           &
       !$acc                   c_1,          &
       !$acc                   c_2,          &
+      !$acc                   c_mom_flux,   &
       !$acc                   c_mom,        &
       !$acc                   c_therm,      &
       !$acc                   c_mom_w3,     &