diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vmix.F b/components/mpas-ocean/src/shared/mpas_ocn_vmix.F
index 8b7894c2070a..87f23d65aa1f 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_vmix.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_vmix.F
@@ -1069,13 +1069,6 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_loglaw(forcingPool, dt,
                             von_karman_sq / &
                             (log(1.0_RKIND + &
                                  (zMidEdge/bottom_roughness)))**2))
-           ! For some reason we can use k-1 even when N=Nsurf
-           kineticEnergyTopOfEdge = 0.5_RKIND * (kineticEnergyCell(k,cell1) + kineticEnergyCell(k,cell2) + &
-                                                 kineticEnergyCell(k-1,cell1) + kineticEnergyCell(k-1,cell2))
-           vertViscTopOfEdge(k,iEdge) = max(vertViscTopOfEdge(k,iEdge), &
-                                            sqrt(CdTemp(k) * kineticEnergyTopOfEdge) * &
-                                            von_karman * (zMidEdge + bottom_roughness) &
-                                           )
          enddo
 
          ! one active layer
@@ -1084,7 +1077,14 @@ subroutine ocn_vel_vmix_tend_implicit_spatially_variable_loglaw(forcingPool, dt,
                 / (1.0_RKIND + dt*CdTemp(N) &
                    * sqrt(kineticEnergyCell(N,cell1) + kineticEnergyCell(N,cell2)) / layerThickEdgeDrag(N,iEdge) )
          else
-
+           do k = Nsurf+1, N
+              kineticEnergyTopOfEdge = 0.5_RKIND * (kineticEnergyCell(k,cell1) + kineticEnergyCell(k,cell2) + &
+                                                    kineticEnergyCell(k-1,cell1) + kineticEnergyCell(k-1,cell2))
+              vertViscTopOfEdge(k,iEdge) = max(vertViscTopOfEdge(k,iEdge), &
+                                               sqrt(CdTemp(k) * kineticEnergyTopOfEdge) * &
+                                               von_karman * (zMidEdge + bottom_roughness) &
+                                              )
+           enddo
            ! tridiagonal matrix algorithm
            C(Nsurf)     = -2.0_RKIND*dt*vertViscTopOfEdge(Nsurf+1,iEdge) &
                       / (layerThickEdgeMean(Nsurf,iEdge) + layerThickEdgeMean(Nsurf+1,iEdge)) &