Major Update: all modules are now correctly submitted

code/CPFEM.f90

@@ -121,8 +121,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
   integer(pInt), intent(in) ::                        element, &          ! FE element number
                                                       IP, &               ! FE integration point number
                                                       ngens               ! size of stress strain law
-  real(pReal), intent(in) ::                          Temperature, &      ! temperature
+  real(pReal), intent(inout) ::                       Temperature         ! temperature
-                                                      dt                  ! time increment
+  real(pReal), intent(in) ::                          dt                  ! time increment
   real(pReal), dimension (3,3), intent(in) ::         ffn, &              ! deformation gradient for t=t0
                                                       ffn1                ! deformation gradient for t=t1
   integer(pInt), intent(in) ::                        mode                ! computation mode  1: regular computation with aged results
@@ -170,8 +170,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
     call lattice_init()
     call material_init()
     call constitutive_init()
-    call crystallite_init()
+    call crystallite_init(Temperature)         ! (have to) use temperature of first IP for whole model
-    call homogenization_init()
+    call homogenization_init(Temperature)
     call CPFEM_init()
     CPFEM_init_done = .true.
   endif
@@ -262,6 +262,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
     
     ! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+-- 
     case (3)
+	  if (IP==1.AND.cp_en==1) write(6,*) 'Temp from CPFEM', Temperature
       materialpoint_Temperature(IP,cp_en)   = Temperature
       materialpoint_F0(:,:,IP,cp_en)        = ffn
       materialpoint_F(:,:,IP,cp_en)         = ffn1
@@ -287,7 +288,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
   cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en)
   jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,1:ngens,IP,cp_en)
   ! return temperature
-  Temperature = materialpoint_Temperature(IP,cp_en)
+  if (theInc > 0_pInt) Temperature = materialpoint_Temperature(IP,cp_en)  ! homogenized result except for potentially non-isothermal starting condition.
   
   return
 

code/constitutive_dislobased.f90

@@ -373,7 +373,6 @@ subroutine constitutive_dislobased_microstructure(Temperature,state,ipc,ip,el)
  real(pReal) Temperature
  type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
 
- Temperature = 298.0
  matID = phase_constitutionInstance(material_phase(ipc,ip,el))
  n = constitutive_dislobased_Nslip(matID)
  !* Quantities derived from state - slip
@@ -448,7 +447,6 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera
  real(pReal), dimension(constitutive_dislobased_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
    gdot_slip,dgdot_dtauslip,tau_slip
 
- Temperature = 298.0
  matID = phase_constitutionInstance(material_phase(ipc,ip,el))
  n = constitutive_dislobased_Nslip(matID)
 
@@ -512,7 +510,6 @@ function constitutive_dislobased_dotState(Tstar_v,Temperature,state,ipc,ip,el)
  real(pReal), dimension(constitutive_dislobased_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
    constitutive_dislobased_dotState
 
- Temperature = 298.0
  matID = phase_constitutionInstance(material_phase(ipc,ip,el))
  n = constitutive_dislobased_Nslip(matID)
 
@@ -605,7 +602,7 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i
    select case(constitutive_dislobased_output(o,matID))
 
      case ('dislodensity')
-       constitutive_dislobased_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n)
+       constitutive_dislobased_postResults(c+1:c+n) = state(ipc,ip,el)%p(6*n+1:7*n)
        c = c + n
 
      case ('rateofshear')
@@ -613,7 +610,7 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i
          tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID)))
 		 if ((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))>0) then
             constitutive_dislobased_postResults(c+i) = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip)*&
-            sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*298.0))
+            sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature))
 	     else
 		    constitutive_dislobased_postResults(c+i) = 0.0_pReal
 		 endif

code/constitutive_j2.f90

@@ -417,6 +417,36 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
 endfunction
 
 
+!****************************************************************
+!* calculates the rate of change of microstructure              *
+!****************************************************************
+pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, el)
+
+  !*** variables and functions from other modules ***!
+  use prec,     only: pReal,pInt,p_vec
+  use lattice,  only: lattice_Sslip_v
+  use mesh,     only: mesh_NcpElems,mesh_maxNips
+  use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance  
+  implicit none
+
+  !*** input variables ***!
+  real(pReal), dimension(6), intent(in) ::  Tstar_v                   ! 2nd Piola Kirchhoff stress tensor in Mandel notation
+  real(pReal), intent(in) ::                Temperature
+  integer(pInt), intent(in)::               g, &                      ! grain number
+                                            ip, &                     ! integration point number
+                                            el                        ! element number
+  type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure
+  
+  !*** output variables ***!
+  real(pReal) constitutive_j2_dotTemperature  ! evolution of state variable
+  
+  ! calculate dotTemperature
+  constitutive_j2_dotTemperature = 0.0_pReal
+
+  return
+endfunction
+
+
 !*********************************************************************
 !* return array of constitutive results                              *
 !*********************************************************************

