subroutine umat(stress,statev,ddsdde,sse,spd,scd,rpl,ddsddt, 
     1 drplde,drpldt,stran,dstran,time,dtime,temp,dtemp,predef,dpred,
     2 cmname,ndi,nshr,ntens,nstatv,props,nprops,coords,drot,
     3 pnewdt,celent,dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)
c
      include 'aba_param.inc'
c
      character*80 cmname
c
      dimension stress(ntens),statev(nstatv),ddsdde(ntens,ntens),
     1 ddsddt(ntens),drplde(ntens),stran(ntens),dstran(ntens),
     2 time(2),predef(1),dpred(1),props(nprops),coords(3),drot(3,3),
     3 dfgrd0(3,3),dfgrd1(3,3)
c ----------------------------------------------------------------------
      real*8 strain(3),eelas(3),eelas_n(3),eplas(3),eplas_n(3),
     &  aa(2,2),unitsecondorder(3),unitfourorder(3,3),stress_n(3),
     &  trialstress(3),devtristress(3),prindamage(2,2),
     &  backstress(3),backstress_n(3),damage_inc(3),damage(3),
     &  damage_n(3),seconddamage(2,2),prindameigen(2),hh(3,3),
     &  fourdamage(3,3),fourdamage_inv(3,3),fourdamage_inv_t(3,3),
     &  damagedir(2,2),damagedir_t(2,2),devstress(3), dnetstress(3),
     &  hh_damage(3,3),hh_damdam(3,3),mphi(3),fourorder_j(3,3),
     &  thermfor(3),thermfor_inc(3),plastnorm(3),fourordermult(3,3),
     &  devunitfourorder(3,3),plastnormmult1(3,3),plastnormmult2(3,3),
     &  plastnormmult3(3,3),temp1(2,2),temp2(3,3),temp3(3,3),temp4(3),
     &  temp5(3,3),temp6(3,3),temp7(3),temp8(3,3),temp9(3),temp10(3),
     &  temp11(3),temp12(3,3),temp13(3,3),temp14(3,3),temp15(3,3),
     &  temp16(3,3),temp17(3,3),temp18(3,3),fourdamage_t(3,3),
     &  plastnorm2(3)
c
       real*8 emod,enu,ebulk3,eg2,eg,elam,ebulk,friangle,dilangle,
     &  yieldo,dkesa1,dkesa11,sigmat,sigmac,fmiu,dkd,hb,dkesa2,alfa0,
     &  dm0,eqplas,eqplas_n,detdamage,hydtristress,equivstress,
     &  equivalentstress,damagenorm,hardstress,yieldcondition,
     &  yieldcondition1,yieldcondition2,trace1,trace2,datalamta,dlamta1,
     &  dlamta2,dlamta3,dlamta4,phy,hp,fyrhs,dlamta11,dlamta22,dlamta33,
     &  dlamta44,hydstress,yeq,datalamta_d1,datalamta_d2,datalamta_d,
     &  dlamta444,toler,limit,tempa1,tempa2,tempa3,equivalentstress2
c       
       integer mm,nn,k1,k2,newton,kewton

       parameter (zero=0.d0,one=1.d0,two=2.d0,three=3.d0,four=4.d0,
     &            nine=9.d0,newton=1000,toler=1.0d-6,mm=3,nn=2,
     &            limit=1.0d-4)
c ----------------------------------------------------------------------------------------
c  UMAT for Finite-deformation Anisotropic Drucker-prager Viscoplasticity/Damage
c  Coupled Thermodynamic Model with Isotropic Hardening/Softening for Plane Stress Problem
c       NTENS = 3 and NDI = 2 for Plane Stress Problem
c       props(1)----Elastic Modulus emod = 229
c       props(2)----Poisson Ratio enu = 0.45
c       props(3)----Frictional Angle of Adhesive Materials friangle = 0.5233
c       props(4)----Dilatant Angle of Adhesive Materials dilangle = 0.314
c       props(5)----Initial Yielding Stress of Adhesive Materials yieldo = 15.155
c       props(6)----Coefficient in Yielding Stress dkesa1 = 0.58466
c       props(7)----Coefficient in Yielding Stress dkesa11 = 11.534
c       props(8)----Tensile Strength of Adhesive Materials sigmat = 13.9
c       props(9)----Compressive Strength of Adhesive Materials sigmac = 22
c       props(10)----Anisotropic Coefficient in Equivalent Stress fmiu = 0.5 (-0.5 < fmiu < 1)
c       props(11)----Damage Coefficient in Damage Potential Function dkd = 1 (dkd > 0)
c       props(12)----Coefficient in Kinematic Hardening Four-order Modulus Tensor hb = 0.01 (hb > 0)
c       props(13)----Coefficient in Damage Potential Function dkesa2 = 0.01 (dkesa2 > 0)
c       props(14)----Viscous Coefficient in Viscoplastic Model alfa0 = 100 (alfa0 > 0)
c       props(15)----Viscous Coefficient in Viscoplastic Model dm0 = 0.2 (dm0 > 0)
c -----------------------------------------------------------------------------------------
c     Get User-defined Parameters
c
      emod     =   props(1)
      enu      =   props(2)   
      friangle =   props(3)   
      dilangle =   props(4)   
      yieldo   =   props(5)           
      dkesa1   =   props(6)  
      dkesa11  =   props(7)     
      sigmat   =   props(8)           
      sigmac   =   props(9)           
      fmiu     =   props(10)    
      dkd      =   props(11)     
      hb       =   props(12)     
      dkesa2   =   props(13)  
      alfa0    =   props(14)   
      dm0      =   props(15)     
