Merge branch 'Fortran-cleaning' into 'development'

consistent variable name

See merge request damask/DAMASK!372

CMakeLists.txt

@@ -4,7 +4,7 @@ include (FindPkgConfig REQUIRED)
 # Dummy project to determine compiler names and version
 project (Prerequisites LANGUAGES)
 set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig")
-pkg_check_modules (PETSC REQUIRED PETSc>=3.12.0 PETSc<3.15.0)
+pkg_check_modules (PETSC REQUIRED PETSc>=3.12.0 PETSc<3.16.0)
 pkg_get_variable (CMAKE_Fortran_COMPILER PETSc fcompiler)
 pkg_get_variable (CMAKE_C_COMPILER PETSc ccompiler)
 

src/DAMASK_interface.f90

@@ -51,9 +51,7 @@ subroutine DAMASK_interface_init
 #include <petsc/finclude/petscsys.h>
 
 #if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR<PETSC_MINOR_MIN || PETSC_VERSION_MINOR>PETSC_MINOR_MAX
-===================================================================================================
+--  UNSUPPORTED PETSc VERSION --- UNSUPPORTED PETSc VERSION --- UNSUPPORTED PETSc VERSION ---
---  WRONG PETSc VERSION --- WRONG PETSc VERSION --- WRONG PETSc VERSION ---  WRONG PETSc VERSION --
-===================================================================================================
 #endif
 
   character(len=pPathLen*3+pStringLen) :: &

src/phase.f90

@@ -161,13 +161,13 @@ module phase
       real(pReal), dimension(3,3) :: P
     end function phase_P
 
-    module function thermal_T(ph,me) result(T)
+    module function thermal_T(ph,en) result(T)
-      integer, intent(in) :: ph,me
+      integer, intent(in) :: ph,en
       real(pReal) :: T
     end function thermal_T
 
-    module function thermal_dot_T(ph,me) result(dot_T)
+    module function thermal_dot_T(ph,en) result(dot_T)
-      integer, intent(in) :: ph,me
+      integer, intent(in) :: ph,en
       real(pReal) :: dot_T
     end function thermal_dot_T
 
@@ -216,10 +216,10 @@ module phase
 
 ! == cleaned:end ===================================================================================
 
-    module function thermal_stress(Delta_t,ph,me) result(converged_)
+    module function thermal_stress(Delta_t,ph,en) result(converged_)
 
       real(pReal), intent(in) :: Delta_t
-      integer, intent(in) :: ph, me
+      integer, intent(in) :: ph, en
       logical :: converged_
 
     end function thermal_stress
@@ -253,8 +253,8 @@ module phase
         f
     end function phase_f_phi
 
-    module function phase_f_T(ph,me) result(f)
+    module function phase_f_T(ph,en) result(f)
-      integer, intent(in) :: ph, me
+      integer, intent(in) :: ph, en
       real(pReal) :: f
     end function phase_f_T
 

src/phase_mechanical_plastic.f90

@@ -37,7 +37,7 @@ submodule(phase:mechanical) plastic
         myPlasticity
     end function plastic_nonlocal_init
 
-    module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+    module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
@@ -46,10 +46,10 @@ submodule(phase:mechanical) plastic
         Mp
       integer,                         intent(in) :: &
         ph, &
-        me
+        en
     end subroutine isotropic_LpAndItsTangent
 
-    pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+    pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
@@ -58,10 +58,10 @@ submodule(phase:mechanical) plastic
         Mp
       integer,                         intent(in) :: &
         ph, &
-        me
+        en
     end subroutine phenopowerlaw_LpAndItsTangent
 
-    pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+    pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
@@ -70,10 +70,10 @@ submodule(phase:mechanical) plastic
         Mp
       integer,                         intent(in) :: &
         ph, &
-        me
+        en
     end subroutine kinehardening_LpAndItsTangent
 
-    module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+    module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
@@ -85,10 +85,10 @@ submodule(phase:mechanical) plastic
         T
       integer,                         intent(in) :: &
         ph, &
-        me
+        en
     end subroutine dislotwin_LpAndItsTangent
 
-    pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+    pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
@@ -100,11 +100,11 @@ submodule(phase:mechanical) plastic
         T
       integer,                         intent(in) :: &
         ph, &
-        me
+        en
     end subroutine dislotungsten_LpAndItsTangent
 
     module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
-                                                       Mp,Temperature,ph,me)
+                                                       Mp,Temperature,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Lp
       real(pReal), dimension(3,3,3,3), intent(out) :: &
@@ -116,55 +116,55 @@ submodule(phase:mechanical) plastic
         Temperature
       integer,                         intent(in) :: &
         ph, &
-        me
+        en
     end subroutine nonlocal_LpAndItsTangent
 
 
-    module subroutine isotropic_dotState(Mp,ph,me)
+    module subroutine isotropic_dotState(Mp,ph,en)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       integer,                      intent(in) :: &
         ph, &
-        me
+        en
     end subroutine isotropic_dotState
 
-    module subroutine phenopowerlaw_dotState(Mp,ph,me)
+    module subroutine phenopowerlaw_dotState(Mp,ph,en)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       integer,                      intent(in) :: &
         ph, &
-        me
+        en
     end subroutine phenopowerlaw_dotState
 
-    module subroutine plastic_kinehardening_dotState(Mp,ph,me)
+    module subroutine plastic_kinehardening_dotState(Mp,ph,en)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       integer,                      intent(in) :: &
         ph, &
-        me
+        en
     end subroutine plastic_kinehardening_dotState
 
-    module subroutine dislotwin_dotState(Mp,T,ph,me)
+    module subroutine dislotwin_dotState(Mp,T,ph,en)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       real(pReal),                  intent(in) :: &
         T
       integer,                      intent(in) :: &
         ph, &
-        me
+        en
     end subroutine dislotwin_dotState
 
-    module subroutine dislotungsten_dotState(Mp,T,ph,me)
+    module subroutine dislotungsten_dotState(Mp,T,ph,en)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       real(pReal),                  intent(in) :: &
         T
       integer,                      intent(in) :: &
         ph, &
-        me
+        en
     end subroutine dislotungsten_dotState
 
-    module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,me,ip,el)
+    module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,en,ip,el)
       real(pReal), dimension(3,3), intent(in) :: &
         Mp                                                                                          !< MandelStress
       real(pReal), intent(in) :: &
@@ -172,47 +172,47 @@ submodule(phase:mechanical) plastic
         timestep                                                                                    !< substepped crystallite time increment
       integer, intent(in) :: &
         ph, &
-        me, &
+        en, &
         ip, &                                                                                       !< current integration point
         el                                                                                          !< current element number
     end subroutine nonlocal_dotState
 
-    module subroutine dislotwin_dependentState(T,ph,me)
+    module subroutine dislotwin_dependentState(T,ph,en)
       integer,       intent(in) :: &
         ph, &
-        me
+        en
       real(pReal),   intent(in) :: &
         T
     end subroutine dislotwin_dependentState
 
-    module subroutine dislotungsten_dependentState(ph,me)
+    module subroutine dislotungsten_dependentState(ph,en)
       integer,       intent(in) :: &
         ph, &
-        me
+        en
     end subroutine dislotungsten_dependentState
 
-    module subroutine nonlocal_dependentState(ph, me, ip, el)
+    module subroutine nonlocal_dependentState(ph, en, ip, el)
       integer, intent(in) :: &
         ph, &
-        me, &
+        en, &
         ip, &                                                                                       !< current integration point
         el                                                                                          !< current element number
     end subroutine nonlocal_dependentState
 
-    module subroutine plastic_kinehardening_deltaState(Mp,ph,me)
+    module subroutine plastic_kinehardening_deltaState(Mp,ph,en)
       real(pReal), dimension(3,3),  intent(in) :: &
         Mp                                                                                          !< Mandel stress
       integer,                      intent(in) :: &
         ph, &
-        me
+        en
     end subroutine plastic_kinehardening_deltaState
 
-    module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
+    module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
       real(pReal), dimension(3,3), intent(in) :: &
         Mp
       integer, intent(in) :: &
         ph, &
-        me
+        en
     end subroutine plastic_nonlocal_deltaState
 
   end interface
@@ -357,22 +357,22 @@ module subroutine plastic_dependentState(co, ip, el)
 
   integer :: &
     ph, &
-    me
+    en
 
 
   ph = material_phaseAt(co,el)
-  me = material_phasememberAt(co,ip,el)
+  en = material_phasememberAt(co,ip,el)
 
   plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
 
     case (PLASTICITY_DISLOTWIN_ID) plasticType
-      call dislotwin_dependentState(thermal_T(ph,me),ph,me)
+      call dislotwin_dependentState(thermal_T(ph,en),ph,en)
 
     case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
