diff --git a/src/gray_equil.f90 b/src/gray_equil.f90
index c1f1c5a..39e7101 100644
--- a/src/gray_equil.f90
+++ b/src/gray_equil.f90
@@ -35,7 +35,7 @@ module gray_equil
   end type
 
   abstract interface
-    subroutine pol_flux_sub(self, R, z, psi_n, dpsidr, dpsidz, &
+    pure subroutine pol_flux_sub(self, R, z, psi_n, dpsidr, dpsidz, &
                             ddpsidrr, ddpsidzz, ddpsidrz)
       ! Computes the normalised poloidal flux ψ_n and its
       ! derivatives wrt (R, z) up to the second order.
@@ -48,7 +48,7 @@ module gray_equil
         psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz
     end subroutine pol_flux_sub
 
-    subroutine pol_curr_sub(self, psi_n, fpol, dfpol)
+    pure subroutine pol_curr_sub(self, psi_n, fpol, dfpol)
       ! Computes the poloidal current function F(ψ_n)
       ! and (optionally) its derivative dF/dψ_n given ψ_n.
       import :: abstract_equil, wp_
@@ -58,7 +58,7 @@ module gray_equil
       real(wp_),             intent(out), optional :: dfpol  ! derivative
     end subroutine pol_curr_sub
 
-    function safety_fun(self, psi_n) result(q)
+    pure function safety_fun(self, psi_n) result(q)
       ! Computes the safety factor q as a function of the
       ! normalised poloidal flux ψ_n.
       !
@@ -69,7 +69,7 @@ module gray_equil
       real(wp_)                         :: q
     end function safety_fun
 
-    function rho_conv_fun(self, rho_in) result(rho_out)
+    pure function rho_conv_fun(self, rho_in) result(rho_out)
       ! Converts between poloidal (ρ_p) and toroidal (ρ_t) normalised radius
       import :: abstract_equil, wp_
       class(abstract_equil), intent(in) :: self
@@ -176,7 +176,7 @@ module gray_equil
 
 contains
 
-  subroutine b_field(self, R, z, B_R, B_z, B_phi)
+  pure subroutine b_field(self, R, z, B_R, B_z, B_phi)
     ! Computes the magnetic field as a function of
     ! (R, z) in cylindrical coordinates
     !
@@ -211,7 +211,7 @@ contains
   end subroutine b_field
 
 
-  function tor_curr(self, R, z) result(J_phi)
+  pure function tor_curr(self, R, z) result(J_phi)
     ! Computes the toroidal current J_φ as a function of (R, z)
     use const_and_precisions, only : mu0_
 
@@ -251,7 +251,7 @@ contains
   ! Analytical model
   !
 
-  subroutine analytic_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
+  pure subroutine analytic_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
                           ddpsidrr, ddpsidzz, ddpsidrz)
     use const_and_precisions, only : pi
 
@@ -376,7 +376,7 @@ contains
   end subroutine analytic_pol_flux
 
 
-  subroutine analytic_pol_curr(self, psi_n, fpol, dfpol)
+  pure subroutine analytic_pol_curr(self, psi_n, fpol, dfpol)
     ! subroutine arguments
     class(analytic_equil), intent(in)            :: self
     real(wp_),             intent(in)            :: psi_n
@@ -394,7 +394,7 @@ contains
   end subroutine analytic_pol_curr
 
 
-  function analytic_safety(self, psi_n) result(q)
+  pure function analytic_safety(self, psi_n) result(q)
     ! function arguments
     class(analytic_equil), intent(in) :: self
     real(wp_),             intent(in) :: psi_n
@@ -412,7 +412,7 @@ contains
   end function analytic_safety
 
 
-  function analytic_pol2tor(self, rho_in) result(rho_out)
+  pure function analytic_pol2tor(self, rho_in) result(rho_out)
     ! function arguments
     class(analytic_equil), intent(in) :: self
     real(wp_),             intent(in) :: rho_in
@@ -438,7 +438,7 @@ contains
   end function analytic_pol2tor
 
 
-  function analytic_tor2pol(self, rho_in) result(rho_out)
+  pure function analytic_tor2pol(self, rho_in) result(rho_out)
     ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius
     use const_and_precisions, only : comp_eps
     use minpack,              only : hybrj1
