C
C  User subroutine to define damage-parameter eta for Mullins effect

      subroutine umullins(numProps,props,uMaxNew,uMaxOld,sedDev,
     *     eta,detadW,dmgDissOld,dmgDissNew,senerNew,numStateV,stateV,
     *     temp,dtemp,numFieldV,fieldV,fieldVInc,cmname,linper)
C
      include 'aba_param.inc'
C
      character*80 cmname
      dimension props(*),stateV(*),fieldV(*),fieldVInc(*)
C
      parameter ( c1 = -1.26551223d0, c2 = 1.00002368d0,
     *     c3 = 0.37409196d0, c4 = 0.09678418d0,
     *     c5 =-0.18628806d0, c6 = 0.27886807d0,
     *     c7 =-1.13520398d0, c8 = 1.48851587d0,
     *     c9 =-0.82215223d0, c10= 0.17087277d0,
     *     one=1.d0,two=2.d0,half=0.5d0,pi=3.14159265358979d0,
     *     zero=0.d0,tol=1.d-15)   
C
      if (cmname(1:6) .eq. 'YEOH1') then
         continue
      else if ((cmname(1:6) .eq. 'YEOH2').and.numStateV.gt.0) then
         uMaxNew = max(uMaxNew,stateV(1))
      end if
C
      rMullins=props(1)
      amMullins=props(2)
      betaMullins=props(3)
      rInv = one/rMullins
C
      oSqrtPi=one/sqrt(pi)
      term=(amMullins+betaMullins*uMaxNew)
      if (uMaxNew.lt.zero) then
         arg1=zero
      else
         arg1=(uMaxNew-sedDev)/term
      end if
C
      if (arg1.gt.tol) then
         z = abs(arg1)
         t = one/(one+half*z)
         erfc = t*exp(-z*z + c1 + t * (c2 + t * (c3 + t * (c4
     *        + t * (c5 + t * (c6 + t * (c7 + t * (c8
     *        + t * (c9 + t * c10)))))))))
         erf = one - erfc
         eta = one - rInv*erf
         if (eta.gt.one) eta=one
      else
         eta = one
      end if
C
      detadW = oSqrtPi*rInv*two*exp(-arg1*arg1)/term
C
C     Damage dissipation calculations
      if (uMaxNew.gt.zero) then
         arg2=uMaxNew/term
      else
         arg2=zero
      end if
      z = abs(arg2)
      t = one/(one+half*z)
      erfc = t*exp(-z*z + c1 + t * (c2 + t * (c3 + t * (c4
     *           + t * (c5 + t * (c6 + t * (c7 + t * (c8
     *           + t * (c9 + t * c10)))))))))
      erf = one - erfc
      term1=exp(-arg2*arg2)
      dmgDissNew=erf*uMaxNew
     *     - oSqrtPi*term*(one-term1)
      dmgDissNew = rInv*dmgDissNew     
C
C    Compute phi(eta)
      term2=exp(-arg1*arg1)
      phiOfEta=oSqrtPi*rInv*term*(term2-one) + (one-eta)*uMaxNew      
C
      senerNew=eta*sedDev+phiOfEta-dmgDissNew
C
      return
      end