      subroutine umat49v(cm,d1,d2,d3,d4,d5,d6,sig1,sig2,
     .sig3,sig4,sig5,sig6,ef,hisv,lft,llt,dt1,capa,etype,
     .time,ipt,temp)
c
c     user material model for belytschko-tsay shell
c     ACC damage model 53
c
c
c     developed by:
c
c     Srdan Simunovic
c     Yongyi Zhu
c     Oak Ridge National Laboratory
c     p.o. box 2008, ms-6025
c     phone: (865) 241-3863
c     fax: (865) 574-7463
c     email: simunovics@ornl.gov
c
c     ACC damage model 53
c
c     modified by:
c
c     Venkatesh Agaram
c     Chrysler Corporation
c     CIMS 483-05-10, 800 Chrysler Dr East
c     Auburn Hills,MI-48326 2757
c     phone: (248) 576 5655
c     fax  : (248) 576 7936
c     email: VA2@CHRYSLER.COM
c
c
c     variables
c
c     cm(1) - young's modulus
c     cm(2) - poisson's ratio
c     cm(3) - tensile strength
c     cm(4) - compressive strength
c     cm(5) - tensile m (6-8, tensile stress-strain curve fit)
c     cm(6) - compressive m (0.5-1, lowered for delamination)
c     cm(7) - maximum damage
c     cm(8) - minimum time step for shell
c     cm(9) - damage viscosity parameter
c
c     d1(nlq) - local x  strain increment
c     d2(nlq) - local y  strain increment
c     d3(nlq) - local z  strain increment
c     d4(nlq) - local xy strain increment
c     d5(nlq) - local yz strain increment
c     d6(nlq) - local zx strain increment
c
c     sig1(nlq) - local x  stress
c     sig2(nlq) - local y  stress
c     sig3(nlq) - local z  stress
c     sig4(nlq) - local xy stress
c     sig5(nlq) - local yz stress
c     sig6(nlq) - local zx stress
c
c     ef(nlq)     - fiber damage parameter wt+wc
c     hisv(nlq,i) - ith history variable
c
c     lft - first element in vectorized form
c     llt - last element in vectorized form
c
c     dt1   - current time step size
c     capa  - reduction factor for transverse shear
c     etype - element type
c     time  - current problem time
c     ipt   - Gauss point number 
c
c
      include 'nlqparm'
c      parameter (nlq=256)	
cdir$ taskcommon failuloc, soundloc, aux35loc	
c
      character*(*) etype
c
      dimension cm(1),d1(1),d2(1),d3(1),d4(1),d5(1),d6(1)
      dimension sig1(1),sig2(1),sig3(1),sig4(1),sig5(1),sig6(1)
      dimension ef(1),hisv(nlq,*),dt1(1)
c
      dimension efail(nlq),dtck(nlq)
c
      dimension dfl(20)
c
      logical failur,failgl
      common/failcm/failur,failgl
      common/failcmloc/ifail(nlq)
c$$$      common/failuloc/sieu(nlq),fail(nlq)

      common/failuloc/sieu(nlq),failu(nlq)

      common/soundloc/sndspd(nlq),sndsp(nlq),diagm(nlq),sarea(nlq),
     .                dxl(nlq)
      common/aux35loc/rhoa(nlq),cb(nlq),davg(nlq),p(nlq)
c
      common /mptcom/ nip,lavv,nmtcon,nnm1,iop


      data zero /0.d0/, one /1.d0/, two /2.d0/, half /0.5d0/,
     .  three/3.d0/,third /0.333333333d0/,twotrd/.66666666666667d0/,
     .  four/4.d0/,oneoe/0.36787944d0/,coep/2.7182818d0/,proc/1.d-6/
c
      if (etype.ne.'shell') then
         write(0,*) 'error: material 49v is for shells only'
         return
      end if
c
c     material parameters
c
      et   = cm(1)
      anu  = cm(2)
      xt   = cm(3)
      xc   = cm(4)
      dmt  = cm(5)
      dmc  = cm(6)
      dmmax= cm(7)
      dtmin= cm(8)
      visco= cm(11)
      twog = et/(one+anu)
      ranu = one/(one-anu)
      canu = anu*ranu
      proce = proc*(xt+xc)/et

