diff --git a/PRIVATE b/PRIVATE
index 226d7b627..b14f78e96 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 226d7b627e44247b800bce0d9eb7bef1aac6f537
+Subproject commit b14f78e96a8e2986aaf6845b98ea77fec92bc997
diff --git a/examples/config/phase/mechanical/plastic/dislotwin_alpha-Brass-shearbanding.yaml b/examples/config/phase/mechanical/plastic/dislotwin_alpha-Brass-shearbanding.yaml
new file mode 100644
index 000000000..414be0b85
--- /dev/null
+++ b/examples/config/phase/mechanical/plastic/dislotwin_alpha-Brass-shearbanding.yaml
@@ -0,0 +1,15 @@
+type: dislotwin
+
+references:
+  - N. Jia et al.,
+    Acta Materialia 60(3):1099-1115, 2012,
+    https://doi.org/10.1016/j.actamat.2011.10.047
+  - N. Jia et al.,
+    Acta Materialia 60:3415-3434, 2012,
+    https://doi.org/10.1016/j.actamat.2012.03.005
+
+gamma_0_sb: 0.0001
+tau_sb: 180.0e6 # tau_hat_sb
+Q_sb: 4.0e-19 # Q_0
+p_sb: 1.15
+q_sb: 1.0
diff --git a/src/CLI.f90 b/src/CLI.f90
index ff3d7349f..cc73365d6 100644
--- a/src/CLI.f90
+++ b/src/CLI.f90
@@ -156,15 +156,15 @@ subroutine CLI_init
         if (CLI_restartInc < 0 .or. stat /=0) then
           print'(/,a)', ' ERROR: Could not parse restart increment: '//trim(arg)
           call quit(1)
-        endif
+        end if
     end select
     if (err /= 0) call quit(1)
-  enddo
+  end do
 
   if (len_trim(loadcaseArg) == 0 .or. len_trim(geometryArg) == 0) then
     print'(/,a)', ' ERROR: Please specify geometry AND load case (-h for help)'
     call quit(1)
-  endif
+  end if
 
   if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg))
   CLI_geomFile = getGeometryFile(geometryArg)
@@ -205,14 +205,14 @@ subroutine setWorkingDirectory(workingDirectoryArg)
   else absolutePath
     workingDirectory = getCWD()
     workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg
-  endif absolutePath
+  end if absolutePath
 
   workingDirectory = trim(rectifyPath(workingDirectory))
   error = setCWD(trim(workingDirectory))
   if(error) then
     print*, 'ERROR: Invalid Working directory: '//trim(workingDirectory)
     call quit(1)
-  endif
+  end if
 
 end subroutine setWorkingDirectory
 
@@ -256,7 +256,7 @@ function getGeometryFile(geometryParameter)
   if (.not. file_exists) then
     print*, 'ERROR: Geometry file does not exists: '//trim(getGeometryFile)
     call quit(1)
-  endif
+  end if
 
 end function getGeometryFile
 
@@ -279,7 +279,7 @@ function getLoadCaseFile(loadCaseParameter)
   if (.not. file_exists) then
     print*, 'ERROR: Load case file does not exists: '//trim(getLoadCaseFile)
     call quit(1)
-  endif
+  end if
 
 end function getLoadCaseFile
 
@@ -300,14 +300,14 @@ function rectifyPath(path)
   l = len_trim(rectifyPath)
   do i = l,3,-1
     if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//'  '
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! remove // from path
   l = len_trim(rectifyPath)
   do i = l,2,-1
     if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' '
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! remove ../ and corresponding directory from rectifyPath
@@ -321,9 +321,9 @@ function rectifyPath(path)
        k = len_trim(rectifyPath)
        rectifyPath(j+1:k-1) = rectifyPath(j+2:k)
        rectifyPath(k:k) = ' '
-     endif
+     end if
      i = j+index(rectifyPath(j+1:l),'../')
-  enddo
+  end do
   if(len_trim(rectifyPath) == 0) rectifyPath = '/'
 
   rectifyPath = trim(rectifyPath)
@@ -349,10 +349,10 @@ function makeRelativePath(a,b)
   do i = 1, min(len_trim(a_cleaned),len_trim(rectifyPath(b_cleaned)))
     if (a_cleaned(i:i) /= b_cleaned(i:i)) exit
     if (a_cleaned(i:i) == '/') posLastCommonSlash = i
-  enddo
+  end do
   do i = posLastCommonSlash+1,len_trim(a_cleaned)
     if (a_cleaned(i:i) == '/') remainingSlashes = remainingSlashes + 1
