following style guide

src/Marc/DAMASK_Marc.f90

@@ -70,7 +70,7 @@ subroutine DAMASK_interface_init
   if (ierr /= 0) then
     print*, 'working directory "'//trim(wd)//'" does not exist'
     call quit(1)
-  endif
+  end if
   symmetricSolver = solverIsSymmetric()
 
 end subroutine DAMASK_interface_init
@@ -111,8 +111,8 @@ logical function solverIsSymmetric()
         s = s + verify(line(s+1:),' ')                                                              ! start of second chunk
         e = s + scan  (line(s+1:),' ')                                                              ! end of second chunk
       solverIsSymmetric = line(s:e) /= '1'
-    endif
+    end if
-  enddo
+  end do
 100 close(fileUnit)
   contains
 
@@ -134,7 +134,7 @@ logical function solverIsSymmetric()
       lc(i:i) = string(i:i)
       n = index(UPPER,lc(i:i))
       if (n/=0) lc(i:i) = LOWER(n:n)
-    enddo
+    end do
   end function lc
 
 end function solverIsSymmetric
@@ -299,7 +299,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
                                   transpose(ffn)
     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
                                   transpose(ffn1)
-  endif
+  end if
 
   defaultNumThreadsInt = omp_get_num_threads()                                                      ! remember number of threads set by Marc
   call omp_set_num_threads(1_pI32)                                                                  ! no openMP
@@ -309,7 +309,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
     call materialpoint_initAll
     debug_Marc => config_debug%get_list('Marc',defaultVal=emptyList)
     debug_basic = debug_Marc%contains('basic')
-  endif
+  end if
 
   computationMode = 0                                                                               ! save initialization value, since it does not result in any calculation
   if (lovl == 4 ) then                                                                              ! jacobian requested by marc
@@ -333,35 +333,35 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
         lastIncConverged = .true.
         outdatedByNewInc = .true.
         print'(a,i6,1x,i2)', '<< HYPELA2 >> new increment..! ',m(1),nn
-      endif
+      end if
     else if ( timinc < theDelta ) then                                                              ! >> cutBack <<
       lastIncConverged = .false.
       outdatedByNewInc = .false.
       terminallyIll = .false.
       cycleCounter = -1                                                                             ! first calc step increments this to cycle = 0
       print'(a,i6,1x,i2)', '<< HYPELA2 >> cutback detected..! ',m(1),nn
-    endif                                                                                           ! convergence treatment end
+    end if                                                                                          ! convergence treatment end
     flush(6)
 
     if (lastLovl /= lovl) then
       cycleCounter  = cycleCounter + 1
       !mesh_cellnode = mesh_build_cellnodes()                                                       ! update cell node coordinates
       !call mesh_build_ipCoordinates()                                                              ! update ip coordinates
-    endif
+    end if
     if (outdatedByNewInc) then
       computationMode = ior(computationMode,materialpoint_AGERESULTS)
       outdatedByNewInc = .false.
-    endif
+    end if
     if (lastIncConverged) then
       computationMode = ior(computationMode,materialpoint_BACKUPJACOBIAN)
       lastIncConverged = .false.
-    endif
+    end if
 
     theTime  = cptim
     theDelta = timinc
     theInc   = inc
 
-  endif
+  end if
   lastLovl = lovl
 
   call materialpoint_general(computationMode,ffn,ffn1,t(1),timinc,int(m(1)),int(nn),stress,ddsdde)
@@ -429,13 +429,13 @@ subroutine uedinc(inc,incsub)
       if (discretization_Marc_FEM2DAMASK_node(n) /= -1) then
         call nodvar(1,n,d_n(1:3,discretization_Marc_FEM2DAMASK_node(n)),nqncomp,nqdatatype)
         if(nqncomp == 2) d_n(3,discretization_Marc_FEM2DAMASK_node(n)) = 0.0_pReal
-      endif
+      end if
-    enddo
+    end do
 
     call discretization_Marc_UpdateNodeAndIpCoords(d_n)
     call materialpoint_results(int(inc),cptim)
 
     inc_written = int(inc)
-  endif
+  end if
 
 end subroutine uedinc

src/Marc/discretization_Marc.f90

@@ -275,8 +275,8 @@ subroutine inputRead_fileFormat(fileFormat,fileContent)
     if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'version') then
       fileFormat = IO_intValue(fileContent(l),chunkPos,2)
       exit
-    endif
+    end if
-  enddo
+  end do
 
 end subroutine inputRead_fileFormat
 
@@ -302,8 +302,8 @@ subroutine inputRead_tableStyles(initialcond,hypoelastic,fileContent)
       initialcond = IO_intValue(fileContent(l),chunkPos,4)
       hypoelastic = IO_intValue(fileContent(l),chunkPos,5)
       exit
-    endif
+    end if
-  enddo
+  end do
 
 end subroutine inputRead_tableStyles
 