c
c     eg----¦Ì,lame constants
c     ebulk----k = lamta + 2/3*¦Ì,lame constants
c
      ebulk3 = emod/(one-two*enu)
      eg2 = emod/(one+enu)
      eg = eg2/two
      elam = (ebulk3-eg2)/three
      ebulk = ebulk3/three 
c     
c  eelas - elastic strain in voigt form
c  eplas - plastic strain in voigt form
c  backstress - back stress in voigt form at the start of new iteration
c
      if((kstep.eq.1).and.(kinc.eq.0)) then   
         do k1 = 1,nstatv
             statev(k1) = zero
         end do
      end if  
c      
      do k1=1,3
        eelas_n(k1) = statev(k1)
        eplas_n(k1) = statev(k1+1*3)
        stress_n(k1) = statev(k1+2*3)
        backstress_n(k1) = statev(k1+3*3)
        eplas(k1) = statev(k1+1*3)
        backstress(k1) = statev(k1+3*3)
      end do
c
       ddsdde = ddsdde * zero

c      
       do k1 = 1,2
         do k2 = 1,2
           ddsdde(k1,k2) = elam
         end do
           ddsdde(k1,k1) = eg2+elam
       end do
       do k1 = 3,3
           ddsdde(k1,k1) = eg
      end do

c
c  engineering strain is used, so "eg" instead of "2*eg" used
c  backup for undamaged elastic matrix
c      
c**********************************************************************
c                          Damage Part
c**********************************************************************
c      
c  calculate principal damage eigenvalues and directions
c  seconddamage-second order damage tensor at the start of new iteration (end of last iteration)
c  restore second-order damage tensor at the end of last iteration
c
c  here,damage is given the number of last literation to prevent errors 
c      
       damage(1) = statev(14)
       damage(2) = statev(15)
       damage(3) = statev(16)
       damage_n(1) = statev(14)
       damage_n(2) = statev(15)
       damage_n(3) = statev(16)

c
c  only two principal damage directions, thus 2*2 damage tensors
c
       seconddamage(1,1) = damage_n(1)
       seconddamage(2,2) = damage_n(2)
       seconddamage(1,2) = damage_n(3)
       seconddamage(2,1) = damage_n(3)
c             
c  prindameigen-principal eigenvalues of second-order damage tensor by spectral decomposition
c  damagedir-eigenvectors of second-order damage tensor by spectral decomposition
c  damagedir_t-transpose of damagedir tensor
c  initial values for matrix fourdamageinv for later initial calculation without damage
c
       fourdamage_inv = zero
       fourdamage_inv(1,1) = one
       fourdamage_inv(2,2) = one
       fourdamage_inv(3,3) = one
       fourdamage_inv_t = transpose(fourdamage_inv)
       fourdamage = fourdamage_inv
       fourdamage_t = fourdamage_inv_t
c
       detdamage = zero
       call mdet(seconddamage,detdamage,nn)
       
c
       if (detdamage.gt.toler) then
          prindameigen = zero
          damagedir = zero
          damagedir_t = zero
          prindamage = zero
          aa = zero
          call spectral(seconddamage,prindameigen,damagedir,nn) 
c
c  damage space transforms to stress space  
c  by refering to "voyiadjis and park, 1997, eq.(29) 
c  anisotropic damage effect tensors for the symmetrization of the effective stress tensor"
c
c  prindamage second-order principal damage tensor
c      
       do k1 = 1,2
          prindamage(k1,k1) = one/sqrt(one-prindameigen(k1))
       end do
c  
c  damage space transforms to stress space  
c  by refering to "voyiadjis and park, 1997, eq.(32) 
c  anisotropic damage effect tensors for the symmetrization of the effective stress tensor"
c 
c  temp1-19: temporary array
c  aa-dyadic of two damagedir tensors in eq.(32):a=b*b,
c  b=damagedir,eigenvector of second-order damage tensor
c      
       temp1 = zero
       damagedir_t = transpose(damagedir)  
       temp1 = matmul(damagedir_t,prindamage)
       aa = matmul(temp1,damagedir)
c
c  plane stress problem which select,1,2,4 in transformation matrix 
c  by refering to "voyiadjis and park, 1997, eq.(36) 
c  anisotropic damage effect tensors for the symmetrization of the effective stress tensor"
c  for 3d problem, complete transformation matrix is required for change in eq.(36) 
c
c  fourdamageinv-inverse of four-order damage tensor, m(-1)for plane stress problem
c
       fourdamage(1,1) = aa(1,1)*aa(1,1)
       fourdamage(1,2) = aa(1,2)*aa(1,2)
       fourdamage(2,1) = aa(2,1)*aa(2,1)
       fourdamage(2,2) = aa(2,2)*aa(2,2)
       fourdamage(1,3) = two*aa(1,1)*aa(1,2)
       fourdamage(3,1) = aa(1,1)*aa(2,1)
       fourdamage(2,3) = two*aa(2,1)*aa(2,2)
       fourdamage(3,2) = aa(1,2)*aa(2,2)
       fourdamage(3,3) = aa(1,1)*aa(2,2)+aa(1,2)*aa(2,1)
       fourdamage_t = transpose(fourdamage)
       call m2inv(fourdamage,fourdamage_inv,mm)       
       fourdamage_inv_t = transpose(fourdamage_inv)
