C
C User subroutine vumullins
      subroutine vumullins ( 
C Read only (unmodifiable) variables -
     1     nblock,
     2     jElem, kIntPt, kLayer, kSecPt,
     3     cmname,
     4     nstatev, nfieldv, nprops,
     5     props, tempOld, tempNew, fieldOld, fieldNew,
     6     stateOld, enerDamageOld,
     7     uMaxOld, uMaxNew, uDev, 
C Write only (modifiable) variables -
     8     eta, detaDuDev,
     9     stateNew, enerDamageNew )
C
      include 'vaba_param.inc' 
C
      dimension props(nprops),
     1     tempOld(nblock),
     2     fieldOld(nblock,nfieldv),
     3     stateOld(nblock,nstatev),
     4     tempNew(nblock),
     5     fieldNew(nblock,nfieldv),
     6     enerDamageOld(nblock),
     7     uMaxOld(nblock), uMaxNew(nblock),
     8     uDev(nblock),
     9     eta(nblock), detaDuDev(nblock),
     1     stateNew(nblock,nstatev),
     2     enerDamageNew(nblock)
C
      character*80 cmname 
C
C Parameters for the evaluation of the Ogden-Roxburgh function
C for *MULLINS EFFECT
C
      parameter ( r1 = -1.26551223d0, r2 = 1.00002368d0,
     *     r3 = 0.37409196d0, r4 = 0.09678418d0,
     *     r5 =-0.18628806d0, r6 = 0.27886807d0,
     *     r7 =-1.13520398d0, r8 = 1.48851587d0,
     *     r9 =-0.82215223d0, r10= 0.17087277d0,
     *     sqrtPi = 1.772453851d0, twoSqrtPiInv = 2.d0/sqrtPi,
     *     rlogTwo=0.693147180559945d0 )
C
      parameter ( zero = 0.d0, one = 1.d0, half = 0.5d0 )
C
      parameter ( rMinExp = -126*rlogTwo )
      parameter ( small = 2**(-126) )
C
      R  = props(1)
      xM = props(2)
      xB = props(3)
      Rinv = one / R
*
*  Evaluation of the Ogden-Roxburgh eta function and derivatives
*  for *MULLINS EFFECT
*
      do k=1, nblock
*
         if ( uDev(k) .lt. uMaxOld(k) ) then
            uMax = uMaxOld(k)
            tmp =  xM + xB * uMax
            test = half + sign( half, small - abs( tmp ) )
            xMinv = ( one - test ) / ( test + tmp )
            arg = xMinv * ( uMax - uDev(k) )
            arg2 = arg * arg
*     -- Evaluate the error function (arg .ge. zero)
            z = abs(arg)
            t = one/(one+half*z)
            arg1 = -z*z + r1 + t * (r2 + t * (r3 + t * (r4
     *           + t * (r5 + t * (r6 + t * (r7 + t * (r8
     *           + t * (r9 + t * r10))))))))
            arg1 = max ( arg1, rMinExp )
            ans = t*exp(arg1)
            erfc = ans
            erf = one - erfc
            etaT = one - Rinv * erf
         else
* --	Assume unloading modulus on primary curve
            etaT = one
            uMax = uDev(k)
            arg2 = zero
            tmp =  xM + xB * uMax
            test = half + sign( half,  small - abs( tmp ) )
            xMinv = ( one - test ) / ( test + tmp )
         end if
* -- Derivative of eta
         eta(k) = etaT
         detaDuDev(k) = twoSqrtPiInv*Rinv*xMinv*exp(max(-arg2,rMinExp))
      end do
*
      do k = 1, nblock
         denom = xM + xB * uMaxNew(k)
         test = half + sign( half, small - abs( denom ) )
         xMinv = ( one - test ) / ( test + denom ) 
         arg = xMinv * uMaxNew(k)
         arg2 = arg * arg
         arg2 = - max ( -arg2, rMinExp )
         z = abs(arg)
         t = one/(one+half*z)
         arg1 = -z*z + r1 + t * (r2 + t * (r3 + t * (r4
     *        + t * (r5 + t * (r6 + t * (r7 + t * (r8
     *        + t * (r9 + t * r10))))))))
         arg1 = max ( arg1, rMinExp )
         ans = t*exp(arg1)
         erfc = ans
         erf = one - erfc
         factor = Rinv * ( erf * uMaxNew(k)
     *        - denom * ( one - exp(- arg2) ) / sqrtPi  )
         enerDamageNew(k) = factor 
      end do
*
      return
      end