-  enddo
+  end do
 
   makeRelativePath = repeat('..'//'/',remainingSlashes)//b_cleaned(posLastCommonSlash+1:len_trim(b_cleaned))
 
diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90
index fcd2c189a..c905164c2 100644
--- a/src/HDF5_utilities.f90
+++ b/src/HDF5_utilities.f90
@@ -509,7 +509,7 @@ subroutine HDF5_addAttribute_str_array(loc_id,attrLabel,attrValue,path)
   do i=1,size(attrValue)
     attrValue_(i) = attrValue(i)//C_NULL_CHAR
     ptr(i) = c_loc(attrValue_(i))
-  enddo
+  end do
 
   call H5Screate_simple_f(1,shape(attrValue_,kind=HSIZE_T),space_id,hdferr,shape(attrValue_,kind=HSIZE_T))
   if(hdferr < 0) error stop 'HDF5 error'
diff --git a/src/IO.f90 b/src/IO.f90
index a6ecb1491..9b9b9d563 100644
--- a/src/IO.f90
+++ b/src/IO.f90
@@ -92,15 +92,15 @@ function IO_readlines(fileName) result(fileContent)
       if (.not. warned) then
         call IO_warning(207,trim(fileName),label1='line',ID1=l)
         warned = .true.
-      endif
+      end if
     else
       line = rawData(startPos:endpos)
-    endif
+    end if
     startPos = endPos + 2                                                                           ! jump to next line start
 
     fileContent(l) = trim(line)//''
     l = l + 1
-  enddo
+  end do
 
 end function IO_readlines
 
@@ -129,7 +129,7 @@ function IO_read(fileName) result(fileContent)
   if (fileLength==0) then
     close(fileUnit)
     return
-  endif
+  end if
 
   read(fileUnit,iostat=myStat) fileContent
   if (myStat /= 0) call IO_error(102,trim(fileName))
@@ -183,8 +183,8 @@ pure function IO_stringPos(string)
     endOfString: if (right < left) then
       IO_stringPos(IO_stringPos(1)*2+1) = len_trim(string)
       exit
-    endif endOfString
-  enddo
+    end if endOfString
+  end do
 
 end function IO_stringPos
 
@@ -204,7 +204,7 @@ function IO_stringValue(string,chunkPos,myChunk)
     call IO_error(110,'IO_stringValue: "'//trim(string)//'"',label1='chunk',ID1=myChunk)
   else validChunk
     IO_stringValue = string(chunkPos(myChunk*2):chunkPos(myChunk*2+1))
-  endif validChunk
+  end if validChunk
 
 end function IO_stringValue
 
@@ -257,8 +257,8 @@ pure function IO_lc(string)
       IO_lc(i:i) = LOWER(n:n)
     else
       IO_lc(i:i) = string(i:i)
-    endif
-  enddo
+    end if
+  end do
 
 end function IO_lc
 
@@ -280,7 +280,7 @@ function IO_rmComment(line)
     IO_rmComment = trim(line)
   else
     IO_rmComment = trim(line(:split-1))
-  endif
+  end if
 
 end function IO_rmComment
 
@@ -302,7 +302,7 @@ integer function IO_stringAsInt(string)
   else valid
     IO_stringAsInt = 0
     call IO_error(111,string)
-  endif valid
+  end if valid
 
 end function IO_stringAsInt
 
@@ -324,7 +324,7 @@ real(pReal) function IO_stringAsFloat(string)
   else valid
     IO_stringAsFloat = 0.0_pReal
     call IO_error(112,string)
-  endif valid
+  end if valid
 
 end function IO_stringAsFloat
 
@@ -344,7 +344,7 @@ logical function IO_stringAsBool(string)
   else
     IO_stringAsBool = .false.
     call IO_error(113,string)
-  endif
+  end if
 
 end function IO_stringAsBool
 
@@ -598,7 +598,7 @@ pure function CRLF2LF(string)
     CRLF2LF(c-n:c-n) = string(c:c)
     if (c == len_trim(string)) exit
     if (string(c:c+1) == CR//LF) n = n + 1
-  enddo
+  end do
 
   CRLF2LF = CRLF2LF(:c-n)
 
@@ -639,17 +639,17 @@ subroutine panel(paneltype,ID,msg,ext_msg,label1,ID1,label2,ID2)
     write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
                                                        max(1,panelwidth+3-len_trim(ext_msg)-4),'x,a)'
     write(IO_STDERR,formatString)          '│ ',trim(ext_msg),                                      '│'
-  endif
+  end if
   if (present(label1)) then
     write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label1)),',i9,',&
                                                        max(1,panelwidth+3-len_trim(label1)-9-7),'x,a)'
     write(IO_STDERR,formatString)          '│ at ',trim(label1),ID1,                                '│'
-  endif
+  end if
   if (present(label2)) then
     write(formatString,'(a,i3.3,a,i3.3,a)') '(1x,a7,a',max(1,len_trim(label2)),',i9,',&
                                                        max(1,panelwidth+3-len_trim(label2)-9-7),'x,a)'
     write(IO_STDERR,formatString)          '│ at ',trim(label2),ID2,                                '│'
-  endif
+  end if
   write(formatString,'(a,i2.2,a)') '(a,',max(1,panelwidth),'x,a)'
   write(IO_STDERR,formatString)          ' │',                                                     '│'
   write(IO_STDERR,'(a)')                  ' └'//DIVIDER//'┘'
diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90
index a887476a9..67abfdd90 100644
--- a/src/YAML_parse.f90
+++ b/src/YAML_parse.f90
@@ -103,7 +103,7 @@ recursive function parse_flow(YAML_flow) result(node)
         class is (tDict)
           call node%set(key,myVal)
       end select
-    enddo
+    end do
   elseif (flow_string(1:1) == '[') then                                                             ! start of a list
     e = 1
     allocate(tList::node)
@@ -116,7 +116,7 @@ recursive function parse_flow(YAML_flow) result(node)
         class is (tList)
           call node%append(myVal)
       end select
-    enddo
+    end do
   else                                                                                              ! scalar value
     allocate(tScalar::node)
       select type (node)
@@ -156,7 +156,7 @@ integer function find_end(str,e_char)
     N_sq = N_sq - merge(1,0,str(i:i) == ']')
     N_cu = N_cu - merge(1,0,str(i:i) == '}')
     i = i + 1
-  enddo
+  end do
   find_end = i
 
 end function find_end
@@ -333,7 +333,7 @@ subroutine skip_empty_lines(blck,s_blck)
   do while(empty .and. len_trim(blck(s_blck:)) /= 0)
     empty = len_trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == 0
     if(empty) s_blck = s_blck + index(blck(s_blck:),IO_EOL)
-  enddo
+  end do
 
 end subroutine skip_empty_lines
 
@@ -387,7 +387,7 @@ logical function flow_is_closed(str,e_char)
     N_cu = N_cu + merge(1,0,line(i:i) == '{')
     N_sq = N_sq - merge(1,0,line(i:i) == ']')
     N_cu = N_cu - merge(1,0,line(i:i) == '}')
-  enddo
+  end do
 
 end function flow_is_closed
 
@@ -410,7 +410,7 @@ subroutine remove_line_break(blck,s_blck,e_char,flow_line)
     flow_line = flow_line//IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))//' '
     line_end  = flow_is_closed(flow_line,e_char)
     s_blck    = s_blck + index(blck(s_blck:),IO_EOL)
-  enddo
+  end do
 
 end subroutine remove_line_break
 
@@ -439,7 +439,7 @@ subroutine list_item_inline(blck,s_blck,inline,offset)
     inline = inline//' '//trim(adjustl(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))))
     s_blck = s_blck + index(blck(s_blck:),IO_EOL)
     indent_next = indentDepth(blck(s_blck:))
-  enddo
+  end do
 
   if(scan(inline,",") > 0) inline = '"'//inline//'"'
 
@@ -481,7 +481,7 @@ recursive subroutine line_isFlow(flow,s_flow,line)
       flow(s_flow:s_flow+1) = ', '
       s_flow = s_flow +2
       s = s + find_end(line(s+1:),']')
-    enddo
+    end do
     s_flow = s_flow - 1
     if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
     flow(s_flow:s_flow) = ']'
@@ -498,7 +498,7 @@ recursive subroutine line_isFlow(flow,s_flow,line)
       flow(s_flow:s_flow+1) = ', '
       s_flow = s_flow +2
       s = s + find_end(line(s+1:),'}')
-    enddo
+    end do
     s_flow = s_flow -1
     if(flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow -1
     flow(s_flow:s_flow) = '}'
@@ -646,7 +646,7 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
       s_flow = s_flow + 2
     end if
 
-  enddo
+  end do
 
   s_flow = s_flow - 1
   if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
@@ -733,7 +733,7 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset)
     flow(s_flow:s_flow) = ' '
     s_flow = s_flow + 1
     offset = 0
-  enddo
+  end do
 
   s_flow = s_flow - 1
   if (flow(s_flow-1:s_flow-1) == ',') s_flow = s_flow - 1
diff --git a/src/YAML_types.f90 b/src/YAML_types.f90
index 857b24bd0..e685209e6 100644
--- a/src/YAML_types.f90
+++ b/src/YAML_types.f90
@@ -411,7 +411,7 @@ function tNode_get_byIndex(self,i) result(node)
 
   do j = 2,i
     item => item%next
-  enddo
+  end do
   node => item%node
 
 end function tNode_get_byIndex
@@ -681,7 +681,7 @@ function tNode_contains(self,k)  result(exists)
           exists = .true.
           return
         end if
-      enddo
+      end do
     class is(tList)
       list => self%asList()
       do j=1, list%length
@@ -689,7 +689,7 @@ function tNode_contains(self,k)  result(exists)
           exists = .true.
           return
         end if
-      enddo
+      end do
     class default
       call IO_error(706,ext_msg='Expected list or dict')
   end select
@@ -731,7 +731,7 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
     end if
     item => item%next
     j = j + 1
-  enddo
+  end do
 
   if (.not. found) then
     call IO_error(143,ext_msg=k)
@@ -1333,7 +1333,7 @@ function tList_as1dString(self)
     scalar => item%node%asScalar()
     tList_as1dString(i) = scalar%asString()
     item => item%next
-  enddo
+  end do
 
 end function tList_as1dString
 
@@ -1384,7 +1384,7 @@ subroutine tDict_set(self,key,node)
     searchExisting: do while (associated(item%next))
       if (item%key == key) exit
       item => item%next
-    enddo searchExisting
+    end do searchExisting
     if (item%key /= key) then
       allocate(item%next)
       item => item%next
diff --git a/src/discretization.f90 b/src/discretization.f90
index 5f4d1c5a6..95fd39fd5 100644
--- a/src/discretization.f90
+++ b/src/discretization.f90
@@ -68,7 +68,7 @@ subroutine discretization_init(materialAt,&
     discretization_sharedNodesBegin = sharedNodesBegin
   else
     discretization_sharedNodesBegin = size(discretization_NodeCoords0,2)
-  endif
+  end if
 
 end subroutine discretization_init
 
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 203c7bb1c..290f3de9b 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -365,7 +365,7 @@ subroutine homogenization_results
       call thermal_results(ho,group)
     end if
 
- enddo
+ end do
 
 end subroutine homogenization_results
 
@@ -383,7 +383,7 @@ subroutine homogenization_forward
     homogState (ho)%state0 = homogState (ho)%state
     if(damageState_h(ho)%sizeState > 0) &
       damageState_h(ho)%state0 = damageState_h(ho)%state
-  enddo
+  end do
 
 end subroutine homogenization_forward
 
@@ -408,7 +408,7 @@ subroutine homogenization_restartWrite(fileHandle)
 
     call HDF5_closeGroup(groupHandle(2))
 
-  enddo
+  end do
 
   call HDF5_closeGroup(groupHandle(1))
 
@@ -435,7 +435,7 @@ subroutine homogenization_restartRead(fileHandle)
 
     call HDF5_closeGroup(groupHandle(2))
 
-  enddo
+  end do
 
   call HDF5_closeGroup(groupHandle(1))
 
@@ -476,7 +476,7 @@ subroutine parseHomogenization
           case default
             call IO_error(500,ext_msg=homogThermal%get_asString('type'))
         end select
-    endif
+    end if
 
     if (homog%contains('damage')) then
       homogDamage => homog%get('damage')
@@ -486,8 +486,8 @@ subroutine parseHomogenization
           case default
             call IO_error(500,ext_msg=homogDamage%get_asString('type'))
         end select
-    endif
-  enddo
+    end if
+  end do
 
 end subroutine parseHomogenization
 
diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index 53d6622ac..2e66da591 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -65,9 +65,9 @@ module subroutine damage_init()
         allocate(damageState_h(ho)%state (1,Nmembers), source=1.0_pReal)
       else
         prm%output = emptyStringArray
-      endif
+      end if
     end associate
-  enddo
+  end do
 
   call pass_init()
 
@@ -174,7 +174,7 @@ module subroutine damage_results(ho,group)
             call results_writeDataset(damagestate_h(ho)%state(1,:),group,prm%output(o),&
                                       'damage indicator','-')
         end select
-      enddo outputsLoop
+      end do outputsLoop
   end associate
 
 end subroutine damage_results
diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90
index c74b9fb6c..911def035 100644
--- a/src/homogenization_mechanical_RGC.f90
+++ b/src/homogenization_mechanical_RGC.f90
@@ -561,7 +561,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
                                                 *cosh(prm%c_alpha*nDefNorm) &
                                                 *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
                                                 *tanh(nDefNorm/num%xSmoo)
-       end do; end do;enddo; end do
+       end do; end do;end do; end do
      end do interfaceLoop
 
 
diff --git a/src/lattice.f90 b/src/lattice.f90
index 1496deb78..2c2da345d 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -484,8 +484,8 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
         case default
           call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice))
       end select