code/constitutive_phenomenological.f90

@@ -355,8 +355,10 @@ subroutine constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,T
  Lp = 0.0_pReal
  do i = 1,constitutive_phenomenological_Nslip(matID)
    tau_slip(i)  = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
+
    gdot_slip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**&
                   constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip(i))
+   
    Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_phenomenological_structure(matID))
  enddo
 
@@ -364,9 +366,11 @@ subroutine constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,T
  dLp_dTstar3333 = 0.0_pReal
  dLp_dTstar = 0.0_pReal
  do i = 1,constitutive_phenomenological_Nslip(matID)
+
    dgdot_dtauslip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**&
                       (constitutive_phenomenological_n_slip(matID)-1.0_pReal)*&
                       constitutive_phenomenological_n_slip(matID)/state(ipc,ip,el)%p(i)
+
    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
           dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
              dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_phenomenological_structure(matID))* &
@@ -409,9 +413,12 @@ function constitutive_phenomenological_dotState(Tstar_v,Temperature,state,ipc,ip
 
 !* Self-Hardening of each system
  do i = 1,n
+
    tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
+
    gdot_slip = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip)/state(ipc,ip,el)%p(i))**&
                constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip)
+
    self_hardening(i) = constitutive_phenomenological_h0(matID)*(1.0_pReal-state(ipc,ip,el)%p(i)/&
                        constitutive_phenomenological_s_sat(matID))**constitutive_phenomenological_w0(matID)*abs(gdot_slip)
  enddo
@@ -424,6 +431,37 @@ function constitutive_phenomenological_dotState(Tstar_v,Temperature,state,ipc,ip
 end function
 
 
+function constitutive_phenomenological_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el)
+!*********************************************************************
+!* rate of change of microstructure                                  *
+!* INPUT:                                                            *
+!*  - Tstar_v         : 2nd Piola Kirchhoff stress tensor (Mandel)   *
+!*  - ipc             : component-ID at current integration point    *
+!*  - ip              : current integration point                    *
+!*  - el              : current element                              *
+!* OUTPUT:                                                           *
+!*  - constitutive_dotTemperature : evolution of temperature         *
+!*********************************************************************
+ use prec, only: pReal,pInt,p_vec
+ use lattice, only: lattice_Sslip_v
+ use mesh, only: mesh_NcpElems,mesh_maxNips
+ use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
+ implicit none
+
+!* Definition of variables
+ integer(pInt) ipc,ip,el
+ integer(pInt) matID,i,n
+ real(pReal) Temperature
+ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
+ real(pReal), dimension(6) :: Tstar_v
+ real(pReal) constitutive_phenomenological_dotTemperature
+
+ constitutive_phenomenological_dotTemperature = 0.0_pReal
+
+ return
+end function
+
+
 pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
 !*********************************************************************
 !* return array of constitutive results                              *
@@ -457,9 +495,11 @@ pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,s
 
  do o = 1,phase_Noutput(material_phase(ipc,ip,el))
    select case(constitutive_phenomenological_output(o,matID))
+
      case ('slipresistance')
        constitutive_phenomenological_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n)
        c = c + n
+
      case ('rateofshear')
        do i = 1,n
          tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
@@ -468,6 +508,7 @@ pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,s
                                                            constitutive_phenomenological_n_slip(matID)
        enddo
        c = c + n
+
    end select
  enddo
  

code/crystallite.f90

@@ -26,7 +26,6 @@ real(pReal), dimension (:,:,:), allocatable ::          crystallite_dt, &
                                                         crystallite_subFrac, &               ! already calculated fraction of increment
                                                         crystallite_subStep, &               ! size of next integration step
                                                         crystallite_Temperature, &           ! Temp of each grain 