@@ -331,16 +331,16 @@ subroutine inputRead_matNumber(matNumber, &
         data_blocks = IO_intValue(fileContent(l+1),chunkPos,1)
       else
         data_blocks = 1
-      endif
+      end if
       allocate(matNumber(data_blocks), source = 0)
       do i = 0, data_blocks - 1
         j = i*(2+tableStyle) + 1
         chunkPos = IO_stringPos(fileContent(l+1+j))
         matNumber(i+1) = IO_intValue(fileContent(l+1+j),chunkPos,1)
-      enddo
+      end do
       exit
-    endif
+    end if
-  enddo
+  end do
 
 end subroutine inputRead_matNumber
 
@@ -368,8 +368,8 @@ subroutine inputRead_NnodesAndElements(nNodes,nElems,&
     elseif(IO_lc(IO_StringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
       chunkPos = IO_stringPos(fileContent(l+1))
       nNodes = IO_IntValue (fileContent(l+1),chunkPos,2)
-    endif
+    end if
-  enddo
+  end do
 
 end subroutine inputRead_NnodesAndElements
 
@@ -411,12 +411,12 @@ subroutine inputRead_NelemSets(nElemSets,maxNelemInSet,&
           if(IO_lc(IO_stringValue(fileContent(l+i),chunkPos,chunkPos(1))) /= 'c') then              ! line finished, read last value
             elemInCurrentSet = elemInCurrentSet + 1                                                 ! data ended
             exit
-          endif
+          end if
-        enddo
+        end do
-      endif
+      end if
       maxNelemInSet = max(maxNelemInSet, elemInCurrentSet)
-    endif
+    end if
-  enddo
+  end do
 
 end subroutine inputRead_NelemSets
 
@@ -448,8 +448,8 @@ subroutine inputRead_mapElemSets(nameElemSet,mapElemSet,&
        elemSet = elemSet+1
        nameElemSet(elemSet)  = trim(IO_stringValue(fileContent(l),chunkPos,4))
        mapElemSet(:,elemSet) = continuousIntValues(fileContent(l+1:),size(mapElemSet,1)-1,nameElemSet,mapElemSet,size(nameElemSet))
-    endif
+    end if
-  enddo
+  end do
 
 end subroutine inputRead_mapElemSets
 
@@ -484,17 +484,17 @@ subroutine inputRead_mapElems(FEM2DAMASK, &
           j = j + 1
           chunkPos = IO_stringPos(fileContent(l+1+i+j))
           nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
-        enddo
+        end do
-      enddo
+      end do
       exit
-    endif
+    end if
-  enddo
+  end do
 
   call math_sort(map_unsorted)
   allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1)
   do i = 1, nElems
     FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i)
-  enddo
+  end do
 
 end subroutine inputRead_mapElems
 
@@ -522,16 +522,16 @@ subroutine inputRead_mapNodes(FEM2DAMASK, &
       chunkPos = [1,1,10]
       do i = 1,nNodes
         map_unsorted(:,i) = [IO_intValue(fileContent(l+1+i),chunkPos,1),i]
-      enddo
+      end do
       exit
-    endif
+    end if
-  enddo
+  end do
 
   call math_sort(map_unsorted)
   allocate(FEM2DAMASK(minval(map_unsorted(1,:)):maxval(map_unsorted(1,:))),source=-1)
   do i = 1, nNodes
     FEM2DAMASK(map_unsorted(1,i)) = map_unsorted(2,i)
-  enddo
+  end do
 
 end subroutine inputRead_mapNodes
 
@@ -560,10 +560,10 @@ subroutine inputRead_elemNodes(nodes, &
       do i=1,nNode
         m = discretization_Marc_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1))
         nodes(1:3,m) = [(mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1),j=1,3)]
-      enddo
+      end do
       exit
-    endif
+    end if
-  enddo
+  end do
 
 end subroutine inputRead_elemNodes
 
@@ -596,17 +596,17 @@ subroutine inputRead_elemType(elem, &
         else
           t_ = mapElemtype(IO_stringValue(fileContent(l+1+i+j),chunkPos,2))
           if (t /= t_) call IO_error(191,IO_stringValue(fileContent(l+1+i+j),chunkPos,2),label1='type',ID1=t)
-        endif
+        end if
         remainingChunks = elem%nNodes - (chunkPos(1) - 2)
         do while(remainingChunks > 0)
           j = j + 1
           chunkPos = IO_stringPos(fileContent(l+1+i+j))
           remainingChunks = remainingChunks - chunkPos(1)
-        enddo
+        end do
-      enddo
+      end do
       exit
-    endif
+    end if
-  enddo
+  end do
 
   contains
 
@@ -686,7 +686,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
           do k = 1,chunkPos(1)-2
             inputRead_connectivityElem(k,e) = &
               discretization_Marc_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2))
-          enddo
+          end do
           nNodesAlreadyRead = chunkPos(1) - 2
           do while(nNodesAlreadyRead < nNodes)                                                      ! read on if not all nodes in one line
             j = j + 1
@@ -694,14 +694,14 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
             do k = 1,chunkPos(1)
               inputRead_connectivityElem(nNodesAlreadyRead+k,e) = &
                 discretization_Marc_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k))
-            enddo
+            end do
             nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
-          enddo
+          end do
-        endif
+        end if
-      enddo
+      end do
       exit
-    endif
+    end if
-  enddo
+  end do
 
 end function inputRead_connectivityElem
 
@@ -749,12 +749,12 @@ subroutine inputRead_material(materialAt,&
           do i = 1,contInts(1)
             e = discretization_Marc_FEM2DAMASK_elem(contInts(1+i))
             materialAt(e) = ID + 1
-          enddo
+          end do
           if (initialcondTableStyle == 0) m = m + 1
-        enddo
+        end do
-      endif
+      end if
-    endif
+    end if
-  enddo
+  end do
 
   if(any(materialAt < 1)) call IO_error(180)
 
@@ -791,9 +791,9 @@ pure subroutine buildCells(connectivity,definition, &
     do c = 1, elem%NcellNodes
       realNode: if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == 1) then
         where(connectivity(:,:,e) == -c) connectivity(:,:,e) = connectivity_elem(c,e)
-      endif realNode
+      end if realNode
-    enddo
+    end do
-  enddo
+  end do
 
   nCellNode = maxval(connectivity_elem)
 
@@ -806,7 +806,7 @@ pure subroutine buildCells(connectivity,definition, &
     do c = 1, elem%NcellNodes
       if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == nParentNodes) &
         candidates_local = [candidates_local,c]
-    enddo
+    end do
     s = size(candidates_local)
 
     if (allocated(candidates_global)) deallocate(candidates_global)
@@ -822,8 +822,8 @@ pure subroutine buildCells(connectivity,definition, &
           if (elem%cellNodeParentNodeWeights(j,c) /= 0) then                                        ! real node 'j' partly defines cell node 'c'
             p = p + 1
             parentsAndWeights(p,1:2) = [connectivity_elem(j,e),elem%cellNodeParentNodeWeights(j,c)]
-          endif
+          end if
-        enddo
+        end do
         ! store (and order) real node IDs and their weights together with the element number and local ID
         do p = 1, nParentNodes
           m = maxloc(parentsAndWeights(:,1),1)
@@ -833,9 +833,9 @@ pure subroutine buildCells(connectivity,definition, &
           candidates_global(nParentNodes*2+1:nParentNodes*2+2,candidateID) = [e,c]
 
           parentsAndWeights(m,1) = -huge(parentsAndWeights(m,1))                                    ! out of the competition
-        enddo
+        end do
-      enddo
+      end do
-    enddo
+    end do
 
     ! sort according to real node IDs + weight (from left to right)
     call math_sort(candidates_global,sortDim=1)                                                     ! sort according to first column
@@ -847,13 +847,13 @@ pure subroutine buildCells(connectivity,definition, &
         do while (n+j<= size(candidates_local)*Nelem)
           if (candidates_global(p-1,n+j)/=candidates_global(p-1,n)) exit
           j = j + 1
-        enddo
+        end do
         e = n+j-1
         if (any(candidates_global(p,n:e)/=candidates_global(p,n))) &
           call math_sort(candidates_global(:,n:e),sortDim=p)
         n = e+1
-      enddo
+      end do
-    enddo
+    end do
 
     i = uniqueRows(candidates_global(1:2*nParentNodes,:))
     allocate(definition(nParentNodes-1)%parents(i,nParentNodes))
@@ -876,15 +876,15 @@ pure subroutine buildCells(connectivity,definition, &
         end where
 
         j = j+1
-      enddo
+      end do
       nCellNode = nCellNode + 1
       definition(nParentNodes-1)%parents(i,:) = parentsAndWeights(:,1)
       definition(nParentNodes-1)%weights(i,:) = parentsAndWeights(:,2)
       i = i + 1
       n = n+j
-    enddo
+    end do
 
-  enddo
+  end do
 
   contains
   !------------------------------------------------------------------------------------------------
@@ -906,10 +906,10 @@ pure subroutine buildCells(connectivity,definition, &
       do while (r+d<= size(A,2))
         if (any(A(:,r)/=A(:,r+d))) exit
         d = d+1
-      enddo
+      end do
       u = u+1
       r = r+d
-    enddo
+    end do
 
   end function uniqueRows
 
@@ -939,10 +939,10 @@ pure function buildCellNodes(node_elem)
         buildCellNodes(:,n) = buildCellNodes(:,n) &
                             + buildCellNodes(:,cellNodeDefinition(i)%parents(j,k)) &
                             * real(cellNodeDefinition(i)%weights(j,k),pReal)
-      enddo
+      end do
       buildCellNodes(:,n) = buildCellNodes(:,n)/real(sum(cellNodeDefinition(i)%weights(j,:)),pReal)
-    enddo
+    end do
-  enddo
+  end do
 
 end function buildCellNodes
 
@@ -970,9 +970,9 @@ pure function buildIPcoordinates(node_cell)
     do n = 1, size(connectivity_cell_reshaped,1)
       buildIPcoordinates(:,i) = buildIPcoordinates(:,i) &
                               + node_cell(:,connectivity_cell_reshaped(n,i))
-    enddo
+    end do
     buildIPcoordinates(:,i) = buildIPcoordinates(:,i)/real(size(connectivity_cell_reshaped,1),pReal)
-  enddo
+  end do
 
 end function buildIPcoordinates
 
@@ -1031,8 +1031,8 @@ pure function IPvolume(elem,node)
                         + dot_product((x7-x1),        math_cross((x5-x0),        (x7-x4)+(x3-x0)))
           IPvolume(i,e) = IPvolume(i,e)/12.0_pReal
       end select
-    enddo
+    end do
-  enddo
+  end do
 
 end function IPvolume
 
@@ -1075,11 +1075,11 @@ pure function IPareaNormal(elem,nElem,node)
               IPareaNormal(1:3,f,i,e) = IPareaNormal(1:3,f,i,e) &
                                       + math_cross(nodePos(1:3,mod(n+0,m)+1) - nodePos(1:3,n), &
                                                    nodePos(1:3,mod(n+1,m)+1) - nodePos(1:3,n)) * 0.5_pReal
-            enddo
+            end do
         end select
-      enddo
+      end do
-    enddo
+    end do
-  enddo
+  end do
 
 end function IPareaNormal
 
@@ -1109,10 +1109,10 @@ function IPneighborhood(elem)
         do n = 1, size(face_unordered)
           face(n,c) = minval(face_unordered)
           face_unordered(minloc(face_unordered)) = huge(face_unordered)
-        enddo
+        end do
         face(n:n+3,c) = [e,i,f]
-      enddo
+      end do
-  enddo; enddo
+  end do; end do
 
 !--------------------------------------------------------------------------------------------------
 ! sort face definitions
@@ -1125,17 +1125,17 @@ function IPneighborhood(elem)
       if(any(face(:c,s) /= face(:c,e))) then
         if(e-1/=s) call math_sort(face(:,s:e-1),sortDim=c)
         s = e
-      endif
+      end if
-    enddo
+    end do
-  enddo
+  end do
 
   IPneighborhood = 0
   do c=1, size(face,2) - 1
     if(all(face(:n-1,c) == face(:n-1,c+1))) then
       IPneighborhood(:,face(n+2,c+1),face(n+1,c+1),face(n+0,c+1)) = face(n:n+3,c+0)
       IPneighborhood(:,face(n+2,c+0),face(n+1,c+0),face(n+0,c+0)) = face(n:n+3,c+1)
-    endif
+    end if
-  enddo
+  end do
 
 end function IPneighborhood
 
@@ -1171,8 +1171,8 @@ function continuousIntValues(fileContent,maxN,lookupName,lookupMap,lookupMaxN)
         if (IO_stringValue(fileContent(l),chunkPos,1) == lookupName(i)) then                        ! found matching name
           continuousIntValues = lookupMap(:,i)                                                      ! return resp. entity list
           exit
-        endif
+        end if
-      enddo
+      end do
       exit
     elseif(containsRange(fileContent(l),chunkPos)) then
       first = IO_intValue(fileContent(l),chunkPos,1)
@@ -1180,20 +1180,20 @@ function continuousIntValues(fileContent,maxN,lookupName,lookupMap,lookupMaxN)
       do i = first, last, sign(1,last-first)
         continuousIntValues(1) = continuousIntValues(1) + 1
         continuousIntValues(1+continuousIntValues(1)) = i
-      enddo
+      end do
       exit
     else
       do i = 1,chunkPos(1)-1                                                                        ! interpret up to second to last value
         continuousIntValues(1) = continuousIntValues(1) + 1
         continuousIntValues(1+continuousIntValues(1)) = IO_intValue(fileContent(l),chunkPos,i)
-      enddo
+      end do
       if ( IO_lc(IO_stringValue(fileContent(l),chunkPos,chunkPos(1))) /= 'c' ) then                 ! line finished, read last value
         continuousIntValues(1) = continuousIntValues(1) + 1
         continuousIntValues(1+continuousIntValues(1)) = IO_intValue(fileContent(l),chunkPos,chunkPos(1))
         exit
-      endif
+      end if
-    endif
+    end if
-  enddo
+  end do
 
 end function continuousIntValues
 
@@ -1210,7 +1210,7 @@ logical function containsRange(str,chunkPos)
   containsRange = .False.
   if(chunkPos(1) == 3) then
     if(IO_lc(IO_stringValue(str,chunkPos,2)) == 'to') containsRange = .True.
-  endif
+  end if
 
 end function containsRange
 

src/Marc/materialpoint_Marc.f90

@@ -126,7 +126,7 @@ subroutine materialpoint_init
     print'(a32,1x,6(i8,1x))',   'materialpoint_dcsdE:           ', shape(materialpoint_dcsdE)
     print'(a32,1x,6(i8,1x),/)', 'materialpoint_dcsdE_knownGood: ', shape(materialpoint_dcsdE_knownGood)
     flush(IO_STDOUT)
-  endif
+  end if
 
 end subroutine materialpoint_init
 
@@ -171,7 +171,7 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
     if (terminallyIll) &
     print'(a,/)', '#           --- terminallyIll ---           #'
     print'(a,/)', '#############################################'; flush (6)
-  endif
+  end if
 
   if (iand(mode, materialpoint_BACKUPJACOBIAN) /= 0) &
     materialpoint_dcsde_knownGood = materialpoint_dcsde
@@ -220,15 +220,15 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
                      -  math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) &
                      +  0.5_pReal * (  Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
                                      + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
-        enddo; enddo; enddo; enddo; enddo; enddo
+        end do; end do; end do; end do; end do; end do
 
         forall(i=1:3, j=1:3,k=1:3,l=1:3) &
           H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k))
 
         materialpoint_dcsde(1:6,1:6,ip,elCP) = math_sym3333to66(J_inverse * H_sym,weighted=.false.)
 
-      endif terminalIllness
+      end if terminalIllness
-    endif validCalculation
+    end if validCalculation
 
     if (debugmaterialpoint%extensive &
         .and. ((debugmaterialpoint%element == elCP .and. debugmaterialpoint%ip == ip) .or. .not. debugmaterialpoint%selective)) then
@@ -237,9 +237,9 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
         print'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))', &
           '<< materialpoint >> Jacobian/GPa at elFE ip ', elFE, ip, transpose(materialpoint_dcsdE(1:6,1:6,ip,elCP))*1.0e-9_pReal
         flush(IO_STDOUT)
-    endif
+    end if
 
-  endif
+  end if
 
   if (all(abs(materialpoint_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) &
     call IO_warning(601,label1='element (CP)',ID1=elCP,label2='IP',ID2=ip)

src/YAML_parse.f90

@@ -764,7 +764,7 @@ end subroutine dct
 
 
 !--------------------------------------------------------------------------------------------------
-! @brief decide whether next block is list or dict
+! @brief Decide whether next block is list or dict.
 !--------------------------------------------------------------------------------------------------
 recursive subroutine decide(blck,flow,s_blck,s_flow,offset)
 
@@ -811,7 +811,7 @@ recursive subroutine decide(blck,flow,s_blck,s_flow,offset)
     end if
   end if
 
-end subroutine
+end subroutine decide
 
 
 !--------------------------------------------------------------------------------------------------

src/grid/DAMASK_grid.f90

@@ -148,7 +148,7 @@ program DAMASK_grid
     call results_openJobFile(parallel=.false.)
     call results_writeDataset_str(fileContent,'setup',fname,'load case definition (grid solver)')
     call results_closeJobFile
-  endif
+  end if
 
   call parallelization_bcast_str(fileContent)
   config_load => YAML_parse_str_asDict(fileContent)
@@ -198,11 +198,11 @@ program DAMASK_grid
   thermalActive: if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') then
     field = field + 1
     ID(field) = FIELD_THERMAL_ID
-  endif thermalActive
+  end if thermalActive
   damageActive: if (solver%get_asString('damage',defaultVal = 'n/a') == 'spectral') then
     field = field + 1
     ID(field) = FIELD_DAMAGE_ID
-  endif damageActive
+  end if damageActive
 
 
 !--------------------------------------------------------------------------------------------------
@@ -235,7 +235,7 @@ program DAMASK_grid
 #endif
       end select
       call loadCases(l)%rot%fromAxisAngle(step_mech%get_as1dFloat('R',defaultVal = real([0.0,0.0,1.0,0.0],pReal)),degrees=.true.)
-    enddo readMech
+    end do readMech
     if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing')
 
     step_discretization => load_step%get_dict('discretization')
@@ -264,9 +264,9 @@ program DAMASK_grid
           write(IO_STDOUT,'(2x,12a)',advance='no') '     x      '
         else
           write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j)
-        endif
+        end if
-        enddo; write(IO_STDOUT,'(/)',advance='no')
+        end do; write(IO_STDOUT,'(/)',advance='no')
-      enddo
+      end do
       if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831
       if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) &
         errorID = 838                                                                               ! no rotation is allowed by stress BC
@@ -280,10 +280,10 @@ program DAMASK_grid
             write(IO_STDOUT,'(2x,12a)',advance='no') '     x      '
           else
             write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal
-          endif
+          end if
-          enddo; write(IO_STDOUT,'(/)',advance='no')
+          end do; write(IO_STDOUT,'(/)',advance='no')
-        enddo
+        end do
-      endif
+      end if
       if (any(dNeq(loadCases(l)%rot%asMatrix(), math_I3))) &
         write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',&
                  transpose(loadCases(l)%rot%asMatrix())
@@ -298,7 +298,7 @@ program DAMASK_grid
         print'(2x,a)', 'r: 1 (constant step width)'
       else
         print'(2x,a,1x,f0.3)', 'r:', loadCases(l)%r
-      endif
+      end if
       print'(2x,a,1x,f0.3)',   't:', loadCases(l)%t
       print'(2x,a,1x,i0)',     'N:', loadCases(l)%N
       if (loadCases(l)%f_out < huge(0)) &
@@ -308,8 +308,8 @@ program DAMASK_grid
 
       if (errorID > 0) call IO_error(errorID,label1='line',ID1=l)
 
-    endif reportAndCheck
+    end if reportAndCheck
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! doing initialization depending on active solvers
@@ -337,14 +337,14 @@ program DAMASK_grid
     else writeHeader
       open(newunit=statUnit,file=trim(getSolverJobName())//&
                                   '.sta',form='FORMATTED', position='APPEND', status='OLD')
-    endif writeHeader
+    end if writeHeader
-  endif
+  end if
 
   writeUndeformed: if (CLI_restartInc < 1) then
     print'(/,1x,a)', '... writing initial configuration to file .................................'
     flush(IO_STDOUT)
     call materialpoint_results(0,0.0_pReal)
-  endif writeUndeformed
+  end if writeUndeformed
 
   loadCaseLooping: do l = 1, size(loadCases)
     t_0 = t                                                                                         ! load case start time
@@ -361,7 +361,7 @@ program DAMASK_grid
       else
         Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) &
                                  / (1.0_pReal-loadCases(l)%r**loadCases(l)%N)
-      endif
+      end if
       Delta_t = Delta_t * real(subStepFactor,pReal)**real(-cutBackLevel,pReal)                      ! depending on cut back level, decrease time step
 
       skipping: if (totalIncsCounter <= CLI_restartInc) then                                  ! not yet at restart inc?
@@ -402,7 +402,7 @@ program DAMASK_grid
               case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
               case(FIELD_DAMAGE_ID);  call grid_damage_spectral_forward(cutBack)
             end select
-          enddo
+          end do
           if (.not. cutBack) call materialpoint_forward
 
 !--------------------------------------------------------------------------------------------------
@@ -422,12 +422,12 @@ program DAMASK_grid
 
               if (.not. solres(field)%converged) exit                                               ! no solution found
 
-            enddo
+            end do
             stagIter = stagIter + 1
             stagIterate =            stagIter < stagItMax &
                          .and.       all(solres(:)%converged) &
                          .and. .not. all(solres(:)%stagConverged)                                   ! stationary with respect to staggered iteration
-          enddo
+          end do
 
 !--------------------------------------------------------------------------------------------------
 ! check solution for either advance or retry
@@ -442,7 +442,7 @@ program DAMASK_grid
               write(statUnit,*) totalIncsCounter, t, cutBackLevel, &
                                 solres(1)%converged, solres(1)%iterationsNeeded
               flush(statUnit)
-            endif
+            end if
           elseif (cutBackLevel < maxCutBack) then                                                   ! further cutbacking tolerated?
             cutBack = .true.
             stepFraction = (stepFraction - 1) * subStepFactor                                       ! adjust to new denominator
@@ -453,9 +453,9 @@ program DAMASK_grid
           else                                                                                      ! no more options to continue
             if (worldrank == 0) close(statUnit)
             call IO_error(950)
-          endif
+          end if
 
-        enddo subStepLooping
+        end do subStepLooping
 
         cutBackLevel = max(0, cutBackLevel - 1)                                                     ! try half number of subincs next inc
 
@@ -463,7 +463,7 @@ program DAMASK_grid
           print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged'
         else
           print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged'
-        endif; flush(IO_STDOUT)
+        end if; flush(IO_STDOUT)
 
         call MPI_Allreduce(signals_SIGUSR1,signal,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
         if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@@ -471,7 +471,7 @@ program DAMASK_grid
           print'(/,1x,a)', '... writing results to file ...............................................'
           flush(IO_STDOUT)
           call materialpoint_results(totalIncsCounter,t)
-        endif
+        end if
         if (signal) call signals_setSIGUSR1(.false.)
         call MPI_Allreduce(signals_SIGUSR2,signal,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
         if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
@@ -485,16 +485,16 @@ program DAMASK_grid
             end select
           end do
           call materialpoint_restartWrite
-        endif
+        end if
         if (signal) call signals_setSIGUSR2(.false.)
         call MPI_Allreduce(signals_SIGINT,signal,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
         if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
         if (signal) exit loadCaseLooping
-      endif skipping
+      end if skipping
 
-    enddo incLooping
+    end do incLooping
 
-  enddo loadCaseLooping
+  end do loadCaseLooping
 
 
 !--------------------------------------------------------------------------------------------------
@@ -523,9 +523,9 @@ subroutine getMaskedTensor(values,mask,tensor)
     do j = 1,3
       mask(i,j) = row%get_asString(j) == 'x'
       if (.not. mask(i,j)) values(i,j) = row%get_asFloat(j)
-    enddo
+    end do
-  enddo
+  end do
 
-end subroutine
+end subroutine getMaskedTensor
 
 end program DAMASK_grid

src/grid/VTI.f90

@@ -222,7 +222,7 @@ subroutine cellsSizeOrigin(c,s,o,header)
   temp = getXMLValue(header,'Origin')
   o = [(IO_floatValue(temp,IO_stringPos(temp),i),i=1,3)]
 
-end subroutine
+end subroutine cellsSizeOrigin
 
 
 !--------------------------------------------------------------------------------------------------
@@ -421,7 +421,7 @@ pure function getXMLValue(line,key)
     end if
   end if
 
-end function
+end function getXMLValue
 
 
 !--------------------------------------------------------------------------------------------------

src/grid/base64.f90

@@ -82,7 +82,7 @@ function base64_to_bytes(base64_str,s,e) result(bytes)
   else
     s_str = 1_pI64
     s_bytes = 1_pI64
-  endif
+  end if
 
   if(present(e)) then
     if(e>base64_nByte(len(base64_str,kind=pI64))) call IO_error(114, ext_msg='e out of range')
@@ -93,7 +93,7 @@ function base64_to_bytes(base64_str,s,e) result(bytes)
     e_bytes = base64_nByte(len(base64_str,kind=pI64)) - base64_nByte(s_str)
     if(base64_str(e_str-0_pI64:e_str-0_pI64) == '=') e_bytes = e_bytes - 1_pI64
     if(base64_str(e_str-1_pI64:e_str-1_pI64) == '=') e_bytes = e_bytes - 1_pI64
-  endif
+  end if
 
   bytes = decodeBase64(base64_str(s_str:e_str))
   bytes = bytes(s_bytes:e_bytes)
@@ -122,8 +122,8 @@ pure function decodeBase64(base64_str) result(bytes)
         charPos(p) = int(index(base64_encoding,base64_str(c+p:c+p))-1,C_SIGNED_CHAR)
       else
         charPos(p) = 0_C_SIGNED_CHAR
-      endif
+      end if
-    enddo
+    end do
 
     call mvbits(charPos(0),0,6,bytes(b+0),2)
     call mvbits(charPos(1),4,2,bytes(b+0),0)
@@ -133,7 +133,7 @@ pure function decodeBase64(base64_str) result(bytes)
     call mvbits(charPos(3),0,6,bytes(b+2),0)
     b = b+3_pI64
     c = c+4_pI64
-  enddo
+  end do
 
 end function decodeBase64
 

src/grid/discretization_grid.f90

@@ -334,7 +334,7 @@ function discretization_grid_getInitialCondition(label) result(ic)
     ic_global = VTI_readDataset_real(IO_read(CLI_geomFile),label)
   else
     allocate(ic_global(0))                                                                          ! needed for IntelMPI
-  endif
+  end if
 
   call MPI_Gather(product(cells(1:2))*cells3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,&
                   1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)

src/grid/grid_mech_FEM.f90

@@ -274,7 +274,7 @@ subroutine grid_mechanical_FEM_init
   elseif (CLI_restartInc == 0) then restartRead
     F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3)                      ! initialize to identity
     F         = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3)
-  endif restartRead
+  end if restartRead
 
   homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3])                          ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(F)
@@ -298,7 +298,7 @@ subroutine grid_mechanical_FEM_init
     call HDF5_closeGroup(groupHandle)
     call HDF5_closeFile(fileHandle)
 
-  endif restartRead2
+  end if restartRead2
 
 end subroutine grid_mechanical_FEM_init
 
@@ -387,7 +387,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
     elseif (deformation_BC%myType=='F') then                                                        ! aim at end of load case is prescribed
       F_aimDot = F_aimDot &
                + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
-    endif
+    end if
 
     if (guess) then
       call VecWAXPY(solution_rate,-1.0_pReal,solution_lastInc,solution_current,err_PETSc)
@@ -397,14 +397,14 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
     else
       call VecSet(solution_rate,0.0_pReal,err_PETSc)
       CHKERRQ(err_PETSc)
-    endif
+    end if
     call VecCopy(solution_current,solution_lastInc,err_PETSc)
     CHKERRQ(err_PETSc)
 
     F_lastInc = F
 
     homogenization_F0 = reshape(F, [3,3,product(cells(1:2))*cells3])
-  endif
+  end if
 
 !--------------------------------------------------------------------------------------------------
 ! update average and local deformation gradients
@@ -477,7 +477,7 @@ subroutine grid_mechanical_FEM_restartWrite
     call HDF5_write(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.)
     call HDF5_closeGroup(groupHandle)
     call HDF5_closeFile(fileHandle)
-  endif
+  end if
 
   call DMDAVecRestoreArrayF90(mechanical_grid,solution_current,u,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -517,7 +517,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,e
     reason = -1
   else
     reason = 0
-  endif
+  end if
 
   print'(/,1x,a)', '... reporting .............................................................'
   print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', &
@@ -567,7 +567,7 @@ subroutine formResidual(da_local,x_local, &
     print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
       'deformation gradient aim       =', transpose(F_aim)
     flush(IO_STDOUT)
-  endif newIteration
+  end if newIteration
 
 !--------------------------------------------------------------------------------------------------
 ! get deformation gradient
@@ -578,9 +578,9 @@ subroutine formResidual(da_local,x_local, &
     do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
       ctr = ctr + 1
       x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
-    enddo; enddo; enddo
+    end do; end do; end do
     F(1:3,1:3,i,j,k-cells3Offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
-  enddo; enddo; enddo
+  end do; end do; end do
   call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc)
   CHKERRQ(err_PETSc)
 
@@ -611,7 +611,7 @@ subroutine formResidual(da_local,x_local, &
     do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
       ctr = ctr + 1
       x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
-    enddo; enddo; enddo
+    end do; end do; end do
     ele = ele + 1
     f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-cells3Offset)))*detJ + &
              matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + &
@@ -621,8 +621,8 @@ subroutine formResidual(da_local,x_local, &
     do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
       ctr = ctr + 1
       r(0:2,i+ii,j+jj,k+kk) = r(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3)
-    enddo; enddo; enddo
+    end do; end do; end do
-  enddo; enddo; enddo
+  end do; end do; end do
   call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc)
   CHKERRQ(err_PETSc)
   call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
@@ -696,7 +696,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
       col(MatStencil_j,ctr+16) = j+jj
       col(MatStencil_k,ctr+16) = k+kk
       col(MatStencil_c,ctr+16) = 2
-    enddo; enddo; enddo
+    end do; end do; end do
     row = col
     ce = ce + 1
     K_ele = 0.0_pReal
@@ -715,7 +715,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
                                           shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ
     call MatSetValuesStencil(Jac,24_pPETScInt,row,24_pPetscInt,col,K_ele,ADD_VALUES,err_PETSc)
     CHKERRQ(err_PETSc)
-  enddo; enddo; enddo
+  end do; end do; end do
   call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
   CHKERRQ(err_PETSc)
   call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
@@ -739,7 +739,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
   do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1)
     ce = ce + 1
     x_scal(0:2,i-1,j-1,k-1) = discretization_IPcoords(1:3,ce)
-  enddo; enddo; enddo
+  end do; end do; end do
   call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,err_PETSc)
   CHKERRQ(err_PETSc)                                                                                ! initialize to undeformed coordinates (ToDo: use ip coordinates)
   call MatNullSpaceCreateRigidBody(coordinates,matnull,err_PETSc)

src/materialpoint.f90

@@ -95,7 +95,7 @@ subroutine materialpoint_init
     call phase_restartRead(fileHandle)
 
     call HDF5_closeFile(fileHandle)
-  endif
+  end if
 
 end subroutine materialpoint_init
 

src/mesh/DAMASK_mesh.f90

@@ -311,7 +311,7 @@ program DAMASK_mesh
           write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
                             solres%converged, solres%iterationsNeeded                               ! write statistics about accepted solution
           flush(statUnit)
-        endif
+        end if
       end do subStepLooping
 
       cutBackLevel = max(0, cutBackLevel - 1)                                                       ! try half number of subincs next inc

src/mesh/FEM_quadrature.f90

@@ -365,16 +365,16 @@ subroutine selfTest
     do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1)
       if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) &
         error stop 'quadrature weights'
-    enddo
+    end do
-  enddo
+  end do
 
   do d = lbound(FEM_quadrature_points,1), ubound(FEM_quadrature_points,1)
     do o = lbound(FEM_quadrature_points(d,:),1), ubound(FEM_quadrature_points(d,:),1)
       n = size(FEM_quadrature_points(d,o)%p,1)/d
       if (any(dNeq(sum(reshape(FEM_quadrature_points(d,o)%p,[d,n]),2),-real(n,pReal)/w(d),1.e-14_pReal))) &
         error stop 'quadrature points'
-    enddo
+    end do
-  enddo
+  end do
 
 end subroutine selfTest
 

src/mesh/discretization_mesh.f90

@@ -140,7 +140,7 @@ subroutine discretization_mesh_init(restart)
     call DMClone(globalMesh,geomMesh,err_PETSc)
   else
     call DMPlexDistribute(globalMesh,0_pPETSCINT,sf,geomMesh,err_PETSc)
-  endif
+  end if
   CHKERRQ(err_PETSc)
 
   allocate(mesh_boundaries(mesh_Nboundaries), source = 0_pPETSCINT)
@@ -154,7 +154,7 @@ subroutine discretization_mesh_init(restart)
     mesh_boundaries(1:nFaceSets) = pFaceSets
     CHKERRQ(err_PETSc)
     call ISRestoreIndicesF90(faceSetIS,pFaceSets,err_PETSc)