@@ -462,7 +462,7 @@ contains
 
     contains
 
-      subroutine equation(n, x, f, df, ldf, flag)
+      pure subroutine equation(n, x, f, df, ldf, flag)
         ! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0
 
         ! optimal step size
@@ -490,7 +490,7 @@ contains
   end function analytic_tor2pol
 
 
-  subroutine analytic_flux_contour(self, psi0, R_min, R, z, &
+  pure subroutine analytic_flux_contour(self, psi0, R_min, R, z, &
                                   R_hi, z_hi, R_lo, z_lo)
     use const_and_precisions, only : pi
 
@@ -525,7 +525,7 @@ contains
   ! Numerical equilibrium
   !
 
-  subroutine numeric_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
+  pure subroutine numeric_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
                           ddpsidrr, ddpsidzz, ddpsidrz)
     ! function arguments
     class(numeric_equil), intent(in)            :: self
@@ -554,7 +554,7 @@ contains
   end subroutine numeric_pol_flux
 
 
-  subroutine numeric_pol_curr(self, psi_n, fpol, dfpol)
+  pure subroutine numeric_pol_curr(self, psi_n, fpol, dfpol)
     ! subroutine arguments
     class(numeric_equil), intent(in)            :: self
     real(wp_),            intent(in)            :: psi_n
@@ -573,7 +573,7 @@ contains
   end subroutine numeric_pol_curr
   
 
-  function numeric_safety(self, psi_n) result(q)
+  pure function numeric_safety(self, psi_n) result(q)
     ! function arguments
     class(numeric_equil), intent(in) :: self
     real(wp_),            intent(in) :: psi_n
@@ -589,7 +589,7 @@ contains
   end function numeric_safety
 
 
-  function numeric_pol2tor(self, rho_in) result(rho_out)
+  pure function numeric_pol2tor(self, rho_in) result(rho_out)
     ! function arguments
     class(numeric_equil), intent(in) :: self
     real(wp_),            intent(in) :: rho_in
@@ -599,7 +599,7 @@ contains
   end function numeric_pol2tor
 
 
-  function numeric_tor2pol(self, rho_in) result(rho_out)
+  pure function numeric_tor2pol(self, rho_in) result(rho_out)
     ! function arguments
     class(numeric_equil), intent(in) :: self
     real(wp_),            intent(in) :: rho_in
@@ -956,7 +956,7 @@ contains
   end subroutine numeric_init
 
 
-  subroutine rescale_boundary(cont, psi, O, t0)
+  pure subroutine rescale_boundary(cont, psi, O, t0)
     ! Given the plasma boundary contour `cont`, the position of the
     ! magnetic axis `O`, and a scaling factor `t0`; this subroutine
     ! rescales the contour by `t0` about `O` while ensuring the
@@ -994,7 +994,7 @@ contains
 
   contains
 
-    function s(t)
+    pure function s(t)
       ! Rescriction of ψ(R, z) on the line Q(t) = O + tN
 
       real(wp_), intent(in) :: t
@@ -1047,7 +1047,7 @@ contains
 
   contains
 
-    subroutine equation(n, x, f, df, ldf, flag)
+    pure subroutine equation(n, x, f, df, ldf, flag)
       ! The equation to solve: f(R,z) = ∇ψ(R,z) = 0
 
       ! subroutine arguments
@@ -1108,7 +1108,7 @@ contains
     z1 = sol(2)
   contains
 
-    subroutine equation(n, x, f, df, ldf, flag)
+    pure subroutine equation(n, x, f, df, ldf, flag)
       ! The equation to solve: f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R] = 0
 
       ! subroutine arguments
@@ -1134,7 +1134,7 @@ contains
   ! Vacuum
   !
 
-  subroutine vacuum_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
+  pure subroutine vacuum_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
                           ddpsidrr, ddpsidzz, ddpsidrz)
     ! function arguments
     class(vacuum), intent(in)            :: self
@@ -1154,7 +1154,7 @@ contains
   end subroutine vacuum_pol_flux
 
 