c       
c  update damage elastic modulus "ddsdde" symmetre matrix
c     
       temp2 = zero
       temp2 = matmul(fourdamage_inv_t,ddsdde)
       ddsdde = matmul(temp2,fourdamage_inv)
      end if
c**********************************************************************      
c       
c  trial stress at the start of new iteration
c
       trialstress = zero
       do k1 = 1,3
             trialstress(k1) = stress_n(k1)
          do k2 = 1,3
             trialstress(k1) = trialstress(k1)+
     &                      ddsdde(k1,k2)*dstran(k2)
          end do
       end do
c
c  trial elastic strain
c
       do k1 = 1,3
         eelas(k1) = eelas_n(k1)+dstran(k1)
      end do
c
c  update elastic strain and stress if no plasticity
c     
       do k1 = 1,3
         strain(k1) = eelas(k1)
         stress(k1) = trialstress(k1)
       end do
c
c	 trial hydrostatic strain and trial deviatoric strain for plasticity      
c
       hydtristress = (trialstress(1)+trialstress(2))/three
       do k1 = 1,2
         devtristress(k1) = trialstress(k1)-hydtristress
       end do
       devtristress(3) = trialstress(3)
c
       do k1 = 1,3
         dnetstress(k1) = devtristress(k1)-backstress_n(k1)
       end do
c
c
c  calculate the trail equivalent stress by anisotropic plasticity 
c  parameters for plane stress problem
c
       do k1 = 1,3
           do k2 = 1,3
               hh(k1,k2) = zero
           end do
      end do
c      
       !hh(1,1) = two*one 
       !hh(1,2) = two*fmiu 
       !hh(2,1) = two*fmiu 
       !hh(2,2) = two*one   
       !hh(3,3) = two*two*(one-fmiu)
       hh(1,1) = one 
       hh(2,2) = one   
       hh(3,3) = one
c
c  hh_damage=m^(-1)*hh*m^(-1),symmetre,four-order plasticity modulus for "equivalentstress"
c
       temp3 = zero
       hh_damage = zero
c       
       temp3 = matmul(fourdamage_t,hh)
       hh_damage = matmul(temp3,fourdamage)
c       temp3 = matmul(fourdamage_inv_t,hh)
c       hh_damage = matmul(temp3,fourdamage_inv)

       do k1 = 1,3
            temp4(k1) = zero
         do k2 = 1,3
            temp4(k1) = temp4(k1)+hh_damage(k1,k2)*dnetstress(k2)
         end do
       end do
c
       equivstress = zero
       do k1 = 1,3
           equivstress = equivstress+temp4(k1)*dnetstress(k1)
       end do
c 
       equivalentstress = sqrt(three/two*equivstress)  
c       
       damagenorm = sqrt(seconddamage(1,1)**2+seconddamage(1,2)**2+
     &                 seconddamage(2,1)**2+seconddamage(2,2)**2)    
c       
c**********************************************************************
c                         Plasticity Part
c**********************************************************************
c  calculate damage hardening stress in yielding function and
c  equivalent plastic strain "eqplas" at the end of last iteration
c  for plastic hardening stress   
c      yield condition checked!
       eqplas = statev(13)
       eqplas_n = statev(13)
       hardstress = yieldo*(one-dkesa1*exp(-dkesa11*eqplas_n))
c       
       yieldcondition1 = cos(friangle)*equivalentstress
c       
       yieldcondition2 = sin(friangle)*(one-damagenorm)*(sigmac-sigmat)/
     &                 sigmac/sigmat*sqrt(three)*hydtristress
c       
       yieldcondition = yieldcondition1+yieldcondition2-
     &                (one-damagenorm)*hardstress
c
c  judge yielding condition
c  solve plastic consistency factor "datalamta" for  
c  equivalent plastic strain increment using newton iterations
c
       temp12 = zero
       temp13 = zero
       temp14 = zero
       temp15 = zero
       temp16 = zero
       temp17 = zero
       temp18 = zero
       
      if (yieldcondition.gt.toler) then
c          
c  trace1 = tr((m-1)*(m-1)*hh_damage:1),hh_damdam=m^(-1)*hh_damage*m^(-1),symmetre  
c  trace1 checked! (squence)
c          
           temp5 = zero
           hh_damdam = zero
           temp5 = matmul(fourdamage_inv_t,fourdamage_inv)
           hh_damdam = matmul(temp5,hh_damage)
           trace1 = hh_damdam(1,1)+hh_damdam(2,2)
c           trace1 = one - damagenorm
           
c              
c  trace2 = tr((m-1)*(m-1):1)   
c  trace2 checked!  
c           
           temp6 = zero
           temp6 = matmul(fourdamage_inv_t,fourdamage_inv)
           trace2 = temp6(1,1)+temp6(2,2)
c           trace2 = one - damagenorm
c    
           datalamta = zero
c          
          do kewton = 1,newton
c              write(7,*) 'datalamta',datalamta
c
c   eq.(49) in our paper checked!
c              
             dlamta1 = cos(friangle)*(equivalentstress-datalamta*three*
     &               eg*cos(dilangle)*trace1)