-    enddo mySystems
-  enddo myFamilies
+    end do mySystems
+  end do myFamilies
 
 end function lattice_characteristicShear_Twin
 
@@ -523,7 +523,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
   do i = 1, sum(Ntwin)
     call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1)                                        ! ToDo: Why always 180 deg?
     lattice_C66_twin(1:6,1:6,i) = R%rotStiffness(C66)
-  enddo
+  end do
 
 end function lattice_C66_twin
 
@@ -572,19 +572,19 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
     C_target_unrotated66 = C_parent66
   else
     call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target))
-  endif
+  end if
 
   do i = 1,6
     if (abs(C_target_unrotated66(i,i))<tol_math_check) &
     call IO_error(135,'matrix diagonal in transformation',label1='entry',ID1=i)
-  enddo
+  end do
 
   call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_cF,a_cI)
 
   do i = 1,sum(Ntrans)
     call R%fromMatrix(Q(1:3,1:3,i))
     lattice_C66_trans(1:6,1:6,i) = R%rotStiffness(C_target_unrotated66)
-  enddo
+  end do
 
  end function lattice_C66_trans
 
@@ -632,7 +632,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
                                               math_cross(normal, direction))
     if (size(nonSchmidCoefficients)>5) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
       + nonSchmidCoefficients(6) * math_outer(direction, direction)
-  enddo
+  end do
 
 end function lattice_nonSchmidMatrix
 
@@ -1432,7 +1432,7 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
     SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
     if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) &
       error stop 'dilatational Schmid matrix for slip'
-  enddo
+  end do
 
 end function lattice_SchmidMatrix_slip
 
@@ -1479,7 +1479,7 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
     SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
     if (abs(math_trace33(SchmidMatrix(1:3,1:3,i))) > tol_math_check) &
       error stop 'dilatational Schmid matrix for twin'
-  enddo
+  end do
 
 end function lattice_SchmidMatrix_twin
 
@@ -1552,7 +1552,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMa
     SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
     SchmidMatrix(1:3,1:3,2,i) = math_outer(coordinateSystem(1:3,3,i),coordinateSystem(1:3,2,i))
     SchmidMatrix(1:3,1:3,3,i) = math_outer(coordinateSystem(1:3,2,i),coordinateSystem(1:3,2,i))
-  enddo
+  end do
 
 end function lattice_SchmidMatrix_cleavage
 
@@ -1719,8 +1719,8 @@ pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
    do i = 1, 6
      do j = i+1, 6
        C66_sym(j,i) = C66_sym(i,j)
-     enddo
-   enddo
+     end do
+   end do
 
 end function lattice_symmetrize_C66
 
@@ -1782,7 +1782,7 @@ function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
 
   do i=1, sum(Nslip); do j=1, sum(Nslip)
     projection(i,j) = abs(math_inner(n(:,i),t(:,j)))
-  enddo; enddo
+  end do; end do
 
 end function slipProjection_transverse
 
@@ -1806,7 +1806,7 @@ function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
 
   do i=1, sum(Nslip); do j=1, sum(Nslip)
     projection(i,j) = abs(math_inner(n(:,i),d(:,j)))
-  enddo; enddo
+  end do; end do
 
 end function slipProjection_direction
 
@@ -1890,8 +1890,8 @@ function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,valu
 
           buildInteraction(l,k) = values(matrix(i,j))
 
-      enddo; enddo
-  enddo; enddo
+      end do; end do
+  end do; end do
 
 end function buildInteraction
 
@@ -1957,8 +1957,8 @@ function buildCoordinateSystem(active,potential,system,lattice,cOverA)
       buildCoordinateSystem(1:3,3,a) = math_cross(direction/norm2(direction),&
                                                   normal   /norm2(normal))
 
-    enddo activeSystems
-  enddo activeFamilies
+    end do activeSystems
+  end do activeFamilies
 
 end function buildCoordinateSystem
 
@@ -2066,7 +2066,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
       U = (a_cI/a_cF) * (math_outer(x,x) + (math_outer(y,y)+math_outer(z,z)) * sqrt(2.0_pReal))
       Q(1:3,1:3,i) = matmul(R%asMatrix(),B%asMatrix())
       S(1:3,1:3,i) = matmul(R%asMatrix(),U) - MATH_I3
-    enddo
+    end do
   else if (present(cOverA)) then
     ss      = MATH_I3
     sd      = MATH_I3
@@ -2125,7 +2125,7 @@ function getlabels(active,potential,system) result(labels)
         write(label(i+1:i+2),'(I2.1)') int(system(j,p))
         label(i+3:i+3) = ' '
         i = i + 3
-      enddo direction
+      end do direction
       label(i:i) = ']'
 
       i = i +1
@@ -2134,13 +2134,13 @@ function getlabels(active,potential,system) result(labels)
         write(label(i+1:i+2),'(I2.1)') int(system(j,p))
         label(i+3:i+3) = ' '
         i = i + 3
-      enddo normal
+      end do normal
       label(i:i) = ')'
 
       labels(a) = label
 
-    enddo activeSystems
-  enddo activeFamilies
+    end do activeSystems
+  end do activeFamilies
 
 end function getlabels
 
@@ -2170,7 +2170,7 @@ pure function lattice_equivalent_nu(C,assumption) result(nu)
       / (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3)))
   else
     error stop 'invalid assumption'
-  endif
+  end if
 
   mu = lattice_equivalent_mu(C,assumption)
   nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu)
@@ -2202,7 +2202,7 @@ pure function lattice_equivalent_mu(C,assumption) result(mu)
        / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6)))
   else
     error stop 'invalid assumption'
-  endif
+  end if
 
 end function lattice_equivalent_mu
 
@@ -2266,7 +2266,7 @@ subroutine selfTest
     if (any(dNeq(T(1,1),[T_hP(1,1),T_hP(2,2)])))           error stop 'Symmetry33_11-22/hP'
     if (any(dNeq(T(1,1),[T_tI(1,1),T_tI(2,2)])))           error stop 'Symmetry33_11-22/tI'
 
