diff --git a/pkg/CVMix-src b/pkg/CVMix-src
index 5e3d0cee52..87c3c0cb3d 160000
--- a/pkg/CVMix-src
+++ b/pkg/CVMix-src
@@ -1 +1 @@
-Subproject commit 5e3d0cee52c6ee66a62277abf2fee8956ed85890
+Subproject commit 87c3c0cb3d0b5ac606f938d413c64a88ec781682
diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90
index 71dcbb9170..b0c9d815f4 100644
--- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90
+++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90
@@ -114,6 +114,7 @@ module MOM_CVMix_KPP
   logical :: LT_K_Enhancement          !< Flags if enhancing mixing coefficients due to LT
   integer :: LT_K_Shape                !< Integer for constant or shape function enhancement
   integer :: LT_K_Method               !< Integer for mixing coefficients LT method
+  real    :: KPP_CVt2                  !< Parameter for Stokes MOST convection entrainment
   real    :: KPP_K_ENH_FAC             !< Factor to multiply by K if Method is CONSTANT [nondim]
   logical :: LT_Vt2_Enhancement        !< Flags if enhancing Vt2 due to LT
   integer :: LT_VT2_METHOD             !< Integer for Vt2 LT method
@@ -436,6 +437,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
       case default ; call MOM_error(FATAL,"KPP_init: "//&
                     "Unrecognized KPP_LT_K_SHAPE option: "//trim(string))
     end select
+    call get_param(paramFile, mdl, "KPP_CVt2", CS%KPP_CVt2, &
+                   'Parameter for Stokes MOST convection entrainment', &
+                   units="nondim", default=1.6)
     call get_param(paramFile, mdl, "KPP_LT_K_METHOD", string ,                   &
                    'Method to enhance mixing coefficient in KPP. '//             &
                    'Valid options are: \n'//                                     &
@@ -529,6 +533,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
                        minVtsqr=US%L_T_to_m_s**2*CS%minVtsqr, &
                        vonKarman=CS%vonKarman,             &
                        surf_layer_ext=CS%surf_layer_ext,   &
+                       CVt2=CS%KPP_CVt2,                   &
                        interp_type=CS%interpType,          &
                        interp_type2=CS%interpType2,        &
                        lEkman=CS%computeEkman,             &
@@ -1162,8 +1167,10 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
         enddo
 
         if (CS%StokesMOST) then
+          ! if k=1, want buoyFlux(i,j,1) - buoyFlux(i,j,2), otherwise
+          ! subtract average of buoyFlux(i,j,k) and buoyFlux(i,j,k+1)
           surfBuoyFlux   = buoy_scale * &
-                          (buoyFlux(i,j,1) - 0.5*(buoyFlux(i,j,k)+buoyFlux(i,j,k+1)) )
+                          (buoyFlux(i,j,1) - 0.5*(buoyFlux(i,j,max(2,k))+buoyFlux(i,j,k+1)) )
           surfBuoyFlux2(k) = surfBuoyFlux
           call Compute_StokesDrift(i,j, iFaceHeight(k),iFaceHeight(k+1), &
               uS_Hi(k+1), vS_Hi(k+1), uS_H(k), vS_H(k), uSbar_H(k), vSbar_H(k), Waves)