subroutine ucure(m,n,nn,kcus,matus,dt,ak,density,volumi,curedat,
     1   tempbeg,tempend,deltime,time,curerate)
#ifdef _IMPLICITNONE
      implicit none
#else
          implicit logical (a-z)
#endif
#include "mrcparm.cmn"
        include 'concom'
        include 'creeps'
        dimension curedat(*),matus(2)
        real*8 ak, curedat, curerate, deltime, density, dt
        integer kcus, m, matus, n, nn
        real*8 tempbeg, tempend, time, volumi
        real*8 K1,K2,K1c,K2c,Kd,A1
        real*8 A2,Ad,n1,n2,mi,E1,E2
        real*8 Ed, H, R, vf, l, ba
        real*8 Tginf, Tgo, g, f,w,ako
        real*8 ank1(200,1000,1),ank2(200,16,1000)
        real*8 usp1(200,1000,1),usp2(200,16,1000)
        
c       temp(m,inc,1),dg(m,nn,inc) 

        common /ankon/ ank1,ank2
        common /uspch/ usp1,usp2
        common /initconv/ ako
        common /kin1/ A1,A2,Ad,E1,E2,Ed,n1,n2,mi,ba,g,w,H,R
c
C--VOLUME FRACTION
        vf=0.6d0
C--INITIAL CONVECTION
        ako=0.038642259   
C--CURE KINETICS       
        A1=17580.23346
        A2=21525.43198  
        Ad=6.48E+18
        E1=70500
        E2=59050
	Ed=136800
	n1=1.8
          n2=1.32  
        mi=1.397913015
        l=0.435d0
	ba=0.47
        Tginf=206d0
        Tgo=-11d0
        H=465000.0
	R=8.3144
        g=0.024810815
	w=0.00048
      f=w*(tempbeg-273.15-Tgo-(Tginf-Tgo)*l*dt/(1-(1-l)*dt))+g
        K1c=A1*exp(-E1/(R*tempbeg))
        K2c=A2*exp(-E2/(R*tempbeg))
	Kd=Ad*exp(-Ed/(R*tempbeg))*exp(-ba/f)
        K1=1/(1/Kd+1/K1c)
        K2=1/(1/Kd+1/K2c)
        curerate=K1*(1-dt)**n1+K2*dt**mi*(1-dt)**n2
	curedat(9)=H*(1-vf)
	volumi=1.0
	ak=dt+deltime*curerate
	if (ak>0.999) then
		curerate=0.0
		ak=1.0
     	 end if

      return
      end
      subroutine uspcht	(spheat,m,nn,kcus,inc1,ncycle1,matus,nstats,
     * temp0,dtemp,time,dtime)
c     m(1)        user element number
c     m(2)        is the internal element number
c     nn          integration point number
c     kcus(1)     layer number
c     kcus(2)     internal layer number
c     temp0       temperature at beginning of increment
c     dtemp       estimated temperature increment
c     time        time at begining of increment
c     dtime       time increment
#ifdef _IMPLICITNONE
      implicit none
#else
          implicit logical (a-z)
#endif
c     ** Start of generated type statements **
      real*8 dtemp, dtime
      integer inc1, matus, ncycle1, nn, nstats
      integer layers, intel, icode, kcus,m
      real*8 spheat, temp0, time,ak
      real*8 Tgo, Tginf,l, sigma, Crcp, Afcp, Bfcp
      real*8 Arcp, Brcp, Drcp, vf, p1, p2, pc, wf, icod
      real*8 cpr, cpf, Tg, tempend, dt, dtdl,ako
      real*8 A1,A2,Ad,E1,E2,Ed,n1,n2,mi,ba,H,R,g
      real*8 f,w,K1c,K2c,K1,K2,Kd,curerate,cv,t
      real*8 dg(200,16,1000)
      real*8 temp(200,1000,1)

      common /uspch/ temp,dg
      common /initconv/ ako
      common /kin1/ A1,A2,Ad,E1,E2,Ed,n1,n2,mi,ba,g,w,H,R
      
c      
       dimension m(2),matus(2),kcus(2)      

c      dg(m(1),nn,inc), temp(m(1),inc,1)

       Tginf=206d0
       Tgo=-11d0
       vf=0.6d0
       l=0.435d0
       sigma=16.5d0
       Crcp=1.1d0
       Afcp=0.0023E3
       Bfcp=0.765E3
       Arcp=0.0025E3
       Brcp=1.80E3
       Drcp=-0.25E3
       p1=1760 !initial fiber density
       p2=1110 !initial resin density
       pc=p2+vf*(p1-p2)!composite density
       wf=p1/pc*vf 
