diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90
index 20fcb3bd5..8acbc4315 100644
--- a/src/encounter/encounter_check.f90
+++ b/src/encounter/encounter_check.f90
@@ -865,7 +865,7 @@ subroutine encounter_check_remove_duplicates(n, nenc, index1, index2, lvdotr)
    end subroutine encounter_check_remove_duplicates
 
 
-   module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr)
+   pure module subroutine encounter_check_sort_aabb_1D(self, n, extent_arr)
       !! author: David A. Minton
       !!
       !! Sorts the bounding box extents along a single dimension prior to the sweep phase. 
diff --git a/src/gr/gr.f90 b/src/gr/gr.f90
index b05465e98..73a182729 100644
--- a/src/gr/gr.f90
+++ b/src/gr/gr.f90
@@ -2,7 +2,7 @@
    use swiftest
 contains
 
-   module pure subroutine gr_kick_getaccb_ns_body(self, system, param) 
+   pure module subroutine gr_kick_getaccb_ns_body(self, system, param) 
       !! author: David A. Minton
       !!
       !! Add relativistic correction acceleration for non-symplectic integrators.
@@ -45,7 +45,7 @@ module pure subroutine gr_kick_getaccb_ns_body(self, system, param)
    end subroutine gr_kick_getaccb_ns_body
 
 
-   module pure subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) 
+   pure module subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) 
       !! author: David A. Minton
       !!
       !! Compute relativisitic accelerations of massive bodies
@@ -75,7 +75,7 @@ module pure subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr)
    end subroutine gr_kick_getacch
 
 
-   module pure subroutine gr_p4_pos_kick(param, x, v, dt)
+   pure module subroutine gr_p4_pos_kick(param, x, v, dt)
       !! author: David A. Minton
       !!
       !! Position kick due to p**4 term in the post-Newtonian correction
@@ -104,7 +104,7 @@ module pure subroutine gr_p4_pos_kick(param, x, v, dt)
    end subroutine gr_p4_pos_kick
 
 
-   module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh) 
+   pure module subroutine gr_pseudovel2vel(param, mu, xh, pv, vh) 
       !! author: David A. Minton
       !!
       !! Converts the relativistic pseudovelocity back into a veliocentric velocity
@@ -136,7 +136,7 @@ module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh)
    end subroutine gr_pseudovel2vel
 
 
-   module pure subroutine gr_pv2vh_body(self, param)
+   pure module subroutine gr_pv2vh_body(self, param)
       !! author: David A. Minton
       !!
       !! Wrapper function that converts from pseudovelocity to heliocentric velocity for swiftest bodies
@@ -161,7 +161,7 @@ module pure subroutine gr_pv2vh_body(self, param)
    end subroutine gr_pv2vh_body
 
 
-   module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
+   pure module subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
       !! author: David A. Minton
       !!
       !! Converts the heliocentric velocity into a pseudovelocity with relativistic corrections. 
@@ -238,7 +238,7 @@ module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
    end subroutine gr_vel2pseudovel
 
 
-   module pure subroutine gr_vh2pv_body(self, param)
+   pure module subroutine gr_vh2pv_body(self, param)
       !! author: David A. Minton
       !!
       !! Wrapper function that converts from heliocentric velocity to pseudovelocity for Swiftest bodies
diff --git a/src/helio/helio_gr.f90 b/src/helio/helio_gr.f90
index 1ff8eb5d5..2c99d0016 100644
--- a/src/helio/helio_gr.f90
+++ b/src/helio/helio_gr.f90
@@ -2,7 +2,7 @@
    use swiftest
 contains
 
-   module pure subroutine helio_gr_kick_getacch_pl(self, param) 
+   pure module subroutine helio_gr_kick_getacch_pl(self, param) 
       !! author: David A. Minton
       !!
       !! Compute relativisitic accelerations of massive bodies
@@ -30,7 +30,7 @@ module pure subroutine helio_gr_kick_getacch_pl(self, param)
    end subroutine helio_gr_kick_getacch_pl
 
 
-   module pure subroutine helio_gr_kick_getacch_tp(self, param)
+   pure module subroutine helio_gr_kick_getacch_tp(self, param)
       !! author: David A. Minton
       !!
       !! Compute relativisitic accelerations of test particles
@@ -56,7 +56,7 @@ module pure subroutine helio_gr_kick_getacch_tp(self, param)
    end subroutine helio_gr_kick_getacch_tp
    
 
-   module pure subroutine helio_gr_p4_pl(self, system, param, dt)
+   pure module subroutine helio_gr_p4_pl(self, system, param, dt)
       !! author: David A. Minton
       !!
       !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction
@@ -84,7 +84,7 @@ module pure subroutine helio_gr_p4_pl(self, system, param, dt)
    end subroutine helio_gr_p4_pl
 
 
-   module pure subroutine helio_gr_p4_tp(self, system, param, dt)
+   pure module subroutine helio_gr_p4_tp(self, system, param, dt)
       !! author: David A. Minton
       !!
       !! Position kick to test particles due to p**4 term in the post-Newtonian correction
diff --git a/src/kick/kick.f90 b/src/kick/kick.f90
index e680fde98..a9a5e8d7a 100644
--- a/src/kick/kick.f90
+++ b/src/kick/kick.f90
@@ -255,7 +255,7 @@ module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
    end subroutine kick_getacch_int_all_tp
 
 
-   module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
+   pure module subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
       !! author: David A. Minton
       !!
       !! Compute direct cross (third) term heliocentric accelerations for a single pair of massive bodies
@@ -286,7 +286,7 @@ module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi,
    end subroutine kick_getacch_int_one_pl
 
 
-   module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, az)
+   pure module subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, GMpl, ax, ay, az)
       !! author: David A. Minton
       !!
       !! Compute direct cross (third) term heliocentric accelerations of a single test particle massive body pair.
diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90
index f7758ef54..39b88ae05 100644
--- a/src/modules/encounter_classes.f90
+++ b/src/modules/encounter_classes.f90
@@ -128,7 +128,7 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in
          logical,              dimension(:), allocatable, intent(out), optional :: lvdotr      !! Array indicating which bodies are approaching
       end subroutine encounter_check_collapse_ragged_list
 