-  enddo
+  end do
 
   call random_number(C)
   C(1,1) = C(1,1) + C(1,2) + 0.1_pReal
diff --git a/src/math.f90 b/src/math.f90
index 882706ea7..a814d7b99 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -135,25 +135,25 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
     s = istart
   else
     s = lbound(a,2)
-  endif
+  end if
 
   if (present(iend)) then
     e = iend
   else
     e = ubound(a,2)
-  endif
+  end if
 
   if (present(sortDim)) then
     d = sortDim
   else
     d = 1
-  endif
+  end if
 
   if (s < e) then
     call qsort_partition(a,ipivot, s,e, d)
     call math_sort(a, s, ipivot-1, d)
     call math_sort(a, ipivot+1, e, d)
-  endif
+  end if
 
 
   contains
@@ -175,11 +175,11 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
       ! find the first element on the right side less than or equal to the pivot point
       do j = iend, istart, -1
         if (a(sort,j) <= a(sort,istart)) exit
-      enddo
+      end do
       ! find the first element on the left side greater than the pivot point
       do i = istart, iend
         if (a(sort,i) > a(sort,istart)) exit
-      enddo
+      end do
       cross: if (i >= j) then ! exchange left value with pivot and return with the partition index
         tmp         = a(:,istart)
         a(:,istart) = a(:,j)
@@ -190,8 +190,8 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
         tmp    = a(:,i)
         a(:,i) = a(:,j)
         a(:,j) = tmp
-      endif cross
-    enddo
+      end if cross
+    end do
 
   end subroutine qsort_partition
 
@@ -216,7 +216,7 @@ pure function math_expand(what,how)
 
   do i = 1, size(how)
     math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1)
-  enddo
+  end do
 
 end function math_expand
 
@@ -251,7 +251,7 @@ pure function math_eye(d)
   math_eye = 0.0_pReal
   do i=1,d
     math_eye(i,i) = 1.0_pReal
-  enddo
+  end do
 
 end function math_eye
 
@@ -270,7 +270,7 @@ pure function math_identity4th()
 #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
+  end do
 #else
   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))
@@ -298,7 +298,7 @@ real(pReal) pure function math_LeviCivita(i,j,k)
     math_LeviCivita = -1.0_pReal
   else
     math_LeviCivita =  0.0_pReal
-  endif
+  end if
 
 end function math_LeviCivita
 
@@ -348,7 +348,7 @@ pure function math_outer(A,B)
 #ifndef __INTEL_COMPILER
   do concurrent(i=1:size(A,1), j=1:size(B,1))
     math_outer(i,j) = A(i)*B(j)
-  enddo
+  end do
 #else
   forall(i=1:size(A,1), j=1:size(B,1)) math_outer(i,j) = A(i)*B(j)
 #endif
@@ -398,7 +398,7 @@ pure function math_mul3333xx33(A,B)
 #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
+  end do
 #else
   forall (i=1:3, j=1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
 #endif
@@ -421,7 +421,7 @@ pure function math_mul3333xx3333(A,B)
 #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
+  end do
 #else
   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))
 #endif
@@ -446,7 +446,7 @@ pure function math_exp33(A,n)
     n_ = n
   else
     n_ = 5
-  endif
+  end if
 
   invFac     = 1.0_pReal                                                                            ! 0!
   B          = math_I3
@@ -456,7 +456,7 @@ pure function math_exp33(A,n)
     invFac = invFac/real(i,pReal)                                                                   ! invfac = 1/(i!)
     B = matmul(B,A)
     math_exp33 = math_exp33 + invFac*B                                                              ! exp = SUM (A^i)/(i!)
-  enddo
+  end do
 
 end function math_exp33
 
@@ -514,7 +514,7 @@ pure subroutine math_invert33(InvA, DetA, error, A)
 
     InvA = InvA/DetA
     error = .false.
-  endif
+  end if
 
 end subroutine math_invert33
 
@@ -541,7 +541,7 @@ pure function math_invSym3333(A)
     error stop 'matrix inversion error'
   else
     math_invSym3333 = math_66toSym3333(temp66)
-  endif
+  end if
 
 end function math_invSym3333
 
@@ -696,7 +696,7 @@ pure function math_9to33(v9)
 
   do i = 1, 9
     math_9to33(MAPPLAIN(1,i),MAPPLAIN(2,i)) = v9(i)
-  enddo
+  end do
 
 end function math_9to33
 
@@ -721,7 +721,7 @@ pure function math_sym33to6(m33,weighted)
     w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
     w = NRMMANDEL
-  endif
+  end if
 
   math_sym33to6 = [(w(i)*m33(MAPNYE(1,i),MAPNYE(2,i)),i=1,6)]
 
@@ -748,12 +748,12 @@ pure function math_6toSym33(v6,weighted)
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
     w = INVNRMMANDEL
-  endif
+  end if
 
   do i=1,6
     math_6toSym33(MAPNYE(1,i),MAPNYE(2,i)) = w(i)*v6(i)
     math_6toSym33(MAPNYE(2,i),MAPNYE(1,i)) = w(i)*v6(i)
-  enddo
+  end do
 
 end function math_6toSym33
 
@@ -772,7 +772,7 @@ pure function math_3333to99(m3333)
 #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
+  end do
 #else
   forall(i=1:9, j=1:9) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
 #endif
@@ -794,7 +794,7 @@ pure function math_99to3333(m99)
 #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
+  end do
 #else
   forall(i=1:9, j=1:9) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
 #endif
@@ -827,7 +827,7 @@ pure function math_sym3333to66(m3333,weighted)
 #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
+  end do
 #else
   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))
 #endif
@@ -1080,8 +1080,8 @@ pure subroutine math_eigh33(w,v,m)
     else fallback2
       v(1:3,2) = v(1:3, 2) / norm
       v(1:3,3) = math_cross(v(1:3,1),v(1:3,2))
-     endif fallback2
-  endif fallback1
+     end if fallback2
+  end if fallback1
 
 end subroutine math_eigh33
 
@@ -1120,7 +1120,7 @@ pure function math_rotationalPart(F) result(R)
     lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal)
   else
     lambda = sqrt(I_C(1)/3.0_pReal)
-  endif
+  end if
 
   I_U = [sum(lambda), lambda(1)*lambda(2)+lambda(2)*lambda(3)+lambda(3)*lambda(1), product(lambda)]
 
@@ -1188,7 +1188,7 @@ pure function math_eigvalsh33(m)
                                                              cos((phi+2.0_pReal*TAU)/3.0_pReal) &
                                                             ] &
                     + I(1)/3.0_pReal
-  endif
+  end if
 
 end function math_eigvalsh33
 
@@ -1238,7 +1238,7 @@ integer pure function math_binomial(n,k)
   do i = 1, k_
     math_binomial = (math_binomial * n_)/i
     n_ = n_ -1
-  enddo
+  end do
 
 end function math_binomial
 
@@ -1302,7 +1302,7 @@ real(pReal) pure elemental function math_clip(a, left, right)
   if (present(right)) math_clip = min(right,math_clip)
   if (present(left) .and. present(right)) then
     if (left>right) error stop 'left > right'
-  endif
+  end if
 
 end function math_clip
 
@@ -1402,7 +1402,7 @@ subroutine selfTest()
 
   do while(abs(math_det33(t33))<1.0e-9_pReal)
     call random_number(t33)
-  enddo
+  end do
   if (any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) &
     error stop 'math_inv33'
 
@@ -1418,7 +1418,7 @@ subroutine selfTest()
 
   do while(math_det33(t33)<1.0e-2_pReal)                                                            ! O(det(F)) = 1
     call random_number(t33)
-  enddo
+  end do
   t33_2 = math_rotationalPart(transpose(t33))
   t33   = math_rotationalPart(t33)
   if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) &
diff --git a/src/parallelization.f90 b/src/parallelization.f90
index 44e5c86ae..d631711d8 100644
--- a/src/parallelization.f90
+++ b/src/parallelization.f90
@@ -97,7 +97,7 @@ subroutine parallelization_init
     open(OUTPUT_UNIT,file='/dev/null',status='replace')                                             ! close() alone will leave some temp files in cwd
   else
     open(OUTPUT_UNIT,encoding='UTF-8')                                                              ! for special characters in output