c
        icode=285 
c       layers=4 !number of internal layers
c       intel=4  !number of integration per layer
c        kcus=4

         

         call elmvar (icode,m(1),nn,kcus,ak)
    

        tempend=temp0+dtemp !temperature at the end 

        temp(m(1),inc1,1)=tempend 
                
       if (inc1>1.0)then
           if (ak-ako==0d0) then
             if (inc1==2d0) then
                 cv=ako
                 t=temp(m(1),inc1-1,1)
           Tg=Tgo+(Tginf-Tgo)*l*cv/(1-(1-l)*cv)
           f=w*(t-273.15-Tg)+g
           K1c=A1*exp(-E1/R/(t))
           K2c=A2*exp(-E2/R/(t))
	      Kd=Ad*exp(-Ed/R/(t))*exp(-ba/f)
           K1=1/(1/Kd+1/K1c)
           K2=1/(1/Kd+1/K2c)
           curerate=K1*(1-cv)**n1+K2*cv**mi*(1-cv)**n2
           ak=cv+dtime*curerate
              if (ak>0.999) then
		     curerate=0.0
		     ak=1.0
	         end if
           dg(m(1),nn,inc1)=ak
c             print*,'cv=',cv, curerate
            else
            t=temp(m(1),inc1-1,1)
           cv=dg(m(1),nn,inc1-1)
           Tg=Tgo+(Tginf-Tgo)*l*cv/(1-(1-l)*cv)
           f=w*(t-273.15-Tg)+g
           K1c=A1*exp(-E1/R/(t))
           K2c=A2*exp(-E2/R/(t))
	      Kd=Ad*exp(-Ed/R/(t))*exp(-ba/f)
           K1=1/(1/Kd+1/K1c)
           K2=1/(1/Kd+1/K2c)
           curerate=K1*(1-cv)**n1+K2*cv**mi*(1-cv)**n2
           ak=cv+dtime*curerate
                if (ak>0.999) then
		         curerate=0.0
		         ak=1.0
	          end if
           dg(m(1),nn,inc1)=ak
c           print*,'cv=', cv,t,curerate
           end if
          end if
        end if
c        print*,inc1, m(1),nn, kcus, lovl, ak
c
c
       Tg=Tgo+(Tginf-Tgo)*l*ak/(1-(1-l)*ak)+273.15 
       cpf=Afcp*(tempend-273.15)+Bfcp
       cpr=Arcp*(tempend-273.15)
       cpr=cpr+Brcp+Drcp/(1+exp(Crcp*(tempend-Tg-sigma)))
       spheat=wf*cpf+(1-wf)*cpr

        return
        end
        subroutine ankond(cond,caniso,n,nn,kc,matno,id,t,dt,time,deltim,
     1                  joulht)
c     cond         conductivity matrix
c     n            element number
c     nn           integration point number
c     kcus(1)      layer number
c     kcus(2)      internal layer number
c     matus(1)     user material identification number
c     matus(2)     internal material identification number
c     t           temperature
c     dt          incremental temperature
c     time        time
c     deltim      increment of time
#ifdef _IMPLICITNONE
      implicit none
#else
         implicit logical (a-z)
#endif
      include 'creeps'
      include 'concom'
c     ** Start of generated type statements **
        real*8 caniso, cond, deltim, dt
        integer id, joulht, kc, matno, n, nn
        real*8 t, time,ak
        real*8 Tg,Tgo,Tginf,l,A1,A2,Ad,E1,E2,Ed
        real*8 n1,n2,mi,ba,H,R,g,cv,f,w,tempbeg
        real*8 K1c,K2c,Kd,K1,K2,curerate,ako
        real*8 vf,Alf,Blf,Btf,akr,bkr,ckr,dkr,ekr,fkr
        real*8 K11,K22,K33,Klf,Ktf,Kr,Kr1,Kr2
        real*8 az,bz,bv,c,dv
        real*8 dg(200,16,1000)
        real*8 temp(200,1000,1)
         
        common /ankon/ temp,dg
        common /initconv/ ako
        common /kin1/ A1,A2,Ad,E1,E2,Ed,n1,n2,mi,ba,g,w,H,R
        

