diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90
index 74a0305ce1..69d835ae95 100644
--- a/src/parameterizations/vertical/MOM_CVMix_conv.F90
+++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90
@@ -174,7 +174,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)
   integer :: kOBL                      !< level of ocean boundary layer extent
   real :: g_o_rho0  ! Gravitational acceleration, perhaps divided by density, times unit conversion factors
                     ! [H s-2 R-1 ~> m4 s-2 kg-1 or m s-2]
-  real :: pref      ! Interface pressures [R L2 T-2 ~> Pa]
+  real :: pRef      ! Interface pressures [R L2 T-2 ~> Pa]
   real :: rhok, rhokm1 ! In situ densities of the layers above and below at the interface pressure [R ~> kg m-3]
   real :: dh_int    ! The distance between layer centers [H ~> m or kg m-2]
   real :: dh, hcorr ! Limited thicknesses and a cumulative correction [Z ~> m]
@@ -211,7 +211,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)
       do K=2,GV%ke
 
         ! pRef is pressure at interface between k and km1 [R L2 T-2 ~> Pa].
-        pRef = pRef + (GV%H_to_RZ*GV%g_Earth) * h(i,j,k)
+        pRef = pRef + (GV%H_to_RZ*GV%g_Earth) * h(i,j,k-1)
         call calculate_density(tv%t(i,j,k), tv%s(i,j,k), pRef, rhok, tv%eqn_of_state)
         call calculate_density(tv%t(i,j,k-1), tv%s(i,j,k-1), pRef, rhokm1, tv%eqn_of_state)