-      call dislotungsten_dependentState(ph,me)
+      call dislotungsten_dependentState(ph,en)
 
     case (PLASTICITY_NONLOCAL_ID) plasticType
-      call nonlocal_dependentState(ph,me,ip,el)
+      call nonlocal_dependentState(ph,en,ip,el)
 
   end select plasticType
 

src/phase_mechanical_plastic_dislotungsten.f90

@@ -17,7 +17,7 @@ submodule(phase:plastic) dislotungsten
       D_0  = 1.0_pReal, &                                                                           !< prefactor for self-diffusion coefficient
       Q_cl = 1.0_pReal                                                                              !< activation energy for dislocation climb
     real(pReal),               allocatable, dimension(:) :: &
-      b_sl, &                                                                                       !< magnitude me Burgers vector [m]
+      b_sl, &                                                                                       !< magnitude of Burgers vector [m]
       D_a, &
       i_sl, &                                                                                       !< Adj. parameter for distance between 2 forest dislocations
       f_at, &                                                                                       !< factor to calculate atomic volume
@@ -271,7 +271,7 @@ end function plastic_dislotungsten_init
 !> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
 pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
-                                                         Mp,T,ph,me)
+                                                         Mp,T,ph,en)
   real(pReal), dimension(3,3),     intent(out) :: &
     Lp                                                                                              !< plastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out) :: &
@@ -283,7 +283,7 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
     T                                                                                               !< temperature
   integer,                     intent(in) :: &
     ph, &
-    me
+    en
 
   integer :: &
     i,k,l,m,n
@@ -296,7 +296,7 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
 
   associate(prm => param(ph))
 
-  call kinetics(Mp,T,ph,me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
+  call kinetics(Mp,T,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
   do i = 1, prm%sum_N_sl
     Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -313,7 +313,7 @@ end subroutine dislotungsten_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine dislotungsten_dotState(Mp,T,ph,me)
+module subroutine dislotungsten_dotState(Mp,T,ph,en)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -321,7 +321,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
     T                                                                                               !< temperature
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal) :: &
     VacancyDiffusion
@@ -337,11 +337,11 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
   associate(prm => param(ph), stt => state(ph),&
             dot => dotState(ph), dst => dependentState(ph))
 
-  call kinetics(Mp,T,ph,me,&
+  call kinetics(Mp,T,ph,en,&
                 gdot_pos,gdot_neg, &
                 tau_pos_out = tau_pos,tau_neg_out = tau_neg)
 
-  dot%gamma_sl(:,me) = (gdot_pos+gdot_neg)                                                          ! ToDo: needs to be abs
+  dot%gamma_sl(:,en) = (gdot_pos+gdot_neg)                                                          ! ToDo: needs to be abs
   VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
 
   where(dEq0(tau_pos))                                                                              ! ToDo: use avg of pos and neg
@@ -350,20 +350,20 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
   else where
     dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), &
                              prm%D_a, &                                                             ! lower limit
-                             dst%Lambda_sl(:,me))                                                   ! upper limit
+                             dst%Lambda_sl(:,en))                                                   ! upper limit
-    dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,me)*abs(dot%gamma_sl(:,me))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation
+    dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation
                                   0.0_pReal, &
                                   prm%dipoleformation)
     v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) &
                   * (1.0_pReal/(dip_distance+prm%D_a))
-    dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,me))/(dip_distance-prm%D_a)                   ! ToDo: Discuss with Franz: Stress dependency?
+    dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(dip_distance-prm%D_a)                   ! ToDo: Discuss with Franz: Stress dependency?
   end where
 
-  dot%rho_mob(:,me) = abs(dot%gamma_sl(:,me))/(prm%b_sl*dst%Lambda_sl(:,me)) &                      ! multiplication
+  dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) &                      ! multiplication
                     - dot_rho_dip_formation &
-                    - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,me)*abs(dot%gamma_sl(:,me))        ! Spontaneous annihilation of 2 single edge dislocations
+                    - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))        ! Spontaneous annihilation of 2 single edge dislocations
-  dot%rho_dip(:,me) = dot_rho_dip_formation &
+  dot%rho_dip(:,en) = dot_rho_dip_formation &
-                    - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,me)*abs(dot%gamma_sl(:,me)) &      ! Spontaneous annihilation of a single edge dislocation with a dipole constituent
+                    - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) &      ! Spontaneous annihilation of a single edge dislocation with a dipole constituent
                     - dot_rho_dip_climb
 
   end associate
@@ -374,22 +374,22 @@ end subroutine dislotungsten_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate derived quantities from state.
 !--------------------------------------------------------------------------------------------------
-module subroutine dislotungsten_dependentState(ph,me)
+module subroutine dislotungsten_dependentState(ph,en)
 
   integer,      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dislocationSpacing
 
   associate(prm => param(ph), stt => state(ph),dst => dependentState(ph))
 
-  dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))
+  dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
-  dst%threshold_stress(:,me) = prm%mu*prm%b_sl &
+  dst%threshold_stress(:,en) = prm%mu*prm%b_sl &
-                             * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,me)+stt%rho_dip(:,me)))
+                             * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
-  dst%Lambda_sl(:,me) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
+  dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
 
   end associate
 
@@ -438,7 +438,7 @@ end subroutine plastic_dislotungsten_results
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics(Mp,T,ph,me, &
+pure subroutine kinetics(Mp,T,ph,en, &
                  dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out)
 
   real(pReal), dimension(3,3),  intent(in) :: &
@@ -447,7 +447,7 @@ pure subroutine kinetics(Mp,T,ph,me, &
     T                                                                                               !< temperature
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal),                  intent(out), dimension(param(ph)%sum_N_sl) :: &
     dot_gamma_pos, &
@@ -478,11 +478,11 @@ pure subroutine kinetics(Mp,T,ph,me, &
   if (present(tau_neg_out)) tau_neg_out = tau_neg
 
   associate(BoltzmannRatio  => prm%Q_s/(kB*T), &
-            dot_gamma_0     => stt%rho_mob(:,me)*prm%b_sl*prm%v_0, &
+            dot_gamma_0     => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, &
-            effectiveLength => dst%Lambda_sl(:,me) - prm%w)
+            effectiveLength => dst%Lambda_sl(:,en) - prm%w)
 
-  significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,me) > tol_math_check)
+  significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check)
-    StressRatio = (abs(tau_pos)-dst%threshold_stress(:,me))/prm%tau_Peierls
+    StressRatio = (abs(tau_pos)-dst%threshold_stress(:,en))/prm%tau_Peierls
     StressRatio_p       = StressRatio** prm%p
     StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
     needsGoodName       = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
@@ -498,7 +498,7 @@ pure subroutine kinetics(Mp,T,ph,me, &
   end where significantPositiveTau
 
   if (present(ddot_gamma_dtau_pos)) then
-    significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,me) > tol_math_check)
+    significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check)
       dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
           * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
       dtk = -1.0_pReal * t_k / tau_pos
@@ -511,8 +511,8 @@ pure subroutine kinetics(Mp,T,ph,me, &
     end where significantPositiveTau2
   endif
 
-  significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,me) > tol_math_check)
+  significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check)
-    StressRatio = (abs(tau_neg)-dst%threshold_stress(:,me))/prm%tau_Peierls
+    StressRatio = (abs(tau_neg)-dst%threshold_stress(:,en))/prm%tau_Peierls
     StressRatio_p       = StressRatio** prm%p
     StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
     needsGoodName       = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
@@ -528,7 +528,7 @@ pure subroutine kinetics(Mp,T,ph,me, &
   end where significantNegativeTau
 
   if (present(ddot_gamma_dtau_neg)) then
-    significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,me) > tol_math_check)
+    significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check)
       dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
           * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
       dtk = -1.0_pReal * t_k / tau_neg

src/phase_mechanical_plastic_dislotwin.f90

@@ -519,12 +519,12 @@ end function plastic_dislotwin_homogenizedC
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
-module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
 
   real(pReal), dimension(3,3),     intent(out) :: Lp
   real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp
   real(pReal), dimension(3,3),     intent(in)  :: Mp
-  integer,                         intent(in)  :: ph,me
+  integer,                         intent(in)  :: ph,en
   real(pReal),                     intent(in)  :: T
 
   integer :: i,k,l,m,n
@@ -565,13 +565,13 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
   associate(prm => param(ph), stt => state(ph))
 
   f_unrotated = 1.0_pReal &
-              - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
+              - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-              - sum(stt%f_tr(1:prm%sum_N_tr,me))
+              - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
 
-  call kinetics_slip(Mp,T,ph,me,dot_gamma_sl,ddot_gamma_dtau_slip)
+  call kinetics_slip(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_slip)
   slipContribution: do i = 1, prm%sum_N_sl
     Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -579,7 +579,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
                        + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
   enddo slipContribution
 