c             
             dlamta2 = sqrt(three)*sin(friangle)*(one-damagenorm)*
     &               (sigmac-sigmat)/sigmac/sigmat
c             
             dlamta3 = one/sqrt(three)*ebulk*sin(dilangle)*(sigmac-
     &               sigmat)/sigmac/sigmat*(one-damagenorm)*trace2
c            
             dlamta4 = hydtristress-datalamta*dlamta3
c             
             phy = dlamta1+dlamta2*dlamta4
c
c   load residual judgement  checked!
c             
             eqplas = eqplas_n+(one-damagenorm)*datalamta
             hardstress = yieldo*(one-dkesa1*exp(-dkesa11*eqplas))
c             
c  for plastic hardening stress  checked!
c             
             hp = yieldo*dkesa1*dkesa11*exp(-dkesa11*eqplas)
c
             fyrhs = (dtime/(dtime+alfa0*datalamta))**dm0*phy-
     &             (one-damagenorm)*hardstress
c                
             dlamta11 = -(dtime/(dtime+alfa0*datalamta))**(dm0-one)*
     &                alfa0*dm0*dtime/(dtime+alfa0*datalamta)**two
c   
             dlamta22 = -three*eg*cos(friangle)*cos(dilangle)*trace1-
     &                ebulk*sin(friangle)*sin(dilangle)*
     &                ((sigmac-sigmat)/sigmac/sigmat)**2*
     &                (one-damagenorm)**two*trace2
c             
             dlamta33 = (one-damagenorm)**two*hp            
c             
             dlamta44 = dlamta22*(dtime/(dtime+alfa0*datalamta))**dm0
c             
             datalamta = datalamta-fyrhs/
     &                 (dlamta11*phy+dlamta44-dlamta33) 
c             
             if(abs(fyrhs).lt.toler) goto 10
          end do
c
   10   continue
c        write(7,*) 'fyrhs=',fyrhs
c        
c  update equivalent plastic strain checked!
c    
        eqplas = eqplas_n+(one-damagenorm)*datalamta
c        
c  update equivalentstress
c         
        equivalentstress2 = equivalentstress-datalamta*three*eg*
     &                      cos(dilangle)*trace1
c        
        do k1 = 1,3
             temp7(k1) = zero
           do k2 = 1,3
             temp7(k1) = temp7(k1)+hh_damdam(k1,k2)*dnetstress(k2)/
     &                 equivalentstress
           end do
        end do
c       
        do k1=1,3
           devstress(k1) = devtristress(k1)-datalamta*three*eg*
     &                   cos(dilangle)*temp7(k1)   
        end do   
c        devstress(3) = devtristress(3)
c
c  update hydrostatic stress and deviatoric stress
c  "hh_damdam"=m(-1)*m(-1)*(hh_damage)
c       
        do k1 = 1,3
             plastnorm(k1) = zero
             plastnorm2(k1) = zero
           do k2 = 1,3
             plastnorm(k1) = plastnorm(k1)+hh_damage(k1,k2)*
     &                     dnetstress(k2)/equivalentstress
             plastnorm2(k1) = plastnorm2(k1)+hh_damage(k1,k2)*
     &                     (devstress(k2)-backstress_n(k2))/
     &                      equivalentstress2
           end do
        end do
c
        hydstress = hydtristress-datalamta/sqrt(three)*ebulk*
     &            sin(dilangle)*(sigmac-sigmat)/sigmac/sigmat*
     &            (one-damagenorm)*trace2
c        
c  update stress
c        
        do k1 = 1,2
           stress(k1) = devstress(k1)+hydstress
        end do
        stress(3) = devstress(3)
c        
c  update the elastic strain and plastic strain checked!
c        
        do k1 = 1,2
          eelas(k1) = eelas(k1)-datalamta*
     &                  (three/two*cos(dilangle)*plastnorm2(k1)+
     &                  one/sqrt(three)*sin(dilangle)*
     &                 (one-damagenorm)*(sigmac-sigmat)/sigmac/sigmat)
          eplas(k1) = eplas_n(k1)+datalamta*
     &                  (three/two*cos(dilangle)*plastnorm2(k1)+
     &                  one/sqrt(three)*sin(dilangle)*
     &                 (one-damagenorm)*(sigmac-sigmat)/sigmac/sigmat)
        end do
c         
        eelas(3) = eelas(3)-datalamta*three/two*
     &               cos(dilangle)*plastnorm2(3)
        eplas(3) = eplas_n(3)+datalamta*three/two*
     &               cos(dilangle)*plastnorm2(3)
      
c        
c  strain at the end of new iteration checked!
c      
        strain(1) = eelas(1) + eplas(1)
        strain(2) = eelas(2) + eplas(2)
        strain(3) = eelas(3) + eplas(3)
c        
c  back stress at the end of new iteration suspect
c        
        do k1=1,3
          backstress(k1) = backstress_n(k1)+datalamta*
     &                       hb*cos(dilangle)*plastnorm2(k1)  
        end do
c
        mphi(1) = -one
        mphi(2) = -one
        mphi(3) = zero
c        
c  y=efsita*(dev(m-1/phi)*c*m(-1))*data(efsita) suspect squence    
c  dev(m-1/phi)simplified as [-1,-1,0]
c
        temp8 = zero
        temp8 = matmul(ddsdde,fourdamage_inv)
        tempa1 = zero
c      
        do k1 = 1,3