c
      if (time.eq.0.0) then
        do i=lft,llt
           efail(i)=1.0
           hisv(i,1)=0.
           hisv(i,2)=0.
           hisv(i,3)=0.
           hisv(i,4)=0.
           hisv(i,5)=0.
           hisv(i,6)=0.
           hisv(i,7)=0.
           hisv(i,11)=1.0
        end do
      else
        do i=lft,llt
           efail(i)=hisv(i,11)
        end do
      end if
c
      if (ipt.eq.1.and.time.eq.0.0) then
        do i=lft,llt
          failu(i)=1.0
        end do
      end if
c
      do i=lft,llt
c
c     compute total strain increments
c
      v1 = d1(i)
      v2 = d2(i)
      v12= d4(i)
c
c     v12= 2.0*d4(i)
c
      wt = hisv(i,1)
      wc = hisv(i,2)
      rt = hisv(i,3)
      rc = hisv(i,4)
c
      wt0 = wt
      wc0 = wc
      rt0 = rt
      rc0 = rc
c
      rwt = one - wt
      rwc = one - wc
c
c     identification of stress
c
c     compute new strains
c
      eps1 = hisv(i,5) + v1
      eps2 = hisv(i,6) + v2
      eps12= hisv(i,7) + v12
c
      hisv(i,5) = eps1
      hisv(i,6) = eps2
      hisv(i,7) = eps12
c
c     strain tensor split
c

c$$$      write(80,'(a4,e10.3,1x,e10.3,1x,e10.3,1x,e10.3,1x,e10.3)')
c$$$     $     'str ',time,eps1,eps2,eps12

      coef1 = sqrt((eps1-eps2)*(eps1-eps2)+eps12*eps12)
      rm1 = (eps1+eps2) + coef1
      rm1 = half*rm1
      rm2 = (eps1+eps2) - coef1
      rm2 = half*rm2
      coefc = sqrt((eps1-eps2+coef1)**2+eps12*eps12)
      if (coefc.gt.proce) then
         a1 = (eps1-eps2 + coef1)/coefc
         a2 = eps12/coefc
      else
c        if ((rm1*rm2).lt.0..and.rm2.lt.rm1) then
c           a1 = zero
c           a2 = one
c        else
            a1 = one
            a2 = zero
c        end if
      end if
      if (rm1.lt.zero) then
         crm1p = zero
         crm1n = one
      else
         crm1p = one
         crm1n = zero
      end if
      if (rm2.lt.zero) then
         crm2p = zero
         crm2n = one
      else
         crm2p = one
         crm2n = zero
      end if
c     crmps= canu
c     crmns= canu
      crmp = canu*crm1p*crm2p
      crmn = canu*crm1n*crm2n
c
c     updated positive strain compoments
c
      ep1 = crm1p*rm1*a1*a1 + crm2p*rm2*a2*a2
      ep2 = crm1p*rm1*a2*a2 + crm2p*rm2*a1*a1
      ep12= (crm1p*rm1-crm2p*rm2)*a1*a2
c
c     updated negative strain compoments
c
      en1 = eps1-ep1
      en2 = eps2-ep2
      en12= half*eps12-ep12
      epm = ep1+ep2
      enm = en1+en2

c$$$      write(80,'(a4,e10.3,1x,e10.3,1x,e10.3,1x,e10.3,1x,e10.3)')
c$$$     $     'epm ',time,ep1,ep2,ep12
c$$$      write(80,'(a4,e10.3,1x,e10.3,1x,e10.3,1x,e10.3,1x,e10.3)')
c$$$     $     'enm ',time,en1,en2,en12

c
c     calculate equivalent positive/negative strain
c
      epeq = ep1*ep1+ep2*ep2+two*ep12*ep12+canu*epm*epm
      eneq = en1*en1+en2*en2+two*en12*en12+canu*enm*enm
c
c     calculate damage criteria in strain space
c
      gt = (twog/xt)**2*epeq*(one-anu)
cc    gc = (twog/xc)**2*eneq
      gc = (twog/xc)**2*eneq*(one-anu)
