      subroutine umshl53(cm,ldcmw,capa,dt1,time,
     $     deps,ldeps,
     $     sig,ldsig,
     $     hisv,ldhv,
     $     fail,
     $     tsh,
     $     ndsh)
#ifdef ACC
#ifdef DP
      implicit double precision (a-h,o-z)                                    dp
#endif
c
c     Plane stress version of the brittle damage model
c     for belytschko-tsay shell element.
c        (CRACK GROWTH MODEL)
c     Haeng-Ki Lee &
c     Srdan Simunovic 
c     Oak Ridge National Lab.
c     P.O. box 2008, MS-6359
c     phone: (865) 574-4837, (865) 241-3863
c     email: leeh@ornl.gov, simunovics@ornl.gov
c
c    Inputs:
c    ------
c
c     cm      - User Shell Material property array
c     ldcm    - Leading dimension of cm array
c     capa    - reduction factor for transverse shear
c     dt1     - current time step size
c     time    - current problem time
c     tsh     - time increment for the element
c     ndsh(4) - node array for the element
c
c     deps(1) - local x  strain increment
c     deps(2) - local y  strain increment
c     deps(3) - local z  strain increment
c     deps(4) - local xy strain increment
c     deps(5) - local yz strain increment
c     deps(6) - local zx strain increment
c
c     ldeps   - leading dimension of deps array
c
c   Outputs
c   -------
c
c     sig(1) - local x  stress
c     sig(2) - local y  stress
c     sig(3) - local z  stress
c     sig(4) - local xy stress
c     sig(5) - local yz stress
c     sig(6) - local zx stress
c 
c     ldsig  - leading dimension of sig array
c
c     hisv(i) - ith history variable
c     fail    - element failure parameter
c
      dimension  cm(ldcmw),deps(ldeps),sig(ldsig),hisv(ldhv),ndsh(4),
     x           dfl(20)
c
c    Local Variables
c    ---------------
c
      data itermax /30/
      data Pi /3.1415927d0/
      data toler/1.0e-10/
      data  onethird /0.333333333d0/
c
      dimension dsig(6),devsig(6),devsigtn(6),tsig(6),sigtn(6),
     x epstn(6),cmat(6,6),dcdcb(6,6),work(20)
c
      logical elastic,converge
c
c Common blocks
c -------------
c
      common /temp1/ncyc
c 
c+++++++++++++++++++++++++Start Model Computations +++++++++++++++++++++
c
      delta_T = dt1
c
c --- Copy material props from cm array into local variables...
c
      rho = cm(1)
      shearm = cm(2)
      bulkm = cm(3)
      c0 = cm(4)
      dN0 = cm(5)
      dmu = cm(6)
      Tau0 = cm(7)
      surfe = cm(8)
      vsf = cm(9)
      dmmax = cm(10)
      dtmin = cm(11)
c
c --- Back out a poisson's ratio and youngs modulus from input
c     shear and bulk moduli....
      dnu = 0.5d0*(3.d0*bulkm - 2.d0*shearm)
      dnu = dnu/(3.d0*bulkm + shearm)
      youngm = 3.d0*bulkm*(1.d0-2.d0*dnu)
      e = youngm
      G = shearm
c
c --- Intialize some variables....
      do i = 1, 6
         tsig(i) = 0.d0
         do j = 1, 6
            cmat(i,j) = 0.d0
            dcdcb(i,j) = 0.d0
         enddo
      enddo
c
c --- Copy history varibales at previous time step to local variables...
c
      if (ncyc .eq. 0) then
c        open(unit=120,file='stress.dat')
         cn = c0
      else 
         cn = hisv(1)
      endif
      Prestn = hisv(2)
      do i = 1, 6
         devsigtn(i) = hisv(4+i)
         sigtn(i) = hisv(10+i)
         epstn(i) = hisv(16+i)
      enddo
c
c --- Compute the current moduli...
      f = dN0*cn*cn*cn
      dnucap = dnu
      call dsmod(f,E,dnu,Ecap,dnucap)
c
c --- Compute the current modulus matrix....
c
      c11 = Ecap*(1.d0 - dnucap)/((1.d0+ dnucap)*(1.d0 - 2.d0*dnucap))
      c12 = Ecap*dnucap/((1.d0+ dnucap)*(1.d0 - 2.d0*dnucap))
      c44 = Ecap/(2.d0*(1.d0 + dnucap))
