ntwin = [0] cause problems

math_expand now works for sum(how) == 0

src/math.f90

@@ -386,6 +386,9 @@ pure function math_expand(what,how)
  real(pReal), dimension(sum(how)) ::  math_expand
  integer(pInt) :: i
 
+ if (sum(how) == 0_pInt) &
+   return
+
  do i = 1_pInt, size(how)
    math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1_pInt,size(what))+1_pInt)
  enddo

src/plastic_phenopowerlaw.f90

@@ -202,7 +202,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
 
  allocate(totalNslip(maxNinstance),                    source=0_pInt)
  allocate(totalNtwin(maxNinstance),                    source=0_pInt)
-  allocate(param(maxNinstance))                                                                      ! one container of parameters per instance
+ allocate(param(maxNinstance))                                                                      ! one container of parameters per instance
 
  rewind(fileUnit)
  phase = 0_pInt