-  endif
+  end if
   call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
 
@@ -182,7 +182,7 @@ subroutine discretization_mesh_init(restart)
   do j = 1, mesh_NcpElems
     call DMGetLabelValue(geomMesh,'Cell Sets',j-1,materialAt(j),err_PETSc)
     CHKERRQ(err_PETSc)
-  enddo
+  end do
   materialAt = materialAt + 1_pPETSCINT
 
   if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element')
@@ -222,7 +222,7 @@ subroutine mesh_FEM_build_ipVolumes(dimPlex)
     call  DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,err_PETSc)
     CHKERRQ(err_PETSc)
     mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
-  enddo
+  end do
 
 end subroutine mesh_FEM_build_ipVolumes
 
@@ -258,11 +258,11 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
         do dirJ = 1_pPETSCINT, dimPlex
           mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
                                                 pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0_pReal)
-        enddo
+        end do
-      enddo
+      end do
       qOffset = qOffset + dimPlex
-    enddo
+    end do
-  enddo
+  end do
 
 end subroutine mesh_FEM_build_ipCoordinates
 

src/mesh/mesh_mech_FEM.f90

@@ -199,11 +199,11 @@ subroutine FEM_mechanical_init(fieldBC)
     CHKERRQ(err_PETSc)
     call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),err_PETSc)
     CHKERRQ(err_PETSc)