c
      cmat(1,1) = c11
      cmat(2,2) = c11
      cmat(3,3) = c11
      cmat(1,2) = c12
      cmat(1,3) = c12
      cmat(2,1) = c12
      cmat(2,3) = c12
      cmat(3,1) = c12
      cmat(3,2) = c12
      cmat(4,4) = c44
      cmat(5,5) = c44
      cmat(6,6) = c44
c
c --- Compute the 3rd direction incremental strain
c     assuming no crack growth....
c
      r8tmp = -cmat(3,1)*deps(1) - cmat(3,2)*deps(2) 
      deps(3) = r8tmp/cmat(3,3)
c
c --- Compute incremental stresses.
      do i = 1, 6
         dsig(i) = 0.d0
         do j = 1, 6
            xx = deps(j)
c     already engineering strain
c            if ( j .gt. 3) xx = 2.d0*deps(j)
            dsig(i) = dsig(i) + cmat(i,j)*xx
         enddo
      enddo
c
c --- Compute incremental pressure...
      delP = -onethird*(dsig(1)+dsig(2)+dsig(3))
c
c --- Compute trial pressure and deviatoric stresses....
c
      Pres = Prestn + delP
c
      devsig(1) = devsigtn(1) + (dsig(1) + delP)
      devsig(2) = devsigtn(2) + (dsig(2) + delP)
      devsig(3) = devsigtn(3) + (dsig(3) + delP)
      devsig(4) = devsigtn(4) + dsig(4)
      devsig(5) = devsigtn(5) + dsig(5)
      devsig(6) = devsigtn(6) + dsig(6)
c
c    --- compute q = (1.5*S:S)^1/2
      q = dip2tens1(3,devsig(1))
      q = sqrt(1.5d0*q)
c
c    --- Check the stress state, (p,q), is within or outside the
c        damage surface...
      BigK = 0.5d0*Pi*((2.d0 - dnu)/(1.d0 - dnu))*surfe*G
c
      pp = Prestn
      if (Prestn .le. 0.d0) then
         Fd = (2.d0 - dnu)*pp*pp - 2.d0*BigK/cn
         Fd = -45.d0*Fd/(4.d0*(5.d0 - dnu))
      else
         Fd = (dmu*pp + Tau0)*(dmu*pp + Tau0 + sqrt(Pi*BigK/cn))
         Fd = Fd + BigK/cn
         Fd = 45.d0*Fd/(2.d0*(3.d0 - 2.d0*dmu*dmu))
      endif
c
c    --- If the stress state is within the damage surface then
c        an elastic response is obtained. Jump to end of routine. 
      elastic = .false.
      ii = 0
      q2 = q*q
      if (q2 .le. Fd) then 
         elastic = .true. 
         ii = 1
         goto 777
      else
         Ds = q2 - Fd
         delta = dtanh(Ds)
      endif
c
c --- The micro-cracks grow in this step and damage accumulates..
c
c --- First evolve the micro-crack's radii.....
c     
      f = dN0*cn*cn*cn
      cbar = cn
      iter = 0
      converge = .false.
c
 100  continue
c
      cbar3 = cbar*cbar*cbar
      cbar2 = cbar*cbar
c
c --- Update the micro-crack volume fraction...
      f = dN0*cbar3

c --- Compute the new moduli....
      call dsmod(f,E,dnu,Ecap,dnucap)

c --- Compute the function...
      vsf1 = 1.d0
      func = cbar - cn - vsf1*vsf*delta_T*delta*sqrt(Ecap/rho)
c
c --- Check the residual...
      resd0 = 0.d0
      if (iter .ne. 0) then
         if (abs(func) .le. toler) then 
            converge = .true.
            goto 200
         endif
      endif
      iter = iter + 1
      if (iter .gt. itermax) goto 200
c
c --- Compute the Jacobian....
      dfdcbar = 3.d0*dN0*cbar2
      dEdnucap = (E/dnu**(10/9))*(10.d0/9.d0)*(dnucap**(1/9))*
     x                 ( ((3.d0-dnu)/(3.d0-dnucap))**(1/9) )
      dEdnucap = dEdnucap + E*((dnucap/dnu)**(10/9))*((3.d0-dnu)**(1/9))
     x        *(-1.d0/9.d0)*( (3.d0-dnucap)**(-10/9) )
c
      c1 = (5.d0/8.d0)
      c2 = (15.d0/64.d0)
      c3 = (45.d0/128.d0)
      c4 = (5.d0/128.d0)