c
c     vector normal to damage surface
c
      gnt1 = ep1*ranu+crmp*ep2
      gnt2 = ep2*ranu+crmp*ep1
      gnt  = gnt1*v1+gnt2*v2+ep12*v12
      gnc1 = en1*ranu+crmn*en2
      gnc2 = en2*ranu+crmn*en1
      gnc  = gnc1*v1+gnc2*v2+en12*v12
c
c     gnt = zero
c     if (crm1p.gt.zero.or.crm2p.gt.zero) then
c        gnt1 = ep1*ranu + crmps*ep2
c        gnt2 = ep2*ranu + crmps*ep1
c        gnt = gnt1*v1 + gnt2*v2 + ep12*v12
c     end if
c     gnc = zero
c     if (crm1n.gt.zero.or.crm2n.gt.zero) then
c        gnc1 = en1*ranu + crmns*en2
c        gnc2 = en2*ranu + crmns*en1
c        gnc  = gnc1*v1 + gnc2*v2 + en2*v12
c     end if
c
c     verify loading/unloading and calculate damage variables
c
      if ((gt-rt).ge.zero .and. gnt.ge.zero) then
c
c     loading, calculate new damage variable
c
         if (wt.lt.dmmax) then
            rt = gt
            coef6 = -oneoe/dmt*gt**(dmt/two)

c$$$            call chk_class(coef6,iclass)
c$$$            if(iclass.eq.4.or.iclass.eq.5.or.
c$$$     $           iclass.eq.8.or.iclass.eq.9)then
c$$$               write(*,*) 'nanq 1',iclass,coef6
c$$$               pause
c$$$            end if

            wt = max(wt,one-exp(coef6))
            if (wt.gt.dmmax) wt = (one+1.d-5)*dmmax
c$$$      write(80,'(a4,e10.3,1x,e10.3,1x,e10.3,1x,e10.3,1x,e10.3)')
c$$$     $     'tdm ',time,gt,rt,gnt,wt

         end if
c
      end if
c
      if ((gc-rc).ge.zero .and. gnc.ge.zero) then
c
c     loading, calculate new damage variable
c
         if (wc.lt.dmmax) then
            rc = gc
            coef6 = -oneoe/dmc*gc**(dmc/two)

c$$$            call chk_class(coef6,iclass)
c$$$            if(iclass.eq.4.or.iclass.eq.5.or.
c$$$     $           iclass.eq.8.or.iclass.eq.9)then
c$$$               write(*,*) 'nanq 2',iclass,coef6
c$$$               pause
c$$$            end if

            wc = max(wc,one-exp(coef6))
            if (wc.gt.dmmax) wc = (one+1.d-5)*dmmax
c$$$      write(80,'(a4,e10.3,1x,e10.3,1x,e10.3,1x,e10.3,1x,e10.3)')
c$$$     $     'cdm ',time,gc,rc,gnc,wc

         end if
c
      end if
c
c
c     viscous regularization
c      
      if(visco.gt.zero) then
       dtvis1 = dt1(i)/visco
c         dtvis1 = 1./visco
         dtvis2 = one + dtvis1
c
         wt  = (wt0 + dtvis1*wt)/dtvis2
         wc  = (wc0 + dtvis1*wc)/dtvis2
         rt = (rt0 + dtvis1*rt)/dtvis2
         rc = (rc0 + dtvis1*rc)/dtvis2
      endif
c
      rwt = (one-wt)*twog
      rwc = (one-wc)*twog
      g12  = twog*anu/(one-anu)
c
c     calculate new stresses
c
c     sig1(i) = rwt*(ep1*ranu+ep2*crmp)+rwc*en1+g12*(en1+crmn*en2)
c     sig2(i) = rwt*(ep2*ranu+ep1*crmp)+rwc*en2+g12*(en2+crmn*en1)
c
      sig1(i) = rwt*(ep1*ranu+ep2*crmp)+rwc*(en1*ranu+en2*crmn)
      sig2(i) = rwt*(ep2*ranu+ep1*crmp)+rwc*(en2*ranu+en1*crmn)
      sig3(i) = zero
      sig4(i) = rwt*ep12+rwc*en12
      wq4 = twog/2.0*capa
      sig5(i) = sig5(i)+wq4*d5(i)
      sig6(i) = sig6(i)+wq4*d6(i)