c       dg(n,nn,inc),temp(n,inc,1)

c     ** End of generated type statements **
        dimension cond(id,id)

        joulht=0
        call elmvar (285,n,nn,kc,ak)

         Tginf=206d0
         Tgo=-11d0
         vf=0.6d0
         l=0.435d0

         temp(n,inc,1)=t

        if (inc>1.0)then
           if (ak-ako==0d0) then
             if (inc==2d0) then
                 cv=ako
                 tempbeg=temp(n,inc-1,1)
           Tg=Tgo+(Tginf-Tgo)*l*cv/(1-(1-l)*cv)
           f=w*(tempbeg-273.15-Tg)+g
           K1c=A1*exp(-E1/R/(tempbeg))
           K2c=A2*exp(-E2/R/(tempbeg))
	      Kd=Ad*exp(-Ed/R/(tempbeg))*exp(-ba/f)
           K1=1/(1/Kd+1/K1c)
           K2=1/(1/Kd+1/K2c)
           curerate=K1*(1-cv)**n1+K2*cv**mi*(1-cv)**n2
           ak=cv+timinc*curerate
                if (ak>0.999) then
		       curerate=0.0
		       ak=1.0
	           end if
           dg(n,nn,inc)=ak
c            print*,'cv=',cv
            else
            tempbeg=temp(n,inc-1,1)
           cv=dg(n,nn,inc-1)
           Tg=Tgo+(Tginf-Tgo)*l*cv/(1-(1-l)*cv)
           f=w*(tempbeg-273.15-Tg)+g
           K1c=A1*exp(-E1/R/(tempbeg))
           K2c=A2*exp(-E2/R/(tempbeg))
	      Kd=Ad*exp(-Ed/R/(tempbeg))*exp(-ba/f)
           K1=1/(1/Kd+1/K1c)
           K2=1/(1/Kd+1/K2c)
           curerate=K1*(1-cv)**n1+K2*cv**mi*(1-cv)**n2
           ak=cv+timinc*curerate
               if (ak>0.999) then
		        curerate=0.0
		        ak=1.0
	          end if
           dg(n,nn,inc)=ak
c            print*,'cv=', cv,tempbeg,curerate
           end if
          end if
         end if

         !required input
       
         Alf=0.0074
         Blf=9.7
         Btf=0.84
         akr=0.0008
         bkr=-0.0011
         ckr= -0.0002
         dkr=-0.0937
         ekr=0.22
         fkr=0.12 !intercept!
         
          !calculations
       
          Kr=akr*(t-273.15)*ak**2+bkr*(t-273.15)*ak
          Kr=Kr+ckr*(t-273.15)+dkr*ak**2+ekr*ak+fkr
          Ktf=Btf
          Klf=Alf*(t-273.15)+Blf
          
          K11=vf*Klf+(1-vf)*Kr

c        farmer model

         az=vf*Kr*(Ktf/Kr-1)+Kr*(0.5-Ktf/(2*Kr)) 
         bz=Kr*(Ktf/Kr-1)
         bz=bz*(vf**2-vf+(Ktf/Kr+1)**2/(2*Ktf/Kr-2)**2)**0.5

         K22=az+bz
        

         K33=K22

        !thermal conductivity matrix
        cond(1,1)=k11
        cond(2,2)=k22
        cond(3,3)=k33
        cond(1,2)=0
        cond(1,3)=0
        cond(2,1)=0
        cond(2,3)=0
        cond(3,1)=0
        cond(3,2)=0
          
        !pause
        return
        end
      subroutine usinc(f,n,ndeg,iflag)
#ifdef _IMPLICITNONE
      implicit none
#else
         implicit logical (a-z)
#endif
c     ** Start of generated type statements **
      real*8 f,inittemp,Ti
      integer iflag, n, ndeg
      common /f/Ti
      iflag=3  !temperature
C--INITTEMP
      f=393    !initial temperaure
      Ti=f
      !print*,'usinc',Ti
      !pause
      return
      end
      subroutine ufilm(uhfilm,utsink,uhnatur,uexpnat,
     *  ueffview,uemiss,uqflux,ts,n,time,inc)
      implicit real*8 (a-h,o-z)
      dimension ts(*),n(*)
