*Crystal plasticity subroutine
*The code is distributed under a BSD license     
      
*If using this code for research or industrial purposes, please cite:
* Tan W and Liu B (2020).
* "A physically-based constitutive model for the shear-dominated response
* and strain rate effect of carbon fibre reinforced composites.",
* Elsevier  Composites Part B: Engineering,193 (2020), 108032

      
* Wei Tan(wei.tan@qmul.ac.uk)
* Queen Mary University London
************************************
      
      
      subroutine vumat(
C Read only (unmodifiable)variables -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2  stepTime, totalTime, dt, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     5  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only (modifiable) variables -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
C
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     1  charLength(nblock), strainInc(nblock,ndir+nshr),
     2  relSpinInc(nblock,nshr), tempOld(nblock),
     3  stretchOld(nblock,ndir+nshr),
     4  defgradOld(nblock,ndir+nshr+nshr),
     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(nblock),
     8  stretchNew(nblock,ndir+nshr),
     8  defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv),
     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock)
C
      character*80 cmname
C
C     define state variables 
C     state(*,1) = fp11
C     state(*,2) = fp22
C     state(*,3) = fp33
C     state(*,4) = fp12
C     state(*,5) = fp23
C     state(*,6) = fp31
C     state(*,7) = fp21
C     state(*,8) = fp32
C     state(*,9) = fp13 
CCCCCCCCCCCCCCCCCCCCCCC
C     state(*,10) = sig11
C     state(*,11) = sig22
C     state(*,12) = sig33
c     state(*,13) = sig12
C     state(*,14) = sig23
C     state(*,15) = sig31
CCCCCCCCCCCCCCCCCCCCCCCCCCC
C     ij, jth component of system i
C     state(*,16) = s11 ij, jth component of system i
C     state(*,17) = s12
C     state(*,18) = s13
C     state(*,19) = m11
C     state(*,20) = m12
C     state(*,21) = m13 
CCCCCCCCCCCCCCCCCCCCCCCCCCC
C     state(*,22) = s21
C     state(*,23) = s22
C     state(*,24) = s23 
C     state(*,25) = m21
C     state(*,26) = m22
C     state(*,27) = m23
CCCCCCCCCCCCCCCCCCCCCCCCC
C     state(*,28) = s31
C     state(*,29) = s32
C     state(*,30) = s33
C     state(*,31) = m31
C     state(*,32) = m32
C     state(*,33) = m33
CCCCCCCCCCCCCCCCCCCCCCCCC
C     state(*,34) = s41
C     state(*,35) = s42
C     state(*,36) = s43
C     state(*,37) = m41
C     state(*,38) = m42
C     state(*,39) = m43
CCCCCCCCCCCCCCCCCCCCCCCC
C     state(*,40) = s51
C     state(*,41) = s52
C     state(*,42) = s53
C     state(*,43) = m51
C     state(*,44) = m52
C     state(*,45) = m53
c     state(*,46) = psig2
C     state(*,47) = psig3 
C     state(*,48) = psig1
C     state(*,49) = plastic dissipation 
C     state(*,50) = so6
C     state(*,53) = mo6

C=========================
C      so11 = props(10)
C      so12 = props(11)
C      so13 = props(12)
C      xmo11 = props(13)
C      xmo12 = props(14)
C      xmo13 = props(15)
C      so21 = props(16)
C      so22 = props(17)
C      so23 = props(18)
C      xmo21 = props(19)
C      xmo22 = props(20)
C      xmo23 = props(21)
C      so31 = props(22)
C      so32 = props(23)
C      so33 = props(24) 
C      xmo31 = props(25)
C      xmo32 = props(26)
C      xmo33 = props(27)
C      so41 = props(28)
C      so42 = props(29)
C      so43 = props(30)
C      xmo41 = props(31)
C      xmo42 = props(32)
C      xmo43 = props(33)
C      so51 = props(34)
C      so52 = props(35)
C      so53 = props(36)
C      xmo51 = props(37)
C      xmo52 = props(38)
C      xmo53 = props(39)
C      so6   = props(47)
C      mo6   = props(50)
C==========================
C     hardening 

C     state(*,56) = gamma(1) (Culumated plastic shear strain system1)
C     state(*,57) = gamma(2)
C     state(*,58) = gamma(3)
C     state(*,59) = gamma(4)
C     state(*,60) = gamma(5)
C     state(*,61) = gamma(6)

CC==================
CC update the strain rate in transverse and shear modulus and shear strain

