diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90
index d2f2072223..6e5baa6fdd 100644
--- a/src/user/MOM_wave_interface.F90
+++ b/src/user/MOM_wave_interface.F90
@@ -798,8 +798,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step)
         MidPoint = 0.0
         do k = 1,GV%ke
           Top = Bottom
-          MidPoint = Bottom - 0.25*(dz(I,j,k)+dz(I-1,j,k))
-          Bottom = Bottom - 0.5*(dz(I,j,k)+dz(I-1,j,k))
+          MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i+1,j,k))
+          Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i+1,j,k))
           CS%Us_x(I,j,k) = CS%TP_STKX0*exp(MidPoint*DecayScale)
         enddo
       enddo
@@ -810,8 +810,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step)
         MidPoint = 0.0
         do k = 1,GV%ke
           Top = Bottom
-          MidPoint = Bottom - 0.25*(dz(i,J,k)+dz(i,J-1,k))
-          Bottom = Bottom - 0.5*(dz(i,J,k)+dz(i,J-1,k))
+          MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i,j+1,k))
+          Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i,j+1,k))
           CS%Us_y(i,J,k) = CS%TP_STKY0*exp(MidPoint*DecayScale)
         enddo
       enddo
@@ -837,7 +837,7 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step)
         bottom = 0.0
         do k = 1,GV%ke
           Top = Bottom
-          level_thick = 0.5*(dz(I,j,k)+dz(I-1,j,k))
+          level_thick = 0.5*(dz(i,j,k)+dz(i+1,j,k))
           MidPoint = Top - 0.5*level_thick
           Bottom = Top - level_thick
 
@@ -894,7 +894,7 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step)
         bottom = 0.0
         do k = 1,GV%ke
           Top = Bottom
-          level_thick = 0.5*(dz(i,J,k)+dz(i,J-1,k))
+          level_thick = 0.5*(dz(i,j,k)+dz(i,j+1,k))
           MidPoint = Top - 0.5*level_thick
           Bottom = Top - level_thick
 
@@ -947,8 +947,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step)
           bottom = 0.0
           do k = 1,GV%ke
             Top = Bottom
-            MidPoint = Top - 0.25*(dz(I,j,k)+dz(I-1,j,k))
-            Bottom = Top - 0.5*(dz(I,j,k)+dz(I-1,j,k))
+            MidPoint = Top - 0.25*(dz(i,j,k)+dz(i+1,j,k))
+            Bottom = Top - 0.5*(dz(i,j,k)+dz(i+1,j,k))
             !bgr note that this is using a u-point I on h-point ustar
             !    this code has only been previous used for uniform
             !    grid cases.  This needs fixed if DHH85 is used for non
@@ -964,8 +964,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, dz, ustar, dt, dynamics_step)
           Bottom = 0.0
           do k = 1,GV%ke
             Top = Bottom
-            MidPoint = Bottom - 0.25*(dz(i,J,k)+dz(i,J-1,k))
-            Bottom = Bottom - 0.5*(dz(i,J,k)+dz(i,J-1,k))
+            MidPoint = Bottom - 0.25*(dz(i,j,k)+dz(i,j+1,k))
+            Bottom = Bottom - 0.5*(dz(i,j,k)+dz(i,j+1,k))
             !bgr note that this is using a v-point J on h-point ustar
             !    this code has only been previous used for uniform
             !    grid cases.  This needs fixed if DHH85 is used for non
@@ -1688,8 +1688,8 @@ subroutine CoriolisStokes(G, GV, dt, h, u, v, Waves)
   do k = 1, GV%ke
     do j = G%jsc, G%jec
       do I = G%iscB, G%iecB
-        DVel = 0.25*((Waves%us_y(i,J+1,k)+Waves%us_y(i-1,J+1,k)) * G%CoriolisBu(I,J+1)) + &
-               0.25*((Waves%us_y(i,J,k)+Waves%us_y(i-1,J,k)) * G%CoriolisBu(I,J))
+        DVel = 0.25*((Waves%us_y(i,J-1,k)+Waves%us_y(i+1,J-1,k)) * G%CoriolisBu(I,J-1)) + &
+               0.25*((Waves%us_y(i,J,k)+Waves%us_y(i+1,J,k)) * G%CoriolisBu(I,J))
         u(I,j,k) = u(I,j,k) + DVEL*dt
       enddo
     enddo
@@ -1698,8 +1698,8 @@ subroutine CoriolisStokes(G, GV, dt, h, u, v, Waves)
   do k = 1, GV%ke
     do J = G%jscB, G%jecB
       do i = G%isc, G%iec
-        DVel = 0.25*((Waves%us_x(I+1,j,k)+Waves%us_x(I+1,j-1,k)) * G%CoriolisBu(I+1,J)) + &
-               0.25*((Waves%us_x(I,j,k)+Waves%us_x(I,j-1,k)) * G%CoriolisBu(I,J))
+        DVel = 0.25*((Waves%us_x(I-1,j,k)+Waves%us_x(I-1,j+1,k)) * G%CoriolisBu(I-1,j)) + &
+               0.25*((Waves%us_x(I,j,k)+Waves%us_x(I,j+1,k)) * G%CoriolisBu(I,J))
         v(i,J,k) = v(i,j,k) - DVEL*dt
       enddo
     enddo