-  subroutine vacuum_pol_curr(self, psi_n, fpol, dfpol)
+  pure subroutine vacuum_pol_curr(self, psi_n, fpol, dfpol)
     ! subroutine arguments
     class(vacuum), intent(in)            :: self
     real(wp_),     intent(in)            :: psi_n
@@ -1169,7 +1169,7 @@ contains
   end subroutine vacuum_pol_curr
 
 
-  function vacuum_safety(self, psi_n) result(q)
+  pure function vacuum_safety(self, psi_n) result(q)
     ! function arguments
     class(vacuum), intent(in) :: self
     real(wp_),     intent(in) :: psi_n
@@ -1182,7 +1182,7 @@ contains
   end function vacuum_safety
 
 
-  function vacuum_conv(self, rho_in) result(rho_out)
+  pure function vacuum_conv(self, rho_in) result(rho_out)
     ! function arguments
     class(vacuum), intent(in) :: self
     real(wp_),     intent(in) :: rho_in
@@ -1195,7 +1195,7 @@ contains
   end function vacuum_conv
 
 
-  subroutine vacuum_flux_contour(self, psi0, R_min, R, z, &
+  pure subroutine vacuum_flux_contour(self, psi0, R_min, R, z, &
                                  R_hi, z_hi, R_lo, z_lo)
     ! subroutine arguments
     class(vacuum), intent(in)    :: self
diff --git a/src/gray_tables.f90 b/src/gray_tables.f90
index 600d876..11a1faa 100644
--- a/src/gray_tables.f90
+++ b/src/gray_tables.f90
@@ -267,7 +267,7 @@ contains
     ! Build a regular (R, z) grid
     dR = (equil%r_range(2) - equil%r_range(1) - comp_eps)/(npsi - 1)
     dz = (equil%z_range(2) - equil%z_range(1))/(npsi - 1)
-    do j = 1, npsi
+    do concurrent (j = 1:npsi)
       R(j) = comp_eps + equil%r_range(1) + dR*(j - 1)
       z(j) = equil%z_range(1) + dz*(j - 1)
     end do
@@ -360,7 +360,7 @@ contains
 
   contains
 
-    subroutine equation(n, x, f, df, ldf, flag)
+    pure subroutine equation(n, x, f, df, ldf, flag)
       ! The equation to solve: f(x) = q(x) - q₀ = 0
       use const_and_precisions, only : comp_eps
 
diff --git a/src/splines.f90 b/src/splines.f90
index a3a83c6..308c7a1 100644
--- a/src/splines.f90
+++ b/src/splines.f90
@@ -108,7 +108,7 @@ contains
   end subroutine spline_simple_deinit
 
 
-  function spline_simple_eval(self, x) result(y)
+  pure function spline_simple_eval(self, x) result(y)
     ! Evaluates the spline at x
     use utils, only : locate
 
@@ -128,7 +128,7 @@ contains
   end function spline_simple_eval
 
 
-  function spline_simple_raw_eval(self, i, dx) result(y)
+  pure function spline_simple_raw_eval(self, i, dx) result(y)
     ! Evaluates the i-th polynomial of the spline at dx
 
     ! subroutine arguments
@@ -142,7 +142,7 @@ contains
   end function spline_simple_raw_eval
 
 
-  function spline_simple_deriv(self, x) result(y)
+  pure function spline_simple_deriv(self, x) result(y)
     ! Computes the derivative of the spline at x
     use utils, only : locate
 
@@ -163,7 +163,7 @@ contains
   end function spline_simple_deriv
 
 
-  subroutine spline_simple_coeffs(x, y, n, c)
+  pure subroutine spline_simple_coeffs(x, y, n, c)
     ! Computes the cubic coefficients of all n polynomials
 
     ! subroutine arguments
@@ -323,7 +323,7 @@ contains
   end subroutine spline_1d_deinit
 
 
-  function spline_1d_eval(self, x) result(y)
+  pure function spline_1d_eval(self, x) result(y)
     ! Evaluates the spline at x
     use dierckx, only : splev
 
@@ -341,7 +341,7 @@ contains
   end function spline_1d_eval
 
 
-  function spline_1d_deriv(self, x, order) result(y)
+  pure function spline_1d_deriv(self, x, order) result(y)
     ! Evaluates the spline n-th order derivative at x
     use dierckx, only : splder
 
@@ -525,7 +525,7 @@ contains
   end subroutine spline_2d_final
 
 
-  function spline_2d_eval(self, x, y) result(z)
+  pure function spline_2d_eval(self, x, y) result(z)
     ! Evaluates the spline at (x, y)
     use dierckx, only : fpbisp
 
@@ -585,7 +585,7 @@ contains
   end subroutine spline_2d_init_deriv
 
 
-  function spline_2d_deriv(self, x, y, n, m) result(z)
+  pure function spline_2d_deriv(self, x, y, n, m) result(z)
     ! Evaluates the spline n-th partial derivative w.r.t x
     ! and m-th partial derivative w.r.t y at (x, y)
     !
@@ -679,7 +679,7 @@ contains
   end subroutine linear_1d_deinit
 
 
-  function linear_1d_eval(self, x) result(y)
+  pure function linear_1d_eval(self, x) result(y)
     ! Evaluates the linear interpolated data at x
     use utils, only : locate
 
@@ -697,7 +697,7 @@ contains
   end function linear_1d_eval
 
 
-  function linear_1d_raw_eval(self, i, x) result(y)
+  pure function linear_1d_raw_eval(self, i, x) result(y)
     ! Performs the linear interpolation in the i-th interval
 
     ! subroutine arguments
diff --git a/src/types.f90 b/src/types.f90
index e6f13ab..7acf18c 100644
--- a/src/types.f90
+++ b/src/types.f90
@@ -70,7 +70,7 @@ module types
 
 contains
 
-  subroutine queue_put(self, val)
+  pure subroutine queue_put(self, val)
     ! Inserts an item of value `val` at the end of the `self` queue
     class(queue), intent(inout) :: self
     class(*),     intent(in)    :: val
@@ -260,7 +260,7 @@ contains
   end subroutine table_save
 
 
-  function contour_init(R, z) result(self)
+  pure function contour_init(R, z) result(self)
     ! Creates a contour
 
     ! functions arguments
@@ -277,7 +277,7 @@ contains
   end function contour_init
 
 
-  function contour_init_rect(Rmin, Rmax, zmin, zmax) result(self)
+  pure function contour_init_rect(Rmin, Rmax, zmin, zmax) result(self)
     ! Given two ranges [Rmin, Rmax], [zmin, zmax] creates a
     ! rectangular contour as follows:
     !
@@ -296,7 +296,7 @@ contains
   end function contour_init_rect
 
 
-  function contour_contains(self, R0, z0) result(inside)
+  pure function contour_contains(self, R0, z0) result(inside)
     ! Tests whether the point (`R`, `z`) lies inside the 2D contour
 
     use utils, only : intlinf, locate_unord
diff --git a/src/utils.f90 b/src/utils.f90
index e56d1c8..d6dfd8e 100644
--- a/src/utils.f90
+++ b/src/utils.f90
@@ -6,7 +6,7 @@ module utils
 
 contains
 
-  function locatef(a,n,x) result(j)
+  pure function locatef(a,n,x) result(j)
 ! Given an array a(n), and a value x, with a(n) monotonic, either
 ! increasing or decreasing, returns a value j such that
 ! a(j) < x <= a(j+1) for a increasing, and such that
@@ -32,7 +32,7 @@ contains
     j=jl
   end function locatef
 
-  subroutine locate(xx,n,x,j)
+  pure subroutine locate(xx,n,x,j)
     integer, intent(in) :: n
     real(wp_), intent(in) :: xx(n), x
     integer, intent(out) :: j
@@ -58,7 +58,7 @@ contains
     j=jl
   end subroutine locate
 
-  subroutine locatex(xx,n,n1,n2,x,j)
+  pure subroutine locatex(xx,n,n1,n2,x,j)
     integer, intent(in) :: n,n1,n2
     real(wp_), intent(in) :: xx(n), x
     integer, intent(out) :: j
@@ -84,7 +84,7 @@ contains
   end subroutine locatex
 
 
-  subroutine locate_unord(array, value, locs, n, nlocs)
+  pure subroutine locate_unord(array, value, locs, n, nlocs)
     ! Given an `array` of size `n` and a `value`, finds at most
     ! `n` locations `locs` such that `value` is between
     ! `array(locs(i))` and `array(locs(i+i))`, in whichever order.
@@ -117,7 +117,7 @@ contains
   end subroutine locate_unord
 
 
-  function intlinf(x1,y1,x2,y2,x) result(y)
+  pure function intlinf(x1,y1,x2,y2,x) result(y)
     !linear interpolation
     !must be x1 != x2
     use const_and_precisions, only : one
@@ -128,7 +128,7 @@ contains
     y=a*y1+(one-a)*y2
   end function intlinf
 
-  subroutine intlin(x1,y1,x2,y2,x,y)
+  pure subroutine intlin(x1,y1,x2,y2,x,y)
     real(wp_), intent(in) :: x1,y1,x2,y2,x
     real(wp_), intent(out) :: y
     real(wp_) :: dx,aa,bb
@@ -142,7 +142,7 @@ contains
     y=aa*y1+bb*y2
   end subroutine intlin
 
-  subroutine vmax(x,n,xmax,imx)
+  pure subroutine vmax(x,n,xmax,imx)
     integer, intent(in) :: n
     real(wp_), intent(in) :: x(n)
     real(wp_), intent(out) :: xmax
@@ -163,7 +163,7 @@ contains
     end do
   end subroutine vmax
 
-  subroutine vmin(x,n,xmin,imn)
+  pure subroutine vmin(x,n,xmin,imn)
     integer, intent(in) :: n
     real(wp_), intent(in) :: x(n)
     real(wp_), intent(out) :: xmin
@@ -184,7 +184,7 @@ contains
     end do
   end subroutine vmin
 
-  subroutine vmaxmini(x,n,xmin,xmax,imn,imx)
+  pure subroutine vmaxmini(x,n,xmin,xmax,imn,imx)
     integer, intent(in) :: n
     real(wp_), intent(in) :: x(n)
     real(wp_), intent(out) :: xmin, xmax
@@ -210,7 +210,7 @@ contains
     end do
   end subroutine vmaxmini
 
-  subroutine vmaxmin(x,n,xmin,xmax)
+  pure subroutine vmaxmin(x,n,xmin,xmax)
     integer, intent(in) :: n
     real(wp_), intent(in) :: x(n)
     real(wp_), intent(out) :: xmin, xmax
@@ -230,7 +230,7 @@ contains
     end do
   end subroutine vmaxmin
 
-  subroutine order(p,q)
+  pure subroutine order(p,q)
 ! returns p,q in ascending order
     real(wp_), intent(inout) :: p,q
     real(wp_) :: temp
@@ -241,7 +241,7 @@ contains
     end if
   end subroutine order
 
-  subroutine bubble(a,n)
+  pure subroutine bubble(a,n)
 ! bubble sorting of array a
     integer, intent(in) :: n
     real(wp_), dimension(n), intent(inout) :: a
@@ -254,7 +254,7 @@ contains
   end subroutine bubble
 
 
-  subroutine range2rect(xmin, xmax, ymin, ymax, x, y)
+  pure subroutine range2rect(xmin, xmax, ymin, ymax, x, y)
     ! Given two ranges [xmin, xmax], [ymin, ymax] builds
     ! the x and y vertices of the following rectangle:
     !
diff --git a/src/vendor/dierckx.f90 b/src/vendor/dierckx.f90
index d339d17..3775b2e 100644
--- a/src/vendor/dierckx.f90
+++ b/src/vendor/dierckx.f90
@@ -1242,7 +1242,7 @@ contains
     end do
   end subroutine fpback
 
-  subroutine fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wx,wy,lx,ly)
+  pure subroutine fpbisp(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wx,wy,lx,ly)
 ! arguments
     integer, intent(in) :: nx, ny, kx, ky, mx, my
     integer, intent(out) :: lx(mx), ly(my)