c
      dnucapdf = c1*(dnucap/dnu)*(-dnu/dnucap**2) + 
     x              c2*((1.d0-dnu)/(1.d0-dnucap))*(-1.d0/(1.d0-dnu)) + 
     x          c3*((1.d0+dnu)/(1.d0+dnucap))*(1.d0/(1.d0+dnu)) + 
     x          c4*((3.d0-dnu)/(3.d0-dnucap))*(-1.d0/(3.d0-dnu))
      dnucapdf = 1.d0/dnucapdf
      
      dEdcbar = dEdnucap*dnucapdf*dfdcbar

      dJac =  1.d0 - vsf1*vsf*delta_T*delta*0.5d0*
     x                (1.d0/sqrt(rho*Ecap))*dEdcbar
c
c --- Compute the correction to cbar....
      delta_cb = -func/dJac

c --- Update cbar...
      cbar = cbar + delta_cb
c
      goto 100

200   if (.not.converge) then
         write(*,*)"**** Error: solution did not converge"
         stop
      endif
c
      dfdcbar = 3.d0*dN0*cbar2
      dEdnucap = (E/dnu**(10/9))*(10.d0/9.d0)*(dnucap**(1/9))*
     x                 ( ((3.d0-dnu)/(3.d0-dnucap))**(1/9) )
      dEdnucap = dEdnucap + E*((dnucap/dnu)**(10/9))*((3.d0-dnu)**(1/9))
     x        *(-1.d0/9.d0)*( (3.d0-dnucap)**(-10/9) )
c
      c1 = (5.d0/8.d0)
      c2 = (15.d0/64.d0)
      c3 = (45.d0/128.d0)
      c4 = (5.d0/128.d0)
c
      dnucapdf = c1*(dnucap/dnu)*(-dnu/dnucap**2) + 
     x              c2*((1.d0-dnu)/(1.d0-dnucap))*(-1.d0/(1.d0-dnu)) + 
     x          c3*((1.d0+dnu)/(1.d0+dnucap))*(1.d0/(1.d0+dnu)) + 
     x          c4*((3.d0-dnu)/(3.d0-dnucap))*(-1.d0/(3.d0-dnu))
      dnucapdf = 1.d0/dnucapdf
c
      dEdcbar = dEdnucap*dnucapdf*dfdcbar
      dnudcbar = dnucapdf*dfdcbar
c
c --- Compute the new degraded mat. modulus matrix...
      c11 = Ecap*(1.d0 - dnucap)/((1.d0+ dnucap)*(1.d0 - 2.d0*dnucap))
      c12 = Ecap*dnucap/((1.d0+ dnucap)*(1.d0 - 2.d0*dnucap))
      c44 = Ecap/(2.d0*(1.d0 + dnucap))
c
      cmat(1,1) = c11
      cmat(2,2) = c11
      cmat(3,3) = c11
      cmat(1,2) = c12
      cmat(1,3) = c12
      cmat(2,1) = c12
      cmat(2,3) = c12
      cmat(3,1) = c12
      cmat(3,2) = c12
      cmat(4,4) = c44
      cmat(5,5) = c44
      cmat(6,6) = c44
c
c --- Compute the derivative of C matrix with cbar....
      fac3 = (1.d0 - dnucap)/((1.d0+dnucap)*(1.d0-2.d0*dnucap))
      fac4 = dnucap/((1.d0+dnucap)*(1.d0-2.d0*dnucap))
      fac5 = ((1.d0+dnucap)*(1.d0-2.d0*dnucap))**2
c
      r8tmp = -dnudcbar*(1.d0+dnucap)*(1.d0-2.d0*dnucap) -
     x          (1.d0-dnucap)*(dnudcbar*(1.d0-2.d0*dnucap) +
     x                (1.d0+dnucap)*(-2.d0*dnudcbar))
      dcdcb(1,1) = dEdcbar*fac3 + Ecap*r8tmp/fac5
      dcdcb(2,2) = dcdcb(1,1)
      dcdcb(3,3) = dcdcb(1,1)

      r8tmp = dnudcbar*(1.d0+dnucap)*(1.d0-2.d0*dnucap) -
     x          dnucap*(dnudcbar*(1.d0-2.d0*dnucap) + 
     x                  (1.d0+dnucap)*(-2.d0*dnudcbar))
      dcdcb(1,2) = dEdcbar*fac4 + Ecap*r8tmp/fac5
      dcdcb(1,3) = dcdcb(1,2)
      dcdcb(2,1) = dcdcb(1,2)
      dcdcb(2,3) = dcdcb(1,2)
      dcdcb(3,1) = dcdcb(1,2)
      dcdcb(3,2) = dcdcb(1,2)
