diff --git a/src/phase.f90 b/src/phase.f90
index 8ea1fa5f7..eb6e61eb8 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -58,12 +58,9 @@ module phase
   type(tDebugOptions) :: debugCrystallite
 
   integer, dimension(:), allocatable, public :: &                                                   !< ToDo: should be protected (bug in Intel compiler)
-    thermal_Nsources, &
-    phase_Nsources, &                                                                               !< number of source mechanisms active in each phase
-    phase_Nkinematics, &                                                                            !< number of kinematic mechanisms active in each phase
+    phase_elasticityInstance, &
     phase_NstiffnessDegradations, &                                                                 !< number of stiffness degradation mechanisms active in each phase
-    phase_plasticInstance, &                                                                     !< instance of particular plasticity of each phase
-    phase_elasticityInstance                                                                        !< instance of particular elasticity of each phase
+    phase_plasticInstance
 
   logical, dimension(:), allocatable, public :: &                                                   ! ToDo: should be protected (bug in Intel Compiler)
     phase_localPlasticity                                                                           !< flags phases with local constitutive law
@@ -351,7 +348,7 @@ subroutine phase_init
 !--------------------------------------------------------------------------------------------------
 ! partition and initialize state
     plasticState(ph)%state = plasticState(ph)%state0
-    if(phase_Nsources(ph) > 0) &
+    if(damageState(ph)%sizeState > 0) &
       damageState(ph)%state  = damageState(ph)%state0
   enddo PhaseLoop2
 
@@ -365,7 +362,7 @@ end subroutine phase_init
 !> @brief Allocate the components of the state structure for a given phase
 !--------------------------------------------------------------------------------------------------
 subroutine phase_allocateState(state, &
-                                  Nconstituents,sizeState,sizeDotState,sizeDeltaState)
+                               Nconstituents,sizeState,sizeDotState,sizeDeltaState)
 
   class(tState), intent(out) :: &
     state
@@ -406,7 +403,7 @@ subroutine phase_restore(ce,includeL)
 
 
   do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce))
-    if (phase_Nsources(material_phaseAt2(co,ce)) > 0) &
+    if (damageState(material_phaseAt2(co,ce))%sizeState > 0) &
     damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = &
       damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce))
   enddo
@@ -429,7 +426,7 @@ subroutine phase_forward()
   call thermal_forward()
 
   do ph = 1, size(damageState)
-    if (phase_Nsources(ph) > 0) &
+    if (damageState(ph)%sizeState > 0) &
       damageState(ph)%state0 = damageState(ph)%state
   enddo
 
@@ -527,7 +524,7 @@ subroutine crystallite_init()
   phases => config_material%get('phase')
 
   do ph = 1, phases%length
-    if (phase_Nsources(ph) > 0) &
+    if (damageState(ph)%sizeState > 0) &
       allocate(damageState(ph)%subState0,source=damageState(ph)%state0)                 ! ToDo: hack
   enddo
 
@@ -574,8 +571,7 @@ subroutine phase_windForward(ip,el)
 
     call mechanical_windForward(ph,me)
 
-    if (phase_Nsources(ph) > 0) &
-      damageState(ph)%state0(:,me) = damageState(ph)%state(:,me)
+    if(damageState(ph)%sizeState > 0) damageState(ph)%state0(:,me) = damageState(ph)%state(:,me)
 
 
   enddo
diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index 2eda9540d..5e4e43017 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -18,6 +18,9 @@ submodule(phase) damagee
   integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:), allocatable :: &
     phase_source                                                                                !< active sources mechanisms of each phase
 
+  integer, dimension(:), allocatable :: &
+    phase_Nsources
+
   type(tDataContainer), dimension(:), allocatable :: current
 
   interface
@@ -156,6 +159,7 @@ module subroutine damage_init
     phase => phases%get(ph)
     sources => phase%get('damage',defaultVal=emptyList)
     if (sources%length > 1) error stop
+    phase_Nsources(ph) = sources%length
 
   enddo
 
@@ -192,7 +196,6 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
   integer :: &
     ph, &
     co, &
-    so, &
     me
 
    phiDot = 0.0_pReal
@@ -201,7 +204,7 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
    do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
      ph = material_phaseAt(co,el)
      me = material_phasememberAt(co,ip,el)
-     do so = 1, phase_Nsources(ph)
+
        select case(phase_source(ph))
          case (DAMAGE_ISOBRITTLE_ID)
            call isobrittle_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
@@ -222,7 +225,6 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
       end select
       phiDot = phiDot + localphiDot
       dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