C     Declare tensor 


      dimension sig(3,3), fp(3,3),dfp(3,3),f(3,3),
     *     fpi(3,3), fs(3,3),e(3,3),xiden(3,3), pk2(3,3),
     *     fsi(3,3), u(3,3),ui(3,3),r(3,3),ep(3,3),
     *     csig(3,3),sm(6,6),vs(6),ve(6),so1(3),so2(3),so3(3),
     *     so4(3),so5(3),xmo1(3),xmo2(3),xmo3(3),xmo4(3),xmo5(3),
     *     s1(3),s2(3),s3(3),s4(3),s5(3),xm1(3),xm2(3),xm3(3),gamma(6),
     *     xm4(3),xm5(3),tao(6),taoh(6),dgam(6),ddfp(3,3),so6(3),xmo6(3),
     *     s6(3),xm6(3)
      parameter(zero=0d0,one=1d0,half=0.5d0,two=2d0,six=6d0)
      
      if (nstatev.lt.55) then
           write(6,*) "not enough state vars"
           stop
      endif
               
      xiden(1,1)=one
      xiden(2,1)=zero
      xiden(3,1)=zero 
      xiden(1,2)=zero
      xiden(2,2)=one
      xiden(3,2)=zero
      xiden(1,3)=zero
      xiden(2,3)=zero
      xiden(3,3)=one       
            
C     define properties 
C     xMu=0.05, dgam0=0.001, taoy=2.0, xn=10.0)   
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      ey1 = props(1)
      ey2 = props(2)
      ey3 = props(46)
      emu12 = props(3)
      emu13 = props(47)
      emu23 = props(4)
      xg12 = props(5)      
      xg31 = props(48)
      xg23 = props(49)
      xmu = props(6)
      dgam0 = props(7)
      taoy = props(8)
      xn = props(9) 
      pfactor0 = props(50)
      pfactor90 = props(51)
      den=props(52)
      hfactor=props(53)      
      eRate=props(54)
      xn2=props(55)
      
      do i = 1,3 
        so1(i) = props(i+9)
        xmo1(i) = props(i+12)
        so2(i) = props(i+15)
        xmo2(i) = props(i+18)
        so3(i) = props(i+21)
        xmo3(i) = props(i+24)
        so4(i) = props(i+27)
        xmo4(i) = props(i+30)
        so5(i) = props(i+33)
        xmo5(i) = props(i+36)
        so6(i) = props(i+39)
        xmo6(i) = props(i+42)        
      enddo 
 
      
C============= Initiation of G12 G13 and G23 =========C 
      xg120 = props(5)
      xg310 = props(48)
      xg230 = props(49)
      ey20  = props(2)  
      ey30  = props(46)        
      refrate=0.0001
  
      taoy0 = props(8)      
      
C=================================
      emu21 = emu12*ey2/ey1
      emu31 = emu13*ey3/ey1
      emu32 = emu23*ey3/ey2


C     derive stiffness matrix 
      
      delta = one-emu12*emu21-emu23*emu32-emu31*emu13
     .        -two*emu12*emu23*emu31
      
      delta = delta/(ey1*ey2*ey3)

      
C     corresponding strain order 11. 22. 33. 12. 23 .31

      sm(1,1) = (one-emu23*emu32)/(ey2*ey3*delta)
      sm(1,2) = (emu21+emu31*emu23)/(ey2*ey3*delta)
      sm(1,3) = (emu31+emu21*emu32)/(ey2*ey3*delta)
      sm(1,4) = zero 
      sm(1,5) = zero
      sm(1,6) = zero
      sm(2,2) = (one-emu31*emu13)/(ey3*ey1*delta)
      sm(2,3) = (emu32+emu31*emu12)/(ey3*ey1*delta)
      sm(2,4) = zero
      sm(2,5) = zero 
      sm(2,6) = zero 
      sm(3,3) = (one-emu12*emu21)/(ey1*ey2*delta)
      sm(3,4) = zero
      sm(3,5) = zero 
      sm(3,6) = zero 
      sm(4,4) = two*xg12
      sm(4,5) = zero 
      sm(4,6) = zero 
      sm(5,5) = two*xg23
      sm(5,6) = zero 
      sm(6,6) = two*xg31
      
      
      sm(2,1) = sm(1,2)
      
      sm(3,1) = sm(1,3)
      sm(3,2) = sm(2,3)
      
      sm(4,1) = sm(1,4)
      sm(4,2) = sm(2,4)
      sm(4,2) = sm(2,4)
      sm(4,3) = sm(3,4) 
      
      sm(5,1) = sm(1,5)
      sm(5,2) = sm(2,5)
      sm(5,3) = sm(3,5)
      sm(5,4) = sm(4,5)
       
      sm(6,1) = sm(1,6)
      sm(6,2) = sm(2,6) 
      sm(6,3) = sm(3,6) 
      sm(6,4) = sm(4,6) 
      sm(6,5) = sm(5,6)  
      
C      print *, 'sm=', sm
         
         
      If (totalTime .eq. zero) then
          
      do km = 1,nblock          

C     define deformation gradient

      f(1,1) = defgradNew(km,1) 
      f(2,2) = defgradNew(km,2) 
      f(3,3) = defgradNew(km,3)
      f(1,2) = defgradNew(km,4)
      f(2,3) = defgradNew(km,5)
      f(3,1) = defgradNew(km,6)
      f(2,1) = defgradNew(km,7)
      f(3,2) = defgradNew(km,8)
      f(1,3) = defgradNew(km,9)


      
      
