diff --git a/src/linalg_inverse.f90 b/src/linalg_inverse.f90
index 954cfe7c..d5b0290d 100644
--- a/src/linalg_inverse.f90
+++ b/src/linalg_inverse.f90
@@ -362,7 +362,7 @@ subroutine mtx_pinverse_dbl(a, ainv, tol, work, olwork, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine mtx_pinverse_cmplx(a, ainv, tol, work, olwork, rwork, err)
+subroutine mtx_pinverse_cmplx(a, ainv, tol, work, olwork, rwork, err)
     !! Computes the Moore-Penrose pseudo-inverse of a M-by-N matrix using the
     !! singular value decomposition of the matrix.
     complex(real64), intent(inout), dimension(:,:) :: a
diff --git a/src/linalg_least_squares.f90 b/src/linalg_least_squares.f90
index 5088d547..00eb9a9c 100644
--- a/src/linalg_least_squares.f90
+++ b/src/linalg_least_squares.f90
@@ -1157,7 +1157,7 @@ subroutine solve_least_squares_mtx_svd_cmplx(a, b, s, arnk, work, &
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine solve_least_squares_vec_svd(a, b, s, arnk, work, olwork, err)
+subroutine solve_least_squares_vec_svd(a, b, s, arnk, work, olwork, err)
     !! Solves the system of equations \(A \vec{x} = \vec{b}\) using a singular 
     !! value decomposition of \(A\).
     real(real64), intent(inout), dimension(:,:) :: a
@@ -1271,7 +1271,7 @@ module subroutine solve_least_squares_vec_svd(a, b, s, arnk, work, olwork, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine solve_least_squares_vec_svd_cmplx(a, b, s, arnk, work, &
+subroutine solve_least_squares_vec_svd_cmplx(a, b, s, arnk, work, &
     olwork, rwork, err)
     !! Solves the system of equations \(A \vec{x} = \vec{b}\) using a singular 
     !! value decomposition of \(A\).
diff --git a/src/linalg_lq.f90 b/src/linalg_lq.f90
index 1d3d6590..4a34290a 100644
--- a/src/linalg_lq.f90
+++ b/src/linalg_lq.f90
@@ -878,7 +878,7 @@ subroutine solve_lq_mtx(a, tau, b, work, olwork, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine solve_lq_mtx_cmplx(a, tau, b, work, olwork, err)
+subroutine solve_lq_mtx_cmplx(a, tau, b, work, olwork, err)
     !! Solves a system of LQ factored equations of the form \(A X = L Q X = B\).
     complex(real64), intent(in), dimension(:,:) :: a
         !! On input, the M-by-N LQ factored matrix as returned by lq_factor.  
diff --git a/src/linalg_qr.f90 b/src/linalg_qr.f90
index 008d55ea..4157375f 100644
--- a/src/linalg_qr.f90
+++ b/src/linalg_qr.f90
@@ -692,7 +692,7 @@ subroutine form_qr_pivot(r, tau, pvt, q, p, work, olwork, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine form_qr_pivot_cmplx(r, tau, pvt, q, p, work, olwork, err)
+subroutine form_qr_pivot_cmplx(r, tau, pvt, q, p, work, olwork, err)
     !! Forms the full M-by-M orthogonal matrix \(Q\) from the elementary
     !! reflectors returned by the base QR factorization algorithm.
     complex(real64), intent(inout), dimension(:,:) :: r
@@ -1271,7 +1271,7 @@ subroutine qr_rank1_update_dbl(q, r, u, v, work, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine qr_rank1_update_cmplx(q, r, u, v, work, rwork, err)
+subroutine qr_rank1_update_cmplx(q, r, u, v, work, rwork, err)
     !! Computes the rank-1 update to an M-by-N QR factored matrix \(A\) where
     !! \(M \ge N\), \(A = Q R\), and \(A_1 = A + \vec{u} \vec{v}^H\) such that
     !! \(A_1 = Q_1 R_1\).
@@ -1663,7 +1663,7 @@ subroutine solve_qr_no_pivot_vec(a, tau, b, work, olwork, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine solve_qr_no_pivot_vec_cmplx(a, tau, b, work, olwork, err)
+subroutine solve_qr_no_pivot_vec_cmplx(a, tau, b, work, olwork, err)
     !! Solves a system of M QR-factored equations of N unknowns.  M must be
     !! greater than or equal to N.
     complex(real64), intent(inout), dimension(:,:) :: a
@@ -2080,7 +2080,7 @@ subroutine solve_qr_pivot_mtx_cmplx(a, tau, jpvt, b, work, olwork, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine solve_qr_pivot_vec(a, tau, jpvt, b, work, olwork, err)
+subroutine solve_qr_pivot_vec(a, tau, jpvt, b, work, olwork, err)
     !! Solves a system of M QR-factored equations of N unknowns.
     real(real64), intent(inout), dimension(:,:) :: a
         !! On input, the M-by-N QR factored matrix as returned by qr_factor.  
@@ -2238,7 +2238,7 @@ module subroutine solve_qr_pivot_vec(a, tau, jpvt, b, work, olwork, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-module subroutine solve_qr_pivot_vec_cmplx(a, tau, jpvt, b, work, olwork, err)
+subroutine solve_qr_pivot_vec_cmplx(a, tau, jpvt, b, work, olwork, err)
     !! Solves a system of M QR-factored equations of N unknowns.
     complex(real64), intent(inout), dimension(:,:) :: a
         !! On input, the M-by-N QR factored matrix as returned by qr_factor.  
diff --git a/src/linalg_sparse.f90 b/src/linalg_sparse.f90
index 4eb42758..9b43a067 100644
--- a/src/linalg_sparse.f90
+++ b/src/linalg_sparse.f90
@@ -839,7 +839,7 @@ function csr_mtx_divide_scalar_1(a, b) result(rst)
 end function
 
 ! ------------------------------------------------------------------------------
-module function csr_transpose(a) result(rst)
+function csr_transpose(a) result(rst)
     !! Transposes a CSR matrix.
     class(csr_matrix), intent(in) :: a
         !! The CSR matrix.
diff --git a/src/linalg_svd.f90 b/src/linalg_svd.f90
index 4ec562e5..49a186d2 100644
--- a/src/linalg_svd.f90
+++ b/src/linalg_svd.f90
@@ -12,7 +12,7 @@ module linalg_svd
     end interface
 contains
 ! ------------------------------------------------------------------------------
-module subroutine svd_dbl(a, s, u, vt, work, olwork, err)
+subroutine svd_dbl(a, s, u, vt, work, olwork, err)
     !! Computes the singular value decomposition of an M-by-N matrix \(A\) such 
     !! that \(A = U S V^T\) where \(U\) is an M-by-M orthogonal matrix, \(S\)
     !! is an M-by-N diagonal matrix containing the singular values, and \(V\)
diff --git a/src/linalg_tri.f90 b/src/linalg_tri.f90
index 11861070..f51c7f2c 100644
--- a/src/linalg_tri.f90
+++ b/src/linalg_tri.f90
@@ -14,7 +14,7 @@ module linalg_tri
     end interface
 contains
 ! ------------------------------------------------------------------------------
-module subroutine solve_tri_mtx(lside, upper, trans, nounit, alpha, a, b, err)
+subroutine solve_tri_mtx(lside, upper, trans, nounit, alpha, a, b, err)
     !! Solves a triangular system of equations of the form 
     !! \(op(A) X = \alpha B\) or \(X op(A) = \alpha B\) where \(A\) is a 
     !! triangular matrix (either upper or lower) for the unknown \(X\).