-  enddo
+  end do
   numBC = 0
   do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
     if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
-  enddo; enddo
+  end do; end do
   allocate(pbcField(numBC), source=0_pPETSCINT)
   allocate(pbcComps(numBC))
   allocate(pbcPoints(numBC))
@@ -229,9 +229,9 @@ subroutine FEM_mechanical_init(fieldBC)
       else
         call ISCreateGeneral(PETSC_COMM_WORLD,0_pPETSCINT,[0_pPETSCINT],PETSC_COPY_VALUES,pbcPoints(numBC),err_PETSc)
         CHKERRQ(err_PETSc)
-      endif
+      end if
-    endif
+    end if
-  enddo; enddo
+  end do; end do
   call DMPlexCreateSection(mechanical_mesh,nolabel,pNumComp,pNumDof, &
                            numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -240,7 +240,7 @@ subroutine FEM_mechanical_init(fieldBC)
   do faceSet = 1, numBC
     call ISDestroy(pbcPoints(faceSet),err_PETSc)
     CHKERRQ(err_PETSc)
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! initialize solver specific parts of PETSc
@@ -299,11 +299,11 @@ subroutine FEM_mechanical_init(fieldBC)
       call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,err_PETSc)
       CHKERRQ(err_PETSc)
       x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal)
-    enddo
+    end do
     px_scal => x_scal
     call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,err_PETSc)
     CHKERRQ(err_PETSc)
