diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90
index 9eef17719..81e1887c6 100644
--- a/src/constitutive_plastic_nonlocal.f90
+++ b/src/constitutive_plastic_nonlocal.f90
@@ -741,10 +741,10 @@ module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
   if (debugConstitutive%extensive &
       .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
              .or. .not. debugConstitutive%selective)) then
-    write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_microstructure at el ip ',el,ip
-    write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', stt%rho_forest(:,of)
-    write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,of)*1e-6
-    write(6,'(a,/,12x,12(f10.5,1x),/)') '<< CONST >> tauBack / MPa', dst%tau_back(:,of)*1e-6
+    print'(/,a,i8,1x,i2,1x,i1,/)', '<< CONST >> nonlocal_microstructure at el ip ',el,ip
+    print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,of)
+    print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,of)*1e-6
+    print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,of)*1e-6
   endif
 #endif
 
@@ -958,8 +958,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,instance,of,ip,el)
   if (debugConstitutive%extensive &
       .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
              .or. .not. debugConstitutive%selective)) then
-    write(6,'(a,/,8(12x,12(e12.5,1x),/))')    '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8)
-    write(6,'(a,/,10(12x,12(e12.5,1x),/),/)') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress
+    print'(a,/,8(12x,12(e12.5,1x),/))',    '<< CONST >> dislocation remobilization', deltaRhoRemobilization(:,1:8)
+    print'(a,/,10(12x,12(e12.5,1x),/),/)', '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress
   endif
 #endif
 
@@ -1047,8 +1047,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
     if (debugConstitutive%basic &
         .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip) &
                .or. .not. debugConstitutive%selective)) then
-      write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
-      write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot
+      print'(a,/,10(12x,12(e12.5,1x),/))', '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
+      print'(a,/,4(12x,12(e12.5,1x),/))', '<< CONST >> gdot / 1/s',gdot
     endif
 #endif
 
@@ -1156,8 +1156,8 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
      .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then
 #ifdef DEBUG
     if (debugConstitutive%extensive) then
-      write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
-      write(6,'(a)') '<< CONST >> enforcing cutback !!!'
+      print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
+      print'(a)', '<< CONST >> enforcing cutback !!!'
     endif
 #endif
     plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
@@ -1268,8 +1268,8 @@ function rhoDotFlux(F,Fp,timestep,  instance,of,ip,el)
                 > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then                                  ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
 #ifdef DEBUG
     if (debugConstitutive%extensive) then
-      write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
-      write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at  ', &
+      print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
+      print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at  ', &
         maxval(abs(v0), abs(gdot) > 0.0_pReal &
                        .and.  prm%CFLfactor * abs(v0) * timestep &
                              > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90
index 341210afc..213322278 100644
--- a/src/grid/DAMASK_grid.f90
+++ b/src/grid/DAMASK_grid.f90
@@ -341,8 +341,7 @@ program DAMASK_grid
     writeHeader: if (interface_restartInc < 1) then
       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')) &
-        write(6,'(/,a)') ' header of statistics file written out'
+      if (debug_grid%contains('basic')) print'(/,a)', ' header of statistics file written out'
       flush(6)
     else writeHeader
       open(newunit=statUnit,file=trim(getSolverJobName())//&
@@ -351,7 +350,7 @@ program DAMASK_grid
   endif
   
   writeUndeformed: if (interface_restartInc < 1) then
-    write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
+    print'(/,a)', ' ... writing initial configuration to file ........................'
     call CPFEM_results(0,0.0_pReal)
   endif writeUndeformed
   
diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90
index 764e8adb6..5e53fe2a5 100644
--- a/src/grid/spectral_utilities.f90
+++ b/src/grid/spectral_utilities.f90
@@ -700,7 +700,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
     temp99_real = math_3333to99(rot_BC%rotate(C))
 
     if(debugGeneral) then
-      write(6,'(/,a)') ' ... updating masked compliance ............................................'
+      print'(/,a)', ' ... updating masked compliance ............................................'
       write(6,'(/,a,/,9(9(2x,f12.7,1x)/))',advance='no') ' Stiffness C (load) / GPa =',&
                                                    transpose(temp99_Real)*1.0e-9_pReal
       flush(6)