Merge branch 'pretty-print-init' into 'development'

Pretty print init

See merge request damask/DAMASK!454

src/CPFEM.f90

@@ -103,7 +103,7 @@ subroutine CPFEM_init
   class(tNode), pointer :: &
     debug_CPFEM
 
-  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
 
   allocate(CPFEM_cs(               6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
   allocate(CPFEM_dcsdE(          6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)

src/CPFEM2.f90

@@ -81,7 +81,7 @@ subroutine CPFEM_init
   integer(HID_T) :: fileHandle
 
 
-  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
 
 
   if (interface_restartInc > 0) then

src/DAMASK_Marc.f90

@@ -46,7 +46,7 @@ subroutine DAMASK_interface_init
   integer                 :: ierr
   character(len=pPathLen) :: wd
 
-  print'(/,a)', ' <<<+-  DAMASK_marc init -+>>>'
+  print'(/,1x,a)', '<<<+-  DAMASK_marc init -+>>>'
 
   print*, 'Roters et al., Computational Materials Science 158:420–478, 2019'
   print*, 'https://doi.org/10.1016/j.commatsci.2018.04.030'

src/DAMASK_interface.f90

@@ -70,7 +70,7 @@ subroutine DAMASK_interface_init
   external :: &
     quit
 
-  print'(/,a)', ' <<<+-  DAMASK_interface init  -+>>>'
+  print'(/,1x,a)', '<<<+-  DAMASK_interface init  -+>>>'
 
   if(worldrank == 0) open(OUTPUT_UNIT, encoding='UTF-8')                                            ! for special characters in output
 

src/HDF5_utilities.f90

@@ -109,7 +109,7 @@ subroutine HDF5_utilities_init
   integer(SIZE_T) :: typeSize
 
 
-  print'(/,a)', ' <<<+-  HDF5_Utilities init  -+>>>'
+  print'(/,1x,a)', '<<<+-  HDF5_Utilities init  -+>>>'
 
 
   call h5open_f(hdferr)

src/IO.f90

@@ -56,7 +56,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine IO_init
 
-  print'(/,a)', ' <<<+-  IO init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  IO init  -+>>>'; flush(IO_STDOUT)
 
   call selfTest
 

src/YAML_types.f90

@@ -187,7 +187,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine YAML_types_init
 
-  print'(/,a)', ' <<<+-  YAML_types init  -+>>>'
+  print'(/,1x,a)', '<<<+-  YAML_types init  -+>>>'
 
   call selfTest
 

src/base64.f90

@@ -27,7 +27,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine base64_init
 
-  print'(/,a)', ' <<<+-  base64 init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  base64 init  -+>>>'; flush(IO_STDOUT)
 
   call selfTest
 

src/config.f90

@@ -30,8 +30,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine config_init
 
-  print'(/,a)', ' <<<+-  config init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  config init  -+>>>'; flush(IO_STDOUT)
-
 
   call parse_material
   call parse_numerics
@@ -50,15 +49,15 @@ subroutine parse_material()
 
 
   inquire(file='material.yaml',exist=fileExists)
-  if(.not. fileExists) call IO_error(100,ext_msg='material.yaml')
+  if (.not. fileExists) call IO_error(100,ext_msg='material.yaml')
 
   if (worldrank == 0) then
-    print*, 'reading material.yaml'; flush(IO_STDOUT)
+    print'(/,1x,a)', 'reading material.yaml'; flush(IO_STDOUT)
     fileContent = IO_read('material.yaml')
     call results_openJobFile(parallel=.false.)
     call results_writeDataset_str(fileContent,'setup','material.yaml','main configuration')
     call results_closeJobFile
-  endif
+  end if
   call parallelization_bcast_str(fileContent)
 
   config_material => YAML_parse_str(fileContent)
@@ -81,19 +80,19 @@ subroutine parse_numerics()
   if (fileExists) then
 
     if (worldrank == 0) then
-      print*, 'reading numerics.yaml'; flush(IO_STDOUT)
+      print'(1x,a)', 'reading numerics.yaml'; flush(IO_STDOUT)
       fileContent = IO_read('numerics.yaml')
       if (len(fileContent) > 0) then
         call results_openJobFile(parallel=.false.)
         call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration')
         call results_closeJobFile
-      endif
+      end if
-    endif
+    end if
     call parallelization_bcast_str(fileContent)
 
     config_numerics => YAML_parse_str(fileContent)
 
-  endif
+  end if
 
 end subroutine parse_numerics
 
@@ -113,19 +112,19 @@ subroutine parse_debug()
   if (fileExists) then
 
     if (worldrank == 0) then
-      print*, 'reading debug.yaml'; flush(IO_STDOUT)
+      print'(1x,a)', 'reading debug.yaml'; flush(IO_STDOUT)
       fileContent = IO_read('debug.yaml')
       if (len(fileContent) > 0) then
         call results_openJobFile(parallel=.false.)
         call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration')
         call results_closeJobFile
-      endif
+      end if
-    endif
+    end if
     call parallelization_bcast_str(fileContent)
 
     config_debug => YAML_parse_str(fileContent)
 
-  endif
+  end if
 
 end subroutine parse_debug
 

src/discretization.f90

@@ -49,7 +49,7 @@ subroutine discretization_init(materialAt,&
   integer, optional,           intent(in) :: &
     sharedNodesBegin                                                                                !< index of first node shared among different processes (MPI)
 
-  print'(/,a)', ' <<<+-  discretization init  -+>>>'; flush(6)
+  print'(/,1x,a)', '<<<+-  discretization init  -+>>>'; flush(6)
 
   discretization_Nelems = size(materialAt,1)
   discretization_nIPs   = size(IPcoords0,2)/discretization_Nelems

src/element.f90

@@ -923,7 +923,7 @@ subroutine tElement_init(self,elemType)
 
   self%nIPneighbors = size(self%IPneighbor,1)
 
-  print'(/,a)', ' <<<+-  element_init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  element_init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'element type:      ',self%elemType
   print*, '  geom type:       ',self%geomType

src/grid/DAMASK_grid.f90

@@ -113,10 +113,10 @@ program DAMASK_grid
 ! init DAMASK (all modules)
 
   call CPFEM_initAll
-  print'(/,a)',   ' <<<+-  DAMASK_grid init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  DAMASK_grid init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
+  print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
-  print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
+  print'(  1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
 
 !-------------------------------------------------------------------------------------------------
@@ -227,12 +227,12 @@ program DAMASK_grid
     loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
 
     reportAndCheck: if (worldrank == 0) then
-      print'(/,a,i0)', ' load case: ', l
+      print'(/,1x,a,1x,i0)', 'load case:', l
-      print*, ' estimate_rate:', loadCases(l)%estimate_rate
+      print'(2x,a,1x,l1)', 'estimate_rate:', loadCases(l)%estimate_rate
       if (loadCases(l)%deformation%myType == 'F') then
-        print*, ' F:'
+        print'(2x,a)', 'F:'
       else
-        print*, ' '//loadCases(l)%deformation%myType//' / 1/s:'
+        print'(2x,a)', loadCases(l)%deformation%myType//' / 1/s:'
       endif
       do i = 1, 3; do j = 1, 3
         if (loadCases(l)%deformation%mask(i,j)) then
@@ -246,8 +246,8 @@ program DAMASK_grid
       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
 
-      if (loadCases(l)%stress%myType == 'P')     print*, ' P / MPa:'
+      if (loadCases(l)%stress%myType == 'P')     print'(2x,a)', 'P / MPa:'
-      if (loadCases(l)%stress%myType == 'dot_P') print*, ' dot_P / MPa/s:'
+      if (loadCases(l)%stress%myType == 'dot_P') print'(2x,a)', 'dot_P / MPa/s:'
 
       if (loadCases(l)%stress%myType /= '') then
         do i = 1, 3; do j = 1, 3
@@ -270,15 +270,15 @@ program DAMASK_grid
       if (loadCases(l)%f_restart < 1)  errorID = 839
 
       if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then
-        print'(a)', '  r: 1 (constant step width)'
+        print'(2x,a)', 'r: 1 (constant step width)'
       else
-        print'(a,f0.3)', '  r: ', loadCases(l)%r
+        print'(2x,a,1x,f0.3)', 'r:', loadCases(l)%r
       endif
-      print'(a,f0.3)',   '  t: ', loadCases(l)%t
+      print'(2x,a,1x,f0.3)',   't:', loadCases(l)%t
-      print'(a,i0)',     '  N: ', loadCases(l)%N
+      print'(2x,a,1x,i0)',     'N:', loadCases(l)%N
-      print'(a,i0)',     '  f_out: ', loadCases(l)%f_out
+      print'(2x,a,1x,i0)',     'f_out:', loadCases(l)%f_out
       if (loadCases(l)%f_restart < huge(0)) &
-        print'(a,i0)',   '  f_restart: ', loadCases(l)%f_restart
+        print'(2x,a,1x,i0)',   'f_restart:', loadCases(l)%f_restart
 
       if (errorID > 0) call IO_error(error_ID = errorID, el = l)
 
@@ -317,7 +317,7 @@ program DAMASK_grid
   endif
 
   writeUndeformed: if (interface_restartInc < 1) then
-    print'(/,a)', ' ... writing initial configuration to file ........................'
+    print'(/,1x,a)', '... writing initial configuration to file .................................'
     flush(IO_STDOUT)
     call CPFEM_results(0,0.0_pReal)
   endif writeUndeformed
@@ -353,8 +353,8 @@ program DAMASK_grid
 
 !--------------------------------------------------------------------------------------------------
 ! report begin of new step
-          print'(/,a)', ' ###########################################################################'
+          print'(/,1x,a)', '###########################################################################'
-          print'(1x,a,es12.5,6(a,i0))', &
+          print'(1x,a,1x,es12.5,6(a,i0))', &
                   'Time', t, &
                   's: Increment ', inc,'/',loadCases(l)%N,&
                   '-', stepFraction,'/',subStepFactor**cutBackLevel,&
@@ -379,7 +379,7 @@ program DAMASK_grid
               case(FIELD_DAMAGE_ID);  call grid_damage_spectral_forward(cutBack)
             end select
           enddo
-          if(.not. cutBack) call CPFEM_forward
+          if (.not. cutBack) call CPFEM_forward
 
 !--------------------------------------------------------------------------------------------------
 ! solve fields
@@ -425,7 +425,7 @@ program DAMASK_grid
             cutBackLevel = cutBackLevel + 1
             t = t - Delta_t
             Delta_t = Delta_t/real(subStepFactor,pReal)                                             ! cut timestep
-            print'(/,a)', ' cutting back '
+            print'(/,1x,a)', 'cutting back '
           else                                                                                      ! no more options to continue
             if (worldrank == 0) close(statUnit)
             call IO_error(950)
@@ -436,26 +436,26 @@ program DAMASK_grid
         cutBackLevel = max(0, cutBackLevel - 1)                                                     ! try half number of subincs next inc
 
         if (all(solres(:)%converged)) then
-          print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged'
+          print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged'
         else
-          print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged'
+          print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged'
         endif; flush(IO_STDOUT)
 
         call MPI_Allreduce(interface_SIGUSR1,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
         if (ierr /= 0) error stop 'MPI error'
         if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then
-          print'(1/,a)', ' ... writing results to file ......................................'
+          print'(/,1x,a)', '... writing results to file ...............................................'
           flush(IO_STDOUT)
           call CPFEM_results(totalIncsCounter,t)
         endif
-        if(signal) call interface_setSIGUSR1(.false.)
+        if (signal) call interface_setSIGUSR1(.false.)
         call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
         if (ierr /= 0) error stop 'MPI error'
         if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then
           call mechanical_restartWrite
           call CPFEM_restartWrite
         endif
-        if(signal) call interface_setSIGUSR2(.false.)
+        if (signal) call interface_setSIGUSR2(.false.)
         call MPI_Allreduce(interface_SIGTERM,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
         if (ierr /= 0) error stop 'MPI error'
         if (signal) exit loadCaseLooping
@@ -468,7 +468,7 @@ program DAMASK_grid
 
 !--------------------------------------------------------------------------------------------------
 ! report summary of whole calculation
-  print'(/,a)', ' ###########################################################################'
+  print'(/,1x,a)', '###########################################################################'
   if (worldrank == 0) close(statUnit)
 
   call quit(0)                                                                                      ! no complains ;)

src/grid/discretization_grid.f90

@@ -72,7 +72,7 @@ subroutine discretization_grid_init(restart)
     fileContent, fname
 
 
-  print'(/,a)', ' <<<+-  discretization_grid init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  discretization_grid init  -+>>>'; flush(IO_STDOUT)
 
 
   if (worldrank == 0) then
@@ -96,9 +96,9 @@ subroutine discretization_grid_init(restart)
   call MPI_Bcast(origin,3,MPI_DOUBLE,0,MPI_COMM_WORLD, ierr)
   if (ierr /= 0) error stop 'MPI error'
 
-  print'(/,a,3(i12  ))',  ' cells  a b c: ', grid
+  print'(/,1x,a,3(i12,1x))',    'cells  a b c: ', grid
-  print'(a,3(es12.5))',   ' size   x y z: ', geomSize
+  print  '(1x,a,3(es12.5,1x))', 'size   x y z: ', geomSize
-  print'(a,3(es12.5))',   ' origin x y z: ', origin
+  print  '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin
 
   if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
 

src/grid/grid_damage_spectral.f90

@@ -76,10 +76,10 @@ subroutine grid_damage_spectral_init()
   character(len=pStringLen) :: &
     snes_type
 
-  print'(/,a)', ' <<<+-  grid_spectral_damage init  -+>>>'
+  print'(/,1x,a)', '<<<+-  grid_spectral_damage init  -+>>>'
 
-  print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
+  print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
-  print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
+  print'(  1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
@@ -137,7 +137,7 @@ subroutine grid_damage_spectral_init()
     call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr)                                    ! variable bounds for variational inequalities like contact mechanics, damage etc.
     call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
     call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
-  endif
+  end if
 
 !--------------------------------------------------------------------------------------------------
 ! init fields
@@ -187,7 +187,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
   else
     solution%converged = .true.
     solution%iterationsNeeded = totalIter
-  endif
+  end if
   stagNorm = maxval(abs(phi_current - phi_stagInc))
   solnNorm = maxval(abs(phi_current))
   call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr)
@@ -201,14 +201,14 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
     call homogenization_set_phi(phi_current(i,j,k),ce)
-  enddo; enddo; enddo
+  end do; end do; end do
 
   call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
   call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
   if (solution%converged) &
-    print'(/,a)', ' ... nonlocal damage converged .....................................'
+    print'(/,1x,a)', '... nonlocal damage converged .....................................'
-  print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage      = ', phi_min, phi_max, stagNorm
+  print'(/,1x,a,f8.6,2x,f8.6,2x,e11.4)', 'Minimum|Maximum|Delta Damage      = ', phi_min, phi_max, stagNorm
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end function grid_damage_spectral_solution
@@ -238,11 +238,11 @@ subroutine grid_damage_spectral_forward(cutBack)
     do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
       ce = ce + 1
       call homogenization_set_phi(phi_current(i,j,k),ce)
-    enddo; enddo; enddo
+    end do; end do; end do
   else
     phi_lastInc = phi_current
     call updateReference
-  endif
+  end if
 
 end subroutine grid_damage_spectral_forward
 
@@ -277,7 +277,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
     vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k))
-  enddo; enddo; enddo
+  end do; end do; end do
   call utilities_FFTvectorForward
   call utilities_fourierVectorDivergence                                                            !< calculate damage divergence in fourier field
   call utilities_FFTscalarBackward
@@ -287,7 +287,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
     scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
                             + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
                             + mu_ref*phi_current(i,j,k)
-  enddo; enddo; enddo
+  end do; end do; end do
 
 !--------------------------------------------------------------------------------------------------
 ! convolution of damage field with green operator
@@ -320,7 +320,7 @@ subroutine updateReference()
   do ce = 1, product(grid(1:2))*grid3
     K_ref  = K_ref  + homogenization_K_phi(ce)
     mu_ref = mu_ref + homogenization_mu_phi(ce)
-  enddo
+  end do
 
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)

src/grid/grid_mech_FEM.f90

@@ -116,7 +116,7 @@ subroutine grid_mechanical_FEM_init
     num_grid, &
     debug_grid
 
