don' hardcode output unit

replaced remaining write(6,..), flush(6) with write(OUTPUT_UNIT,...),
flush(OUTPUT_UNIT)

src/CPFEM.f90

@@ -106,7 +106,7 @@ subroutine CPFEM_init
     num_commercialFEM, &
     debug_CPFEM
 
-  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(6)
+  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(OUTPUT_UNIT)
 
   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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+        flush(OUTPUT_UNIT)
     endif
 
   endif

src/CPFEM2.f90

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

src/DAMASK_marc.f90

@@ -265,9 +265,9 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
     print'(a,i2)',         ' Coordinates:                   ', ncrd
     print'(a,i12)',        ' Nodes:                         ', nnode
     print'(a,i1)',         ' Deformation gradient:          ', itel
-    write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
+    write(OUTPUT_UNIT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
                                   transpose(ffn)
-    write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
+    write(OUTPUT_UNIT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
                                   transpose(ffn1)
   endif
 
@@ -312,7 +312,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
       cycleCounter = -1                                                                             ! first calc step increments this to cycle = 0
       print'(a,i6,1x,i2)', '<< HYPELA2 >> cutback detected..! ',m(1),nn
     endif                                                                                           ! convergence treatment end
-    flush(6)
+    flush(OUTPUT_UNIT)
 
     if (lastLovl /= lovl) then
       cycleCounter  = cycleCounter + 1

src/IO.f90

@@ -52,7 +52,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine IO_init
 
-  print'(/,a)', ' <<<+-  IO init  -+>>>'; flush(6)
+  print'(/,a)', ' <<<+-  IO init  -+>>>'; flush(OUTPUT_UNIT)
 
   call selfTest
 

src/base64.f90

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

src/config.f90

@@ -35,7 +35,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine config_init
 
-  print'(/,a)', ' <<<+-  config init  -+>>>'; flush(6)
+  print'(/,a)', ' <<<+-  config init  -+>>>'; flush(OUTPUT_UNIT)
 
   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(6)
+  print*, 'reading '//fname; flush(OUTPUT_UNIT)
   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(6)
+    print*, 'reading numerics.yaml'; flush(OUTPUT_UNIT)
     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(6)
+    print*, 'reading debug.yaml'; flush(OUTPUT_UNIT)
     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(6)
+  print'(/,a)', ' <<<+-  constitutive init  -+>>>'; flush(OUTPUT_UNIT)
 
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+  print*, ' writing field and constitutive data required for restart to file';flush(OUTPUT_UNIT)
 
   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(6)
+  print'(/,a)', ' <<<+-  damage_local init  -+>>>'; flush(OUTPUT_UNIT)
 
 !----------------------------------------------------------------------------------------------
 ! 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(6)
+  print'(/,a)', ' <<<+-  element_init  -+>>>'; flush(OUTPUT_UNIT)
 
   print*, 'element type:      ',self%elemType
   print*, '  geom type:       ',self%geomType

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(6)
+  print'(/,a)', ' <<<+-  homogenization init  -+>>>'; flush(OUTPUT_UNIT)
 
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
 
   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(6)
+      flush(OUTPUT_UNIT)
     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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+      print*, '... done and happy'; flush(OUTPUT_UNIT)
 #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(6)
+      flush(OUTPUT_UNIT)
     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(6)
+      print'(a,/)', ' ... broken'; flush(OUTPUT_UNIT)
 #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(6)
+     print'(a,/)', ' ... not yet done'; flush(OUTPUT_UNIT)
 #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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+    flush(OUTPUT_UNIT)
     !$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(6)
+    flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
 
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
 
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(/,a)', ' <<<+-  lattice init  -+>>>'; flush(OUTPUT_UNIT)
 
   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(6)
+  print'(/,a)', ' <<<+-  material init  -+>>>'; flush(OUTPUT_UNIT)
 
   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(6)
+  print'(/,a)', ' <<<+-  math init  -+>>>'; flush(OUTPUT_UNIT)
 
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
   randomSeed  = num_generic%get_asInt('random_seed', defaultVal = 0)

src/results.f90

@@ -65,7 +65,7 @@ subroutine results_init(restart)
 
   character(len=pStringLen) :: commandLine
 
-  print'(/,a)', ' <<<+-  results init  -+>>>'; flush(6)
+  print'(/,a)', ' <<<+-  results init  -+>>>'; flush(OUTPUT_UNIT)
 
   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(6)
+  print'(/,a)', ' <<<+-  rotations init  -+>>>'; flush(OUTPUT_UNIT)
 
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   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(6)
+  print'(a,i2)', ' # instances: ',Ninstance; flush(OUTPUT_UNIT)
   if(Ninstance == 0) return
 
   phases => config_material%get('phase')