sort_ascending Subroutine

   subroutine sort_ascending(n, dx, dy)

      integer, intent(in) :: n
      real(wp), dimension(*), intent(inout) :: dx !! array of values to be sorted
      real(wp), dimension(*), intent(inout), optional :: dy !! array to be (optionally) carried along

      logical :: carry_dy !! if `dy` is to be also sorted

      integer, parameter :: max_size_for_insertion_sort = 20 !! max size for using insertion sort.
                                                             !! (otherwise, use quicksort)


      carry_dy = present(dy)
      call quicksort(1, n)

   contains

      recursive subroutine quicksort(ilow, ihigh)

         !! Sort the array (ascending order).

         integer, intent(in) :: ilow
         integer, intent(in) :: ihigh

         integer :: ipivot !! pivot element
         integer :: i      !! counter
         integer :: j      !! counter

         if (ihigh - ilow <= max_size_for_insertion_sort .and. ihigh > ilow) then

            ! do insertion sort:
            do i = ilow + 1, ihigh
               do j = i, ilow + 1, -1
                  if (dx(j) < dx(j - 1)) then
                     call swap(dx(j), dx(j - 1))
                     if (carry_dy) call swap(dy(j), dy(j - 1))
                  else
                     exit
                  end if
               end do
            end do

         else if (ihigh - ilow > max_size_for_insertion_sort) then

            ! do the normal quicksort:
            call partition(ilow, ihigh, ipivot)
            call quicksort(ilow, ipivot - 1)
            call quicksort(ipivot + 1, ihigh)

         end if

      end subroutine quicksort

      subroutine partition(ilow, ihigh, ipivot)

         !! Partition the array

         integer, intent(in)  :: ilow
         integer, intent(in)  :: ihigh
         integer, intent(out) :: ipivot

         integer :: i, ip, im

         im = (ilow + ihigh)/2
         call swap(dx(ilow), dx(im))
         if (carry_dy) call swap(dy(ilow), dy(im))
         ip = ilow
         do i = ilow + 1, ihigh
            if (dx(i) < dx(ilow)) then
               ip = ip + 1
               call swap(dx(ip), dx(i))
               if (carry_dy) call swap(dy(ip), dy(i))
            end if
         end do
         call swap(dx(ilow), dx(ip))
         if (carry_dy) call swap(dy(ilow), dy(ip))
         ipivot = ip

      end subroutine partition

      subroutine swap(v1, v2)
         !! swap two real values
         real(wp), intent(inout) :: v1
         real(wp), intent(inout) :: v2
         real(wp) :: tmp
         tmp = v1
         v1 = v2
         v2 = tmp
      end subroutine swap

   end subroutine sort_ascending