-      module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr)
+      pure module subroutine encounter_check_sort_aabb_1D(self, n, extent_arr)
          implicit none
          class(encounter_bounding_box_1D), intent(inout) :: self       !! Bounding box structure along a single dimension
          integer(I4B),                     intent(in)    :: n          !! Number of bodies with extents
diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90
index 78a4bdc34..81d600c48 100644
--- a/src/modules/helio_classes.f90
+++ b/src/modules/helio_classes.f90
@@ -108,21 +108,21 @@ module subroutine helio_drift_linear_tp(self, cb, dt, lbeg)
          logical,                       intent(in)    :: lbeg !! Argument that determines whether or not this is the beginning or end of the step
       end subroutine helio_drift_linear_tp
 
-      module pure subroutine helio_gr_kick_getacch_pl(self, param) 
+      pure module subroutine helio_gr_kick_getacch_pl(self, param) 
          use swiftest_classes, only : swiftest_parameters
          implicit none
          class(helio_pl),            intent(inout) :: self   !! Helio massive body particle data structure
          class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
       end subroutine helio_gr_kick_getacch_pl
    
-      module pure subroutine helio_gr_kick_getacch_tp(self, param)
+      pure module subroutine helio_gr_kick_getacch_tp(self, param)
          use swiftest_classes, only : swiftest_parameters
          implicit none
          class(helio_tp),            intent(inout) :: self   !! Helio massive body particle data structure
          class(swiftest_parameters), intent(in)    :: param  !! Current run configuration parameters 
       end subroutine helio_gr_kick_getacch_tp
       
-      module pure subroutine helio_gr_p4_pl(self, system, param, dt)
+      pure module subroutine helio_gr_p4_pl(self, system, param, dt)
          use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
          implicit none
          class(helio_pl),              intent(inout) :: self   !! Swiftest particle object
@@ -131,7 +131,7 @@ module pure subroutine helio_gr_p4_pl(self, system, param, dt)
          real(DP),                     intent(in)    :: dt     !! Step size
       end subroutine helio_gr_p4_pl
    
-      module pure subroutine helio_gr_p4_tp(self, system, param, dt)
+      pure module subroutine helio_gr_p4_tp(self, system, param, dt)
          use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
          implicit none
          class(helio_tp),              intent(inout) :: self   !! Swiftest particle object
diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90
index ca2d97347..d2f1adeef 100644
--- a/src/modules/swiftest_classes.f90
+++ b/src/modules/swiftest_classes.f90
@@ -518,14 +518,14 @@ module pure elemental subroutine drift_one(mu, px, py, pz, vx, vy, vz, dt, iflag
          integer(I4B), intent(out)      :: iflag !! iflag : error status flag for Danby drift (0 = OK, nonzero = ERROR)
       end subroutine drift_one
 
-      module pure subroutine gr_kick_getaccb_ns_body(self, system, param)
+      pure module subroutine gr_kick_getaccb_ns_body(self, system, param)
          implicit none
          class(swiftest_body),         intent(inout) :: self   !! Swiftest generic body object
          class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object
          class(swiftest_parameters),   intent(in)    :: param  !! Current run configuration parameters 
       end subroutine gr_kick_getaccb_ns_body
 
-      module pure subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) 
+      pure module subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr) 
          implicit none
          real(DP), dimension(:),     intent(in)  :: mu     !! Gravitational constant
          real(DP), dimension(:,:),   intent(in)  :: x      !! Position vectors
@@ -535,7 +535,7 @@ module pure subroutine gr_kick_getacch(mu, x, lmask, n, inv_c2, agr)
          real(DP), dimension(:,:),   intent(out) :: agr    !! Accelerations
       end subroutine gr_kick_getacch
 
-      module pure subroutine gr_p4_pos_kick(param, x, v, dt)
+      pure module subroutine gr_p4_pos_kick(param, x, v, dt)
          implicit none
          class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
          real(DP), dimension(:),     intent(inout) :: x     !! Position vector
@@ -543,7 +543,7 @@ module pure subroutine gr_p4_pos_kick(param, x, v, dt)
          real(DP),                   intent(in)    :: dt    !! Step size
       end subroutine gr_p4_pos_kick
 
-      module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh) 
+      pure module subroutine gr_pseudovel2vel(param, mu, xh, pv, vh) 
          implicit none
          class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
          real(DP),                   intent(in)  :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
@@ -552,13 +552,13 @@ module pure subroutine gr_pseudovel2vel(param, mu, xh, pv, vh)
          real(DP), dimension(:),     intent(out) :: vh    !! Swiftestcentric velocity vector 
       end subroutine gr_pseudovel2vel
 
-      module pure subroutine gr_pv2vh_body(self, param)
+      pure module subroutine gr_pv2vh_body(self, param)
          implicit none
          class(swiftest_body),       intent(inout) :: self  !! Swiftest particle object
          class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
       end subroutine gr_pv2vh_body
 
-      module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
+      pure module subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
          implicit none
          class(swiftest_parameters), intent(in)  :: param !! Current run configuration parameters 
          real(DP),                   intent(in)  :: mu    !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body
@@ -567,7 +567,7 @@ module pure subroutine gr_vel2pseudovel(param, mu, xh, vh, pv)
          real(DP), dimension(:),     intent(out) :: pv    !! Pseudovelocity vector - see Saha & Tremain (1994), eq. (32)
       end subroutine gr_vel2pseudovel
 
-      module pure subroutine gr_vh2pv_body(self, param)
+      pure module subroutine gr_vh2pv_body(self, param)
          implicit none
          class(swiftest_body),       intent(inout) :: self  !! Swiftest particle object
          class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
@@ -869,7 +869,7 @@ module subroutine kick_getacch_int_all_tp(ntp, npl, xtp, xpl, GMpl, lmask, acc)
          real(DP),     dimension(:,:), intent(inout) :: acc    !! Acceleration vector array 
       end subroutine kick_getacch_int_all_tp
 
-      module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
+      pure module subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi, ayi, azi, axj, ayj, azj)
          !$omp declare simd(kick_getacch_int_one_pl)
          implicit none
          real(DP), intent(in)  :: rji2            !! Square of distance between the two bodies