C     calculate Green strain 
      
C     calculate e 
      
      do i = 1,3
          do j = 1,3
              e(i,j) = zero
          enddo
      enddo 
      
      do i = 1,3
          do j = 1,3
              do k = 1,3
                 e(i,j) = e(i,j)+half*(f(k,i)*f(k,j))
              enddo
          enddo
      enddo
      
      do i=1,3
        do j=1,3
           e(i,j)=e(i,j)-half*xiden(i,j)
        enddo
      enddo
      
C     Caluclate Pk2 stress 
CCCCCC     initialize PK2 stress 
CCCCCC

CCCCCC
cCCCCCvecterize e

      do i = 1,6
         vs(i) = zero
      enddo

      ve(1) = e(1,1)
      ve(2) = e(2,2)
      ve(3) = e(3,3)
      ve(4) = e(1,2)
      ve(5) = e(2,3)
      ve(6) = e(3,1) 
      
      do i = 1,6 
         do j = 1,6
            vs(i) = vs(i)+sm(i,j)*ve(j)
         enddo
      enddo
     
CCCCCC store vs 
      
      pk2(1,1) = vs(1)
      pk2(2,2) = vs(2)
      pk2(3,3) = vs(3)
      pk2(1,2) = vs(4)
      pk2(2,3) = vs(5)
      pk2(3,1) = vs(6) 
      pk2(1,3) = pk2(3,1) 
      pk2(2,1) = pk2(1,2) 
      pk2(3,2) = pk2(2,3) 

      
C     rotate to update Cauchy stress 

C     initialize sig, clear the Old value
      do i = 1,3
          do j = 1,3
              sig(i,j) = zero
          enddo
      enddo

C     assign new values 
      do i = 1,3
          do j = 1,3
              do k = 1,3
                  do l = 1,3 
                     sig(i,j) = sig(i,j)+f(i,k)*pk2(k,l)*f(j,l)
                  enddo                  
              enddo
          enddo
      enddo
      
      detf = f(1,1)*(f(2,2)*f(3,3)-f(3,2)*f(2,3))
      detf = detf - f(1,2)*(f(2,1)*f(3,3)-f(2,3)*f(3,1))
      detf = detf + f(1,3)*(f(2,1)*f(3,2)-f(2,2)*f(3,1))
      
c     divide by J             
      do i = 1,3
         do j = 1,3
            sig(i,j) = sig(i,j)/detf
         enddo 
      enddo 
          
C     Inverse u   

      u(1,1) = stretchNew(km,1)
      u(2,2) = stretchNew(km,2) 
      u(3,3) = stretchNew(km,3)
      u(1,2) = stretchNew(km,4)
      u(2,3) = stretchNew(km,5)
      u(3,1) = stretchNew(km,6)
      u(2,1) = u(1,2)
      u(1,3) = u(3,1)
      u(3,2) = u(2,3)
      
      detu = u(1,1)*(u(2,2)*u(3,3)-u(3,2)*u(2,3))
      detu = detu-u(1,2)*(u(2,1)*u(3,3)-u(3,1)*u(2,3))
      detu = detu+u(1,3)*(u(2,1)*u(3,2)-u(3,1)*u(2,2))
      
      ui(1,1) = (u(2,2)*u(3,3)-u(2,3)*u(3,2))/detu
      ui(1,2) = -(u(1,2)*u(3,3)-u(1,3)*u(3,2))/detu
      ui(1,3) = (u(1,2)*u(2,3)-u(1,3)*u(2,2))/detu
      ui(2,1) = -(u(2,1)*u(3,3)-u(2,3)*u(3,1))/detu
      ui(2,2) = (u(1,1)*u(3,3)-u(1,3)*u(3,1))/detu
      ui(2,3) = -(u(1,1)*u(2,3)-u(1,3)*u(2,1))/detu
      ui(3,1) = (u(2,1)*u(3,2)-u(2,2)*u(3,1))/detu
      ui(3,2) = -(u(1,1)*u(3,2)-u(1,2)*u(3,1))/detu
      ui(3,3) = (u(1,1)*u(2,2)-u(1,2)*u(2,1))/detu
      
C     calculate r 
      do i = 1,3
         do j = 1,3
            r(i,j) = zero
         enddo
      enddo
      
      do i = 1,3
         do j = 1,3
            do k = 1,3
               r(i,j) = r(i,j)+f(i,k)*ui(k,j)
            enddo
         enddo
      enddo     
      
C     cauchy stress to corotational system      

      do i = 1,3
         do j = 1,3
            csig(i,j) = zero 
         enddo
      enddo
      
      do i = 1,3
         do j = 1,3
            do k = 1,3
               do l = 1,3
                  csig(i,j) = csig(i,j)+r(k,i)*sig(k,l)*r(l,j)
               enddo
            enddo
         enddo
      enddo

C     Get rid of the rotation problem 

      do i = 1,6 
         do j = 1,6
            vs(i) = vs(i)+sm(i,j)*strainInc(km,j)
         enddo
      enddo
      
C     Use the PK2/material stress to initialize 
C     the corotational stress
      csig(1,1) = vs(1)
      csig(2,2) = vs(2)
      csig(3,3) = vs(3)
      csig(1,2) = vs(4)
      csig(2,3) = vs(5)
      csig(1,3) = vs(6)
      csig(2,1) = csig(1,2)
      csig(3,2) = csig(2,3)
      csig(3,1) = csig(1,3)      
      
C     Update StressNew and StateNew for Stress 
    
      stressNew(km,1) = csig(1,1)
      stressNew(km,2) = csig(2,2)
      stressNew(km,3) = csig(3,3) 
      stressNew(km,4) = csig(1,2)
      stressNew(km,5) = csig(2,3)
      stressNew(km,6) = csig(3,1) 
      


      
      stateNew(km,10) = sig(1,1)
      stateNew(km,11) = sig(2,2)
      stateNew(km,12) = sig(3,3) 
      stateNew(km,13) = sig(1,2)
      stateNew(km,14) = sig(2,3)
      stateNew(km,15) = sig(3,1) 
      
      enddo
      
      
C================================            
      else 

      do 100 km = 1,nblock
C     define plastic deformation gradient tensor 

      fp(1,1) = stateOld(km,1)
      fp(2,2) = stateOld(km,2)
      fp(3,3) = stateOld(km,3)
      fp(1,2) = stateOld(km,4)
      fp(2,3) = stateOld(km,5)
      fp(3,1) = stateOld(km,6)
      fp(2,1) = stateOld(km,7)
      fp(3,2) = stateOld(km,8)
      fp(1,3) = stateOld(km,9)

C     derive stiffness matrix 
      
C    strain rate effect of shear modulus
C       do i=1,6
C          sr(i)=strainInc(km,i)/dt
C       enddo
C      print *, 'strain rate', strainrate 
C      print *, 'Totaltime=', Totaltime  
C      print *, 'sr=', sr 

C     define culumated plastic shear strain gamma 
      
      gamma(1) = stateOld(km,56)
      gamma(2) = stateOld(km,57)
      gamma(3) = stateOld(km,58)
      gamma(4) = stateOld(km,59)
      gamma(5) = stateOld(km,60)
      gamma(6) = stateOld(km,61)

      xg12=xg120*(xn2*log10(eRate/refrate)+1)      
      xg23=xg310*(xn2*log10(eRate/refrate)+1)      
      xg31=xg230*(xn2*log10(eRate/refrate)+1)
      ey2=ey20*(xn2*log10(eRate/refrate)+1)      
      ey3=ey30*(xn2*log10(eRate/refrate)+1) 
      
      taoy=taoy0*(xn*log10(eRate/refrate)+1)    
C       print *,'xg12','=',xg12    
C       print *,'taoy','=',taoy 
      
C     corresponding strain order 11. 22. 33. 12. 23 .31

      sm(1,1) = (one-emu23*emu32)/(ey2*ey3*delta)
      sm(1,2) = (emu21+emu31*emu23)/(ey2*ey3*delta)
      sm(1,3) = (emu31+emu21*emu32)/(ey2*ey3*delta)
      sm(1,4) = zero 
      sm(1,5) = zero
      sm(1,6) = zero
      sm(2,2) = (one-emu31*emu13)/(ey3*ey1*delta)
      sm(2,3) = (emu32+emu31*emu12)/(ey3*ey1*delta)
      sm(2,4) = zero
      sm(2,5) = zero 
      sm(2,6) = zero 
      sm(3,3) = (one-emu12*emu21)/(ey1*ey2*delta)
      sm(3,4) = zero
      sm(3,5) = zero 
      sm(3,6) = zero 
      sm(4,4) = two*xg12
      sm(4,5) = zero 
      sm(4,6) = zero 
      sm(5,5) = two*xg23
      sm(5,6) = zero 
      sm(6,6) = two*xg31
      
      
      sm(2,1) = sm(1,2)
      
      sm(3,1) = sm(1,3)
      sm(3,2) = sm(2,3)
      
      sm(4,1) = sm(1,4)
      sm(4,2) = sm(2,4)
      sm(4,2) = sm(2,4)
      sm(4,3) = sm(3,4) 
      
      sm(5,1) = sm(1,5)
      sm(5,2) = sm(2,5)
      sm(5,3) = sm(3,5)
      sm(5,4) = sm(4,5)
       
      sm(6,1) = sm(1,6)
      sm(6,2) = sm(2,6) 
      sm(6,3) = sm(3,6) 
      sm(6,4) = sm(4,6) 
      sm(6,5) = sm(5,6)  
     
      
C     define slipDirection s and slipNormal m
  
      do i = 1,3
          s1(i)  = stateOld(km,i+15)
          xm1(i) = stateOld(km,i+18)
          s2(i)  = stateOld(km,i+21)
          xm2(i) = stateOld(km,i+24)
          s3(i)  = stateOld(km,i+27)
          xm3(i) = stateOld(km,i+30)    
          s4(i)  = stateOld(km,i+33)
          xm4(i) = stateOld(km,i+36)
          s5(i)  = stateOld(km,i+39)
          xm5(i) = stateOld(km,i+42)     
          s6(i)  = stateOld(km,i+49)
          xm6(i) = stateOld(km,i+52)         
      end do 
    
C     calculate shear stress 
C     initialize tao (clear the effect of previous step) 
C     tao(i),p(i) corresponding to ith slip system 

      do i = 1,6
         tao(i) = zero
      enddo 
         
C     define stress tensor 
          
      sig(1,1) = stateOld(km,10)
      sig(2,2) = stateOld(km,11)
      sig(3,3) = stateOld(km,12)
      sig(1,2) = stateOld(km,13)
      sig(2,3) = stateOld(km,14)
      sig(3,1) = stateOld(km,15)    
      
      sig(2,1) = sig(1,2)
      sig(1,3) = sig(3,1)
      sig(3,2) = sig(2,3)    
      
      psig2 = stateOld(km,46)
      psig3 = stateOld(km,47) 
      psig1 = stateOld(km,48)

      do i = 1,3
          do j = 1,3
              tao(1) = tao(1) + s1(i)*sig(i,j)*xm1(j)
              tao(2) = tao(2) + s2(i)*sig(i,j)*xm2(j)
              tao(3) = tao(3) + s3(i)*sig(i,j)*xm3(j)
              tao(4) = tao(4) + s4(i)*sig(i,j)*xm4(j)
              tao(5) = tao(5) + s5(i)*sig(i,j)*xm5(j)
              tao(6) = tao(6) + s6(i)*sig(i,j)*xm6(j)
          end do
      end do   

C      print *, 'psig2=', psig2
C      print *, 'psig3=', psig3


C     calcualte pressure 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      p =-half*pfactor0*(psig2+psig3)-half*pfactor90*(psig1+psig3)

C      print *, 'p=', p
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
C     define hardening and softening

      do i = 1,6 
        if (p .GE. zero) then 
             taoh(i) = max(abs(tao(i))-xmu*p,zero)
         else 
            taoh(i) = abs(tao(i))
        endif 
      enddo    

C      print *, 'taoh=', taoh      
      
C===== define slip hardening ==========================
      
      gammatol = zero 
	  
      do i = 1,6
	     gammatol = gammatol + gamma(i) 
      enddo
	  
      gammatol = one/six*gammatol
      taoyhard    = taoy + hfactor*gammatol
C      taoy0=taoyhard
C      print *, 'taoyhard=', taoyhard      
     
C===== finish =========================================
    
      
C     calculate slip rate dgam              
C     reference slip rate dgam0, strength taoy, power xn
C      dgam0=dgam0/charLength(km)
      do i = 1,6         
         dgam(i) = dgam0*sign(one,tao(i))*(taoh(i)/taoyhard)**(2.0)
C      dgam(i)=dgam0*10**((taoh(i)/taoy-1)**(xn))
C      dgam(i)=dgam0*10**(0.5*(taoh(i)/taoy)**(xn))
C         dgam(i)=dgam0*10**((taoh(i)/taoyhard-one)*(one/xn)) 
C         dgam(i)=dgam(i)*sign(one,tao(i))
C      dgamax=max(dgamax,abs(dgam(i)))      
      enddo 

C C     Strain rate tensor Dij
C       do i = 1,3
C           do j = 1,3
C               erate(1) = erate(1) + 0.5*dgam(1)*(so1(i)*xmo1(j)
C      .                 +so1(j)*xmo1(i))        
C               erate(2) = erate(2) + 0.5*dgam(2)*(so2(i)*xmo2(j)
C      .                 +so2(j)*xmo2(i))    
C               erate(3) = erate(3) + 0.5*dgam(3)*(so3(i)*xmo3(j)
C      .                 +so3(j)*xmo3(i))    
C               erate(4) = erate(4) + 0.5*dgam(4)*(so4(i)*xmo4(j)
C      .                 +so4(j)*xmo4(i))    
C               erate(5) = erate(5) + 0.5*dgam(5)*(so5(i)*xmo5(j)
C      .                 +so5(j)*xmo5(i))    
C               erate(6) = erate(6) + 0.5*dgam(6)*(so6(i)*xmo6(j)
C      .                 +so6(j)*xmo6(i))    
C           end do
C       end do 
      
C      xg12=xg120*(0.035*log10(dgamax/dgam0)+1)  

C      print *, 'xg12=', xg12     
C       print *, 'dgam=', dgam
C       print *, 'erate=', erate       
C       print *, 'dgam(3)=', dgam(3)
      
C     define time derivative of plastic defgrad dFp(i,j) 
C     initialize dFp
CCCCCCC Wrong in the previous code,can't directly use dfp
CCCCCCC otherwise dfp will be at least doubled 

      do i = 1,3
         do j = 1,3
             dfp(i,j) = zero
             ddfp(i,j) = zero
         enddo 
      enddo 

      do i = 1,3
          do j = 1,3
                  ddfp(i,j) = ddfp(i,j)+dgam(1)*so1(i)*xmo1(j)
     .                       +dgam(2)*so2(i)*xmo2(j)
     .                       +dgam(3)*so3(i)*xmo3(j)
     .                       +dgam(4)*so4(i)*xmo4(j)
     .                       +dgam(5)*so5(i)*xmo5(j)
     .                       +dgam(6)*so6(i)*xmo6(j)     
          enddo
      enddo 
      
       

      do i=1,3
         do j=1,3
            do k=1,3
               dfp(i,j) = dfp(i,j) + ddfp(i,k)*fp(k,j)
            enddo
         enddo
      enddo



C     update Fp, Forward Euler 
C     don't need to initalize Fp       
      do i = 1,3
          do j = 1,3
              fp(i,j) = fp(i,j) + dfp(i,j)*dt
          enddo 
      enddo  

C ========= store cumulated strain ===================C   
C      
	   
      do i = 1,6
         gamma(i) = gamma(i) + dgam(i)*dt	   
      enddo 
	   
C      print *, 'gamma=', gamma
      
      stateNew(km,56)= gamma(1)
      stateNew(km,57)= gamma(2)
      stateNew(km,58)= gamma(3)
      stateNew(km,59)= gamma(4)
      stateNew(km,60)= gamma(5)	  
      stateNew(km,61)= gamma(6)	
C ========= finish ===================================C 

      

C     Update state variable for Fp  fp  
      
      stateNew(km,1) = fp(1,1)
      stateNew(km,2) = fp(2,2)
      stateNew(km,3) = fp(3,3)
      stateNew(km,4) = fp(1,2)
      stateNew(km,5) = fp(2,3)
      stateNew(km,6) = fp(3,1)
      stateNew(km,7) = fp(2,1)
      stateNew(km,8) = fp(3,2)
      stateNew(km,9) = fp(1,3)
      
C     Calculate elastic part of defgrad
C     Store total Defgrad in Array F
   
      f(1,1) = defgradNew(km,1) 
      f(2,2) = defgradNew(km,2) 
      f(3,3) = defgradNew(km,3)
      f(1,2) = defgradNew(km,4)
      f(2,3) = defgradNew(km,5)
      f(3,1) = defgradNew(km,6)
      f(2,1) = defgradNew(km,7)
      f(3,2) = defgradNew(km,8)
      f(1,3) = defgradNew(km,9)
     

      
      
C     calculate inverse of updated Fp named Fpi
      
C     1 calculate the determinate of Fp named detFp
      detfp = fp(1,1)*(fp(2,2)*fp(3,3)-fp(3,2)*fp(2,3))
      detfp = detfp - fp(1,2)*(fp(2,1)*fp(3,3)-fp(2,3)*fp(3,1))
      detfp = detfp + fp(1,3)*(fp(2,1)*fp(3,2)-fp(2,2)*fp(3,1))
c

C     2 calculate the Fpi 


      fpi(1,1) = (fp(2,2)*fp(3,3)-fp(2,3)*fp(3,2))/detfp
      fpi(1,2) = -(fp(1,2)*fp(3,3)-fp(1,3)*fp(3,2))/detfp
      fpi(1,3) = (fp(1,2)*fp(2,3)-fp(1,3)*fp(2,2))/detfp
      fpi(2,1) = -(fp(2,1)*fp(3,3)-fp(2,3)*fp(3,1))/detfp
      fpi(2,2) = (fp(1,1)*fp(3,3)-fp(1,3)*fp(3,1))/detfp
      fpi(2,3) = -(fp(1,1)*fp(2,3)-fp(1,3)*fp(2,1))/detfp
      fpi(3,1) = (fp(2,1)*fp(3,2)-fp(2,2)*fp(3,1))/detfp
      fpi(3,2) = -(fp(1,1)*fp(3,2)-fp(1,2)*fp(3,1))/detfp
      fpi(3,3) = (fp(1,1)*fp(2,2)-fp(1,2)*fp(2,1))/detfp

 
C     calculate Fs (Fe) 
C     initialize Fs

      do i = 1,3
          do j = 1,3
              fs(i,j) = zero
          enddo
      enddo 

      do i = 1,3
          do j = 1,3
              do k = 1,3
                  fs(i,j) = fs(i,j)+f(i,k)*fpi(k,j)
              enddo 
          enddo
      enddo 
      
      
C     calculate E 
C     initialize E 
      
      do i = 1,3
          do j = 1,3
              e(i,j) = zero
          enddo
      enddo 
   

C  
      
      do i = 1,3
          do j = 1,3
              do k = 1,3
                 e(i,j)=e(i,j)+half*(fs(k,i)*fs(k,j))
              enddo
          enddo
      enddo


C  

      do i = 1,3
          do j = 1,3
              e(i,j) =e(i,j)-half*xiden(i,j)
          enddo
      enddo
      

C  

C     calculate 2nd PK stress Pk2

C     Lame's constant Lamda and Nu is used 
C     sigma = lamda*tr(E)*I+2*Nu*E
C     initialize PK2 stress 

cCCCCCvecterize e

      do i = 1,6
         vs(i) = zero
      enddo 
      
      ve(1) = e(1,1)
      ve(2) = e(2,2)
      ve(3) = e(3,3)
      ve(4) = e(1,2)
      ve(5) = e(2,3) 
      ve(6) = e(3,1) 
      
      do i = 1,6 
         do j = 1,6
            vs(i) = vs(i)+sm(i,j)*ve(j)
         enddo
      enddo
      

     
CCCCCC store vs 
      
      pk2(1,1) = vs(1)
      pk2(2,2) = vs(2)
      pk2(3,3) = vs(3)
      pk2(1,2) = vs(4)
      pk2(2,3) = vs(5)
      pk2(3,1) = vs(6) 
      pk2(1,3) = pk2(3,1) 
      pk2(2,1) = pk2(1,2) 
      pk2(3,2) = pk2(2,3) 
      
      stateNew(km,46) = pk2(2,2)
      stateNew(km,47) = pk2(3,3) 
      stateNew(km,48) = pk2(1,1)
    



C     rotate to update Cauchy stress 

C     initialize sig, clear the Old value
      do i = 1,3
          do j = 1,3
              sig(i,j) = zero
          enddo
      enddo

C     assign new values 
      do i = 1,3
          do j = 1,3
              do k = 1,3
                  do l = 1,3 
                     sig(i,j) = sig(i,j)+fs(i,k)*pk2(k,l)*fs(j,l)
                  enddo                  
              enddo
          enddo
      enddo
 

      
      detfs = fs(1,1)*(fs(2,2)*fs(3,3)-fs(3,2)*fs(2,3))
      detfs = detfs - fs(1,2)*(fs(2,1)*fs(3,3)-fs(2,3)*fs(3,1))
      detfs = detfs + fs(1,3)*(fs(2,1)*fs(3,2)-fs(2,2)*fs(3,1))
      
      
 

      
c     divide by J             
      do i = 1,3
         do j = 1,3
            sig(i,j) = sig(i,j)/detfs
         enddo 
      enddo 

C     Inverse u   

C     Inverse u   

      u(1,1) = stretchNew(km,1)
      u(2,2) = stretchNew(km,2) 
      u(3,3) = stretchNew(km,3)
      u(1,2) = stretchNew(km,4)
      u(2,3) = stretchNew(km,5)
      u(3,1) = stretchNew(km,6)
      u(2,1) = u(1,2)
      u(1,3) = u(3,1)
      u(3,2) = u(2,3)
      
      detu = u(1,1)*(u(2,2)*u(3,3)-u(3,2)*u(2,3))
      detu = detu-u(1,2)*(u(2,1)*u(3,3)-u(3,1)*u(2,3))
      detu = detu+u(1,3)*(u(2,1)*u(3,2)-u(3,1)*u(2,2))
      
      ui(1,1) = (u(2,2)*u(3,3)-u(2,3)*u(3,2))/detu
      ui(1,2) = -(u(1,2)*u(3,3)-u(1,3)*u(3,2))/detu
      ui(1,3) = (u(1,2)*u(2,3)-u(1,3)*u(2,2))/detu
      ui(2,1) = -(u(2,1)*u(3,3)-u(2,3)*u(3,1))/detu
      ui(2,2) = (u(1,1)*u(3,3)-u(1,3)*u(3,1))/detu
      ui(2,3) = -(u(1,1)*u(2,3)-u(1,3)*u(2,1))/detu
      ui(3,1) = (u(2,1)*u(3,2)-u(2,2)*u(3,1))/detu
      ui(3,2) = -(u(1,1)*u(3,2)-u(1,2)*u(3,1))/detu
      ui(3,3) = (u(1,1)*u(2,2)-u(1,2)*u(2,1))/detu
      
C     calculate r 
      do i = 1,3
         do j = 1,3
            r(i,j) = zero
         enddo
      enddo
      
      do i = 1,3
         do j = 1,3
            do k = 1,3
               r(i,j) = r(i,j)+f(i,k)*ui(k,j)
            enddo
         enddo
      enddo     
      
C     cauchy stress to corotational system      

      do i = 1,3
         do j = 1,3
            csig(i,j) = zero 
         enddo
      enddo
      
      do i = 1,3
         do j = 1,3
            do k = 1,3
               do l = 1,3
                  csig(i,j) = csig(i,j)+r(k,i)*sig(k,l)*r(l,j)
               enddo
            enddo
         enddo
      enddo
      
C     Update StressNew and StateNew for Stress 
C     Not very clear how to deal with sig(3,3) 
 
 
      stressNew(km,1) = csig(1,1)
      stressNew(km,2) = csig(2,2)
      stressNew(km,3) = csig(3,3) 
      stressNew(km,4) = csig(1,2)
      stressNew(km,5) = csig(2,3)
      stressNew(km,6) = csig(3,1) 
      
      
      stateNew(km,10) = sig(1,1)
      stateNew(km,11) = sig(2,2)
      stateNew(km,12) = sig(3,3) 
      stateNew(km,13) = sig(1,2)
      stateNew(km,14) = sig(2,3)
      stateNew(km,15) = sig(3,1) 
      
C     Update slip sets s m 
C     clear the previous value 
      
      do i = 1,3
         s1(i)  = zero
         xm1(i) = zero
         s2(i)  = zero
         xm2(i) = zero
         s3(i)  = zero
         xm3(i) = zero
         s4(i)  = zero
         xm4(i) = zero
         s5(i)  = zero
         xm5(i) = zero 
         s6(i)  = zero 
         xm6(i) = zero 
      enddo
           
      
      do i = 1,3
          do j = 1,3
              s1(i) = s1(i) + fs(i,j)*so1(j)
              s2(i) = s2(i) + fs(i,j)*so2(j)
              s3(i) = s3(i) + fs(i,j)*so3(j)
              s4(i) = s4(i) + fs(i,j)*so4(j)
              s5(i) = s5(i) + fs(i,j)*so5(j)
              s6(i) = s6(i) + fs(i,j)*so6(j)
          enddo
      enddo
      

C     2 calculate Fsi 
      
      fsi(1,1) = (fs(2,2)*fs(3,3)-fs(2,3)*fs(3,2))/detfs
      fsi(1,2) = -(fs(1,2)*fs(3,3)-fs(1,3)*fs(3,2))/detfs
      fsi(1,3) = (fs(1,2)*fs(2,3)-fs(1,3)*fs(2,2))/detfs
      fsi(2,1) = -(fs(2,1)*fs(3,3)-fs(2,3)*fs(3,1))/detfs
      fsi(2,2) = (fs(1,1)*fs(3,3)-fs(1,3)*fs(3,1))/detfs
      fsi(2,3) = -(fs(1,1)*fs(2,3)-fs(1,3)*fs(2,1))/detfs
      fsi(3,1) = (fs(2,1)*fs(3,2)-fs(2,2)*fs(3,1))/detfs
      fsi(3,2) = -(fs(1,1)*fs(3,2)-fs(1,2)*fs(3,1))/detfs
      fsi(3,3) = (fs(1,1)*fs(2,2)-fs(1,2)*fs(2,1))/detfs
 

      
      do i = 1,3
          do j = 1,3
              xm1(i) = xm1(i) + xmo1(j)*Fsi(j,i)
              xm2(i) = xm2(i) + xmo2(j)*Fsi(j,i)
              xm3(i) = xm3(i) + xmo3(j)*Fsi(j,i)
              xm4(i) = xm4(i) + xmo4(j)*Fsi(j,i)
              xm5(i) = xm5(i) + xmo5(j)*Fsi(j,i)
              xm6(i) = xm6(i) + xmo6(j)*Fsi(j,i)              
          enddo
      enddo 


C     update StateV s,m 

      do i = 1,3
          stateNew(km,i+15) = s1(i)
          stateNew(km,i+18) = xm1(i) 
          stateNew(km,i+21) = s2(i) 
          stateNew(km,i+24) = xm2(i) 
          stateNew(km,i+27) = s3(i)
          stateNew(km,i+30) = xm3(i)     
          stateNew(km,i+33) = s4(i) 
          stateNew(km,i+36) = xm4(i)
          stateNew(km,i+39) = s5(i) 
          stateNew(km,i+42) = xm5(i)   
          stateNew(km,i+49) = s6(i)
          stateNew(km,i+52) = xm6(i)          
      end do 
      
C     calculte total energy 

      do i = 1,3
         do j = 1,3
             ep(i,j) = zero
         enddo 
      enddo 

      do i = 1,3
         do j = 1,3
            do k = 1,3
                do l = 1,3
                     ep(i,j) = ep(i,j) + fs(i,k)*ddfp(k,l)*fsi(l,j)
                 enddo
             enddo
          enddo
       enddo

      ep(1,2) = half*(ep(1,2) + ep(2,1))
      ep(2,1) = ep(1,2)
      
      ep(2,3) = half*(ep(2,3) + ep(3,2))
      ep(3,2) = ep(2,3)
      
      ep(3,1) = half*(ep(1,3) + ep(3,1))
      ep(1,3) = ep(3,1)
      
      plas = stateOld(km,49)
        do i = 1,3
          do j = 1,3
                plas = plas + sig(i,j)*ep(i,j)*dt
          enddo
        enddo 
        
      stateNew(km,49) = plas

      elas=zero      
      do i = 1,3
          do j = 1,3
                elas   = elas + half*pk2(i,j)*e(i,j)
          enddo
       enddo 
      elas = elas + half*pk233*e33
      
      enerInternNew(km) = elas/den + plas/density(km)

   
  100 continue
 
      endif
      return
      end