MAPL_Sort.F90 Source File

!------------------------------------------------------------------------------
!               Global Modeling and Assimilation Office (GMAO)                !
!                    Goddard Earth Observing System (GEOS)                    !
!                                 MAPL Component                              !
!------------------------------------------------------------------------------
!
#include "MAPL_ErrLog.h"
!
!>
!### MODULE: `MAPL_SortMod`
!
! Author: GMAO SI-Team
!
! `MAPL_SortMod` -- A utility to sort integers
! 
! `MAPL_Sort` is a utility to do a quicksort on integers. 
! The general interface is:
!
!```fortran       
!       subroutine MAPL_Sort(A[,B])
!         integer(kind=[4,8]),           intent(INOUT) :: A(:)
!         integer(kind=4),     optional, intent(INOUT) :: B(size(A))
!         real   (kind=[4,8]), optional, intent(INOUT) :: B(size(A))
!
!       subroutine MAPL_Sort(A[,DIM])
!         integer(kind=[4,8]),           intent(INOUT) :: A(:,:)
!         integer(kind=4),     optional, intent(IN   ) :: DIM
!
!       subroutine MAPL_Sort(A,B[,DIM])
!         integer(kind=[4,8]),       intent(INOUT) :: A(:)
!         integer(kind=4),           intent(INOUT) :: B(:,:)
!         integer(kind=4), optional, intent(IN   ) :: DIM
!```
!
! `MAPL_Sort` sorts the key (contained in a row or column of A)
! in ascending order and reorders the data in B or in non-key rows or columns of A
! in the same order; i.e., it does the same exchanges as were done 
! to the key in sorting it.  If, for example, on input B(:) contains the ordered integers
! from 1 to size(A), on output it will contain the positions of the elements of
! the sorted A in the unsorted A. In the last two signatures, DIM is the dimension
! of A or B being reordered. In the last signature, for example, DIM=1 corresponds
! to a B ordered as B(size(A),:), whereas DIM=2 corresponds to B(:,size(A)).
! The default is DIM=2. The quicksort is coded in C and does not appear here.
!
module MAPL_SortMod
   use, intrinsic :: ISO_FORTRAN_ENV, only: REAL32, REAL64
   use, intrinsic :: ISO_FORTRAN_ENV, only: INT32, INT64

   use MAPL_ExceptionHandling

  implicit none
  private

! !PUBLIC ROUTINES:

  public MAPL_Sort


!=============================================================================

interface MAPL_Sort
   module procedure SORT1SS
   module procedure SORT1SR
   module procedure SORT1SD
   module procedure SORT1LS
   module procedure SORT1LR
   module procedure SORT1LD

   module procedure SORT2LS
   module procedure SORT2LR
   module procedure SORT2LD
   module procedure SORT2SS
   module procedure SORT2SR
   module procedure SORT2SD

   module procedure SORT2AS
   module procedure SORT2AL
end interface MAPL_Sort

interface 
   subroutine qsorts(a, b, r, n) bind(C,name="QSORTS")
      use, intrinsic :: iso_fortran_env, only: INT32
      integer(kind=INT32), intent(inout) :: a(*)
      type(*), intent(inout) :: b(*)
      integer, value, intent(in) :: r
      integer, value, intent(in) :: n
   end subroutine qsorts
   subroutine qsortl(a, b, r, n) bind(C,name="QSORTL")
      use, intrinsic :: iso_fortran_env, only: INT64
      integer(kind=INT64), intent(inout) :: a(*)
      type(*), intent(inout) :: b(*)
      integer, value, intent(in) :: r
      integer, value, intent(in) :: n
   end subroutine qsortl
end interface

contains

subroutine SORT1SS(A,B)
  integer(kind=INT32),           intent(INOUT) :: A(:)
  integer(kind=INT32), optional, intent(INOUT) :: B(:)
  if(present(B)) then
     call QSORTS(A,B,size(A),1)
  else
     call QSORTS(A,A,size(A),0)
  endif