-  call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw,ddot_gamma_dtau_tw)
+  call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
   twinContibution: do i = 1, prm%sum_N_tw
     Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -587,7 +587,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
                        + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
   enddo twinContibution
 
-  call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_tr)
+  call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
   transContibution: do i = 1, prm%sum_N_tr
     Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -632,7 +632,7 @@ end subroutine dislotwin_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine dislotwin_dotState(Mp,T,ph,me)
+module subroutine dislotwin_dotState(Mp,T,ph,en)
 
   real(pReal), dimension(3,3),  intent(in):: &
     Mp                                                                                              !< Mandel stress
@@ -640,7 +640,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
     T                                                                                               !< temperature at integration point
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   integer :: i
   real(pReal) :: &
@@ -664,11 +664,11 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
             dot => dotState(ph), dst => dependentState(ph))
 
   f_unrotated = 1.0_pReal &
-              - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
+              - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-              - sum(stt%f_tr(1:prm%sum_N_tr,me))
+              - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-  call kinetics_slip(Mp,T,ph,me,dot_gamma_sl)
+  call kinetics_slip(Mp,T,ph,en,dot_gamma_sl)
-  dot%gamma_sl(:,me) = abs(dot_gamma_sl)
+  dot%gamma_sl(:,en) = abs(dot_gamma_sl)
 
   rho_dip_distance_min = prm%D_a*prm%b_sl
 
@@ -680,11 +680,11 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
       dot_rho_dip_climb(i) = 0.0_pReal
     else significantSlipStress
       rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
-      rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me))
+      rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,en))
       rho_dip_distance = math_clip(rho_dip_distance, left  = rho_dip_distance_min(i))
 
       dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) &
-                               * stt%rho_mob(i,me)*abs(dot_gamma_sl(i))
+                               * stt%rho_mob(i,en)*abs(dot_gamma_sl(i))
 
       if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then
         dot_rho_dip_climb(i) = 0.0_pReal
@@ -698,25 +698,25 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
         v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
              * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
 
-        dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,me) &
+        dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
                              / (rho_dip_distance-rho_dip_distance_min(i))
       endif
     endif significantSlipStress
   enddo slipState
 
-  dot%rho_mob(:,me) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,me)) &
+  dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) &
                     - dot_rho_dip_formation &
-                    - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,me)*abs(dot_gamma_sl)
+                    - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl)
 
-  dot%rho_dip(:,me) = dot_rho_dip_formation &
+  dot%rho_dip(:,en) = dot_rho_dip_formation &
-                    - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) &
+                    - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) &
                     - dot_rho_dip_climb
 
-  call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw)
+  call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
-  dot%f_tw(:,me) = f_unrotated*dot_gamma_tw/prm%gamma_char
+  dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
 
-  call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr)
+  call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
-  dot%f_tr(:,me) = f_unrotated*dot_gamma_tr
+  dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
 
   end associate
 
@@ -726,11 +726,11 @@ end subroutine dislotwin_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate derived quantities from state.
 !--------------------------------------------------------------------------------------------------
-module subroutine dislotwin_dependentState(T,ph,me)
+module subroutine dislotwin_dependentState(T,ph,en)
 
   integer,       intent(in) :: &
     ph, &
-    me
+    en
   real(pReal),   intent(in) :: &
     T
 
@@ -752,50 +752,50 @@ module subroutine dislotwin_dependentState(T,ph,me)
             stt => state(ph),&
             dst => dependentState(ph))
 
-  sumf_tw  = sum(stt%f_tw(1:prm%sum_N_tw,me))
+  sumf_tw  = sum(stt%f_tw(1:prm%sum_N_tw,en))
-  sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,me))
+  sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
 
   Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
 
   !* rescaled volume fraction for topology
-  f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw                                                ! this is per system ...
+  f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw                                                ! this is per system ...
   f_over_t_tr = sumf_tr/prm%t_tr                                                                    ! but this not
                                                                                                     ! ToDo ...Physically correct, but naming could be adjusted
 
-  inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl
+  inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl
   if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
     inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
   if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
     inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
-  dst%Lambda_sl(:,me) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
+  dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
 
   inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
-  dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
+  dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
 
   inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
-  dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
+  dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
 
   !* threshold stress for dislocation motion
-  dst%tau_pass(:,me) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,me)+stt%rho_dip(:,me)))
+  dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
   !* threshold stress for growing twin/martensite
   if(prm%sum_N_tw == prm%sum_N_sl) &
-    dst%tau_hat_tw(:,me) = Gamma/(3.0_pReal*prm%b_tw) &
+    dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
                          + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
   if(prm%sum_N_tr == prm%sum_N_sl) &
-    dst%tau_hat_tr(:,me) = Gamma/(3.0_pReal*prm%b_tr) &
+    dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
                          + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
                          + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
 
-  dst%V_tw(:,me) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,me)**2.0_pReal
+  dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
-  dst%V_tr(:,me) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,me)**2.0_pReal
+  dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
 
 
   x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)        ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
-  dst%tau_r_tw(:,me) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
+  dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
 
   x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)        ! ToDo: In the paper, this is the Burgers vector for slip
-  dst%tau_r_tr(:,me) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
+  dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
 
   end associate
 
@@ -860,7 +860,7 @@ end subroutine plastic_dislotwin_results
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_slip(Mp,T,ph,me, &
+pure subroutine kinetics_slip(Mp,T,ph,en, &
                               dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip)
 
   real(pReal), dimension(3,3),  intent(in) :: &
@@ -869,7 +869,7 @@ pure subroutine kinetics_slip(Mp,T,ph,me, &
     T                                                                                               !< temperature
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: &
     dot_gamma_sl
@@ -898,7 +898,7 @@ pure subroutine kinetics_slip(Mp,T,ph,me, &
     tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
   enddo
 
-  tau_eff = abs(tau)-dst%tau_pass(:,me)
+  tau_eff = abs(tau)-dst%tau_pass(:,en)
 
   significantStress: where(tau_eff > tol_math_check)
     stressRatio    = tau_eff/prm%tau_0
@@ -907,7 +907,7 @@ pure subroutine kinetics_slip(Mp,T,ph,me, &
     v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q)
     v_run_inverse  = prm%B/(tau_eff*prm%b_sl)
 
-    dot_gamma_sl = sign(stt%rho_mob(:,me)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)
+    dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)
 
     dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio &
                          * (stressRatio**(prm%p-1.0_pReal)) &
@@ -916,7 +916,7 @@ pure subroutine kinetics_slip(Mp,T,ph,me, &
     dV_run_inverse_dTau  = -1.0_pReal * v_run_inverse/tau_eff
     dV_dTau              = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) &
                          / (v_wait_inverse+v_run_inverse)**2.0_pReal
-    ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,me)*prm%b_sl
+    ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl
   else where significantStress
     dot_gamma_sl    = 0.0_pReal
     ddot_gamma_dtau = 0.0_pReal
@@ -937,7 +937,7 @@ end subroutine kinetics_slip
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
+pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,&
                               dot_gamma_tw,ddot_gamma_dtau_tw)
 
   real(pReal), dimension(3,3),  intent(in) :: &
@@ -946,7 +946,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
     T                                                                                               !< temperature
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
   real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
     dot_gamma_sl
 
@@ -970,11 +970,11 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
     isFCC: if (prm%fccTwinTransNucleation) then
       s1=prm%fcc_twinNucleationSlipPair(1,i)
       s2=prm%fcc_twinNucleationSlipPair(2,i)
-      if (tau(i) < dst%tau_r_tw(i,me)) then                                                         ! ToDo: correct?
+      if (tau(i) < dst%tau_r_tw(i,en)) then                                                         ! ToDo: correct?
-        Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,me)+stt%rho_dip(s2,me))+&
+        Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
-               abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,me)+stt%rho_dip(s1,me)))/&                     ! ToDo: MD: it would be more consistent to use shearrates from state
+               abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&                     ! ToDo: MD: it would be more consistent to use shearrates from state
                 (prm%L_tw*prm%b_sl(i))*&
-                (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,me)-tau(i))))                       ! P_ncs
+                (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i))))                       ! P_ncs
       else
         Ndot0=0.0_pReal
       end if
@@ -984,8 +984,8 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
   enddo
 
   significantStress: where(tau > tol_math_check)
-    StressRatio_r   = (dst%tau_hat_tw(:,me)/tau)**prm%r
+    StressRatio_r   = (dst%tau_hat_tw(:,en)/tau)**prm%r
-    dot_gamma_tw    = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r)
+    dot_gamma_tw    = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r)
     ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
   else where significantStress
     dot_gamma_tw    = 0.0_pReal
@@ -1006,7 +1006,7 @@ end subroutine kinetics_twin
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
+pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,&
                               dot_gamma_tr,ddot_gamma_dtau_tr)
 
   real(pReal), dimension(3,3),  intent(in) :: &
