anchor.diff:        real(rs),intent(in) :: VF_in(:), C_in(:,:,:), rh_in(:)
         real(rs),intent(out) :: Cave(6,6), rhave
         integer :: i,n
  -      real(rs),allocatable :: VF(:)
  +!      real(rs),allocatable :: VF(:)
         real(rs) :: S(6,6), S_in(6,6), C_temp(6,6)
         
   !  Get size of arrays and check they're consistent
         n = size(VF_in)
         if (size(C_in,1) /= n .or. size(rh_in) /= n) then
  -         write(0,'(a)') 'anisotropy_ajn: CIJ_Voigt_av: input VF, C and rh must be same length.'
  +         write(0,'(a)') 'anisotropy_ajn: CIJ_Reuss_av: input VF, C and rh must be same length.'
            stop
         elseif (size(C_in,2) /= 6 .or. size(C_in,3) /= 6) then
  -         write(0,'(a)') 'anisotropy_ajn: CIJ_Voigt_av: C must be nx6x6 array.'
  +         write(0,'(a)') 'anisotropy_ajn: CIJ_Reuss_av: C must be nx6x6 array.'
            stop
         endif
         
   !  Allocate space for arrays
  -      allocate(VF(n))
  +!      allocate(VF(n))
         
   !  Normalise the volume fractions to sum to unity
  -      VF = VF_in / sum(VF_in)
  +!      VF = VF_in / sum(VF_in)
         
   !  Construct Reuss average
         Cave = 0.  ;  S = 0.;  rhave = 0.
  
positive.code:          !  Find compliance from input stiffness
           C_temp = C_in(i,:,:)
           call inverse(6,6,C_temp,S_in)
           S = S + VF(i)*S_in
           rhave = rhave + VF(i)*rh_in(i)
        enddo
        
  !  Find stiffness from compliance
        call inverse(6,6,S,Cave)
        
        deallocate(VF)
        
     end subroutine CIJ_Reuss_av
  !-------------------------------------------------------------------------------
  
negative.code:       endif
        
  !  Allocate space for arrays
  !      allocate(VF(n))
        
  !  Normalise the volume fractions to sum to unity
  !      VF = VF_in / sum(VF_in)
        
  !  Construct Voigt average
        Cave = 0.  ;  rhave = 0.
        do i=1,n
           Cave = Cave + VF_in(i)*C_in(i,:,:)/sum(VF_in)
           rhave = rhave + VF_in(i)*rh_in(i)/sum(VF_in)
        enddo
        
  !      deallocate(VF)
        
     end subroutine CIJ_Voigt_av
  !-------------------------------------------------------------------------------
  