GridSorting.f90

Go to the documentation of this file.

module gridsorting
  use kindmodule, only: i4b, dp, lgp
  use constantsmodule, only: dhalf
  use cellwithnbrsmodule, only: globalcelltype
  use mathutilmodule, only: is_close
  use basedismodule, only: dis_transform_xy
  implicit none
  private
 
  public :: quicksortgrid
 
contains
  ! Sort an array of integers
  subroutine quicksortgrid(array, arraySize, idxToGlobal)
    integer, intent(inout), dimension(:) :: array
    integer, intent(in) :: arraysize
    type(globalcelltype), dimension(:), pointer :: idxtoglobal
    ! local
    integer :: qsort_threshold = 8
    include "qsort_inline.inc"
 
  contains
    subroutine init()
    end subroutine init
 
    ! Compare two grid cells, this doesn't work as
    ! smooth for staggered discretizations though...
    function lessthan(n, m) result(isLess)
      integer(I4B), intent(in) :: n
      integer(I4B), intent(in) :: m
      logical(LGP) :: isless
      ! local
      type(globalcelltype), pointer :: gcn, gcm
      real(dp) :: xn, yn, zn, xm, ym, zm
      real(dp), dimension(:), pointer, contiguous :: dis_top_n, dis_bot_n, &
                                                     dis_top_m, dis_bot_m
      real(dp), dimension(:), pointer, contiguous :: dis_xc_n, dis_yc_n, &
                                                     dis_xc_m, dis_yc_m
      real(dp) :: xorigin_n, yorigin_n, angrot_n, &
                  xorigin_m, yorigin_m, angrot_m
 
      ! get coordinates
      gcn => idxtoglobal(array(n))
      gcm => idxtoglobal(array(m))
 
      ! get model data
      ! for n:
      dis_top_n => gcn%v_model%dis_top%get_array()
      dis_bot_n => gcn%v_model%dis_bot%get_array()
      dis_xc_n => gcn%v_model%dis_xc%get_array()
      dis_yc_n => gcn%v_model%dis_yc%get_array()
      xorigin_n = gcn%v_model%dis_xorigin%get()
      yorigin_n = gcn%v_model%dis_yorigin%get()
      angrot_n = gcn%v_model%dis_angrot%get()
      ! for m:
      dis_top_m => gcm%v_model%dis_top%get_array()
      dis_bot_m => gcm%v_model%dis_bot%get_array()
      dis_xc_m => gcm%v_model%dis_xc%get_array()
      dis_yc_m => gcm%v_model%dis_yc%get_array()
      xorigin_m = gcm%v_model%dis_xorigin%get()
      yorigin_m = gcm%v_model%dis_yorigin%get()
      angrot_m = gcm%v_model%dis_angrot%get()
 
      ! convert coordinates
      call dis_transform_xy(dis_xc_n(gcn%index), dis_yc_n(gcn%index), &
                            xorigin_n, yorigin_n, angrot_n, &
                            xn, yn)
      zn = dhalf * (dis_top_n(gcn%index) + &
                    dis_bot_n(gcn%index))
 
      call dis_transform_xy(dis_xc_m(gcm%index), dis_yc_m(gcm%index), &
                            xorigin_m, yorigin_m, angrot_m, &
                            xm, ym)
      zm = dhalf * (dis_top_m(gcm%index) + &
                    dis_bot_m(gcm%index))
 
      ! compare
      if (.not. is_close(zn, zm, 10 * epsilon(zn))) then
        isless = zn > zm
      else if (.not. is_close(yn, ym, 10 * epsilon(yn))) then
        isless = yn > ym
      else if (.not. is_close(xn, xm, 10 * epsilon(xn))) then
        isless = xn < xm
      else
        isless = .false.
      end if
 
    end function lessthan
 
    ! swap indices
    subroutine swap(a, b)
      integer, intent(in) :: a, b
      integer :: hold
 
      hold = array(a)
      array(a) = array(b)
      array(b) = hold
 
    end subroutine swap
 
    ! circular shift-right by one
    subroutine rshift(left, right)
      integer, intent(in) :: left, right
      integer :: hold
 
      hold = array(right)
      array(left + 1:right) = array(left:right - 1)
      array(left) = hold
 
    end subroutine rshift
  end subroutine quicksortgrid
end module gridsorting