@@ -1015,7 +1015,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
     T                                                                                               !< temperature
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
   real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
     dot_gamma_sl
 
@@ -1038,11 +1038,11 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
     isFCC: if (prm%fccTwinTransNucleation) then
       s1=prm%fcc_twinNucleationSlipPair(1,i)
       s2=prm%fcc_twinNucleationSlipPair(2,i)
-      if (tau(i) < dst%tau_r_tr(i,me)) then                                                         ! ToDo: correct?
+      if (tau(i) < dst%tau_r_tr(i,en)) then                                                         ! ToDo: correct?
-        Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,me)+stt%rho_dip(s2,me))+&
+        Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
-               abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,me)+stt%rho_dip(s1,me)))/&                     ! ToDo: MD: it would be more consistent to use shearrates from state
+               abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&                     ! ToDo: MD: it would be more consistent to use shearrates from state
                 (prm%L_tr*prm%b_sl(i))*&
-                (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,me)-tau(i))))                       ! P_ncs
+                (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i))))                       ! P_ncs
       else
         Ndot0=0.0_pReal
       end if
@@ -1052,8 +1052,8 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
   enddo
 
   significantStress: where(tau > tol_math_check)
-    StressRatio_s   = (dst%tau_hat_tr(:,me)/tau)**prm%s
+    StressRatio_s   = (dst%tau_hat_tr(:,en)/tau)**prm%s
-    dot_gamma_tr    = dst%V_tr(:,me) * Ndot0*exp(-StressRatio_s)
+    dot_gamma_tr    = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s)
     ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
   else where significantStress
     dot_gamma_tr  = 0.0_pReal

src/phase_mechanical_plastic_isotropic.f90

@@ -156,7 +156,7 @@ end function plastic_isotropic_init
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
-module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
   real(pReal), dimension(3,3),     intent(out) :: &
     Lp                                                                                              !< plastic velocity gradient
@@ -167,7 +167,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
     Mp                                                                                              !< Mandel stress
   integer,                     intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal), dimension(3,3) :: &
     Mp_dev                                                                                          !< deviatoric part of the Mandel stress
@@ -185,7 +185,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
   norm_Mp_dev = sqrt(squarenorm_Mp_dev)
 
   if (norm_Mp_dev > 0.0_pReal) then
-    dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(me))) **prm%n
+    dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
 
     Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -248,13 +248,13 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine isotropic_dotState(Mp,ph,me)
+module subroutine isotropic_dotState(Mp,ph,en)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal) :: &
     dot_gamma, &                                                                                    !< strainrate
@@ -270,7 +270,7 @@ module subroutine isotropic_dotState(Mp,ph,me)
     norm_Mp = sqrt(math_tensordot(math_deviatoric33(Mp),math_deviatoric33(Mp)))
   endif
 
-  dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(me))) **prm%n
+  dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(en))) **prm%n
 
   if (dot_gamma > 1e-12_pReal) then
     if (dEq0(prm%c_1)) then
@@ -280,15 +280,15 @@ module subroutine isotropic_dotState(Mp,ph,me)
                   + 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
-    dot%xi(me) = dot_gamma &
+    dot%xi(en) = dot_gamma &
                * ( prm%h_0 + prm%h_ln * log(dot_gamma) ) &
-               * abs( 1.0_pReal - stt%xi(me)/xi_inf_star )**prm%a &
+               * abs( 1.0_pReal - stt%xi(en)/xi_inf_star )**prm%a &
-               * sign(1.0_pReal, 1.0_pReal - stt%xi(me)/xi_inf_star)
+               * sign(1.0_pReal, 1.0_pReal - stt%xi(en)/xi_inf_star)
   else
-    dot%xi(me) = 0.0_pReal
+    dot%xi(en) = 0.0_pReal
   endif
 
-  dot%gamma(me) = dot_gamma                                                                         ! ToDo: not really used
+  dot%gamma(en) = dot_gamma                                                                         ! ToDo: not really used
 
   end associate
 

src/phase_mechanical_plastic_kinehardening.f90

@@ -230,7 +230,7 @@ end function plastic_kinehardening_init
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
-pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
   real(pReal), dimension(3,3),     intent(out) :: &
     Lp                                                                                              !< plastic velocity gradient
@@ -241,7 +241,7 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
     Mp                                                                                              !< Mandel stress
   integer,               intent(in) :: &
     ph, &
-    me
+    en
 
   integer :: &
     i,k,l,m,n
@@ -254,7 +254,7 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
 
   associate(prm => param(ph))
 
-  call kinetics(Mp,ph,me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
+  call kinetics(Mp,ph,en,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
 
   do i = 1, prm%sum_N_sl
     Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i)
@@ -272,13 +272,13 @@ end subroutine kinehardening_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_kinehardening_dotState(Mp,ph,me)
+module subroutine plastic_kinehardening_dotState(Mp,ph,en)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal) :: &
     sumGamma
@@ -289,22 +289,22 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,me)
   associate(prm => param(ph), stt => state(ph),&
             dot => dotState(ph))
 
-  call kinetics(Mp,ph,me,gdot_pos,gdot_neg)
+  call kinetics(Mp,ph,en,gdot_pos,gdot_neg)
-  dot%accshear(:,me) = abs(gdot_pos+gdot_neg)
+  dot%accshear(:,en) = abs(gdot_pos+gdot_neg)
-  sumGamma = sum(stt%accshear(:,me))
+  sumGamma = sum(stt%accshear(:,en))
 
 
-  dot%crss(:,me) = matmul(prm%interaction_SlipSlip,dot%accshear(:,me)) &
+  dot%crss(:,en) = matmul(prm%interaction_SlipSlip,dot%accshear(:,en)) &
                  * (  prm%h_inf_f &
                      + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
                      * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
                    )
 
-  dot%crss_back(:,me) = stt%sense(:,me)*dot%accshear(:,me) * &
+  dot%crss_back(:,en) = stt%sense(:,en)*dot%accshear(:,en) * &
            ( prm%h_inf_b + &
              (prm%h_0_b - prm%h_inf_b &
-               + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,me))*(stt%accshear(:,me)-stt%gamma0(:,me))&
+               + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,en))*(stt%accshear(:,en)-stt%gamma0(:,en))&
-             ) *exp(-(stt%accshear(:,me)-stt%gamma0(:,me)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,me))) &
+             ) *exp(-(stt%accshear(:,en)-stt%gamma0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,en))) &
            )
 
   end associate
@@ -315,13 +315,13 @@ end subroutine plastic_kinehardening_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate (instantaneous) incremental change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_kinehardening_deltaState(Mp,ph,me)
+module subroutine plastic_kinehardening_deltaState(Mp,ph,en)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     gdot_pos,gdot_neg, &
@@ -329,22 +329,22 @@ module subroutine plastic_kinehardening_deltaState(Mp,ph,me)
 
   associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph))
 
-  call kinetics(Mp,ph,me,gdot_pos,gdot_neg)
+  call kinetics(Mp,ph,en,gdot_pos,gdot_neg)
-  sense = merge(state(ph)%sense(:,me), &                                                            ! keep existing...
+  sense = merge(state(ph)%sense(:,en), &                                                            ! keep existing...
                 sign(1.0_pReal,gdot_pos+gdot_neg), &                                                ! ...or have a defined
                 dEq0(gdot_pos+gdot_neg,1e-10_pReal))                                                ! current sense of shear direction
 
 
 !--------------------------------------------------------------------------------------------------
-! switch in sense me shear?
+! switch in sense of shear?
-  where(dNeq(sense,stt%sense(:,me),0.1_pReal))
+  where(dNeq(sense,stt%sense(:,en),0.1_pReal))
-    dlt%sense (:,me) = sense - stt%sense(:,me)                                                      ! switch sense
+    dlt%sense (:,en) = sense - stt%sense(:,en)                                                      ! switch sense
-    dlt%chi0  (:,me) = abs(stt%crss_back(:,me)) - stt%chi0(:,me)                                    ! remember current backstress magnitude
+    dlt%chi0  (:,en) = abs(stt%crss_back(:,en)) - stt%chi0(:,en)                                    ! remember current backstress magnitude
-    dlt%gamma0(:,me) = stt%accshear(:,me) - stt%gamma0(:,me)                                        ! remember current accumulated shear
+    dlt%gamma0(:,en) = stt%accshear(:,en) - stt%gamma0(:,en)                                        ! remember current accumulated shear
   else where
-    dlt%sense (:,me) = 0.0_pReal
+    dlt%sense (:,en) = 0.0_pReal
-    dlt%chi0  (:,me) = 0.0_pReal
+    dlt%chi0  (:,en) = 0.0_pReal
-    dlt%gamma0(:,me) = 0.0_pReal
+    dlt%gamma0(:,en) = 0.0_pReal
   end where
 
   end associate
@@ -397,14 +397,14 @@ end subroutine plastic_kinehardening_results
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics(Mp,ph,me, &
+pure subroutine kinetics(Mp,ph,en, &
                          gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal),                  intent(out), dimension(param(ph)%sum_N_sl) :: &
     gdot_pos, &
@@ -421,21 +421,21 @@ pure subroutine kinetics(Mp,ph,me, &
   associate(prm => param(ph), stt => state(ph))
 
   do i = 1, prm%sum_N_sl
-    tau_pos(i) =       math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,me)
+    tau_pos(i) =       math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,en)
-    tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,me), &
+    tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,en), &
                        0.0_pReal, prm%nonSchmidActive)
   enddo
 
   where(dNeq0(tau_pos))
     gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &                  ! 1/2 if non-Schmid active
