broken Intel compiler (do concurrent)

errors occur under unclear conditions, better avoid at the moment:

https://community.intel.com/t5/Intel-Fortran-Compiler/do-concurrent-broken-with-openmp/m-p/1301528#M156942
https://community.intel.com/t5/Intel-Fortran-Compiler/Bug-with-do-concurrent-and-openmp/m-p/1173473m

src/element.f90

@@ -151,6 +151,7 @@ module element
   integer, dimension(NIPNEIGHBOR(CELLTYPE(1)),NIP(1)),   parameter :: IPNEIGHBOR1 = &
     reshape([&
       -2,-3,-1  &
+! Note: This fix is for gfortran 9 only. gfortran 8 supports neither, gfortran > 9 both variants
 #if !defined(__GFORTRAN__)
     ],shape(IPNEIGHBOR1))
 #else

src/math.f90

@@ -90,12 +90,12 @@ subroutine math_init
   integer, dimension(:), allocatable :: randInit
   class(tNode), pointer :: &
     num_generic
- 
+
   print'(/,a)', ' <<<+-  math init  -+>>>'; flush(IO_STDOUT)
 
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
   randomSeed  = num_generic%get_asInt('random_seed', defaultVal = 0)
- 
+
   call random_seed(size=randSize)
   allocate(randInit(randSize))
   if (randomSeed > 0) then
@@ -260,9 +260,15 @@ pure function math_identity4th()
   integer :: i,j,k,l
 
 
+#ifndef __INTEL_COMPILER
   do concurrent(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
+#else
+  do i=1,3; do j=1,3; do k=1,3; do 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
 
@@ -330,9 +336,15 @@ pure function math_outer(A,B)
   integer :: i,j
 
 
+#ifndef __INTEL_COMPILER
   do concurrent(i=1:size(A,1), j=1:size(B,1))
     math_outer(i,j) = A(i)*B(j)
   enddo
+#else
+  do i=1,size(A,1); do j=1,size(B,1)
+    math_outer(i,j) = A(i)*B(j)
+  enddo; enddo
+#endif
 
 end function math_outer
 
@@ -373,9 +385,15 @@ pure function math_mul3333xx33(A,B)
   integer :: i,j
 
 
+#ifndef __INTEL_COMPILER
   do concurrent(i=1:3, j=1:3)
     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
+    math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
+  enddo; enddo
+#endif
 
 end function math_mul3333xx33
 
@@ -390,9 +408,16 @@ pure function math_mul3333xx3333(A,B)
   real(pReal), dimension(3,3,3,3), intent(in) :: B
   real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333
 
+
+#ifndef __INTEL_COMPILER
   do concurrent(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))
   enddo
+#else
+  do i=1,3; do j=1,3; do k=1,3; do l=1,3
+    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
 
@@ -725,9 +750,15 @@ pure function math_3333to99(m3333)
   integer :: i,j
 
 
+#ifndef __INTEL_COMPILER
   do concurrent(i=1:9, j=1:9)
     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
+    math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
+  enddo; enddo
+#endif
 
 end function math_3333to99
 
@@ -742,10 +773,15 @@ pure function math_99to3333(m99)
 
   integer :: i,j
 
-
+#ifndef __INTEL_COMPILER
   do concurrent(i=1:9, j=1:9)
     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
+    math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
+  enddo; enddo
+#endif
 
 end function math_99to3333
 
@@ -772,9 +808,15 @@ pure function math_sym3333to66(m3333,weighted)
     w = NRMMANDEL
   endif
 
+#ifndef __INTEL_COMPILER
   do concurrent(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))
   enddo
+#else
+  do i=1,6; do j=1,6
+    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
 
@@ -874,10 +916,11 @@ subroutine math_eigh(w,v,error,m)
   real(pReal), dimension(size(m,1)),            intent(out) :: w                                    !< eigenvalues
   real(pReal), dimension(size(m,1),size(m,1)),  intent(out) :: v                                    !< eigenvectors
   logical,                                      intent(out) :: error
-  
+
   integer :: ierr
   real(pReal), dimension(size(m,1)**2) :: work
 
+
   v = m                                                                                             ! copy matrix to input (doubles as output) array
   call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr)
   error = (ierr /= 0)