@@ -1321,7 +1321,7 @@ contains
     end do
   end subroutine fpbisp
 
-  subroutine fpbspl(t,n,k,x,l,h)
+  pure subroutine fpbspl(t,n,k,x,l,h)
 !  subroutine fpbspl evaluates the (k+1) non-zero b-splines of
 !  degree k at t(l) <= x < t(l+1) using the stable recurrence
 !  relation of de boor and cox.
@@ -3942,7 +3942,7 @@ contains
     end do
   end subroutine fpcurf
 
-  subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
+  pure subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
 !  subroutine splder evaluates in a number of points x(i),i=1,2,...,m
 !  the derivative of order nu of a spline s(x) of degree k,given in
 !  its b-spline representation.
@@ -4081,7 +4081,7 @@ contains
     end do
   end subroutine splder
 
-  subroutine splev(t,n,c,k,x,y,m,ier)
+  pure subroutine splev(t,n,c,k,x,y,m,ier)
 !  subroutine splev evaluates in a number of points x(i),i=1,2,...,m
 !  a spline s(x) of degree k, given in its b-spline representation.
 !
diff --git a/src/vendor/minpack.f90 b/src/vendor/minpack.f90
index c230b48..e5a65d8 100644
--- a/src/vendor/minpack.f90
+++ b/src/vendor/minpack.f90
@@ -47,7 +47,7 @@ module minpack
 
 contains
 