c           tempa1 = tempa1+mphi(k1)*abs(eelas(k1))
           tempa1 = tempa1+mphi(k1)*abs(eelas(k1))
        end do
c   
        do k1 = 1,3
           temp9(k1) = zero
           do k2 = 1,3
              temp9(k1) = temp9(k1)+temp8(k1,k2)*
     &                  abs(eelas(k2)-eelas_n(k2))
           end do
        end do
c   
c       thermfor_inc means deltaY, thermfor means Y
c        
        do k1 = 1,3
           thermfor_inc(k1) = -two*tempa1*temp9(k1)
        end do
c
        do k1 = 1,3
           temp10(k1) = zero
           do k2 = 1,3
              temp10(k1) = temp10(k1)+temp8(k1,k2)*abs(eelas(k2))    
c              temp10(k1) = temp8(k1,k2)*abs(eelas(k2))
           end do
        end do
c      
        do k1 = 1,3
           thermfor(k1) = -tempa1*temp10(k1)
        end do
c  
        fourorder_j = zero
        do k1 = 1,2
            do k2 = 1,2
              fourorder_j(k1,k2) = one/two           
           end do
        end do
c      
        do k1 = 1,3
           fourorder_j(k1,k1) = one          
        end do
c      
        do k1 = 1,3
              temp11(k1) = zero
           do k2 = 1,3
              temp11(k1) = temp11(k1)+fourorder_j(k1,k2)*
     &                   thermfor(k2)             
           end do
      end do            
c       
        yeq = sqrt(three/two*(thermfor(1)*temp11(1)+
     &                      thermfor(2)*temp11(2)+
     &                      thermfor(3)*temp11(3))) 
c   
        if (yeq.lt.limit) then
             yeq = limit
        end if
c       
c       datalamta_d suspect
c      
        datalamta_d1 = three/two*(temp11(1)*thermfor_inc(1)+
     &               temp11(2)*thermfor_inc(2)+
     &               temp11(3)*thermfor_inc(3))/yeq/dkd
c 
        datalamta_d2 = dkesa2*(sigmac-sigmat)/sigmac/sigmat/
     &               sqrt(three)*(thermfor_inc(1)+thermfor_inc(2)+
     &                thermfor_inc(3))/dkd
!        datalamta_d1 = three/two*temp11(3)*thermfor_inc(3)/yeq/dkd     
!c 
!        datalamta_d2 = dkesa2*(sigmac-sigmat)/sigmac/sigmat/
!     &               sqrt(three)*thermfor_inc(3)/dkd
c
        datalamta_d = datalamta_d1+datalamta_d2
c        
        unitsecondorder = zero
        do k1 = 1,2
             unitsecondorder(k1) = one
        end do
c
        do k1 = 1,3
            damage_inc(k1) = datalamta_d*(three/two*temp11(k1)/yeq
     &                       +dkesa2/sqrt(three)*(sigmac-sigmat)/sigmac
     &                       /sigmat*unitsecondorder(k1))
        end do
        do k1 = 1,2
            damage_inc(k1) = damage_inc(k1)/100.0
        end do
        do k1 = 1,3
            if(damage_inc(k1) .lt. zero) then
                damage_inc(k1) = zero
            end if 
        end do
c        
        do k1 = 1,3
           damage(k1) = damage_n(k1)+damage_inc(k1)
           if(damage(k1) .gt. 0.20) then
               damage(k1) = 0.25
               end if
      end do
      
      !do k1 = 3,3
      !     damage(k1) = damage_n(k1)+damage_inc(k1)
      !     if(damage(k1) .gt. 0.25) then
      !         damage(k1) = 0.25
      !         end if
      !end do
      
c       damage(2) = enu*damage(1)
c
c  formulate the jacobian (material tangent) problem
c   
        unitfourorder = zero      
c
        do k1 = 1,3
             unitfourorder(k1,k1) = one
        end do  
c      
      do k1 = 1,3
           do k2 = 1,3
              fourordermult(k1,k2) = unitsecondorder(k1)*
     &                               unitsecondorder(k2)
           end do
      end do
c   n*n 
        do k1 = 1,3
           do k2 = 1,3
              plastnormmult1(k1,k2) = plastnorm(k1)*plastnorm(k2)
           end do
        end do
c   n*1
        do k1 = 1,3
           do k2 = 1,3
              plastnormmult2(k1,k2) = plastnorm(k1)*unitsecondorder(k2)
           end do
        end do      
c
        do k1 = 1,3
           do k2 = 1,3
              plastnormmult3(k1,k2) = unitsecondorder(k1)*plastnorm(k2)
           end do
        end do  
c
        do k1 = 1,3
           do k2 = 1,3
              devunitfourorder(k1,k2) = unitfourorder(k1,k2)-one/three
           end do
        end do
c        
c   plastic hardening modulus at the end of new iteration checked!
c      
        hp = yieldo*dkesa1*dkesa11*exp(-dkesa11*eqplas)
c
        dlamta444 = zero
        dlamta444 = -dlamta22+(one-damagenorm)**2*
     &          (dtime/(dtime+alfa0*datalamta))**(-dm0)*hp+
     &           alfa0*dm0/(dtime+alfa0*datalamta)*phy