-             * sign(abs(tau_pos/stt%crss(:,me))**prm%n,  tau_pos)
+             * sign(abs(tau_pos/stt%crss(:,en))**prm%n,  tau_pos)
   else where
     gdot_pos = 0.0_pReal
   end where
 
   where(dNeq0(tau_neg))
     gdot_neg = prm%dot_gamma_0 * 0.5_pReal &                                                        ! only used if non-Schmid active, always 1/2
-             * sign(abs(tau_neg/stt%crss(:,me))**prm%n,  tau_neg)
+             * sign(abs(tau_neg/stt%crss(:,en))**prm%n,  tau_neg)
   else where
     gdot_neg = 0.0_pReal
   end where

src/phase_mechanical_plastic_nonlocal.f90

@@ -49,7 +49,7 @@ submodule(phase:plastic) nonlocal
   !END DEPRECATED
 
   real(pReal), dimension(:,:,:,:,:,:), allocatable :: &
-    compatibility                                                                                   !< slip system compatibility between me and my neighbors
+    compatibility                                                                                   !< slip system compatibility between en and my neighbors
 
   type :: tInitialParameters                                                                        !< container type for internal constitutive parameters
     real(pReal) :: &
@@ -558,11 +558,11 @@ end function plastic_nonlocal_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates quantities characterizing the microstructure
 !--------------------------------------------------------------------------------------------------
-module subroutine nonlocal_dependentState(ph, me, ip, el)
+module subroutine nonlocal_dependentState(ph, en, ip, el)
 
   integer, intent(in) :: &
     ph, &
-    me, &
+    en, &
     ip, &
     el
 
@@ -615,9 +615,9 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
 
   associate(prm => param(ph),dst => microstructure(ph), stt => state(ph))
 
-  rho = getRho(ph,me)
+  rho = getRho(ph,en)
 
-  stt%rho_forest(:,me) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
+  stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
                        + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
 
 
@@ -627,13 +627,13 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
     myInteractionMatrix = prm%h_sl_sl &
                         * spread((  1.0_pReal - prm%f_F &
                                    + prm%f_F &
-                                   * log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,me),prm%rho_significant))) &
+                                   * log(0.35_pReal * prm%b_sl * sqrt(max(stt%rho_forest(:,en),prm%rho_significant))) &
                                    / log(0.35_pReal * prm%b_sl * 1e6_pReal))** 2.0_pReal,2,prm%sum_N_sl)
   else
     myInteractionMatrix = prm%h_sl_sl
   endif
 
-  dst%tau_pass(:,me) = prm%mu * prm%b_sl &
+  dst%tau_pass(:,en) = prm%mu * prm%b_sl &
                      * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
 
 !*** calculate the dislocation stress of the neighboring excess dislocation densities
@@ -643,10 +643,10 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
   ! ToDo: MD: this is most likely only correct for F_i = I
   !#################################################################################################
 
-  rho0 = getRho0(ph,me)
+  rho0 = getRho0(ph,en)
   if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
-    invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))
+    invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
-    invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me))
+    invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en))
 
     rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
     rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
@@ -654,7 +654,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
     rhoExcess(1,:) = rho_edg_delta
     rhoExcess(2,:) = rho_scr_delta
 
-    FVsize = geom(ph)%V_0(me) ** (1.0_pReal/3.0_pReal)
+    FVsize = geom(ph)%V_0(en) ** (1.0_pReal/3.0_pReal)
 
     !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
 
@@ -734,7 +734,7 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
       where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal
 
         ! ... gives the local stress correction when multiplied with a factor
-      dst%tau_back(s,me) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
+      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
@@ -745,9 +745,9 @@ module subroutine nonlocal_dependentState(ph, me, ip, el)
       .and. ((debugConstitutive%element == el .and. debugConstitutive%ip == ip)&
              .or. .not. debugConstitutive%selective)) then
     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(:,me)
+    print'(a,/,12x,12(e10.3,1x))', '<< CONST >> rhoForest', stt%rho_forest(:,en)
-    print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,me)*1e-6
+    print'(a,/,12x,12(f10.5,1x))', '<< CONST >> tauThreshold / MPa', dst%tau_pass(:,en)*1e-6
-    print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,me)*1e-6
+    print'(a,/,12x,12(f10.5,1x),/)', '<< CONST >> tauBack / MPa', dst%tau_back(:,en)*1e-6
   endif
 #endif
 
@@ -760,14 +760,14 @@ end subroutine nonlocal_dependentState
 !> @brief calculates plastic velocity gradient and its tangent
 !--------------------------------------------------------------------------------------------------
 module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
-                                                   Mp,Temperature,ph,me)
+                                                   Mp,Temperature,ph,en)
   real(pReal), dimension(3,3), intent(out) :: &
     Lp                                                                                              !< plastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out) :: &
     dLp_dMp
   integer, intent(in) :: &
     ph, &
-    me
+    en
   real(pReal), intent(in) :: &
     Temperature                                                                                     !< temperature
 
@@ -800,7 +800,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
   ns = prm%sum_N_sl
 
   !*** shortcut to state variables
-  rho = getRho(ph,me)
+  rho = getRho(ph,en)
   rhoSgl = rho(:,sgl)
 
   do s = 1,ns
@@ -815,12 +815,12 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
       tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_pos(1:3,1:3,s))
     endif
   enddo
-  tauNS = tauNS + spread(dst%tau_back(:,me),2,4)
+  tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
-  tau   = tau   + dst%tau_back(:,me)
+  tau   = tau   + dst%tau_back(:,en)
 
   ! edges
   call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
-                tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, ph)
+                tau, tauNS(:,1), dst%tau_pass(:,en),1,Temperature, ph)
   v(:,2)         = v(:,1)
   dv_dtau(:,2)   = dv_dtau(:,1)
   dv_dtauNS(:,2) = dv_dtauNS(:,1)
@@ -829,7 +829,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
   if (prm%nonSchmidActive) then
     do t = 3,4
       call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
-                    tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, ph)
+                    tau, tauNS(:,t), dst%tau_pass(:,en),2,Temperature, ph)
     enddo
   else
     v(:,3:4)         = spread(v(:,1),2,2)
@@ -837,7 +837,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
     dv_dtauNS(:,3:4) = spread(dv_dtauNS(:,1),2,2)
   endif
 
-  stt%v(:,me) = pack(v,.true.)
+  stt%v(:,en) = pack(v,.true.)
 
   !*** Bauschinger effect
   forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
@@ -866,13 +866,13 @@ end subroutine nonlocal_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief (instantaneous) incremental change of microstructure
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
+module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
 
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< MandelStress
   integer, intent(in) :: &
     ph, &
-    me
+    en
 
   integer :: &
     ns, &                                                                                           ! short notation for the total number of active slip systems
@@ -898,10 +898,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
   ns = prm%sum_N_sl
 
   !*** shortcut to state variables
-  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me)
+  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
-  forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),me)
+  forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
 
-  rho =  getRho(ph,me)
+  rho =  getRho(ph,en)
   rhoDip = rho(:,dip)
 
   !****************************************************************************
@@ -922,7 +922,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
 
   !*** calculate limits for stable dipole height
   do s = 1,prm%sum_N_sl
-    tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,me)
+    tau(s) = math_tensordot(Mp, prm%Schmid(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
 
@@ -946,10 +946,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,me)
                                        / (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
 
   forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
-  forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),me) = dUpper(s,c)
+  forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c)
 
-  plasticState(ph)%deltaState(:,me) = 0.0_pReal
+  plasticState(ph)%deltaState(:,en) = 0.0_pReal
-  del%rho(:,me) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
+  del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
 
   end associate
 
@@ -960,7 +960,7 @@ end subroutine plastic_nonlocal_deltaState
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
 module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
-                                            ph,me,ip,el)
+                                            ph,en,ip,el)
 
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< MandelStress
@@ -969,7 +969,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
     timestep                                                                                        !< substepped crystallite time increment
   integer, intent(in) :: &
     ph, &
-    me, &
+    en, &
     ip, &                                                                                           !< current integration point
     el                                                                                              !< current element number
 