-  endif
+  end if
 #endif
 
   print'(/,1x,a)', '<<<+-  parallelization init  -+>>>'
@@ -142,8 +142,8 @@ subroutine parallelization_init
 !$   if (OMP_NUM_THREADS < 1_pI32) then
 !$     print'(1x,a)', 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default'
 !$     OMP_NUM_THREADS = 4_pI32
-!$   endif
-!$ endif
+!$   end if
+!$ end if
 !$ print'(1x,a,i0)',   'OMP_NUM_THREADS: ',OMP_NUM_THREADS
 !$ call omp_set_num_threads(OMP_NUM_THREADS)
 
diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index afde549fc..93357bee5 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -255,7 +255,7 @@ module subroutine mechanical_init(phases)
 #else
     output_mechanical(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
-  enddo
+  end do
 
   do ce = 1, size(material_phaseID,2)
     ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
@@ -267,14 +267,14 @@ module subroutine mechanical_init(phases)
       phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), &
                                                         transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
       phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co)))
-    enddo
-  enddo
+    end do
+  end do
 
   do ph = 1, phases%length
     phase_mechanical_F0(ph)%data  = phase_mechanical_F(ph)%data
     phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data
     phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(ph)%data
-  enddo
+  end do
 
 
   call elastic_init(phases)
@@ -284,7 +284,7 @@ module subroutine mechanical_init(phases)
   call plastic_init()
   do ph = 1,phases%length
     plasticState(ph)%state0 = plasticState(ph)%state
-  enddo
+  end do
 
   num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
 
@@ -473,25 +473,25 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
         Lpguess      = Lpguess_old &
                      + deltaLp * stepLengthLp
         cycle LpLoop
-      endif
+      end if
 
       calculateJacobiLp: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then
         jacoCounterLp = jacoCounterLp + 1
 
         do o=1,3; do p=1,3
           dFe_dLp(o,1:3,p,1:3) = - Delta_t * A(o,p)*transpose(invFi_new)                            ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j)
-        enddo; enddo
+        end do; end do
         dRLp_dLp = math_eye(9) &
                  - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp))
         temp_9 = math_33to9(residuumLp)
         call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr)                                          ! solve dRLp/dLp * delta Lp = -res for delta Lp
         if (ierr /= 0) return ! error
         deltaLp = - math_9to33(temp_9)
-      endif calculateJacobiLp
+      end if calculateJacobiLp
 
       Lpguess = Lpguess &
               + deltaLp * steplengthLp
-    enddo LpLoop
+    end do LpLoop
 
     call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
                                 S, Fi_new, ph,en)
@@ -513,7 +513,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
       Liguess      = Liguess_old &
                    + deltaLi * steplengthLi
       cycle LiLoop
-    endif
+    end if
 
     calculateJacobiLi: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then
       jacoCounterLi = jacoCounterLi + 1
@@ -522,10 +522,10 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
       do o=1,3; do p=1,3
         dFe_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*temp_33                                        ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j)
         dFi_dLi(1:3,o,1:3,p) = -Delta_t*math_I3(o,p)*invFi_current
-      enddo; enddo
+      end do; end do
       do o=1,3; do p=1,3
         dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
-      enddo; enddo
+      end do; end do
       dRLi_dLi  = math_eye(9) &
                 - math_3333to99(math_mul3333xx3333(dLi_dS,  math_mul3333xx3333(dS_dFe, dFe_dLi) &
                                                           + math_mul3333xx3333(dS_dFi, dFi_dLi)))  &
@@ -534,11 +534,11 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
       call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr)                                            ! solve dRLi/dLp * delta Li = -res for delta Li
       if (ierr /= 0) return ! error
       deltaLi = - math_9to33(temp_9)
-    endif calculateJacobiLi
+    end if calculateJacobiLi
 
     Liguess = Liguess &
             + deltaLi * steplengthLi
-  enddo LiLoop
+  end do LiLoop
 
   invFp_new = matmul(invFp_current,B)
   call math_invert33(Fp_new,devNull,error,invFp_new)
@@ -613,9 +613,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
     if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then
       broken = plastic_deltaState(ph,en)
       exit iteration
-    endif
+    end if
 
-  enddo iteration
+  end do iteration
 
 
   contains
@@ -638,7 +638,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(b
       damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
     else
       damper = 1.0_pReal
-    endif
+    end if
 
   end function damper
 
@@ -844,7 +844,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
 #else
       dotState = IEEE_FMA(A(n,stage),plastic_RKdotState(1:sizeDotState,n),dotState)
 #endif
-    enddo
+    end do
 
 #ifndef __INTEL_LLVM_COMPILER
     plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState*Delta_t
@@ -858,7 +858,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
     dotState = plastic_dotState(Delta_t*C(stage), ph,en)
     if (any(IEEE_is_NaN(dotState))) exit
 
-  enddo
+  end do
   if(broken) return
 
 
@@ -950,7 +950,7 @@ subroutine results(group,ph)
 
     do i = 1, size(dataset,1)
       to_quaternion(:,i) = dataset(i)%asQuaternion()
-    enddo
+    end do
 
  end function to_quaternion
 
@@ -974,7 +974,7 @@ module subroutine mechanical_forward()
     phase_mechanical_Lp0(ph) = phase_mechanical_Lp(ph)
     phase_mechanical_S0(ph)  = phase_mechanical_S(ph)
     plasticState(ph)%state0 = plasticState(ph)%state
-  enddo
+  end do
 
 end subroutine mechanical_forward
 
@@ -1037,7 +1037,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
         subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en)
         subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en)
         subState0 = plasticState(ph)%state(:,en)
-      endif
+      end if
 !--------------------------------------------------------------------------------------------------
 !  cut back (reduced time and restore)
     else
@@ -1048,10 +1048,10 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
       if (subStep < 1.0_pReal) then                                                                 ! actual (not initial) cutback
         phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
         phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
-      endif
+      end if
       plasticState(ph)%state(:,en) = subState0
       todo = subStep > num%subStepMinCryst                                                          ! still on track or already done (beyond repair)
-    endif
+    end if
 
 !--------------------------------------------------------------------------------------------------
 !  prepare for integration
@@ -1060,9 +1060,9 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
       subF = subF0 &
            + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
       converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en)
-    endif
+    end if
 
-  enddo cutbackLooping
+  end do cutbackLooping
 
 end function phase_mechanical_constitutive
 
@@ -1086,14 +1086,14 @@ module subroutine mechanical_restore(ce,includeL)
     if (includeL) then
       phase_mechanical_Lp(ph)%data(1:3,1:3,en) = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
       phase_mechanical_Li(ph)%data(1:3,1:3,en) = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
-    endif                                                                                           ! maybe protecting everything from overwriting makes more sense
+    end if                                                                                           ! maybe protecting everything from overwriting makes more sense
 
     phase_mechanical_Fp(ph)%data(1:3,1:3,en)   = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
     phase_mechanical_Fi(ph)%data(1:3,1:3,en)   = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
     phase_mechanical_S(ph)%data(1:3,1:3,en)    = phase_mechanical_S0(ph)%data(1:3,1:3,en)
 
     plasticState(ph)%state(:,en) = plasticState(ph)%State0(:,en)
-  enddo
+  end do
 
 end subroutine mechanical_restore
 
@@ -1164,7 +1164,7 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
       lhs_3333(1:3,o,1:3,p) = IEEE_FMA(invFi,invFi(p,o),lhs_3333(1:3,o,1:3,p))
       rhs_3333(1:3,1:3,o,p) = IEEE_FMA(matmul(invSubFi0,dLidS(1:3,1:3,o,p)),-Delta_t,rhs_3333(1:3,1:3,o,p))
 #endif
