C Infinitesimal damage model 
C by Ravi Kiran Yellavajjala
C Date: 
      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 )
      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     for the ferrite phase     
      IF (cmname(1:4) .EQ. 'FERR') THEN
	      Call vumat_Dam(
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     for the martensite phase     
      ELSE IF(cmname(1:4) .EQ. 'MART') THEN
	      Call vumat_Dam1(
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     for the retained austenite phase
      ELSE IF(cmname(1:6) .EQ. 'RA_MAT') THEN
	      Call vumat_Pha(
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 )
	 
      END IF
     
      RETURN
      END
C
      subroutine vumat_Dam(
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 )
      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
      REAL(8) :: E_mod, v, sig_0, eps_e_tr(1,6), zz, H, Q, b
      REAL(8) :: Ra, p_strain, fac, dev_str(1,6), s_tri(1,6)
      REAL(8) :: sig_y, phi, pre, tot, kb, mu, vol_strain
      REAL(8) :: cons, K_tp1, K_tp2, sig_v, K_tp, Dt1, dg
      REAL(8) :: p_strain1, q_tri, q_tri1, dir(1,6)
      REAL(8) :: s_new(1,6), cr_strain, K_sig, plas_strain(1,6)
      REAL(8) :: trans_strain(1,6), C2, red_flag, hh
      REAL(8) :: H_a, H_m, Q_a, Q_m, A_con, B_con, sig_max, p_max 
      REAL(8) :: vol_exp, K_vol, tot_strain(1,6), sig_a, sig_m
      REAL(8) :: dil_strain(1,3), K_tv, tr_dil, aa_con, stress_thresold
      REAL(8) :: flag_stress, mmm, cons1, phidot, cons3, m, n1, dam
      REAL(8) :: d_dam, A1, B1, C1, D1, F1, F2, der(2,2), res(2,1), dam1
      REAL(8) :: max_iter, T_sig, det_der, der_inv(2,2), Xa(2,1), conv
      REAL(8) :: P1, P2, cri_dam, T1, lambda, lambda_cr, num, qtri_new 
      REAL(8) :: T_avg, E1, G1, dcon_dg, dcons1_dg, cc1, cc2, df1_dg
      REAL(8) :: dF1_dam, dF2_dg
      COMPLEX(8) :: dgc, damc, F1_dg, F1_dam, F2_dg, F2_dam, cons_c


      Parameter (Zero = 0.d0, One = 1.d0, Two = 2.d0, Three = 3.d0,
     1 Tole = 1.0D-6)
      hh=10.0d0**-30.0d0

      max_iter=40.0
C------material properties---------
C E_mod - Youngs modulus
C v - Poissons Ratio
      fac=sqrt(Three/Two)
      E_mod=200000.0d0
      v=0.30d0
	  

C--sig_a,sig_m - yield stress of austenite, martensite  
      sig_0=props(1)

C-- H, Q, b - isotropic hardening parameters
      H=props(2)
      Q=props(3)
      b=props(4)
C--- power hard mod
      n1=props(5)
C----Strain rate parameters
      phidot=props(6)
      m=props(7)
C-- cr_strain-strain at which damage initiates	  	  
      cr_strain=props(8)
      cri_dam=props(9)
      lambda_cr=props(10)
C----kb-bulk mod. and mu-shear mod.	  
      kb=E_mod/(3.0d0*(1.0d0-(2.0d0*v))) 
      mu=E_mod/(2.0d0*(1.0d0+v))

	  
      do i=1,nblock   
C------initializing state variables---------
C----stateOld(i,1:6)-elastic strain
C----stateOld(i,7)-p strain
C----stateOld(i,8)-martensitic vol fraction
C----stateOld(i,9)-old Von Mises stress
         tot=totalTime
	     if (tot.eq.0.0d0) then
              stateOld(i,1)=0.0d0
              stateOld(i,2)=0.0d0
              stateOld(i,3)=0.0d0
              stateOld(i,4)=0.0d0
              stateOld(i,5)=0.0d0
              stateOld(i,6)=0.0d0
			  
              stateOld(i,7)=0.0d0
              stateOld(i,8)=0.0d0
              stateOld(i,9)=0.0d0
              stateOld(i,10)=0.0d0
              stateOld(i,11)=0.0d0
              stateOld(i,12)=0.0d0
			  
              stateOld(i,13)=0.0d0
              stateOld(i,14)=0.0d0
			  
              stateOld(i,15)=1.0d0
              stateOld(i,16)=1
              stateOld(i,17)=0.0d0



         end if 
C----finding total strain-----

         tot_strain(1,1)=stateOld(i,1)+strainInc(i,1)
         tot_strain(1,2)=stateOld(i,2)+strainInc(i,2)
         tot_strain(1,3)=stateOld(i,3)+strainInc(i,3)
         tot_strain(1,4)=stateOld(i,4)+strainInc(i,4)
         tot_strain(1,5)=stateOld(i,5)+strainInc(i,5)
         tot_strain(1,6)=stateOld(i,6)+strainInc(i,6)
		 
C----finding elastic strain-----
         eps_e_tr(1,1)=tot_strain(1,1)-stateOld(i,7)
         eps_e_tr(1,2)=tot_strain(1,2)-stateOld(i,8)
         eps_e_tr(1,3)=tot_strain(1,3)-stateOld(i,9)
         eps_e_tr(1,4)=tot_strain(1,4)-stateOld(i,10)
         eps_e_tr(1,5)=tot_strain(1,5)-stateOld(i,11)
         eps_e_tr(1,6)=tot_strain(1,6)-stateOld(i,12)
		 
         p_strain=stateOld(i,13)
		 

		 
         dam=stateOld(i,14)
         lambda=stateOld(i,15)
         T_sig=stateOld(i,17)
         num=stateOld(i,18)

C-----volumetric component of eps_e_tr
         vol_strain=(eps_e_tr(1,1)+eps_e_tr(1,2)
     1 +eps_e_tr(1,3))
C-----deviatoric component of eps_e_tr--------------------
         dev_str(1,1)=eps_e_tr(1,1)-(One/Three)*vol_strain
         dev_str(1,2)=eps_e_tr(1,2)-(One/Three)*vol_strain
         dev_str(1,3)=eps_e_tr(1,3)-(One/Three)*vol_strain
         dev_str(1,4)=eps_e_tr(1,4)
         dev_str(1,5)=eps_e_tr(1,5)
         dev_str(1,6)=eps_e_tr(1,6)
		 
C-----Determination of pressure------------------------------
         pre=kb*(tot_strain(1,1)+tot_strain(1,2)+
     1 tot_strain(1,3))
C-----Determination of trial deviatoric stress (s_tri)--------
         s_tri(1,1)=Two*mu*dev_str(1,1)
         s_tri(1,2)=Two*mu*dev_str(1,2)
         s_tri(1,3)=Two*mu*dev_str(1,3)
         s_tri(1,4)=Two*mu*dev_str(1,4)
         s_tri(1,5)=Two*mu*dev_str(1,5)
         s_tri(1,6)=Two*mu*dev_str(1,6)
		 
C-----Determination of trial J2 stress-----------------
         q_tri=fac*sqrt(s_tri(1,1)**2.0d0+s_tri(1,2)**2.0d0
     1 +s_tri(1,3)**2.0d0+Two*(s_tri(1,4)**2.0d0+s_tri(1,5)**2.0d0
     2 +s_tri(1,6)**2.0d0))
         q_tri1=q_tri/(fac)
C-----direction vector
         dir=s_tri/((q_tri1))

C-----isotropic hardening law
         Ra=(H*p_strain**n1)+Q*(One-exp(-b*p_strain))
C         call mprint(Ra,1,1,1,'Ra',6)
C-----Determination of current yield strength
         sig_y=sig_0+Ra
         phi=((q_tri)-sig_y)*(1.0d0-dam)
C         call mprint(q_tri,1,1,1,'q_tri',6)
C         call mprint(phi,1,1,1,'phi',6)
      if (phi.LT.1.0D-6) then
C-----Update the state variables
             stateNew(i,1)=tot_strain(1,1)
             stateNew(i,2)=tot_strain(1,2)
             stateNew(i,3)=tot_strain(1,3)
             stateNew(i,4)=tot_strain(1,4)
             stateNew(i,5)=tot_strain(1,5)
             stateNew(i,6)=tot_strain(1,6)	
			 
             stateNew(i,7)=stateOld(i,7)
             stateNew(i,8)=stateOld(i,8)
             stateNew(i,9)=stateOld(i,9)
             stateNew(i,10)=stateOld(i,10)
             stateNew(i,11)=stateOld(i,11)
             stateNew(i,12)=stateOld(i,12)
			 
             stateNew(i,13)=stateOld(i,13)
             stateNew(i,14)=stateOld(i,14)
			 
             stateNew(i,15)=stateOld(i,15)
             stateNew(i,16)=stateOld(i,16)
			 stateNew(i,17)=(dabs(pre)/q_tri)

			 
             



C-----Update the stress
             stressNew(i,1)=(s_tri(1,1)+pre)
             stressNew(i,2)=(s_tri(1,2)+pre)
             stressNew(i,3)=(s_tri(1,3)+pre)
             stressNew(i,4)=s_tri(1,4)
             stressNew(i,5)=s_tri(1,5)
             stressNew(i,6)=s_tri(1,6)
         end if
		 
C------transformation & Plasticity step--------
         if (phi.GE.1.0D-6) then

		     red_flag=0.0d0
		     dg=0.0d0
             d_dam=0.0d0
             dam1=dam
             Dt1=dt

             conv=1.0d0
             
            
			
             A1=0.313d0
             B1=2.013d0
             G1=0.0000082887d0
             C1=-1.46d0
             D1=3.58d0
             P1=0.250d0
             P2=0.250d0
			 
             
             T1=0.0d0

			 

             do while ((conv.EQ.1.0d0).AND.(red_flag.LT.max_iter))
			 
			     red_flag=red_flag+1.0d0
				 
                 cons1=sig_0+(H*(p_strain+dg)**n1)+Q*(One-
     1 exp(-b*(p_strain+dg)))
				 
                 dcon_dg=Dt1*phidot*(1.0d0/m)*(((q_tri-Three*mu*dg/
     1 (1.0d0-dam1))/(cons1)))**((1.0d0/m)-1.0d0)
	             dcons1_dg=H*n1*(p_strain+dg)**(n1-1.0d0)+b*Q*
     1 exp(-b*(p_strain+dg))
	             cc1=(-3.0d0*cons1*mu)/(1.0d0-dam1)
                 cc2=(q_tri-Three*mu*dg/
     1 (1.0d0-dam1))
	 
	             dF1_dg=1.0d0-dcon_dg*(cc1-dcons1_dg*cc2)/cons1**2.0d0
                 dF1_dam=dcon_dg*((3.0d0*mu*dg)/cons1)*
     1 (1.0d0-dam1)**(-2.0d0)
                 dF2_dg=-G1*E1*EXP(E1*T_sig)-
     1 (P1*exp(P2*(p_strain+dg))*T1)-dg*P1*T1*P2*exp(P2*(p_strain+dg))
				 
				 
				 

                 F1=dg-Dt1*phidot*(((q_tri-Three*mu*dg/(1.0d0-dam1)
     1 )/(cons1)))**(1.0d0/m)
				 
                 E1=(A1*EXP(B1*T_sig))
                 F2=dam1-dam-dg*(0.0006*T_sig**2.0d0-0.0003*T_sig
     1 -5.0d0**-5.0d0)-dg*(P1*exp(P2*(p_strain+dg))*T1)
                 dgc=CMPLX(dg,hh)
                 damc=CMPLX(dam1,hh)
				 
                 cons_c=sig_0+(H*(p_strain+dgc)**n1)+Q*(One-
     1 exp(-b*(p_strain+dgc)))
                 F1_dg=dgc-Dt1*phidot*(((q_tri-Three*mu*dgc/(1.0d0-dam1)
     1 )/(cons_c)))**(1.0d0/m)
                 F1_dam=dg-Dt1*phidot*(((q_tri-Three*mu*dg/(1.0d0-damc)
     1 )/(cons1)))**(1.0d0/m)


                 F2_dg=dam1-dam-dgc*(0.0006*T_sig**2.0d0-0.0003*T_sig
     1 -5.0d0**-5.0d0)-dgc*(P1*exp(P2*(p_strain+dgc))*T1)
                 F2_dam=damc-dam-dg*(0.0006*T_sig**2.0d0-0.0003*T_sig
     1 -5.0d0**-5.0d0)-dg*(P1*exp(P2*(p_strain+dg))*T1)
	 

                 der(1,1)=AIMAG(F1_dg)/hh
                 der(1,2)=AIMAG(F1_dam)/hh
                 der(2,1)=AIMAG(F2_dg)/hh
                 der(2,2)=AIMAG(F2_dam)/hh

				 


                 res(1,1)=F1
                 res(2,1)=F2
				 
                 call two_inverse(der,der_inv)
                 det_der=der(1,1)*der(2,2)-der(1,2)*der(2,1)
                 if ((isnan(det_der)).OR.(red_flag.GE.max_iter)) then
C			         call mprint(iter,1,1,1,'isnan encountered',6)
					 dg=0.0d0
					 dam1=dam
				     EXIT
                 end if
                 Xa=matmul(der_inv, res)
                 dg=dg-Xa(1,1)
                 dam1=dam1-Xa(2,1)

				 
                 if ((dabs(Xa(1,1)).lt.tole*dabs(dg)).and.
     1 (dabs(F1).lt.tole*sig_0)) then
                     if ((dabs(Xa(2,1)).lt.tole*dabs(dam1)).and.
     1 (dabs(F2).lt.tole*sig_0)) then
	                     conv=0.0d0
                     end if
                 end if

             end do
             if (red_flag.GE.max_iter) then
				dg=0.0d0
                dam1=dam
             end if
             lambda=lambda+(-1.195d0*T_sig+3.83d0)*dg
             if (dam1.LT.0.0d0) then
				
                dam1=0.0d0
             end if


C------update eff plastic strain--------
             p_strain=p_strain+dg
			 
             if ((p_strain).LT.(cr_strain)) then
		 
                 dam1=0.0d0
             end if
			 


C-----updating plastic strain
C             plas_strain=(fac*dg*dir)/(1.0d0-dam1)
             plas_strain=(fac*dg*dir)
	 
C----updating the deviatoric stress
             s_new=(s_tri-(Two*mu*fac*dg*dir)*(1.0d0-dam1))
             q_tri_new=fac*sqrt(s_new(1,1)**2.0d0+s_new(1,2)**2.0d0
     1 +s_new(1,3)**2.0d0+Two*(s_new(1,4)**2.0d0+s_new(1,5)**2.0d0
     2 +s_new(1,6)**2.0d0))
C             T_avg=(T_sig*num+(dabs(pre)/q_tri_new))/(num+1.0d0)
             T_sig=(dabs(pre)/q_tri_new)
C-----Update the elastic strain

             stateNew(i,1)=tot_strain(1,1)
             stateNew(i,2)=tot_strain(1,2)
             stateNew(i,3)=tot_strain(1,3)
             stateNew(i,4)=tot_strain(1,4)
             stateNew(i,5)=tot_strain(1,5)
             stateNew(i,6)=tot_strain(1,6)	
			 
             stateNew(i,7)=stateOld(i,7)+plas_strain(1,1)
             stateNew(i,8)=stateOld(i,8)+plas_strain(1,2)
             stateNew(i,9)=stateOld(i,9)+plas_strain(1,3)
             stateNew(i,10)=stateOld(i,10)+plas_strain(1,4)
             stateNew(i,11)=stateOld(i,11)+plas_strain(1,5)
             stateNew(i,12)=stateOld(i,12)+plas_strain(1,6)
			 
             stateNew(i,13)=p_strain
             stateNew(i,14)=dam1
			 
             stateNew(i,15)=lambda
             stateNew(i,16)=stateOld(i,16)
			 
             if (((dam1).GT.(cri_dam)).or.((lambda).GT.(lambda_cr))) 
     1 then
		 
                 stateNew(i,16)=0
             end if

			 stateNew(i,17)=T_sig



C-----updating the stress-------
             stressNew(i,1)=s_new(1,1)+pre*(1.0d0-dam1)
             stressNew(i,2)=s_new(1,2)+pre*(1.0d0-dam1)
             stressNew(i,3)=s_new(1,3)+pre*(1.0d0-dam1)
             stressNew(i,4)=s_new(1,4)
             stressNew(i,5)=s_new(1,5)
             stressNew(i,6)=s_new(1,6)

			 
         end if 
      end do 
      return
      end subroutine vumat_Dam
	  
      subroutine vumat_Dam1(
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 )
      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
      REAL(8) :: E_mod, v, sig_0, eps_e_tr(1,6), zz, H, Q, b
      REAL(8) :: Ra, p_strain, fac, dev_str(1,6), s_tri(1,6)
      REAL(8) :: sig_y, phi, pre, tot, kb, mu, vol_strain
      REAL(8) :: cons, K_tp1, K_tp2, sig_v, K_tp, Dt1, dg
      REAL(8) :: p_strain1, q_tri, q_tri1, dir(1,6)
      REAL(8) :: s_new(1,6), cr_strain, K_sig, plas_strain(1,6)
      REAL(8) :: trans_strain(1,6), C2, red_flag, hh
      REAL(8) :: H_a, H_m, Q_a, Q_m, A_con, B_con, sig_max, p_max 
      REAL(8) :: vol_exp, K_vol, tot_strain(1,6), sig_a, sig_m
      REAL(8) :: dil_strain(1,3), K_tv, tr_dil, aa_con, stress_thresold
      REAL(8) :: flag_stress, mmm, cons1, phidot, cons3, m, n1, dam
      REAL(8) :: d_dam, A1, B1, C1, D1, F1, F2, der(2,2), res(2,1), dam1
      REAL(8) :: max_iter, T_sig, det_der, der_inv(2,2), Xa(2,1), conv
      REAL(8) :: P1, P2, cri_dam, T1, lambda, lambda_cr, num, qtri_new 
      REAL(8) :: T_avg, E1, G1, dcon_dg, dcons1_dg, cc1, cc2, df1_dg
      REAL(8) :: dF1_dam, dF2_dg
      COMPLEX(8) :: dgc, damc, F1_dg, F1_dam, F2_dg, F2_dam, cons_c


      Parameter (Zero = 0.d0, One = 1.d0, Two = 2.d0, Three = 3.d0,
     1 Tole = 1.0D-6)
      hh=10.0d0**-30.0d0

      max_iter=40.0
C------material properties---------
C E_mod - Youngs modulus
C v - Poissons Ratio
      fac=sqrt(Three/Two)
      E_mod=200000.0d0
      v=0.30d0
	  

C--sig_a,sig_m - yield stress of austenite, martensite  
      sig_0=props(1)

C-- H, Q, b - isotropic hardening parameters
      H=props(2)
      Q=props(3)
      b=props(4)
C--- power hard mod
      n1=props(5)
C----Strain rate parameters
      phidot=props(6)
      m=props(7)
C-- cr_strain-strain at which damage initiates	  	  
      cr_strain=props(8)
      cri_dam=props(9)
      lambda_cr=props(10)
C----kb-bulk mod. and mu-shear mod.	  
      kb=E_mod/(3.0d0*(1.0d0-(2.0d0*v))) 
      mu=E_mod/(2.0d0*(1.0d0+v))

	  
      do i=1,nblock   
C------initializing state variables---------
C----stateOld(i,1:6)-elastic strain
C----stateOld(i,7)-p strain
C----stateOld(i,8)-martensitic vol fraction
C----stateOld(i,9)-old Von Mises stress
         tot=totalTime
	     if (tot.eq.0.0d0) then
              stateOld(i,1)=0.0d0
              stateOld(i,2)=0.0d0
              stateOld(i,3)=0.0d0
              stateOld(i,4)=0.0d0
              stateOld(i,5)=0.0d0
              stateOld(i,6)=0.0d0
			  
              stateOld(i,7)=0.0d0
              stateOld(i,8)=0.0d0
              stateOld(i,9)=0.0d0
              stateOld(i,10)=0.0d0
              stateOld(i,11)=0.0d0
              stateOld(i,12)=0.0d0
			  
              stateOld(i,13)=0.0d0
              stateOld(i,14)=0.0d0
			  
              stateOld(i,15)=1.0d0
              stateOld(i,16)=1
              stateOld(i,17)=0.0d0



         end if 
C----finding total strain-----

         tot_strain(1,1)=stateOld(i,1)+strainInc(i,1)
         tot_strain(1,2)=stateOld(i,2)+strainInc(i,2)
         tot_strain(1,3)=stateOld(i,3)+strainInc(i,3)
         tot_strain(1,4)=stateOld(i,4)+strainInc(i,4)
         tot_strain(1,5)=stateOld(i,5)+strainInc(i,5)
         tot_strain(1,6)=stateOld(i,6)+strainInc(i,6)
		 
C----finding elastic strain-----
         eps_e_tr(1,1)=tot_strain(1,1)-stateOld(i,7)
         eps_e_tr(1,2)=tot_strain(1,2)-stateOld(i,8)
         eps_e_tr(1,3)=tot_strain(1,3)-stateOld(i,9)
         eps_e_tr(1,4)=tot_strain(1,4)-stateOld(i,10)
         eps_e_tr(1,5)=tot_strain(1,5)-stateOld(i,11)
         eps_e_tr(1,6)=tot_strain(1,6)-stateOld(i,12)
		 
         p_strain=stateOld(i,13)
		 

		 
         dam=stateOld(i,14)
         lambda=stateOld(i,15)
         T_sig=stateOld(i,17)
         num=stateOld(i,18)

C-----volumetric component of eps_e_tr
         vol_strain=(eps_e_tr(1,1)+eps_e_tr(1,2)
     1 +eps_e_tr(1,3))
C-----deviatoric component of eps_e_tr--------------------
         dev_str(1,1)=eps_e_tr(1,1)-(One/Three)*vol_strain
         dev_str(1,2)=eps_e_tr(1,2)-(One/Three)*vol_strain
         dev_str(1,3)=eps_e_tr(1,3)-(One/Three)*vol_strain
         dev_str(1,4)=eps_e_tr(1,4)
         dev_str(1,5)=eps_e_tr(1,5)
         dev_str(1,6)=eps_e_tr(1,6)
		 
C-----Determination of pressure------------------------------
         pre=kb*(tot_strain(1,1)+tot_strain(1,2)+
     1 tot_strain(1,3))
C-----Determination of trial deviatoric stress (s_tri)--------
         s_tri(1,1)=Two*mu*dev_str(1,1)
         s_tri(1,2)=Two*mu*dev_str(1,2)
         s_tri(1,3)=Two*mu*dev_str(1,3)
         s_tri(1,4)=Two*mu*dev_str(1,4)
         s_tri(1,5)=Two*mu*dev_str(1,5)
         s_tri(1,6)=Two*mu*dev_str(1,6)
		 
C-----Determination of trial J2 stress-----------------
         q_tri=fac*sqrt(s_tri(1,1)**2.0d0+s_tri(1,2)**2.0d0
     1 +s_tri(1,3)**2.0d0+Two*(s_tri(1,4)**2.0d0+s_tri(1,5)**2.0d0
     2 +s_tri(1,6)**2.0d0))
         q_tri1=q_tri/(fac)
C-----direction vector
         dir=s_tri/((q_tri1))

C-----isotropic hardening law
         Ra=(H*p_strain**n1)+Q*(One-exp(-b*p_strain))
C         call mprint(Ra,1,1,1,'Ra',6)
C-----Determination of current yield strength
         sig_y=sig_0+Ra
         phi=((q_tri)-sig_y)*(1.0d0-dam)
C         call mprint(q_tri,1,1,1,'q_tri',6)
C         call mprint(phi,1,1,1,'phi',6)
      if (phi.LT.1.0D-6) then
C-----Update the state variables
             stateNew(i,1)=tot_strain(1,1)
             stateNew(i,2)=tot_strain(1,2)
             stateNew(i,3)=tot_strain(1,3)
             stateNew(i,4)=tot_strain(1,4)
             stateNew(i,5)=tot_strain(1,5)
             stateNew(i,6)=tot_strain(1,6)	
			 
             stateNew(i,7)=stateOld(i,7)
             stateNew(i,8)=stateOld(i,8)
             stateNew(i,9)=stateOld(i,9)
             stateNew(i,10)=stateOld(i,10)
             stateNew(i,11)=stateOld(i,11)
             stateNew(i,12)=stateOld(i,12)
			 
             stateNew(i,13)=stateOld(i,13)
             stateNew(i,14)=stateOld(i,14)
			 
             stateNew(i,15)=stateOld(i,15)
             stateNew(i,16)=stateOld(i,16)
			 stateNew(i,17)=(dabs(pre)/q_tri)

			 
             



C-----Update the stress
             stressNew(i,1)=(s_tri(1,1)+pre)
             stressNew(i,2)=(s_tri(1,2)+pre)
             stressNew(i,3)=(s_tri(1,3)+pre)
             stressNew(i,4)=s_tri(1,4)
             stressNew(i,5)=s_tri(1,5)
             stressNew(i,6)=s_tri(1,6)
         end if
		 
C------transformation & Plasticity step--------
         if (phi.GE.1.0D-6) then

		     red_flag=0.0d0
		     dg=0.0d0
             d_dam=0.0d0
             dam1=dam
             Dt1=dt

             conv=1.0d0
             
            
			
             A1=0.313d0
             B1=2.013d0
             G1=0.0000082887d0
             C1=-1.46d0
             D1=3.58d0
             P1=0.250d0
             P2=0.250d0
			 
             
             T1=0.0d0

			 

             do while ((conv.EQ.1.0d0).AND.(red_flag.LT.max_iter))
			 
			     red_flag=red_flag+1.0d0
				 
                 cons1=sig_0+(H*(p_strain+dg)**n1)+Q*(One-
     1 exp(-b*(p_strain+dg)))
				 
                 dcon_dg=Dt1*phidot*(1.0d0/m)*(((q_tri-Three*mu*dg/
     1 (1.0d0-dam1))/(cons1)))**((1.0d0/m)-1.0d0)
	             dcons1_dg=H*n1*(p_strain+dg)**(n1-1.0d0)+b*Q*
     1 exp(-b*(p_strain+dg))
	             cc1=(-3.0d0*cons1*mu)/(1.0d0-dam1)
                 cc2=(q_tri-Three*mu*dg/
     1 (1.0d0-dam1))
	 
	             dF1_dg=1.0d0-dcon_dg*(cc1-dcons1_dg*cc2)/cons1**2.0d0
                 dF1_dam=dcon_dg*((3.0d0*mu*dg)/cons1)*
     1 (1.0d0-dam1)**(-2.0d0)
                 dF2_dg=-G1*E1*EXP(E1*T_sig)-
     1 (P1*exp(P2*(p_strain+dg))*T1)-dg*P1*T1*P2*exp(P2*(p_strain+dg))
				 
				 
				 

                 F1=dg-Dt1*phidot*(((q_tri-Three*mu*dg/(1.0d0-dam1)
     1 )/(cons1)))**(1.0d0/m)
				 
                 E1=(A1*EXP(B1*T_sig))
                 F2=dam1-dam-dg*(0.0007*T_sig**2.0d0-0.0002*T_sig
     1 -3.0d0**-6.0d0)-dg*(P1*exp(P2*(p_strain+dg))*T1)
                 dgc=CMPLX(dg,hh)
                 damc=CMPLX(dam1,hh)
				 
                 cons_c=sig_0+(H*(p_strain+dgc)**n1)+Q*(One-
     1 exp(-b*(p_strain+dgc)))
                 F1_dg=dgc-Dt1*phidot*(((q_tri-Three*mu*dgc/(1.0d0-dam1)
     1 )/(cons_c)))**(1.0d0/m)
                 F1_dam=dg-Dt1*phidot*(((q_tri-Three*mu*dg/(1.0d0-damc)
     1 )/(cons1)))**(1.0d0/m)


                 F2_dg=dam1-dam-dgc*(0.0007*T_sig**2.0d0-0.0002*T_sig
     1 -3.0d0**-6.0d0)-dgc*(P1*exp(P2*(p_strain+dgc))*T1)
                 F2_dam=damc-dam-dg*(0.0007*T_sig**2.0d0-0.0002*T_sig
     1 -3.0d0**-6.0d0)-dg*(P1*exp(P2*(p_strain+dg))*T1)
	 

                 der(1,1)=AIMAG(F1_dg)/hh
                 der(1,2)=AIMAG(F1_dam)/hh
                 der(2,1)=AIMAG(F2_dg)/hh
                 der(2,2)=AIMAG(F2_dam)/hh

				 


                 res(1,1)=F1
                 res(2,1)=F2
				 
                 call two_inverse(der,der_inv)
                 det_der=der(1,1)*der(2,2)-der(1,2)*der(2,1)
                 if ((isnan(det_der)).OR.(red_flag.GE.max_iter)) then
C			         call mprint(iter,1,1,1,'isnan encountered',6)
					 dg=0.0d0
					 dam1=dam
				     EXIT
                 end if
                 Xa=matmul(der_inv, res)
                 dg=dg-Xa(1,1)
                 dam1=dam1-Xa(2,1)

				 
                 if ((dabs(Xa(1,1)).lt.tole*dabs(dg)).and.
     1 (dabs(F1).lt.tole*sig_0)) then
                     if ((dabs(Xa(2,1)).lt.tole*dabs(dam1)).and.
     1 (dabs(F2).lt.tole*sig_0)) then
	                     conv=0.0d0
                     end if
                 end if

             end do
             if (red_flag.GE.max_iter) then
				dg=0.0d0
                dam1=dam
             end if
             if (dam1.LT.0.0d0) then
				
                dam1=0.0d0
             end if
             lambda=lambda+(-1.195d0*T_sig+3.83d0)*dg



C------update eff plastic strain--------
             p_strain=p_strain+dg
			 
             if ((p_strain).LT.(cr_strain)) then
		 
                 dam1=0.0d0
             end if
			 


C-----updating plastic strain
C             plas_strain=(fac*dg*dir)/(1.0d0-dam1)
             plas_strain=(fac*dg*dir)
	 
C----updating the deviatoric stress
             s_new=(s_tri*-(Two*mu*fac*dg*dir)(1.0d0-dam1))
             q_tri_new=fac*sqrt(s_new(1,1)**2.0d0+s_new(1,2)**2.0d0
     1 +s_new(1,3)**2.0d0+Two*(s_new(1,4)**2.0d0+s_new(1,5)**2.0d0
     2 +s_new(1,6)**2.0d0))
C             T_avg=(T_sig*num+(dabs(pre)/q_tri_new))/(num+1.0d0)
             T_sig=(dabs(pre)/q_tri_new)
C-----Update the elastic strain

             stateNew(i,1)=tot_strain(1,1)
             stateNew(i,2)=tot_strain(1,2)
             stateNew(i,3)=tot_strain(1,3)
             stateNew(i,4)=tot_strain(1,4)
             stateNew(i,5)=tot_strain(1,5)
             stateNew(i,6)=tot_strain(1,6)	
			 
             stateNew(i,7)=stateOld(i,7)+plas_strain(1,1)
             stateNew(i,8)=stateOld(i,8)+plas_strain(1,2)
             stateNew(i,9)=stateOld(i,9)+plas_strain(1,3)
             stateNew(i,10)=stateOld(i,10)+plas_strain(1,4)
             stateNew(i,11)=stateOld(i,11)+plas_strain(1,5)
             stateNew(i,12)=stateOld(i,12)+plas_strain(1,6)
			 
             stateNew(i,13)=p_strain
             stateNew(i,14)=dam1
			 
             stateNew(i,15)=lambda
             stateNew(i,16)=stateOld(i,16)
			 
             if (((dam1).GT.(cri_dam)).or.((lambda).GT.(lambda_cr))) 
     1 then
		 
                 stateNew(i,16)=0
             end if

			 stateNew(i,17)=T_sig



C-----updating the stress-------
             stressNew(i,1)=s_new(1,1)+pre*(1.0d0-dam1)
             stressNew(i,2)=s_new(1,2)+pre*(1.0d0-dam1)
             stressNew(i,3)=s_new(1,3)+pre*(1.0d0-dam1)
             stressNew(i,4)=s_new(1,4)
             stressNew(i,5)=s_new(1,5)
             stressNew(i,6)=s_new(1,6)

			 
         end if 
      end do 
      return
      end subroutine vumat_Dam1
	  
      subroutine vumat_Pha(
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 )
      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
      REAL(8) :: E_mod, v, sig_0, eps_e_tr(1,6), zz, H, Q, b
      REAL(8) :: Ra, p_strain, fac, dev_str(1,6), s_tri(1,6)
      REAL(8) :: sig_y, phi, pre, tot, kb, mu, vol_strain
      REAL(8) :: cons, K_tp1, K_tp2, sig_v, K_tp, Dt1, dg
      REAL(8) :: der, p_strain1, q_tri, q_tri1, dir(1,6)
      REAL(8) :: s_new(1,6), cr_strain, K_sig, plas_strain(1,6)
      REAL(8) :: trans_strain(1,6), C1, C2, res, red_flag, hh
      REAL(8) :: H_a, H_m, Q_a, Q_m, A_con, B_con, sig_max, p_max 
      REAL(8) :: vol_exp, K_vol, tot_strain(1,6), sig_a, sig_m
      REAL(8) :: dil_strain(1,3), K_tv, tr_dil, aa_con, stress_thresold
      REAL(8) :: flag_stress, mmm, cons1, phidot, cons3, m, n1
      REAL(8) :: cons2_real, dcons_dg, dcons2_dg, dres_dg
      COMPLEX(8) :: ii, dgc, cons2, resc


      Parameter (Zero = 0.d0, One = 1.d0, Two = 2.d0, Three = 3.d0,
     1 Toler = 1.0D-6)
      hh=10.0d0**-30.0d0
      ii=CMPLX(0.0d0,1.0d0)

C------material properties---------
C E_mod - Youngs modulus
C v - Poissons Ratio
      fac=sqrt(Three/Two)
      E_mod=200000.0d0
      v=0.30d0
C--sig_a,sig_m - yield stress of austenite, martensite  
      sig_a=props(1)
      sig_m=props(2)
C-- H, Q, b - isotropic hardening parameters
      H_a=props(3)
      H_m=props(4)
      Q_a=props(5)
      Q_m=props(6)
      b=props(7)
C-- K_tp1,K_tp2,K_tv - martensite vol evol parameters
      K_tp1=props(8)
      K_tp2=props(9)
      K_tv=props(10)
C-- cr_strain,sig_max,p_max - martensite vol evol parameters	  	  
      cr_strain=props(11)
      sig_max=props(12)
      p_max=props(13)
C--- power hard mod
      n1=props(14)
C----Strain rate parameters
      phidot=props(15)
      m=props(16)
C----kb-bulk mod. and mu-shear mod.	  
      kb=E_mod/(3.0d0*(1.0d0-(2.0d0*v))) 
      mu=E_mod/(2.0d0*(1.0d0+v))


	  
      do i=1,nblock   
C------initializing state variables---------
C----stateOld(i,1:6)-elastic strain
C----stateOld(i,7)-p strain
C----stateOld(i,8)-martensitic vol fraction
C----stateOld(i,9)-old Von Mises stress
         tot=totalTime
	     if (tot.eq.0.0d0) then
              stateOld(i,1)=0.0d0
              stateOld(i,2)=0.0d0
              stateOld(i,3)=0.0d0
              stateOld(i,4)=0.0d0
              stateOld(i,5)=0.0d0
              stateOld(i,6)=0.0d0
			  
              stateOld(i,7)=0.0d0
              stateOld(i,8)=0.0d0
              stateOld(i,9)=0.0d0
              stateOld(i,10)=0.0d0
              stateOld(i,11)=0.0d0
              stateOld(i,12)=0.0d0
			  
              stateOld(i,13)=0.0d0
              stateOld(i,14)=0.0d0
              stateOld(i,15)=0.0d0
              stateOld(i,16)=0.0d0
              stateOld(i,17)=0.0d0
              stateOld(i,18)=0.0d0
			  
              stateOld(i,19)=0.0d0
			  
              stateOld(i,20)=0.0d0
			  
              stateOld(i,21)=0.0d0
			  
              stateOld(i,22)=0.0d0
              stateOld(i,23)=0.0d0
              stateOld(i,24)=0.0d0
			  stateOld(i,25)=0.0d0


         end if 
C----finding total strain-----

         tot_strain(1,1)=stateOld(i,1)+strainInc(i,1)
         tot_strain(1,2)=stateOld(i,2)+strainInc(i,2)
         tot_strain(1,3)=stateOld(i,3)+strainInc(i,3)
         tot_strain(1,4)=stateOld(i,4)+strainInc(i,4)
         tot_strain(1,5)=stateOld(i,5)+strainInc(i,5)
         tot_strain(1,6)=stateOld(i,6)+strainInc(i,6)
		 
C----finding elastic strain-----
         eps_e_tr(1,1)=tot_strain(1,1)-stateOld(i,7)-stateOld(i,13)
         eps_e_tr(1,2)=tot_strain(1,2)-stateOld(i,8)-stateOld(i,14)
         eps_e_tr(1,3)=tot_strain(1,3)-stateOld(i,9)-stateOld(i,15)
         eps_e_tr(1,4)=tot_strain(1,4)-stateOld(i,10)-stateOld(i,16)
         eps_e_tr(1,5)=tot_strain(1,5)-stateOld(i,11)-stateOld(i,17)
         eps_e_tr(1,6)=tot_strain(1,6)-stateOld(i,12)-stateOld(i,18)
		 
         p_strain=stateOld(i,19)
         zz=stateOld(i,20)
         dil_strain(1,1)=stateOld(i,22)
         dil_strain(1,2)=stateOld(i,23)
         dil_strain(1,3)=stateOld(i,24)
C-------trace of transformation dilational strain---------
         tr_dil=dil_strain(1,1)+dil_strain(1,2)+dil_strain(1,3)
C-----volumetric component of eps_e_tr
         vol_strain=(eps_e_tr(1,1)+eps_e_tr(1,2)
     1 +eps_e_tr(1,3))
C-----deviatoric component of eps_e_tr--------------------
         dev_str(1,1)=eps_e_tr(1,1)-(One/Three)*vol_strain
         dev_str(1,2)=eps_e_tr(1,2)-(One/Three)*vol_strain
         dev_str(1,3)=eps_e_tr(1,3)-(One/Three)*vol_strain
         dev_str(1,4)=eps_e_tr(1,4)
         dev_str(1,5)=eps_e_tr(1,5)
         dev_str(1,6)=eps_e_tr(1,6)
		 
C-----Determination of pressure------------------------------
         pre=kb*((tot_strain(1,1)+tot_strain(1,2)+tot_strain(1,3))-
     1 (One/Three)*tr_dil)
C-----Determination of trial deviatoric stress (s_tri)--------
         s_tri(1,1)=Two*mu*dev_str(1,1)
         s_tri(1,2)=Two*mu*dev_str(1,2)
         s_tri(1,3)=Two*mu*dev_str(1,3)
         s_tri(1,4)=Two*mu*dev_str(1,4)
         s_tri(1,5)=Two*mu*dev_str(1,5)
         s_tri(1,6)=Two*mu*dev_str(1,6)
		 
C-----Determination of trial J2 stress-----------------
         q_tri=fac*sqrt(s_tri(1,1)**2.0d0+s_tri(1,2)**2.0d0
     1 +s_tri(1,3)**2.0d0+Two*(s_tri(1,4)**2.0d0+s_tri(1,5)**2.0d0
     2 +s_tri(1,6)**2.0d0))
         q_tri1=q_tri/fac
C-----direction vector
         dir=s_tri/q_tri1
C-----Homogenization of hardening properties---------
         Q=(1.0d0-zz)*Q_a+zz*Q_m
         H=(1.0d0-zz)*H_a+zz*H_m
         sig_0=(1.0d0-zz)*sig_a+zz*sig_m
C-----isotropic hardening law
         Ra=(H*p_strain**n1)+Q*(One-exp(-b*p_strain))
C         call mprint(Ra,1,1,1,'Ra',6)
C-----Determination of current yield strength
         sig_y=sig_0+Ra
         phi=q_tri-sig_y
C         call mprint(q_tri,1,1,1,'q_tri',6)

      if (phi.LT.1.0D-6) then
C-----Update the state variables
             stateNew(i,1)=tot_strain(1,1)
             stateNew(i,2)=tot_strain(1,2)
             stateNew(i,3)=tot_strain(1,3)
             stateNew(i,4)=tot_strain(1,4)
             stateNew(i,5)=tot_strain(1,5)
             stateNew(i,6)=tot_strain(1,6)	
			 
             stateNew(i,7)=stateOld(i,7)
             stateNew(i,8)=stateOld(i,8)
             stateNew(i,9)=stateOld(i,9)
             stateNew(i,10)=stateOld(i,10)
             stateNew(i,11)=stateOld(i,11)
             stateNew(i,12)=stateOld(i,12)
			 
             stateNew(i,13)=stateOld(i,13)
             stateNew(i,14)=stateOld(i,14)
             stateNew(i,15)=stateOld(i,15)
             stateNew(i,16)=stateOld(i,16)
             stateNew(i,17)=stateOld(i,17)
             stateNew(i,18)=stateOld(i,18)
			 
             stateNew(i,19)=stateOld(i,19)
             stateNew(i,20)=stateOld(i,20)
             stateNew(i,21)=q_tri
			 
             stateNew(i,22)=stateOld(i,22)
             stateNew(i,23)=stateOld(i,23)
             stateNew(i,24)=stateOld(i,24)




C-----Update the stress
             stressNew(i,1)=s_tri(1,1)+pre
             stressNew(i,2)=s_tri(1,2)+pre
             stressNew(i,3)=s_tri(1,3)+pre
             stressNew(i,4)=s_tri(1,4)
             stressNew(i,5)=s_tri(1,5)
             stressNew(i,6)=s_tri(1,6)
         end if
		 
C------transformation & Plasticity step--------
         if (phi.GE.1.0D-6) then

		     red_flag=0.0d0
		     dg=0.0d0
		     sig_v=stateOld(i,21)

             Dt1=dt

			 
			 K_tp=K_tp1+K_tp2*sig_v
C             call mprint(K_tp,1,1,1,'K_tp',6)
		     cons=(1.0d0+Three*K_tp*mu*(Two*(One-zz))*zz*Dt1)
             cons1=sig_0+(H*p_strain**n1+Q*(One-exp
     1 (-b*p_strain)))
	         

             res=dg-Dt1*phidot*(((q_tri-Three*mu*dg)/(cons*cons1))
     1 )**(1.0d0/m)

             do while ((dabs(res).GT.Toler).AND.(red_flag.LT.30.0d0))
			     red_flag=red_flag+1.0d0
				 
                 dgc=dg+hh*ii
                 cons2=sig_0+(H*(p_strain+dgc)**n1+Q*(One-exp
     1 (-b*(p_strain+dgc))))
                 cons2_real=sig_0+(H*(p_strain+dg)**n1+Q*(One-exp
     1 (-b*(p_strain+dg))))
	             dcons_dg=Dt1*phidot*(1.0d0/m)*(((q_tri-Three*mu*dgc)/
     1 (cons*cons2)))**(1.0d0/m-1.0d0)
	             dcons2_dg=H*n1*(p_strain+dgc)**(n1-1.0d0)+
     1 b*Q*exp(-b*(p_strain+dg))
	             dres_dg=1.0d0-dcons_dg*(1.0d0/cons)*(-cons2_real*3.0d0*mu
     1 -dcons2_dg*(q_tri-Three*mu*dg))/cons2_real**2.0d0

                 resc=dgc-Dt1*phidot*(((q_tri-Three*mu*dgc)/(cons*cons2)
     1 ))**(1.0d0/m)

                 der=dres_dg
            

                 dg1=dg-res/der
                 p_strain1=p_strain+dg1
                 cons3=sig_0+(H*(p_strain1)**n1+Q*(One-exp
     1 (-b*(p_strain1))))

                 res=dg1-Dt1*phidot*(((q_tri-Three*mu*dg1)/(cons*cons3))
     1 )**(1.0d0/m)
                 dg=dg1

             end do
             if (red_flag.GE.30.0d0) then 
                 dg=0.0d0
             end if
			 

C------update eff plastic strain--------
             p_strain=p_strain+dg
C------Update martensite volume fraction
C------Evolution of Martensite volume fraction-------  
             A_con=sig_v/sig_max
             B_con=(p_strain-cr_strain)/(p_max-cr_strain)			 
             if (p_strain.GE.cr_strain) then
                 if (zz.LT.1.0d0) then
                     if (A_con.GT.3.0d0) then
				         A_con=3.0d0
                     end if
                     if (B_con.GT.3.0d0) then
				         B_con=3.0d0
                     end if

                     aa_con=((0.30d0*
     1 ((1.0d0/A_con)+(B_con))**1.86d0))-0.042d0
                     mmm=((0.5d0*
     1 ((1.0d0/A_con)+(B_con))**2.0d0))-0.188d0
					 zz=aa_con

					 dil_strain(1,1)=(One/Three)*K_tv*zz
					 dil_strain(1,2)=(One/Three)*K_tv*zz
					 dil_strain(1,3)=(One/Three)*K_tv*zz

                 end if 
                 if (zz.LT.0.0d0) then
                     zz=0.0d0
                 end if
                 if (zz.GE.1.0d0) then
				     zz=1.0d0
                 end if
				 
		     endif
             if (p_strain.LT.cr_strain) then
		         zz=0.0d0
             endif
C             call mprint(zz,1,1,1,'zz',6)
C-----updating plastic strain
             plas_strain=fac*dg*dir
	 
C----updating the deviatoric stress
             s_new=(s_tri-Two*mu*fac*dg*dir)/(One+Three*mu*K_tp*Two*
     1 (One-zz)*zz*Dt1)
	 
C-----updating transformation strain
             trans_strain=Three*K_tp*s_new*(One-zz)*zz*Dt1
C-----Update the elastic strain

             stateNew(i,1)=tot_strain(1,1)
             stateNew(i,2)=tot_strain(1,2)
             stateNew(i,3)=tot_strain(1,3)
             stateNew(i,4)=tot_strain(1,4)
             stateNew(i,5)=tot_strain(1,5)
             stateNew(i,6)=tot_strain(1,6)	
			 
             stateNew(i,7)=stateOld(i,7)+plas_strain(1,1)
             stateNew(i,8)=stateOld(i,8)+plas_strain(1,2)
             stateNew(i,9)=stateOld(i,9)+plas_strain(1,3)
             stateNew(i,10)=stateOld(i,10)+plas_strain(1,4)
             stateNew(i,11)=stateOld(i,11)+plas_strain(1,5)
             stateNew(i,12)=stateOld(i,12)+plas_strain(1,6)
			 
             stateNew(i,13)=stateOld(i,13)+trans_strain(1,1)
             stateNew(i,14)=stateOld(i,14)+trans_strain(1,2)
             stateNew(i,15)=stateOld(i,15)+trans_strain(1,3)
             stateNew(i,16)=stateOld(i,16)+trans_strain(1,4)
             stateNew(i,17)=stateOld(i,17)+trans_strain(1,5)
             stateNew(i,18)=stateOld(i,18)+trans_strain(1,6)
			 
             stateNew(i,22)=dil_strain(1,1)
             stateNew(i,23)=dil_strain(1,2)
             stateNew(i,24)=dil_strain(1,3)



C-----updating the stress-------
             stressNew(i,1)=s_new(1,1)+pre
             stressNew(i,2)=s_new(1,2)+pre
             stressNew(i,3)=s_new(1,3)+pre
             stressNew(i,4)=s_new(1,4)
             stressNew(i,5)=s_new(1,5)
             stressNew(i,6)=s_new(1,6)
C-----sig_v-------
             sig_v=fac*sqrt(s_new(1,1)**2.0d0+s_new(1,2)**2.0d0
     1 +s_new(1,3)**2.0d0+Two*(s_new(1,4)**2.0d0+s_new(1,5)
     2 **2.0d0+s_new(1,6)**2.0d0))
C----update other state variables--------
             stateNew(i,19)=p_strain
			 stateNew(i,20)=zz
			 stateNew(i,21)=sig_v
			 stateNew(i,25)=mmm
			 
         end if 
      end do 
      return
      end subroutine vumat_Pha
      
      
      SUBROUTINE two_inverse(jac,jac_inv)
      REAL(8) :: jac(2,2), jac1(2,2), jac_inv(2,2)	  
	  
      jac1(1,1)=jac(2,2)
      jac1(1,2)=-jac(1,2)
      jac1(2,1)=-jac(2,1)
      jac1(2,2)=jac(1,1)			
      jac_inv=(1.0d0/(jac(1,1)*jac(2,2)-jac(1,2)*jac(2,1)))*jac1
      end SUBROUTINE two_inverse
	  
      subroutine mprint(a,ii,jj,mm,name,iow)
c-----[--.----+----.----+----.-----------------------------------------]
c      Purpose: Output array of integer values

c      Inputs:
c         a(mm,*)  - Array to output
c         ii       - Number of rows to output
c         jj       - Number of columns to output
c         mm       - Dimension of array
c         name     - Name of array to appear with outputs

c      Outputs:
c         none
c-----[--.----+----.----+----.-----------------------------------------]
C      include 'vaba_param.inc'

      character name*(*), aname*30
      integer   ii,jj,mm,nn, ja,jb,j,i,n,irow,iow
      real*8    a(mm,*)

      
c     Print an ii x jj array whose dimension is mm for first subscript

      aname = name
      nn    = (jj+5)/6
      jb    = 0

      do n = 1,nn

        ja = jb + 1
        jb = min(jj,ja + 5)

c       Output to O file

        if(n.eq.1 .or. ii.gt.1) then
          write(iow,2000) aname,(j,j=ja,jb)
        endif
        do i = 1,ii
          if(ii.eq.1) then
            irow = n
          else
            irow = i
          endif
          write(iow,2001) irow,(a(i,j),j=ja,jb)
        end do ! i

      end do ! n

c     Flush buffer to force outputs (fix later)

c     Formats

 2000 format(/4x,'Matrix: ',a30/4x,'row/col',i6,5i12)
 2001 format(i8,1p,6e12.4)

      end
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! vuel Dan Gerbig
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine vuel(nblock,rhs,amass,dtimestable,svars,nsvars,
     *                energy,
     *                nnode,ndofel,props,nprops,jprops,njprops,
     *                coords,mcrd,u,du,v,a,
     *                jtype,jelem,
     *                time,period,dtimecur,dtimeprev,kstep,kinc,
     *                lflags,
     *                dmassscalefactor,
     *                predef,npredef,
     *                jdltyp, adlmag)

      include 'vaba_param.inc'

      ! operational code keys
      parameter ( jmasscalc            = 1,
     *            jintforceanddtstable = 2,
     *            jexternforce         = 3)

      ! flag indices
      parameter (iprocedure = 1,
     *           inlgeom    = 2,
     *           iopcode    = 3,
     *           nflags     = 3)

      ! energy array indices
      parameter ( ielpd = 1,
     *            ielcd = 2,
     *            ielie = 3,
     *            ielts = 4,
     *            ieldd = 5,
     *            ielbv = 6,
     *            ielde = 7,
     *            ielhe = 8,
     *            ielke = 9,
     *            ielth = 10,
     *            ieldmd = 11,
     *            ieldc = 12,
     *            nelenergy = 12)

      ! predefined variables indices
      parameter ( ipredvaluenew = 1,
     *            ipredvalueold = 2,
     *            npred         = 2)

      ! time indices
      parameter (isteptime  = 1,
     *           itotaltime = 2,
     *           ntime      = 2)

      ! procedure flags
      parameter ( jdynexplicit = 17)

      dimension rhs(nblock,ndofel),amass(nblock,ndofel,ndofel),
     *          dtimestable(nblock),
     *          svars(nblock,nsvars),energy(nblock,nelenergy),
     *          props(nprops),jprops(njprops),
     *          jelem(nblock),time(ntime),lflags(nflags),
     *          coords(nblock,nnode,mcrd),
     *          u(nblock,ndofel), du(nblock,ndofel),
     *          v(nblock,ndofel), a(nblock, ndofel),
     *          dmassscalefactor(nblock),  
     *          predef(nblock,nnode,npredef,npred),
     *          adlmag(nblock)

      ! local variables
      ! constants
      parameter ( zero = 0.d0, half = 0.5d0, one = 1.d0 )
      ! deformed filaments
      dimension dx(3), dy1(3), dy2(3), dy3(3), volAvF(3, 3)

      if ((jtype.eq.1).and.(lflags(iprocedure).eq.jdynexplicit)) then
        ! user spring for proportional loading

        ! properties
        spring_kx = props(1)
        prop_load_ratio = props(2)
        uspring_rho = props(3)
        uspring_a0 = props(4)
        i_acting_dof = jprops(1)

        i_m_offset = (i_acting_dof - 1)*3
        i_d_offset = 12

        if (lflags(iopcode).eq.jmasscalc) then
          ! define mass matrix

          do kblock = 1,nblock
            ! compute mass
            ! use original length
            alenx0 = coords(kblock, 5, 1) 
     *        - coords(kblock, i_acting_dof, 1)
            aleny0 = coords(kblock, 5, 2)
     *        - coords(kblock, i_acting_dof, 2)
            alenz0 = coords(kblock, 5, 3)
     *        - coords(kblock, i_acting_dof, 3)
            alen0 = sqrt(alenX0**2 + alenY0**2 + alenZ0**2)

            ! half mass
            am0 = alen0*uspring_a0*uspring_rho*half

            ! master node (on RVE)
            amass(kblock, i_m_offset + 1, i_m_offset + 1) = am0
            amass(kblock, i_m_offset + 2, i_m_offset + 2) = am0
            amass(kblock, i_m_offset + 3, i_m_offset + 3) = am0
            ! dummy node (on other end of spring)
            amass(kblock, i_d_offset + 1, i_d_offset + 1) = am0
            amass(kblock, i_d_offset + 2, i_d_offset + 2) = am0
            amass(kblock, i_d_offset + 3, i_d_offset + 3) = am0
          end do

        else if (lflags(iopcode).eq.jintforceanddtstable) then
          ! define force

          do kblock = 1,nblock
            ! compute mass
            ! use original length
            alenx0 = coords(kblock, 5, 1) 
     *        - coords(kblock, i_acting_dof, 1)
            aleny0 = coords(kblock, 5, 2)
     *        - coords(kblock, i_acting_dof, 2)
            alenz0 = coords(kblock, 5, 3)
     *        - coords(kblock, i_acting_dof, 3)
            alen0 = sqrt(alenx0**2 + aleny0**2 + alenz0**2)

            ! full mass
            am0 = alen0*uspring_a0*uspring_rho

            ! proportional loading

            ! calculate volume averaged deformation gradient
            ! reference filaments, from origin to master nodes
            ! are aligned with coordinate system
            dx(1) = coords(kblock, 1, 1)
            dx(2) = coords(kblock, 2, 2)
            dx(3) = coords(kblock, 3, 3)
            ! deformed filaments, from orgin to master nodes
            dy1(1) = u(kblock, 1) + dx(1)
            dy1(2) = u(kblock, 2)
            dy1(3) = u(kblock, 3)
            dy2(1) = u(kblock, 4)
            dy2(2) = u(kblock, 5) + dx(2)
            dy2(3) = u(kblock, 6)
            dy3(1) = u(kblock, 7)
            dy3(2) = u(kblock, 8)
            dy3(3) = u(kblock, 9) + dx(3)
            ! calc F
            call CalcDefGrad(dx, dy1, dy2, dy3, volAvF)

            ! e1 displacement of nmx, ndx
            u_nmx_1 = u(kblock, 1)
            u_ndx_1 = u(kblock, 10)
            ! stretch delta in spring x
            spring_d = u_ndx_1 - u_nmx_1
            
            ! internal force in spring x
            spring_fx = spring_d*spring_kx

            ! e{2, 3} displacement of nm{y, z}, nd{y, z}
            u_nmuser = u(kblock, i_m_offset + i_acting_dof)
            u_nduser = u(kblock, i_d_offset + i_acting_dof)
            ! stretch delta in uspring
            uspring_d = u_nduser - u_nmuser

            ! calculate user spring constant uspring_k
            !
            ! want to satisfy
            !   sigma_11*prop_load_ratio = sigma_aa
            ! where
            !   sigma = true stress tensor
            !   no sum on subscript aa
            ! transform to nominal (1st PK stress) S
            !   sigma = 1/J * F * S
            ! result is
            !   S_11*F_11*prop_load_ratio = S_aa*F_aa
            ! use nominal dimensions of RVE in dx() to get areas
            ! for stresses, and solve for force
            uspring_f = (prop_load_ratio*volAvF(1, 1)*spring_fx
     *        *dx(1))/(volAvF(i_acting_dof, i_acting_dof)
     *        *dx(i_acting_dof))
            if (uspring_d.eq.zero) then
              uspring_k = spring_kx*prop_load_ratio
            else
              uspring_k = uspring_f/uspring_d
            end if

            uspring_k = abs(uspring_k)

            ! stable time increment in the element
            dtimestable(kblock) = sqrt(am0/uspring_k)

            ! assemble internal load in RHS
            ! master node (on RVE)
            rhs(kblock, i_m_offset + i_acting_dof) = -uspring_f
            ! dummy node (on other end spring)
            rhs(kblock, i_d_offset + i_acting_dof) =  uspring_f

            ! internal energy calculation
            uspring_d_old = svars(kblock, 1)
            uspring_f_old = svars(kblock, 2)

            ! average force * displacement increment
            energy(kblock, ielie) = energy(kblock, ielie) +
     *        half*(uspring_f + uspring_f_old)*
     *        (uspring_d - uspring_d_old)

            svars(kblock, 1) = uspring_d
            svars(kblock, 2) = uspring_f
          end do

        end if

      end if
      
      return
      end subroutine vuel 

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! CalcDefGrad
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      subroutine CalcDefGrad(dx, dy1, dy2, dy3, F)
      implicit none

      ! transferred
      real*8, intent(in) :: dx(3)
      real*8, intent(in) :: dy1(3), dy2(3), dy3(3)
      real*8, intent(out) :: F(3, 3)

      ! calculate volume average F
      F(1, 1) = dy1(1)/dx(1)
      F(2, 1) = dy1(2)/dx(1)
      F(3, 1) = dy1(3)/dx(1)

      F(1, 2) = dy2(1)/dx(2)
      F(2, 2) = dy2(2)/dx(2)
      F(3, 2) = dy2(3)/dx(2)

      F(1, 3) = dy3(1)/dx(3)
      F(2, 3) = dy3(2)/dx(3)
      F(3, 3) = dy3(3)/dx(3)

      return
      end subroutine CalcDefGrad
	  