-  print'(/,a)', ' <<<+-  grid_mechanical_FEM init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  grid_mechanical_FEM init  -+>>>'; flush(IO_STDOUT)
 
 !-------------------------------------------------------------------------------------------------
 ! debugging options
@@ -234,7 +234,7 @@ subroutine grid_mechanical_FEM_init
 !--------------------------------------------------------------------------------------------------
 ! init fields
   restartRead: if (interface_restartInc > 0) then
-    print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
+    print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
 
     fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
     groupHandle = HDF5_openGroup(fileHandle,'solver')
@@ -272,7 +272,7 @@ subroutine grid_mechanical_FEM_init
   CHKERRQ(ierr)
 
   restartRead2: if (interface_restartInc > 0) then
-    print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
+    print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
     call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
     call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
     if(ierr /=0) error stop 'MPI error'
@@ -442,7 +442,7 @@ subroutine grid_mechanical_FEM_restartWrite
   call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
   CHKERRQ(ierr)
 
-  print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
+  print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
 
   fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
@@ -506,12 +506,12 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
     reason = 0
   endif
 
-  print'(1/,a)', ' ... reporting .............................................................'
+  print'(/,1x,a)', '... reporting .............................................................'
-  print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
+  print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', &
           err_div/divTol,  ' (',err_div,' / m, tol = ',divTol,')'
-  print'(a,f12.2,a,es8.2,a,es9.2,a)',    ' error stress BC  = ', &
+  print'(1x,a,f12.2,a,es8.2,a,es9.2,a)',    'error stress BC  = ', &
           err_BC/BCTol,    ' (',err_BC, ' Pa,  tol = ',BCTol,')'
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end subroutine converged
@@ -547,10 +547,10 @@ subroutine formResidual(da_local,x_local, &
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax
-    if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
-      ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
+      'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
-      ' deformation gradient aim       =', transpose(F_aim)
+      'deformation gradient aim       =', transpose(F_aim)
     flush(IO_STDOUT)
   endif newIteration
 

src/grid/grid_mech_spectral_basic.f90

@@ -111,13 +111,13 @@ subroutine grid_mechanical_spectral_basic_init
     num_grid, &
     debug_grid
 
-  print'(/,a)', ' <<<+-  grid_mechanical_spectral_basic init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  grid_mechanical_spectral_basic init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
+  print'(/,1x,a)', 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
-  print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
 
-  print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
+  print'(  1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
-  print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
 
 !-------------------------------------------------------------------------------------------------
 ! debugging options
@@ -186,30 +186,30 @@ subroutine grid_mechanical_spectral_basic_init
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                   ! places pointer on PETSc data
 
   restartRead: if (interface_restartInc > 0) then
-    print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
+    print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
 
     fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
     groupHandle = HDF5_openGroup(fileHandle,'solver')
 
     call HDF5_read(P_aim,groupHandle,'P_aim',.false.)
     call MPI_Bcast(P_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(F_aim,groupHandle,'F_aim',.false.)
     call MPI_Bcast(F_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.)
     call MPI_Bcast(F_aim_lastInc,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.)
     call MPI_Bcast(F_aimDot,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(F,groupHandle,'F')
     call HDF5_read(F_lastInc,groupHandle,'F_lastInc')
 
   elseif (interface_restartInc == 0) then restartRead
     F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)                         ! initialize to identity
     F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
-  endif restartRead
+  end if restartRead
 
   homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
@@ -219,13 +219,13 @@ subroutine grid_mechanical_spectral_basic_init
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                ! deassociate pointer
 
   restartRead2: if (interface_restartInc > 0) then
-    print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
+    print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
     call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
     call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.)
     call MPI_Bcast(C_volAvgLastInc,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
 
     call HDF5_closeGroup(groupHandle)
     call HDF5_closeFile(fileHandle)
@@ -234,7 +234,7 @@ subroutine grid_mechanical_spectral_basic_init
                        MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr)
     call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr)
     call MPI_File_close(fileUnit,ierr)
-  endif restartRead2
+  end if restartRead2
 
   call utilities_updateGamma(C_minMaxAvg)
   call utilities_saveReferenceStiffness
@@ -263,7 +263,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
 !--------------------------------------------------------------------------------------------------
 ! update stiffness (and gamma operator)
   S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
-  if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
+  if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
 
 !--------------------------------------------------------------------------------------------------
 ! solve BVP
@@ -328,7 +328,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
     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
 
     Fdot = utilities_calculateRate(guess, &
                                    F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
@@ -336,7 +336,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
     F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
 
     homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
-  endif
+  end if
 
 !--------------------------------------------------------------------------------------------------
 ! update average and local deformation gradients
@@ -385,7 +385,7 @@ subroutine grid_mechanical_spectral_basic_restartWrite
 
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
 
-  print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
+  print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
 
   fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
@@ -406,7 +406,7 @@ subroutine grid_mechanical_spectral_basic_restartWrite
     call HDF5_write(C_minMaxAvg,groupHandle,'C_minMaxAvg',.false.)
     call HDF5_closeGroup(groupHandle)
     call HDF5_closeFile(fileHandle)
-  endif
+  end if
 
   if (num%update_gamma) call utilities_saveReferenceStiffness
 
@@ -443,14 +443,14 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     reason = -1
   else
     reason = 0
-  endif
+  end if
 
-  print'(1/,a)', ' ... reporting .............................................................'
+  print'(/,1x,a)', '... reporting .............................................................'
-  print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
+  print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', &
           err_div/divTol,  ' (',err_div,' / m, tol = ',divTol,')'
-  print'(a,f12.2,a,es8.2,a,es9.2,a)',    ' error stress BC  = ', &
+  print'(1x,a,f12.2,a,es8.2,a,es9.2,a)',    'error stress BC  = ', &
           err_BC/BCTol,    ' (',err_BC, ' Pa,  tol = ',BCTol,')'
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end subroutine converged
@@ -485,12 +485,12 @@ subroutine formResidual(in, F, &
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
-    if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
-      ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
+      'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
-      ' deformation gradient aim       =', transpose(F_aim)
+      'deformation gradient aim       =', transpose(F_aim)
     flush(IO_STDOUT)
-  endif newIteration
+  end if newIteration
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate constitutive response

src/grid/grid_mech_spectral_polarisation.f90

@@ -124,10 +124,10 @@ subroutine grid_mechanical_spectral_polarisation_init
     num_grid, &
     debug_grid
 
-  print'(/,a)', ' <<<+-  grid_mechanical_spectral_polarization init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  grid_mechanical_spectral_polarization init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
+  print'(/,1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
-  print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
 
 !-------------------------------------------------------------------------------------------------
 ! debugging options
@@ -208,23 +208,23 @@ subroutine grid_mechanical_spectral_polarisation_init
   F_tau => FandF_tau(9:17,:,:,:)
 
   restartRead: if (interface_restartInc > 0) then
-    print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
+    print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
 
     fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
     groupHandle = HDF5_openGroup(fileHandle,'solver')
 
     call HDF5_read(P_aim,groupHandle,'P_aim',.false.)
     call MPI_Bcast(P_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(F_aim,groupHandle,'F_aim',.false.)
     call MPI_Bcast(F_aim,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(F_aim_lastInc,groupHandle,'F_aim_lastInc',.false.)
     call MPI_Bcast(F_aim_lastInc,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(F_aimDot,groupHandle,'F_aimDot',.false.)
     call MPI_Bcast(F_aimDot,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(F,groupHandle,'F')
     call HDF5_read(F_lastInc,groupHandle,'F_lastInc')
     call HDF5_read(F_tau,groupHandle,'F_tau')
@@ -235,7 +235,7 @@ subroutine grid_mechanical_spectral_polarisation_init
     F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
     F_tau = 2.0_pReal*F
     F_tau_lastInc = 2.0_pReal*F_lastInc
-  endif restartRead
+  end if restartRead
 
   homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
@@ -245,13 +245,13 @@ subroutine grid_mechanical_spectral_polarisation_init
   call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)                        ! deassociate pointer
 
   restartRead2: if (interface_restartInc > 0) then
-    print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
+    print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
     call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
     call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
     call HDF5_read(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.)
     call MPI_Bcast(C_volAvgLastInc,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-    if(ierr /=0) error stop 'MPI error'
+    if (ierr /=0) error stop 'MPI error'
 
     call HDF5_closeGroup(groupHandle)
     call HDF5_closeFile(fileHandle)
@@ -260,7 +260,7 @@ subroutine grid_mechanical_spectral_polarisation_init
                        MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr)
     call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr)
     call MPI_File_close(fileUnit,ierr)
-  endif restartRead2
+  end if restartRead2
 
   call utilities_updateGamma(C_minMaxAvg)
   call utilities_saveReferenceStiffness
@@ -291,11 +291,11 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
 !--------------------------------------------------------------------------------------------------
 ! update stiffness (and gamma operator)
   S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
-  if(num%update_gamma) then
+  if (num%update_gamma) then
     call utilities_updateGamma(C_minMaxAvg)
     C_scale = C_minMaxAvg
     S_scale = math_invSym3333(C_minMaxAvg)
-  endif
+  end if
 
 !--------------------------------------------------------------------------------------------------
 ! solve BVP
@@ -364,7 +364,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
     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
 
     Fdot     = utilities_calculateRate(guess, &
                                        F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
@@ -376,14 +376,14 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
     F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
 
     homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
-  endif
+  end if
 
 !--------------------------------------------------------------------------------------------------
 ! update average and local deformation gradients
   F_aim = F_aim_lastInc + F_aimDot * Delta_t
-  if(stress_BC%myType=='P')     P_aim = P_aim &
+  if (stress_BC%myType=='P')     P_aim = P_aim &
                                       + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t
-  if(stress_BC%myType=='dot_P') P_aim = P_aim &
+  if (stress_BC%myType=='dot_P') P_aim = P_aim &
                                       + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
 
   F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, &                                      ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
@@ -399,8 +399,8 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
                   + math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, &
                     math_mul3333xx33(C_scale,matmul(transpose(F_lambda33),F_lambda33)-math_I3)))
        F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
-    enddo; enddo; enddo
+    end do; end do; end do
-  endif
+  end if
 
   call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
 
@@ -441,7 +441,7 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite
   F     => FandF_tau(0: 8,:,:,:)
   F_tau => FandF_tau(9:17,:,:,:)
 
-  print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
+  print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
 
   fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
@@ -463,9 +463,9 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite
     call HDF5_write(C_volAvgLastInc,groupHandle,'C_volAvgLastInc',.false.)
     call HDF5_closeGroup(groupHandle)
     call HDF5_closeFile(fileHandle)
-  endif
+  end if
 
-  if(num%update_gamma) call utilities_saveReferenceStiffness
+  if (num%update_gamma) call utilities_saveReferenceStiffness
 
   call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
 
@@ -502,16 +502,16 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     reason = -1
   else
     reason = 0
-  endif
+  end if
 
-  print'(1/,a)', ' ... reporting .............................................................'
+  print'(/,1x,a)', '... reporting .............................................................'
-  print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
+  print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', &
             err_div/divTol,  ' (',err_div, ' / m, tol = ',divTol,')'
-  print  '(a,f12.2,a,es8.2,a,es9.2,a)', ' error curl       = ', &
+  print  '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error curl       = ', &
             err_curl/curlTol,' (',err_curl,' -,   tol = ',curlTol,')'
-  print  '(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC  = ', &
+  print  '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC  = ', &
             err_BC/BCTol,    ' (',err_BC,  ' Pa,  tol = ',BCTol,')'
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end subroutine converged
@@ -565,12 +565,12 @@ subroutine formResidual(in, FandF_tau, &
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
-    if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
-      ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
+      'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
-      ' deformation gradient aim       =', transpose(F_aim)
+      'deformation gradient aim       =', transpose(F_aim)
     flush(IO_STDOUT)
-  endif newIteration
+  end if newIteration
 
 !--------------------------------------------------------------------------------------------------
 !
@@ -580,7 +580,7 @@ subroutine formResidual(in, FandF_tau, &
       num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
       num%alpha*matmul(F(1:3,1:3,i,j,k), &
                          math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
-  enddo; enddo; enddo
+  end do; end do; end do
 
 !--------------------------------------------------------------------------------------------------
 ! doing convolution in Fourier space
@@ -621,7 +621,7 @@ subroutine formResidual(in, FandF_tau, &
                        residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
                        math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
                        + residual_F_tau(1:3,1:3,i,j,k)
-  enddo; enddo; enddo
+  end do; end do; end do
 
 !--------------------------------------------------------------------------------------------------
 ! calculating curl

src/grid/grid_thermal_spectral.f90

@@ -75,10 +75,10 @@ subroutine grid_thermal_spectral_init(T_0)
   class(tNode), pointer :: &
     num_grid
 
-  print'(/,a)', ' <<<+-  grid_thermal_spectral init  -+>>>'
+  print'(/,1x,a)', '<<<+-  grid_thermal_spectral init  -+>>>'
 
-  print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
+  print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
-  print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
+  print'(  1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
@@ -141,7 +141,7 @@ subroutine grid_thermal_spectral_init(T_0)
     T_lastInc(i,j,k) = T_current(i,j,k)
     T_stagInc(i,j,k) = T_current(i,j,k)
     call homogenization_thermal_setField(T_0,0.0_pReal,ce)
-  enddo; enddo; enddo
+  end do; end do; end do
 
   call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr)
   T_PETSc(xstart:xend,ystart:yend,zstart:zend) = T_current
@@ -182,7 +182,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
   else
     solution%converged = .true.
     solution%iterationsNeeded = totalIter
-  endif
+  end if
   stagNorm = maxval(abs(T_current - T_stagInc))
   solnNorm = maxval(abs(T_current))
   call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD,ierr)
@@ -196,14 +196,14 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
     call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
-  enddo; enddo; enddo
+  end do; end do; end do
 
   call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
   call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr)
   if (solution%converged) &
-    print'(/,a)', ' ... thermal conduction converged ..................................'
+    print'(/,1x,a)', '... thermal conduction converged ..................................'
-  print'(/,a,f8.4,2x,f8.4,2x,f8.4)', ' Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm
+  print'(/,1x,a,f8.4,2x,f8.4,2x,f8.4)', 'Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end function grid_thermal_spectral_solution
@@ -234,11 +234,11 @@ subroutine grid_thermal_spectral_forward(cutBack)
     do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
       ce = ce + 1
       call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
-    enddo; enddo; enddo
+    end do; end do; end do
   else
     T_lastInc = T_current
     call updateReference
-  endif
+  end if
 
 end subroutine grid_thermal_spectral_forward
 
@@ -272,7 +272,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
     vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k))
-  enddo; enddo; enddo
+  end do; end do; end do
   call utilities_FFTvectorForward
   call utilities_fourierVectorDivergence                                                            !< calculate temperature divergence in fourier field
   call utilities_FFTscalarBackward
@@ -282,7 +282,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
     scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
                             + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
                             + mu_ref*T_current(i,j,k)
-  enddo; enddo; enddo
+  end do; end do; end do
 
 !--------------------------------------------------------------------------------------------------
 ! convolution of temperature field with green operator
@@ -310,7 +310,7 @@ subroutine updateReference()
   do ce = 1, product(grid(1:2))*grid3
     K_ref  = K_ref  + homogenization_K_T(ce)
     mu_ref = mu_ref + homogenization_mu_T(ce)
-  enddo
+  end do
 
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)

src/grid/spectral_utilities.f90

@@ -177,19 +177,19 @@ subroutine spectral_utilities_init
     num_grid, &
     debug_grid                                                                                      ! pointer to grid  debug options
 
-  print'(/,a)', ' <<<+-  spectral_utilities init  -+>>>'
+  print'(/,1x,a)', '<<<+-  spectral_utilities init  -+>>>'
 
-  print*, 'M. Diehl, Diploma Thesis TU München, 2010'
+  print'(/,1x,a)', 'M. Diehl, Diploma Thesis TU München, 2010'
-  print*, 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL
 
-  print*, 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
+  print'(  1x,a)', 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
-  print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
 
-  print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
+  print'(  1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
-  print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL
 
-  print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
+  print'(  1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
-  print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
+  print'(  1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
 !--------------------------------------------------------------------------------------------------
 ! set debugging parameters
@@ -200,15 +200,15 @@ subroutine spectral_utilities_init
   debugRotation   =  debug_grid%contains('rotation')
   debugPETSc      =  debug_grid%contains('PETSc')
 
-  if(debugPETSc) print'(3(/,a),/)', &
+  if (debugPETSc) print'(3(/,1x,a),/)', &
-                 ' Initializing PETSc with debug options: ', &
+                 'Initializing PETSc with debug options: ', &
                  trim(PETScDebug), &
-                 ' add more using the "PETSc_options" keyword in numerics.yaml'
+                 'add more using the "PETSc_options" keyword in numerics.yaml'
   flush(IO_STDOUT)
 
   call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
   CHKERRQ(ierr)
-  if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
+  if (debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
   CHKERRQ(ierr)
   call PetscOptionsInsertString(PETSC_NULL_OPTIONS,&
                                 num_grid%get_asString('PETSc_options',defaultVal=''),ierr)
@@ -271,7 +271,7 @@ subroutine spectral_utilities_init
   if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match'
   call fftw_set_timelimit(num_grid%get_asFloat('fftw_timelimit',defaultVal=-1.0_pReal))
 
-  print*, 'FFTW initialized'; flush(IO_STDOUT)
+  print'(/,1x,a)', 'FFTW initialized'; flush(IO_STDOUT)
 
 !--------------------------------------------------------------------------------------------------
 ! MPI allocation
@@ -342,10 +342,10 @@ subroutine spectral_utilities_init
 ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
   do k = grid3Offset+1, grid3Offset+grid3
     k_s(3) = k - 1
-    if(k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3)                                                 ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
+    if (k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3)                                                 ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
       do j = 1, grid(2)
         k_s(2) = j - 1
-        if(j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2)                                             ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
+        if (j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2)                                             ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
           do i = 1, grid1Red
             k_s(1) = i - 1                                                                          ! symmetry, junst running from 0,1,...,N/2,N/2+1
             xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s)
@@ -357,7 +357,7 @@ subroutine spectral_utilities_init
             endwhere
   enddo; enddo; enddo
 
-  if(num%memory_efficient) then                                                                     ! allocate just single fourth order tensor
+  if (num%memory_efficient) then                                                                     ! allocate just single fourth order tensor
     allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal))
   else                                                                                              ! precalculation of gamma_hat field
     allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