c
      r8tmp = 0.5d0/((1.d0+dnucap)*(1.d0+dnucap))
      dcdcb(4,4) = dEdcbar*r8tmp*(1.d0+dnucap) - Ecap*r8tmp*dnudcbar
      dcdcb(5,5) = dcdcb(4,4)
      dcdcb(6,6) = dcdcb(4,4)
c
c --- Compute the 3rd direction incremental strain...
c
      r8tmp = - cmat(3,1)*deps(1) - cmat(3,2)*deps(2) 
     x           - dcdcb(3,1)*(cbar-cn)*epstn(1)
     x           - dcdcb(3,2)*(cbar-cn)*epstn(2)
     x           - dcdcb(3,3)*(cbar-cn)*epstn(3)
     x           - dcdcb(3,1)*(cbar-cn)*deps(1) - 
     x             dcdcb(3,2)*(cbar-cn)*deps(2)
      deps3o=deps(3)
      deps(3) = r8tmp/(cmat(3,3) + dcdcb(3,3)*(cbar-cn))
      deps3d=abs((deps3o-deps(3))/deps3o)
      if(deps3d .gt. 0.1)then
c        write(*,*)'deps3',deps3o,deps(3),deps3d
      endif
c
c --- Compute the new stresses....
      do i = 1, 6
         work(i) = 0.d0
         do j = 1, 6
            xx = epstn(j)+deps(j)
c$$$            if (j .gt. 3) xx = 2.d0*(epstn(j)+deps(j))
            work(i) = work(i) + dcdcb(i,j)*(cbar-cn)*xx
c            write(7,*)time, work(2)
         enddo
      enddo
      do i = 1, 6
         do j = 1, 6
            xx = deps(j)
c$$$            if (j .gt. 3) xx = 2.d0*deps(j)
            work(i) = work(i) + cmat(i,j)*xx
c            write(7,*)time, work(2)            
         enddo
         tsig(i) = sigtn(i) + work(i)
c         write(7,*)time,tsig(2)
      enddo
c
c --- Move the pressure and other values at time, tn.
      Prestn = -onethird*(tsig(1)+tsig(2)+tsig(3))
      cn = cbar
      do i = 1, 6
         fac1 = 0.d0
         if (i .eq. 1 .or. i .eq. 2 .or. i .eq. 3) fac1 = 1.d0
         devsigtn(i) = tsig(i) + fac1*Prestn
         sigtn(i) = tsig(i)
      enddo
c
777   if (elastic) then
         Prestn = Pres
         do i = 1, 6
            devsigtn(i) = devsig(i)
            fac1 = 0.d0
            if (i .eq. 1 .or. i .eq. 2 .or. i .eq. 3) fac1 = 1.d0
            sigtn(i) = devsigtn(i) - fac1*Prestn
         enddo
      endif
c
c --- Copy stresses into input sig array
      do i = 1, 6
         sig(i) = sigtn(i)
      enddo
c
      effs = 1./sqrt(2.)*sqrt((sig(1) - sig(2))**2 + 
     x       (sig(2) - sig(3))**2 + (sig(3) - sig(1))**2 + 
     x       6.*(sig(4)**2 + sig(5)**2 + sig(6)**2))
c
c  --- Update the total strain...
      do i = 1, 6
         epstn(i) = epstn(i) + deps(i)
      enddo
c
      dfl(1)=2.00001
      dfl(2)=cn
      dfl(3)=0.
      dfl(4)=dmmax
      dfl(5)=1.7
      dfl(6)=dtmin
      dfl(7)=tsh
      dfl(8)=time
      dfl(9)=hisv(4)
      dfl(10)=0.01
      dfl(11)=hisv(3)
      dfl(12)=0.
      dfl(13)=hisv(23)
      dfl(14)=effs

      call umfail(1,dfl,14)

      hisv(4)=dfl(9)
      hisv(3)=dfl(11)
      hisv(23)=dfl(14)
c      write(8,*)time,dfl(11),dfl(12)
      if(abs(dfl(12)) .ne. 0.0)then
         fail=hisv(3)
         write(7,*)time,fail
      else
         fail=1.
      endif
      
      hisv(3)=min(fail,hisv(3))
c
      do i=1,6
         sig(i)=fail*sig(i)
      enddo
c      write(9,*)time,fail
      
c
c --- Update history varibales....
      hisv(1) = cn
      hisv(2) = Prestn