-  enddo
+  end do
   call utilities_constitutiveResponse(0.0_pReal,devNull,.true.)
 
 end subroutine FEM_mechanical_init
@@ -348,7 +348,7 @@ type(tSolutionState) function FEM_mechanical_solution( &
     FEM_mechanical_solution%converged = .true.
     call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,err_PETSc)
     CHKERRQ(err_PETSc)
-  endif
+  end if
 
   print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
@@ -409,9 +409,9 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
                                        0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
         call ISDestroy(bcPoints,err_PETSc)
         CHKERRQ(err_PETSc)
-      endif
+      end if
-    endif
+    end if
-  enddo; enddo
+  end do; end do
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate field derivatives
@@ -433,10 +433,10 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
           i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
           BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
             matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
-        enddo
+        end do
-      enddo
+      end do
       homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
-    enddo
+    end do
     if (num%BBarStabilisation) then
       detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature,pReal))
       do qPt = 0, nQuadrature-1
@@ -444,11 +444,11 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
         homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) &
                                                 * (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0_pReal/real(dimPlex,pReal))
 
-      enddo
+      end do
-    endif
+    end if
     call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc)
     CHKERRQ(err_PETSc)
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate constitutive response
@@ -475,20 +475,20 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
           i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
           BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
             matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