@@ -384,7 +384,7 @@ subroutine utilities_updateGamma(C)
 
   C_ref = C
 
-  if(.not. num%memory_efficient) then
+  if (.not. num%memory_efficient) then
     gamma_hat =  cmplx(0.0_pReal,0.0_pReal,pReal)                                                   ! for the singular point and any non invertible A
     do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red
       if (any([i,j,k] /= 1)) then                                                                   ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
@@ -497,12 +497,12 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
   logical :: err
 
 
-  print'(/,a)', ' ... doing gamma convolution ...............................................'
+  print'(/,1x,a)', '... doing gamma convolution ...............................................'
   flush(IO_STDOUT)
 
 !--------------------------------------------------------------------------------------------------
 ! do the actual spectral method calculation (mechanical equilibrium)
-  memoryEfficient: if(num%memory_efficient) then
+  memoryEfficient: if (num%memory_efficient) then
     do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
       if (any([i,j,k+grid3Offset] /= 1)) then                                                       ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
         forall(l = 1:3, m = 1:3) &
@@ -567,7 +567,7 @@ real(pReal) function utilities_divergenceRMS()
   integer :: i, j, k, ierr
   complex(pReal), dimension(3)   :: rescaledGeom
 
-  print'(/,a)', ' ... calculating divergence ................................................'
+  print'(/,1x,a)', '... calculating divergence ................................................'
   flush(IO_STDOUT)
 
   rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
@@ -593,9 +593,9 @@ real(pReal) function utilities_divergenceRMS()
                + sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
                                   conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal)
   enddo; enddo
-  if(grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal                    ! counted twice in case of grid(1) == 1
+  if (grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal                    ! counted twice in case of grid(1) == 1
   call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
-  if(ierr /=0) error stop 'MPI error'
+  if (ierr /=0) error stop 'MPI error'
   utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt                                     ! RMS in real space calculated with Parsevals theorem from Fourier space
 
 end function utilities_divergenceRMS
@@ -610,7 +610,7 @@ real(pReal) function utilities_curlRMS()
   complex(pReal), dimension(3,3) :: curl_fourier
   complex(pReal), dimension(3)   :: rescaledGeom
 
-  print'(/,a)', ' ... calculating curl ......................................................'
+  print'(/,1x,a)', '... calculating curl ......................................................'
   flush(IO_STDOUT)
 
   rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
@@ -655,9 +655,9 @@ real(pReal) function utilities_curlRMS()
   enddo; enddo
 
   call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
-  if(ierr /=0) error stop 'MPI error'
+  if (ierr /=0) error stop 'MPI error'
   utilities_curlRMS = sqrt(utilities_curlRMS) * wgt
-  if(grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal                                ! counted twice in case of grid(1) == 1
+  if (grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal                                ! counted twice in case of grid(1) == 1
 
 end function utilities_curlRMS
 
@@ -686,13 +686,13 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
 
   mask_stressVector = .not. reshape(transpose(mask_stress), [9])
   size_reduced = count(mask_stressVector)
-  if(size_reduced > 0) then
+  if (size_reduced > 0) then
     temp99_real = math_3333to99(rot_BC%rotate(C))
 
-    if(debugGeneral) then
+    if (debugGeneral) then
-      print'(/,a)', ' ... updating masked compliance ............................................'
+      print'(/,1x,a)', '... updating masked compliance ............................................'
-      print'(/,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', &
+      print'(/,1x,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', &
-        ' Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal
+        'Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal
       flush(IO_STDOUT)
     endif
 
@@ -711,10 +711,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
     errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pReal))
     if (debugGeneral .or. errmatinv) then
       write(formatString, '(i2)') size_reduced
-      formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
+      formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
-      print trim(formatString), ' C * S (load) ', transpose(matmul(c_reduced,s_reduced))
+      print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced))
-      print trim(formatString), ' S (load) ', transpose(s_reduced)
+      print trim(formatString), 'S (load) ', transpose(s_reduced)
-      if(errmatinv) error stop 'matrix inversion error'
+      if (errmatinv) error stop 'matrix inversion error'
     endif
     temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9])
   else
@@ -723,9 +723,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
 
   utilities_maskedCompliance = math_99to3333(temp99_Real)
 
-  if(debugGeneral) then
+  if (debugGeneral) then
-    print'(/,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', &
+    print'(/,1x,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', &
-      ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
+      'Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
     flush(IO_STDOUT)
   endif
 
@@ -810,7 +810,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   real(pReal)                     :: dPdF_norm_max, dPdF_norm_min
   real(pReal), dimension(2) :: valueAndRank                                                         !< pair of min/max norm of dPdF to synchronize min/max of dPdF
 
-  print'(/,a)', ' ... evaluating constitutive response ......................................'
+  print'(/,1x,a)', '... evaluating constitutive response ......................................'
   flush(IO_STDOUT)
 
   homogenization_F  = reshape(F,[3,3,product(grid(1:2))*grid3])                                     ! set materialpoint target F to estimated field
@@ -824,11 +824,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
   P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
-  if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
+  if (debugRotation) print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
-    ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal
+    'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal
-  if(present(rotation_BC)) P_av = rotation_BC%rotate(P_av)
+  if (present(rotation_BC)) P_av = rotation_BC%rotate(P_av)
-  print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
+  print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
-    ' Piola--Kirchhoff stress       / MPa =', transpose(P_av)*1.e-6_pReal
+    'Piola--Kirchhoff stress       / MPa =', transpose(P_av)*1.e-6_pReal
   flush(IO_STDOUT)
 
   dPdF_max = 0.0_pReal
@@ -1017,7 +1017,7 @@ subroutine utilities_updateCoords(F)
   call utilities_FFTtensorForward()
 
   do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
-    if(any([i,j,k+grid3Offset] /= 1)) then
+    if (any([i,j,k+grid3Offset] /= 1)) then
       vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) &
                                      / sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal)
     else
@@ -1031,7 +1031,7 @@ subroutine utilities_updateCoords(F)
  ! average F
   if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
   call MPI_Bcast(Favg,9,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
-  if(ierr /=0) error stop 'MPI error'
+  if (ierr /=0) error stop 'MPI error'
 
  !--------------------------------------------------------------------------------------------------
  ! pad cell center fluctuations along z-direction (needed when running MPI simulation)
@@ -1042,22 +1042,22 @@ subroutine utilities_updateCoords(F)
 
   ! send bottom layer to process below
   call MPI_Isend(IPfluct_padded(:,:,:,2),      c,MPI_DOUBLE,rank_b,0,MPI_COMM_WORLD,request(1),ierr)
-  if(ierr /=0) error stop 'MPI error'
+  if (ierr /=0) error stop 'MPI error'
   call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0,MPI_COMM_WORLD,request(2),ierr)
-  if(ierr /=0) error stop 'MPI error'
+  if (ierr /=0) error stop 'MPI error'
 
   ! send top layer to process above
   call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,1,MPI_COMM_WORLD,request(3),ierr)
-  if(ierr /=0) error stop 'MPI error'
+  if (ierr /=0) error stop 'MPI error'
   call MPI_Irecv(IPfluct_padded(:,:,:,1),      c,MPI_DOUBLE,rank_b,1,MPI_COMM_WORLD,request(4),ierr)
-  if(ierr /=0) error stop 'MPI error'
+  if (ierr /=0) error stop 'MPI error'
 
   call MPI_Waitall(4,request,status,ierr)
-  if(ierr /=0) error stop 'MPI error'
+  if (ierr /=0) error stop 'MPI error'
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
   ! ToDo
 #else
-  if(any(status(MPI_ERROR,:) /= 0)) error stop 'MPI error'
+  if (any(status(MPI_ERROR,:) /= 0)) error stop 'MPI error'
 #endif
 
  !--------------------------------------------------------------------------------------------------
@@ -1094,10 +1094,10 @@ subroutine utilities_saveReferenceStiffness
     fileUnit,ierr
 
   if (worldrank == 0) then
-    print'(a)', ' writing reference stiffness data required for restart to file'; flush(IO_STDOUT)
+    print'(/,1x,a)', '... writing reference stiffness data required for restart to file .........'; flush(IO_STDOUT)
     open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',&
          status='replace',access='stream',action='write',iostat=ierr)
-    if(ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref')
+    if (ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref')
     write(fileUnit) C_ref
     close(fileUnit)
   endif

src/homogenization.f90

@@ -199,7 +199,7 @@ subroutine homogenization_init()
     num_homog, &
     num_homogGeneric
 
-  print'(/,a)', ' <<<+-  homogenization init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  homogenization init  -+>>>'; flush(IO_STDOUT)
 
 
   allocate(homogState      (size(material_name_homogenization)))

src/homogenization_damage.f90

@@ -41,7 +41,7 @@ module subroutine damage_init()
   integer :: ho,Nmembers
 
 
-  print'(/,a)', ' <<<+-  homogenization:damage init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:damage init  -+>>>'
 
 
   configHomogenizations => config_material%get('homogenization')

src/homogenization_damage_pass.f90

@@ -8,7 +8,7 @@ contains
 
 module subroutine pass_init()
   
-  print'(/,a)', ' <<<+-  homogenization:damage:pass init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:damage:pass init  -+>>>'
 
 end subroutine pass_init
 

src/homogenization_mechanical.f90

@@ -68,7 +68,7 @@ module subroutine mechanical_init(num_homog)
   class(tNode), pointer :: &
     num_homogMech
 
-  print'(/,a)', ' <<<+-  homogenization:mechanical init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:mechanical init  -+>>>'
 
   call material_parseHomogenization2()
 
@@ -114,7 +114,7 @@ module subroutine mechanical_partition(subF,ce)
 
   do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
     call phase_set_F(Fs(1:3,1:3,co),co,ce)
-  enddo
+  end do
 
 
 end subroutine mechanical_partition
@@ -138,7 +138,7 @@ module subroutine mechanical_homogenize(Delta_t,ce)
                                             + phase_P(co,ce)
     homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
                                             + phase_mechanical_dPdF(Delta_t,co,ce)
-  enddo
+  end do
 
   homogenization_P(1:3,1:3,ce)            = homogenization_P(1:3,1:3,ce) &
                                           / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
@@ -173,11 +173,11 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
         dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce)
         Fs(:,:,co)        = phase_F(co,ce)
         Ps(:,:,co)        = phase_P(co,ce)
-      enddo
+      end do
       doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce)
   else
     doneAndHappy = .true.
-  endif
+  end if
 
 end function mechanical_updateState
 
@@ -241,7 +241,7 @@ subroutine material_parseHomogenization2()
       case default
         call IO_error(500,ext_msg=homogMech%get_asString('type'))
     end select
-  enddo
+  end do
 
 end subroutine material_parseHomogenization2
 

src/homogenization_mechanical_RGC.f90

@@ -87,16 +87,16 @@ module subroutine RGC_init(num_homogMech)
     homog, &
     homogMech
 
-  print'(/,a)', ' <<<+-  homogenization:mechanical:RGC init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:mechanical:RGC init  -+>>>'
 
-  print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
+  print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
   flush(IO_STDOUT)
 
-  print*, 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
+  print'(/,1x,a)', 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
-  print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
 
-  print*, 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
+  print'(/,1x,a)', 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
-  print*, 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL
 
 
   material_homogenization => config_material%get('homogenization')
@@ -186,7 +186,7 @@ module subroutine RGC_init(num_homogMech)
 
     end associate
 
-  enddo
+  end do
 
 end subroutine RGC_init
 
@@ -222,9 +222,9 @@ module subroutine RGC_partitionDeformation(F,avgF,ce)
       nVect = interfaceNormal(intFace,ho,en)
       forall (i=1:3,j=1:3) &
         F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                           ! calculating deformation relaxations due to interface relaxation
-    enddo
+    end do
     F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF                                                    ! resulting relaxed deformation gradient
-  enddo
+  end do
 
   end associate
 
@@ -257,10 +257,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
   real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
   real(pReal), dimension(:), allocatable   :: resid,relax,p_relax,p_resid,drelax
 
-  zeroTimeStep: if(dEq0(dt)) then
+  zeroTimeStep: if (dEq0(dt)) then
     doneAndHappy = .true.                                                                   ! pretend everything is fine and return
     return
-  endif zeroTimeStep
+  end if zeroTimeStep
 
   ho  = material_homogenizationID(ce)
   en = material_homogenizationEntry(ce)
@@ -319,10 +319,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
         tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) &        ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface
                                       + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j)
         resid(i+3*(iNum-1)) = tract(iNum,i)                                                         ! translate the local residual into global 1-dimensional residual array
-      enddo
+      end do
-    enddo
+    end do
 
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! convergence check for stress residual
@@ -347,7 +347,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
   elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then                           ! try to restart when residual blows up exceeding maximum bound
     doneAndHappy = [.true.,.false.]                                                         ! with direct cut-back
    return
-  endif
+  end if
 
 !---------------------------------------------------------------------------------------------------
 ! construct the global Jacobian matrix for updating the global relaxation vector array when
@@ -373,11 +373,11 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
         do i=1,3; do j=1,3; do k=1,3; do l=1,3
           smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
                                              + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l)
-        enddo;enddo;enddo;enddo
+        end do;end do;end do;end do
 ! projecting the material tangent dPdF into the interface
 ! to obtain the Jacobian matrix contribution of dPdF
-      endif
+      end if
-    enddo
+    end do
 
 !--------------------------------------------------------------------------------------------------
 ! identify the right/up/front grain (+|P)
@@ -394,10 +394,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
         do i=1,3; do j=1,3; do k=1,3; do l=1,3
           smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
                                              + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l)