-    enddo; enddo
+    end do; end do
     call math_invert(temp_99,error,math_3333to99(lhs_3333))
     if (error) then
       call IO_warning(600,'inversion error in analytic tangent calculation', &
@@ -1172,9 +1172,9 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
       dFidS = 0.0_pReal
     else
       dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
-    endif
+    end if
     dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
-  endif
+  end if
 
   call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
                                              phase_mechanical_S(ph)%data(1:3,1:3,en), &
@@ -1191,7 +1191,7 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
     rhs_3333(p,o,1:3,1:3)  = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
     temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
                            + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
-  enddo; enddo
+  end do; end do
 #ifndef __INTEL_LLVM_COMPILER
   lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t &
            + math_mul3333xx3333(dSdFi,dFidS)
@@ -1206,14 +1206,14 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
     dSdF = rhs_3333
   else
     dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
-  endif
+  end if
 
 !--------------------------------------------------------------------------------------------------
 ! calculate dFpinvdF
   temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
   do o=1,3; do p=1,3
     dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t
-  enddo; enddo
+  end do; end do
 
 !--------------------------------------------------------------------------------------------------
 ! assemble dPdF
@@ -1224,13 +1224,13 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
   dPdF = 0.0_pReal
   do p=1,3
     dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1))
-  enddo
+  end do
   do o=1,3; do p=1,3
     dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
                       + matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,en),dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
                       + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) &
                       + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
-  enddo; enddo
+  end do; end do
 
 end function phase_mechanical_dPdF
 
diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90
index f0b86319c..4cd02d090 100644
--- a/src/phase_mechanical_eigen.f90
+++ b/src/phase_mechanical_eigen.f90
@@ -61,7 +61,7 @@ module subroutine eigen_init(phases)
 
   if (maxval(Nmodels) /= 0) then
     where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID
-  endif
+  end if
 
   allocate(model_damage(phases%length),  source = EIGEN_UNDEFINED_ID)
 
diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90
index 15a8ae2a7..6469c6e45 100644
--- a/src/phase_mechanical_plastic.f90
+++ b/src/phase_mechanical_plastic.f90
@@ -435,7 +435,7 @@ function plastic_active(plastic_label) result(active_plastic)
     mech  => phase%get('mechanical')
     pl    => mech%get('plastic',defaultVal = emptyDict)
     active_plastic(ph) = pl%get_asString('type',defaultVal='none') == plastic_label
-  enddo
+  end do
 
 end function plastic_active
 
diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90
index a29efeb56..dfd84ffc9 100644
--- a/src/phase_mechanical_plastic_dislotwin.f90
+++ b/src/phase_mechanical_plastic_dislotwin.f90
@@ -9,28 +9,28 @@
 !--------------------------------------------------------------------------------------------------
 submodule(phase:plastic) dislotwin
 
+  real(pReal), parameter :: gamma_char_tr = sqrt(0.125_pReal)                                       !< Characteristic shear for transformation
   type :: tParameters
     real(pReal) :: &
-      Q_cl                = 1.0_pReal, &                                                            !< activation energy for dislocation climb
-      omega               = 1.0_pReal, &                                                            !< frequency factor for dislocation climb
-      D                   = 1.0_pReal, &                                                            !< grain size
-      p_sb                = 1.0_pReal, &                                                            !< p-exponent in shear band velocity
-      q_sb                = 1.0_pReal, &                                                            !< q-exponent in shear band velocity
-      i_tw                = 1.0_pReal, &                                                            !< adjustment parameter to calculate MFP for twinning
-      i_tr                = 1.0_pReal, &                                                            !< adjustment parameter to calculate MFP for transformation
-      L_tw                = 1.0_pReal, &                                                            !< length of twin nuclei
-      L_tr                = 1.0_pReal, &                                                            !< length of trans nuclei
-      x_c                 = 1.0_pReal, &                                                            !< critical distance for formation of twin/trans nucleus
-      V_cs                = 1.0_pReal, &                                                            !< cross slip volume
-      xi_sb               = 1.0_pReal, &                                                            !< value for shearband resistance
-      v_sb                = 1.0_pReal, &                                                            !< value for shearband velocity_0
-      E_sb                = 1.0_pReal, &                                                            !< activation energy for shear bands
-      h                   = 1.0_pReal, &                                                            !< stack height of hex nucleus
-      gamma_char_tr       = sqrt(0.125_pReal), &                                                    !< Characteristic shear for transformation
-      a_cF                = 1.0_pReal, &
-      cOverA_hP           = 1.0_pReal, &
-      V_mol               = 1.0_pReal, &
-      rho                 = 1.0_pReal
+      Q_cl       = 1.0_pReal, &                                                                     !< activation energy for dislocation climb
+      omega      = 1.0_pReal, &                                                                     !< frequency factor for dislocation climb
+      D          = 1.0_pReal, &                                                                     !< grain size
+      p_sb       = 1.0_pReal, &                                                                     !< p-exponent in shear band velocity
+      q_sb       = 1.0_pReal, &                                                                     !< q-exponent in shear band velocity
+      i_tw       = 1.0_pReal, &                                                                     !< adjustment parameter to calculate MFP for twinning
+      i_tr       = 1.0_pReal, &                                                                     !< adjustment parameter to calculate MFP for transformation
+      L_tw       = 1.0_pReal, &                                                                     !< length of twin nuclei
+      L_tr       = 1.0_pReal, &                                                                     !< length of trans nuclei
+      x_c        = 1.0_pReal, &                                                                     !< critical distance for formation of twin/trans nucleus
+      V_cs       = 1.0_pReal, &                                                                     !< cross slip volume
+      tau_sb     = 1.0_pReal, &                                                                     !< value for shearband resistance
+      gamma_0_sb = 1.0_pReal, &                                                                     !< value for shearband velocity_0
+      E_sb       = 1.0_pReal, &                                                                     !< activation energy for shear bands
+      h          = 1.0_pReal, &                                                                     !< stack height of hex nucleus
+      a_cF       = 1.0_pReal, &
+      cOverA_hP  = 1.0_pReal, &
+      V_mol      = 1.0_pReal, &
+      rho        = 1.0_pReal
     type(tPolynomial) :: &
       Gamma_sf, &                                                                                   !< stacking fault energy
       Delta_G                                                                                       !< free energy difference between austensite and martensite
@@ -331,18 +331,18 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! shearband related parameters
-    prm%v_sb = pl%get_asFloat('v_sb',defaultVal=0.0_pReal)
-    if (prm%v_sb > 0.0_pReal) then
-      prm%xi_sb        = pl%get_asFloat('xi_sb')
+    prm%gamma_0_sb = pl%get_asFloat('gamma_0_sb',defaultVal=0.0_pReal)
+    if (prm%gamma_0_sb > 0.0_pReal) then
+      prm%tau_sb       = pl%get_asFloat('tau_sb')
       prm%E_sb         = pl%get_asFloat('Q_sb')
       prm%p_sb         = pl%get_asFloat('p_sb')
       prm%q_sb         = pl%get_asFloat('q_sb')
 
       ! sanity checks
-      if (prm%xi_sb         <  0.0_pReal) extmsg = trim(extmsg)//' xi_sb'
-      if (prm%E_sb          <  0.0_pReal) extmsg = trim(extmsg)//' Q_sb'
-      if (prm%p_sb          <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb'
-      if (prm%q_sb          <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb'
+      if (prm%tau_sb <  0.0_pReal) extmsg = trim(extmsg)//' tau_sb'
+      if (prm%E_sb   <  0.0_pReal) extmsg = trim(extmsg)//' Q_sb'
+      if (prm%p_sb   <= 0.0_pReal) extmsg = trim(extmsg)//' p_sb'
+      if (prm%q_sb   <= 0.0_pReal) extmsg = trim(extmsg)//' q_sb'
     end if
 
 !--------------------------------------------------------------------------------------------------
@@ -364,13 +364,13 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), &
                                                    phase_lattice(ph))
       if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
-    endif slipAndTwinActive
+    end if slipAndTwinActive
 
     slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
       prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dFloat('h_sl-tr'), &
                                                     phase_lattice(ph))
       if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
-    endif slipAndTransActive
+    end if slipAndTransActive
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -569,7 +569,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
     Lp      = Lp      * f_matrix
     dLp_dMp = dLp_dMp * f_matrix
 