@@ -880,7 +880,7 @@ module pure subroutine kick_getacch_int_one_pl(rji2, xr, yr, zr, Gmi, Gmj, axi,
          real(DP), intent(inout) :: axj, ayj, azj !! Acceleration vector components of body j
       end subroutine kick_getacch_int_one_pl
 
-      module pure subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, ax, ay, az)
+      pure module subroutine kick_getacch_int_one_tp(rji2, xr, yr, zr, Gmpl, ax, ay, az)
          !$omp declare simd(kick_getacch_int_one_tp)
          implicit none
          real(DP), intent(in)  :: rji2         !! Square of distance between the test particle and massive body
@@ -1006,14 +1006,14 @@ module subroutine orbel_el2xv_vec(self, cb)
          class(swiftest_cb),           intent(inout) :: cb   !! Swiftest central body object
       end subroutine orbel_el2xv_vec
 
-      module pure subroutine orbel_scget(angle, sx, cx)
+      pure module subroutine orbel_scget(angle, sx, cx)
          !$omp declare simd(orbel_scget)
          implicit none
          real(DP), intent(in)  :: angle
          real(DP), intent(out) :: sx, cx
       end subroutine orbel_scget
 
-      module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q)
+      pure module subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q)
          !$omp declare simd(orbel_xv2aeq)
          implicit none
          real(DP), intent(in)  :: mu !! Gravitational constant
@@ -1024,7 +1024,7 @@ module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q)
          real(DP), intent(out) :: q  !! periapsis
       end subroutine orbel_xv2aeq
 
-      module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi)
+      pure module subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi)
          !$omp declare simd(orbel_xv2aqt)
          implicit none
          real(DP), intent(in)  :: mu    !! Gravitational constant
@@ -1036,7 +1036,7 @@ module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tper
          real(DP), intent(out) :: tperi !! time of pericenter passage
       end subroutine orbel_xv2aqt
 
-      module pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
+      pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
          implicit none
          real(DP), intent(in)  :: mu    !! Gravitational constant
          real(DP), intent(in)  :: px,py,pz    !! Position vector
@@ -1361,7 +1361,7 @@ end subroutine util_fill_arr_logical
    end interface
 
    interface
-      module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
+      pure module subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
          !$omp declare simd(util_flatten_eucl_ij_to_k)
          implicit none
          integer(I4B), intent(in)  :: n !! Number of bodies
@@ -1370,7 +1370,7 @@ module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
          integer(I8B), intent(out) :: k !! Index of the flattened matrix
       end subroutine util_flatten_eucl_ij_to_k
 
-      module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j)
+      pure module subroutine util_flatten_eucl_k_to_ij(n, k, i, j)
          implicit none
          integer(I4B), intent(in)  :: n !! Number of bodies
          integer(I8B), intent(in)  :: k !! Index of the flattened matrix
@@ -1596,46 +1596,46 @@ end function util_solve_rkf45
    end interface
 
    interface util_sort
-      module pure subroutine util_sort_i4b(arr)
+      pure module subroutine util_sort_i4b(arr)
          implicit none
          integer(I4B), dimension(:), intent(inout) :: arr
       end subroutine util_sort_i4b
 
-      module pure subroutine util_sort_index_i4b(arr,ind)
+      pure module subroutine util_sort_index_i4b(arr,ind)
          implicit none
          integer(I4B), dimension(:), intent(in)  :: arr
          integer(I4B), dimension(:), allocatable, intent(inout) :: ind
       end subroutine util_sort_index_i4b
 
-      module pure subroutine util_sort_index_I4B_I8Bind(arr,ind)
+      pure module subroutine util_sort_index_I4B_I8Bind(arr,ind)
          implicit none
          integer(I4B), dimension(:), intent(in)  :: arr
          integer(I8B), dimension(:), allocatable, intent(inout) :: ind
       end subroutine util_sort_index_I4b_I8Bind
 
-      module pure subroutine util_sort_index_I8B_I8Bind(arr,ind)
+      pure module subroutine util_sort_index_I8B_I8Bind(arr,ind)
          implicit none
          integer(I8B), dimension(:), intent(in)  :: arr
          integer(I8B), dimension(:), allocatable, intent(inout) :: ind
       end subroutine util_sort_index_I8B_I8Bind
 
-      module pure subroutine util_sort_sp(arr)
+      pure module subroutine util_sort_sp(arr)
          implicit none
          real(SP), dimension(:), intent(inout) :: arr
       end subroutine util_sort_sp
 
-      module pure subroutine util_sort_index_sp(arr,ind)
+      pure module subroutine util_sort_index_sp(arr,ind)
          implicit none
          real(SP), dimension(:), intent(in)  :: arr
          integer(I4B), dimension(:), allocatable, intent(inout) :: ind
       end subroutine util_sort_index_sp
 
-      module pure subroutine util_sort_dp(arr)
+      pure module subroutine util_sort_dp(arr)
          implicit none
          real(DP), dimension(:), intent(inout) :: arr
       end subroutine util_sort_dp
 
-      module pure subroutine util_sort_index_dp(arr,ind)
+      pure module subroutine util_sort_index_dp(arr,ind)
          implicit none
          real(DP), dimension(:), intent(in)  :: arr
          integer(I4B), dimension(:), allocatable, intent(inout) :: ind
@@ -1643,35 +1643,35 @@ end subroutine util_sort_index_dp
    end interface util_sort
 
    interface util_sort_rearrange
-      module pure subroutine util_sort_rearrange_arr_char_string(arr, ind, n)
+      pure module subroutine util_sort_rearrange_arr_char_string(arr, ind, n)
          implicit none
          character(len=STRMAX), dimension(:), allocatable, intent(inout) :: arr !! Destination array 
          integer(I4B),          dimension(:),              intent(in)    :: ind !! Index to rearrange against
          integer(I4B),                                     intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       end subroutine util_sort_rearrange_arr_char_string
 
-      module pure subroutine util_sort_rearrange_arr_DP(arr, ind, n)
+      pure module subroutine util_sort_rearrange_arr_DP(arr, ind, n)
          implicit none
          real(DP),     dimension(:), allocatable, intent(inout) :: arr !! Destination array 
          integer(I4B), dimension(:),              intent(in)  :: ind !! Index to rearrange against
          integer(I4B),                            intent(in)  :: n   !! Number of elements in arr and ind to rearrange
       end subroutine util_sort_rearrange_arr_DP
 
