standard conforming name

src/CPFEM.f90

@@ -106,7 +106,7 @@ subroutine CPFEM_init
     num_commercialFEM, &
     debug_CPFEM
 
-  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
 
   allocate(CPFEM_cs(               6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
   allocate(CPFEM_dcsdE(          6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
@@ -132,7 +132,7 @@ subroutine CPFEM_init
     print'(a32,1x,6(i8,1x))',   'CPFEM_cs:              ', shape(CPFEM_cs)
     print'(a32,1x,6(i8,1x))',   'CPFEM_dcsdE:           ', shape(CPFEM_dcsdE)
     print'(a32,1x,6(i8,1x),/)', 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 
 end subroutine CPFEM_init
@@ -250,7 +250,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
           '<< CPFEM >> stress/MPa at elFE ip ',   elFE, ip, CPFEM_cs(1:6,ip,elCP)*1.0e-6_pReal
         print'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))', &
           '<< CPFEM >> Jacobian/GPa at elFE ip ', elFE, ip, transpose(CPFEM_dcsdE(1:6,1:6,ip,elCP))*1.0e-9_pReal
-        flush(OUTPUT_UNIT)
+        flush(IO_STDOUT)
     endif
 
   endif

src/CPFEM2.f90

@@ -76,7 +76,7 @@ end subroutine CPFEM_initAll
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_init
 
-  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
 
   if (interface_restartInc > 0) call crystallite_restartRead
 

src/IO.f90

@@ -7,8 +7,8 @@
 !--------------------------------------------------------------------------------------------------
 module IO
   use, intrinsic :: ISO_fortran_env, only: &
-    OUTPUT_UNIT, &
+    IO_STDOUT => OUTPUT_UNIT, &
-    ERROR_UNIT
+    IO_STDERR => ERROR_UNIT
 
   use prec
 
@@ -20,7 +20,7 @@ module IO
   character, parameter, public :: &
     IO_EOL = new_line('DAMASK'), &                                                                  !< end of line character
     IO_COMMENT = '#'
-  character(len=*), parameter, private :: &
+  character(len=*), parameter :: &
     IO_DIVIDER = '───────────────────'//&
                  '───────────────────'//&
                  '───────────────────'//&
@@ -42,7 +42,7 @@ module IO
     IO_stringAsBool, &
     IO_error, &
     IO_warning, &
-    OUTPUT_UNIT
+    IO_STDOUT
 
 contains
 
@@ -52,7 +52,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine IO_init
 
-  print'(/,a)', ' <<<+-  IO init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  IO init  -+>>>'; flush(IO_STDOUT)
 
   call selfTest
 
@@ -543,29 +543,29 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
   end select
 
   !$OMP CRITICAL (write2out)
-  write(ERROR_UNIT,'(/,a)')                ' ┌'//IO_DIVIDER//'┐'
+  write(IO_STDERR,'(/,a)')                ' ┌'//IO_DIVIDER//'┐'
-  write(ERROR_UNIT,'(a,24x,a,40x,a)')      ' │','error',                                             '│'
+  write(IO_STDERR,'(a,24x,a,40x,a)')      ' │','error',                                             '│'
-  write(ERROR_UNIT,'(a,24x,i3,42x,a)')     ' │',error_ID,                                            '│'
+  write(IO_STDERR,'(a,24x,i3,42x,a)')     ' │',error_ID,                                            '│'
-  write(ERROR_UNIT,'(a)')                  ' ├'//IO_DIVIDER//'┤'
+  write(IO_STDERR,'(a)')                  ' ├'//IO_DIVIDER//'┤'
   write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(msg)),',',&
                                                      max(1,72-len_trim(msg)-4),'x,a)'
-  write(ERROR_UNIT,formatString)            '│ ',trim(msg),                                          '│'
+  write(IO_STDERR,formatString)            '│ ',trim(msg),                                          '│'
   if (present(ext_msg)) then
     write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
                                                        max(1,72-len_trim(ext_msg)-4),'x,a)'
-    write(ERROR_UNIT,formatString)          '│ ',trim(ext_msg),                                      '│'
+    write(IO_STDERR,formatString)          '│ ',trim(ext_msg),                                      '│'
   endif
   if (present(el)) &
-    write(ERROR_UNIT,'(a19,1x,i9,44x,a3)') ' │ at element    ',el,                                   '│'
+    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at element    ',el,                                   '│'
   if (present(ip)) &
-    write(ERROR_UNIT,'(a19,1x,i9,44x,a3)') ' │ at IP         ',ip,                                   '│'
+    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at IP         ',ip,                                   '│'
   if (present(g)) &
-    write(ERROR_UNIT,'(a19,1x,i9,44x,a3)') ' │ at constituent',g,                                    '│'
+    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at constituent',g,                                    '│'
   if (present(instance)) &
-    write(ERROR_UNIT,'(a19,1x,i9,44x,a3)') ' │ at instance   ',instance,                             '│'
+    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at instance   ',instance,                             '│'
-  write(ERROR_UNIT,'(a,69x,a)')            ' │',                                                     '│'
+  write(IO_STDERR,'(a,69x,a)')            ' │',                                                     '│'
-  write(ERROR_UNIT,'(a)')                  ' └'//IO_DIVIDER//'┘'
+  write(IO_STDERR,'(a)')                  ' └'//IO_DIVIDER//'┘'
-  flush(ERROR_UNIT)
+  flush(IO_STDERR)
   call quit(9000+error_ID)
   !$OMP END CRITICAL (write2out)
 
@@ -628,27 +628,27 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
     end select
 
   !$OMP CRITICAL (write2out)
-  write(ERROR_UNIT,'(/,a)')                ' ┌'//IO_DIVIDER//'┐'
+  write(IO_STDERR,'(/,a)')                ' ┌'//IO_DIVIDER//'┐'
-  write(ERROR_UNIT,'(a,24x,a,38x,a)')      ' │','warning',                                           '│'
+  write(IO_STDERR,'(a,24x,a,38x,a)')      ' │','warning',                                           '│'
-  write(ERROR_UNIT,'(a,24x,i3,42x,a)')     ' │',warning_ID,                                          '│'
+  write(IO_STDERR,'(a,24x,i3,42x,a)')     ' │',warning_ID,                                          '│'
-  write(ERROR_UNIT,'(a)')                  ' ├'//IO_DIVIDER//'┤'
+  write(IO_STDERR,'(a)')                  ' ├'//IO_DIVIDER//'┤'
   write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(msg)),',',&
                                                      max(1,72-len_trim(msg)-4),'x,a)'
-  write(ERROR_UNIT,formatString)            '│ ',trim(msg),                                          '│'
+  write(IO_STDERR,formatString)            '│ ',trim(msg),                                          '│'
   if (present(ext_msg)) then
     write(formatString,'(a,i6.6,a,i6.6,a)') '(1x,a4,a',max(1,len_trim(ext_msg)),',',&
                                                        max(1,72-len_trim(ext_msg)-4),'x,a)'
-    write(ERROR_UNIT,formatString)          '│ ',trim(ext_msg),                                      '│'
+    write(IO_STDERR,formatString)          '│ ',trim(ext_msg),                                      '│'
   endif
   if (present(el)) &
-    write(ERROR_UNIT,'(a19,1x,i9,44x,a3)') ' │ at element    ',el,                                   '│'
+    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at element    ',el,                                   '│'
   if (present(ip)) &
-    write(ERROR_UNIT,'(a19,1x,i9,44x,a3)') ' │ at IP         ',ip,                                   '│'
+    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at IP         ',ip,                                   '│'
   if (present(g)) &
-    write(ERROR_UNIT,'(a19,1x,i9,44x,a3)') ' │ at constituent',g,                                    '│'
+    write(IO_STDERR,'(a19,1x,i9,44x,a3)') ' │ at constituent',g,                                    '│'
-  write(ERROR_UNIT,'(a,69x,a)')            ' │',                                                     '│'
+  write(IO_STDERR,'(a,69x,a)')            ' │',                                                     '│'
-  write(ERROR_UNIT,'(a)')                  ' └'//IO_DIVIDER//'┘'
+  write(IO_STDERR,'(a)')                  ' └'//IO_DIVIDER//'┘'
-  flush(ERROR_UNIT)
+  flush(IO_STDERR)
   !$OMP END CRITICAL (write2out)
 
 end subroutine IO_warning

src/base64.f90

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

src/config.f90

@@ -35,7 +35,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine config_init
 
-  print'(/,a)', ' <<<+-  config init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  config init  -+>>>'; flush(IO_STDOUT)
 
   call parse_material
   call parse_numerics
@@ -59,7 +59,7 @@ subroutine parse_material
     inquire(file=fname,exist=fileExists)
     if(.not. fileExists) call IO_error(100,ext_msg=fname)
   endif
-  print*, 'reading '//fname; flush(OUTPUT_UNIT)
+  print*, 'reading '//fname; flush(IO_STDOUT)
   config_material => YAML_parse_file(fname)
 
 end subroutine parse_material
@@ -75,7 +75,7 @@ subroutine parse_numerics
   config_numerics => emptyDict
   inquire(file='numerics.yaml', exist=fexist)
   if (fexist) then
-    print*, 'reading numerics.yaml'; flush(OUTPUT_UNIT)
+    print*, 'reading numerics.yaml'; flush(IO_STDOUT)
     config_numerics => YAML_parse_file('numerics.yaml')
   endif
 
@@ -92,7 +92,7 @@ subroutine parse_debug
   config_debug => emptyDict
   inquire(file='debug.yaml', exist=fexist)
   fileExists: if (fexist) then
-    print*, 'reading debug.yaml'; flush(OUTPUT_UNIT)
+    print*, 'reading debug.yaml'; flush(IO_STDOUT)
     config_debug => YAML_parse_file('debug.yaml')
   endif fileExists
 

src/constitutive.f90

@@ -446,7 +446,7 @@ subroutine constitutive_init
   call damage_init
   call thermal_init
 
-  print'(/,a)', ' <<<+-  constitutive init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  constitutive init  -+>>>'; flush(IO_STDOUT)
 
   constitutive_source_maxSizeDotState = 0
   PhaseLoop2:do p = 1,phases%length

src/constitutive_plastic_disloTungsten.f90

@@ -100,7 +100,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
 
   myPlasticity = plastic_active('disloTungsten')
   Ninstance = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
   
   print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016'

src/constitutive_plastic_dislotwin.f90

@@ -147,7 +147,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
   myPlasticity = plastic_active('dislotwin')
   Ninstance = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
   
   print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'

src/constitutive_plastic_isotropic.f90

@@ -71,7 +71,7 @@ module function plastic_isotropic_init() result(myPlasticity)
 
   myPlasticity = plastic_active('isotropic')
   Ninstance = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'

src/constitutive_plastic_kinehardening.f90

@@ -83,7 +83,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
 
   myPlasticity = plastic_active('kinehardening')
   Ninstance = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   allocate(param(Ninstance))

src/constitutive_plastic_none.f90

@@ -35,7 +35,7 @@ module function plastic_none_init() result(myPlasticity)
   enddo
 
   Ninstance = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   do p = 1, phases%length

src/constitutive_plastic_nonlocal.f90

@@ -189,7 +189,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
 
   myPlasticity = plastic_active('nonlocal')
   Ninstance = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) then
     call geometry_plastic_nonlocal_disable
     return

src/constitutive_plastic_phenopowerlaw.f90

@@ -92,7 +92,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 
   myPlasticity = plastic_active('phenopowerlaw')
   Ninstance = count(myPlasticity)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   allocate(param(Ninstance))

src/crystallite.f90

@@ -294,7 +294,7 @@ subroutine crystallite_init
     print'(a42,1x,i10)', '    # of elements:                       ', eMax
     print'(a42,1x,i10)', '    # of integration points/element:     ', iMax
     print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 
 #endif
@@ -1561,7 +1561,7 @@ subroutine crystallite_restartWrite
   integer(HID_T) :: fileHandle, groupHandle
   character(len=pStringLen) :: fileName, datasetName
 
-  print*, ' writing field and constitutive data required for restart to file';flush(OUTPUT_UNIT)
+  print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT)
 
   write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
   fileHandle = HDF5_openFile(fileName,'a')

