anchor.diff:  
         if(mtype==VISCOELASTIC+2) then   ! Arrhenius
   
  -        hsn1 = mvar(2)*( 1.d0/(tn1-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
  -        hsn = mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
  -        asn1 = dexp(-hsn1)
  -        asn = dexp(-hsn)
  +        if( tn1==tn ) then
  +          hsn = mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
  +          trsinc = dexp(hsn)
  +        else
  +          hsn1 = mvar(2)*( 1.d0/(tn1-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
  +          hsn = mvar(2)*( 1.d0/(tn-mvar(3))-1.d0/(mvar(1)-mvar(3)) )/mvar(4)
  +          asn1 = dexp(-hsn1)
  +          asn = dexp(-hsn)
  +          trsinc = (asn1-asn)/(hsn1-hsn)
  +        endif
   
         else   ! WLF
   
  -        hsn1 = mvar(2)*( tn1-mvar(1) )/ (mvar(3)+tn1-mvar(1) )*dlog(10.d0)
  -        hsn = mvar(2)*( tn-mvar(1) )/ (mvar(3)+tn-mvar(1) )*dlog(10.d0)
  -        asn1 = dexp(hsn1)
  -        asn = dexp(hsn)
  -
  -      endif
  -
  -      if(dabs(hsn1-hsn) < 1.d-04) then
  -        trsinc = dexp( (hsn+hsn1)/2.d0 )
  -      else
  -        trsinc = (asn1-asn)/(hsn1-hsn)
  +        if( mvar(3)+tn1-mvar(1)<=0.d0 ) then
  +           trsinc= -1.d0; return
  +        endif
  +        if( mvar(3)+tn-mvar(1)<=0.d0 ) then
  +           trsinc= -1.d0; return
  +        endif
  +        if( tn1==tn ) then
  +          hsn = mvar(2)*( tn-mvar(1) )/ (mvar(3)+tn-mvar(1) )*dlog(10.d0)
  +          trsinc = dexp(hsn)
  +        else
  +          hsn1 = mvar(2)*( tn1-mvar(1) )/ (mvar(3)+tn1-mvar(1) )*dlog(10.d0)
  +          hsn = mvar(2)*( tn-mvar(1) )/ (mvar(3)+tn-mvar(1) )*dlog(10.d0)
  +          asn1 = dexp(hsn1)
  +          asn = dexp(hsn)
  +          trsinc = (asn1-asn)/(hsn1-hsn)
  +        endif
  +
         endif
   
     end function
  
positive.code:       real(kind=kreal) :: G,Gg,K,Kg,Kth, exp_n,mu_0,mu_n,dq_n,dtau, theta
        real(kind=kreal) :: ina(1), outa(2), EE, PP, ddt
  
        real(kind=kreal) :: devstrain(6), en(6)
        type(tTable), pointer  :: dicval
        logical :: ierr
  
  
negative.code:           PP = outa(2)
          endif
          if( matl%mtype>VISCOELASTIC ) ddt=trsinc(temp,tempn, matl%mtype, matl%variables)*dt
          if( ddt<=0.d0 ) then
             call calElasticMatrix( matl, sectType, D, temp  )
             sig = matmul( D(:,:), eps(:) )
             vsig = 0.d0
             return
          endif
        else
          EE = matl%variables(M_YOUNGS)
          PP = matl%variables(M_POISSON)
  