-        enddo
+        end do
-      enddo
+      end do
       f_scal = f_scal &
              +  matmul(transpose(BMat), &
                        reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), &
                                shape=[dimPlex*dimPlex]))*qWeights(qPt+1_pPETSCINT)
-    enddo
+    end do
     f_scal = f_scal*abs(detJ)
     pf_scal => f_scal
     call DMPlexVecSetClosure(dm_local,section,f_local,cell,pf_scal,ADD_VALUES,err_PETSc)
     CHKERRQ(err_PETSc)
     call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc)
     CHKERRQ(err_PETSc)
-  enddo
+  end do
   call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
   CHKERRQ(err_PETSc)
 
@@ -559,9 +559,9 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
                                        0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
         call ISDestroy(bcPoints,err_PETSc)
         CHKERRQ(err_PETSc)
-      endif
+      end if
-    endif
+    end if
-  enddo; enddo
+  end do; end do
   call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc)
   CHKERRQ(err_PETSc)
   do cell = cellStart, cellEnd-1                                                                    !< loop over all elements
@@ -583,8 +583,8 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
           i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp
           BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = &
             matmul(reshape(pInvcellJ,[dimPlex,dimPlex]),basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex))
-        enddo
+        end do
-      enddo
+      end do
       MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), &
                                     shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
                             shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT)