@@ -1017,13 +1017,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   tau = 0.0_pReal
   gdot = 0.0_pReal
 
-  rho  = getRho(ph,me)
+  rho  = getRho(ph,en)
   rhoSgl = rho(:,sgl)
   rhoDip = rho(:,dip)
-  rho0 = getRho0(ph,me)
+  rho0 = getRho0(ph,en)
   my_rhoSgl0 = rho0(:,sgl)
 
-  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me)
+  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
   gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
 
 #ifdef DEBUG
@@ -1038,7 +1038,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   !****************************************************************************
   !*** limits for stable dipole height
   do s = 1,ns
-    tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,me)
+    tau(s) = math_tensordot(Mp, prm%Schmid(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
 
@@ -1059,20 +1059,20 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then
     forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal)
       rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) &                           ! assuming double-cross-slip of screws to be decisive for multiplication
-                                  * sqrt(stt%rho_forest(s,me)) / prm%i_sl(s) ! &                    ! mean free path
+                                  * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! &                    ! mean free path
                                   ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4)))           ! ratio of screw to overall velocity determines edge generation
       rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) &                            ! assuming double-cross-slip of screws to be decisive for multiplication
-                                  * sqrt(stt%rho_forest(s,me)) / prm%i_sl(s) ! &                    ! mean free path
+                                  * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! &                    ! mean free path
                                   ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4)))           ! ratio of edge to overall velocity determines screw generation
     endforall
 
   else isBCC
     rhoDotMultiplication(:,1:4) = spread( &
           (sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) &
-        * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4)
+        * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4)
   endif isBCC
 
-  forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me)
+  forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
 
 
   !****************************************************************************
@@ -1115,10 +1115,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   if (lattice_structure(ph) == LATTICE_fcc_ID) &
     forall (s = 1:ns, prm%colinearSystem(s) > 0) &
       rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) &
-        * 0.25_pReal * sqrt(stt%rho_forest(s,me)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
+        * 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
 
 
-  !*** thermally activated annihilation me edge dipoles by climb
+  !*** thermally activated annihilation of edge dipoles by climb
   rhoDotThermalAnnihilation = 0.0_pReal
   selfDiffusion = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature))
   vClimb = prm%V_at * selfDiffusion * prm%mu &
@@ -1128,7 +1128,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
                                          - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
                                                                   - rhoDotSingle2DipoleGlide(s,9))  ! make sure that we do not annihilate more dipoles than we have
 
-  rhoDot = rhoDotFlux(timestep, ph,me,ip,el) &
+  rhoDot = rhoDotFlux(timestep, ph,en,ip,el) &
          + rhoDotMultiplication &
          + rhoDotSingle2DipoleGlide &
          + rhoDotAthermalAnnihilation &
@@ -1145,8 +1145,8 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
 #endif
     plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
   else
-    dot%rho(:,me) = pack(rhoDot,.true.)
+    dot%rho(:,en) = pack(rhoDot,.true.)
-    dot%gamma(:,me) = sum(gdot,2)
+    dot%gamma(:,en) = sum(gdot,2)
   endif
 
   end associate
@@ -1157,13 +1157,13 @@ end subroutine nonlocal_dotState
 !---------------------------------------------------------------------------------------------------
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
-function rhoDotFlux(timestep,ph,me,ip,el)
+function rhoDotFlux(timestep,ph,en,ip,el)
 
   real(pReal), intent(in) :: &
     timestep                                                                                        !< substepped crystallite time increment
   integer, intent(in) :: &
     ph, &
-    me, &
+    en, &
     ip, &                                                                                           !< current integration point
     el                                                                                              !< current element number
 
@@ -1174,11 +1174,11 @@ function rhoDotFlux(timestep,ph,me,ip,el)
     n, &                                                                                            !< index of my current neighbor
     neighbor_el, &                                                                                  !< element number of my neighbor
     neighbor_ip, &                                                                                  !< integration point of my neighbor
-    neighbor_n, &                                                                                   !< neighbor index pointing to me when looking from my neighbor
+    neighbor_n, &                                                                                   !< neighbor index pointing to en when looking from my neighbor
     opposite_neighbor, &                                                                            !< index of my opposite neighbor
     opposite_ip, &                                                                                  !< ip of my opposite neighbor
     opposite_el, &                                                                                  !< element index of my opposite neighbor
-    opposite_n, &                                                                                   !< neighbor index pointing to me when looking from my opposite neighbor
+    opposite_n, &                                                                                   !< neighbor index pointing to en when looking from my opposite neighbor
     t, &                                                                                            !< type of dislocation
     no,&                                                                                            !< neighbor offset shortcut
     np,&                                                                                            !< neighbor phase shortcut
@@ -1204,12 +1204,12 @@ function rhoDotFlux(timestep,ph,me,ip,el)
     neighbor_F, &                                                                                   !< total deformation gradient of my neighbor
     my_Fe, &                                                                                        !< my elastic deformation gradient
     neighbor_Fe, &                                                                                  !< elastic deformation gradient of my neighbor
-    Favg                                                                                            !< average total deformation gradient of me and my neighbor
+    Favg                                                                                            !< average total deformation gradient of en and my neighbor
   real(pReal), dimension(3) :: &
-    normal_neighbor2me, &                                                                           !< interface normal pointing from my neighbor to me in neighbor's lattice configuration
+    normal_neighbor2me, &                                                                           !< interface normal pointing from my neighbor to en in neighbor's lattice configuration
-    normal_neighbor2me_defConf, &                                                                   !< interface normal pointing from my neighbor to me in shared deformed configuration
+    normal_neighbor2me_defConf, &                                                                   !< interface normal pointing from my neighbor to en in shared deformed configuration
-    normal_me2neighbor, &                                                                           !< interface normal pointing from me to my neighbor in my lattice configuration
+    normal_me2neighbor, &                                                                           !< interface normal pointing from en to my neighbor in my lattice configuration
-    normal_me2neighbor_defConf                                                                      !< interface normal pointing from me to my neighbor in shared deformed configuration
+    normal_me2neighbor_defConf                                                                      !< interface normal pointing from en to my neighbor in shared deformed configuration
   real(pReal) :: &
     area, &                                                                                         !< area of the current interface
     transmissivity, &                                                                               !< overall transmissivity of dislocation flux to neighboring material point
@@ -1224,16 +1224,16 @@ function rhoDotFlux(timestep,ph,me,ip,el)
 
   gdot = 0.0_pReal
 
-  rho  = getRho(ph,me)
+  rho  = getRho(ph,en)
   rhoSgl = rho(:,sgl)
-  rho0 = getRho0(ph,me)
+  rho0 = getRho0(ph,en)
   my_rhoSgl0 = rho0(:,sgl)
 
-  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),me)                   !ToDo: MD: I think we should use state0 here
+  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)                   !ToDo: MD: I think we should use state0 here
   gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
 
 
-  forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),me)
+  forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
 
   !****************************************************************************
   !*** calculate dislocation fluxes (only for nonlocal plasticity)
@@ -1268,8 +1268,8 @@ function rhoDotFlux(timestep,ph,me,ip,el)
     m(1:3,:,3) = -prm%slip_transverse
     m(1:3,:,4) =  prm%slip_transverse
 
-    my_F = phase_mechanical_F(ph)%data(1:3,1:3,me)
+    my_F = phase_mechanical_F(ph)%data(1:3,1:3,en)
-    my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)))
+    my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
 
     neighbors: do n = 1,nIPneighbors
 
@@ -1316,7 +1316,7 @@ function rhoDotFlux(timestep,ph,me,ip,el)
           .or. neighbor_rhoSgl0 < prm%rho_significant) &
           neighbor_rhoSgl0 = 0.0_pReal
         normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
-                                     IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el))          ! normal of the interface in (average) deformed configuration (pointing neighbor => me)
+                                     IPareaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el))          ! normal of the interface in (average) deformed configuration (pointing neighbor => en)
         normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
                            / math_det33(neighbor_Fe)                                                ! interface normal in the lattice configuration of my neighbor
         area = IParea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me)
@@ -1325,7 +1325,7 @@ function rhoDotFlux(timestep,ph,me,ip,el)
           do t = 1,4
             c = (t + 1) / 2
             topp = t + mod(t,2) - mod(t+1,2)
-            if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal &         ! flux from my neighbor to me == entering flux for me
+            if (neighbor_v0(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal &         ! flux from my neighbor to en == entering flux for en
                 .and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then                                ! ... only if no sign change in flux density
               lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) &
                          * math_inner(m(1:3,s,t), normal_neighbor2me) * area                        ! positive line length that wants to enter through this interface
@@ -1344,7 +1344,7 @@ function rhoDotFlux(timestep,ph,me,ip,el)
 
       !* FLUX FROM ME TO MY NEIGHBOR
       !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with local properties).