-      module pure subroutine util_sort_rearrange_arr_DPvec(arr, ind, n)
+      pure module subroutine util_sort_rearrange_arr_DPvec(arr, ind, n)
          implicit none
          real(DP),     dimension(:,:), allocatable, intent(inout) :: arr !! Destination array 
          integer(I4B), dimension(:),                intent(in)    :: ind !! Index to rearrange against
          integer(I4B),                              intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       end subroutine util_sort_rearrange_arr_DPvec
 
-      module pure subroutine util_sort_rearrange_arr_I4B(arr, ind, n)
+      pure module subroutine util_sort_rearrange_arr_I4B(arr, ind, n)
          implicit none
          integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array 
          integer(I4B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
          integer(I4B),                             intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       end subroutine util_sort_rearrange_arr_I4B
 
-      module pure subroutine util_sort_rearrange_arr_I4B_I8Bind(arr, ind, n)
+      pure module subroutine util_sort_rearrange_arr_I4B_I8Bind(arr, ind, n)
          implicit none
          integer(I4B), dimension(:), allocatable, intent(inout) :: arr !! Destination array 
          integer(I8B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
@@ -1685,14 +1685,14 @@ module subroutine util_sort_rearrange_arr_info(arr, ind, n)
          integer(I4B),                                   intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       end subroutine util_sort_rearrange_arr_info
 
-      module pure subroutine util_sort_rearrange_arr_logical(arr, ind, n)
+      pure module subroutine util_sort_rearrange_arr_logical(arr, ind, n)
          implicit none
          logical,      dimension(:), allocatable, intent(inout) :: arr !! Destination array 
          integer(I4B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
          integer(I4B),                            intent(in)    :: n   !! Number of elements in arr and ind to rearrange
       end subroutine util_sort_rearrange_arr_logical
 
-      module pure subroutine util_sort_rearrange_arr_logical_I8Bind(arr, ind, n)
+      pure module subroutine util_sort_rearrange_arr_logical_I8Bind(arr, ind, n)
          implicit none
          logical,      dimension(:), allocatable, intent(inout) :: arr !! Destination array 
          integer(I8B), dimension(:),              intent(in)    :: ind !! Index to rearrange against
diff --git a/src/modules/swiftest_operators.f90 b/src/modules/swiftest_operators.f90
index 13cb57839..770bc0047 100644
--- a/src/modules/swiftest_operators.f90
+++ b/src/modules/swiftest_operators.f90
@@ -14,92 +14,92 @@ module swiftest_operators
    !********************************************************************************************************************************
 
    interface operator(.cross.)
-      module pure function operator_cross_sp(A, B) result(C)
+      pure module function operator_cross_sp(A, B) result(C)
          !$omp declare simd(operator_cross_sp)
          implicit none
          real(SP), dimension(:), intent(in) :: A, B
          real(SP), dimension(3) :: C
       end function operator_cross_sp
 
-      module pure function operator_cross_dp(A, B) result(C)
+      pure module function operator_cross_dp(A, B) result(C)
          !$omp declare simd(operator_cross_dp)
          implicit none
          real(DP), dimension(:), intent(in) :: A, B
          real(DP), dimension(3) :: C
       end function operator_cross_dp
 
-      module pure function operator_cross_qp(A, B) result(C)
+      pure module function operator_cross_qp(A, B) result(C)
          !$omp declare simd(operator_cross_qp)
          implicit none
          real(QP), dimension(:), intent(in) :: A, B
          real(QP), dimension(3) :: C
       end function operator_cross_qp
 
-      module pure function operator_cross_i1b(A, B) result(C)
+      pure module function operator_cross_i1b(A, B) result(C)
          !$omp declare simd(operator_cross_i1b)
          implicit none
          integer(I1B), dimension(:), intent(in) :: A, B
          integer(I1B), dimension(3) :: C
       end function operator_cross_i1b
 
-      module pure function operator_cross_i2b(A, B) result(C)
+      pure module function operator_cross_i2b(A, B) result(C)
          !$omp declare simd(operator_cross_i2b)
          implicit none
          integer(I2B), dimension(:), intent(in) :: A, B
          integer(I2B), dimension(3) :: C
       end function operator_cross_i2b
 
-      module pure function operator_cross_i4b(A, B) result(C)
+      pure module function operator_cross_i4b(A, B) result(C)
          !$omp declare simd(operator_cross_i4b)
          implicit none
          integer(I4B), dimension(:), intent(in) :: A, B
          integer(I4B), dimension(3) :: C
       end function operator_cross_i4b
 
-      module pure function operator_cross_i8b(A, B) result(C)
+      pure module function operator_cross_i8b(A, B) result(C)
          !$omp declare simd(operator_cross_i8b)
          implicit none
          integer(I8B), dimension(:), intent(in) :: A, B
          integer(I8B), dimension(3) :: C
       end function operator_cross_i8b
 
-      module pure function operator_cross_el_sp(A, B) result(C)
+      pure module function operator_cross_el_sp(A, B) result(C)
          implicit none
          real(SP), dimension(:,:), intent(in) :: A, B
          real(SP), dimension(:,:), allocatable :: C
       end function operator_cross_el_sp
 
-      module pure function operator_cross_el_dp(A, B) result(C)
+      pure module function operator_cross_el_dp(A, B) result(C)
          implicit none
          real(DP), dimension(:,:), intent(in) :: A, B
          real(DP), dimension(:,:), allocatable :: C
       end function operator_cross_el_dp
 
-      module pure function operator_cross_el_qp(A, B) result(C)
+      pure module function operator_cross_el_qp(A, B) result(C)
          implicit none
          real(QP), dimension(:,:), intent(in) :: A, B
          real(QP), dimension(:,:), allocatable :: C
       end function operator_cross_el_qp
 
-      module pure function operator_cross_el_i1b(A, B) result(C)
+      pure module function operator_cross_el_i1b(A, B) result(C)
          implicit none
          integer(I1B), dimension(:,:), intent(in) :: A, B
          integer(I1B), dimension(:,:), allocatable :: C
       end function operator_cross_el_i1b
 
-      module pure function operator_cross_el_i2b(A, B) result(C)
+      pure module function operator_cross_el_i2b(A, B) result(C)
          implicit none
          integer(I2B), dimension(:,:), intent(in) :: A, B
          integer(I2B), dimension(:,:), allocatable :: C
       end function operator_cross_el_i2b
 
-      module pure function operator_cross_el_i4b(A, B) result(C)
+      pure module function operator_cross_el_i4b(A, B) result(C)
          implicit none
          integer(I4B), dimension(:,:), intent(in) :: A, B
          integer(I4B), dimension(:,:), allocatable :: C
       end function operator_cross_el_i4b
 
-      module pure function operator_cross_el_i8b(A, B) result(C)
+      pure module function operator_cross_el_i8b(A, B) result(C)
          implicit none
          integer(I8B), dimension(:,:), intent(in) :: A, B
          integer(I8B), dimension(:,:), allocatable :: C
@@ -111,40 +111,40 @@ end function operator_cross_el_i8b
    !********************************************************************************************************************************
 
    interface operator(.mag.)
-      module pure function operator_mag_sp(A) result(B)
+      pure module function operator_mag_sp(A) result(B)
          !$omp declare simd(operator_mag_sp)
          implicit none
          real(SP), dimension(:), intent(in) :: A
          real(SP)                           :: B
       end function operator_mag_sp
 
-      module pure function operator_mag_dp(A) result(B)
+      pure module function operator_mag_dp(A) result(B)
          !$omp declare simd(operator_mag_dp)
          implicit none
          real(DP), dimension(:), intent(in) :: A
          real(DP)                           :: B
       end function operator_mag_dp
 
-      module pure function operator_mag_qp(A) result(B)
+      pure module function operator_mag_qp(A) result(B)
          !$omp declare simd(operator_mag_qp)
          implicit none
          real(QP), dimension(:), intent(in) :: A
          real(QP)                           :: B
       end function operator_mag_qp
 
-      module pure function operator_mag_el_sp(A) result(B)
+      pure module function operator_mag_el_sp(A) result(B)
          implicit none
          real(SP), dimension(:,:), intent(in) :: A
          real(SP), dimension(:), allocatable  :: B
       end function operator_mag_el_sp
 
-      module pure function operator_mag_el_dp(A) result(B)
+      pure module function operator_mag_el_dp(A) result(B)
          implicit none
          real(DP), dimension(:,:), intent(in) :: A
          real(DP), dimension(:), allocatable  :: B
       end function operator_mag_el_dp
 
-      module pure function operator_mag_el_qp(A) result(B)
+      pure module function operator_mag_el_qp(A) result(B)
          implicit none
          real(QP), dimension(:,:), intent(in) :: A
          real(QP), dimension(:), allocatable  :: B
diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90
index a01d60e6a..f501c1f4f 100644
--- a/src/modules/symba_classes.f90
+++ b/src/modules/symba_classes.f90
@@ -306,7 +306,7 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l
          logical                                   :: lany_encounter !! Returns true if there is at least one close encounter      
       end function symba_encounter_check_tp
 
-      module pure subroutine symba_gr_p4_pl(self, system, param, dt)
+      pure module subroutine symba_gr_p4_pl(self, system, param, dt)
          use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
          implicit none
          class(symba_pl),              intent(inout) :: self   !! SyMBA massive body object
@@ -315,7 +315,7 @@ module pure subroutine symba_gr_p4_pl(self, system, param, dt)
          real(DP),                     intent(in)    :: dt     !! Step size
       end subroutine symba_gr_p4_pl
    
-      module pure subroutine symba_gr_p4_tp(self, system, param, dt)
+      pure module subroutine symba_gr_p4_tp(self, system, param, dt)
          use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system
          implicit none
          class(symba_tp),              intent(inout) :: self   !! SyMBA test particle object
diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90
index a7cd2f49c..2069564e9 100644
--- a/src/modules/whm_classes.f90
+++ b/src/modules/whm_classes.f90
@@ -156,21 +156,21 @@ module subroutine whm_kick_vh_tp(self, system, param, t, dt, lbeg)
          logical,                      intent(in)    :: lbeg   !! Logical flag indicating whether this is the beginning of the half step or not. 
       end subroutine whm_kick_vh_tp
 
-      module pure subroutine whm_gr_kick_getacch_pl(self, param)
+      pure module subroutine whm_gr_kick_getacch_pl(self, param)
          use swiftest_classes, only : swiftest_cb, swiftest_parameters
          implicit none
          class(whm_pl),              intent(inout) :: self  !! WHM massive body particle data structure
          class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters 
       end subroutine whm_gr_kick_getacch_pl
 
-      module pure subroutine whm_gr_kick_getacch_tp(self, param)
+      pure module subroutine whm_gr_kick_getacch_tp(self, param)
          use swiftest_classes, only : swiftest_cb, swiftest_parameters
          implicit none
          class(whm_tp),              intent(inout) :: self  !! WHM test particle data structure
          class(swiftest_parameters), intent(in)    :: param !! Current run configuration parameters
       end subroutine whm_gr_kick_getacch_tp
 
-      module pure subroutine whm_gr_p4_pl(self, system, param, dt)
+      pure module subroutine whm_gr_p4_pl(self, system, param, dt)
          use swiftest_classes, only : swiftest_parameters
          implicit none
          class(whm_pl),                intent(inout) :: self  !! WHM massive body object
@@ -179,7 +179,7 @@ module pure subroutine whm_gr_p4_pl(self, system, param, dt)
          real(DP),                     intent(in)    :: dt    !! Step size
       end subroutine whm_gr_p4_pl
 
-      module pure subroutine whm_gr_p4_tp(self, system, param, dt)
+      pure module subroutine whm_gr_p4_tp(self, system, param, dt)
          use swiftest_classes, only : swiftest_parameters
          implicit none
          class(whm_tp),                intent(inout) :: self   !! WHM test particle object
diff --git a/src/operators/operator_cross.f90 b/src/operators/operator_cross.f90
index 5a9f8fd28..edb802535 100644
--- a/src/operators/operator_cross.f90
+++ b/src/operators/operator_cross.f90
@@ -7,7 +7,7 @@
    !! Vector list implementations:   C(1:3, :) = A(1:3, :) .cross. B(1:3, :)
 contains
 
-   module pure function operator_cross_sp(A, B) result(C)
+   pure module function operator_cross_sp(A, B) result(C)
       implicit none
       real(SP), dimension(:), intent(in) :: A, B
       real(SP), dimension(3) :: C
@@ -17,7 +17,7 @@ module pure function operator_cross_sp(A, B) result(C)
       return
    end function operator_cross_sp
 
-   module pure function operator_cross_dp(A, B) result(C)
+   pure module function operator_cross_dp(A, B) result(C)
       implicit none
       real(DP), dimension(:), intent(in) :: A, B
       real(DP), dimension(3) :: C
@@ -27,7 +27,7 @@ module pure function operator_cross_dp(A, B) result(C)
       return
    end function operator_cross_dp
 
-   module pure function operator_cross_qp(A, B) result(C)
+   pure module function operator_cross_qp(A, B) result(C)
       implicit none
       real(QP), dimension(:), intent(in) :: A, B
       real(QP), dimension(3) :: C
@@ -37,7 +37,7 @@ module pure function operator_cross_qp(A, B) result(C)
       return
    end function operator_cross_qp
 
-   module pure function operator_cross_i1b(A, B) result(C)
+   pure module function operator_cross_i1b(A, B) result(C)
       implicit none
       integer(I1B), dimension(:), intent(in) :: A, B
       integer(I1B), dimension(3) :: C
@@ -47,7 +47,7 @@ module pure function operator_cross_i1b(A, B) result(C)
       return
    end function operator_cross_i1b
 
-   module pure function operator_cross_i2b(A, B) result(C)
+   pure module function operator_cross_i2b(A, B) result(C)
       implicit none
       integer(I2B), dimension(:), intent(in) :: A, B
       integer(I2B), dimension(3) :: C
@@ -57,7 +57,7 @@ module pure function operator_cross_i2b(A, B) result(C)
       return
    end function operator_cross_i2b
 
-   module pure function operator_cross_i4b(A, B) result(C)
+   pure module function operator_cross_i4b(A, B) result(C)
       implicit none
       integer(I4B), dimension(:), intent(in) :: A, B
       integer(I4B), dimension(3) :: C
@@ -67,7 +67,7 @@ module pure function operator_cross_i4b(A, B) result(C)
       return
    end function operator_cross_i4b
 
-   module pure function operator_cross_i8b(A, B) result(C)     
+   pure module function operator_cross_i8b(A, B) result(C)     
       implicit none
       integer(I8B), dimension(:), intent(in) :: A, B
       integer(I8B), dimension(3) :: C
@@ -77,7 +77,7 @@ module pure function operator_cross_i8b(A, B) result(C)
       return
    end function operator_cross_i8b
 
-   module pure function operator_cross_el_sp(A, B) result(C)
+   pure module function operator_cross_el_sp(A, B) result(C)
       implicit none
       real(SP), dimension(:,:), intent(in)  :: A, B
       real(SP), dimension(:,:), allocatable :: C
@@ -91,7 +91,7 @@ module pure function operator_cross_el_sp(A, B) result(C)
       return
    end function operator_cross_el_sp
 
-   module pure function operator_cross_el_dp(A, B) result(C)
+   pure module function operator_cross_el_dp(A, B) result(C)
       implicit none
       real(DP), dimension(:,:), intent(in)  :: A, B
       real(DP), dimension(:,:), allocatable :: C
@@ -105,7 +105,7 @@ module pure function operator_cross_el_dp(A, B) result(C)
       return
    end function operator_cross_el_dp
 
-   module pure function operator_cross_el_qp(A, B) result(C)
+   pure module function operator_cross_el_qp(A, B) result(C)
       implicit none
       real(QP), dimension(:,:), intent(in)  :: A, B
       real(QP), dimension(:,:), allocatable :: C
@@ -119,7 +119,7 @@ module pure function operator_cross_el_qp(A, B) result(C)
       return
    end function operator_cross_el_qp
 
-   module pure function operator_cross_el_i1b(A, B) result(C)
+   pure module function operator_cross_el_i1b(A, B) result(C)
       implicit none
       integer(I1B), dimension(:,:), intent(in)  :: A, B
       integer(I1B), dimension(:,:), allocatable :: C
@@ -133,7 +133,7 @@ module pure function operator_cross_el_i1b(A, B) result(C)
       return
    end function operator_cross_el_i1b
 
-   module pure function operator_cross_el_i2b(A, B) result(C)
+   pure module function operator_cross_el_i2b(A, B) result(C)
       implicit none
       integer(I2B), dimension(:,:), intent(in)  :: A, B
       integer(I2B), dimension(:,:), allocatable :: C
@@ -147,7 +147,7 @@ module pure function operator_cross_el_i2b(A, B) result(C)
       return
    end function operator_cross_el_i2b
 
-   module pure function operator_cross_el_i4b(A, B) result(C)
+   pure module function operator_cross_el_i4b(A, B) result(C)
       implicit none
       integer(I4B), dimension(:,:), intent(in)  :: A, B
       integer(I4B), dimension(:,:), allocatable :: C
@@ -161,7 +161,7 @@ module pure function operator_cross_el_i4b(A, B) result(C)
       return
    end function operator_cross_el_i4b
 
-   module pure function operator_cross_el_i8b(A, B) result(C)
+   pure module function operator_cross_el_i8b(A, B) result(C)
       implicit none
       integer(I8B), dimension(:,:), intent(in)  :: A, B
       integer(I8B), dimension(:,:), allocatable :: C
diff --git a/src/operators/operator_mag.f90 b/src/operators/operator_mag.f90
index 5a054d5ce..92c19cb4b 100644
--- a/src/operators/operator_mag.f90
+++ b/src/operators/operator_mag.f90
@@ -6,7 +6,7 @@
    !! Vector list implementations:   B(:) = .mag. A(1:3, :)
    contains
 
-   module pure function operator_mag_sp(A) result(B)
+   pure module function operator_mag_sp(A) result(B)
       implicit none
       real(SP), dimension(:), intent(in) :: A
       real(SP)                           :: B
@@ -14,7 +14,7 @@ module pure function operator_mag_sp(A) result(B)
       return
    end function operator_mag_sp
 
-   module pure function operator_mag_dp(A) result(B)
+   pure module function operator_mag_dp(A) result(B)
       implicit none
       real(DP), dimension(:), intent(in) :: A
       real(DP)                           :: B
@@ -22,7 +22,7 @@ module pure function operator_mag_dp(A) result(B)
       return
    end function operator_mag_dp
 
-   module pure function operator_mag_el_sp(A) result(B)
+   pure module function operator_mag_el_sp(A) result(B)
       implicit none
       real(SP), dimension(:,:), intent(in) :: A
       real(SP), dimension(:), allocatable  :: B
@@ -36,7 +36,7 @@ module pure function operator_mag_el_sp(A) result(B)
       return
    end function operator_mag_el_sp
 
-   module pure function operator_mag_el_dp(A) result(B)
+   pure module function operator_mag_el_dp(A) result(B)
       implicit none
       real(DP), dimension(:,:), intent(in) :: A
       real(DP), dimension(:), allocatable  :: B
@@ -50,7 +50,7 @@ module pure function operator_mag_el_dp(A) result(B)
       return 
    end function operator_mag_el_dp
 
-   module pure function operator_mag_el_qp(A) result(B)
+   pure module function operator_mag_el_qp(A) result(B)
       implicit none
       real(QP), dimension(:,:), intent(in) :: A
       real(QP), dimension(:), allocatable  :: B
diff --git a/src/orbel/orbel.f90 b/src/orbel/orbel.f90
index e35d1e20f..31f7f23c3 100644
--- a/src/orbel/orbel.f90
+++ b/src/orbel/orbel.f90
@@ -128,7 +128,7 @@ pure subroutine orbel_el2xv(mu, a, ie, inc, capom, omega, capm, x, v)
    end subroutine orbel_el2xv
 
 
-   module pure subroutine orbel_scget(angle, sx, cx)
+   pure module subroutine orbel_scget(angle, sx, cx)
       !! author: David A. Minton
       !!
       !! Efficiently compute the sine and cosine of an input angle
@@ -683,7 +683,7 @@ real(DP) pure function orbel_fhybrid(e,n)
    end function orbel_fhybrid
    
 
-   module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q)
+   pure module subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q)
       !! author: David A. Minton
       !!
       !! Compute semimajor axis, eccentricity, and pericentric distance from relative Cartesian position and velocity
@@ -748,7 +748,7 @@ module pure subroutine orbel_xv2aeq(mu, px, py, pz, vx, vy, vz, a, e, q)
    end subroutine orbel_xv2aeq
 
 
-   module pure subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi)
+   pure module subroutine orbel_xv2aqt(mu, px, py, pz, vx, vy, vz, a, q, capm, tperi)
       !! author: David A. Minton
       !!
       !! Compute semimajor axis, pericentric distance, mean anomaly, and time to nearest pericenter passage from