-    enddo
   enddo
 
 end subroutine phase_damage_getRateAndItsTangents
@@ -258,7 +260,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
   ph = material_phaseAt(co,el)
   me = material_phaseMemberAt(co,ip,el)
 
-  if (phase_Nsources(ph) == 0) then
+  if (damageState(ph)%sizeState == 0) then
     broken = .false.
     return
   endif
@@ -377,7 +379,7 @@ function phase_damage_collectDotState(ph,me) result(broken)
 
   broken = .false.
 
-  if (phase_Nsources(ph)==1) then
+  if (damageState(ph)%sizeState > 0) then
 
     sourceType: select case (phase_source(ph))
 
@@ -420,7 +422,7 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
 
   broken = .false.
 
-  if (phase_Nsources(ph) == 0) return
+  if (damageState(ph)%sizeState == 0) return
 
      sourceType: select case (phase_source(ph))
 
@@ -461,7 +463,7 @@ function source_active(source_label)  result(active_source)
     phase => phases%get(ph)
     sources => phase%get('damage',defaultVal=emptyList)
     src => sources%get(1)
-    active_source(ph) = src%get_asString('type') == source_label
+    active_source(ph) = src%get_asString('type',defaultVal = 'x') == source_label
   enddo
 
 
diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index bc4b2cd47..ef962b7e6 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -3,6 +3,7 @@
 !----------------------------------------------------------------------------------------------------
 submodule(phase) mechanics
 
+
   enum, bind(c); enumerator :: &
     ELASTICITY_UNDEFINED_ID, &
     ELASTICITY_HOOKE_ID, &
@@ -22,8 +23,6 @@ submodule(phase) mechanics
     KINEMATICS_THERMAL_EXPANSION_ID
   end enum
 
-  integer(kind(KINEMATICS_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
-    phase_kinematics
   integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:),   allocatable :: &
     phase_elasticity                                                                                !< elasticity of each phase
   integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
@@ -1159,7 +1158,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
   subState0 = plasticState(ph)%State0(:,me)
 
-  if (phase_Nsources(ph) > 0) &
+  if (damageState(ph)%sizeState > 0) &
     damageState(ph)%subState0(:,me) = damageState(ph)%state0(:,me)
 
   subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
@@ -1187,7 +1186,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
         subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me)
         subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me)
         subState0 = plasticState(ph)%state(:,me)
-        if (phase_Nsources(ph) > 0) &
+        if (damageState(ph)%sizeState > 0) &
           damageState(ph)%subState0(:,me) = damageState(ph)%state(:,me)
 
       endif
@@ -1203,7 +1202,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
         phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0
       endif
       plasticState(ph)%state(:,me) = subState0
-      if (phase_Nsources(ph) > 0) &
+      if (damageState(ph)%sizeState > 0) &
         damageState(ph)%state(:,me) = damageState(ph)%subState0(:,me)
 
       todo = subStep > num%subStepMinCryst                          ! still on track or already done (beyond repair)
diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90
index 439f2bef9..d55ff8759 100644
--- a/src/phase_mechanical_eigen.f90
+++ b/src/phase_mechanical_eigen.f90
@@ -1,20 +1,26 @@
 submodule(phase:mechanics) eigendeformation
 
+  integer, dimension(:), allocatable :: &
+    Nmodels
+
+  integer(kind(KINEMATICS_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
+    model
+  integer(kind(KINEMATICS_UNDEFINED_ID)),  dimension(:), allocatable :: &
+    model_damage
+
   interface
-    module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics)
-      integer, intent(in) :: kinematics_length
-      logical, dimension(:,:), allocatable :: myKinematics
+    module function kinematics_cleavage_opening_init() result(myKinematics)
+      logical, dimension(:), allocatable :: myKinematics
     end function kinematics_cleavage_opening_init
 
-    module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics)
-      integer, intent(in) :: kinematics_length
-      logical, dimension(:,:), allocatable :: myKinematics
+    module function kinematics_slipplane_opening_init() result(myKinematics)
+      logical, dimension(:), allocatable :: myKinematics
     end function kinematics_slipplane_opening_init
 
-    module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics)
+    module function thermalexpansion_init(kinematics_length) result(myKinematics)
       integer, intent(in) :: kinematics_length
       logical, dimension(:,:), allocatable :: myKinematics
-    end function kinematics_thermal_expansion_init
+    end function thermalexpansion_init
 
     module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
       integer, intent(in) :: ph, me
@@ -65,28 +71,27 @@ module subroutine eigendeformation_init(phases)
   print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation init  -+>>>'
 
 !--------------------------------------------------------------------------------------------------