c
      d3(i) = -canu*(eps1+eps2)
c
      dfl(1)=2.00001
      dfl(2)=wt
      dfl(3)=wc
      dfl(4)=dmmax
      dfl(5)=1.5
      dfl(6)=dtmin
      if (dtmin.gt.zero) then
         snds=sqrt(sndspd(lft)/rhoa(lft))
         dfl(7)=sarea(i)/(sqrt(diagm(i))*snds)
      else
         dfl(7)=zero
      end if      
      dfl(8)=time
      dfl(9)=hisv(i,12)
      dfl(10)=0.01
      dfl(10)=0.25
      dfl(11)=hisv(i,11)
      dfl(12)=0.
c	
      call umfail(1,dfl,12)
c
      hisv(i,12)=dfl(9)
      hisv(i,11)=dfl(11)
c
      if (dfl(12).ne.0.)then
         efail(i)=hisv(i,11)
      else
         efail(i)=1.
      endif
c
      hisv(i,11)=min(efail(i),hisv(i,11))
c
c$$$      sig1(i) = efail(i) * sig1(i)
c$$$      sig2(i) = efail(i) * sig2(i)
c$$$      sig3(i) = efail(i) * sig3(i)
c$$$      sig4(i) = efail(i) * sig4(i)
c$$$      sig5(i) = efail(i) * sig5(i)
c$$$      sig6(i) = efail(i) * sig6(i)

c
      hisv(i,1) = wt
      hisv(i,2) = wc
      hisv(i,3) = rt
      hisv(i,4) = rc
c 
      if (ipt.eq.1) then
         failu(i)=efail(i)
      else    
         failu(i)=max(failu(i),efail(i))
      end if

      if(ipt.eq.5)then
         if(failu(i).eq.zero)then
            failur=.true.
            ifail(i)=1
c$$$            ieint=nnm1+i
c$$$            ieext=nelmntid(ieint,2)
c$$$      write(80,'(a4,e10.3,1x,i10,1x,e10.3,1x,e10.3,1x,i10,1x,i10)')
c$$$     $     'flr ',time,i,failu(i),efail(i)
         endif
      endif
            

c$$$      write(80,'(a4,e10.3,1x,i10,1x,e10.3,1x,e10.3,1x,e10.3,1x,e10.3)')
c$$$     $     'flu ',time,ipt,failu(i),efail(i)

      ef(i)=wt+wc

c
      end do
c
      return
      end
c
      subroutine umfail(failaw,a,n)
c
c     user material model for belytschko-tsay shell
c
c     srdan simunovic
c     lockheed martin energy research corporation
c     p.o. box 2008, ms-6140
c     phone: (423) 241-3863
c     fax: (423) 574-4839
c     email: qs9@ornl.gov
c
c     disclaimer: your mileage may vary.
c
      implicit none
c
      integer failaw,n,i,ndm,imax
c#ifdef DP
c      double precision a(n),dtot,dmax,zero,one
c#else
      real a(n),dtot,dmax,zero,one
c#endif
c
      data zero /0.d0/, one /1.d0/
c      
      if(failaw.eq.1)then
c........................
c
c
c     a(1)       - # of damage variables (in)
c     a(1+ndm)   - damage variables (in)
c     a(1+ndm+1) - max. allowed damage (in)
c     a(1+ndm+2) - factor for max. damage (in)
c     a(1+ndm+3) - min. time (in)
c     a(1+ndm+4) - el. time (in)
c     a(1+ndm+5) - problem time (in)
c     a(1+ndm+6) - time when deletion starts (in/out)
c     a(1+ndm+7) - increment of failure (in)
c     a(1+ndm+8) - current failure level (in/out)
c     a(1+ndm+9) - failure indicator (-1=dt, 1-ndm= d. v.) (out)
c
         ndm=nint(a(1))
         dtot=0.
         dmax=0.
         imax=0
         a(1+ndm+9)=0.
         do i=1,ndm
            dtot=dtot+a(1+i)
            dmax=max(dmax,a(1+i))
            if(dmax.eq.a(1+i))then
               imax=i
            endif
         enddo
