subroutine tetrahedron_sort(a, b1, b2, t)
!===========================================!
!A simple bubble sort subroutine for size 4 !
!Assuming the size of a is four !
!===========================================!
real(kind=dp), DIMENSION(4), INTENT(INOUT) :: a, b1, b2
real(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: t
real(kind=dp), DIMENSION(4) :: b1_temp, b2_temp
real(kind=dp), DIMENSION(3, 4) :: t_temp
INTEGER :: i, j
INTEGER, DIMENSION(4) :: reference
real(kind=dp) :: temp
integer :: temp2
DO i = 1, 4
reference(i) = i
END DO
b1_temp = b1; b2_temp = b2
DO j = 1, 3
t_temp(j, 1) = 0
END DO
DO i = 2, 4
DO j = 1, 3
t_temp(j, i) = t(j, i - 1)
END DO
END DO
DO i = 4, 1, -1
DO j = 1, i - 1, +1
IF (a(j) > a(j + 1)) THEN
temp = a(j)
a(j) = a(j + 1)
a(j + 1) = temp
temp2 = reference(j)
reference(j) = reference(j + 1)
reference(j + 1) = temp2
END IF
END DO
END DO
! rearrange t(j,i)
DO i = 2, 4
DO j = 1, 3
t(j, i - 1) = t_temp(j, reference(i)) - t_temp(j, reference(1))
END DO
END DO
! rearrange b(j,i)
DO i = 1, 4
b1(i) = b1_temp(reference(i))
b2(i) = b2_temp(reference(i))
END DO
end subroutine tetrahedron_sort