-! initialize kinematic mechanisms
-  allocate(phase_Nkinematics(phases%length),source = 0)
+! explicit eigen mechanisms
+  allocate(Nmodels(phases%length),source = 0)
+
   do ph = 1,phases%length
     phase => phases%get(ph)
     kinematics => phase%get('kinematics',defaultVal=emptyList)
-    phase_Nkinematics(ph) = kinematics%length
-    kinematics => phase%get('damage',defaultVal=emptyList)
-    if(kinematics%length >0) then
-      damage => kinematics%get(1)
-      if(damage%get_asString('type',defaultVal='n/a') == 'anisobrittle')  phase_Nkinematics(ph) =  phase_Nkinematics(ph) +1
-      if(damage%get_asString('type',defaultVal='n/a') == 'isoductile'  )  phase_Nkinematics(ph) =  phase_Nkinematics(ph) +1
-    endif
+    Nmodels(ph) = kinematics%length
   enddo
 
-  allocate(phase_kinematics(maxval(phase_Nkinematics),phases%length), source = KINEMATICS_undefined_ID)
+  allocate(model(maxval(Nmodels),phases%length), source = KINEMATICS_undefined_ID)
 
-  if(maxval(phase_Nkinematics) /= 0) then
-    where(kinematics_cleavage_opening_init(maxval(phase_Nkinematics)))  phase_kinematics = KINEMATICS_cleavage_opening_ID
-    where(kinematics_slipplane_opening_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_slipplane_opening_ID
-    where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) phase_kinematics = KINEMATICS_thermal_expansion_ID
+  if(maxval(Nmodels) /= 0) then
+    where(thermalexpansion_init(maxval(Nmodels))) model = KINEMATICS_thermal_expansion_ID
   endif
 
+  allocate(model_damage(phases%length),  source = KINEMATICS_UNDEFINED_ID)
+
+  where(kinematics_cleavage_opening_init())  model_damage = KINEMATICS_cleavage_opening_ID
+  where(kinematics_slipplane_opening_init()) model_damage = KINEMATICS_slipplane_opening_ID
+
+
 end subroutine eigendeformation_init
 
 
@@ -125,11 +130,10 @@ end function kinematics_active
 !--------------------------------------------------------------------------------------------------
 !> @brief checks if a kinematic mechanism is active or not
 !--------------------------------------------------------------------------------------------------
-function kinematics_active2(kinematics_label,kinematics_length)  result(active_kinematics)
+function kinematics_active2(kinematics_label)  result(active_kinematics)
 
-  character(len=*), intent(in)         :: kinematics_label                                          !< name of kinematic mechanism
-  integer,          intent(in)         :: kinematics_length                                         !< max. number of kinematics in system
-  logical, dimension(:,:), allocatable :: active_kinematics
+  character(len=*), intent(in)       :: kinematics_label                                            !< name of kinematic mechanism
+  logical, dimension(:), allocatable :: active_kinematics
 
   class(tNode), pointer :: &
     phases, &
@@ -139,13 +143,14 @@ function kinematics_active2(kinematics_label,kinematics_length)  result(active_k
   integer :: p
 
   phases => config_material%get('phase')
-  allocate(active_kinematics(kinematics_length,phases%length), source = .false. )
+  allocate(active_kinematics(phases%length), source = .false. )
   do p = 1, phases%length
     phase => phases%get(p)
     kinematics => phase%get('damage',defaultVal=emptyList)
+    if(kinematics%length < 1) return
     kinematics_type => kinematics%get(1)
     if (.not. kinematics_type%contains('type')) continue
-    active_kinematics(1,p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label
+    active_kinematics(p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label
   enddo
 
 
@@ -181,7 +186,9 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
     detFi
   integer :: &
     k, i, j
+  logical :: active
 
+  active = .false.
   Li = 0.0_pReal
   dLi_dS  = 0.0_pReal
   dLi_dFi = 0.0_pReal
@@ -190,30 +197,37 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
   plasticType: select case (phase_plasticity(ph))
     case (PLASTICITY_isotropic_ID) plasticType
       call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,phase_plasticInstance(ph),me)
-    case default plasticType
-      my_Li = 0.0_pReal
-      my_dLi_dS = 0.0_pReal
+      Li = Li + my_Li
+      dLi_dS = dLi_dS + my_dLi_dS
+      active = .true.
   end select plasticType
 
-  Li = Li + my_Li
-  dLi_dS = dLi_dS + my_dLi_dS
 
-  KinematicsLoop: do k = 1, phase_Nkinematics(ph)
-    kinematicsType: select case (phase_kinematics(k,ph))
-      case (KINEMATICS_cleavage_opening_ID) kinematicsType
-        call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me)
-      case (KINEMATICS_slipplane_opening_ID) kinematicsType
-        call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me)
+  KinematicsLoop: do k = 1, Nmodels(ph)
+    kinematicsType: select case (model(k,ph))
       case (KINEMATICS_thermal_expansion_ID) kinematicsType
         call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,me)
