some versions of ifort have problems with 'do concurrent'

src/grid/spectral_utilities.f90

@@ -376,20 +376,32 @@ subroutine utilities_updateGamma(C)
     !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err)
     do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells1Red
       if (any([i,j,k] /= 1)) then                                                                   ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
-        do concurrent (l = 1:3, m = 1:3)
+#ifndef __INTEL_COMPILER
+        do concurrent(l = 1:3, m = 1:3)
           xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
         end do
         do concurrent(l = 1:3, m = 1:3)
           temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
         end do
+#else
+        forall(l = 1:3, m = 1:3) &
+          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset)
+        forall(l = 1:3, m = 1:3) &
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+#endif
         A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) =  temp33_cmplx%re
         A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
         if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
           call math_invert(A_inv, err, A)
           temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
+#ifndef __INTEL_COMPILER
           do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
             gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
           end do
+#else
+          forall(l=1:3, m=1:3, n=1:3, o=1:3) &
+            gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_cmplx(l,n) * xiDyad_cmplx(o,m)
+#endif
         end if
       end if
     end do; end do; end do
@@ -497,23 +509,37 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
     !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat)
     do k = 1, cells3; do j = 1, cells(2); do i = 1, cells1Red
       if (any([i,j,k+cells3Offset] /= 1)) then                                                      ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
+#ifndef __INTEL_COMPILER
         do concurrent(l = 1:3, m = 1:3)
           xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
         end do
         do concurrent(l = 1:3, m = 1:3)
           temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
         end do
+#else
+        forall(l = 1:3, m = 1:3) &
+          xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k)
+        forall(l = 1:3, m = 1:3) &
+          temp33_cmplx(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx)
+#endif
         A(1:3,1:3) = temp33_cmplx%re; A(4:6,4:6) =  temp33_cmplx%re
         A(1:3,4:6) = temp33_cmplx%im; A(4:6,1:3) = -temp33_cmplx%im
         if (abs(math_det33(A(1:3,1:3))) > 1e-16) then
           call math_invert(A_inv, err, A)
           temp33_cmplx = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal)
+#ifndef __INTEL_COMPILER
           do concurrent(l=1:3, m=1:3, n=1:3, o=1:3)
             gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m)
           end do
           do concurrent(l = 1:3, m = 1:3)
             temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
           end do
+#else
+          forall(l=1:3, m=1:3, n=1:3, o=1:3) &
+            gamma_hat(l,m,n,o,1,1,1) = temp33_cmplx(l,n)*xiDyad_cmplx(o,m)
+          forall(l = 1:3, m = 1:3) &
+            temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,1,1,1)*tensorField_fourier(1:3,1:3,i,j,k))
+#endif
           tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx
         else
           tensorField_fourier(1:3,1:3,i,j,k) = cmplx(0.0_pReal,0.0_pReal,pReal)
@@ -524,9 +550,14 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
   else memoryEfficient
     !$OMP PARALLEL DO PRIVATE(l,m,temp33_cmplx)
     do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells1Red
+#ifndef __INTEL_COMPILER
       do concurrent(l = 1:3, m = 1:3)
         temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k))
       end do
+#else
+      forall(l = 1:3, m = 1:3) &
+        temp33_cmplx(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k)*tensorField_fourier(1:3,1:3,i,j,k))
+#endif
       tensorField_fourier(1:3,1:3,i,j,k) = temp33_cmplx
     end do; end do; end do
     !$OMP END PARALLEL DO
@@ -758,7 +789,7 @@ subroutine utilities_fourierVectorDivergence()
 
 
   scalarField_fourier(1:cells1Red,1:cells(2),1:cells3) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(2),1:cells3) &
-                                                             *conjg(-xi1st))
+                                                             *conjg(-xi1st),1)
 
 end subroutine utilities_fourierVectorDivergence
 

src/math.f90

@@ -262,9 +262,8 @@ pure function math_identity4th()
     math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
   enddo
 #else
-  do i=1,3; do j=1,3; do k=1,3; do l=1,3
+  forall(i=1:3, j=1:3, k=1:3, l=1:3) &
     math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
-  enddo; enddo; enddo; enddo
 #endif
 
 end function math_identity4th
@@ -338,9 +337,7 @@ pure function math_outer(A,B)
     math_outer(i,j) = A(i)*B(j)
   enddo
 #else
-  do i=1,size(A,1); do j=1,size(B,1)
+  forall(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j)
-    math_outer(i,j) = A(i)*B(j)
-  enddo; enddo
 #endif
 
 end function math_outer
@@ -387,9 +384,7 @@ pure function math_mul3333xx33(A,B)
     math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
   enddo
 #else
-  do i=1,3; do j=1,3
+  forall (i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
-    math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
-  enddo; enddo
 #endif
 
 end function math_mul3333xx33
@@ -411,9 +406,7 @@ pure function math_mul3333xx3333(A,B)
     math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
   enddo
 #else
-  do i=1,3; do j=1,3; do k=1,3; do l=1,3
+  forall(i=1:3, j=1:3, k=1:3, l=1:3) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
-    math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
-  enddo; enddo; enddo; enddo
 #endif
 
 end function math_mul3333xx3333
@@ -752,9 +745,7 @@ pure function math_3333to99(m3333)
     math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
   enddo
 #else
-  do i=1,9; do j=1,9
+  forall(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
-    math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
-  enddo; enddo
 #endif
 
 end function math_3333to99
@@ -775,9 +766,7 @@ pure function math_99to3333(m99)
     math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
   enddo
 #else
-  do i=1,9; do j=1,9
+  forall(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
-    math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
-  enddo; enddo
 #endif
 
 end function math_99to3333
@@ -810,9 +799,7 @@ pure function math_sym3333to66(m3333,weighted)
     math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
   enddo
 #else
-  do i=1,6; do j=1,6
+  forall(i=1:6, j=1:6) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
-    math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
-  enddo; enddo
 #endif
 
 end function math_sym3333to66
@@ -950,9 +937,7 @@ pure function math_3333toVoigt66_stiffness(C) result(C_tilde)
     C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
   end do
 #else
-  do i=1,6; do j=1,6
+  forall(i=1:6, j=1:6) C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
-    C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
-  end do; end do
 #endif
 
 end function math_3333toVoigt66_stiffness

src/phase_mechanical_plastic.f90

@@ -379,10 +379,9 @@ module function plastic_deltaState(ph, en) result(broken)
     en
   logical :: broken
 
-  real(pReal),               dimension(3,3) :: &
+  real(pReal), dimension(3,3) :: &
     Mp
   integer :: &
-    myOffset, &
     mySize