-    shearBandingContribution: if (dNeq0(prm%v_sb)) then
+    shearBandingContribution: if (dNeq0(prm%gamma_0_sb)) then
 
       E_kB_T = prm%E_sb/(K_B*T)
       call math_eigh33(eigValues,eigVectors,Mp)                                                     ! is Mp symmetric by design?
@@ -580,10 +580,10 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
         tau = math_tensordot(Mp,P_sb)
 
         significantShearBandStress: if (abs(tau) > tol_math_check) then
-          StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
-          dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
-          ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
-                          * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
+          StressRatio_p = (abs(tau)/prm%tau_sb)**prm%p_sb
+          dot_gamma_sb = sign(prm%gamma_0_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
+          ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%tau_sb &
+                          * (abs(tau)/prm%tau_sb)**(prm%p_sb-1.0_pReal) &
                           * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
 
           Lp = Lp + dot_gamma_sb * P_sb
@@ -697,7 +697,7 @@ module function dislotwin_dotState(Mp,ph,en) result(dotState)
     dot_f_tw = f_matrix*dot_gamma_tw/prm%gamma_char_tw
 
     if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,dot_gamma_tr)
-    dot_f_tr = f_matrix*dot_gamma_tr/prm%gamma_char_tr
+    dot_f_tr = f_matrix*dot_gamma_tr/gamma_char_tr
 
   end associate
 
@@ -1026,9 +1026,9 @@ pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,&
         dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
 
         V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i)
-        dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*prm%gamma_char_tr
+        dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*gamma_char_tr
         if (present(ddot_gamma_dtau_tr)) &
-          ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*prm%gamma_char_tr
+          ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*gamma_char_tr
       else
         dot_gamma_tr(i) = 0.0_pReal
         if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index 2aab37529..0f5836df9 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -122,25 +122,25 @@ module subroutine thermal_init(phases)
 
     allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
 
-  enddo
+  end do
 
   allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
 
   if (maxval(thermal_Nsources) /= 0) then
     where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
     where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
-  endif
+  end if
 
   thermal_source_maxSizeDotState = 0
   do ph = 1,phases%length
 
     do so = 1,thermal_Nsources(ph)
       thermalState(ph)%p(so)%state  = thermalState(ph)%p(so)%state0
-    enddo
+    end do
 
     thermal_source_maxSizeDotState  = max(thermal_source_maxSizeDotState, &
                                           maxval(thermalState(ph)%p%sizeDotState))
-  enddo
+  end do
 
 end subroutine thermal_init
 
@@ -170,7 +170,7 @@ module function phase_f_T(ph,en) result(f)
 
    end select
 
-  enddo
+  end do
 
 end function phase_f_T
 
@@ -195,7 +195,7 @@ function phase_thermal_collectDotState(ph,en) result(broken)
 
     broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))
 
-  enddo SourceLoop
+  end do SourceLoop
 
 end function phase_thermal_collectDotState
 
@@ -262,7 +262,7 @@ function integrateThermalState(Delta_t, ph,en) result(broken)
     sizeDotState = thermalState(ph)%p(so)%sizeDotState
     thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) &
                                                     + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t
-  enddo
+  end do
 
 end function integrateThermalState
 
@@ -277,7 +277,7 @@ module subroutine thermal_restartWrite(groupHandle,ph)
 
   do so = 1,thermal_Nsources(ph)
     call HDF5_write(thermalState(ph)%p(so)%state,groupHandle,'omega_thermal')
-  enddo
+  end do
 
 end subroutine thermal_restartWrite
 
@@ -292,7 +292,7 @@ module subroutine thermal_restartRead(groupHandle,ph)
 
   do so = 1,thermal_Nsources(ph)
     call HDF5_read(thermalState(ph)%p(so)%state0,groupHandle,'omega_thermal')
-  enddo
+  end do
 
 end subroutine thermal_restartRead
 
@@ -305,8 +305,8 @@ module subroutine thermal_forward()
   do ph = 1, size(thermalState)
     do so = 1, size(thermalState(ph)%p)
       thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state
-    enddo
-  enddo
+    end do
+  end do
 
 end subroutine thermal_forward
 
@@ -380,8 +380,8 @@ function thermal_active(source_label,src_length)  result(active_source)
     do s = 1, sources%length
       src => sources%get(s)
       active_source(s,p) = src%get_asString('type') == source_label
-    enddo
-  enddo
+    end do
+  end do
 
 
 end function thermal_active
diff --git a/src/polynomials.f90 b/src/polynomials.f90
index c883528ca..e04db16f5 100644
--- a/src/polynomials.f90
+++ b/src/polynomials.f90
@@ -112,7 +112,7 @@ pure function eval(self,x) result(y)
 #else
     y = IEEE_FMA(y,x-self%x_ref,self%coef(o))
 #endif
-  enddo
+  end do
 
 end function eval
 