-  subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa)
+  pure subroutine hybrj1(fcn,n,x,fvec,fjac,ldfjac,tol,info,wa,lwa)
     use const_and_precisions, only : zero, one
 !   arguments
     integer, intent(in) :: n, ldfjac, lwa
@@ -154,7 +154,7 @@ contains
     real(wp_), parameter :: factor=1.0e2_wp_
     
     interface
-      subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
+      pure subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
         use const_and_precisions, only : wp_
         integer, intent(in) :: n,ldfjac,iflag
         real(wp_), intent(in) :: x(n)
@@ -185,7 +185,7 @@ contains
     if (info == 5) info = 4
   end subroutine hybrj1
 
-  subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, &
+  pure subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, &
                      factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, &
                      wa3,wa4)
     use const_and_precisions,          only : zero, one, epsmch=>comp_eps
@@ -353,7 +353,7 @@ contains
                p001 = 1.0e-3_wp_, p0001 = 1.0e-4_wp_
     
     interface
-      subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
+      pure subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
         use const_and_precisions, only : wp_
         integer, intent(in) :: n,ldfjac,iflag
         real(wp_), intent(in) :: x(n)
@@ -637,7 +637,7 @@ contains
 
   end subroutine hybrj
 
-  subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
+  pure subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
     use const_and_precisions, only : zero, one, epsmch=>comp_eps
 !   arguments
     integer, intent(in) :: n, lr