src/damage_local.f90

@@ -49,7 +49,7 @@ subroutine damage_local_init
     homog, &
     homogDamage
 
-  print'(/,a)', ' <<<+-  damage_local init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  damage_local init  -+>>>'; flush(IO_STDOUT)
 
 !----------------------------------------------------------------------------------------------
 ! read numerics parameter and do sanity check

src/element.f90

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

src/grid/DAMASK_grid.f90

@@ -99,7 +99,7 @@ program DAMASK_grid
 ! init DAMASK (all modules)
  
   call CPFEM_initAll 
-  print'(/,a)',   ' <<<+-  DAMASK_spectral init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)',   ' <<<+-  DAMASK_spectral init  -+>>>'; flush(IO_STDOUT)
   
   print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
   print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
@@ -279,11 +279,11 @@ program DAMASK_grid
       endif
       do i = 1, 3; do j = 1, 3
         if(newLoadCase%deformation%maskLogical(i,j)) then
-          write(OUTPUT_UNIT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
+          write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j)
         else
-          write(OUTPUT_UNIT,'(2x,12a)',advance='no') '     *      '
+          write(IO_STDOUT,'(2x,12a)',advance='no') '     *      '
           endif
-        enddo; write(OUTPUT_UNIT,'(/)',advance='no')
+        enddo; write(IO_STDOUT,'(/)',advance='no')
       enddo
       if (any(newLoadCase%stress%maskLogical .eqv. &
               newLoadCase%deformation%maskLogical)) errorID = 831                                   ! exclusive or masking only