end subroutine SORT1SS

subroutine SORT1SR(A,B)
  integer(kind=INT32),           intent(INOUT) :: A(:)
  real   (kind=REAL32),           intent(INOUT) :: B(:)
  call QSORTS(A,B,size(A),1)
end subroutine SORT1SR

subroutine SORT1SD(A,B)
  integer(kind=INT32),           intent(INOUT) :: A(:)
  real   (kind=REAL64),           intent(INOUT) :: B(:)
  call QSORTS(A,B,size(A),2)
end subroutine SORT1SD

subroutine SORT1LS(A,B)
  integer(kind=INT64), intent(INOUT) :: A(:)
  integer(kind=INT32), optional, intent(INOUT) :: B(:)
  if(present(B)) then
     call QSORTL(A,B,size(A),1)
  else
     call QSORTL(A,A,size(A),0)
  endif
end subroutine SORT1LS

subroutine SORT1LR(A,B)
  integer(kind=INT64),           intent(INOUT) :: A(:)
  real   (kind=REAL32),           intent(INOUT) :: B(:)
  call QSORTL(A,B,size(A),1)
end subroutine SORT1LR

subroutine SORT1LD(A,B)
  integer(kind=INT64),           intent(INOUT) :: A(:)
  real   (kind=REAL64),           intent(INOUT) :: B(:)
  call QSORTL(A,B,size(A),2)
end subroutine SORT1LD


subroutine SORT2SS(A,B,DIM)
  integer(kind=INT32),   intent(INOUT) :: A(:)
  integer(kind=INT32),   intent(INOUT) :: B(:,:)
  integer, optional, intent(IN   ) :: DIM

  integer :: uDIM, rc

  if(present(DIM)) then
     uDIM = DIM
  else
     uDIM = 2
  end if
  _ASSERT(uDIM>0, 'uDim needs to be greater than 0')
  _ASSERT(uDIM<3, 'uDim needs to be less than 3')
  _ASSERT(size(A)==size(B,uDIM),'needs informative message')
  if(uDIM==1) then
     call QSORTS(A,B,size(A),-size(B,2))
  else
     call QSORTS(A,B,size(A), size(B,1))
  end if

end subroutine SORT2SS

subroutine SORT2SR(A,B,DIM)
  integer(kind=INT32),   intent(INOUT) :: A(:)
  real   (kind=REAL32),   intent(INOUT) :: B(:,:)
  integer, optional, intent(IN   ) :: DIM

  integer :: uDIM, rc

  if(present(DIM)) then
     uDIM = DIM
  else
     uDIM = 2
  end if
  _ASSERT(uDIM>0, 'uDim needs to be greater than 0')
  _ASSERT(uDIM<3, 'uDim needs to be less than 3')  
  _ASSERT(size(A)==size(B,uDIM),'needs informative message')
  if(uDIM==1) then
     call QSORTS(A,B,size(A),-size(B,2))
  else
     call QSORTS(A,B,size(A), size(B,1))
  end if

end subroutine SORT2SR

subroutine SORT2SD(A,B,DIM)
  integer(kind=INT32),   intent(INOUT) :: A(:)
  real   (kind=REAL64),   intent(INOUT) :: B(:,:)
  integer, optional, intent(IN   ) :: DIM

  integer :: uDIM, rc

  if(present(DIM)) then
     uDIM = DIM
  else
     uDIM = 2
  end if
  _ASSERT(uDIM>0, 'uDim needs to be greater than 0')
  _ASSERT(uDIM<3, 'uDim needs to be less than 3')
  _ASSERT(size(A)==size(B,uDIM),'needs informative message')
  if(uDIM==1) then
     call QSORTS(A,B,size(A),-size(B,2)*2)
  else
     call QSORTS(A,B,size(A), size(B,1)*2)
  end if

end subroutine SORT2SD

