diff --git a/grint_module.f90 b/grint_module.f90
index 2a868ad..57601c1 100644
--- a/grint_module.f90
+++ b/grint_module.f90
@@ -155,7 +155,9 @@ END FUNCTION func_derivs
        IF ( change > 0.0 ) THEN ! Better fit
 
           c = c_new
-          FORALL ( t = 1:nterms ) sigma(t) = SQRT ( array(t,t) / alpha(t,t) )
+          DO t = 1, nterms
+             sigma(t) = SQRT ( array(t,t) / alpha(t,t) )
+          END DO
 
           IF ( verbose ) THEN
              WRITE ( unit=output_unit, fmt='(i5)', advance='no' ) iter
diff --git a/initialize_module.f90 b/initialize_module.f90
index 7199992..cb0feba 100644
--- a/initialize_module.f90
+++ b/initialize_module.f90
@@ -302,12 +302,12 @@ SUBROUTINE chain_velocities ( temperature, constraints )
     ! in which case both these quantities are conserved
     ! NB there is at present no check for a singular inertia tensor in the angular momentum fix!
     ! We assume centre of mass is already at the origin
-    ! We assume unit molecular mass and employ Lennard-Jones units
+    ! We assume unit atom mass and employ Lennard-Jones units
 
-    REAL                 :: temp
+    REAL                 :: temp, r_sq
     REAL, DIMENSION(3)   :: v_cm, r_cm, ang_mom, ang_vel
     REAL, DIMENSION(3,3) :: inertia
-    INTEGER              :: i, xyz
+    INTEGER              :: i
     REAL, PARAMETER      :: tol = 1.e-6
 
     WRITE ( unit=output_unit, fmt='(a,t40,f15.6)' ) 'Chain velocities at temperature', temperature
@@ -331,7 +331,10 @@ SUBROUTINE chain_velocities ( temperature, constraints )
     DO i = 1, n
        ang_mom = ang_mom + cross_product ( r(:,i), v(:,i) )
        inertia = inertia - outer_product ( r(:,i), r(:,i) )
-       FORALL ( xyz=1:3 ) inertia(xyz,xyz) = inertia(xyz,xyz) + DOT_PRODUCT ( r(:,i), r(:,i) )
+       r_sq = DOT_PRODUCT ( r(:,i), r(:,i) )
+       inertia(1,1) = inertia(1,1) + r_sq
+       inertia(2,2) = inertia(2,2) + r_sq
+       inertia(3,3) = inertia(3,3) + r_sq
     END DO
 
     ! Solve linear system to get angular velocity
diff --git a/maths_module.f90 b/maths_module.f90
index 4d58e51..437102c 100644
--- a/maths_module.f90
+++ b/maths_module.f90
@@ -802,7 +802,7 @@ FUNCTION nematic_order ( e ) RESULT ( order )
 
     ! Note that this is not the same as the order parameter characterizing a crystal
 
-    INTEGER              :: i, n
+    INTEGER              :: n
     REAL, DIMENSION(3,3) :: q         ! order tensor
     REAL                 :: h, g, psi ! used in eigenvalue calculation
 
@@ -813,9 +813,12 @@ FUNCTION nematic_order ( e ) RESULT ( order )
     n = SIZE(e,dim=2)
 
     ! Order tensor: outer product of each orientation vector, summed over n molecules
+    ! Then normalized and made traceless
     q = SUM ( SPREAD ( e, dim=2, ncopies=3) * SPREAD ( e, dim=1, ncopies=3 ), dim = 3 )
-    q = 1.5 * q / REAL(n)                ! Normalize
-    FORALL (i=1:3) q(i,i) = q(i,i) - 0.5 ! Make traceless
+    q = 1.5 * q / REAL(n)
+    q(1,1) = q(1,1) - 0.5
+    q(2,2) = q(2,2) - 0.5
+    q(3,3) = q(3,3) - 0.5
 
     ! Trigonometric solution of characteristic cubic equation, assuming real roots
 
diff --git a/t_tensor.f90 b/t_tensor.f90
index ab4300a..d2236d8 100644
--- a/t_tensor.f90
+++ b/t_tensor.f90
@@ -55,7 +55,7 @@ PROGRAM t_tensor
   REAL, DIMENSION(3,3,3,3,3) :: tt5
 
   REAL    :: r12_mag, c1, c2, c12, v12t, v12e
-  INTEGER :: i, ioerr
+  INTEGER :: ioerr
   REAL    :: d_min, d_max, mu1_mag, mu2_mag, quad1_mag, quad2_mag
 
   NAMELIST /nml/ d_min, d_max, mu1_mag, mu2_mag, quad1_mag, quad2_mag
@@ -116,11 +116,15 @@ PROGRAM t_tensor
   mu2 = mu2_mag * e2
 
   ! Quadrupole tensors in space-fixed frame (traceless)
-  quad1 = 1.5 * outer_product ( e1, e1 )        
-  FORALL (i=1:3) quad1(i,i) = quad1(i,i) - 0.5
+  quad1 = 1.5 * outer_product ( e1, e1 )
+  quad1(1,1) = quad1(1,1) - 0.5
+  quad1(2,2) = quad1(2,2) - 0.5
+  quad1(3,3) = quad1(3,3) - 0.5
   quad1 = quad1_mag * quad1
   quad2 = 1.5 * outer_product ( e2, e2 )
-  FORALL (i=1:3) quad2(i,i) = quad2(i,i) - 0.5
+  quad2(1,1) = quad2(1,1) - 0.5
+  quad2(2,2) = quad2(2,2) - 0.5
+  quad2(3,3) = quad2(3,3) - 0.5
   quad2 = quad2_mag * quad2
 
   ! The T tensors of each rank: T2, T3, T4, T5
@@ -259,11 +263,12 @@ FUNCTION t2_tensor ( r, r3 ) RESULT ( t2 )
     REAL, DIMENSION(3),   INTENT(in)  :: r  ! Unit vector from 2 to 1
     REAL,                 INTENT(in)  :: r3 ! Third power of the modulus of r12
 
-    INTEGER :: i
-
     t2 = 3.0 * outer_product ( r, r ) ! Starting point
 
-    FORALL (i=1:3) t2(i,i) = t2(i,i) - 1.0 ! Make traceless
+    ! Make traceless
+    t2(1,1) = t2(1,1) - 1.0
+    t2(2,2) = t2(2,2) - 1.0
+    t2(3,3) = t2(3,3) - 1.0
 
     t2 = t2 / r3 ! Scale by third power of distance