diff --git a/src/gray_plasma.f90 b/src/gray_plasma.f90
index a3f79e4..c44cb19 100644
--- a/src/gray_plasma.f90
+++ b/src/gray_plasma.f90
@@ -20,6 +20,39 @@ module gray_plasma
       procedure(zeff_fun),    deferred :: zeff
   end type
 
+  abstract interface
+    subroutine density_sub(self, psin, dens, ddens)
+      ! Computes the density its first derivative as a function of
+      ! normalised poloidal flux.
+      !
+      ! Note: density has units of 10¹⁹ m⁻³.
+      import :: abstract_plasma, wp_
+      class(abstract_plasma), intent(in)  :: self
+      real(wp_),              intent(in)  :: psin         ! normalised poloidal flux
+      real(wp_),              intent(out) :: dens, ddens  ! density and first derivative
+    end subroutine density_sub
+
+    function temp_fun(self, psin) result(temp)
+      ! Computes the temperature as a function of the
+      ! normalised poloidal flux.
+      !
+      ! Note: temperature has units of keV.
+      import :: abstract_plasma, wp_
+      class(abstract_plasma), intent(in) :: self
+      real(wp_),              intent(in) :: psin
+      real(wp_)                          :: temp
+    end function temp_fun
+
+    function zeff_fun(self, psin) result(zeff)
+      ! Computes the effective charge Z_eff as a
+      ! function of the normalised poloidal flux.
+      import :: abstract_plasma, wp_
+      class(abstract_plasma), intent(in) :: self
+      real(wp_),              intent(in) :: psin
+      real(wp_)                          :: zeff
+    end function zeff_fun
+  end interface
+
   ! Analytical plasma description
   type, extends(abstract_plasma) :: analytic_plasma
     private
@@ -57,48 +90,6 @@ module gray_plasma
     procedure :: zeff    => numeric_zeff
   end type
 
-  abstract interface
-    subroutine density_sub(self, psin, dens, ddens)
-      ! Computes the density its first derivative as a function of
-      ! normalised poloidal flux.
-      !
-      ! Note: density has units of 10¹⁹ m⁻³.
-      import :: abstract_plasma, wp_
-      class(abstract_plasma), intent(in)  :: self
-      real(wp_),              intent(in)  :: psin         ! normalised poloidal flux
-      real(wp_),              intent(out) :: dens, ddens  ! density and first derivative
-    end subroutine density_sub
-
-    function temp_fun(self, psin) result(temp)
-      ! Computes the temperature as a function of the
-      ! normalised poloidal flux.
-      !
-      ! Note: temperature has units of keV.
-      import :: abstract_plasma, wp_
-      class(abstract_plasma), intent(in) :: self
-      real(wp_),              intent(in) :: psin
-      real(wp_)                          :: temp
-    end function temp_fun
-
-    function zeff_fun(self, psin) result(zeff)
-      ! Computes the effective charge Z_eff as a
-      ! function of the normalised poloidal flux.
-      import :: abstract_plasma, wp_
-      class(abstract_plasma), intent(in) :: self
-      real(wp_),              intent(in) :: psin
-      real(wp_)                          :: zeff
-    end function zeff_fun
-  end interface
-
-  ! Parameters of the analytical profiles model
-  type analytic_model
-    real(wp_) :: dens0     ! Density scaling factor
-    real(wp_) :: n1, n2    ! Density exponents
-    real(wp_) :: te0, te1  ! Temperature at ψ=0, ψ=1
-    real(wp_) :: t1, t2    ! Temperature exponents
-    real(wp_) :: zeff      ! Effective charge
-  end type
-
   private
   public abstract_plasma                  ! The abstract plasma object
   public analytic_plasma, numeric_plasma  ! Implementations
@@ -168,7 +159,7 @@ contains
     use logger, only : log_error
 
     ! subroutine arguments
-    class(numeric_plasma), intent(in) :: self
+    class(numeric_plasma), intent(in)  :: self
     real(wp_),             intent(in)  :: psin
     real(wp_),             intent(out) :: dens, ddens
 
diff --git a/src/splines.f90 b/src/splines.f90
index 7d0b431..94f3825 100644
--- a/src/splines.f90
+++ b/src/splines.f90
@@ -559,12 +559,10 @@ contains
       ! first m elements contain the spline coefficients
       m = (self%nknots_x - 3 - 1) * (self%nknots_y - 3 - 1)
 
-      do i = 0, size(self%partial, dim=1) - 1
-      do j = 0, size(self%partial, dim=2) - 1
-        if (associated(self%partial(i, j)%ptr)) then
-          self%partial(i, j)%ptr(1:m) = a * self%partial(i, j)%ptr(1:m)
-        end if
-      end do
+      do concurrent (i = 0:size(self%partial, dim=1) - 1, &
+                     j = 0:size(self%partial, dim=2) - 1, &
+                     associated(self%partial(i, j)%ptr))
+        self%partial(i, j)%ptr(1:m) = a * self%partial(i, j)%ptr(1:m)
       end do
 
     end if