-      case default kinematicsType
-        my_Li = 0.0_pReal
-        my_dLi_dS = 0.0_pReal
+        Li = Li + my_Li
+        dLi_dS = dLi_dS + my_dLi_dS
+        active = .true.
     end select kinematicsType
-    Li = Li + my_Li
-    dLi_dS = dLi_dS + my_dLi_dS
   enddo KinematicsLoop
 
+  select case (model_damage(ph))
+    case (KINEMATICS_cleavage_opening_ID)
+      call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me)
+      Li = Li + my_Li
+      dLi_dS = dLi_dS + my_dLi_dS
+      active = .true.
+    case (KINEMATICS_slipplane_opening_ID)
+      call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me)
+      Li = Li + my_Li
+      dLi_dS = dLi_dS + my_dLi_dS
+      active = .true.
+  end select
+
+  if(.not. active) return
+
   FiInv = math_inv33(Fi)
   detFi = math_det33(Fi)
   Li = matmul(matmul(Fi,Li),FiInv)*detFi                                                            !< push forward to intermediate configuration
diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90
index 5bb9847f4..aaa1beb97 100644
--- a/src/phase_mechanical_eigen_cleavageopening.f90
+++ b/src/phase_mechanical_eigen_cleavageopening.f90
@@ -28,12 +28,11 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function kinematics_cleavage_opening_init(kinematics_length) result(myKinematics)
+module function kinematics_cleavage_opening_init() result(myKinematics)
 
-  integer, intent(in)                  :: kinematics_length
-  logical, dimension(:,:), allocatable :: myKinematics
+  logical, dimension(:), allocatable :: myKinematics
 
-  integer :: Ninstances,p,k
+  integer :: p
   integer, dimension(:), allocatable :: N_cl                                                        !< active number of cleavage systems per family
   character(len=pStringLen) :: extmsg = ''
   class(tNode), pointer :: &
@@ -42,24 +41,24 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
     kinematics, &
     kinematic_type
 
-  print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation:cleavageopening init  -+>>>'
 
-  myKinematics = kinematics_active2('anisobrittle',kinematics_length)
-  Ninstances = count(myKinematics)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
+  myKinematics = kinematics_active2('anisobrittle')
+  if(count(myKinematics) == 0) return
+
+  print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation:cleavageopening init  -+>>>'
+  print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
+
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
 
   do p = 1, phases%length
-    phase => phases%get(p)
-    if(count(myKinematics(:,p)) == 0) cycle
-    kinematics => phase%get('damage')
-    do k = 1, kinematics%length
-      if(myKinematics(k,p)) then
-        associate(prm  => param(p))
-        kinematic_type => kinematics%get(k)
+    if(myKinematics(p)) then
+      phase => phases%get(p)
+      kinematics => phase%get('damage')
+
+      associate(prm  => param(p))
+        kinematic_type => kinematics%get(1)
 
         N_cl = kinematic_type%get_asInts('N_cl')
         prm%sum_N_cl = sum(abs(N_cl))
@@ -83,9 +82,8 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
         if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(cleavage_opening)')
-        end associate
-      endif
-    enddo
+      end associate
+    endif
   enddo
 
 
diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90
index 64988d87f..ccba98b14 100644
--- a/src/phase_mechanical_eigen_slipplaneopening.f90
+++ b/src/phase_mechanical_eigen_slipplaneopening.f90
@@ -32,12 +32,11 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function kinematics_slipplane_opening_init(kinematics_length) result(myKinematics)
+module function kinematics_slipplane_opening_init() result(myKinematics)
 
-  integer, intent(in)                  :: kinematics_length
-  logical, dimension(:,:), allocatable :: myKinematics
+  logical, dimension(:), allocatable :: myKinematics
 
-  integer :: Ninstances,p,i,k
+  integer :: p,i
   character(len=pStringLen) :: extmsg = ''
   integer,     dimension(:),   allocatable :: N_sl
   real(pReal), dimension(:,:), allocatable :: d,n,t
@@ -49,26 +48,26 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
     kinematics, &
     kinematic_type
 
-  print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation:slipplaneopening init  -+>>>'
 
