anchor.diff:        endif
   
   ! Eshelby's Sijkl tensor (appendix of Tandon and Weng 1984).
  -     s11 = (t3 + (t4-1.0)/t1 - (t3 + t4/t1)*g)/(2.0*t2)              
  -     s22 = (t4/(t1*2.0) + (t3 - 9.0/(4.0*t1))*g)/(4.0*t2)            
  -     s33 = s22                                                       
  -     s23 = (del_in**2/(2.0*t1) - (t3 + 3.0/(4.0*t1))*g)/(4.0*t2)       
  -     s32 = s23                                                       
  -     s21 = (-2.0*del_in*del_in/t1 + (t4/t1 - t3)*g)/(4.0*t2)               
  -     s31 = s21                                                       
  -     s12 = (-1.0*(t3 + 1.0/t1) + (t3 + 3.0/(2.0*t1))*g)/(2.0*t2)     
  -     s13 = s12                                                       
  -     s44 = (del_in*del_in/(2.0*t1) + (t3 - 3.0/(4.0*t1))*g)/(4.0*t2)       
  +     s11 = (t3 + (t4-1.0)/t1 - (t3 + t4/t1)*g)/(2.0*t2)
  +     s22 = (t4/(t1*2.0) + (t3 - 9.0/(4.0*t1))*g)/(4.0*t2)
  +     s33 = s22
  +     s23 = (del_in**2/(2.0*t1) - (t3 + 3.0/(4.0*t1))*g)/(4.0*t2)
  +     s32 = s23
  +     s21 = (-2.0*del_in*del_in/t1 + (t4/t1 - t3)*g)/(4.0*t2)
  +     s31 = s21
  +     s12 = (-1.0*(t3 + 1.0/t1) + (t3 + 3.0/(2.0*t1))*g)/(2.0*t2)
  +     s13 = s12
  +     s44 = (del_in*del_in/(2.0*t1) + (t3 - 3.0/(4.0*t1))*g)/(4.0*t2)
        s66 = (t3 - (t1+2.0)/t1 - (t3 - 3.0*(t1+2.0)/t1)*g/2.0)/(4.0*t2)
  -     s55 = s66                                                       
  -     
  +     s55 = s66
  +
   ! Tandon and Weng's B terms (after equation 17).
  -     B1 = c_in*D1 + D2 + (1.0-c_in)*(D1*s11 + 2.0*s21)    
  -     B2 = c_in + D3 + (1.0-c_in)*(D1*s12 + s22 + s23)     
  -     B3 = c_in + D3 + (1.0-c_in)*(s11 + (1.0+D1)*s21)     
  -     B4 = c_in*D1 + D2 + (1.0-c_in)*(s12 + D1*s22 + s23)  
  -     B5 = c_in + D3 + (1.0-c_in)*(s12 + s22 + D1*s23)     
  -     
  +     B1 = c_in*D1 + D2 + (1.0-c_in)*(D1*s11 + 2.0*s21)
  +     B2 = c_in + D3 + (1.0-c_in)*(D1*s12 + s22 + s23)
  +     B3 = c_in + D3 + (1.0-c_in)*(s11 + (1.0+D1)*s21)
  +     B4 = c_in*D1 + D2 + (1.0-c_in)*(s12 + D1*s22 + s23)
  +     B5 = c_in + D3 + (1.0-c_in)*(s12 + s22 + D1*s23)
  +
   ! Tandon and Weng's A terms (after equation 20).
  -     A1 = D1*(B4 + B5) - 2.0*B2    
  +     A1 = D1*(B4 + B5) - 2.0*B2
        A2 = (1.0 + D1)*B2 - (B4 + B5)
  -     A3 = B1 - D1*B3               
  -     A4 = (1.0 + D1)*B1 - 2.0*B3   
  -     A5 = (1.0 - D1)/(B4 - B5)     
  -     A = 2.0*B2*B3 - B1*(B4+B5)    
  -     
  -! Tandon and Weng (1984) equations (25) (28) (31) (32) 
  +     A3 = B1 - D1*B3
  +     A4 = (1.0 + D1)*B1 - 2.0*B3
  +     A5 = (1.0 - D1)/(B4 - B5)
  +     A = 2.0*B2*B3 - B1*(B4+B5)
  +
  +! Tandon and Weng (1984) equations (25) (28) (31) (32)
        E11 = E0 /(1.0+c_in*(A1+2.0*anu*A2)/A)
        E22 = E0 &
            /(1.0+c_in*(-2.0*anu*A3 + (1.0-anu)*A4 + (1.0+anu)*A5*A)/(2.0*A))
        amu12 = amu*(1.0 + c_in/(amu/(amui-amu) + 2.0*(1.0-c_in)*s66))
        amu23 = amu*(1.0 + c_in/(amu/(amui-amu) + 2.0*(1.0-c_in)*s44))
  -     
  +
   ! Sayers equation (36)
        anu31 = anu - c_in*(anu*(A1+2.0*anu*A2)+(A3-anu*A4)) &
                   /(A + c_in*(A1+2.0*anu*A2))
  -     
  +
   ! T&W equation (36)
   !     aK12 term; bmps=plane strain bulk modulus
        anum = (1.0+anu)*(1.0-2.0*anu)
  
positive.code:                          'Require TI with hexad // x3.'
           stop
        endif
        
        eps = (C(1,1) - C(3,3))/(2._rs*C(3,3))
        gam = (C(6,6) - C(4,4))/(2._rs*C(4,4))
        del = ((C(1,3)+C(4,4))**2-(C(3,3)-C(4,4))**2)/(2._rs*C(3,3)*(C(3,3)-C(4,4)))
     
     end subroutine CIJ_VTI2thom
  !-------------------------------------------------------------------------------
  
  
negative.code:        end do ! j
        end do ! i
        ABS_PIVOT = abs(PIVOT)
  
        ! HAVE PIVOT, INTERCHANGE ROW, COL POINTERS
        HOLD = ROW(k)
        ROW(k) = ROW(I_PIVOT)
        ROW(I_PIVOT) = HOLD
        HOLD = COL(k)
        COL(k) = COL(J_PIVOT)
        COL(J_PIVOT) = HOLD
  
        ! CHECK FOR NEAR SINGULAR
        if( ABS_PIVOT < 1.0D-52*NORM1 ) then
         do j=1,n
           AI(ROW(k),j) = 0.0D0
  