diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90
index 9bc71dd15f..99c2ac9799 100644
--- a/src/ALE/MOM_regridding.F90
+++ b/src/ALE/MOM_regridding.F90
@@ -1923,6 +1923,90 @@ subroutine convective_adjustment(G, GV, h, tv)
 
 end subroutine convective_adjustment
 
+!------------------------------------------------------------------------------
+!> Return the index of a sorted array of scalar values
+subroutine sort_scalar_k_1d(G, GV, phi, ksort)
+  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure
+  type(verticalGrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure
+  real, dimension(SZK_(GV)), &
+                           intent(in)    :: phi  !< Array of scalar quantity to be sorted
+  integer, dimension(SZK_(GV)), &
+                           intent(out) :: ksort !< An array of indicies for a 
+                                                  !! monotonically increasing scalar
+!------------------------------------------------------------------------------
+! Check each water column to see if a given scalar is monotonically increasing.
+! If not, return an array of the sorted indices (bubble sort algorithm).
+! No need to return the sorted scalar array itself.
+!------------------------------------------------------------------------------
+
+  ! Local variables
+  integer   :: k
+  real      :: P0, P1       ! temperatures
+  logical   :: monotonic
+
+  ! Repeat swapping of indices until complete
+  do
+    monotonic = .true.
+    do k = 1,GV%ke-1
+      ! Gather information of scalar value in current and next cells
+      P0 = phi(k)  ; P1 = phi(k+1)
+      ! If the scalar value of the current cell is larger than the scalar
+      ! below it, we swap the cell indices
+      if ( P0 > P1 ) then
+        ksort(k) = k+1 ; ksort(k+1) = k
+        monotonic = .false.
+      endif
+    enddo  ! k
+
+    if ( monotonic ) exit
+  enddo
+
+end subroutine sort_scalar_k_1d
+
+!------------------------------------------------------------------------------
+!> Return the index of a sorted array of scalar values
+subroutine sort_scalar_k(G, GV, phi, ksort)
+  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure
+  type(verticalGrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(in)    :: phi  !< Array of scalar quantity to be sorted
+  integer, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(out) :: ksort !< An array of indicies for a 
+                                                  !! monotonically increasing scalar
+!------------------------------------------------------------------------------
+! Check each water column to see if a given scalar is monotonically increasing.
+! If not, return an array of the sorted indices (bubble sort algorithm).
+! No need to return the sorted scalar array itself.
+!------------------------------------------------------------------------------
+
+  ! Local variables
+  integer   :: i, j, k
+  real      :: P0, P1       ! temperatures
+  logical   :: monotonic
+
+  ! Loop on columns
+  do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
+
+    ! Repeat swapping of indices until complete
+    do
+      monotonic = .true.
+      do k = 1,GV%ke-1
+        ! Gather information of scalar value in current and next cells
+        P0 = phi(i,j,k)  ; P1 = phi(i,j,k+1)
+        ! If the scalar value of the current cell is larger than the scalar
+        ! below it, we swap the cell indices
+        if ( P0 > P1 ) then
+          ksort(i,j,k) = k+1 ; ksort(i,j,k+1) = k
+          monotonic = .false.
+        endif
+      enddo  ! k
+
+      if ( monotonic ) exit
+    enddo
+
+  enddo ; enddo  ! i & j
+
+end subroutine sort_scalar_k
 
 !------------------------------------------------------------------------------
 !> Return a uniform resolution vector in the units of the coordinate