      module mergesort_class
      use face_abs_class, only : face_abs
      implicit none

      private

      public :: mergesort

      interface mergesort
        module procedure ms_int64, ms_face
      end interface

      contains

      subroutine ms_int64(array,n,idx)
      implicit none

      integer          :: n,currsize,l,r,mid
      integer          :: idx(*)
      integer (kind=8) :: array(*)

      currsize = 1
      do while (currsize .le. n)
        l = 1
        do while (l .le. n)
          mid = min(l + currsize - 1,n)
          r = min(l+2*currsize-1,n)
          ! Merge
          call mergea_int64(array,l,mid,r,idx)
          l = l + 2*currsize
        end do
        currsize = 2*currsize
      end do

      end subroutine ms_int64

      subroutine mergea_int64(array,l,m,r,idx)
      implicit none

      integer :: l,m,r,i,j,k,n1,n2
      integer :: idx(*)
      integer :: lidx(m-l+1),ridx(r-m)
      integer (kind=8) :: array(*)
      integer (kind=8), dimension(:), allocatable :: larr, rarr

      n1 = m-l+1
      n2 = r-m

      ! Allocate and Fill Work Arrays
      allocate(larr(1:n1),source=array(l:m))
      allocate(rarr(1:n2),source=array(m+1:r))

      lidx(:) = idx(l:m)
      ridx(:) = idx(m+1:r)

      i = 1
      j = 1
      k = l
      do while (i .le. n1 .and. j .le. n2)
        if (larr(i) .le. rarr(j)) then
          array(k) = larr(i)
          idx(k)   = lidx(i)
          i = i + 1
        else
          array(k) = rarr(j)
          idx(k)   = ridx(j)
          j = j + 1
        end if
        k = k + 1
      end do
      do while (i .le. n1)
        array(k) = larr(i)
        idx(k)   = lidx(i)
        i = i + 1
        k = k + 1
      end do
      do while (j .le. n2)
        array(k) = rarr(j)
        idx(k)   = ridx(j)
        j = j + 1
        k = k + 1
      end do
      end subroutine mergea_int64

      subroutine ms_face(array,n,idx)
      implicit none

      integer          :: i
      integer          :: n,currsize,l,r,mid
      integer          :: idx(*)
      class(face_abs)  :: array(*)

      currsize = 1
      do while (currsize .le. n)
        l = 1
        do while (l .le. n)
          mid = min(l + currsize - 1,n)
          r = min(l+2*currsize-1,n)
          ! Merge
          call mergea_face(array,l,mid,r,idx)
          l = l + 2*currsize
        end do
        currsize = 2*currsize
      end do

      end subroutine ms_face

      subroutine mergea_face(array,l,m,r,idx)
      implicit none

      integer :: l,m,r,i,j,k,n1,n2
      integer :: idx(*)
      integer :: lidx(m-l+1),ridx(r-m)

      class(face_abs)                            :: array(*)
      class(face_abs), dimension(:), allocatable :: larr, rarr 
 

      n1 = m-l+1
      n2 = r-m

      ! Allocate and Fill Work Arrays
      allocate(larr(1:n1),source=array(l:m))
      allocate(rarr(1:n2),source=array(m+1:r))
      ! Manually Fill Work Arrays to correct 'source' problem
      do i = 1,n1
        larr(i) = array(l+i-1)
      end do
      do i = 1,n2
        rarr(i) = array(m+i)
      end do

      lidx(:) = idx(l:m)
      ridx(:) = idx(m+1:r)

      i = 1
      j = 1
      k = l
      do while (i .le. n1 .and. j .le. n2)
        if (larr(i) .le. rarr(j)) then
          array(k) = larr(i)
          idx(k)   = lidx(i)
          i = i + 1
        else
          array(k) = rarr(j)
          idx(k)   = ridx(j)
          j = j + 1
        end if
        k = k + 1
      end do
      do while (i .le. n1)
        array(k) = larr(i)
        idx(k)   = lidx(i)
        i = i + 1
        k = k + 1
      end do
      do while (j .le. n2)
        array(k) = rarr(j)
        idx(k)   = ridx(j)
        j = j + 1
        k = k + 1
      end do
      end subroutine mergea_face

      end module mergesort_class