-                                                        crystallite_Temperature0, &         ! Temp of each grain at start of FE inc
                                                         crystallite_partionedTemperature0, & ! Temp of each grain at start of homog inc
                                                         crystallite_subTemperature0          ! Temp of each grain at start of crystallite inc
 real(pReal), dimension (:,:,:,:), allocatable ::        crystallite_Tstar_v, &               ! current 2nd Piola-Kirchhoff stress vector (end of converged time step)
@@ -62,7 +61,7 @@ CONTAINS
 !********************************************************************
 ! allocate and initialize per grain variables
 !********************************************************************
-subroutine crystallite_init()
+subroutine crystallite_init(Temperature)
 
   !*** variables and functions from other modules ***!
   use prec, only:             pInt, &
@@ -84,6 +83,7 @@ subroutine crystallite_init()
   implicit none
 
   !*** input variables ***!
+  real(pReal) Temperature
  
   !*** output variables ***!
  
@@ -100,13 +100,12 @@ subroutine crystallite_init()
   iMax = mesh_maxNips
   eMax = mesh_NcpElems
 
-  allocate(crystallite_Temperature(gMax,iMax,eMax));                crystallite_Temperature = 0.0_pReal
+  allocate(crystallite_Temperature(gMax,iMax,eMax));                crystallite_Temperature = Temperature
   allocate(crystallite_P(3,3,gMax,iMax,eMax));                                crystallite_P = 0.0_pReal
   allocate(crystallite_Fe(3,3,gMax,iMax,eMax));                              crystallite_Fe = 0.0_pReal
   allocate(crystallite_Fp(3,3,gMax,iMax,eMax));                              crystallite_Fp = 0.0_pReal
   allocate(crystallite_Lp(3,3,gMax,iMax,eMax));                              crystallite_Lp = 0.0_pReal
   allocate(crystallite_Tstar_v(6,gMax,iMax,eMax));                      crystallite_Tstar_v = 0.0_pReal
-  allocate(crystallite_Temperature0(gMax,iMax,eMax));              crystallite_Temperature0 = 0.0_pReal
   allocate(crystallite_F0(3,3,gMax,iMax,eMax));                              crystallite_F0 = 0.0_pReal
   allocate(crystallite_Fp0(3,3,gMax,iMax,eMax));                            crystallite_Fp0 = 0.0_pReal
   allocate(crystallite_Lp0(3,3,gMax,iMax,eMax));                            crystallite_Lp0 = 0.0_pReal
@@ -139,7 +138,7 @@ subroutine crystallite_init()
       myNgrains = homogenization_Ngrains(mesh_element(3,e))
       do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)         ! iterate over IPs of this element
         do g = 1,myNgrains
-          crystallite_partionedTemperature0(g,i,e) = crystallite_Temperature0(g,i,e)
+          crystallite_partionedTemperature0(g,i,e) = Temperature                   ! isothermal assumption
           crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e))    ! plastic def gradient reflects init orientation
           crystallite_F0(:,:,g,i,e)  = math_I3
           crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e)
@@ -166,7 +165,6 @@ subroutine crystallite_init()
     write(6,'(a32,x,7(i5,x))') 'crystallite_Fe:                ', shape(crystallite_Fe)
     write(6,'(a32,x,7(i5,x))') 'crystallite_Fp:                ', shape(crystallite_Fp)
     write(6,'(a32,x,7(i5,x))') 'crystallite_Lp:                ', shape(crystallite_Lp)
-    write(6,'(a32,x,7(i5,x))') 'crystallite_Temperature0:      ', shape(crystallite_Temperature0)
     write(6,'(a32,x,7(i5,x))') 'crystallite_F0:                ', shape(crystallite_F0)
     write(6,'(a32,x,7(i5,x))') 'crystallite_Fp0:               ', shape(crystallite_Fp0)
     write(6,'(a32,x,7(i5,x))') 'crystallite_Lp0:               ', shape(crystallite_Lp0)
@@ -281,11 +279,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
 
   write (6,*)
   write (6,*) 'Crystallite request from Materialpoint'