c      hisv(3) = min(fail,hisv(3))
c      hisv(4) = dfl(11)
      do i = 1, 6
         hisv(4+i) = devsigtn(i)
         hisv(10+i) = sigtn(i)
         hisv(16+i) = epstn(i)
      enddo
c      write(9,*)time,fail
      
c
      q2 = dip2tens1(3,devsigtn(1))
      q2 = sqrt(1.5d0*q2)
c
c      if (f .ge. 1.5d0) goto 999
c
999   continue

      if(f .ge. 1.5d0)then
         write(*,*) 'bummer f',f
      endif
c      write(9,*)time,fail
c      fail=1.0

c      write(7,*)time,ecap
c     write(120,333)(deps(i),i=1,6)
c     write(120,333)f,cn,Ecap/E,dnucap/dnu,Prestn
333   format(6(e10.4,3x))
c     write(120,333)(epstn(i),i=1,6)
c     write(120,333)(sigtn(i),i=1,6)
c     write(120,333)(devsigtn(i),i=1,6)
c     write(120,*)'use',epstn(2),sigtn(2)
c     write(120,*)'ncyc',ncyc
#endif
      return
      end
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      function dip2tens1(if,s)
c
#ifdef ACC
#ifdef DP
      implicit double precision (a-h,o-z)                                    dp
#endif
c
      dimension s(*)
c
      if (if .eq. 1) then
         q = s(1)*s(1) + s(2)*s(2) + s(4)*s(4) + 2.d0*s(3)*s(3)
      elseif (if .eq. 3) then
         q = s(1)*s(1) + s(2)*s(2) + s(3)*s(3) + 2.d0*s(4)*s(4) + 
     x           2.d0*s(5)*s(5) + 2.d0*s(6)*s(6)
      endif
      dip2tens1 = q
c
#endif
      return
      end
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      function dip2tens(if,a,b)
c
#ifdef ACC
#ifdef DP
      implicit double precision (a-h,o-z)                                    dp
#endif
c
      dimension a(*),b(*)
c
      if (if .eq. 1) then
         q = a(1)*b(1) + a(2)*b(2) + a(4)*b(4) + 
     x          2.d0*a(3)*b(3)
      elseif (if .eq. 3) then
         q = a(1)*b(1) + a(2)*b(2) + a(3)*b(3) + 
     x          2.d0*a(4)*b(4) + 2.d0*a(5)*b(5) + 2.d0*a(6)*b(6)
      endif
      dip2tens = q
c
#endif
      return
      end
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      subroutine dsmod(f,e,dnu,ecap,dnucap)
c
#ifdef ACC
#ifdef DP
      implicit double precision (a-h,o-z)                                    dp
#endif
c
      data toler/1.0e-10/
      data itermax /30/
c
      logical converge
c
      converge = .false.
      iter = 0
c
100   continue
c
      func = (5.d0/8.d0)*log(dnu/dnucap) +
     x (15.d0/64.d0)*log((1.d0-dnucap)/(1.d0-dnu)) +
     x (45.d0/128.d0)*log((1.d0 + dnucap)/(1.d0+dnu)) + 
     x (5.d0/128.d0)*log((3.d0-dnucap)/(3.d0-dnu)) - f
c
      resd0 = 0.d0
      if (iter .ne. 0) then
         if (abs(func) .le. toler) then
            converge = .true.
            goto 200
         endif
      endif
      iter = iter + 1
      if (iter .gt. itermax) goto 200

      
      dJac = -(5.d0/8.d0)*(dnucap/dnu)*(dnu/dnucap**2) -
     x   (15.d0/64.d0)*((1.d0-dnu)/(1.d0-dnucap))*(1.d0/(1.d0-dnu)) + 
     x   (45.d0/128.d0)*((1.d0+dnu)/(1.d0+dnucap))*(1.d0/(1.d0+dnu)) -
     x   (5.d0/128.d0)*((3.d0-dnu)/(3.d0-dnucap))*(1.d0/(3.d0-dnu))


      delta_nu = -func/dJac
      dnucap = dnucap + delta_nu

      goto 100

200   if (.not.converge) then
         write(*,*)"**** Error: solution did not converge (dsmod)"
         stop
      endif
c
      Ecap = (dnucap/dnu)**(10.d0/9.d0)
      Ecap = Ecap*(((3.d0-dnu)/(3.d0-dnucap))**(1.d0/9.d0))
      Ecap = Ecap*E
      
c
#endif
      return
      end