-      !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me.
+      !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to en.
       !* So the net flux in the direction of my neighbor is equal to zero:
       !*    leaving flux to neighbor == entering flux from opposite neighbor
       !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
@@ -1353,7 +1353,7 @@ function rhoDotFlux(timestep,ph,me,ip,el)
       if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTICITY_NONLOCAL_ID) then
 
         normal_me2neighbor_defConf = math_det33(Favg) &
-                                   * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el))  ! normal of the interface in (average) deformed configuration (pointing me => neighbor)
+                                   * matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el))  ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
         normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
                            / math_det33(my_Fe)                                                      ! interface normal in my lattice configuration
         area = IParea(n,ip,el) * norm2(normal_me2neighbor)
@@ -1361,7 +1361,7 @@ function rhoDotFlux(timestep,ph,me,ip,el)
         do s = 1,ns
           do t = 1,4
             c = (t + 1) / 2
-            if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then             ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
+            if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then             ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
               if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then                                     ! no sign change in flux density
                 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
@@ -1403,13 +1403,13 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
 
   integer :: &
     n, &                                                                                            ! neighbor index
-    me, &
+    en, &
     neighbor_e, &                                                                                   ! element index of my neighbor
     neighbor_i, &                                                                                   ! integration point index of my neighbor
     neighbor_me, &
     neighbor_phase, &
     ns, &                                                                                           ! number of active slip systems
-    s1, &                                                                                           ! slip system index (me)
+    s1, &                                                                                           ! slip system index (en)
     s2                                                                                              ! slip system index (my neighbor)
   real(pReal), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: &
     my_compatibility                                                                                ! my_compatibility for current element and ip
@@ -1424,7 +1424,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
   associate(prm => param(ph))
   ns = prm%sum_N_sl
 
-  me = material_phaseMemberAt(1,i,e)
+  en = material_phaseMemberAt(1,i,e)
   !*** start out fully compatible
   my_compatibility = 0.0_pReal
   forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
@@ -1450,7 +1450,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
     elseif (prm%chi_GB >= 0.0_pReal) then
       !* GRAIN BOUNDARY !
       !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
-      if (any(dNeq(material_orientation0(1,ph,me)%asQuaternion(), &
+      if (any(dNeq(material_orientation0(1,ph,en)%asQuaternion(), &
                    material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. &
           (.not. phase_localPlasticity(neighbor_phase))) &
         forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
@@ -1767,21 +1767,21 @@ end subroutine kinetics
 !> @brief returns copy of current dislocation densities from state
 !> @details raw values is rectified
 !--------------------------------------------------------------------------------------------------
-pure function getRho(ph,me)
+pure function getRho(ph,en)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho
 
 
   associate(prm => param(ph))
 
-    getRho = reshape(state(ph)%rho(:,me),[prm%sum_N_sl,10])
+    getRho = reshape(state(ph)%rho(:,en),[prm%sum_N_sl,10])
 
     ! ensure positive densities (not for imm, they have a sign)
     getRho(:,mob) = max(getRho(:,mob),0.0_pReal)
     getRho(:,dip) = max(getRho(:,dip),0.0_pReal)
 
-    where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
+    where(abs(getRho) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
       getRho = 0.0_pReal
 
   end associate
@@ -1793,21 +1793,21 @@ end function getRho
 !> @brief returns copy of current dislocation densities from state
 !> @details raw values is rectified
 !--------------------------------------------------------------------------------------------------
-pure function getRho0(ph,me)
+pure function getRho0(ph,en)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   real(pReal), dimension(param(ph)%sum_N_sl,10) :: getRho0
 
 
   associate(prm => param(ph))
 
-    getRho0 = reshape(state0(ph)%rho(:,me),[prm%sum_N_sl,10])
+    getRho0 = reshape(state0(ph)%rho(:,en),[prm%sum_N_sl,10])
 
     ! ensure positive densities (not for imm, they have a sign)
     getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
     getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
 
-    where (abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
+    where (abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(en)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
       getRho0 = 0.0_pReal
 
   end associate

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -283,7 +283,7 @@ end function plastic_phenopowerlaw_init
 !> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
 !  equally (Taylor assumption). Twinning happens only in untwinned volume
 !--------------------------------------------------------------------------------------------------
-pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
   real(pReal), dimension(3,3),     intent(out) :: &
     Lp                                                                                              !< plastic velocity gradient
@@ -294,7 +294,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
     Mp                                                                                              !< Mandel stress
   integer,               intent(in) :: &
     ph, &
-    me
+    en
 
   integer :: &
     i,k,l,m,n
@@ -309,7 +309,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
 
   associate(prm => param(ph))
 
-  call kinetics_slip(Mp,ph,me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
+  call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
   slipSystems: do i = 1, prm%sum_N_sl
     Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -318,7 +318,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
                        + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
   enddo slipSystems
 
-  call kinetics_twin(Mp,ph,me,gdot_twin,dgdot_dtautwin)
+  call kinetics_twin(Mp,ph,en,gdot_twin,dgdot_dtautwin)
   twinSystems: do i = 1, prm%sum_N_tw
     Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -334,13 +334,13 @@ end subroutine phenopowerlaw_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
-module subroutine phenopowerlaw_dotState(Mp,ph,me)
+module subroutine phenopowerlaw_dotState(Mp,ph,en)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal) :: &
     c_SlipSlip,c_TwinSlip,c_TwinTwin, &
@@ -353,8 +353,8 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me)
   associate(prm => param(ph), stt => state(ph), &
   dot => dotState(ph))
 
-  sumGamma = sum(stt%gamma_slip(:,me))
+  sumGamma = sum(stt%gamma_slip(:,en))
-  sumF     = sum(stt%gamma_twin(:,me)/prm%gamma_char)
+  sumF     = sum(stt%gamma_twin(:,en)/prm%gamma_char)
 
 !--------------------------------------------------------------------------------------------------
 ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
@@ -366,23 +366,23 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me)
 !  calculate left and right vectors
   left_SlipSlip  = 1.0_pReal + prm%h_int
   xi_slip_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
-  right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,me) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl &
+  right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl &
-                 * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,me) / (prm%xi_inf_sl+xi_slip_sat_offset))
+                 * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset))
 
 !--------------------------------------------------------------------------------------------------
 ! shear rates
-  call kinetics_slip(Mp,ph,me,gdot_slip_pos,gdot_slip_neg)
+  call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg)
-  dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg)
+  dot%gamma_slip(:,en) = abs(gdot_slip_pos+gdot_slip_neg)
-  call kinetics_twin(Mp,ph,me,dot%gamma_twin(:,me))
+  call kinetics_twin(Mp,ph,en,dot%gamma_twin(:,en))
 
 !--------------------------------------------------------------------------------------------------
 ! hardening
-  dot%xi_slip(:,me) = c_SlipSlip * left_SlipSlip * &
+  dot%xi_slip(:,en) = c_SlipSlip * left_SlipSlip * &
-                      matmul(prm%h_sl_sl,dot%gamma_slip(:,me)*right_SlipSlip) &
+                      matmul(prm%h_sl_sl,dot%gamma_slip(:,en)*right_SlipSlip) &
-                    + matmul(prm%h_sl_tw,dot%gamma_twin(:,me))
+                    + matmul(prm%h_sl_tw,dot%gamma_twin(:,en))
 
-  dot%xi_twin(:,me) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,me)) &
+  dot%xi_twin(:,en) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,en)) &
-                    + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,me))
+                    + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,en))
   end associate
 
 end subroutine phenopowerlaw_dotState
@@ -430,14 +430,14 @@ end subroutine plastic_phenopowerlaw_results
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_slip(Mp,ph,me, &
+pure subroutine kinetics_slip(Mp,ph,en, &
                               gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal),                  intent(out), dimension(param(ph)%sum_N_sl) :: &
     gdot_slip_pos, &
@@ -461,14 +461,14 @@ pure subroutine kinetics_slip(Mp,ph,me, &
 
   where(dNeq0(tau_slip_pos))
     gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &          ! 1/2 if non-Schmid active
-                  * sign(abs(tau_slip_pos/stt%xi_slip(:,me))**prm%n_sl,  tau_slip_pos)
+                  * sign(abs(tau_slip_pos/stt%xi_slip(:,en))**prm%n_sl,  tau_slip_pos)
   else where
     gdot_slip_pos = 0.0_pReal
   end where
 
   where(dNeq0(tau_slip_neg))
     gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal &                                                ! only used if non-Schmid active, always 1/2
-                  * sign(abs(tau_slip_neg/stt%xi_slip(:,me))**prm%n_sl,  tau_slip_neg)
+                  * sign(abs(tau_slip_neg/stt%xi_slip(:,en))**prm%n_sl,  tau_slip_neg)
   else where
     gdot_slip_neg = 0.0_pReal
   end where
@@ -499,14 +499,14 @@ end subroutine kinetics_slip
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_twin(Mp,ph,me,&
+pure subroutine kinetics_twin(Mp,ph,en,&
                               gdot_twin,dgdot_dtau_twin)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     ph, &
-    me
+    en
 
   real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
     gdot_twin
@@ -524,8 +524,8 @@ pure subroutine kinetics_twin(Mp,ph,me,&
   enddo
 
   where(tau_twin > 0.0_pReal)
-    gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,me)/prm%gamma_char)) &                              ! only twin in untwinned volume fraction
+    gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,en)/prm%gamma_char)) &                              ! only twin in untwinned volume fraction
-              * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,me))**prm%n_tw
+              * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,en))**prm%n_tw
   else where
     gdot_twin = 0.0_pReal
   end where

