subroutine SortQuick( im, jm, km, arr, arr1, arr2 )

    integer , intent(in   )           :: im
    integer , intent(in   )           :: jm
    integer , intent(in   )           :: km
    real(DP), intent(inout)           :: arr ( im, jm, km )
    real(DP), intent(inout), optional :: arr1( im, jm, km )
    real(DP), intent(inout), optional :: arr2( im, jm, km )


    !
    ! local varialbes
    !
    integer :: i
    integer :: j
    integer :: ks, ke


    do j = 1, jm
      do i = 1, im

        ks = 1
        ke = km
        call sort_quick0( im, jm, km, i, j, ks, ke, arr, arr1, arr2 )

      end do
    end do


  end subroutine SortQuick

  recursive subroutine sort_quick0( im, jm, km, i, j, ks, ke, arr, arr1, arr2 )

    integer , intent(in   )           :: im
    integer , intent(in   )           :: jm
    integer , intent(in   )           :: km
    integer , intent(in   )           :: i
    integer , intent(in   )           :: j
    integer , intent(in   )           :: ks, ke
    real(DP), intent(inout)           :: arr ( im, jm, km )
    real(DP), intent(inout), optional :: arr1( im, jm, km )
    real(DP), intent(inout), optional :: arr2( im, jm, km )


    !
    ! local varialbes
    !
    real(DP) :: pivot
    integer  :: k1, k2, ks1, ke1, ks2, ke2

!!$        logical      :: sw_find1, sw_find2


    if( ks == ke ) return


    !
    ! select pivot
    !
    do k1 = ks+1, ke
      if( arr( i, j, ks ) .ne. arr( i, j, k1 ) ) exit
    end do
    if( k1 > ke ) return
    if( arr( i, j, ks ) .ge. arr( i, j, k1 ) ) then
      pivot = arr( i, j, ks )
    else
      pivot = arr( i, j, k1 )
    end if


    k1 = ks
    k2 = ke

    do

      do
        if( arr( i, j, k1 ) >= pivot ) exit
        k1 = k1 + 1
        ! MEMO:
        ! I never fail to find a value greater than or equal to pivot.
!!$              if( i1 .ge. ie ) exit
      end do
      do
        if( arr( i, j, k2 ) < pivot ) exit
        k2 = k2 - 1
        ! MEMO:
        ! I never fail to find a value less than or equal to pivot.
        ! This is due to a way to determine pivot.
        ! A value different from pivot MUST be included in array.
!!$              if( i2 .le. is ) exit
      end do

      ! MEMO:
      ! The values i1 and i2 must not be the same. 
      ! This is due to a way to determine pivot.
      ! A value different from pivot MUST be included in array.
!!$           if( i1 .ge. i2 ) then
      if( k1 .gt. k2 ) then

        ! MEMO:
        ! The value arr( i1 ) must not be greater than or equal to pivot.
        ! This is due to a way to determine pivot.
        ! A value different from pivot MUST be included in array.
!!$              if( arr( i1 ) .ge. pivot ) then
!!$                 ie1 = i1 - 1
!!$              else
!!$                 ie1 = i1
!!$              end if

        ke1 = k1 - 1

        exit
      end if

      call swap( im, jm, km, arr, i, j, k1, k2 )
      if( present( arr1 ) ) call swap( im, jm, km, arr1, i, j, k1, k2 )
      if( present( arr2 ) ) call swap( im, jm, km, arr2, i, j, k1, k2 )


    end do


    ks1 = ks
!!$        ie1 = i1 - 1
    call sort_quick0( im, jm, km, i, j, ks1, ke1, arr, arr1, arr2 )
    ks2 = ke1 + 1
    ke2 = ke
    call sort_quick0( im, jm, km, i, j, ks2, ke2, arr, arr1, arr2 )


  end subroutine sort_quick0

  subroutine swap( im, jm, km, arr, i, j, k1, k2 )

    integer , intent(in   ) :: im
    integer , intent(in   ) :: jm
    integer , intent(in   ) :: km
    real(DP), intent(inout) :: arr( im, jm, km )
    integer , intent(in   ) :: i
    integer , intent(in   ) :: j
    integer , intent(in   ) :: k1, k2


    !
    ! local variables
    !
    real(DP) :: rarr


    rarr            = arr( i, j, k1 )
    arr( i, j, k1 ) = arr( i, j, k2 )
    arr( i, j, k2 ) = rarr


  end subroutine swap