      program testsort
*****************************************************************
*                                                               *
*     Copyright 2011 Sourendu Gupta.                            *
*     This program is distributed without warranty under the    *
*     GNU General Public Licence.                               *
*                                                               *
*****************************************************************
      integer, parameter :: mboot=100000
      real, dimension(mboot):: array, scratch
*
      call random_seed (size=n)
*
      do call=1,10000
         call random_number(array)
         call mergesort(1,mboot)
      end do
*
 6001 format(e16.8)
*
      contains
      recursive subroutine mergesort(lbgn,lend)
*****************************************************************
*                                                               *
*     Divide and conquer algorithm for sorting an array         *
*                                                               *
*****************************************************************
      implicit none
      integer, parameter :: magic=4             ! magic number
      integer, intent (in) :: lbgn, lend
      integer :: lmid, i, j, now1, now2, ndum
*****************************************************************
*                                                               *
*     The algorithm specification                               *
*                                                               *
*****************************************************************
      lmid=lend-lbgn+1
      if (lmid.lt.magic) then
*****************************************************************
*                                                               *
*        Array size is smaller than the magic number: simple    *
*                                                               *
*****************************************************************
         if (lmid.eq.2) then            ! exchange if necessary
            if (array(lbgn).gt.array(lend)) then
               call swap(lbgn,lend)
            end if
         else if (lmid.gt.2) then       ! bubble sort small array
            do i=lbgn,lend-1
               do j=i+1,lend
                  if (array(i).gt.array(j)) then
                     call swap(i,j)
                  end if
               end do
            end do
         end if                         ! do nothing if single element
      else      ! large array, first divide and then join
*****************************************************************
*                                                               *
*        Array size is larger than the magic number:            *
*        Divide into smaller pieces and then join the results   *
*                                                               *
*****************************************************************
         lmid=lbgn+lmid/2
*****************************************************************
*        Divide into smaller pieces                             *
*****************************************************************
         call mergesort(lbgn,lmid)      ! divide the problem ...
         call mergesort(lmid+1,lend)    ! ... into two pieces
*****************************************************************
*        Join two sorted lists                                  *
*****************************************************************
         now1=lbgn
         now2=lmid+1
         ndum=lbgn
 10      if ((now1.gt.lmid).and.(now2.le.lend)) then
c           Trailing part of list 2 is already in the correct place
c           Just copy the merged list into position before this
            array(lbgn:ndum-1)=scratch(lbgn:ndum-1)
         else if((now2.gt.lend).and.(now1.le.lmid)) then        
            do i=now1,lmid              ! append to end of merged list
               scratch(ndum)=array(i)
               ndum=ndum+1
            end do
            array(lbgn:lend)=scratch(lbgn:lend)
         else if (array(now1).gt.array(now2)) then
            scratch(ndum)=array(now2)   ! write the smaller element
            now2=now2+1
            ndum=ndum+1
            go to 10
         else 
            scratch(ndum)=array(now1)   ! write the smaller element
            now1=now1+1
            ndum=ndum+1
            go to 10
         end if 
      end if
*
      end subroutine mergesort
      subroutine swap(i,j)
*****************************************************************
*                                                               *
*     Swap two elements of the array                            *
*                                                               *
*****************************************************************
      implicit none
      integer, intent (in) :: i, j
      real :: tmp
*
      tmp = array(i)
      array(i) = array(j)
      array(j) = tmp
*
      end subroutine swap
      end program testsort
