diff --git a/code/math.f90 b/code/math.f90
index 4fbc4fd38..ac3f54322 100644
--- a/code/math.f90
+++ b/code/math.f90
@@ -144,6 +144,7 @@ real(pReal), dimension(4,36), parameter, private :: &
 
  public :: &
    math_init, &
+   math_qsort, &
    math_range, &
    math_identity2nd, &
    math_identity4th, &
@@ -235,7 +236,6 @@ real(pReal), dimension(4,36), parameter, private :: &
 #endif
             
  private :: &
-   qsort, &
    math_partition, &
    halton, &
    halton_memory, &
@@ -346,7 +346,7 @@ end subroutine math_init
 ! Sorting is done with respect to array(1,:)
 ! and keeps array(2:N,:) linked to it.
 !--------------------------------------------------------------------------------------------------
-recursive subroutine qsort(a, istart, iend)
+recursive subroutine math_qsort(a, istart, iend)
 
  implicit none
  integer(pInt), dimension(:,:), intent(inout) :: a
@@ -355,11 +355,11 @@ recursive subroutine qsort(a, istart, iend)
 
  if (istart < iend) then
    ipivot = math_partition(a,istart, iend)
-   call qsort(a, istart, ipivot-1_pInt)
-   call qsort(a, ipivot+1_pInt, iend)
+   call math_qsort(a, istart, ipivot-1_pInt)
+   call math_qsort(a, ipivot+1_pInt, iend)
  endif
 
-end subroutine qsort
+end subroutine math_qsort
 
 
 !--------------------------------------------------------------------------------------------------
diff --git a/code/mesh.f90 b/code/mesh.f90
index 6d5150fcb..15fccacec 100644
--- a/code/mesh.f90
+++ b/code/mesh.f90
@@ -2383,7 +2383,7 @@ subroutine mesh_marc_count_cpElements(myUnit)
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_marc_map_elements(myUnit)
 
- use math, only: qsort
+ use math, only: math_qsort
  use IO,   only: IO_lc, &
                  IO_stringValue, &
                  IO_stringPos, &
@@ -2421,7 +2421,7 @@ subroutine mesh_marc_map_elements(myUnit)
    endif
  enddo
 
-660 call qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))                 ! should be mesh_NcpElems
+660 call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))            ! should be mesh_NcpElems
 
 end subroutine mesh_marc_map_elements
 
@@ -2432,7 +2432,7 @@ end subroutine mesh_marc_map_elements
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_marc_map_nodes(myUnit)
 
- use math, only: qsort
+ use math, only: math_qsort
  use IO,   only: IO_lc, &
                  IO_stringValue, &
                  IO_stringPos, &
@@ -2469,7 +2469,7 @@ subroutine mesh_marc_map_nodes(myUnit)
    endif
  enddo
 
-650 call qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
+650 call math_qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
  
 end subroutine mesh_marc_map_nodes
 
@@ -2995,7 +2995,7 @@ end subroutine mesh_abaqus_count_cpElements
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_abaqus_map_elements(myUnit)
 
- use math, only: qsort
+ use math, only: math_qsort
  use IO,   only: IO_lc, &
                  IO_stringValue, &
                  IO_stringPos, &
@@ -3045,7 +3045,7 @@ subroutine mesh_abaqus_map_elements(myUnit)
    endselect
  enddo
 
-660 call qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))                  ! should be mesh_NcpElems
+660 call math_qsort(mesh_mapFEtoCPelem,1_pInt,int(size(mesh_mapFEtoCPelem,2_pInt),pInt))                  ! should be mesh_NcpElems
 
  if (int(size(mesh_mapFEtoCPelem),pInt) < 2_pInt) call IO_error(error_ID=907_pInt)
 
@@ -3058,7 +3058,7 @@ end subroutine mesh_abaqus_map_elements
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_abaqus_map_nodes(myUnit)
 
- use math, only: qsort
+ use math, only: math_qsort
  use IO,   only: IO_lc, &
                  IO_stringValue, &
                  IO_stringPos, &
@@ -3109,7 +3109,7 @@ subroutine mesh_abaqus_map_nodes(myUnit)
    endif
  enddo
 
-650 call qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
+650 call math_qsort(mesh_mapFEtoCPnode,1_pInt,int(size(mesh_mapFEtoCPnode,2_pInt),pInt))
 
  if (int(size(mesh_mapFEtoCPnode),pInt) == 0_pInt) call IO_error(error_ID=908_pInt)