src/phase_thermal.f90

@@ -26,7 +26,7 @@ submodule(phase) thermal
   integer(kind(THERMAL_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
     thermal_source
 
-  type(tDataContainer), dimension(:), allocatable :: current          ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(me) reads quite good
+  type(tDataContainer), dimension(:), allocatable :: current          ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(en) reads quite good
 
   type(tThermalParameters), dimension(:), allocatable :: param
 
@@ -46,23 +46,23 @@ submodule(phase) thermal
     end function externalheat_init
 
 
-    module subroutine externalheat_dotState(ph, me)
+    module subroutine externalheat_dotState(ph, en)
       integer, intent(in) :: &
         ph, &
-        me
+        en
     end subroutine externalheat_dotState
 
-    module function dissipation_f_T(ph,me) result(f_T)
+    module function dissipation_f_T(ph,en) result(f_T)
       integer, intent(in) :: &
         ph, &
-        me
+        en
       real(pReal) :: f_T
     end function dissipation_f_T
 
-    module function externalheat_f_T(ph,me)  result(f_T)
+    module function externalheat_f_T(ph,en)  result(f_T)
       integer, intent(in) :: &
         ph, &
-        me
+        en
       real(pReal) :: f_T
     end function externalheat_f_T
 
@@ -137,9 +137,9 @@ end subroutine thermal_init
 !----------------------------------------------------------------------------------------------
 !< @brief calculates thermal dissipation rate
 !----------------------------------------------------------------------------------------------
-module function phase_f_T(ph,me) result(f)
+module function phase_f_T(ph,en) result(f)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   real(pReal) :: f
 
 
@@ -152,10 +152,10 @@ module function phase_f_T(ph,me) result(f)
    select case(thermal_source(so,ph))
 
      case (THERMAL_DISSIPATION_ID)
-       f = f + dissipation_f_T(ph,me)
+       f = f + dissipation_f_T(ph,en)
 
      case (THERMAL_EXTERNALHEAT_ID)
-       f = f + externalheat_f_T(ph,me)
+       f = f + externalheat_f_T(ph,en)
 
    end select
 
@@ -167,9 +167,9 @@ end function phase_f_T
 !--------------------------------------------------------------------------------------------------
 !> @brief contains the constitutive equation for calculating the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
-function phase_thermal_collectDotState(ph,me) result(broken)
+function phase_thermal_collectDotState(ph,en) result(broken)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   logical :: broken
 
   integer :: i
@@ -180,9 +180,9 @@ function phase_thermal_collectDotState(ph,me) result(broken)
   SourceLoop: do i = 1, thermal_Nsources(ph)
 
     if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) &
-      call externalheat_dotState(ph,me)
+      call externalheat_dotState(ph,en)
 
-    broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me)))
+    broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))
 
   enddo SourceLoop
 
@@ -218,14 +218,14 @@ module function phase_K_T(co,ce) result(K)
 end function phase_K_T
 
 
-module function thermal_stress(Delta_t,ph,me) result(converged_)           ! ?? why is this called "stress" when it seems closer to "updateState" ??
+module function thermal_stress(Delta_t,ph,en) result(converged_)           ! ?? why is this called "stress" when it seems closer to "updateState" ??
 
   real(pReal), intent(in) :: Delta_t
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   logical :: converged_
 
 
-  converged_ = .not. integrateThermalState(Delta_t,ph,me)
+  converged_ = .not. integrateThermalState(Delta_t,ph,en)
 
 end function thermal_stress
 
@@ -233,10 +233,10 @@ end function thermal_stress
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate state with 1st order explicit Euler method
 !--------------------------------------------------------------------------------------------------
-function integrateThermalState(Delta_t, ph,me) result(broken)
+function integrateThermalState(Delta_t, ph,en) result(broken)
 
   real(pReal), intent(in) :: Delta_t
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   logical :: &
     broken
 
@@ -244,13 +244,13 @@ function integrateThermalState(Delta_t, ph,me) result(broken)
     so, &
     sizeDotState
 
-  broken = phase_thermal_collectDotState(ph,me)
+  broken = phase_thermal_collectDotState(ph,en)
   if (broken) return
 
   do so = 1, thermal_Nsources(ph)
     sizeDotState = thermalState(ph)%p(so)%sizeDotState
-    thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%state0(1:sizeDotState,me) &
+    thermalState(ph)%p(so)%state(1:sizeDotState,en) = thermalState(ph)%p(so)%state0(1:sizeDotState,en) &
-                                                    + thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t
+                                                    + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t
   enddo
 
 end function integrateThermalState
@@ -273,13 +273,13 @@ end subroutine thermal_forward
 !----------------------------------------------------------------------------------------------
 !< @brief Get temperature (for use by non-thermal physics)
 !----------------------------------------------------------------------------------------------
-module function thermal_T(ph,me) result(T)
+module function thermal_T(ph,en) result(T)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   real(pReal) :: T
 
 
-  T = current(ph)%T(me)
+  T = current(ph)%T(en)
 
 end function thermal_T
 
@@ -287,13 +287,13 @@ end function thermal_T
 !----------------------------------------------------------------------------------------------
 !< @brief Get rate of temperature (for use by non-thermal physics)
 !----------------------------------------------------------------------------------------------
-module function thermal_dot_T(ph,me) result(dot_T)
+module function thermal_dot_T(ph,en) result(dot_T)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   real(pReal) :: dot_T
 
 
-  dot_T = current(ph)%dot_T(me)
+  dot_T = current(ph)%dot_T(en)
 
 end function thermal_dot_T
 

src/phase_thermal_dissipation.f90

@@ -69,15 +69,15 @@ end function dissipation_init
 !--------------------------------------------------------------------------------------------------
 !> @brief Ninstancess dissipation rate
 !--------------------------------------------------------------------------------------------------
-module function dissipation_f_T(ph,me) result(f_T)
+module function dissipation_f_T(ph,en) result(f_T)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   real(pReal) :: &
     f_T
 
 
   associate(prm => param(ph))
-    f_T = prm%kappa*sum(abs(mechanical_S(ph,me)*mechanical_L_p(ph,me)))
+    f_T = prm%kappa*sum(abs(mechanical_S(ph,en)*mechanical_L_p(ph,en)))
   end associate
 
 end function dissipation_f_T

src/phase_thermal_externalheat.f90

@@ -81,18 +81,18 @@ end function externalheat_init
 !> @brief rate of change of state
 !> @details state only contains current time to linearly interpolate given heat powers
 !--------------------------------------------------------------------------------------------------
-module subroutine externalheat_dotState(ph, me)
+module subroutine externalheat_dotState(ph, en)
 
   integer, intent(in) :: &
     ph, &
-    me
+    en
 
   integer :: &
     so
 
   so = source_thermal_externalheat_offset(ph)
 
-  thermalState(ph)%p(so)%dotState(1,me) = 1.0_pReal                                                 ! state is current time
+  thermalState(ph)%p(so)%dotState(1,en) = 1.0_pReal                                                 ! state is current time
 
 end subroutine externalheat_dotState
 
@@ -100,11 +100,11 @@ end subroutine externalheat_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local heat generation rate
 !--------------------------------------------------------------------------------------------------
-module function externalheat_f_T(ph,me) result(f_T)
+module function externalheat_f_T(ph,en) result(f_T)
 
   integer, intent(in) :: &
     ph, &
-    me
+    en
   real(pReal) :: &
     f_T
 
@@ -117,14 +117,14 @@ module function externalheat_f_T(ph,me) result(f_T)
 
   associate(prm => param(ph))
     do interval = 1, prm%nIntervals                                                                 ! scan through all rate segments
-      frac_time = (thermalState(ph)%p(so)%state(1,me) - prm%t_n(interval)) &
+      frac_time = (thermalState(ph)%p(so)%state(1,en) - prm%t_n(interval)) &
                 / (prm%t_n(interval+1) - prm%t_n(interval))                                         ! fractional time within segment
       if (     (frac_time <  0.0_pReal .and. interval == 1) &
           .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
           .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
         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 me bounds
+                                                                                                   ! ...or extrapolate if outside of bounds
     enddo
   end associate