C--FILM      
      uhfilm=0    !film coefficient
      !utsink=298  !ambient temperature
      uexpnat=1   !exponent
      ueffview=1  ! effective view factor
      uhnatur=5
      T1=428.1
      T2=488
      Dt=10300
      r=2
      r=r/60
      if (TIME<=(T1-393)/r) then
          UTSINK=393+r*TIME
      end if
      if (TIME>(T1-393)/r.AND.TIME<=Dt+(T1-393)/r) then
        UTSINK=T1
       end if
      if (TIME>Dt+(T1-393)/r.AND.TIME<=(T2-T1)/r+Dt+(T1-393)/r) then
         UTSINK=T1+r*(TIME-(Dt+(T1-393)/r))
      end if
      if (TIME>(T2-T1)/r+Dt+(T1-393)/r) then
	UTSINK=T2
      end if
      uemiss=0
      uqflux=0
c      pause
      return
      end
      subroutine forcdt(x1,x2,x3,f,t,time,dtime,ndeg,node,
     *                  x4,xord,ncrd,iacflg,inc,ipass)
      implicit real*8 (a-h,o-z)
      include 'creeps'
      dimension f(ndeg),t(ndeg),xord(ncrd)
      common /f/Ti
      COMMON StoredINC,Tbound,DTmax,Zeroflag,writeflag,tcure,Temp
c	1 TIME,Tempstore
c      common /curetime/amin
       T1=428.1
       T2=488
       Dt=10300
       r=2
       r=r/60
        if (TIME<=(T1-393)/r) then
	T=393+r*TIME
         end if
	if (TIME>(T1-393)/r.AND.TIME<=Dt+(T1-393)/r) then
 		T=T1
	 end if
	if (TIME>Dt+(T1-393)/r.AND.TIME<=(T2-T1)/r+Dt+(T1-393)/r) then
 	    T=T1+r*(TIME-(Dt+(T1-393)/r))
	end if
	if (TIME>(T2-T1)/r+Dt+(T1-393)/r) then
            T=T2
      end if
      Tbound=T(1)
      return
      end
          SUBROUTINE PLOTV (V,S,SP,ETOT,EPLAS,ECREEP,T,M,NN,KCUS,NDI,
     + NSHEAR,JPLTCD)
	IMPLICIT REAL *8 (A-H, O-Z)
	DIMENSION S(*),SP(*),ETOT(*),EPLAS(*),ECREEP(*),T(*),
     + M(2),KCUS(2)
	COMMON StoredINC, Tbound, DTmax,amin,a1,Zeroflag,writeflag,tcure,
     1 Temp,TIME,tc1
	include '../common/concom'
	include '../common/creeps'
	if (INC==0.AND.Zeroflag==0) then
		DTmax=-8065014
		amin=8065014
		Zeroflag=1
	end if
	if (INC>StoredInc) then
		amin=8065014
		StoredInc=INC
	end if
	call elmvar(285,M,NN,KCUS,conv)
	if (conv<amin .AND. conv>0) then
		amin=conv
      end if         
	RETURN
       END
      SUBROUTINE UPSTNO (NQCODE,NODEID,VALNO,NQNCOMP,NQTYPE,
     1	NQAVER,NQCOMPTYPE,NQDATATYPE,NQCOMPNAME)
		IMPLICIT REAL*8 (A-H,O-Z)
          DIMENSION  Y78(500),x21(1),y21(1)
          DIMENSION A78(3,3), B78(3,1),AINV(3,3),COFACTOR(3,3),COEF(3,1)
		COMMON StoredINC, Tbound, DTmax,amin,Zeroflag,writeflag,tcure,Temp,
	1	TIME,Tempstore,a1,x01,flag,i,j,step,x211,NINC,DET
          COMMON x1,x2,x3,x4,p1,p2,p3,p4,A,B,C,D,y1,y3,y4,y2,Y78,A78,B78
     1    ,COFACTOR,AINV,COEF,x21,y21
		include '../common/concom'
		include '../common/creeps'
		include '../common/blnk'
		include '../common/lass'
		DIMENSION VALNO(*)
     		CHARACTER*24 NQCOMPNAME(*)
		call nodvar(14,NODEID,Temp,NQNCOMP,NQDATATYPE)
		VALNO(1)=Temp
          step=50
          NINC=500
          If ((Temp-Tbound>Y78(INC)) .AND. INC>0 ) then
              Y78(INC)=Temp-Tbound
              !Print*, Y78(INC),Temp,Tbound, INC
              !pause
          else if (Y78(INC)>100) then
              Y78(INC)=0
          end if
          if (INC==NINC .AND. flag==0 ) then