diff --git a/src/results.f90 b/src/results.f90
index bd7f7d34d..cf1b4de5c 100644
--- a/src/results.f90
+++ b/src/results.f90
@@ -95,11 +95,11 @@ subroutine results_init(restart)
     call results_openJobFile
     call get_command(commandLine)
     call results_addAttribute('call (restart at '//date//')',trim(commandLine))
-    call h5gmove_f(resultsFile,'setup','tmp',hdferr)
+    call H5Gmove_f(resultsFile,'setup','tmp',hdferr)
     call results_addAttribute('description','input data used to run the simulation up to restart at '//date,'tmp')
     call results_closeGroup(results_addGroup('setup'))
     call results_addAttribute('description','input data used to run the simulation','setup')
-    call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr)
+    call H5Gmove_f(resultsFile,'tmp','setup/previous',hdferr)
   end if
 
   call results_closeJobFile
@@ -333,8 +333,8 @@ subroutine results_removeLink(link)
   integer                      :: hdferr
 
 
-  call h5ldelete_f(resultsFile,link, hdferr)
-  if (hdferr < 0) call IO_error(1,ext_msg = 'results_removeLink: h5ldelete_soft_f ('//trim(link)//')')
+  call H5Ldelete_f(resultsFile,link, hdferr)
+  if (hdferr < 0) call IO_error(1,ext_msg = 'results_removeLink: H5Ldelete_soft_f ('//trim(link)//')')
 
 end subroutine results_removeLink
 
@@ -522,7 +522,7 @@ subroutine results_mapping_phase(ID,entry,label)
   writeSize = 0
   writeSize(worldrank) = size(entry(1,:))                                                           ! total number of entries of this process
 
-  call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+  call H5Pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 #ifndef PETSC
@@ -530,7 +530,7 @@ subroutine results_mapping_phase(ID,entry,label)
 #else
 !--------------------------------------------------------------------------------------------------
 ! MPI settings and communication
-  call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+  call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)   ! get output at each process
@@ -558,82 +558,82 @@ subroutine results_mapping_phase(ID,entry,label)
 
 !---------------------------------------------------------------------------------------------------
 ! compound type: label(ID) + entry
-  call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
+  call H5Tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
+  call H5Tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tget_size_f(dt_id, type_size_string, hdferr)
+  call H5Tget_size_f(dt_id, type_size_string, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   pI64_t = h5kind_to_type(kind(entryGlobal),H5_INTEGER_KIND)
-  call h5tget_size_f(pI64_t, type_size_int, hdferr)
+  call H5Tget_size_f(pI64_t, type_size_int, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
+  call H5Tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
+  call H5Tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! create memory types for each component of the compound type
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
+  call H5Tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
+  call H5Tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tclose_f(dt_id, hdferr)
+  call H5Tclose_f(dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! create dataspace in memory (local shape = hyperslab) and in file (global shape)
-  call h5screate_simple_f(2,myShape,memspace_id,hdferr,myShape)
+  call H5Screate_simple_f(2,myShape,memspace_id,hdferr,myShape)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5screate_simple_f(2,totalShape,filespace_id,hdferr,totalShape)
+  call H5Screate_simple_f(2,totalShape,filespace_id,hdferr,totalShape)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
+  call H5Sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! write the components of the compound type individually
-  call h5pset_preserve_f(plist_id, .true., hdferr)
+  call H5Pset_preserve_f(plist_id, .true., hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   loc_id = results_openGroup('/cell_to')
-  call h5dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr)
+  call H5Dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
+  call H5Dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
+  call H5Dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! close all
   call HDF5_closeGroup(loc_id)
-  call h5pclose_f(plist_id, hdferr)
+  call H5Pclose_f(plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5sclose_f(filespace_id, hdferr)
+  call H5Sclose_f(filespace_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5sclose_f(memspace_id, hdferr)
+  call H5Sclose_f(memspace_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dclose_f(dset_id, hdferr)
+  call H5Dclose_f(dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(dtype_id, hdferr)
+  call H5Tclose_f(dtype_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(label_id, hdferr)
+  call H5Tclose_f(label_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(entry_id, hdferr)
+  call H5Tclose_f(entry_id, hdferr)
 
   call executionStamp('cell_to/phase','cell ID and constituent ID to phase results')
 
@@ -678,7 +678,7 @@ subroutine results_mapping_homogenization(ID,entry,label)
   writeSize = 0
   writeSize(worldrank) = size(entry)                                                                ! total number of entries of this process
 
-  call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
+  call H5Pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 #ifndef PETSC
@@ -686,7 +686,7 @@ subroutine results_mapping_homogenization(ID,entry,label)
 #else
 !--------------------------------------------------------------------------------------------------
 ! MPI settings and communication
-  call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+  call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)   ! get output at each process
@@ -710,82 +710,82 @@ subroutine results_mapping_homogenization(ID,entry,label)
 
 !---------------------------------------------------------------------------------------------------
 ! compound type: label(ID) + entry
-  call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
+  call H5Tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
+  call H5Tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tget_size_f(dt_id, type_size_string, hdferr)
+  call H5Tget_size_f(dt_id, type_size_string, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   pI64_t = h5kind_to_type(kind(entryGlobal),H5_INTEGER_KIND)
-  call h5tget_size_f(pI64_t, type_size_int, hdferr)
+  call H5Tget_size_f(pI64_t, type_size_int, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_string + type_size_int, dtype_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
+  call H5Tinsert_f(dtype_id, 'label', 0_SIZE_T, dt_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
+  call H5Tinsert_f(dtype_id, 'entry', type_size_string, pI64_t, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! create memory types for each component of the compound type
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_string, label_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
+  call H5Tinsert_f(label_id, 'label', 0_SIZE_T, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
+  call H5Tcreate_f(H5T_COMPOUND_F, type_size_int, entry_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
+  call H5Tinsert_f(entry_id, 'entry', 0_SIZE_T, pI64_t, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5tclose_f(dt_id, hdferr)
+  call H5Tclose_f(dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! create dataspace in memory (local shape = hyperslab) and in file (global shape)
-  call h5screate_simple_f(1,myShape,memspace_id,hdferr,myShape)
+  call H5Screate_simple_f(1,myShape,memspace_id,hdferr,myShape)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5screate_simple_f(1,totalShape,filespace_id,hdferr,totalShape)
+  call H5Screate_simple_f(1,totalShape,filespace_id,hdferr,totalShape)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
+  call H5Sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! write the components of the compound type individually
-  call h5pset_preserve_f(plist_id, .true., hdferr)
+  call H5Pset_preserve_f(plist_id, .true., hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   loc_id = results_openGroup('/cell_to')
-  call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr)
+  call H5Dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
+  call H5Dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
+  call H5Dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
 
 !--------------------------------------------------------------------------------------------------
 ! close all
   call HDF5_closeGroup(loc_id)
-  call h5pclose_f(plist_id, hdferr)
+  call H5Pclose_f(plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5sclose_f(filespace_id, hdferr)
+  call H5Sclose_f(filespace_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5sclose_f(memspace_id, hdferr)
+  call H5Sclose_f(memspace_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dclose_f(dset_id, hdferr)
+  call H5Dclose_f(dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(dtype_id, hdferr)
+  call H5Tclose_f(dtype_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(label_id, hdferr)
+  call H5Tclose_f(label_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5tclose_f(entry_id, hdferr)
+  call H5Tclose_f(entry_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   call executionStamp('cell_to/homogenization','cell ID to homogenization results')
diff --git a/src/rotations.f90 b/src/rotations.f90
index d548404c9..f6ef57783 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -182,7 +182,7 @@ subroutine fromEulers(self,eu,degrees)
     Eulers = eu
   else
     Eulers = merge(eu*INRAD,eu,degrees)
-  endif
+  end if
 
   if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) &
     call IO_error(402,ext_msg='fromEulers')
@@ -206,14 +206,14 @@ subroutine fromAxisAngle(self,ax,degrees,P)
     angle = ax(4)
   else
     angle = merge(ax(4)*INRAD,ax(4),degrees)
-  endif
+  end if
 
   if (.not. present(P)) then
     axis = ax(1:3)
   else
     axis = ax(1:3) * merge(-1.0_pReal,1.0_pReal,P == 1)
     if(abs(P) /= 1) call IO_error(402,ext_msg='fromAxisAngle (P)')
-  endif
+  end if
 
   if(dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) &
     call IO_error(402,ext_msg='fromAxisAngle')
@@ -284,7 +284,7 @@ pure function rotVector(self,v,active) result(vRot)
     passive = .not. active
   else
     passive = .true.
-  endif
+  end if
 
   if (dEq0(norm2(v))) then
     vRot = v
@@ -294,7 +294,7 @@ pure function rotVector(self,v,active) result(vRot)
               multiplyQuaternion(conjugateQuaternion(self%q), multiplyQuaternion(v_normed, self%q)), &
               passive)
     vRot = q(2:4)*norm2(v)
-  endif
+  end if
 
 end function rotVector
 
@@ -318,7 +318,7 @@ pure function rotTensor2(self,T,active) result(tRot)
     passive = .not. active
   else
     passive = .true.
-  endif
+  end if
 
   tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
                matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
@@ -347,14 +347,14 @@ pure function rotTensor4(self,T,active) result(tRot)
     R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
   else
     R = self%asMatrix()
-  endif
+  end if
 
   tRot = 0.0_pReal
   do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3
   do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3
     tRot(i,j,k,l) = tRot(i,j,k,l) &
                   + R(i,m) * R(j,n) * R(k,o) * R(l,p) * T(m,n,o,p)
-  enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
+  end do; end do; end do; end do; end do; end do; end do; end do
 
 end function rotTensor4
 
@@ -379,7 +379,7 @@ pure function rotStiffness(self,C,active) result(cRot)
     R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
   else
     R = self%asMatrix()
-  endif
+  end if
 
   M = reshape([R(1,1)**2,                   R(2,1)**2,                   R(3,1)**2, &
                R(2,1)*R(3,1),               R(1,1)*R(3,1),               R(1,1)*R(2,1), &
@@ -468,7 +468,7 @@ pure function qu2eu(qu) result(eu)
     eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), &
           atan2( 2.0_pReal*chi, q03-q12 ), &
           atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
-  endif degenerated
+  end if degenerated
   where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
 
 end function qu2eu
@@ -526,8 +526,8 @@ pure function om2qu(om) result(qu)
       else
           s = 2.0_pReal * sqrt( 1.0_pReal + om(3,3) - om(1,1) - om(2,2) )
           qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
-      endif
-  endif
+      end if
+  end if
   if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
   qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4)))
   qu = qu/norm2(qu)
@@ -593,7 +593,7 @@ function om2ax(om) result(ax)
     ax(1:3) = VR(1:3,i)
     where (                dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
       ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])
-  endif
+  end if
 
 end function om2ax
 
@@ -826,7 +826,7 @@ subroutine selfTest()
             cos(TAU*x(2))*B,&
             sin(TAU*x(1))*A]
       if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
-    endif
+    end if
 
 
     if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu2om'
diff --git a/src/system_routines.f90 b/src/system_routines.f90
index 209013be4..a74c2e417 100644
--- a/src/system_routines.f90
+++ b/src/system_routines.f90
@@ -123,7 +123,7 @@ function getCWD()
     getCWD = c_f_string(getCWD_Cstring)
   else
     error stop 'invalid working directory'
-  endif
+  end if
 
 end function getCWD
 
@@ -145,7 +145,7 @@ function getHostName()
     getHostName = c_f_string(getHostName_Cstring)
   else
     getHostName = 'n/a (Error!)'
-  endif
+  end if
 
 end function getHostName