-        enddo;enddo;enddo;enddo
+        end do;end do;end do;end do
-      endif
+      end if
-    enddo
+    end do
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
@@ -443,10 +443,10 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
                                                       + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) &
                                                       + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) &
                                                       + (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
-      enddo; enddo
+      end do; end do
-    enddo
+    end do
     pmatrix(:,ipert) = p_resid/num%pPert
-  enddo
+  end do
 
 
 !--------------------------------------------------------------------------------------------------
@@ -455,7 +455,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
   do i=1,3*nIntFaceTot
     rmatrix(i,i) = num%viscModus*num%viscPower/(num%refRelaxRate*dt)* &                             ! tangent due to numerical viscosity traction appears
                    (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_pReal)              ! only in the main diagonal term
-  enddo
+  end do
 
 
 !--------------------------------------------------------------------------------------------------
@@ -472,7 +472,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
   drelax = 0.0_pReal
   do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot
     drelax(i) = drelax(i) - jnverse(i,j)*resid(j)                                                   ! Calculate the correction for the state variable
-  enddo; enddo
+  end do; end do
   stt%relaxationVector(:,en) = relax + drelax                                                       ! Updateing the state variable for the next iteration
   if (any(abs(drelax) > num%maxdRelax)) then                                                        ! Forcing cutback when the incremental change of relaxation vector becomes too large
     doneAndHappy = [.true.,.false.]
@@ -481,7 +481,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
     print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
     flush(IO_STDOUT)
     !$OMP END CRITICAL (write2out)
-  endif
+  end if
 
   end associate
 
@@ -545,9 +545,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
        do i = 1,3; do j = 1,3
          do k = 1,3; do l = 1,3
            nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l)                        ! compute the interface mismatch tensor from the jump of deformation gradient
-         enddo; enddo
+         end do; end do
          nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal                                                 ! compute the norm of the mismatch tensor
-       enddo; enddo
+       end do; end do
        nDefNorm = max(nDefToler,sqrt(nDefNorm))                                                     ! approximation to zero mismatch if mismatch is zero (singularity)
        nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm                       ! total amount of mismatch experienced by the grain (at all six interfaces)
 
@@ -560,11 +560,11 @@ 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)
-       enddo; enddo;enddo; enddo
+       end do; end do;enddo; end do
-     enddo interfaceLoop
+     end do interfaceLoop
 
 
-   enddo grainLoop
+   end do grainLoop
 
    end associate
 
@@ -594,7 +594,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
       gVol(i) = math_det33(fDef(1:3,1:3,i))                                                         ! compute the volume of individual grains
       vDiscrep     = vDiscrep - gVol(i)/real(nGrain,pReal)                                          ! calculate the difference/dicrepancy between
                                                                                                     ! the volume of the cluster and the the total volume of grains
-    enddo
+    end do
 
     !----------------------------------------------------------------------------------------------
     ! calculate the stress and penalty due to volume discrepancy
@@ -603,7 +603,7 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
       vPen(:,:,i) = -1.0_pReal/real(nGrain,pReal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
                          sign((abs(vDiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vDiscrep)* &
                          gVol(i)*transpose(math_inv33(fDef(:,:,i)))
-    enddo
+    end do
 
   end subroutine volumePenalty
 
@@ -633,9 +633,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
       nVect = interfaceNormal([iBase,1,1,1],ho,en)
       do i = 1,3; do j = 1,3
         surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j)           ! compute the component of (the inverse of) the stretch in the direction of the normal
-      enddo; enddo
+      end do; end do
       surfaceCorrection(iBase) = sqrt(surfaceCorrection(iBase))*detF                                ! get the surface correction factor (area contraction/enlargement)
-    enddo
+    end do
 
   end function surfaceCorrection
 
@@ -690,9 +690,9 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
         nVect   = interfaceNormal(intFace,ho,en)
         forall (i=1:3,j=1:3) &
           F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                         ! effective relaxations
-      enddo
+      end do
       F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF                                                  ! relaxed deformation gradient
-    enddo
+    end do
 
     end associate
 
@@ -727,7 +727,7 @@ module subroutine RGC_results(ho,group)
         call results_writeDataset(dst%relaxationrate_avg,group,trim(prm%output(o)), &
                                   'average relaxation rate','m/s')
     end select
-  enddo outputsLoop
+  end do outputsLoop
   end associate
 
 end subroutine RGC_results
@@ -756,7 +756,7 @@ pure function relaxationVector(intFace,ho,en)
     relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en)
   else
     relaxationVector = 0.0_pReal
-  endif
+  end if
 
   end associate
 
@@ -855,7 +855,7 @@ integer pure function interface4to1(iFace4D, nGDim)
       else
         interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) &
                       + nGDim(2)*nGDim(3)*(iFace4D(2)-1)
-      endif
+      end if
 
     case(2)
       if ((iFace4D(3) == 0) .or. (iFace4D(3) == nGDim(2))) then
@@ -864,7 +864,7 @@ integer pure function interface4to1(iFace4D, nGDim)
         interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) &
                       + nGDim(3)*nGDim(1)*(iFace4D(3)-1) &
                       + (nGDim(1)-1)*nGDim(2)*nGDim(3)                                              ! total # of interfaces normal || e1
-      endif
+      end if
 
     case(3)
       if ((iFace4D(4) == 0) .or. (iFace4D(4) == nGDim(3))) then
@@ -874,7 +874,7 @@ integer pure function interface4to1(iFace4D, nGDim)
                       + nGDim(1)*nGDim(2)*(iFace4D(4)-1) &
                       + (nGDim(1)-1)*nGDim(2)*nGDim(3) &                                            ! total # of interfaces normal || e1
                       + nGDim(1)*(nGDim(2)-1)*nGDim(3)                                              ! total # of interfaces normal || e2
-      endif
+      end if
 
     case default
       interface4to1 = -1
@@ -918,7 +918,7 @@ pure function interface1to4(iFace1D, nGDim)
     interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1
     interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)),nGDim(2))+1
     interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/real(nGDim(1),pReal)/real(nGDim(2),pReal))+1
-  endif
+  end if
 
 end function interface1to4
 

src/homogenization_mechanical_isostrain.f90

@@ -17,9 +17,9 @@ module subroutine isostrain_init
     ho, &
     Nmembers
 
-  print'(/,a)', ' <<<+-  homogenization:mechanical:isostrain init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:mechanical:isostrain init  -+>>>'
 
-  print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
+  print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
   flush(IO_STDOUT)
 
   do ho = 1, size(homogenization_type)
@@ -30,7 +30,7 @@ module subroutine isostrain_init
     allocate(homogState(ho)%state0(0,Nmembers))
     allocate(homogState(ho)%state (0,Nmembers))
 
-  enddo
+  end do
 
 end subroutine isostrain_init
 

src/homogenization_mechanical_pass.f90

@@ -17,15 +17,15 @@ module subroutine pass_init
     ho, &
     Nmembers
 
-  print'(/,a)', ' <<<+-  homogenization:mechanical:pass init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:mechanical:pass init  -+>>>'
 
-  print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_NONE_ID)
+  print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_NONE_ID)
   flush(IO_STDOUT)
 
   do ho = 1, size(homogenization_type)
-    if(homogenization_type(ho) /= HOMOGENIZATION_NONE_ID) cycle
+    if (homogenization_type(ho) /= HOMOGENIZATION_NONE_ID) cycle
 
-    if(homogenization_Nconstituents(ho) /= 1) &
+    if (homogenization_Nconstituents(ho) /= 1) &
       call IO_error(211,ext_msg='N_constituents (pass)')
 
     Nmembers = count(material_homogenizationID == ho)
@@ -33,7 +33,7 @@ module subroutine pass_init
     allocate(homogState(ho)%state0(0,Nmembers))
     allocate(homogState(ho)%state (0,Nmembers))
 
-  enddo
+  end do
 
 end subroutine pass_init
 

src/homogenization_thermal.f90

@@ -44,7 +44,7 @@ module subroutine thermal_init()
   integer :: ho
 
 
-  print'(/,a)', ' <<<+-  homogenization:thermal init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:thermal init  -+>>>'
 
 
   configHomogenizations => config_material%get('homogenization')
@@ -65,9 +65,9 @@ module subroutine thermal_init()
 #endif
       else
         prm%output = emptyStringArray
-      endif
+      end if
     end associate
-  enddo
+  end do
 
   call pass_init()
 
@@ -89,7 +89,7 @@ module subroutine thermal_partition(ce)
   dot_T = current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce))
   do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
     call phase_thermal_setField(T,dot_T,co,ce)
-  enddo
+  end do
 
 end subroutine thermal_partition
 
@@ -108,7 +108,7 @@ module function homogenization_mu_T(ce) result(mu)
   mu = phase_mu_T(1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     mu = mu + phase_mu_T(co,ce)
-  enddo
+  end do
 
   mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
@@ -129,7 +129,7 @@ module function homogenization_K_T(ce) result(K)
   K = phase_K_T(1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     K = K + phase_K_T(co,ce)
-  enddo
+  end do
 
   K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
@@ -150,7 +150,7 @@ module function homogenization_f_T(ce) result(f)
   f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
-  enddo
+  end do
 
   f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
@@ -189,7 +189,7 @@ module subroutine thermal_results(ho,group)
         case('T')
           call results_writeDataset(current(ho)%T,group,'T','temperature','K')
       end select
-    enddo outputsLoop
+    end do outputsLoop
   end associate
 
 end subroutine thermal_results

src/homogenization_thermal_pass.f90

@@ -8,7 +8,7 @@ contains
 
 module subroutine pass_init()
   
-  print'(/,a)', ' <<<+-  homogenization:thermal:pass init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:thermal:pass init  -+>>>'
 
 end subroutine pass_init
 

src/lattice.f90

@@ -409,7 +409,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine lattice_init
 
-  print'(/,a)', ' <<<+-  lattice init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  lattice init  -+>>>'; flush(IO_STDOUT)
 
   call selfTest
 

src/material.f90

@@ -63,11 +63,11 @@ subroutine material_init(restart)
   logical, intent(in) :: restart
 
 
-  print'(/,a)', ' <<<+-  material init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  material init  -+>>>'; flush(IO_STDOUT)
 
 
   call parse
-  print*, 'parsed material.yaml'
+  print'(/,1x,a)', 'parsed material.yaml'
 
 
   if (.not. restart) then
@@ -75,7 +75,7 @@ subroutine material_init(restart)
     call results_mapping_phase(material_phaseID,material_phaseEntry,material_name_phase)
     call results_mapping_homogenization(material_homogenizationID,material_homogenizationEntry,material_name_homogenization)
     call results_closeJobFile
-  endif
+  end if
 
 end subroutine material_init
 
@@ -115,7 +115,7 @@ subroutine parse()
   do h=1, homogenizations%length
     homogenization => homogenizations%get(h)
     homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents')
-  enddo
+  end do
   homogenization_maxNconstituents = maxval(homogenization_Nconstituents)
 
   allocate(counterPhase(phases%length),source=0)
@@ -139,7 +139,7 @@ subroutine parse()
       material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization'))
       counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1
       material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce))
-    enddo
+    end do
 
     frac = 0.0_pReal
     do co = 1, constituents%length
@@ -153,12 +153,12 @@ subroutine parse()
         material_phaseMemberAt(co,ip,el)      = counterPhase(material_phaseAt(co,el))
         material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el))
         material_phaseID(co,ce) = material_phaseAt(co,el)
-      enddo
+      end do
 
-    enddo
+    end do
     if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent')
 
-  enddo
+  end do
 
   allocate(material_O_0(materials%length))
   allocate(material_F_i_0(materials%length))
@@ -191,15 +191,15 @@ subroutine sanityCheck(materials,homogenizations)
                            constituents
   integer :: m
 
-  if(maxval(discretization_materialAt) > materials%length) &
+  if (maxval(discretization_materialAt) > materials%length) &
     call IO_error(155,ext_msg='More materials requested than found in material.yaml')
 
   do m = 1, materials%length
     material => materials%get(m)
     constituents   => material%get('constituents')
     homogenization => homogenizations%get(material%get_asString('homogenization'))
-    if(constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
+    if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
-  enddo
+  end do
 
 end subroutine sanityCheck
 
@@ -220,12 +220,12 @@ function getKeys(dict)
   do i=1, dict%length
     temp(i) = dict%getKey(i)
     l = max(len_trim(temp(i)),l)
-  enddo
+  end do
 
   allocate(character(l)::getKeys(dict%length))
   do i=1, dict%length
     getKeys(i) = trim(temp(i))
-  enddo
+  end do
 
 end function getKeys
 

src/math.f90

@@ -91,7 +91,7 @@ subroutine math_init
   class(tNode), pointer :: &
     num_generic
 
-  print'(/,a)', ' <<<+-  math init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  math init  -+>>>'; flush(IO_STDOUT)
 
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
   randomSeed  = num_generic%get_asInt('random_seed', defaultVal = 0)
@@ -109,9 +109,9 @@ subroutine math_init
   call random_seed(put = randInit)
   call random_number(randTest)
 
-  print'(a,i2)',                ' size  of random seed:     ', randSize
+  print'(/,a,i2)',                ' size  of random seed:     ', randSize
-  print'(a,i0)',                ' value of random seed:     ', randInit(1)
+  print'(  a,i0)',                ' value of random seed:     ', randInit(1)
-  print'(a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest
+  print'(  a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest
 
   call random_seed(put = randInit)
 

src/mesh/DAMASK_mesh.f90

@@ -86,7 +86,7 @@ program DAMASK_mesh
 !--------------------------------------------------------------------------------------------------
 ! init DAMASK (all modules)
   call CPFEM_initAll
-  print'(/,a)', ' <<<+-  DAMASK_mesh init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  DAMASK_mesh init  -+>>>'; flush(IO_STDOUT)
 
 !---------------------------------------------------------------------
 ! reading field information from numerics file and do sanity checks
@@ -115,16 +115,16 @@ program DAMASK_mesh
         case('$Loadcase')
           N_def = N_def + 1
       end select
-    enddo                                                                                           ! count all identifiers to allocate memory and do sanity check
+    end do                                                                                           ! count all identifiers to allocate memory and do sanity check
-  enddo
+  end do
 
-  if(N_def < 1) call IO_error(error_ID = 837)
+  if (N_def < 1) call IO_error(error_ID = 837)
   allocate(loadCases(N_def))
 
   do i = 1, size(loadCases)
     allocate(loadCases(i)%fieldBC(1))
     loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID
-  enddo
+  end do
 
   do i = 1, size(loadCases)
     loadCases(i)%fieldBC(1)%nComponents = dimPlex                                                   !< X, Y (, Z) displacements
@@ -138,12 +138,12 @@ program DAMASK_mesh
         case (3)
           loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
       end select
-    enddo
+    end do
     do component = 1, loadCases(i)%fieldBC(1)%nComponents
       allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
       allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
-    enddo
+    end do
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! reading the load case and assign values to the allocated data structure
@@ -163,7 +163,7 @@ program DAMASK_mesh
           currentFaceSet = -1
           do faceSet = 1, mesh_Nboundaries
             if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
-          enddo
+          end do
           if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
         case('t')
           loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1)
@@ -192,11 +192,11 @@ program DAMASK_mesh
                    .true.
                loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = &
                    IO_floatValue(line,chunkPos,i+1)
-             endif
+             end if
-           enddo
+           end do
       end select
-    enddo
+    end do
-  enddo
+  end do
 
 !--------------------------------------------------------------------------------------------------
 ! consistency checks and output of load case
@@ -204,28 +204,28 @@ program DAMASK_mesh
   errorID = 0
   checkLoadcases: do currentLoadCase = 1, size(loadCases)
     write (loadcase_string, '(i0)' ) currentLoadCase
-    print'(a,i0)', ' load case: ', currentLoadCase
+    print'(/,1x,a,i0)', 'load case: ', currentLoadCase
     if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
-      print'(a)', '  drop guessing along trajectory'
+      print'(2x,a)', 'drop guessing along trajectory'
-    print'(a)', '  Field '//trim(FIELD_MECH_label)
+    print'(2x,a)', 'Field '//trim(FIELD_MECH_label)
 
     do faceSet = 1, mesh_Nboundaries
        do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents
          if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) &
-           print'(a,i2,a,i2,a,f12.7)', '    Face  ', mesh_boundaries(faceSet), &
+           print'(a,i2,a,i2,a,f12.7)', &
+           '    Face ', mesh_boundaries(faceSet), &
            ' Component ', component, &
-                                       ' Value ', loadCases(currentLoadCase)%fieldBC(1)% &
+           ' Value ', loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(faceSet)
-                                                          componentBC(component)%Value(faceSet)
+      end do
-      enddo
+    end do
-    enddo
+    print'(2x,a,f12.6)', 'time:       ', loadCases(currentLoadCase)%time
-    print'(a,f12.6)', '  time:       ', loadCases(currentLoadCase)%time
     if (loadCases(currentLoadCase)%incs < 1)             errorID = 835                            ! non-positive incs count
-    print'(a,i5)',    '  increments: ', loadCases(currentLoadCase)%incs
+    print'(2x,a,i5)',    'increments: ', loadCases(currentLoadCase)%incs
     if (loadCases(currentLoadCase)%outputfrequency < 1)  errorID = 836                            ! non-positive result frequency