c
        do k1 = 1,3
          do k2 = 1,3
           temp12(k1,k2) = 2*eg*(unitfourorder(k1,k2)-datalamta*
     &       three*eg*cos(dilangle)*hh_damdam(k1,k2)
     &       /equivalentstress)*devunitfourorder(k2,k1)    
c     
           temp13(k1,k2) = nine*eg*eg*cos(dilangle)*datalamta*
     &         temp6(k1,k2)/equivalentstress*plastnormmult1(k2,k1)    
c                
           temp14(k1,k2) = nine*eg*eg*cos(dilangle)*cos(friangle)*
     &         temp6(k1,k2)*plastnormmult1(k2,k1)/dlamta444
c           
           temp15(k1,k2) = three*sqrt(three)*eg*ebulk*cos(dilangle)*
     &         sin(friangle)*(one-damagenorm)*(sigmac-sigmat)/sigmac/
     &         sigmat*temp6(k1,k2)/dlamta444*plastnormmult2(k2,k1)
c
           temp16(k1,k2) = sqrt(three)*eg*ebulk*sin(dilangle)*
     &         cos(friangle)*(one-damagenorm)*(sigmac-sigmat)/sigmac/
     &         sigmat*trace2/dlamta444*plastnormmult3(k1,k2)
c
           temp17(k1,k2) = ebulk*fourordermult(k1,k2)
c   
           temp18(k1,k2) = ebulk*ebulk*sin(dilangle)*sin(friangle)
     &         *((sigmac-sigmat)/sigmac/sigmat)**2*(one-damagenorm)**2*
     &         trace2/dlamta444*fourordermult(k1,k2)
c     
           ddsdde(k1,k2) = temp12(k1,k2)+temp13(k1,k2)-temp14(k1,k2)-
     &                    temp15(k1,k2)-temp16(k1,k2)+temp17(k1,k2)-  
     &                    temp18(k1,k2)

          end do
       end do       
       !ddsdde(3,1) = -62.0
       !ddsdde(3,2) = -62.0
       !ddsdde(1,3) = -62.0
       !ddsdde(2,3) = -62.0
       !ddsdde(1,2) = ddsdde(1,2)-60.0
       !ddsdde(2,1) = ddsdde(2,1)-60.0
       
       end if
c        
        statev(1) = eelas(1)
        statev(2) = eelas(2)
        statev(3) = eelas(3)
        statev(4) = eplas(1)
        statev(5) = eplas(2)
        statev(6) = eplas(3)
        statev(7) = stress(1)
        statev(8) = stress(2)
        statev(9) = stress(3)
        statev(10) = backstress(1)
        statev(11) = backstress(2)
        statev(12) = backstress(3)
        statev(13) = eqplas
        statev(14) = damage(1)
        statev(15) = damage(2)
        statev(16) = damage(3)
c        
c       statev below only used for drawing the plot
c        
        statev(17) = strain(1)
        statev(18) = strain(2)
        statev(19) = strain(3)
        statev(20) = hydstress
c        
c       key to make/cancel explanatory note ctrl+k+c/ctrl+k+u        
c       symmetre matrice check
c        
      if(noel.eq.4261) then
         !write(7,*) 'noel',noel
         !write(7,*) 'ddsdde',ddsdde
         !write(7,*) 'stress',stress
         !write(7,*) 'devstress',devstress
         ! write(7,*) 'temp12',temp12
         !write(7,*) 'temp13',temp13
         !write(7,*) 'temp14',temp14
         !write(7,*) 'temp15',temp15
         !write(7,*) 'temp16',temp16
         !write(7,*) 'temp17',temp17
         !write(7,*) 'temp18',temp18
      end if  
!c      
!      if((kstep.eq.1).and.(kinc.eq.250)) then   
!         write(7,*) 'temp18',temp18
!         write(7,*) 'hydstress',hydstress
!         write(7,*) 'datalamta_d1',datalamta_d1
!         write(7,*) 'datalamta_d2',datalamta_d2
!         write(7,*) 'yeq',yeq
!         write(7,*) 'thermfor',thermfor
!      end if  
!c      
      !if((kstep.eq.1).and.(kinc.eq.5000)) then   
      !   write(7,*) 'eqplas_n',eqplas_n
      !end if  
      !if((ddsdde(3,3).lt.zero)) then   
         !write(7,*) 'temp6(3,3)',temp6(3,3)
         !write(7,*) 'dlamta444',dlamta444
         !write(7,*) 'plastnormmult1(3,3)',plastnormmult1(3,3)
         !write(7,*) 'temp13',temp13
         !write(7,*) 'temp14',temp14
         !write(7,*) 'temp15',temp15
         !write(7,*) 'temp16',temp16
         !write(7,*) 'temp17',temp17
         !write(7,*) 'temp18',temp18
         !write(7,*) 'hydstress',hydstress
         !write(7,*) 'datalamta_d1',datalamta_d1
         !write(7,*) 'datalamta_d2',datalamta_d2
         !write(7,*) 'yeq',yeq
         !write(7,*) 'thermfor',thermfor
      !end if  
c      
c        write(7,*) 'datalamta_d',datalamta_d
        !write(7,*) 'aa',aa
        !write(7,*) 'hh_damage',hh_damage
        !write(7,*) 'dnetstress',dnetstress         
        !write(7,*) 'hh_damdam',hh_damdam
c        
c       plasiticity check
c
        if(yieldcondition.gt.zero) then
        !write(7,*) 'hardstress',hardstress
        !write(7,*) 'yieldcondition',yieldcondition
        !write(7,*) 'yieldcondition1',yieldcondition1
        !write(7,*) 'yieldcondition2',yieldcondition2
        !write(7,*) 'eqplas',eqplas
        !write(7,*) 'damagenorm',damagenorm
        !write(7,*) 'yeq',yeq
        !write(7,*) 'ddsdde',ddsdde
        !write(7,*) 'dlamta444',dlamta444
        !write(7,*) 'thermfor',thermfor
        !write(7,*) 'thermfor_inc',thermfor_inc     
        !write(7,*) 'damage',damage
        !write(7,*) 'damage_inc',damage_inc    
        !write(7,*) 'temp12',temp12
        !write(7,*) 'temp13',temp13
        !write(7,*) 'temp14',temp14
        !write(7,*) 'temp15',temp15
        !write(7,*) 'temp16',temp16
        !write(7,*) 'temp17',temp17
        !write(7,*) 'temp18',temp18
        !write(7,*) 'temp6(3,3)',temp6(3,3)
        !write(7,*) 'dlamta444',dlamta444
        !write(7,*) 'plastnormmult1(3,3)',plastnormmult1(3,3)
        !write(7,*) 'dnetstress(3)',dnetstress(3)
        !write(7,*) 'equivalentstress',equivalentstress
        end if
c 

       return
       
       end
c
c**********************************************************************
	     subroutine m2inv(a,c,n)
c  get inverse matrix
      ! a(n,n) - array of coefficients for matrix a
      ! n      - dimension
      ! c(n,n) - inverse matrix of a
c**********************************************************************
       integer n,i, j, k  
       real*8 a(n,n),c(n,n),al(n,n),au(n,n),b(n),d(n),x(n),tem_a(n,n)
       real*8 coeff
       parameter (zero=0.d0,one=1.d0)

       do i=1,n
          do j=1,n
            tem_a(i,j)=a(i,j)
          end do
       end do
       
       coeff=zero
       
       do k=1, n-1
          do i=k+1, n
              coeff=tem_a(i,k)/tem_a(k,k)
              al(i,k) = coeff
              do j=k+1,n
                  tem_a(i,j) = tem_a(i,j)-coeff*tem_a(k,j)         
              end do
           end do
       end do
      
       do i=1,n
          al(i,i) = one
       end do
      
       do j=1,n
          do i=1,j
              au(i,j) = tem_a(i,j)
          end do
       end do
      
       do k=1,n
          b(k) = one
          d(1) = b(1)
          do i=2,n
              d(i)=b(i)
              do j=1,i-1
                 d(i) = d(i) - al(i,j)*d(j)
              end do
          end do
          
          x(n)=d(n)/au(n,n)
          do i = n-1,1,-1
              x(i) = d(i)
              do j=n,i+1,-1
                  x(i)=x(i)-au(i,j)*x(j)
              end do
            x(i) = x(i)/au(i,i)
          end do
          
          do i=1,n
              c(i,k) = x(i)
          end do
            b(k)=zero
        end do    
      
      return
      end
c**********************************************************************
	    subroutine mdet(a,det,m)
c 	this subroutine calculates the determinant
c 	of a m by m matrix [a].
c**********************************************************************
    	 real*8 a(m,m),det
    	 integer m
  
       det=0.d0
       if (m.eq.3) then
	        det =	a(1,1)*a(2,2)*a(3,3)
     &	        + a(1,2)*a(2,3)*a(3,1)
     &	        + a(1,3)*a(2,1)*a(3,2)
     &		      - a(3,1)*a(2,2)*a(1,3)
     &		      - a(3,2)*a(2,3)*a(1,1)
     &		      - a(3,3)*a(2,1)*a(1,2)
     
       else if (m.eq.2) then
	        det = a(1,1)*a(2,2)-a(2,1)*a(1,2)

       end if
c
	    return
	    end
c**********************************************************************
c	the following subroutines calculate the spectral
c	decomposition of a symmetric three by three matrix
c**********************************************************************
	    subroutine spectral(a,d,v,m)
c	this subroutine calculates the eigenvalues and eigenvectors of
c	a symmetric m by m matrix [a].
c
c	the output consists of a vector d containing the three
c	eigenvalues in ascending order, and
c	a matrix [v] whose columns contain the correspondin	eigenvectors.
c**********************************************************************
	    real*8 d(m),v(m,m),a(m,m),e(m,m)
	    integer m,i,j,np,nrot
	    
	    np=m
  	   do i = 1,m
          do j= 1,m
            e(i,j) = a(i,j)
          end do
       end do   

	     call jacobi(e,m,np,d,v,nrot)
	     call eigsrt(d,v,m,np)

	    return
	    end
c
c**********************************************************************
	    subroutine jacobi(a,n,np,d,v,nrot)
c	computes all eigenvalues and eigenvectors of a symmetric
c	matrix [a], which is of size n by n, stored in a physical
c	np by bp array. on output, elements of [a] above the diagonal
c	are destroyed, but the diagonal and sub-diagonal are unchanged
c	and give full information about the original symmetric matrix.
c	vector d returns the eigenvalues of [a] in its first n elements.
c	[v] is a matrix with the same logical and physical dimensions as
c	[a] whose columns contain, on output, the normalized
c	eigenvectorsof [a]. nrot returns the number of jacobi rotations
c	which were required.