-  myKinematics = kinematics_active2('isoductile',kinematics_length)
-  Ninstances = count(myKinematics)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
+  myKinematics = kinematics_active2('isoductile')
+  if(count(myKinematics) == 0) return
+  print'(/,a)', ' <<<+-  phase:mechanics:eigendeformation:slipplaneopening init  -+>>>'
+  print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
+
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
 
   do p = 1, phases%length
-    phase => phases%get(p)
-    mech  => phase%get('mechanics')
-    pl    => mech%get('plasticity')
-    if(count(myKinematics(:,p)) == 0) cycle
-    kinematics => phase%get('damage')
-    do k = 1, kinematics%length
-      if(myKinematics(k,p)) then
-        associate(prm  => param(p))
-        kinematic_type => kinematics%get(k)
+    if(myKinematics(p)) then
+      phase => phases%get(p)
+      mech  => phase%get('mechanics')
+      pl    => mech%get('plasticity')
+
+      kinematics => phase%get('damage')
+
+      associate(prm  => param(p))
+        kinematic_type => kinematics%get(1)
 
         prm%dot_o    = kinematic_type%get_asFloat('dot_o')
         prm%q        = kinematic_type%get_asFloat('q')
@@ -103,9 +102,8 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
 !  exit if any parameter is out of range
         if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(slipplane_opening)')
 
-        end associate
-      endif
-    enddo
+      end associate
+    endif
   enddo
 
 
diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90
index 6630f0eb7..1041e0957 100644
--- a/src/phase_mechanical_eigen_thermalexpansion.f90
+++ b/src/phase_mechanical_eigen_thermalexpansion.f90
@@ -23,7 +23,7 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function kinematics_thermal_expansion_init(kinematics_length) result(myKinematics)
+module function thermalexpansion_init(kinematics_length) result(myKinematics)
 
   integer, intent(in)                  :: kinematics_length
   logical, dimension(:,:), allocatable :: myKinematics
@@ -77,7 +77,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
   enddo
 
 
-end function kinematics_thermal_expansion_init
+end function thermalexpansion_init
 
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index 5d3d89ccb..012554aef 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -3,6 +3,9 @@
 !----------------------------------------------------------------------------------------------------
 submodule(phase) thermal
 
+  integer, dimension(:), allocatable :: &
+    thermal_Nsources
+
   type(tSourceState),  allocatable, dimension(:) :: &
     thermalState
 
@@ -36,8 +39,6 @@ submodule(phase) thermal
     end function externalheat_init
 
 
-
-
     module subroutine externalheat_dotState(ph, me)
       integer, intent(in) :: &
         ph, &
diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90
index 3c0095401..9c75763dc 100644
--- a/src/phase_thermal_dissipation.f90
+++ b/src/phase_thermal_dissipation.f90
@@ -31,15 +31,14 @@ module function dissipation_init(source_length) result(mySources)
     phase, &
     sources, thermal, &
     src
-  integer :: Ninstances,so,Nconstituents,ph
+  integer :: so,Nconstituents,ph
 
-  print'(/,a)', ' <<<+-  phase:thermal:dissipation init  -+>>>'
 
   mySources = thermal_active('dissipation',source_length)
+  if(count(mySources) == 0) return
+  print'(/,a)', ' <<<+-  phase:thermal:dissipation init  -+>>>'
+  print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
-  Ninstances = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90
index d2874ded0..c7f894215 100644
--- a/src/phase_thermal_externalheat.f90
+++ b/src/phase_thermal_externalheat.f90
@@ -8,7 +8,7 @@ submodule(phase:thermal) externalheat
 
 
   integer,           dimension(:),   allocatable :: &
-    source_thermal_externalheat_offset                                                           !< which source is my current thermal dissipation mechanism?
+    source_thermal_externalheat_offset                                                              !< which source is my current thermal dissipation mechanism?
 
   type :: tParameters                                                                               !< container type for internal constitutive parameters
     real(pReal), dimension(:), allocatable :: &
@@ -38,15 +38,14 @@ module function externalheat_init(source_length) result(mySources)
     phase, &
     sources, thermal, &
     src
-  integer :: Ninstances,so,Nconstituents,ph
+  integer :: so,Nconstituents,ph
 
-  print'(/,a)', ' <<<+-  phase:thermal:externalheat init  -+>>>'
 
   mySources = thermal_active('externalheat',source_length)
+  if(count(mySources) == 0) return
+  print'(/,a)', ' <<<+-  phase:thermal:externalheat init  -+>>>'
+  print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
-  Ninstances = count(mySources)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-  if(Ninstances == 0) return
 
   phases => config_material%get('phase')
   allocate(param(phases%length))