-    print'(a,i5)',    '  output  frequency:  ', &
+    print'(2x,a,i5)',    'output frequency: ', &
                loadCases(currentLoadCase)%outputfrequency
     if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string)                 ! exit with error message
-  enddo checkLoadcases
+  end do checkLoadcases
 
 !--------------------------------------------------------------------------------------------------
 ! doing initialization depending on active solvers
@@ -235,9 +235,9 @@ program DAMASK_mesh
   if (worldrank == 0) then
     open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
     write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded'                  ! statistics file
-  endif
+  end if
 
-  print'(/,a)', ' ... writing initial configuration to file ........................'
+  print'(/,1x,a)', '... writing initial configuration to file .................................'
   flush(IO_STDOUT)
   call CPFEM_results(0,0.0_pReal)
 
@@ -262,7 +262,7 @@ program DAMASK_mesh
 
 !--------------------------------------------------------------------------------------------------
 ! report begin of new step
-        print'(/,a)', ' ###########################################################################'
+        print'(/,1x,a)', '###########################################################################'
         print'(1x,a,es12.5,6(a,i0))',&
                 'Time', time, &
                 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
@@ -281,60 +281,60 @@ program DAMASK_mesh
         stagIterate = .true.
         do while (stagIterate)
           solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
-          if(.not. solres(1)%converged) exit
+          if (.not. solres(1)%converged) exit
 
           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
         cutBack = .False.
-        if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then                       ! no solution found
+        if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then                       ! no solution found
           if (cutBackLevel < maxCutBack) then                                                       ! do cut back
             cutBack = .True.
             stepFraction = (stepFraction - 1) * subStepFactor                                       ! adjust to new denominator
             cutBackLevel = cutBackLevel + 1
             time    = time - timeinc                                                                ! rewind time
             timeinc = timeinc/2.0_pReal
-            print'(/,a)', ' cutting back'
+            print'(/,1x,a)', 'cutting back'
           else                                                                                      ! default behavior, exit if spectral solver does not converge
             if (worldrank == 0) close(statUnit)
             call IO_error(950)
-          endif
+          end if
         else
           guess = .true.                                                                            ! start guessing after first converged (sub)inc
           timeIncOld = timeinc
-        endif
+        end if
         if (.not. cutBack .and. worldrank == 0) &
           write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
                             solres%converged, solres%iterationsNeeded                               ! write statistics about accepted solution
-      enddo subStepLooping
+      end do subStepLooping
 
       cutBackLevel = max(0, cutBackLevel - 1)                                                       ! try half number of subincs next inc
 
       if (all(solres(:)%converged)) then
-        print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged'
+        print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged'
       else
-        print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged'
+        print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged'
-      endif; flush(IO_STDOUT)
+      end if; flush(IO_STDOUT)
 
       if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then                            ! at output frequency
-        print'(/,a)', ' ... writing results to file ......................................'
+        print'(/,1x,a)', '... writing results to file ...............................................'
         call FEM_mechanical_updateCoords
         call CPFEM_results(totalIncsCounter,time)
-      endif
+      end if
 
 
-    enddo incLooping
+    end do incLooping
 
-  enddo loadCaseLooping
+  end do loadCaseLooping
 
 
 !--------------------------------------------------------------------------------------------------
 ! report summary of whole calculation
-  print'(/,a)', ' ###########################################################################'
+  print'(/,1x,a)', '###########################################################################'
   if (worldrank == 0) close(statUnit)
 
   call quit(0)                                                                                      ! no complains ;)

src/mesh/FEM_quadrature.f90

@@ -41,10 +41,10 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine FEM_quadrature_init()
 
-  print'(/,a)', ' <<<+-  FEM_quadrature init  -+>>>'; flush(6)
+  print'(/,1x,a)', '<<<+-  FEM_quadrature init  -+>>>'; flush(6)
 
-  print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
+  print'(/,1x,a)', 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
-  print*, 'https://www.jstor.org/stable/43693493'
+  print'(  1x,a)', 'https://www.jstor.org/stable/43693493'
 
 !--------------------------------------------------------------------------------------------------
 ! 2D linear

src/mesh/FEM_utilities.f90

@@ -96,7 +96,7 @@ subroutine FEM_utilities_init
   logical :: debugPETSc                                                                             !< use some in debug defined options for more verbose PETSc solution
 
 
-  print'(/,a)',   ' <<<+-  FEM_utilities init  -+>>>'
+  print'(/,1x,a)',   '<<<+-  FEM_utilities init  -+>>>'
 
   num_mesh    => config_numerics%get('mesh',defaultVal=emptyDict)
 
@@ -111,10 +111,10 @@ subroutine FEM_utilities_init
   debug_mesh  => config_debug%get('mesh',defaultVal=emptyList)
   debugPETSc  =  debug_mesh%contains('PETSc')
 
-  if(debugPETSc) print'(3(/,a),/)', &
+  if(debugPETSc) print'(3(/,1x,a),/)', &
-                 ' Initializing PETSc with debug options: ', &
+                 'Initializing PETSc with debug options: ', &
                  trim(PETScDebug), &
-                 ' add more using the "PETSc_options" keyword in numerics.yaml'
+                 'add more using the "PETSc_options" keyword in numerics.yaml'
   flush(IO_STDOUT)
   call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
   CHKERRQ(ierr)
@@ -149,8 +149,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
 
   PetscErrorCode :: ierr
 
-
+  print'(/,1x,a)', '... evaluating constitutive response ......................................'
-  print'(/,a)', ' ... evaluating constitutive response ......................................'
 
   call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
   if (.not. terminallyIll) &

src/mesh/discretization_mesh.f90

@@ -88,7 +88,7 @@ subroutine discretization_mesh_init(restart)
   integer :: p_i                                                                                    !< integration order (quadrature rule)
   type(tvec) :: coords_node0
 
-  print'(/,a)',   ' <<<+-  discretization_mesh init  -+>>>'
+  print'(/,1x,a)',   '<<<+-  discretization_mesh init  -+>>>'
 
 !--------------------------------------------------------------------------------
 ! read numerics parameter
@@ -106,6 +106,7 @@ subroutine discretization_mesh_init(restart)
   CHKERRQ(ierr)
   call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
   CHKERRQ(ierr)
+  print'()'
   call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,ierr)
   CHKERRQ(ierr)
 

src/mesh/mesh_mech_FEM.f90

@@ -113,7 +113,7 @@ subroutine FEM_mechanical_init(fieldBC)
   class(tNode), pointer :: &
     num_mesh
 
-  print'(/,a)', ' <<<+-  FEM_mech init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  FEM_mech init  -+>>>'; flush(IO_STDOUT)
 
 !-----------------------------------------------------------------------------
 ! read numerical parametes and do sanity checks
@@ -302,7 +302,7 @@ type(tSolutionState) function FEM_mechanical_solution( &
     CHKERRQ(ierr)
   endif
 
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end function FEM_mechanical_solution
@@ -663,8 +663,8 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso
   print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), &
                   ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
                                                   int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
-  print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
+  print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
-    ' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
+    'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
   flush(IO_STDOUT)
 
 end subroutine FEM_mechanical_converged

src/parallelization.f90

@@ -76,11 +76,11 @@ subroutine parallelization_init
   call MPI_Comm_rank(MPI_COMM_WORLD,worldrank,err)
   if (err /= 0)                              error stop 'Could not determine worldrank'
 
-  if (worldrank == 0) print'(/,a)',  ' <<<+-  parallelization init  -+>>>'
+  if (worldrank == 0) print'(/,1x,a)', '<<<+-  parallelization init  -+>>>'
 
   call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err)
   if (err /= 0)                              error stop 'Could not determine worldsize'
-  if (worldrank == 0) print'(a,i3)', ' MPI processes: ',worldsize
+  if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize
 
   call MPI_Type_size(MPI_INTEGER,typeSize,err)
   if (err /= 0)                              error stop 'Could not determine MPI integer size'
@@ -97,16 +97,16 @@ subroutine parallelization_init
 
 !$ call get_environment_variable(name='OMP_NUM_THREADS',value=NumThreadsString,STATUS=got_env)
 !$ if(got_env /= 0) then
-!$   print*, 'Could not get $OMP_NUM_THREADS, using default'
+!$   print'(1x,a)', 'Could not get $OMP_NUM_THREADS, using default'
 !$   OMP_NUM_THREADS = 4_pI32
 !$ else
 !$   read(NumThreadsString,'(i6)') OMP_NUM_THREADS
 !$   if (OMP_NUM_THREADS < 1_pI32) then
-!$     print*, 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default'
+!$     print'(1x,a)', 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default'
 !$     OMP_NUM_THREADS = 4_pI32
 !$   endif
 !$ endif
-!$ print'(a,i2)',   ' OMP_NUM_THREADS: ',OMP_NUM_THREADS
+!$ print'(1x,a,1x,i2)',   'OMP_NUM_THREADS:',OMP_NUM_THREADS
 !$ call omp_set_num_threads(OMP_NUM_THREADS)
 
 end subroutine parallelization_init

src/phase.f90

@@ -343,7 +343,7 @@ subroutine phase_init
     phase
 
 
-  print'(/,a)', ' <<<+-  phase init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  phase init  -+>>>'; flush(IO_STDOUT)
 
   debug_constitutive => config_debug%get('phase', defaultVal=emptyList)
   debugConstitutive%basic     = debug_constitutive%contains('basic')
@@ -371,20 +371,20 @@ subroutine phase_init
       phase_cOverA(ph) = phase%get_asFloat('c/a')
     phase_rho(ph) = phase%get_asFloat('rho',defaultVal=0.0_pReal)
     allocate(phase_O_0(ph)%data(count(material_phaseID==ph)))
-  enddo
+  end do
 
   do ce = 1, size(material_phaseID,2)
     ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
     do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
       ph = material_phaseID(co,ce)
       phase_O_0(ph)%data(material_phaseEntry(co,ce)) = material_O_0(ma)%data(co)
-    enddo
+    end do
-  enddo
+  end do
 
   allocate(phase_O(phases%length))
   do ph = 1,phases%length
     phase_O(ph)%data = phase_O_0(ph)%data
-  enddo
+  end do
 
   call mechanical_init(phases)
   call damage_init
@@ -471,7 +471,7 @@ subroutine phase_results()
     call mechanical_results(group,ph)
     call damage_results(group,ph)
 
-  enddo
+  end do
 
 end subroutine phase_results
 
@@ -495,7 +495,7 @@ subroutine crystallite_init()
     phases
 
 
-  print'(/,a)', ' <<<+-  crystallite init  -+>>>'
+  print'(/,1x,a)', '<<<+-  crystallite init  -+>>>'
 
   cMax = homogenization_maxNconstituents
   iMax = discretization_nIPs
@@ -515,28 +515,28 @@ subroutine crystallite_init()
   num%nState                 = num_crystallite%get_asInt   ('nState',           defaultVal=20)
   num%nStress                = num_crystallite%get_asInt   ('nStress',          defaultVal=40)
 
-  if(num%subStepMinCryst   <= 0.0_pReal)      call IO_error(301,ext_msg='subStepMinCryst')
+  if (num%subStepMinCryst   <= 0.0_pReal)      call IO_error(301,ext_msg='subStepMinCryst')
-  if(num%subStepSizeCryst  <= 0.0_pReal)      call IO_error(301,ext_msg='subStepSizeCryst')
+  if (num%subStepSizeCryst  <= 0.0_pReal)      call IO_error(301,ext_msg='subStepSizeCryst')
-  if(num%stepIncreaseCryst <= 0.0_pReal)      call IO_error(301,ext_msg='stepIncreaseCryst')
+  if (num%stepIncreaseCryst <= 0.0_pReal)      call IO_error(301,ext_msg='stepIncreaseCryst')
 
-  if(num%subStepSizeLp <= 0.0_pReal)          call IO_error(301,ext_msg='subStepSizeLp')
+  if (num%subStepSizeLp <= 0.0_pReal)          call IO_error(301,ext_msg='subStepSizeLp')
-  if(num%subStepSizeLi <= 0.0_pReal)          call IO_error(301,ext_msg='subStepSizeLi')
+  if (num%subStepSizeLi <= 0.0_pReal)          call IO_error(301,ext_msg='subStepSizeLi')
 
-  if(num%rtol_crystalliteState  <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState')
+  if (num%rtol_crystalliteState  <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState')
-  if(num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress')
+  if (num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress')
-  if(num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress')
+  if (num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress')
 
-  if(num%iJacoLpresiduum < 1)                 call IO_error(301,ext_msg='iJacoLpresiduum')
+  if (num%iJacoLpresiduum < 1)                 call IO_error(301,ext_msg='iJacoLpresiduum')
 
-  if(num%nState < 1)                          call IO_error(301,ext_msg='nState')
+  if (num%nState < 1)                          call IO_error(301,ext_msg='nState')
-  if(num%nStress< 1)                          call IO_error(301,ext_msg='nStress')
+  if (num%nStress< 1)                          call IO_error(301,ext_msg='nStress')
 
 
   phases => config_material%get('phase')
 
-  print'(a42,1x,i10)', '    # of elements:                       ', eMax
+  print'(/,a42,1x,i10)', '    # of elements:                       ', eMax
-  print'(a42,1x,i10)', '    # of integration points/element:     ', iMax
+  print'(  a42,1x,i10)', '    # of integration points/element:     ', iMax
-  print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax
+  print'(  a42,1x,i10)', 'max # of constituents/integration point: ', cMax
   flush(IO_STDOUT)
 
 
@@ -547,9 +547,9 @@ subroutine crystallite_init()
       do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
         call crystallite_orientations(co,ip,el)
         call plastic_dependentState(co,ip,el)                                                       ! update dependent state variables to be consistent with basic states
-     enddo
+     end do
-    enddo
+    end do
-  enddo
+  end do
   !$OMP END PARALLEL DO
 
 
@@ -642,7 +642,7 @@ subroutine phase_restartWrite(fileHandle)
 
     call HDF5_closeGroup(groupHandle(2))
 
-  enddo
+  end do
 
   call HDF5_closeGroup(groupHandle(1))
 
@@ -670,7 +670,7 @@ subroutine phase_restartRead(fileHandle)
 
     call HDF5_closeGroup(groupHandle(2))
 
-  enddo
+  end do
 
   call HDF5_closeGroup(groupHandle(1))
 

src/phase_damage.f90

@@ -83,7 +83,7 @@ module subroutine damage_init
    source
   logical:: damage_active
 
-  print'(/,a)', ' <<<+-  phase:damage init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:damage init  -+>>>'
 
   phases => config_material%get('phase')
 
@@ -108,16 +108,16 @@ module subroutine damage_init
       param(ph)%D(1,1) = source%get_asFloat('D_11')
       if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33')
       param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph))
-    endif
+    end if
 
-  enddo
+  end do
 
   allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID)
 
   if (damage_active) then
     where(isobrittle_init()  ) phase_damage = DAMAGE_ISOBRITTLE_ID
     where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID
-  endif
+  end if
 
   phase_damage_maxSizeDotState     = maxval(damageState%sizeDotState)
 
@@ -181,7 +181,7 @@ module subroutine damage_restore(ce)
     if (damageState(material_phaseID(co,ce))%sizeState > 0) &
     damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = &
       damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce))
-  enddo
+  end do
 
 end subroutine damage_restore
 
@@ -242,11 +242,11 @@ function integrateDamageState(Delta_t,ph,en) result(broken)
   if (damageState(ph)%sizeState == 0) then
     broken = .false.
     return
-  endif
+  end if
 
   converged_ = .true.
   broken = phase_damage_collectDotState(ph,en)
-  if(broken) return
+  if (broken) return
 
   size_so = damageState(ph)%sizeDotState
   damageState(ph)%state(1:size_so,en) = damageState(ph)%state0  (1:size_so,en) &
@@ -255,11 +255,11 @@ function integrateDamageState(Delta_t,ph,en) result(broken)
 
   iteration: do NiterationState = 1, num%nState
 
-    if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1)
+    if (nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1)
     source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en)
 
     broken = phase_damage_collectDotState(ph,en)
-    if(broken) exit iteration
+    if (broken) exit iteration
 
 
       zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
@@ -274,12 +274,12 @@ function integrateDamageState(Delta_t,ph,en) result(broken)
                                                damageState(ph)%atol(1:size_so))
 
 
