anchor.diff:      !---------------------------------------------------------------------
   
       real(kind=kreal) :: ALP, alp0, D(6, 6), B(6, ndof*nn)
  -    real(kind=kreal) :: B4, B6, B8, det, ecoord(3, nn)
  +    real(kind=kreal) :: Z1(3), det, ecoord(3, nn)
       integer(kind=kint) :: j, LX, serr
       real(kind=kreal) :: estrain(6), SGM(6), H(nn)
       real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
  -    real(kind=kreal) :: wg, outa(1), ina(1), Bbar(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
  -    real(kind=kreal) :: TEMPC, TEMP0, TEMPC0, TEMP00, THERMAL_EPS, tm(6,6)
  +    real(kind=kreal) :: wg, outa(1), ina(1), gderiv1_ave(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
  +    real(kind=kreal) :: TEMPC, TEMP0, V0, jacob_ave, THERMAL_EPS, tm(6,6)
       logical :: ierr, matlaniso
   
       !---------------------------------------------------------------------
  
positive.code:         dstrain(4) = dot_product(Fbar(1:3,1),Fbar(1:3,2))
          dstrain(5) = dot_product(Fbar(1:3,2),Fbar(1:3,3))
          dstrain(6) = dot_product(Fbar(1:3,3),Fbar(1:3,1))
  
          if( mtype == NEOHOOKE .OR. mtype == MOONEYRIVLIN .OR.  mtype == ARRUDABOYCE  .OR.         &
              mtype == USERELASTIC .OR. mtype == USERHYPERELASTIC .OR. mtype == USERMATERIAL ) then
  
negative.code:       !**
        if( matlaniso ) then
          do j = 1, 3
            estrain(j) = ALPO(j)*(TEMPC-ref_temp)-alpo0(j)*(TEMP0-ref_temp)
          end do
          estrain(4:6) = 0.0D0
          call transformation(coordsys, tm)
          estrain(:) = matmul( estrain(:), tm  )      ! to global coord
          estrain(4:6) = estrain(4:6)*2.0D0
        else
          THERMAL_EPS = ALP*(TEMPC-ref_temp)-alp0*(TEMP0-ref_temp)
          estrain(1:3) = THERMAL_EPS
          estrain(4:6) = 0.0D0
        end if
  