@@ -292,17 +292,17 @@ program DAMASK_grid
       print*, ' stress / GPa:'
       do i = 1, 3; do j = 1, 3
         if(newLoadCase%stress%maskLogical(i,j)) then
-          write(OUTPUT_UNIT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
+          write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal
         else
-          write(OUTPUT_UNIT,'(2x,12a)',advance='no') '     *      '
+          write(IO_STDOUT,'(2x,12a)',advance='no') '     *      '
         endif
-        enddo; write(OUTPUT_UNIT,'(/)',advance='no')
+        enddo; write(IO_STDOUT,'(/)',advance='no')
       enddo
       if (any(abs(matmul(newLoadCase%rot%asMatrix(), &
                          transpose(newLoadCase%rot%asMatrix()))-math_I3) > &
                  reshape(spread(tol_math_check,1,9),[ 3,3]))) errorID = 846                         ! given rotation matrix contains strain
       if (any(dNeq(newLoadCase%rot%asMatrix(), math_I3))) &
-        write(OUTPUT_UNIT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
+        write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
                  transpose(newLoadCase%rot%asMatrix())
       if (newLoadCase%time < 0.0_pReal) errorID = 834                                               ! negative time increment
       print'(a,f0.3)', ' time: ', newLoadCase%time
@@ -342,7 +342,7 @@ program DAMASK_grid
       open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
       write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded'                ! statistics file
       if (debug_grid%contains('basic')) print'(/,a)', ' header of statistics file written out'
-      flush(OUTPUT_UNIT)
+      flush(IO_STDOUT)
     else writeHeader
       open(newunit=statUnit,file=trim(getSolverJobName())//&
                                   '.sta',form='FORMATTED', position='APPEND', status='OLD')
@@ -405,7 +405,7 @@ program DAMASK_grid
           write(incInfo,'(4(a,i0))') &
                   'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
                   '-', stepFraction,'/',subStepFactor**cutBackLevel
-          flush(OUTPUT_UNIT)
+          flush(IO_STDOUT)
 
 !--------------------------------------------------------------------------------------------------
 ! forward fields
@@ -489,11 +489,11 @@ program DAMASK_grid
           print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged'
         else
           print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged'
-        endif; flush(OUTPUT_UNIT)
+        endif; flush(IO_STDOUT)
   
         if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then                          ! at output frequency
           print'(1/,a)', ' ... writing results to file ......................................'
-          flush(OUTPUT_UNIT)
+          flush(IO_STDOUT)
           call CPFEM_results(totalIncsCounter,time)
         endif
         if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then

src/grid/discretization_grid.f90

@@ -65,7 +65,7 @@ subroutine discretization_grid_init(restart)
   integer(C_INTPTR_T) :: &
     devNull, z, z_offset
 
-  print'(/,a)', ' <<<+-  discretization_grid init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  discretization_grid init  -+>>>'; flush(IO_STDOUT)
 
   if(index(interface_geomFile,'.vtr') /= 0) then
     call readVTR(grid,geomSize,origin,microstructureAt)

src/grid/grid_damage_spectral.f90

@@ -205,10 +205,10 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution)
   call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
   if (solution%converged) &
     print'(/,a)', ' ... nonlocal damage converged .....................................'
-  write(OUTPUT_UNIT,'(/,a,f8.6,2x,f8.6,2x,e11.4,/)',advance='no') ' Minimum|Maximum|Delta Damage      = ',&
+  write(IO_STDOUT,'(/,a,f8.6,2x,f8.6,2x,e11.4,/)',advance='no') ' Minimum|Maximum|Delta Damage      = ',&
                                                           phi_min, phi_max, stagNorm
   print'(/,a)', ' ==========================================================================='
-  flush(OUTPUT_UNIT) 
+  flush(IO_STDOUT) 
 
 end function grid_damage_spectral_solution
 

src/grid/grid_mech_FEM.f90

@@ -122,7 +122,7 @@ subroutine grid_mech_FEM_init
   PetscScalar, pointer, dimension(:,:,:,:) :: &
   u_current,u_lastInc
 
-  print'(/,a)', ' <<<+-  grid_mech_FEM init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  grid_mech_FEM init  -+>>>'; flush(IO_STDOUT)
 
 !-----------------------------------------------------------------------------------------------
 ! debugging options
@@ -413,7 +413,7 @@ subroutine grid_mech_FEM_restartWrite
   call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
   call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
 
-  print*, 'writing solver data required for restart to file'; flush(OUTPUT_UNIT)
+  print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
 
   write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
   fileHandle  = HDF5_openFile(fileName,'w')
@@ -481,7 +481,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
   print'(a,f12.2,a,es8.2,a,es9.2,a)',    ' error stress BC  = ', &
           err_BC/BCTol,    ' (',err_BC, ' Pa,  tol = ',BCTol,')'
   print'(/,a)', ' ==========================================================================='
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
 end subroutine converged
 
@@ -517,11 +517,11 @@ subroutine formResidual(da_local,x_local, &
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax
     if (debugRotation) &
-      write(OUTPUT_UNIT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+      write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    write(OUTPUT_UNIT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+    write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif newIteration
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_mech_spectral_basic.f90

@@ -109,7 +109,7 @@ subroutine grid_mech_spectral_basic_init
   character(len=pStringLen) :: &
     fileName
 
-  print'(/,a)', ' <<<+-  grid_mech_spectral_basic init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  grid_mech_spectral_basic init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
   print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
@@ -375,7 +375,7 @@ subroutine grid_mech_spectral_basic_restartWrite
 
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
 
-  print'(a)', ' writing solver data required for restart to file'; flush(OUTPUT_UNIT)
+  print'(a)', ' writing solver data required for restart to file'; flush(IO_STDOUT)
 
   write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
   fileHandle  = HDF5_openFile(fileName,'w')
@@ -441,7 +441,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
   print'(a,f12.2,a,es8.2,a,es9.2,a)',    ' error stress BC  = ', &
           err_BC/BCTol,    ' (',err_BC, ' Pa,  tol = ',BCTol,')'
   print'(/,a)', ' ==========================================================================='
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
 end subroutine converged
 
@@ -475,11 +475,11 @@ subroutine formResidual(in, F, &
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
     if (debugRotation) &
-      write(OUTPUT_UNIT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+      write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    write(OUTPUT_UNIT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+    write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif newIteration
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_mech_spectral_polarisation.f90

@@ -123,7 +123,7 @@ subroutine grid_mech_spectral_polarisation_init
   character(len=pStringLen) :: &
     fileName
 
-  print'(/,a)', ' <<<+-  grid_mech_spectral_polarisation init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  grid_mech_spectral_polarisation init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
   print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
@@ -433,7 +433,7 @@ subroutine grid_mech_spectral_polarisation_restartWrite
   F     => FandF_tau(0: 8,:,:,:)
   F_tau => FandF_tau(9:17,:,:,:)
 
-  print*, 'writing solver data required for restart to file'; flush(OUTPUT_UNIT)
+  print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
 
   write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
   fileHandle  = HDF5_openFile(fileName,'w')
@@ -505,7 +505,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
   print  '(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC         = ', &
             err_BC/BCTol,    ' (',err_BC,  ' Pa,  tol = ',BCTol,')'
   print'(/,a)', ' ==========================================================================='
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
 end subroutine converged
 
@@ -559,11 +559,11 @@ subroutine formResidual(in, FandF_tau, &
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
     if(debugRotation) &
-      write(OUTPUT_UNIT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+      write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    write(OUTPUT_UNIT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
+    write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif newIteration
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_thermal_spectral.f90

@@ -199,10 +199,10 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution)
   call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr)
   if (solution%converged) &
     print'(/,a)', ' ... thermal conduction converged ..................................'
-  write(OUTPUT_UNIT,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',&
+  write(IO_STDOUT,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',&
                                                         T_min, T_max, stagNorm
   print'(/,a)', ' ==========================================================================='
-  flush(OUTPUT_UNIT) 
+  flush(IO_STDOUT) 
 
 end function grid_thermal_spectral_solution
 

src/grid/spectral_utilities.f90

@@ -211,7 +211,7 @@ subroutine spectral_utilities_init
   if(debugPETSc) print'(3(/,a),/)', &
                  ' Initializing PETSc with debug options: ', &
                  trim(PETScDebug), &
-                 ' add more using the PETSc_Options keyword in numerics.yaml '; flush(OUTPUT_UNIT)
+                 ' add more using the PETSc_Options keyword in numerics.yaml '; flush(IO_STDOUT)
 
   num_grid => config_numerics%get('grid',defaultVal=emptyDict)
 
@@ -280,7 +280,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(OUTPUT_UNIT)
+  print*, 'FFTW initialized'; flush(IO_STDOUT)
 
 !--------------------------------------------------------------------------------------------------
 ! MPI allocation
@@ -507,7 +507,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
 
 
   print'(/,a)', ' ... doing gamma convolution ...............................................'
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
 !--------------------------------------------------------------------------------------------------
 ! do the actual spectral method calculation (mechanical equilibrium)
@@ -577,7 +577,7 @@ real(pReal) function utilities_divergenceRMS()
   complex(pReal), dimension(3)   :: rescaledGeom
 
   print'(/,a)', ' ... calculating divergence ................................................'
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
   rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
 
@@ -621,7 +621,7 @@ real(pReal) function utilities_curlRMS()
   complex(pReal), dimension(3)   :: rescaledGeom
 
   print'(/,a)', ' ... calculating curl ......................................................'
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
   rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
 
@@ -701,9 +701,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
 
     if(debugGeneral) then
       print'(/,a)', ' ... updating masked compliance ............................................'
-      write(OUTPUT_UNIT,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',&
+      write(IO_STDOUT,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',&
                                                    transpose(temp99_Real)*1.0e-9_pReal
-      flush(OUTPUT_UNIT)
+      flush(IO_STDOUT)
     endif
 
     do i = 1,9; do j = 1,9
@@ -723,9 +723,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
     if (debugGeneral .or. errmatinv) then
       write(formatString, '(i2)') size_reduced
       formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
-      write(OUTPUT_UNIT,trim(formatString),advance='no') ' C * S (load) ', &
+      write(IO_STDOUT,trim(formatString),advance='no') ' C * S (load) ', &
                                                              transpose(matmul(c_reduced,s_reduced))
-      write(OUTPUT_UNIT,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced)
+      write(IO_STDOUT,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced)
       if(errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
     endif
     temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9])
@@ -736,9 +736,9 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
   utilities_maskedCompliance = math_99to3333(temp99_Real)
 
   if(debugGeneral) then
-    write(OUTPUT_UNIT,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
+    write(IO_STDOUT,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
       ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 
 end function utilities_maskedCompliance
@@ -823,7 +823,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   real(pReal), dimension(2) :: valueAndRank                                                         !< pair of min/max norm of dPdF to synchronize min/max of dPdF
 
   print'(/,a)', ' ... evaluating constitutive response ......................................'
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
   materialpoint_F  = reshape(F,[3,3,1,product(grid(1:2))*grid3])                                    ! set materialpoint target F to estimated field
 
@@ -833,13 +833,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt                                                   ! average of P
   call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
   if (debugRotation) &
-  write(OUTPUT_UNIT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
+  write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
                                                       transpose(P_av)*1.e-6_pReal
   if(present(rotation_BC)) &
     P_av = rotation_BC%rotate(P_av)
-  write(OUTPUT_UNIT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress       / MPa =',&
+  write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress       / MPa =',&
                                                       transpose(P_av)*1.e-6_pReal
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
   dPdF_max = 0.0_pReal
   dPdF_norm_max = 0.0_pReal
@@ -1095,7 +1095,7 @@ subroutine utilities_saveReferenceStiffness
     fileUnit,ierr
 
   if (worldrank == 0) then
-    print'(a)', ' writing reference stiffness data required for restart to file'; flush(OUTPUT_UNIT)
+    print'(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')

src/homogenization.f90

@@ -186,7 +186,7 @@ subroutine homogenization_init
   materialpoint_F = materialpoint_F0                                                                ! initialize to identity
   allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem),              source=0.0_pReal)
 
-  print'(/,a)', ' <<<+-  homogenization init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  homogenization init  -+>>>'; flush(IO_STDOUT)
 
   num%nMPstate          = num_homogGeneric%get_asInt  ('nMPstate',     defaultVal=10)
   num%subStepMinHomog   = num_homogGeneric%get_asFloat('subStepMin',   defaultVal=1.0e-3_pReal)

src/homogenization_mech_RGC.f90

@@ -95,7 +95,7 @@ module subroutine mech_RGC_init(num_homogMech)
   print'(/,a)', ' <<<+-  homogenization_mech_rgc init  -+>>>'
 
   Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
 
   print*, '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
@@ -247,7 +247,7 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
         print'(1x,3(e15.8,1x))',(F(i,j,iGrain), j = 1,3)
       enddo
       print*,' '
-      flush(OUTPUT_UNIT)
+      flush(IO_STDOUT)
     endif
 #endif
   enddo
@@ -376,7 +376,7 @@ module procedure mech_RGC_updateState
                '@ grain ',stresLoc(3),' in component ',stresLoc(1),stresLoc(2)
     print'(a,e15.8,a,i3,a,i2)',' Max residual: ',residMax, &
                ' @ iface ',residLoc(1),' in direction ',residLoc(2)
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 #endif
 
@@ -388,7 +388,7 @@ module procedure mech_RGC_updateState
     mech_RGC_updateState = .true.
 #ifdef DEBUG
     if (debugHomog%extensive .and. prm%of_debug == of) &
-      print*, '... done and happy'; flush(OUTPUT_UNIT)
+      print*, '... done and happy'; flush(IO_STDOUT)
 #endif
 
 !--------------------------------------------------------------------------------------------------
@@ -416,7 +416,7 @@ module procedure mech_RGC_updateState
       print'(a,e15.8,/)', ' Volume discrepancy: ',      dst%volumeDiscrepancy(of)
       print'(a,e15.8)',   ' Maximum relaxation rate: ', dst%relaxationRate_max(of)
       print'(a,e15.8,/)', ' Average relaxation rate: ', dst%relaxationRate_avg(of)
-      flush(OUTPUT_UNIT)
+      flush(IO_STDOUT)
     endif
 #endif
 
@@ -429,7 +429,7 @@ module procedure mech_RGC_updateState
 
 #ifdef DEBUG
     if (debugHomog%extensive .and. prm%of_debug == of) &
-      print'(a,/)', ' ... broken'; flush(OUTPUT_UNIT)
+      print'(a,/)', ' ... broken'; flush(IO_STDOUT)
 #endif
 
    return
@@ -437,7 +437,7 @@ module procedure mech_RGC_updateState
  else                                                                                               ! proceed with computing the Jacobian and state update
 #ifdef DEBUG
    if (debugHomog%extensive .and. prm%of_debug == of) &
-     print'(a,/)', ' ... not yet done'; flush(OUTPUT_UNIT)
+     print'(a,/)', ' ... not yet done'; flush(IO_STDOUT)
 #endif
 
  endif
@@ -499,7 +499,7 @@ module procedure mech_RGC_updateState
       print'(1x,100(e11.4,1x))',(smatrix(i,j), j = 1,3*nIntFaceTot)
     enddo
     print*,' '
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 #endif
 
@@ -559,7 +559,7 @@ module procedure mech_RGC_updateState
       print'(1x,100(e11.4,1x))',(pmatrix(i,j), j = 1,3*nIntFaceTot)
     enddo
     print*,' '
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 #endif
 
@@ -578,7 +578,7 @@ module procedure mech_RGC_updateState
       print'(1x,100(e11.4,1x))',(rmatrix(i,j), j = 1,3*nIntFaceTot)
     enddo
     print*,' '
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 #endif
 
@@ -593,7 +593,7 @@ module procedure mech_RGC_updateState
       print'(1x,100(e11.4,1x))',(jmatrix(i,j), j = 1,3*nIntFaceTot)
     enddo
     print*,' '
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 #endif
 
@@ -609,7 +609,7 @@ module procedure mech_RGC_updateState
       print'(1x,100(e11.4,1x))',(jnverse(i,j), j = 1,3*nIntFaceTot)
     enddo
     print*,' '
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 #endif
 
@@ -625,7 +625,7 @@ module procedure mech_RGC_updateState
     !$OMP CRITICAL (write2out)
     print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback'
     print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
     !$OMP END CRITICAL (write2out)
   endif
 
@@ -636,7 +636,7 @@ module procedure mech_RGC_updateState
       print'(1x,2(e15.8,1x))', stt%relaxationVector(i,of)
     enddo
     print*,' '
-    flush(OUTPUT_UNIT)
+    flush(IO_STDOUT)
   endif
 #endif
 

src/homogenization_mech_isostrain.f90

@@ -40,7 +40,7 @@ module subroutine mech_isostrain_init
   print'(/,a)',   ' <<<+-  homogenization_mech_isostrain init  -+>>>'
 
   Ninstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
 
   allocate(param(Ninstance))                                                                        ! one container of parameters per instance
 

src/homogenization_mech_none.f90

@@ -21,7 +21,7 @@ module subroutine mech_none_init
   print'(/,a)', ' <<<+-  homogenization_mech_none init  -+>>>'
 
   Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
 
   do h = 1, size(homogenization_type)
     if (homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle

src/kinematics_cleavage_opening.f90

@@ -49,7 +49,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
 
   myKinematics = kinematics_active('cleavage_opening',kinematics_length)
   Ninstance = count(myKinematics)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')

src/kinematics_slipplane_opening.f90

@@ -52,7 +52,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
 
   myKinematics = kinematics_active('slipplane_opening',kinematics_length)
   Ninstance = count(myKinematics)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')

src/kinematics_thermal_expansion.f90

@@ -42,7 +42,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
 
   myKinematics = kinematics_active('thermal_expansion',kinematics_length)
   Ninstance = count(myKinematics)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')

src/lattice.f90

@@ -457,7 +457,7 @@ subroutine lattice_init
     phase, &
     elasticity
  
-  print'(/,a)', ' <<<+-  lattice init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  lattice init  -+>>>'; flush(IO_STDOUT)
 
   phases => config_material%get('phase')
   Nphases = phases%length

src/material.f90

@@ -164,7 +164,7 @@ subroutine material_init(restart)
     material_homogenization
   character(len=pStringLen) :: sectionName
 
-  print'(/,a)', ' <<<+-  material init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  material init  -+>>>'; flush(IO_STDOUT)
 
   phases => config_material%get('phase')
   allocate(material_name_phase(phases%length))

src/math.f90

@@ -91,7 +91,7 @@ subroutine math_init
   class(tNode), pointer :: &
     num_generic
  
-  print'(/,a)', ' <<<+-  math init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  math init  -+>>>'; flush(IO_STDOUT)
 
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
   randomSeed  = num_generic%get_asInt('random_seed', defaultVal = 0)

src/mesh/DAMASK_mesh.f90

@@ -78,7 +78,7 @@ program DAMASK_mesh
 !--------------------------------------------------------------------------------------------------
 ! init DAMASK (all modules)
   call CPFEM_initAll
-  print'(/,a)', ' <<<+-  DAMASK_mesh init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  DAMASK_mesh init  -+>>>'; flush(IO_STDOUT)
 
 !--------------------------------------------------------------------- 
 ! reading field information from numerics file and do sanity checks
@@ -299,7 +299,7 @@ program DAMASK_mesh
         write(incInfo,'(4(a,i0))') &
                'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
                '-',stepFraction, '/', subStepFactor**cutBackLevel
-        flush(OUTPUT_UNIT)
+        flush(IO_STDOUT)
 
 !--------------------------------------------------------------------------------------------------
 ! forward fields
@@ -363,7 +363,7 @@ program DAMASK_mesh
         print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged'
       else
         print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged'
-      endif; flush(OUTPUT_UNIT)
+      endif; flush(IO_STDOUT)
 
       if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then                            ! at output frequency
         print'(/,a)', ' ... writing results to file ......................................'

src/mesh/FEM_utilities.f90

@@ -122,7 +122,7 @@ subroutine FEM_utilities_init
                  ' Initializing PETSc with debug options: ', &
                  trim(PETScDebug), &
                  ' add more using the PETSc_Options keyword in numerics.yaml '
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
   call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
   CHKERRQ(ierr)
   if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)

src/mesh/mesh_mech_FEM.f90

@@ -110,7 +110,7 @@ subroutine FEM_mech_init(fieldBC)
   class(tNode), pointer :: &
     num_mesh
    
-  print'(/,a)', ' <<<+-  FEM_mech init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  FEM_mech init  -+>>>'; flush(IO_STDOUT)
 
 !-----------------------------------------------------------------------------
 ! read numerical parametes and do sanity checks
@@ -319,7 +319,7 @@ type(tSolutionState) function FEM_mech_solution( &
   endif
 
   print'(/,a)', ' ==========================================================================='
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
 end function FEM_mech_solution
 
@@ -682,9 +682,9 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm
   print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), &
                   ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
                                                   int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
-  write(OUTPUT_UNIT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
+  write(IO_STDOUT,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
                                                         transpose(P_av)*1.e-6_pReal
-  flush(OUTPUT_UNIT)
+  flush(IO_STDOUT)
 
 end subroutine FEM_mech_converged
 

src/results.f90

@@ -65,7 +65,7 @@ subroutine results_init(restart)
 
   character(len=pStringLen) :: commandLine
 
-  print'(/,a)', ' <<<+-  results init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  results init  -+>>>'; flush(IO_STDOUT)
 
   print*, '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

src/rotations.f90

@@ -104,7 +104,7 @@ contains
 subroutine rotations_init
 
   call quaternions_init
-  print'(/,a)', ' <<<+-  rotations init  -+>>>'; flush(OUTPUT_UNIT)
+  print'(/,a)', ' <<<+-  rotations init  -+>>>'; flush(IO_STDOUT)
 
   print*, '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'

src/source_damage_anisoBrittle.f90

@@ -53,7 +53,7 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
 
   mySources = source_active('damage_anisoBrittle',source_length)
   Ninstance = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')

src/source_damage_anisoDuctile.f90

@@ -47,7 +47,7 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
 
   mySources = source_active('damage_anisoDuctile',source_length)
   Ninstance = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')

src/source_damage_isoBrittle.f90

@@ -43,7 +43,7 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
 
   mySources = source_active('damage_isoBrittle',source_length)
   Ninstance = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')

src/source_damage_isoDuctile.f90

@@ -45,7 +45,7 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
 
   mySources = source_active('damage_isoDuctile',source_length)
   Ninstance = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')

src/source_thermal_dissipation.f90

@@ -41,7 +41,7 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
 
   mySources = source_active('thermal_dissipation',source_length)
   Ninstance = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')

src/source_thermal_externalheat.f90

@@ -45,7 +45,7 @@ module function source_thermal_externalheat_init(source_length) result(mySources
 
   mySources = source_active('thermal_externalheat',source_length)
   Ninstance = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')