-    if(converged_) then
+    if (converged_) then
       broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en)
       exit iteration
-    endif
+    end if
 
-  enddo iteration
+  end do iteration
 
   broken = broken .or. .not. converged_
 
@@ -302,7 +302,7 @@ function integrateDamageState(Delta_t,ph,en) result(broken)
     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
 
@@ -358,7 +358,7 @@ function phase_damage_collectDotState(ph,en) result(broken)
 
     broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en)))
 
-  endif
+  end if
 
 end function phase_damage_collectDotState
 
@@ -419,12 +419,12 @@ function phase_damage_deltaState(Fe, ph, en) result(broken)
     case (DAMAGE_ISOBRITTLE_ID) sourceType
       call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en)
       broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))
-      if(.not. broken) then
+      if (.not. broken) then
         myOffset = damageState(ph)%offsetDeltaState
         mySize   = damageState(ph)%sizeDeltaState
         damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = &
         damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en)
-      endif
+      end if
 
   end select sourceType
 
@@ -454,7 +454,7 @@ function source_active(source_label)  result(active_source)
     sources => phase%get('damage',defaultVal=emptyList)
     src => sources%get(1)
     active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label
-  enddo
+  end do
 
 
 end function source_active
@@ -497,7 +497,7 @@ module subroutine damage_forward()
   do ph = 1, size(damageState)
     if (damageState(ph)%sizeState > 0) &
       damageState(ph)%state0 = damageState(ph)%state
-  enddo
+  end do
 
 end subroutine damage_forward
 

src/phase_damage_anisobrittle.f90

@@ -46,10 +46,10 @@ module function anisobrittle_init() result(mySources)
 
 
   mySources = source_active('anisobrittle')
-  if(count(mySources) == 0) return
+  if (count(mySources) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:damage:anisobrittle init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:damage:anisobrittle init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')
@@ -57,7 +57,7 @@ module function anisobrittle_init() result(mySources)
 
 
   do ph = 1, phases%length
-    if(mySources(ph)) then
+    if (mySources(ph)) then
       phase => phases%get(ph)
       sources => phase%get('damage')
 
@@ -94,14 +94,14 @@ module function anisobrittle_init() result(mySources)
         Nmembers = count(material_phaseID==ph)
         call phase_allocateState(damageState(ph),Nmembers,1,1,0)
         damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
-        if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
+        if (any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
 
       end associate
 
       if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)')
-    endif
+    end if
 
-  enddo
+  end do
 
 end function anisobrittle_init
 
@@ -136,7 +136,7 @@ module subroutine anisobrittle_dotState(S, ph,en)
             * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
                (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
                (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**prm%q)
-    enddo
+    end do
   end associate
 
 end subroutine anisobrittle_dotState
@@ -159,7 +159,7 @@ module subroutine anisobrittle_results(phase,group)
         case ('f_phi')
           call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³')
       end select
-    enddo outputsLoop
+    end do outputsLoop
   end associate
 
 end subroutine anisobrittle_results
@@ -200,7 +200,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
       forall (k=1:3,l=1:3,m=1:3,n=1:3) &
         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
                             + dudotd_dt*prm%cleavage_systems(k,l,1,i) * prm%cleavage_systems(m,n,1,i)
-    endif
+    end if
 
     traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
     if (abs(traction_t) > traction_crit + tol_math_check) then
@@ -210,7 +210,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
       forall (k=1:3,l=1:3,m=1:3,n=1:3) &
         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
                             + dudott_dt*prm%cleavage_systems(k,l,2,i) * prm%cleavage_systems(m,n,2,i)
-    endif
+    end if
 
     traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
     if (abs(traction_n) > traction_crit + tol_math_check) then
@@ -220,8 +220,8 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
       forall (k=1:3,l=1:3,m=1:3,n=1:3) &
         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
                             + dudotn_dt*prm%cleavage_systems(k,l,3,i) * prm%cleavage_systems(m,n,3,i)
-    endif
+    end if
-  enddo
+  end do
   end associate
 
 end subroutine damage_anisobrittle_LiAndItsTangent

src/phase_damage_isobrittle.f90

@@ -44,10 +44,10 @@ module function isobrittle_init() result(mySources)
 
 
   mySources = source_active('isobrittle')
-  if(count(mySources) == 0) return
+  if (count(mySources) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:damage:isobrittle init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:damage:isobrittle init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')
@@ -56,7 +56,7 @@ module function isobrittle_init() result(mySources)
   allocate(deltaState(phases%length))
 
   do ph = 1, phases%length
-    if(mySources(ph)) then
+    if (mySources(ph)) then
       phase => phases%get(ph)
       sources => phase%get('damage')
 
@@ -77,7 +77,7 @@ module function isobrittle_init() result(mySources)
         Nmembers = count(material_phaseID==ph)
         call phase_allocateState(damageState(ph),Nmembers,1,1,1)
         damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
-        if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
+        if (any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
 
         stt%r_W => damageState(ph)%state(1,:)
         dlt%r_W => damageState(ph)%deltaState(1,:)
@@ -86,9 +86,9 @@ module function isobrittle_init() result(mySources)
 
 
       if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)')
-    endif
+    end if
 
-  enddo
+  end do
 
 
 end function isobrittle_init
@@ -141,7 +141,7 @@ module subroutine isobrittle_results(phase,group)
         case ('f_phi')
           call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless
       end select
-    enddo outputsLoop
+    end do outputsLoop
 
   end associate
 

src/phase_mechanical.f90

@@ -217,7 +217,7 @@ module subroutine mechanical_init(phases)
     phase, &
     mech
 
-  print'(/,a)', ' <<<+-  phase:mechanical init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical init  -+>>>'
 
 !-------------------------------------------------------------------------------------------------
   allocate(output_constituent(phases%length))

src/phase_mechanical_eigen.f90

@@ -44,7 +44,7 @@ module subroutine eigen_init(phases)
     kinematics, &
     mechanics
 
-  print'(/,a)', ' <<<+-  phase:mechanical:eigen init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:eigen init  -+>>>'
 
 !--------------------------------------------------------------------------------------------------
 ! explicit eigen mechanisms
@@ -55,11 +55,11 @@ module subroutine eigen_init(phases)
     mechanics => phase%get('mechanical')
     kinematics => mechanics%get('eigen',defaultVal=emptyList)
     Nmodels(ph) = kinematics%length
-  enddo
+  end do
 
   allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID)
 
-  if(maxval(Nmodels) /= 0) then
+  if (maxval(Nmodels) /= 0) then
     where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID
   endif
 
@@ -97,8 +97,8 @@ function kinematics_active(kinematics_label,kinematics_length)  result(active_ki
     do k = 1, kinematics%length
       kinematics_type => kinematics%get(k)
       active_kinematics(k,ph) = kinematics_type%get_asString('type') == kinematics_label
-    enddo
+    end do
-  enddo
+  end do
 
 
 end function kinematics_active
@@ -125,11 +125,11 @@ function kinematics_active2(kinematics_label)  result(active_kinematics)
   do ph = 1, phases%length
     phase => phases%get(ph)
     kinematics => phase%get('damage',defaultVal=emptyList)
-    if(kinematics%length < 1) return
+    if (kinematics%length < 1) return
     kinematics_type => kinematics%get(1)
     if (.not. kinematics_type%contains('type')) continue
     active_kinematics(ph) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label
-  enddo
+  end do
 
 
 end function kinematics_active2
@@ -188,7 +188,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
         dLi_dS = dLi_dS + my_dLi_dS
         active = .true.
     end select kinematicsType
-  enddo KinematicsLoop
+  end do KinematicsLoop
 
   select case (model_damage(ph))
     case (KINEMATICS_cleavage_opening_ID)
@@ -198,7 +198,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
       active = .true.
   end select
 
-  if(.not. active) return
+  if (.not. active) return
 
   FiInv = math_inv33(Fi)
   detFi = math_det33(Fi)
@@ -209,7 +209,7 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
     dLi_dS(1:3,1:3,i,j)  = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
     dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
     dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
-  enddo; enddo
+  end do; end do
 
 end subroutine phase_LiAndItsTangents
 

src/phase_mechanical_eigen_cleavageopening.f90

@@ -21,8 +21,8 @@ module function damage_anisobrittle_init() result(myKinematics)
   myKinematics = kinematics_active2('anisobrittle')
   if(count(myKinematics) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:eigen:cleavageopening init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:eigen:cleavageopening init  -+>>>'
-  print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
 
 end function damage_anisobrittle_init
 

src/phase_mechanical_eigen_thermalexpansion.f90

@@ -36,25 +36,25 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
     kinematics, &
     kinematic_type
 
-  print'(/,a)', ' <<<+-  phase:mechanical:eigen:thermalexpansion init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:eigen:thermalexpansion init  -+>>>'
 
   myKinematics = kinematics_active('thermalexpansion',kinematics_length)
   Ninstances = count(myKinematics)
-  print'(a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
+  if (Ninstances == 0) return
 
   phases => config_material%get('phase')
   allocate(param(Ninstances))
   allocate(kinematics_thermal_expansion_instance(phases%length), source=0)
 
   do p = 1, phases%length
-    if(any(myKinematics(:,p))) kinematics_thermal_expansion_instance(p) = count(myKinematics(:,1:p))
+    if (any(myKinematics(:,p))) kinematics_thermal_expansion_instance(p) = count(myKinematics(:,1:p))
     phase => phases%get(p)
-    if(count(myKinematics(:,p)) == 0) cycle
+    if (count(myKinematics(:,p)) == 0) cycle
     mech => phase%get('mechanical')
     kinematics => mech%get('eigen')
     do k = 1, kinematics%length
-      if(myKinematics(k,p)) then
+      if (myKinematics(k,p)) then
         associate(prm  => param(kinematics_thermal_expansion_instance(p)))
           kinematic_type => kinematics%get(k)
 
@@ -67,14 +67,14 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
             prm%A(3,3,1) = kinematic_type%get_asFloat('A_33')
             prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T',defaultVal=0.0_pReal)
             prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal)
-          endif
+          end if
           do i=1, size(prm%A,3)
             prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p))
-          enddo
+          end do
         end associate
-      endif
+      end if
-    enddo
+    end do
-  enddo
+  end do
 
 end function thermalexpansion_init
 

src/phase_mechanical_elastic.f90

@@ -26,10 +26,10 @@ module subroutine elastic_init(phases)
     elastic
 
 
-  print'(/,a)', ' <<<+-  phase:mechanical:elastic init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:elastic init  -+>>>'
-  print'(/,a)', ' <<<+-  phase:mechanical:elastic:Hooke init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:elastic:Hooke init  -+>>>'
 
-  print'(a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT)
 
   allocate(param(phases%length))
 
@@ -48,7 +48,7 @@ module subroutine elastic_init(phases)
       if (any(phase_lattice(ph) == ['hP','tI'])) then
         prm%C66(1,3) = elastic%get_asFloat('C_13')
         prm%C66(3,3) = elastic%get_asFloat('C_33')
-      endif
+      end if
       if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66')
 
       prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph))
@@ -59,7 +59,7 @@ module subroutine elastic_init(phases)
       prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66))                                       ! Literature data is in Voigt notation
 
     end associate
-  enddo
+  end do
 
 end subroutine elastic_init
 
@@ -98,7 +98,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
   do i =1, 3;do j=1,3
     dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi)))                !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
     dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3))                             !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
-  enddo; enddo
+  end do; end do
 
 end subroutine phase_hooke_SandItsTangents
 

src/phase_mechanical_plastic.f90

@@ -222,7 +222,7 @@ contains
 module subroutine plastic_init
 
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic init  -+>>>'
 
   where(plastic_none_init())              phase_plasticity = PLASTICITY_NONE_ID
   where(plastic_isotropic_init())         phase_plasticity = PLASTICITY_ISOTROPIC_ID

src/phase_mechanical_plastic_dislotungsten.f90

@@ -99,13 +99,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
 
 
   myPlasticity = plastic_active('dislotungsten')
-  if(count(myPlasticity) == 0) return
+  if (count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:dislotungsten init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:dislotungsten init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
-  print*, 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
+  print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
-  print*, 'https://doi.org/10.1016/j.ijplas.2015.09.002'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
 
 
   phases => config_material%get('phase')
@@ -116,7 +116,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
 
 
   do ph = 1, phases%length
-    if(.not. myPlasticity(ph)) cycle
+    if (.not. myPlasticity(ph)) cycle
 
     associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
 
@@ -243,7 +243,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
     dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
-    if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
+    if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
 
     allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal)
     allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),  source=0.0_pReal)

src/phase_mechanical_plastic_dislotwin.f90

@@ -150,17 +150,17 @@ module function plastic_dislotwin_init() result(myPlasticity)
   myPlasticity = plastic_active('dislotwin')
   if(count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:dislotwin init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:dislotwin init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
-  print*, 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004'
+  print'(/,1x,a)', 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004'
-  print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
 
-  print*, 'F. Roters et al., Computational Materials Science 39:91–95, 2007'
+  print'(/,1x,a)', 'F. Roters et al., Computational Materials Science 39:91–95, 2007'
-  print*, 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL
 
-  print*, 'S.L. Wong et al., Acta Materialia 118:140–151, 2016'
+  print'(/,1x,a)', 'S.L. Wong et al., Acta Materialia 118:140–151, 2016'
-  print*, 'https://doi.org/10.1016/j.actamat.2016.07.032'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2016.07.032'
 
 
   phases => config_material%get('phase')
@@ -256,7 +256,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
       rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
       allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
       allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0))
-    endif slipActive
+    end if slipActive
 
 !--------------------------------------------------------------------------------------------------
 ! twin related parameters
@@ -283,7 +283,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
       if (.not. prm%fccTwinTransNucleation) then
         prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw')
         prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,N_tw)
-      endif
+      end if
 
       ! expand: family => system
       prm%b_tw = math_expand(prm%b_tw,N_tw)