subroutine SORT2LS(A,B,DIM)
  integer(kind=INT64),   intent(INOUT) :: A(:)
  integer(kind=INT32),   intent(INOUT) :: B(:,:)
  integer, optional, intent(IN   ) :: DIM

  integer :: uDIM, rc

  if(present(DIM)) then
     uDIM = DIM
  else
     uDIM = 2
  end if
  _ASSERT(uDIM>0, 'uDim needs to be greater than 0')
  _ASSERT(uDIM<3, 'uDim needs to be less than 3')  
  _ASSERT(size(A)==size(B,uDIM),'needs informative message')
  if(uDIM==1) then
     call QSORTL(A,B,size(A),-size(B,2))
  else
     call QSORTL(A,B,size(A), size(B,1))
  end if
end subroutine SORT2LS

subroutine SORT2LR(A,B,DIM)
  integer(kind=INT64),   intent(INOUT) :: A(:)
  real   (kind=REAL32),   intent(INOUT) :: B(:,:)
  integer, optional, intent(IN   ) :: DIM

  integer :: uDIM, rc

  if(present(DIM)) then
     uDIM = DIM
  else
     uDIM = 2
  end if
  _ASSERT(uDIM>0, 'uDim needs to be greater than 0')
  _ASSERT(uDIM<3, 'uDim needs to be less than 3')  
  _ASSERT(size(A)==size(B,uDIM),'needs informative message')
  if(uDIM==1) then
     call QSORTL(A,B,size(A),-size(B,2))
  else
     call QSORTL(A,B,size(A), size(B,1))
  end if
end subroutine SORT2LR

subroutine SORT2LD(A,B,DIM)
  integer(kind=INT64),   intent(INOUT) :: A(:)
  real   (kind=REAL64),   intent(INOUT) :: B(:,:)
  integer, optional, intent(IN   ) :: DIM

  integer :: uDIM, rc

  if(present(DIM)) then
     uDIM = DIM
  else
     uDIM = 2
  end if
  _ASSERT(uDIM>0, 'uDim needs to be greater than 0')
  _ASSERT(uDIM<3, 'uDim needs to be less than 3')  
  _ASSERT(size(A)==size(B,uDIM),'needs informative message')
  if(uDIM==1) then
     call QSORTL(A,B,size(A),-size(B,2)*2)
  else
     call QSORTL(A,B,size(A), size(B,1)*2)
  end if
end subroutine SORT2LD

subroutine SORT2AS(B,DIM)
  integer(kind=INT32),   intent(INOUT) :: B(:,:)
  integer, optional, intent(IN   ) :: DIM

  integer :: uDIM, rc

  if(present(DIM)) then
     uDIM = DIM
  else
     uDIM = 2
  end if

  _ASSERT(uDIM>0, 'uDim needs to be greater than 0')
  _ASSERT(uDIM<3, 'uDim needs to be less than 3')

  if(uDIM==1) then
     call QSORTS(B(:,1),B(:,2:),size(B,1),-(size(B,2)-1))
  else
     call QSORTS(B(1,:),B(2:,:),size(B,2), (size(B,1)-1))
  end if
end subroutine SORT2AS

subroutine SORT2AL(B,DIM)
  integer(kind=INT64),   intent(INOUT) :: B(:,:)
  integer, optional, intent(IN   ) :: DIM

  integer :: uDIM, rc

  if(present(DIM)) then
     uDIM = DIM
  else
     uDIM = 2
  end if

  _ASSERT(uDIM>0, 'uDim needs to be greater than 0')
  _ASSERT(uDIM<3, 'uDim needs to be less than 3')

  if(uDIM==1) then
     call QSORTL(B(:,1),B(:,2:),size(B,1),-(size(B,2)-1)*2)
  else
     call QSORTL(B(1,:),B(2:,:),size(B,2), (size(B,1)-1)*2)
  end if
end subroutine SORT2AL

end module MAPL_SortMod