@@ -602,8 +602,8 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
         BMatAvg = BMatAvg + BMat
       else
         K_eA = K_eA + matmul(transpose(BMat),MatA)
-      endif
+      end if
-    enddo
+    end do
     if (num%BBarStabilisation) then
       FInv = math_inv33(FAvg)
       K_e = K_eA*math_det33(FAvg/real(nQuadrature,pReal))**(1.0_pReal/real(dimPlex,pReal)) + &
@@ -612,7 +612,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
              K_eB)/real(dimPlex,pReal)
     else
       K_e = K_eA
-    endif
+    end if
     K_e = (K_e + eps*math_eye(int(cellDof))) * abs(detJ)
 #ifndef __INTEL_COMPILER
     pK_e(1:cellDOF**2) => K_e
@@ -624,7 +624,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
     CHKERRQ(err_PETSc)
     call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc)
     CHKERRQ(err_PETSc)
-  enddo
+  end do
   call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
   CHKERRQ(err_PETSc)
   call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
@@ -704,9 +704,9 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
                                          0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old)
           call ISDestroy(bcPoints,err_PETSc)
           CHKERRQ(err_PETSc)
-        endif
+        end if
-      endif
+      end if
-    enddo; enddo
+    end do; end do
     call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
     CHKERRQ(err_PETSc)
 
@@ -716,7 +716,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
     CHKERRQ(err_PETSc)
     call VecScale(solution_rate,timeinc_old**(-1),err_PETSc)
     CHKERRQ(err_PETSc)
-  endif
+  end if
   call VecCopy(solution_rate,solution,err_PETSc)
   CHKERRQ(err_PETSc)
   call VecScale(solution,timeinc,err_PETSc)
@@ -800,7 +800,7 @@ subroutine FEM_mechanical_updateCoords()
     call DMPlexGetPointLocal(dm_local, p, s, e, err_PETSc)
     CHKERRQ(err_PETSc)
     nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e)
-  enddo
+  end do
 
   call discretization_setNodeCoords(nodeCoords)
   call VecRestoreArrayF90(x_local,nodeCoords_linear,err_PETSc)
@@ -827,9 +827,9 @@ subroutine FEM_mechanical_updateCoords()
                                      x_scal(nOffset+1:nOffset+dimPlex))
           q = q+dimPlex
           nOffset = nOffset+dimPlex
-        enddo
+        end do
-      enddo
+      end do
-    enddo
+    end do
     call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,err_PETSc)
     CHKERRQ(err_PETSc)
   end do

src/phase_mechanical_plastic_nonlocal.f90

@@ -1783,6 +1783,6 @@ subroutine storeGeometry(ph)
     end do
   end do
 
-end subroutine
+end subroutine storeGeometry
 
 end submodule nonlocal