@@ -299,11 +299,11 @@ module function plastic_dislotwin_init() result(myPlasticity)
       if (any(prm%r             < 0.0_pReal)) extmsg = trim(extmsg)//' p_tw'
       if (.not. prm%fccTwinTransNucleation) then
         if (any(prm%dot_N_0_tw  < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tw'
-      endif
+      end if
     else twinActive
       allocate(prm%gamma_char,prm%b_tw,prm%dot_N_0_tw,prm%t_tw,prm%r,source=emptyRealArray)
       allocate(prm%h_tw_tw(0,0))
-    endif twinActive
+    end if twinActive
 
 !--------------------------------------------------------------------------------------------------
 ! transformation related parameters
@@ -335,7 +335,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
       if (phase_lattice(ph) /= 'cF') then
         prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr')
         prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr)
-      endif
+      end if
       prm%t_tr = pl%get_as1dFloat('t_tr')
       prm%t_tr = math_expand(prm%t_tr,N_tr)
       prm%s    = pl%get_as1dFloat('p_tr',defaultVal=[0.0_pReal])
@@ -349,11 +349,11 @@ module function plastic_dislotwin_init() result(myPlasticity)
       if (any(prm%s             < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr'
       if (phase_lattice(ph) /= 'cF') then
         if (any(prm%dot_N_0_tr  < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr'
-      endif
+      end if
     else transActive
       allocate(prm%s,prm%b_tr,prm%t_tr,prm%dot_N_0_tr,source=emptyRealArray)
       allocate(prm%h_tr_tr(0,0))
-    endif transActive
+    end if transActive
 
 !--------------------------------------------------------------------------------------------------
 ! shearband related parameters
@@ -369,7 +369,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
       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'
-    endif
+    end if
 
 !--------------------------------------------------------------------------------------------------
 ! parameters required for several mechanisms and their interactions
@@ -383,19 +383,19 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%T_ref       = pl%get_asFloat('T_ref')
       prm%Gamma_sf(1) = pl%get_asFloat('Gamma_sf')
       prm%Gamma_sf(2) = pl%get_asFloat('Gamma_sf,T',defaultVal=0.0_pReal)
-    endif
+    end if
 
     slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
       prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
                                                    phase_lattice(ph))
       if (prm%fccTwinTransNucleation .and. size(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,N_tr,pl%get_as1dFloat('h_sl-tr'), &
                                                     phase_lattice(ph))
       if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
-    endif slipAndTransActive
+    end if slipAndTransActive
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -465,7 +465,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
 !  exit if any parameter is out of range
     if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(dislotwin)')
 
-  enddo
+  end do
 
 end function plastic_dislotwin_init
 
@@ -494,11 +494,11 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
     do i=1,prm%sum_N_tw
       homogenizedC = homogenizedC &
                    + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
-    enddo
+    end do
     do i=1,prm%sum_N_tr
       homogenizedC = homogenizedC &
                    + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
-    enddo
+    end do
 
   end associate
 
@@ -566,7 +566,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
                        + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
-  enddo slipContribution
+  end do slipContribution
 
   call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
   twinContibution: do i = 1, prm%sum_N_tw
@@ -574,7 +574,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
                        + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
-  enddo twinContibution
+  end do twinContibution
 
   call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
   transContibution: do i = 1, prm%sum_N_tr
@@ -582,7 +582,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
                        + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
-  enddo transContibution
+  end do transContibution
 
   Lp      = Lp      * f_unrotated
   dLp_dMp = dLp_dMp * f_unrotated
@@ -608,10 +608,10 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
         forall (k=1:3,l=1:3,m=1:3,n=1:3) &
           dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
                            + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n)
-      endif significantShearBandStress
+      end if significantShearBandStress
-    enddo
+    end do
 
-  endif shearBandingContribution
+  end if shearBandingContribution
 
   end associate
 
@@ -686,9 +686,9 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
 
           dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
                                / (d_hat-prm%d_caron(i))
-        endif
+        end if
-      endif significantSlipStress
+      end if significantSlipStress
-    enddo slipState
+    end do slipState
 
     dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) &
                       - dot_rho_dip_formation &
@@ -833,7 +833,7 @@ module subroutine plastic_dislotwin_results(ph,group)
 
       end select
 
-    enddo
+    end do
 
   end associate
 
@@ -968,8 +968,8 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
         end if
       else isFCC
         Ndot0=prm%dot_N_0_tw(i)
-      endif isFCC
+      end if isFCC
-    enddo
+    end do
 
     significantStress: where(tau > tol_math_check)
       StressRatio_r   = (dst%tau_hat_tw(:,en)/tau)**prm%r
@@ -1037,8 +1037,8 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
         end if
       else isFCC
         Ndot0=prm%dot_N_0_tr(i)
-      endif isFCC
+      end if isFCC
-    enddo
+    end do
 
     significantStress: where(tau > tol_math_check)
       StressRatio_s   = (dst%tau_hat_tr(:,en)/tau)**prm%s

src/phase_mechanical_plastic_isotropic.f90

@@ -68,11 +68,11 @@ module function plastic_isotropic_init() result(myPlasticity)
   myPlasticity = plastic_active('isotropic')
   if(count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:isotropic init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:isotropic init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
-  print*, 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
+  print'(/,a)', 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
-  print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
+  print'(/,a)', 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
@@ -140,7 +140,7 @@ module function plastic_isotropic_init() result(myPlasticity)
 !  exit if any parameter is out of range
     if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(isotropic)')
 
-  enddo
+  end do
 
 end function plastic_isotropic_init
 
@@ -232,7 +232,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
     else
       Li      = 0.0_pReal
       dLi_dMi = 0.0_pReal
-    endif
+    end if
 
   end associate
 
@@ -262,7 +262,7 @@ module subroutine isotropic_dotState(Mp,ph,en)
     norm_Mp = sqrt(math_tensordot(Mp,Mp))
   else
     norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp)))
-  endif
+  end if
 
   dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n
 
@@ -273,13 +273,13 @@ module subroutine isotropic_dotState(Mp,ph,en)
       xi_inf_star = prm%xi_inf &
                   + asinh( (dot_gamma / prm%c_1)**(1.0_pReal / prm%c_2))**(1.0_pReal / prm%c_3) &
                   / prm%c_4 * (dot_gamma / prm%dot_gamma_0)**(1.0_pReal / prm%n)
-    endif
+    end if
     dot%xi(en) = dot_gamma &
                * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
                * sign(abs(1.0_pReal - stt%xi(en)/xi_inf_star)**prm%a *prm%h, 1.0_pReal-stt%xi(en)/xi_inf_star)
   else
     dot%xi(en) = 0.0_pReal
-  endif
+  end if
 
   end associate
 
@@ -303,7 +303,7 @@ module subroutine plastic_isotropic_results(ph,group)
         call results_writeDataset(stt%xi,group,trim(prm%output(o)), &
                                     'resistance against plastic flow','Pa')
     end select
-  enddo outputsLoop
+  end do outputsLoop
   end associate
 
 end subroutine plastic_isotropic_results

src/phase_mechanical_plastic_kinehardening.f90

@@ -83,8 +83,8 @@ module function plastic_kinehardening_init() result(myPlasticity)
   myPlasticity = plastic_active('kinehardening')
   if(count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:kinehardening init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:kinehardening init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')
@@ -95,7 +95,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
 
 
   do ph = 1, phases%length
-    if(.not. myPlasticity(ph)) cycle
+    if (.not. myPlasticity(ph)) cycle
 
     associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph))
 
@@ -125,7 +125,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
       else
         prm%P_nS_pos = prm%P
         prm%P_nS_neg = prm%P
-      endif
+      end if
       prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
                                                    phase_lattice(ph))
 
@@ -161,7 +161,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
       xi_0 = emptyRealArray
       allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray)
       allocate(prm%h_sl_sl(0,0))
-    endif slipActive
+    end if slipActive
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -217,7 +217,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
 !  exit if any parameter is out of range
     if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(kinehardening)')
 
-  enddo
+  end do
 
 
 end function plastic_kinehardening_init
@@ -258,7 +258,7 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
                        + ddot_gamma_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) &
                        + ddot_gamma_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i)
-  enddo
+  end do
 
   end associate
 
@@ -382,7 +382,7 @@ module subroutine plastic_kinehardening_results(ph,group)
                                     'plastic shear','1',prm%systems_sl)
       end select
 
-    enddo
+    end do
 
   end associate
 
@@ -424,7 +424,7 @@ pure subroutine kinetics(Mp,ph,en, &
       tau_pos(i) =       math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%chi(i,en)
       tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%chi(i,en), &
                          0.0_pReal, prm%nonSchmidActive)
-    enddo
+    end do
 
     where(dNeq0(tau_pos))
       dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &           ! 1/2 if non-Schmid active
@@ -446,14 +446,14 @@ pure subroutine kinetics(Mp,ph,en, &
       else where
         ddot_gamma_dtau_pos = 0.0_pReal
       end where
-    endif
+    end if
     if (present(ddot_gamma_dtau_neg)) then
       where(dNeq0(dot_gamma_neg))
         ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg
       else where
         ddot_gamma_dtau_neg = 0.0_pReal
       end where
-    endif
+    end if
 
   end associate
 

src/phase_mechanical_plastic_none.f90

@@ -22,16 +22,16 @@ module function plastic_none_init() result(myPlasticity)
 
 
   myPlasticity = plastic_active('none')
-  if(count(myPlasticity) == 0) return
+  if (count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:none init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:none init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
   phases => config_material%get('phase')
   do ph = 1, phases%length
-    if(.not. myPlasticity(ph)) cycle
+    if (.not. myPlasticity(ph)) cycle
     call phase_allocateState(plasticState(ph),count(material_phaseID == ph),0,0,0)
-  enddo
+  end do
 
 end function plastic_none_init
 

src/phase_mechanical_plastic_nonlocal.f90

@@ -197,19 +197,19 @@ module function plastic_nonlocal_init() result(myPlasticity)
 
   myPlasticity = plastic_active('nonlocal')
   Ninstances = count(myPlasticity)
-  if(Ninstances == 0) then
+  if (Ninstances == 0) then
     call geometry_plastic_nonlocal_disable
     return
-  endif
+  end if
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:nonlocal init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:nonlocal init  -+>>>'
-  print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
 
-  print*, 'C. Reuber et al., Acta Materialia 71:333–348, 2014'
+  print'(/,1x,a)', 'C. Reuber et al., Acta Materialia 71:333–348, 2014'
-  print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
 
-  print*, 'C. Kords, Dissertation RWTH Aachen, 2014'
+  print'(/,1x,a)', 'C. Kords, Dissertation RWTH Aachen, 2014'
-  print*, 'http://publications.rwth-aachen.de/record/229993'
+  print'(  1x,a)', 'http://publications.rwth-aachen.de/record/229993'
 
 
   phases => config_material%get('phase')
@@ -224,7 +224,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
   allocate(dependentState(phases%length))
 
   do ph = 1, phases%length
-    if(.not. myPlasticity(ph)) cycle
+    if (.not. myPlasticity(ph)) cycle
 
     associate(prm => param(ph),  dot => dotState(ph),   stt => state(ph), &
               st0 => state0(ph), del => deltaState(ph), dst => dependentState(ph))
@@ -259,7 +259,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
       else
         prm%P_nS_pos = prm%P_sl
         prm%P_nS_neg = prm%P_sl
-      endif
+      end if
 
       prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), &
                                                    phase_lattice(ph))
@@ -280,8 +280,8 @@ module function plastic_nonlocal_init() result(myPlasticity)
            if (all(dEq0 (math_cross(prm%slip_direction(1:3,s1),prm%slip_direction(1:3,s2)))) .and. &
                any(dNeq0(math_cross(prm%slip_normal   (1:3,s1),prm%slip_normal   (1:3,s2))))) &
              prm%colinearSystem(s1) = s2
-        enddo
+        end do
-      enddo
+      end do
 
       ini%rho_u_ed_pos_0  = pl%get_as1dFloat('rho_u_ed_pos_0',   requiredSize=size(ini%N_sl))
       ini%rho_u_ed_neg_0  = pl%get_as1dFloat('rho_u_ed_neg_0',   requiredSize=size(ini%N_sl))
@@ -391,7 +391,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
       if (prm%f_ed_mult    < 0.0_pReal .or. prm%f_ed_mult    > 1.0_pReal) &
                                                   extmsg = trim(extmsg)//' f_ed_mult'
 
-    endif slipActive
+    end if slipActive
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -506,7 +506,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
 !  exit if any parameter is out of range
     if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(nonlocal)')
 
-  enddo
+  end do
 
   allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,&
                          discretization_nIPs,discretization_Nelems), source=0.0_pReal)
@@ -527,24 +527,24 @@ module function plastic_nonlocal_init() result(myPlasticity)
       do s = 1,param(ph)%sum_N_sl
         l = l + 1
         iRhoU(s,t,ph) = l
-      enddo
+      end do
-    enddo
+    end do
     l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest
     do t = 1,4
       do s = 1,param(ph)%sum_N_sl
         l = l + 1
         iV(s,t,ph) = l
-      enddo
+      end do
-    enddo
+    end do
     do t = 1,2
       do s = 1,param(ph)%sum_N_sl
         l = l + 1
         iD(s,t,ph) = l
-      enddo
+      end do
-    enddo
+    end do
     if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) &
       error stop 'state indices not properly set (nonlocal)'
-  enddo
+  end do
 
 end function plastic_nonlocal_init
 
@@ -625,7 +625,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
                                    / log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
   else
     myInteractionMatrix = prm%h_sl_sl
-  endif
+  end if
 
   dst%tau_pass(:,en) = prm%mu * prm%b_sl &
                      * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
@@ -680,14 +680,14 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
           connection_latticeConf(1:3,n) = 0.0_pReal
           rho_edg_delta_neighbor(:,n) = rho_edg_delta
           rho_scr_delta_neighbor(:,n) = rho_scr_delta
-        endif
+        end if
       else
         ! free surface -> use central values instead
         connection_latticeConf(1:3,n) = 0.0_pReal
         rho_edg_delta_neighbor(:,n) = rho_edg_delta
         rho_scr_delta_neighbor(:,n) = rho_scr_delta
-      endif
+      end if
-    enddo
+    end do
 
     neighbor_rhoExcess(1,:,:) = rho_edg_delta_neighbor
     neighbor_rhoExcess(2,:,:) = rho_scr_delta_neighbor
@@ -709,12 +709,12 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
                                - connection_latticeConf(1:3,neighbors(2))
           rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) &
                                     - neighbor_rhoExcess(c,s,neighbors(2))
-        enddo
+        end do
         invConnections = math_inv33(connections)
         if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error')
 
         rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences))
-      enddo
+      end do
 
         ! ... plus gradient from deads ...
       rhoExcessGradient(1) = rhoExcessGradient(1) + sum(rho(s,imm_edg)) / FVsize
@@ -731,8 +731,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
       dst%tau_back(s,en) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
                          * (  rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) &
                             + rhoExcessGradient_over_rho(2))
-    enddo
+    end do
-  endif
+  end if
 
  end associate
 
@@ -790,8 +790,8 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
     else
       tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
       tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
-    endif
+    end if
-  enddo
+  end do
   tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
   tau   = tau   + dst%tau_back(:,en)
 
@@ -807,12 +807,12 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
     do t = 3,4
       call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
                     tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph)
-    enddo
+    end do
   else
     v(:,3:4)         = spread(v(:,1),2,2)
     dv_dtau(:,3:4)   = spread(dv_dtau(:,1),2,2)
     dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2)
-  endif
+  end if
 
   stt%v(:,en) = pack(v,.true.)
 
@@ -833,7 +833,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
           + prm%P_sl(i,j,s) &
           * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
              - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4))  * prm%b_sl(s)
-  enddo
+  end do
 
   end associate
 
@@ -899,7 +899,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
   do s = 1,prm%sum_N_sl
     tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) +dst%tau_back(s,en)
     if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
-  enddo
+  end do
 
   dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
   dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
@@ -980,7 +980,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   if (timestep <= 0.0_pReal) then
     plasticState(ph)%dotState = 0.0_pReal
     return
-  endif
+  end if
 
   associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph))
 
@@ -1002,7 +1002,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   do s = 1,prm%sum_N_sl
     tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en)
     if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
-  enddo
+  end do
 
   dLower = prm%minDipoleHeight
   dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
@@ -1032,7 +1032,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
     rhoDotMultiplication(:,1:4) = spread( &
           (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) &
         * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4)                                   ! eq. 3.26
-  endif isBCC
+  end if isBCC
 
   forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
 
@@ -1062,7 +1062,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
                                     + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) &
                                     - rhoDotSingle2DipoleGlide(:,2*c-1) &
                                     - rhoDotSingle2DipoleGlide(:,2*c)
-  enddo
+  end do
 
 
   ! athermal annihilation
@@ -1103,13 +1103,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
     if (debugConstitutive%extensive) then
       print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
       print'(a)', '<< CONST >> enforcing cutback !!!'
-    endif
+    end if
 #endif
     plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
   else
     dot%rho(:,en) = pack(rhoDot,.true.)
     dot%gamma(:,en) = sum(dot_gamma,2)
-  endif
+  end if
 
   end associate
 
@@ -1217,11 +1217,11 @@ function rhoDotFlux(timestep,ph,en,ip,el)
                              > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
         ' at a timestep of ',timestep
       print*, '<< CONST >> enforcing cutback !!!'
-    endif
+    end if
 #endif
       rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN)                                             ! enforce cutback
       return
-    endif
+    end if
 
 
     !*** be aware of the definition of slip_transverse = slip_direction x slip_normal !!!
@@ -1255,7 +1255,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
         Favg = 0.5_pReal * (my_F + neighbor_F)
       else                                                                                          ! if no neighbor, take my value as average
         Favg = my_F
-      endif
+      end if
 
       neighbor_v0 = 0.0_pReal        ! needed for check of sign change in flux density below
 
@@ -1300,10 +1300,10 @@ function rhoDotFlux(timestep,ph,en,ip,el)
                 rhoDotFlux(:,topp) = rhoDotFlux(:,topp) &
                                    + lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2.0_pReal ! transferring to opposite signed mobile dislocation type
 
-            endif
+            end if
-          enddo
+          end do
-        enddo
+        end do
-      endif; endif
+      end if; end if
 
 
       !* FLUX FROM ME TO MY NEIGHBOR
@@ -1330,20 +1330,20 @@ function rhoDotFlux(timestep,ph,en,ip,el)
                 transmissivity = sum(compatibility(c,:,s,n,ip,el)**2.0_pReal)                       ! overall transmissivity from this slip system to my neighbor
               else                                                                                  ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
                 transmissivity = 0.0_pReal
