explicit arguments instead of global variables

src/constitutive.f90

@@ -42,8 +42,7 @@ module constitutive
     KINEMATICS_SLIPPLANE_OPENING_ID, &
     KINEMATICS_THERMAL_EXPANSION_ID
   end enum
-  real(pReal),               dimension(:,:,:),        allocatable :: &
+
-    crystallite_subdt                                                                               !< substepped time increment of each grain
   type(rotation),            dimension(:,:,:),        allocatable :: &
     crystallite_orientation                                                                         !< current orientation
   real(pReal),               dimension(:,:,:,:,:),    allocatable :: &
@@ -874,7 +873,6 @@ subroutine crystallite_init
            crystallite_Fe,crystallite_Lp, &
            source = crystallite_F)
 
-  allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal)
   allocate(crystallite_orientation(cMax,iMax,eMax))
 
 
@@ -1103,11 +1101,11 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
     lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
     do o=1,3; do p=1,3
       lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
-                            + crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
+                            + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * dt
       lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
                             + invFi*invFi(p,o)
       rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
-                            - crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
+                            - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * dt
     enddo; enddo
     call math_invert(temp_99,error,math_3333to99(lhs_3333))
     if (error) then
@@ -1136,7 +1134,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
     temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
                            + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
   enddo; enddo
-  lhs_3333 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) &
+  lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * dt &
            + math_mul3333xx3333(dSdFi,dFidS)
 
   call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
@@ -1152,8 +1150,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
 ! calculate dFpinvdF
   temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
   do o=1,3; do p=1,3
-    dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) &
+    dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt
-                          * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
   enddo; enddo
 
 !--------------------------------------------------------------------------------------------------

src/constitutive_mech.f90

@@ -1526,7 +1526,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   todo = .true.
   converged_ = .false.                                                      ! pretend failed step of 1/subStepSizeCryst
 
-  crystallite_subdt(co,ip,el) = dt
   todo = .true.
   NiterationCrystallite = 0
   cutbackLooping: do while (todo)

src/homogenization.f90

@@ -63,7 +63,8 @@ module homogenization
         el                                                                                          !< element number
     end subroutine mech_partition
 
-    module subroutine mech_homogenize(ip,el)
+    module subroutine mech_homogenize(dt,ip,el)
+     real(pReal), intent(in) :: dt
      integer, intent(in) :: &
        ip, &                                                                                        !< integration point
        el                                                                                           !< element number
@@ -257,7 +258,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
         do co = 1, myNgrains
           call crystallite_orientations(co,ip,el)
         enddo
-        call mech_homogenize(ip,el)
+        call mech_homogenize(dt,ip,el)
       enddo IpLooping3
     enddo elementLooping3
     !$OMP END PARALLEL DO

src/homogenization_mech.f90

@@ -138,11 +138,13 @@ end subroutine mech_partition
 !--------------------------------------------------------------------------------------------------
 !> @brief Average P and dPdF from the individual constituents.
 !--------------------------------------------------------------------------------------------------
-module subroutine mech_homogenize(ip,el)
+module subroutine mech_homogenize(dt,ip,el)
 
+  real(pReal), intent(in) :: dt
   integer, intent(in) :: &
        ip, &                                                                                        !< integration point
        el                                                                                           !< element number
+
   integer :: co,ce
   real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
 
@@ -152,11 +154,11 @@ module subroutine mech_homogenize(ip,el)
 
     case (HOMOGENIZATION_NONE_ID) chosenHomogenization
         homogenization_P(1:3,1:3,ce)            = crystallite_P(1:3,1:3,1,ip,el)
-        homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el)
+        homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(dt,1,ip,el)
 
     case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
       do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
-        dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
+        dPdFs(:,:,:,:,co) = crystallite_stressTangent(dt,co,ip,el)
       enddo
       call mech_isostrain_averageStressAndItsTangent(&
         homogenization_P(1:3,1:3,ce), &
@@ -167,7 +169,7 @@ module subroutine mech_homogenize(ip,el)
 
     case (HOMOGENIZATION_RGC_ID) chosenHomogenization
       do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
-        dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
+        dPdFs(:,:,:,:,co) = crystallite_stressTangent(dt,co,ip,el)
       enddo
       call mech_RGC_averageStressAndItsTangent(&
         homogenization_P(1:3,1:3,ce), &
@@ -202,7 +204,7 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
 
   if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
       do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
-        dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
+        dPdFs(:,:,:,:,co) = crystallite_stressTangent(subdt,co,ip,el)
       enddo
       doneAndHappy = &
           mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &