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,Tg
        real*8 ank1(3312,1000,1),ank2(3312,4,1000)
        real*8 usp1(3312,1000,1),usp2(3312,4,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--VOLUME FRACTION
      vf=0.65d0
C--INITIAL CONVECTION
      ako=0.017
C--CURE KINETICS       
      A1=19006.41379
      A2=22082.69835  
      Ad=6.7557E+18
      E1=72892.53275 
      E2=57827.01505
      Ed=137561.5952
      n1=1.965580509
      n2=1.525367094  
      mi=1.288659069
      l=0.435d0
      ba=0.452381424
      Tginf=206d0
      Tgo=-11d0
      H=465000.0
      R=8.3144
      g=0.029212404
      w=0.000466142
      Tg=Tgo+(Tginf-Tgo)*l*dt/(1-(1-l)*dt)
      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-dt)**n1+K2*dt**mi*(1-dt)**n2
      curedat(9)=H*(1-vf)
      volumi=1.0
	ak=dt+deltime*curerate
      !print*, 'CURE',dt, deltime, curerate, ak
      !pause
	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.63d0 !change
       l=0.435d0
       sigma=16.5d0
       Crcp=1.1d0
       Afcp=0.0014E3 !change
       Bfcp=0.841E3  !change
       Arcp=0.0025E3
       Brcp=1.80E3
       Drcp=-0.25E3
       p1= 2540!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
       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
      if (matus(1)==3) then
          spheat=0.840E3
      end if
      !print*, matus(1),spheat
      !pause
        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.63d0
         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= 1.03!!0.84 ->carbon fiber
         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=Btf!Alf*(t-273.15)+Blf
          
          K11=vf*Klf+(1-vf)*Kr
          !PRINT*, Kr,K11
          !pause
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
         !PRINT*, K11,K22
         !pause
         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
      
        
        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=120+273.15 !406.43    !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=58.5
      !T1=433 
      !T2=488
      !Dt=7400
      !r=1.5
      !r=r/60
      uemiss=0
      uqflux=0
      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,t1
      if (TIME<5400) then
          T=393
      else if (TIME<=1195+(5880) .AND. TIME>=(5400)) then
           T=0.0238*TIME+265.9!0.0232*TIME+404.536
      else if ((TIME >1195+(5880)) .AND. (TIME<=1260+(5880))) then
           T=432.77
       else if ((TIME >1260+(5880)).AND.(TIME<=1296.594+(5880))) then
           T=432.535
       else if ((TIME >1296.54+(5880)).AND.(TIME<=1320.219+(5880))) then
           T=432.24
       else if ((TIME >1320.219+(5880)).AND.(TIME<=1353.219+5880)) then
           T=431.87
       else if ((TIME >1353.219+(5880)).AND.(TIME<=1387.219+5880)) then
           T=431.51
       else if ((TIME >1387.219+5880).AND.(TIME<=1414.531+5880)) then
           T=431.445
       else if ((TIME >1414.531+(5880)).AND.(TIME<=1460.852+5880)) then
           T=431.495
       else if ((TIME >1460.852+5880).AND.(TIME<=1493.852+(5880))) then
           T=431.62
       else if ((TIME >1493.852+(5880)).AND.(TIME<=1506+(5880))) then
           T=431.78
       else if ((TIME>1506+(5880)) .AND. (TIME<=12480)) then
      t1=(0.0013*TIME)+426.8
      T=(2E-12)*(TIME*TIME*TIME)-9E-08*(TIME*TIME)+t1
       else if ((TIME>12480) .AND. (TIME<=13280)) then
          T=121.15+0.025*TIME
       else if ((TIME>13280) .AND. (TIME<=18600)) then
           T=453.15
      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       
      !print*, amin,conv,'hello'
       !   pause
	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),TN1(350),TN5(350),TN12(350)
          REAL*8  t,tINC,TN1F,TN5F,TN12F
          INTEGER tINCf,i
          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,TN1,TN5,TN12
		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=60
          NINC=310
          t=63.5 !input in minutes
          tINC=(t*60)/step !  increment of input
          tINCf=floor(tINC) ! to increment to kontinotero pros ta kato
          !print*, tINCf
          !pause
          if (NODEID==1) then
              TN1(INC)=Temp
              !print*, NODEID,TN1(INC)
              !pause
          end if
          if (NODEID==11) then
              TN5(INC)=Temp
          end if
          if (NODEID==33) then
              TN12(INC)=Temp
          end if
          if (INC==NINC) then
             
           do i=0,NINC
              if (tINCf==i) then
                  !print*, TN1(i),tINC,i,TN1(i+1)
                  !pause
                  
                  TN1F=TN1(i)+(tINC-i)*(TN1(i+1)-TN1(i))
                  !print*, TN1F,TN5F,TN12F
              ! pause
                open (unit=1,file="C:\\Users\s250171\\Desktop\\TN1.txt"
	1		    ,form="formatted",status="replace",action="write")
			    write (unit=1,fmt=*),TN1F
			    close (unit=1)
                  TN5F=TN5(i)+(tINC-i)*(TN5(i+1)-TN5(i))
                open (unit=1,file="C:\\Users\s250171\\Desktop\\TN5.txt"
	1		,form="formatted",status="replace",action="write")
			write (unit=1,fmt=*),TN5F
			close (unit=1)
                  TN12F=TN12(i)+(tINC-i)*(TN12(i+1)-TN12(i))
               open (unit=1,file="C:\\Users\s250171\\Desktop\\TN12.txt"
	1		 ,form="formatted",status="replace",action="write")
			 write (unit=1,fmt=*),TN12F
			 close (unit=1)
               !print*, TN1F,TN5F,TN12F
               !pause
              end if
           end do
          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-60
       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
              !print*, tcure
              !pause
	end if
	if (amin<0.88 .AND. CPTIM==18600) then
              print*, amin
              pause
			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