-  write (6,'(a,/,(f12.7,x))') 'crystallite_partionedTemperature0  of 1 1 1',crystallite_partionedTemperature0(1,1,1)
+  write (6,'(a,/,(f12.7,x))')      'crystallite_partionedTemp0 of 1 1 1' ,crystallite_partionedTemperature0(1,1,1)
-  write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0  of 1 1 1',crystallite_partionedF0(1:3,:,1,1,1)
+  write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0  of 1 1 1'   ,crystallite_partionedF0(1:3,:,1,1,1)
-  write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1',crystallite_partionedFp0(1:3,:,1,1,1)
+  write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1'   ,crystallite_partionedFp0(1:3,:,1,1,1)
-  write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF   of 1 1 1',crystallite_partionedF(1:3,:,1,1,1)
+  write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF   of 1 1 1'   ,crystallite_partionedF(1:3,:,1,1,1)
-  write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1',crystallite_partionedLp0(1:3,:,1,1,1)
+  write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1'   ,crystallite_partionedLp0(1:3,:,1,1,1)
 
  
   !$OMP PARALLEL DO
@@ -456,8 +454,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
               if (             crystallite_requested(g,i,e) & 
                   .and.       crystallite_onTrack(g,i,e) &
                   .and. .not. crystallite_converged(g,i,e)) then                      ! all undone crystallites
-                crystallite_converged(g,i,e) = crystallite_updateState(g,i,e).AND.&
+                crystallite_converged(g,i,e) = (crystallite_updateState(g,i,e).AND.&
-				                               crystallite_updateTemperature(g,i,e)				                                 
+				                               crystallite_updateTemperature(g,i,e))			                                 
                 if (crystallite_converged(g,i,e)) then
                   !$OMP CRITICAL (distributionState)
                     debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1
@@ -724,8 +722,7 @@ endsubroutine
  crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum
  
  ! setting flag to true if state is below relative Tolerance, otherwise set it to false
- crystallite_updateTemperature = maxval(abs(residuum/crystallite_Temperature(g,i,e)), &
+ crystallite_updateTemperature = abs(residuum/crystallite_Temperature(g,i,e)) < rTol_crystalliteTemperature
-                                  crystallite_Temperature(g,i,e) /= 0.0_pReal) < rTol_crystalliteTemperature
                                   
  if (debugger) write(6,'(a,/,f12.4)') 'updated temperature: ', crystallite_Temperature(g,i,e)
  

code/homogenization.f90

@@ -51,7 +51,7 @@ CONTAINS
 !**************************************
 !*      Module initialization         *
 !**************************************
-subroutine homogenization_init()
+subroutine homogenization_init(Temperature)
  use prec, only: pReal,pInt
  use math, only: math_I3
  use IO, only: IO_error, IO_open_file
@@ -62,6 +62,7 @@ subroutine homogenization_init()
  use homogenization_isostrain
 ! use homogenization_RGC
 
+ real(pReal) Temperature
  integer(pInt), parameter :: fileunit = 200
  integer(pInt) e,i,g,myInstance
 
@@ -83,7 +84,7 @@ subroutine homogenization_init()
  allocate(materialpoint_subF0(3,3,mesh_maxNips,mesh_NcpElems));        materialpoint_subF0   = 0.0_pReal
  allocate(materialpoint_subF(3,3,mesh_maxNips,mesh_NcpElems));         materialpoint_subF    = 0.0_pReal
  allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems));            materialpoint_P       = 0.0_pReal
- allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems));      materialpoint_Temperature = 0.0_pReal
+ allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems));      materialpoint_Temperature = Temperature
  allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems));          materialpoint_subFrac = 0.0_pReal
  allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems));          materialpoint_subStep = 0.0_pReal
  allocate(materialpoint_subdt(mesh_maxNips,mesh_NcpElems));            materialpoint_subdt   = 0.0_pReal