@@ -888,7 +888,7 @@ module subroutine orbel_xv2el_vec(self, cb)
    end subroutine orbel_xv2el_vec 
 
 
-   module pure subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
+   pure module subroutine orbel_xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm)
       !! author: David A. Minton
       !!
       !! Compute osculating orbital elements from relative Cartesian position and velocity
diff --git a/src/symba/symba_gr.f90 b/src/symba/symba_gr.f90
index f743b124a..6340bd429 100644
--- a/src/symba/symba_gr.f90
+++ b/src/symba/symba_gr.f90
@@ -2,7 +2,7 @@
    use swiftest
 contains
 
-   module pure subroutine symba_gr_p4_pl(self, system, param, dt)
+   pure module subroutine symba_gr_p4_pl(self, system, param, dt)
       !! author: David A. Minton
       !!
       !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction
@@ -31,7 +31,7 @@ module pure subroutine symba_gr_p4_pl(self, system, param, dt)
    end subroutine symba_gr_p4_pl
 
 
-   module pure subroutine symba_gr_p4_tp(self, system, param, dt)
+   pure module subroutine symba_gr_p4_tp(self, system, param, dt)
       !! author: David A. Minton
       !!
       !! Position kick to test particles due to p**4 term in the post-Newtonian correction
diff --git a/src/util/util_flatten.f90 b/src/util/util_flatten.f90
index 1f6f60fe6..5305b2fe7 100644
--- a/src/util/util_flatten.f90
+++ b/src/util/util_flatten.f90
@@ -2,7 +2,7 @@
    use swiftest
 contains
 
-   module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
+   pure module subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
       !! author: Jacob R. Elliott and David A. Minton
       !!
       !! Turns i,j indices into k index for use in the Euclidean distance matrix for pl-pl interactions.
@@ -29,7 +29,7 @@ module pure subroutine util_flatten_eucl_ij_to_k(n, i, j, k)
    end subroutine util_flatten_eucl_ij_to_k
 
 
-   module pure subroutine util_flatten_eucl_k_to_ij(n, k, i, j)
+   pure module subroutine util_flatten_eucl_k_to_ij(n, k, i, j)
       !! author: Jacob R. Elliott and David A. Minton
       !!
       !! Turns k index into i,j indices for use in the Euclidean distance matrix for pl-pl interactions.
diff --git a/src/util/util_sort.f90 b/src/util/util_sort.f90
index 453c3a2d3..ebc6223c4 100644
--- a/src/util/util_sort.f90
+++ b/src/util/util_sort.f90
@@ -57,7 +57,7 @@ module subroutine util_sort_body(self, sortby, ascending)
    end subroutine util_sort_body
 
 
-   module pure subroutine util_sort_dp(arr)
+   pure module subroutine util_sort_dp(arr)
       !! author: David A. Minton
       !!
       !! Sort input DP precision array in place into ascending numerical order using quicksort.
