anchor.diff:      real(kind=kreal)   :: D(6,6), B(6,ndof*nn), B1(6,ndof*nn), spfunc(nn), ina(1)
       real(kind=kreal)   :: gderiv(nn,3), gdispderiv(3,3), det, WG, ttc,tt0, ttn,outa(1)
       integer(kind=kint) :: i, j, k, LX, mtype, serr
  +    integer(kind=kint) :: isEp
       real(kind=kreal)   :: naturalCoord(3), rot(3,3), mat(6,6), EPSTH(6)
       real(kind=kreal)   :: totaldisp(3,nn), elem(3,nn), elem1(3,nn), coordsys(3,3), tm(6,6)
       real(kind=kreal)   :: dstrain(6),dstress(6),dumstress(3,3),dum(3,3)
  
positive.code:         gausses(LX)%strain(1:6) = gausses(LX)%strain_bak(1:6)+ dstrain(1:6)+EPSTH(:)
  
          IF( isViscoelastic(mtype) .AND. tincr /= 0.0D0 ) THEN
            gausses(LX)%pMaterial%mtype = mtype
            IF( PRESENT(TT) .AND. PRESENT(T0) ) THEN
              CALL StressUpdate( gausses(LX), D3, dstrain, gausses(LX)%stress, tincr, ttc, tt0 )
            ELSE
  
negative.code:         ttc = DOT_PRODUCT(TT, spfunc)
          tt0 = DOT_PRODUCT(T0, spfunc)
          ttn = DOT_PRODUCT(TN, spfunc)
          CALL MatlMatrix( gausses(LX), D3, D, tincr, coordsys, ttc, isEp )
  
          ina(1) = ttc
          IF( matlaniso ) THEN
  