c
         if((a(1+ndm+4).gt.zero) .and. (a(1+ndm+4).le.a(1+ndm+3)))then
            if(a(1+ndm+6).eq.0.)then
               a(1+ndm+6)=a(1+ndm+5)
               a(1+ndm+8)=1.
            endif
c           a(1+ndm+8)=0.
            a(1+ndm+8)=a(1+ndm+8)-0.1
            a(1+ndm+8)=max(zero,a(1+ndm+8))
            a(1+ndm+9)=-1.
            return
         endif
c
         if(dmax.gt.a(1+ndm+1))then
            if(dtot.gt.(a(1+ndm+2)*a(1+ndm+1)))then
               a(1+ndm+9)=real(imax)
               if(a(1+ndm+6).eq.0.)then
                  a(1+ndm+6)=a(1+ndm+5)
                  a(1+ndm+8)=1.
               endif
               a(1+ndm+8)=a(1+ndm+8)-a(1+ndm+7)
               a(1+ndm+8)=max(zero,a(1+ndm+8))
            endif
         endif
c
      else if(failaw.eq.2)then
c.............................
c
c
c     a(1)       - # of damage variables (in)
c     a(1+ndm)   - damage variables (in)
c     a(1+ndm+1) - max. allowed damage (in)
c     a(1+ndm+2) - factor for max. damage (in)
c     a(1+ndm+3) - min. time (in)
c     a(1+ndm+4) - el. time (in)
c     a(1+ndm+5) - current failure level (in/out)
c     a(1+ndm+6) - failure indicator (-1=dt, 1-ndm= d. v.) (out)
c
         ndm=nint(a(1))
         dtot=0
         dmax=0
         imax=0
         a(1+ndm+6)=0
         do i=1,ndm
            dtot=dtot+a(1+i)
            dmax=max(dmax,a(1+i))
            if(dmax.eq.a(1+i))then
               imax=i
            endif
         enddo
c
         if(a(1+ndm+4).le.a(1+ndm+3))then
            a(1+ndm+6)=-1.
            a(1+ndm+5)=0
            return
         endif
c
         if(dmax.gt.(a(1+ndm+2)*a(1+ndm+1)))then
            a(1+ndm+6)=real(imax)
            a(1+ndm+5)=(a(1+ndm+1)-dmax)/(a(1+ndm+1)-
     $           (a(1+ndm+2)*a(1+ndm+1)))
            a(1+ndm+5)=max(zero,a(1+ndm+5))
         endif
c
      else if(failaw.eq.3)then
c.............................
c
c
c     a(1)       - # of damage variables (in)
c     a(1+ndm)   - damage variables (in)
c     a(1+ndm+1) - max. allowed cumulative damage (in) > 1.
c     a(1+ndm+2) - factor for max. damage (in) < 1.
c     a(1+ndm+3) - min. time (in)
c     a(1+ndm+4) - el. time (in)
c     a(1+ndm+5) - current failure level (in/out)
c     a(1+ndm+6) - failure indicator (-1=dt, 1-ndm= d. v.) (out)
c
         ndm=nint(a(1))
         dtot=0
         dmax=0
         imax=0
         a(1+ndm+6)=0
         do i=1,ndm
            dtot=dtot+a(1+i)
            dmax=max(dmax,a(1+i))
            if(dmax.eq.a(1+i))then
               imax=i
            endif
         enddo
c
         if(a(1+ndm+4).le.a(1+ndm+3))then
            a(1+ndm+6)=-1.
            a(1+ndm+5)=0
            return
         endif
c
         if(dtot.gt.(a(1+ndm+2)*a(1+ndm+1)))then
            a(1+ndm+6)=real(imax)
            a(1+ndm+5)=(a(1+ndm+1)-dtot)/(a(1+ndm+1)-
     $           (a(1+ndm+2)*a(1+ndm+1)))
            a(1+ndm+5)=max(zero,a(1+ndm+5))
         endif
c
      else
         write(*,*) 'incorrect failure law', failaw
      endif
c
      return
      end