@@ -72,7 +72,7 @@ module pure subroutine util_sort_dp(arr)
    end subroutine util_sort_dp
 
 
-   module pure subroutine util_sort_index_dp(arr, ind)
+   pure module subroutine util_sort_index_dp(arr, ind)
       !! author: David A. Minton
       !!
       !! Sort input DP precision array by index in ascending numerical order using quick sort.
@@ -186,7 +186,7 @@ pure subroutine partition_DP(arr, marker, ind)
    end subroutine partition_DP
  
 
-   module pure subroutine util_sort_i4b(arr)
+   pure module subroutine util_sort_i4b(arr)
       !! author: David A. Minton
       !!
       !! Sort input integer array in place into ascending numerical order using quick sort.
@@ -202,7 +202,7 @@ module pure subroutine util_sort_i4b(arr)
    end subroutine util_sort_i4b
 
 
-   module pure subroutine util_sort_index_I4B(arr, ind)
+   pure module subroutine util_sort_index_I4B(arr, ind)
       !! author: David A. Minton
       !!
       !! Sort input integer array by index in ascending numerical order using quicksort.
@@ -230,7 +230,7 @@ module pure subroutine util_sort_index_I4B(arr, ind)
    end subroutine util_sort_index_I4B
 
 
-   module pure subroutine util_sort_index_I4B_I8Bind(arr, ind)
+   pure module subroutine util_sort_index_I4B_I8Bind(arr, ind)
       !! author: David A. Minton
       !!
       !! Sort input integer array by index in ascending numerical order using quicksort.