-              endif
+              end if
               lineLength = my_rhoSgl0(s,t) * v0(s,t) &
                          * math_inner(m(1:3,s,t), normal_me2neighbor) * area                        ! positive line length of mobiles that wants to leave through this interface
               rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el)                      ! subtract dislocation flux from current type
               rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
                                 + lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) &
                                 * sign(1.0_pReal, v0(s,t))                                          ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
-            endif
+            end if
-          enddo
+          end do
-        enddo
+        end do
-      endif; endif
+      end if; end if
 
-    enddo neighbors
+    end do neighbors
-  endif
+  end if
 
   end associate
 
@@ -1433,7 +1433,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
                                                        mis%rotate(prm%slip_normal(1:3,s2)))) &
                                       * abs(math_inner(prm%slip_direction(1:3,s1), &
                                                        mis%rotate(prm%slip_direction(1:3,s2))))
-        enddo neighborSlipSystems
+        end do neighborSlipSystems
 
         my_compatibilitySum = 0.0_pReal
         belowThreshold = .true.
@@ -1446,15 +1446,15 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
               my_compatibility(:,:,s1,n) = sign((1.0_pReal - my_compatibilitySum)/nThresholdValues,&
                                                  my_compatibility(:,:,s1,n))
           my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
-        enddo
+        end do
 
         where(belowThreshold) my_compatibility(1,:,s1,n) = 0.0_pReal
         where(belowThreshold) my_compatibility(2,:,s1,n) = 0.0_pReal
 
-      enddo mySlipSystems
+      end do mySlipSystems
-    endif
+    end if
 
-  enddo neighbors
+  end do neighbors
 
   compatibility(:,:,:,:,i,e) = my_compatibility
 
@@ -1532,7 +1532,7 @@ module subroutine plastic_nonlocal_results(ph,group)
                                     'passing stress for slip','Pa', prm%systems_sl)
       end select
 
-    enddo
+    end do
 
   end associate
 
@@ -1581,7 +1581,7 @@ subroutine stateInit(ini,phase,Nentries)
         s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
         meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume
         stt%rhoSglMobile(s,e) = densityBinning
-      enddo
+      end do
     else                                ! homogeneous distribution with noise
       do e = 1, Nentries
         do f = 1,size(ini%N_sl,1)
@@ -1594,12 +1594,12 @@ subroutine stateInit(ini,phase,Nentries)
             stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1)
             stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2)
             stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2)
-          enddo
+          end do
           stt%rho_dip_edg(from:upto,e)   = ini%rho_d_ed_0(f)
           stt%rho_dip_scr(from:upto,e)   = ini%rho_d_sc_0(f)
-        enddo
+        end do
-      enddo
+      end do
-    endif
+    end if
 
   end associate
 
@@ -1688,8 +1688,8 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
       dv_dtau(s)   = v(s)**2.0_pReal * (dtSolidSolution_dtau / lambda_S + prm%B / (prm%b_sl(s) * tauEff**2.0_pReal))
       dv_dtauNS(s) = v(s)**2.0_pReal * dtPeierls_dtau / lambda_P
 
-    endif
+    end if
-  enddo
+  end do
 
   end associate
 
@@ -1760,8 +1760,8 @@ subroutine storeGeometry(ph)
   do ce = 1, size(material_homogenizationEntry,1)
     do co = 1, homogenization_maxNconstituents
       if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
-    enddo
+    end do
-  enddo
+  end do
 
 end subroutine
 

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -95,8 +95,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
   myPlasticity = plastic_active('phenopowerlaw')
   if(count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:phenopowerlaw init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:phenopowerlaw init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')
@@ -105,7 +105,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
   allocate(dotState(phases%length))
 
   do ph = 1, phases%length
-    if(.not. myPlasticity(ph)) cycle
+    if (.not. myPlasticity(ph)) cycle
 
     associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
 
@@ -129,7 +129,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
       else
         prm%P_nS_pos = prm%P_sl
         prm%P_nS_neg = prm%P_sl
-      endif
+      end if
       prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'),phase_lattice(ph))
 
       xi_0_sl             = pl%get_as1dFloat('xi_0_sl',   requiredSize=size(N_sl))
@@ -158,7 +158,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
       xi_0_sl = emptyRealArray
       allocate(prm%xi_inf_sl,prm%h_int,source=emptyRealArray)
       allocate(prm%h_sl_sl(0,0))
-    endif slipActive
+    end if slipActive
 
 !--------------------------------------------------------------------------------------------------
 ! twin related parameters
@@ -192,7 +192,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
       xi_0_tw = emptyRealArray
       allocate(prm%gamma_char,source=emptyRealArray)
       allocate(prm%h_tw_tw(0,0))
-    endif twinActive
+    end if twinActive
 
 !--------------------------------------------------------------------------------------------------
 ! slip-twin related parameters
@@ -206,7 +206,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
       allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw))                                              ! at least one dimension is 0
       allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl))                                              ! at least one dimension is 0
       prm%h_0_tw_sl = 0.0_pReal
-    endif slipAndTwinActive
+    end if slipAndTwinActive
 
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
@@ -263,7 +263,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 !  exit if any parameter is out of range
     if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(phenopowerlaw)')
 
-  enddo
+  end do
 
 end function plastic_phenopowerlaw_init
 
@@ -306,7 +306,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
                        + ddot_gamma_dtau_sl_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) &
                        + ddot_gamma_dtau_sl_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i)
-  enddo slipSystems
+  end do slipSystems
 
   call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
   twinSystems: do i = 1, prm%sum_N_tw
@@ -314,7 +314,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
                        + ddot_gamma_dtau_tw(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
-  enddo twinSystems
+  end do twinSystems
 
   end associate
 
@@ -397,7 +397,7 @@ module subroutine plastic_phenopowerlaw_results(ph,group)
 
       end select
 
-    enddo
+    end do
 
   end associate
 
@@ -438,7 +438,7 @@ pure subroutine kinetics_sl(Mp,ph,en, &
       tau_sl_pos(i) =       math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i))
       tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), &
                             0.0_pReal, prm%nonSchmidActive)
-    enddo
+    end do
 
     where(dNeq0(tau_sl_pos))
       dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &     ! 1/2 if non-Schmid active
@@ -460,14 +460,14 @@ pure subroutine kinetics_sl(Mp,ph,en, &
       else where
         ddot_gamma_dtau_sl_pos = 0.0_pReal
       end where
-    endif
+    end if
     if (present(ddot_gamma_dtau_sl_neg)) then
       where(dNeq0(dot_gamma_sl_neg))
         ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg
       else where
         ddot_gamma_dtau_sl_neg = 0.0_pReal
       end where
-    endif
+    end if
 
   end associate
 
@@ -517,7 +517,7 @@ pure subroutine kinetics_tw(Mp,ph,en,&
       else where
         ddot_gamma_dtau_tw = 0.0_pReal
       end where
-    endif
+    end if
 
   end associate
 

src/phase_thermal.f90

@@ -86,7 +86,7 @@ module subroutine thermal_init(phases)
     Nmembers
 
 
-  print'(/,a)', ' <<<+-  phase:thermal init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:thermal init  -+>>>'
 
   allocate(current(phases%length))
 

src/phase_thermal_dissipation.f90

@@ -36,8 +36,8 @@ module function dissipation_init(source_length) result(mySources)
 
   mySources = thermal_active('dissipation',source_length)
   if(count(mySources) == 0) return
-  print'(/,a)', ' <<<+-  phase:thermal:dissipation init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:thermal:dissipation init  -+>>>'
-  print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')
@@ -45,11 +45,11 @@ module function dissipation_init(source_length) result(mySources)
 
   do ph = 1, phases%length
     phase => phases%get(ph)
-    if(count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1
+    if (count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1
     thermal => phase%get('thermal')
     sources => thermal%get('source')
     do so = 1, sources%length
-      if(mySources(so,ph)) then
+      if (mySources(so,ph)) then
         associate(prm  => param(ph))
           src => sources%get(so)
 
@@ -58,9 +58,9 @@ module function dissipation_init(source_length) result(mySources)
           call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
 
         end associate
-      endif
+      end if
-    enddo
+    end do
-  enddo
+  end do
 
 
 end function dissipation_init

src/phase_thermal_externalheat.f90

@@ -43,8 +43,8 @@ module function externalheat_init(source_length) result(mySources)
 
   mySources = thermal_active('externalheat',source_length)
   if(count(mySources) == 0) return
-  print'(/,a)', ' <<<+-  phase:thermal:externalheat init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:thermal:externalheat init  -+>>>'
-  print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')
@@ -53,11 +53,11 @@ module function externalheat_init(source_length) result(mySources)
 
   do ph = 1, phases%length
     phase => phases%get(ph)
-    if(count(mySources(:,ph)) == 0) cycle
+    if (count(mySources(:,ph)) == 0) cycle
     thermal => phase%get('thermal')
     sources => thermal%get('source')
     do so = 1, sources%length
-      if(mySources(so,ph)) then
+      if (mySources(so,ph)) then
         source_thermal_externalheat_offset(ph) = so
         associate(prm  => param(ph))
           src => sources%get(so)
@@ -70,9 +70,9 @@ module function externalheat_init(source_length) result(mySources)
           Nmembers = count(material_phaseID == ph)
           call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
         end associate
-      endif
+      end if
-    enddo
+    end do
-  enddo
+  end do
 
 end function externalheat_init
 
@@ -125,7 +125,7 @@ module function externalheat_f_T(ph,en) result(f_T)
         f_T = prm%f_T(interval  ) * (1.0_pReal - frac_time) + &
               prm%f_T(interval+1) * frac_time                                                      ! interpolate heat rate between segment boundaries...
                                                                                                    ! ...or extrapolate if outside of bounds
-    enddo
+    end do
   end associate
 
 end function externalheat_f_T

src/prec.f90

@@ -69,15 +69,15 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine prec_init
 
-  print'(/,a)', ' <<<+-  prec init  -+>>>'
+  print'(/,1x,a)', '<<<+-  prec init  -+>>>'
 
-  print'(a,i3)',    ' Size of integer in bit: ',bit_size(0)
+  print'(/,a,i3)',    ' integer size / bit:   ',bit_size(0)
-  print'(a,i19)',   '   Maximum value:        ',huge(0)
+  print'(  a,i19)',   '   maximum value:      ',huge(0)
-  print'(/,a,i3)',  ' Size of float in bit:   ',storage_size(0.0_pReal)
+  print'(/,a,i3)',    ' float size / bit:     ',storage_size(0.0_pReal)
-  print'(a,e10.3)', '   Maximum value:        ',huge(0.0_pReal)
+  print'(  a,e10.3)', '   maximum value:      ',huge(0.0_pReal)
-  print'(a,e10.3)', '   Minimum value:        ',PREAL_MIN
+  print'(  a,e10.3)', '   minimum value:      ',PREAL_MIN
-  print'(a,e10.3)', '   Epsilon value:        ',PREAL_EPSILON
+  print'(  a,e10.3)', '   epsilon value:      ',PREAL_EPSILON
-  print'(a,i3)',    '   Decimal precision:    ',precision(0.0_pReal)
+  print'(  a,i3)',    '   decimal precision:  ',precision(0.0_pReal)
 
   call selfTest
 

src/results.f90

@@ -70,10 +70,10 @@ subroutine results_init(restart)
   character(len=:), allocatable :: date
 
 
-  print'(/,a)', ' <<<+-  results init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  results init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017'
+  print'(/,1x,a)', 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017'
-  print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1007/s40192-017-0084-5'
 
   if (.not. restart) then
     resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w')
@@ -98,7 +98,7 @@ subroutine results_init(restart)
     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)
-  endif
+  end if
 
   call results_closeJobFile
 
@@ -222,7 +222,7 @@ subroutine results_addAttribute_str(attrLabel,attrValue,path)
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
   else
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
-  endif
+  end if
 
 end subroutine results_addAttribute_str
 
@@ -241,7 +241,7 @@ subroutine results_addAttribute_int(attrLabel,attrValue,path)
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
   else
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
-  endif
+  end if
 
 end subroutine results_addAttribute_int
 
@@ -260,7 +260,7 @@ subroutine results_addAttribute_real(attrLabel,attrValue,path)
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
   else
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
-  endif
+  end if
 
 end subroutine results_addAttribute_real
 
@@ -279,7 +279,7 @@ subroutine results_addAttribute_str_array(attrLabel,attrValue,path)
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
   else
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
-  endif
+  end if
 
 end subroutine results_addAttribute_str_array
 
@@ -298,7 +298,7 @@ subroutine results_addAttribute_int_array(attrLabel,attrValue,path)
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
   else
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
-  endif
+  end if
 
 end subroutine results_addAttribute_int_array
 
@@ -317,7 +317,7 @@ subroutine results_addAttribute_real_array(attrLabel,attrValue,path)
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue, path)
   else
     call HDF5_addAttribute(resultsFile,attrLabel, attrValue)
-  endif
+  end if
 
 end subroutine results_addAttribute_real_array
 
@@ -390,7 +390,7 @@ subroutine results_writeVectorDataset_real(dataset,group,label,description,SIuni
 
   if (present(systems)) then
     if (size(systems)*size(dataset,2) == 0 ) return !ToDo: maybe also implement for other results_write (not sure about scalar)
-  endif
+  end if
 
   groupHandle = results_openGroup(group)
   call HDF5_write(dataset,groupHandle,label)
@@ -422,7 +422,7 @@ subroutine results_writeTensorDataset_real(dataset,group,label,description,SIuni
     transposed_ = transposed
   else
     transposed_ = .true.
-  endif
+  end if
 
   groupHandle = results_openGroup(group)
   if(transposed_) then
@@ -430,11 +430,11 @@ subroutine results_writeTensorDataset_real(dataset,group,label,description,SIuni
     allocate(dataset_transposed,mold=dataset)
     do i=1,size(dataset_transposed,3)
       dataset_transposed(:,:,i) = transpose(dataset(:,:,i))
-    enddo
+    end do
     call HDF5_write(dataset_transposed,groupHandle,label)
   else
     call HDF5_write(dataset,groupHandle,label)
-  endif
+  end if
   call executionStamp(group//'/'//label,description,SIunit)
   call HDF5_closeGroup(groupHandle)
 
@@ -456,7 +456,7 @@ subroutine results_writeVectorDataset_int(dataset,group,label,description,SIunit
 
   if (present(systems)) then
     if (size(systems)*size(dataset,2) == 0 ) return !ToDo: maybe also implement for other results_write (not sure about scalar)
-  endif
+  end if
 
   groupHandle = results_openGroup(group)
   call HDF5_write(dataset,groupHandle,label)
@@ -542,16 +542,16 @@ subroutine results_mapping_phase(ID,entry,label)
   do co = 1, size(ID,1)
     do ce = 1, size(ID,2)
       entryOffset(ID(co,ce),worldrank) = entryOffset(ID(co,ce),worldrank) +1
-    enddo
+    end do
-  enddo
+  end do
   call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,MPI_COMM_WORLD,ierr)! get offset at each process
   if(ierr /= 0) error stop 'MPI error'
   entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2)
   do co = 1, size(ID,1)
     do ce = 1, size(ID,2)
       entryGlobal(co,ce) = entry(co,ce) -1 + entryOffset(ID(co,ce),worldrank)
-    enddo
+    end do
-  enddo
+  end do
 #endif
 
   myShape = int([size(ID,1),writeSize(worldrank)], HSIZE_T)
@@ -694,13 +694,13 @@ subroutine results_mapping_homogenization(ID,entry,label)
   entryOffset = 0
   do ce = 1, size(ID,1)
     entryOffset(ID(ce),worldrank) = entryOffset(ID(ce),worldrank) +1
-  enddo
+  end do
   call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,MPI_COMM_WORLD,ierr)! get offset at each process
   if(ierr /= 0) error stop 'MPI error'
   entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2)
   do ce = 1, size(ID,1)
     entryGlobal(ce) = entry(ce) -1 + entryOffset(ID(ce),worldrank)
-  enddo
+  end do
 #endif
 
   myShape = int([writeSize(worldrank)], HSIZE_T)

src/rotations.f90

@@ -103,10 +103,10 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine rotations_init
 
-  print'(/,a)', ' <<<+-  rotations init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  rotations init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
+  print'(/,1x,a)', 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
-  print*, 'https://doi.org/10.1088/0965-0393/23/8/083501'
+  print'(  1x,a)', 'https://doi.org/10.1088/0965-0393/23/8/083501'
 
   call selfTest