@@ -180,8 +181,7 @@ subroutine materialpoint_stressAndItsTangent(&
  use constitutive, only:  constitutive_state0, &
                           constitutive_partionedState0, &
                           constitutive_state
- use crystallite, only:   crystallite_Temperature0, &
+ use crystallite, only:   crystallite_Temperature, &
-                          crystallite_Temperature, &
                           crystallite_F0, &
                           crystallite_Fp0, &
                           crystallite_Fp, &
@@ -210,6 +210,7 @@ subroutine materialpoint_stressAndItsTangent(&
 
  write (6,*)
  write (6,*) 'Material Point start'
+ write (6,'(a,/,(f12.7,x))') 'Temp0  of   8 1'  ,materialpoint_Temperature(8,1)
  write (6,'(a,/,3(3(f12.7,x)/))') 'F0  of   8 1',materialpoint_F0(1:3,:,8,1)
  write (6,'(a,/,3(3(f12.7,x)/))') 'F   of   8 1',materialpoint_F(1:3,:,8,1)
  write (6,'(a,/,3(3(f12.7,x)/))') 'Fp0 of 1 8 1',crystallite_Fp0(1:3,:,1,8,1)
@@ -222,7 +223,7 @@ subroutine materialpoint_stressAndItsTangent(&
      
      ! initialize restoration points of grain...
      forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p   ! ...microstructures
-     crystallite_partionedTemperature0(1:myNgrains,i,e)   = crystallite_Temperature0(1:myNgrains,i,e)! ...temperatures
+     crystallite_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e)             ! ...temperatures
      crystallite_partionedFp0(:,:,1:myNgrains,i,e)   = crystallite_Fp0(:,:,1:myNgrains,i,e)          ! ...plastic def grads
      crystallite_partionedLp0(:,:,1:myNgrains,i,e)   = crystallite_Lp0(:,:,1:myNgrains,i,e)          ! ...plastic velocity grads
      crystallite_partionedF0(:,:,1:myNgrains,i,e)    = crystallite_F0(:,:,1:myNgrains,i,e)           ! ...def grads

code/homogenization_isostrain.f90

@@ -257,7 +257,7 @@ function homogenization_isostrain_averageTemperature(&
 
 ! homID = homogenization_typeInstance(mesh_element(3,el))
  Ngrains = homogenization_Ngrains(mesh_element(3,el))
- homogenization_isostrain_averageTemperature = sum(Temperature)/Ngrains
+ homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/Ngrains
 
  return
 

code/numerics.f90

@@ -17,6 +17,7 @@ real(pReal)                     relevantStrain, &                       ! strain
                                 pert_Fg, &                              ! strain perturbation for FEM Jacobi
                                 subStepMin, &                           ! minimum (relative) size of sub-step allowed during cutback in crystallite
                                 rTol_crystalliteState, &                ! relative tolerance in crystallite state loop 
+                                rTol_crystalliteTemperature, &          ! relative tolerance in crystallite temperature loop 
                                 rTol_crystalliteStress, &               ! relative tolerance in crystallite stress loop
                                 aTol_crystalliteStress, &               ! absolute tolerance in crystallite stress loop
                                 resToler, &                             ! relative tolerance of residual in GIA iteration
@@ -89,6 +90,7 @@ subroutine numerics_init()
   nStress                 = 40_pInt
   subStepMin              = 1.0e-3_pReal
   rTol_crystalliteState   = 1.0e-6_pReal
+  rTol_crystalliteTemperature = 1.0e-6_pReal
   rTol_crystalliteStress  = 1.0e-6_pReal
   aTol_crystalliteStress  = 1.0e-8_pReal
   resToler                = 1.0e-4_pReal
@@ -130,6 +132,8 @@ subroutine numerics_init()
               subStepMin = IO_floatValue(line,positions,2)
         case ('rtol_crystallitestate')
               rTol_crystalliteState = IO_floatValue(line,positions,2)
+        case ('rtol_crystalliteTemperature')
+              rTol_crystalliteTemperature = IO_floatValue(line,positions,2)
         case ('rtol_crystallitestress')
               rTol_crystalliteStress = IO_floatValue(line,positions,2)
         case ('atol_crystallitestress')
@@ -165,6 +169,7 @@ subroutine numerics_init()
   write(6,'(a24,x,i8)')   'nStress:                ',nStress
   write(6,'(a24,x,e8.1)') 'subStepMin:             ',subStepMin
   write(6,'(a24,x,e8.1)') 'rTol_crystalliteState:  ',rTol_crystalliteState
+  write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp:   ',rTol_crystalliteTemperature
   write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress
   write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress
   write(6,'(a24,x,e8.1)') 'resToler:               ',resToler
@@ -173,7 +178,7 @@ subroutine numerics_init()
   write(6,'(a24,x,i8)')   'NRiterMax:              ',NRiterMax
   write(6,*)
   
-  ! sanity check  
+  ! sanity check  (Temperature check missing!!!!!!!)
   if (relevantStrain <= 0.0_pReal)          call IO_error(260)
   if (iJacoStiffness < 1_pInt)              call IO_error(261)
   if (iJacoLpresiduum < 1_pInt)             call IO_error(262)