@@ -805,7 +805,7 @@ contains
     end do
   end subroutine dogleg
 
-  function enorm(n,x)
+  pure function enorm(n,x)
     use const_and_precisions, only : zero, one
     real(wp_) :: enorm
     integer, intent(in) :: n
@@ -903,7 +903,7 @@ contains
     end if
   end function enorm
 
-  subroutine qform(m,n,q,ldq,wa)
+  pure subroutine qform(m,n,q,ldq,wa)
     use const_and_precisions, only : zero, one
 !   arguments
     integer, intent(in) :: m,n,ldq
@@ -994,7 +994,7 @@ contains
     end do
   end subroutine qform
 
-  subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
+  pure subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
     use const_and_precisions, only : zero, one, epsmch=>comp_eps
 !   arguments
     integer, intent(in) :: m, n, lda, lipvt
@@ -1156,7 +1156,7 @@ contains
     end do
   end subroutine qrfac
 
-  subroutine r1mpyq(m,n,a,lda,v,w)
+  pure subroutine r1mpyq(m,n,a,lda,v,w)
     use const_and_precisions, only : one
 !   arguments
     integer, intent(in) :: m, n, lda
@@ -1247,7 +1247,7 @@ contains
     end do
   end subroutine r1mpyq
 
-  subroutine r1updt(m,n,s,ls,u,v,w,sing)
+  pure subroutine r1updt(m,n,s,ls,u,v,w,sing)
     use const_and_precisions, only : zero, one, giant=>comp_huge
 !   arguments
     integer, intent(in) :: m, n, ls