c          print*, maxval(Y78) ,maxloc(Y78)
c          pause
          x21=maxloc(Y78)
          y21=maxval(Y78)
          A78(2,2)=x21(1)*step
          B78(2,1)=y21(1)
          A78(1,2)=A78(2,2)-step
          B78(1,1)=Y78(x21(1)-1)
          A78(3,2)=A78(2,2)+step
          B78(3,1)=Y78(x21(1)+1)
          A78(1,3)= A78(1,2)**2
          A78(2,3)= A78(2,2)**2
          A78(3,3)= A78(3,2)**2
          A78(1,1)=1
          A78(2,1)=1
          A78(3,1)=1
c          print*, B78(1,1),A78(1,2),B78(2,1),A78(2,2),B78(3,1),A78(3,2)
c          pause
          DET =   A78(1,1)*A78(2,2)*A78(3,3)  
     1       - A78(1,1)*A78(2,3)*A78(3,2)  
     1       - A78(1,2)*A78(2,1)*A78(3,3)  
     1       + A78(1,2)*A78(2,3)*A78(3,1)  
     1       + A78(1,3)*A78(2,1)*A78(3,2)  
     1       - A78(1,3)*A78(2,2)*A78(3,1)
          COFACTOR(1,1) = +(A78(2,2)*A78(3,3)-A78(2,3)*A78(3,2))
          COFACTOR(1,2) = -(A78(2,1)*A78(3,3)-A78(2,3)*A78(3,1))
          COFACTOR(1,3) = +(A78(2,1)*A78(3,2)-A78(2,2)*A78(3,1))
          COFACTOR(2,1) = -(A78(1,2)*A78(3,3)-A78(1,3)*A78(3,2))
          COFACTOR(2,2) = +(A78(1,1)*A78(3,3)-A78(1,3)*A78(3,1))
          COFACTOR(2,3) = -(A78(1,1)*A78(3,2)-A78(1,2)*A78(3,1))
          COFACTOR(3,1) = +(A78(1,2)*A78(2,3)-A78(1,3)*A78(2,2))
          COFACTOR(3,2) = -(A78(1,1)*A78(2,3)-A78(1,3)*A78(2,1))
          COFACTOR(3,3) = +(A78(1,1)*A78(2,2)-A78(1,2)*A78(2,1))
          AINV = TRANSPOSE(COFACTOR) / DET
          COEF=matmul(AINV,B78)
          C=COEF(1,1)
          B=COEF(2,1)
          A=COEF(3,1)
          x01=(-B)/(2*A)
c         print*, x01
c         pause
          if (x01<A78(3,2) .AND. x01 > A78(3,1) ) then
              DTmax=A*(x01**2)+B*x01+C
c              print*,DTmax
c              pause
			open (unit=1,file="C:\\Users\\s250171\\Desktop\\DTmax_2dwells.txt"
	1		,form="formatted",status="replace",action="write")
			write (unit=1,fmt=*),DTmax 
			close (unit=1)
              flag=1
          end if
      end if
	RETURN
      END
       SUBROUTINE UEDINC (INC1,INCSUB1)
       	IMPLICIT REAL *8 (A-H, O-Z)
	COMMON StoredINC, Tbound,DTmax,amin,Zeroflag,writeflag,tcure,Temp,TIME,
	1	Tempstore,tc1,a1,step
	include '../common/concom'
	include '../common/creeps'
           tc1=CPTIM-50
       if (amin<0.88) then
       a1=amin
       end if
      !print*,tc1,CPTIM,a1,amin,writeflag
      !pause
	if (amin>0.88.AND.writeflag==1) then
              tcure=((0.88*(CPTIM-tc1))+(amin*tc1)-(a1*CPTIM))/(amin-a1)
			open (unit=2,file="C:\\Users\\s250171\\Desktop\\tcure.txt"
	1		,form="formatted",status="replace",action="write")
			write (unit=2,fmt=*),tcure 
			close (unit=2)
			writeflag=0
	end if
	if (amin<0.88 .AND. CPTIM==25000) then
			tcure=80001
			open (unit=2,file="C:\\Users\\s250171\\Desktop\\tcure.txt"
	1		,form="formatted",status="replace",action="write")
			write (unit=2,fmt=*),tcure 
			close (unit=2)
	endif
         RETURN
       END