@@ -510,7 +510,7 @@ pure subroutine partition_I8B_I8Bind(arr, marker, ind)
    end subroutine partition_I8B_I8Bind
 
 
-   module pure subroutine util_sort_sp(arr)
+   pure module subroutine util_sort_sp(arr)
       !! author: David A. Minton
       !!
       !! Sort input DP precision array in place into ascending numerical order using quicksort.
@@ -525,7 +525,7 @@ module pure subroutine util_sort_sp(arr)
    end subroutine util_sort_sp
 
 
-   module pure subroutine util_sort_index_sp(arr, ind)
+   pure module subroutine util_sort_index_sp(arr, ind)
       !! author: David A. Minton
       !!
       !! Sort input DP precision array by index in ascending numerical order using quicksort.
@@ -774,7 +774,7 @@ module subroutine util_sort_rearrange_body(self, ind)
    end subroutine util_sort_rearrange_body
 
 
-   module pure subroutine util_sort_rearrange_arr_char_string(arr, ind, n)
+   pure module subroutine util_sort_rearrange_arr_char_string(arr, ind, n)
       !! author: David A. Minton
       !!
       !! Rearrange a single array of character string in-place from an index list.
@@ -795,7 +795,7 @@ module pure subroutine util_sort_rearrange_arr_char_string(arr, ind, n)
    end subroutine util_sort_rearrange_arr_char_string
 
 
-   module pure subroutine util_sort_rearrange_arr_DP(arr, ind, n)
+   pure module subroutine util_sort_rearrange_arr_DP(arr, ind, n)
       !! author: David A. Minton
       !!
       !! Rearrange a single array of DP type in-place from an index list.
@@ -816,7 +816,7 @@ module pure subroutine util_sort_rearrange_arr_DP(arr, ind, n)
    end subroutine util_sort_rearrange_arr_DP
 
 
-   module pure subroutine util_sort_rearrange_arr_DPvec(arr, ind, n)
+   pure module subroutine util_sort_rearrange_arr_DPvec(arr, ind, n)
       !! author: David A. Minton
       !!
       !! Rearrange a single array of (NDIM,n) DP-type vectors in-place from an index list.
@@ -837,7 +837,7 @@ module pure subroutine util_sort_rearrange_arr_DPvec(arr, ind, n)
    end subroutine util_sort_rearrange_arr_DPvec
 
 
-   module pure subroutine util_sort_rearrange_arr_I4B(arr, ind, n)
+   pure module subroutine util_sort_rearrange_arr_I4B(arr, ind, n)
       !! author: David A. Minton
       !!
       !! Rearrange a single array of integers in-place from an index list.
@@ -857,7 +857,7 @@ module pure subroutine util_sort_rearrange_arr_I4B(arr, ind, n)
       return
    end subroutine util_sort_rearrange_arr_I4B
 
-   module pure subroutine util_sort_rearrange_arr_I4B_I8Bind(arr, ind, n)
+   pure module subroutine util_sort_rearrange_arr_I4B_I8Bind(arr, ind, n)
       !! author: David A. Minton
       !!
       !! Rearrange a single array of integers in-place from an index list.
@@ -878,7 +878,7 @@ module pure subroutine util_sort_rearrange_arr_I4B_I8Bind(arr, ind, n)
    end subroutine util_sort_rearrange_arr_I4B_I8Bind
 
 
-   module pure subroutine util_sort_rearrange_arr_logical(arr, ind, n)
+   pure module subroutine util_sort_rearrange_arr_logical(arr, ind, n)
       !! author: David A. Minton
       !!
       !! Rearrange a single array of logicals in-place from an index list.
@@ -899,7 +899,7 @@ module pure subroutine util_sort_rearrange_arr_logical(arr, ind, n)
    end subroutine util_sort_rearrange_arr_logical
 
 
-   module pure subroutine util_sort_rearrange_arr_logical_I8Bind(arr, ind, n)
+   pure module subroutine util_sort_rearrange_arr_logical_I8Bind(arr, ind, n)
       !! author: David A. Minton
       !!
       !! Rearrange a single array of logicals in-place from an index list.
diff --git a/src/whm/whm_gr.f90 b/src/whm/whm_gr.f90
index 12ae82a35..36ca5629e 100644
--- a/src/whm/whm_gr.f90
+++ b/src/whm/whm_gr.f90
@@ -2,7 +2,7 @@
    use swiftest
 contains
 
-   module pure subroutine whm_gr_kick_getacch_pl(self, param) 
+   pure module subroutine whm_gr_kick_getacch_pl(self, param) 
       !! author: David A. Minton
       !!
       !! Compute relativisitic accelerations of massive bodies
@@ -35,7 +35,7 @@ module pure subroutine whm_gr_kick_getacch_pl(self, param)
    end subroutine whm_gr_kick_getacch_pl
 
 
-   module pure subroutine whm_gr_kick_getacch_tp(self, param)
+   pure module subroutine whm_gr_kick_getacch_tp(self, param)
       !! author: David A. Minton
       !!
       !! Compute relativisitic accelerations of test particles
@@ -61,7 +61,7 @@ module pure subroutine whm_gr_kick_getacch_tp(self, param)
    end subroutine whm_gr_kick_getacch_tp
    
 
-   module pure subroutine whm_gr_p4_pl(self, system, param, dt)
+   pure module subroutine whm_gr_p4_pl(self, system, param, dt)
       !! author: David A. Minton
       !!
       !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction
@@ -88,7 +88,7 @@ module pure subroutine whm_gr_p4_pl(self, system, param, dt)
    end subroutine whm_gr_p4_pl
 
 
-   module pure subroutine whm_gr_p4_tp(self, system, param, dt)
+   pure module subroutine whm_gr_p4_tp(self, system, param, dt)
       !! author: David A. Minton
       !!
       !! Position kick to test particles due to p**4 term in the post-Newtonian correction