c	this subroutine is taken from "numerical recipes", page 346.
c**********************************************************************
      real*8 a(np,np),d(np),v(np,np),b(100),z(100)
      real*8 sm,tresh,g,t,c,s,tau,h,theta
      integer ip,n,nrot,i,iq,j
 
      sm=0.d0
      tresh=0.d0
      g=0.d0
      t=0.d0
      c=0.d0
      s=0.d0
      tau=0.d0
      h=0.d0
      theta=0.d0
      
c	initialize [v] to the identity matrix

	    do ip = 1,n
	      do iq = 1,n
	           v(ip,iq) = 0.d0
        end do
             v(ip,ip) = 1.d0
      end do

c	initialize [b] and [d] to the diagonal of [a], and z to zero.
c	the vector z will accumulate terms of the form t*a_pq as
c	in equation (11.1.14)

    	do ip = 1,n
	      b(ip) = a(ip,ip)
	      d(ip) = b(ip)
	      z(ip) = 0.d0
      end do
c
      	nrot = 0
	      do i = 1,50

c	sum off-diagonal elements

          do ip = 1, n-1
             do  iq = ip + 1, n
	              sm = sm + abs(a(ip,iq))
             end do
          end do 

c	if sum = 0., then return. this is the normal return
c	which relies on quadratic convergence to machine
c	underflow.

        if (sm .eq. 0.d0) return
c
c	  if (sm .lt. 1.0d-15) return

c	in the first three sweeps carry out the pq rotation only if
c	|a_pq| > tresh, where tresh is some threshold value,
c	see equation (11.1.25). thereafter tresh = 0.

          if ( i .lt. 4) then
            tresh = 0.2*sm/n**2
          else
            tresh = 0.d0
          end if
c
          do ip = 1, n-1
            do iq = ip+1,n
              g = 100.d0*abs(a(ip,iq))

c	after four sweeps, skip the rotation if the
c	off-diagonal element is small.

	      if ((i .gt. 4) .and. (abs(d(ip))+g .eq. abs(d(ip)))
     &       .and. (abs(d(iq))+g .eq. abs(d(iq)))) then
                a(ip,iq) = 0.d0
        else if ( abs(a(ip,iq)) .gt. tresh) then
                h = d(iq) - d(ip)
          if ((abs(h)+g) .eq. abs(h)) then

c	t = 1./(2.*theta), equation(11.1.10)

	          t =a(ip,iq)/h
	        else
	          theta = 0.5d0*h/a(ip,iq)
	          t =1.0/(abs(theta)+sqrt(1.d0+theta**2))
	          if (theta .lt. 0.d0) t = -t
	        end if
	        
	        c = 1.0/sqrt(1.d0 + t**2)
	        s = t*c
	        tau = s/(1.d0 + c)
	        h = t*a(ip,iq)
	        z(ip) = z(ip) - h
	        z(iq) = z(iq) + h
	        d(ip) = d(ip) - h
	        d(iq) = d(iq) + h
	        a(ip,iq) = 0.d0

c	case of rotations 1 <= j < p

	        do j = 1, ip-1
	          g = a(j,ip)
	          h = a(j,iq)
	          a(j,ip) = g - s*(h + g*tau)
	          a(j,iq) = h + s*(g - h*tau)
          end do

c	case of rotations p < j < q

	        do j = ip+1, iq-1
	          g = a(ip,j)
	          h = a(j,iq)
	          a(ip,j) = g - s*(h + g*tau)
	          a(j,iq) = h + s*(g - h*tau)
          end do

c	case of rotations q < j <= n

	        do j = iq+1, n
            g = a(ip,j)
	          h = a(iq,j)
	          a(ip,j) = g - s*(h + g*tau)
	          a(iq,j) = h + s*(g - h*tau)
          end do
 
	        do j = 1,n
	          g = v(j,ip)
	          h = v(j,iq)
	          v(j,ip) = g - s*(h + g*tau)
	          v(j,iq) = h + s*(g - h*tau)
          end do

	        nrot = nrot + 1
       end if
 
         end do
        end do

c	update d with the sum of t*a_pq, and reinitialize z

	      do ip = 1, n
	       b(ip) = b(ip) + z(ip)
	       d(ip) = b(ip)
	       z(ip) = 0.d0
        end do

      end do

c	if the algorithm has reached this stage, then
c	there are too many sweeps, print a diagnostic
c	and stop.

c	write (*,'(/1x,a/)') '50 iters in jacobi should never happen'

  	  return
	    end

c**********************************************************************
	    subroutine eigsrt(d,v,n,np)
c	given the eigenvalues [d] and eigenvectors [v] as output from
c	jacobi, this routine sorts the eigenvalues into ascending order,
c	and rearranges the columns of [v] accordingly.

c	this subroutine is taken from "numerical recipes", p. 348.
c**********************************************************************
	    real*8 d(np),v(np,np),p
      integer i,n,k,np,j
      
	    do i = 1,n-1
	      k = i
	      p = d(i)
	   
	      do j = i+1,n
	        if (d(j) .ge. p) then
	         k = j
	         p = d(j)
	        end if
        end do

	      if (k .ne. i) then
	        d(k) = d(i)
	        d(i) = p
	          do j = 1,n
	            p = v(j,i)
	            v(j,i) = v(j,k)
	            v(j,k) = p
            end do
         end if

       end do
    
	    return
	    end