cforc3.F File Reference

#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "scr02_c.inc"
#include "scr17_c.inc"
#include "scr18_c.inc"
#include "param_c.inc"
#include "com_xfem1.inc"
#include "parit_c.inc"
#include "timeri_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine

cforc3 (timers, elbuf_str, jft, jlt, pm, ixc, x, f, m, v, vr, failwave, nvc, mtn, geo, tf, npf, bufmat, partsav, dt2t, neltst, ityptst, stifn, stifr, fsky, iadc, itab, d, dr, tani, offset, eani, f11, f12, f13, f14, f21, f22, f23, f24, f31, f32, f33, f34, m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, indxof, ipartc, thke, group_param, mat_elem, nel, istrain, ihbe, ithk, iofc, ipla, nft, ismstr, npt, kfts, fzero, igeo, ipm, ifailure, itask, jthe, temp, fthe, fthesky, iexpan, gresav, grth, xedge4n, igrth, msc, dmelc, jsms, table, iparg, ixfem, knod2elc, sensors, elcutc, inod_crk, iel_crk, ibordnode, nodenr, iadc_crk, nodedge, crknodiad, condn, condnsky, stack, isubstack, xfem_str, ig, crkedge, drape_sh4n, ipri, nloc_dmg, indx_drape, igre, jtur, output, dt, snpc, stf, glob_therm, userl_avail, maxfunc, sbufmat, ipart)

Function/Subroutine Documentation

cforc3()

subroutine cforc3	(	type(timer_), intent(inout)	timers,
		type (elbuf_struct_), target	elbuf_str,
		integer	jft,
		integer	jlt,
			pm,
		integer, dimension(nixc,*)	ixc,
			x,
			f,
			m,
			v,
			vr,
		type (failwave_str_)	failwave,
		integer	nvc,
		integer	mtn,
			geo,
			tf,
		integer, dimension(*)	npf,
			bufmat,
			partsav,
			dt2t,
		integer	neltst,
		integer	ityptst,
			stifn,
			stifr,
			fsky,
		integer, dimension(4,*)	iadc,
		integer, dimension(*)	itab,
			d,
			dr,
			tani,
		integer	offset,
			eani,
			f11,
			f12,
			f13,
			f14,
			f21,
			f22,
			f23,
			f24,
			f31,
			f32,
			f33,
			f34,
			m11,
			m12,
			m13,
			m14,
			m21,
			m22,
			m23,
			m24,
			m31,
			m32,
			m33,
			m34,
		integer, dimension(mvsiz)	indxof,
		integer, dimension(*)	ipartc,
			thke,
		type (group_param_), intent(in)	group_param,
		type (mat_elem_), intent(inout)	mat_elem,
		integer	nel,
		integer	istrain,
		integer	ihbe,
		integer	ithk,
		integer	iofc,
		integer	ipla,
		integer	nft,
		integer	ismstr,
		integer	npt,
		integer	kfts,
			fzero,
		integer, dimension(npropgi,*)	igeo,
		integer, dimension(npropmi,*)	ipm,
		integer	ifailure,
		integer	itask,
		integer	jthe,
			temp,
			fthe,
			fthesky,
		integer	iexpan,
			gresav,
		integer, dimension(*)	grth,
		integer, dimension(4,*)	xedge4n,
		integer, dimension(*)	igrth,
			msc,
			dmelc,
		integer	jsms,
		type (ttable), dimension(*)	table,
		integer, dimension(*)	iparg,
		integer	ixfem,
		integer, dimension(*)	knod2elc,
		type (sensors_), intent(inout)	sensors,
		integer, dimension(2,*)	elcutc,
		integer, dimension(*)	inod_crk,
		integer, dimension(*)	iel_crk,
		integer, dimension(*)	ibordnode,
		integer, dimension(*)	nodenr,
		integer, dimension(4,*)	iadc_crk,
		integer, dimension(2,*)	nodedge,
		integer, dimension(*)	crknodiad,
			condn,
			condnsky,
		type (stack_ply)	stack,
		integer	isubstack,
		type (elbuf_struct_), dimension(ngroup,nxel)	xfem_str,
		integer, intent(in)	ig,
		type (xfem_edge_), dimension(*)	crkedge,
		type (drape_), dimension(numelc_drape)	drape_sh4n,
		integer	ipri,
		type (nlocal_str_), target	nloc_dmg,
		integer, dimension(scdrape)	indx_drape,
		integer, intent(in)	igre,
		integer, intent(in)	jtur,
		type (output_), intent(inout)	output,
		type (dt_), intent(in)	dt,
		integer, intent(in)	snpc,
		integer, intent(in)	stf,
		type (glob_therm_), intent(inout)	glob_therm,
		integer, intent(in)	userl_avail,
		integer, intent(in)	maxfunc,
		integer, intent(in)	sbufmat,
		integer, dimension(lipart1,npart), intent(in)	ipart )

Parameters

[in,out]

output

output structure

Definition at line 82 of file cforc3.F.

C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE timer_mod
      USE table_mod
      USE mat_elem_mod
      USE crackxfem_mod
      USE stack_mod
      USE failwave_mod
      USE nlocal_reg_mod
      USE drape_mod
      USE sensor_mod
      USE output_mod
      USE elbufdef_mod
      USE dt_mod
      use glob_therm_mod
      use dttherm_mod
      use element_mod , only : nixc
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "com04_c.inc"
#include      "scr02_c.inc"
#include      "scr17_c.inc"
#include      "scr18_c.inc"
#include      "param_c.inc"
#include      "com_xfem1.inc"
#include      "parit_c.inc"
#include      "timeri_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      TYPE(TIMER_), INTENT(INOUT) :: TIMERS
      INTEGER, INTENT(IN) :: USERL_AVAIL ! Flag for User libraries availability
      INTEGER,INTENT(IN)    :: MAXFUNC   ! Maximum number of functions
      INTEGER, INTENT(IN) :: SBUFMAT     ! size of Bufmat
      INTEGER, INTENT(IN) :: STF         ! size of STF
      INTEGER, INTENT(IN) :: SNPC        ! size of SNPC
      INTEGER, INTENT(IN) :: JTUR
      INTEGER, INTENT(IN) :: IGRE
      INTEGER, INTENT(IN) :: IG
      INTEGER JFT,JLT,NVC,MTN,NELTST,ITYPTST,OFFSET,
     .   NEL,ISTRAIN,IHBE,ITASK,JTHE,JSMS,
     .   ITHK,IOFC,IPLA,NFT,ISMSTR,NPT,KFTS,IFAILURE,IEXPAN,
     .   ISUBSTACK,IPRI
      INTEGER NPF(*),IXC(NIXC,*),IADC(4,*),INDXOF(MVSIZ),IPARTC(*),
     .   IGEO(NPROPGI,*),IPM(NPROPMI,*),XEDGE4N(4,*),ITAB(*),
     .   GRTH(*),IGRTH(*),IPARG(*),IXFEM,KNOD2ELC(*),
     .   ELCUTC(2,*),INOD_CRK(*),IEL_CRK(*),IBORDNODE(*),
     .   NODENR(*),IADC_CRK(4,*),NODEDGE(2,*),CRKNODIAD(*),INDX_DRAPE(SCDRAPE)
      my_real
     .   f11(mvsiz),f12(mvsiz),f13(mvsiz),f14(mvsiz),
     .   f21(mvsiz),f22(mvsiz),f23(mvsiz),f24(mvsiz),
     .   f31(mvsiz),f32(mvsiz),f33(mvsiz),f34(mvsiz),
     .   m11(mvsiz),m12(mvsiz),m13(mvsiz),m14(mvsiz),
     .   m21(mvsiz),m22(mvsiz),m23(mvsiz),m24(mvsiz),
     .   m31(mvsiz),m32(mvsiz),m33(mvsiz),m34(mvsiz),
     .   pm(npropm,*),x(*),f(*),m(*),v(*),vr(*),
     .   geo(npropg,*),tf(*), bufmat(*),partsav(*),stifn(*),stifr(*),
     .   fsky(*),d(*),dr(*),tani(6,*),eani(*), thke(*),dt2t,
     .   fzero(3,4,*),temp(*), fthe(*),fthesky(*),gresav(*),
     .   msc(*), dmelc(*),condn(*),condnsky(*)
C
      TYPE (TTABLE) TABLE(*)
      TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_STR
      TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP,NXEL) :: XFEM_STR  ! take xfem_str
      TYPE (XFEM_EDGE_)   , DIMENSION(*) :: CRKEDGE
      TYPE (FAILWAVE_STR_) :: FAILWAVE 
      TYPE (NLOCAL_STR_), TARGET :: NLOC_DMG
      TYPE (DRAPE_) , DIMENSION(NUMELC_DRAPE) :: DRAPE_SH4N
      TYPE (OUTPUT_)     ,INTENT(INOUT) :: OUTPUT !< output structure
      TYPE (SENSORS_)    ,INTENT(INOUT) :: SENSORS
      TYPE (MAT_ELEM_)   ,INTENT(INOUT) :: MAT_ELEM
      TYPE (GROUP_PARAM_),INTENT(IN)    :: GROUP_PARAM
      TYPE (DT_), INTENT(IN) :: DT
      type (glob_therm_) ,intent(inout) :: glob_therm
    
      INTEGER, DIMENSION(LIPART1,NPART), INTENT(IN)     :: IPART
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER MAT(MVSIZ),PID(MVSIZ),NDT(MVSIZ),NGL(MVSIZ),FWAVE_EL(NEL),
     . I,J,IGTYP,ICSEN,IUN,NPG,ILAY,NLAY,IXEL,IXLAY,NXLAY,
     .   IBID,NG,IR,IS,L_DIRA,L_DIRB,J1,J2,IGMAT,NPTTOT,IREP,IFAILWAVE,IDRAPE,
     .   NPTT,IT,ACTIFXFEM,SEDRAPE,NUMEL_DRAPE,IMAT
      my_real 
     .   sti(mvsiz),stir(mvsiz),sigy(mvsiz),rho(mvsiz),
     .   x2(mvsiz),x3(mvsiz),x4(mvsiz),y2(mvsiz),y3(mvsiz),
     .   y4(mvsiz),z2(mvsiz),ssp(mvsiz),viscmx(mvsiz),
     .   vx1(mvsiz), vx2(mvsiz), vx3(mvsiz), vx4(mvsiz),
     .   vy1(mvsiz), vy2(mvsiz), vy3(mvsiz), vy4(mvsiz),
     .   vz1(mvsiz), vz2(mvsiz), vz3(mvsiz), vz4(mvsiz),
     .
     . thk02(mvsiz),
     .   x1g(mvsiz), x2g(mvsiz), x3g(mvsiz), x4g(mvsiz),
     .   y1g(mvsiz), y2g(mvsiz), y3g(mvsiz), y4g(mvsiz),
     .   z1g(mvsiz), z2g(mvsiz), z3g(mvsiz), z4g(mvsiz),
     .   e1x(mvsiz), e1y(mvsiz), e1z(mvsiz), e2x(mvsiz),
     .   e2y(mvsiz), e2z(mvsiz), e3x(mvsiz), e3y(mvsiz),e3z(mvsiz)
      my_real 
     .   exx(mvsiz), eyy(mvsiz), exy(mvsiz), exz(mvsiz),eyz(mvsiz),
     .   kxx(mvsiz), kyy(mvsiz), kxy(mvsiz), px1(mvsiz),
     .   px2(mvsiz), py1(mvsiz), py2(mvsiz), thk0(mvsiz),
     .   off(mvsiz), nu(mvsiz) , shf(mvsiz), area(mvsiz),
     .   g(mvsiz)  , ym(mvsiz) , a11(mvsiz), a12(mvsiz),
     .   vl1(mvsiz,3),vl2(mvsiz,3),vl3(mvsiz,3),vl4(mvsiz,3),
     .   vrl1(mvsiz,3),vrl2(mvsiz,3),vrl3(mvsiz,3),vrl4(mvsiz,3), 
     .   dt1c(mvsiz), dt2c(mvsiz),
     .   aldt(mvsiz),alpe(mvsiz),vhx(mvsiz),vhy(mvsiz)
      my_real 
     .   h1(mvsiz),  h2(mvsiz),  h3(mvsiz),  vol0(mvsiz),vol00(mvsiz),
     .   h11(mvsiz), h12(mvsiz), h13(mvsiz), h14(mvsiz),  
     .   h21(mvsiz), h22(mvsiz), h23(mvsiz), h24(mvsiz),  
     .   h31(mvsiz), h32(mvsiz), h33(mvsiz), h34(mvsiz),  
     .   b11(mvsiz), b12(mvsiz), b13(mvsiz), b14(mvsiz),  
     .   b21(mvsiz), b22(mvsiz), b23(mvsiz), b24(mvsiz),  
     .   rx1(mvsiz), rx2(mvsiz), rx3(mvsiz), rx4(mvsiz), 
     .   ry1(mvsiz), ry2(mvsiz), ry3(mvsiz), ry4(mvsiz),
     .   zcfac(mvsiz,2),gs(mvsiz),
     .   srh1(mvsiz),srh2(mvsiz),srh3(mvsiz),a_i(mvsiz),
     .   die(mvsiz),tempel(mvsiz),them(mvsiz,4),bid,
     .   ux1(mvsiz),ux2(mvsiz),ux3(mvsiz),ux4(mvsiz),
     .   uy1(mvsiz),uy2(mvsiz),uy3(mvsiz),uy4(mvsiz),
     .
     . dtcsub(mvsiz),conde(mvsiz),a11r(mvsiz)
 
      my_real 
     .  areat(mvsiz),px1t(mvsiz),px2t(mvsiz),py1t(mvsiz),py2t(mvsiz),
     .  f_def(mvsiz,8), u13x(mvsiz),u24x(mvsiz),u13y(mvsiz),u24y(mvsiz),
     .  wxy(mvsiz)
C
      INTEGER, DIMENSION(NEL) :: OFFLY
      INTEGER, ALLOCATABLE, DIMENSION(:) :: ELCRKINI
      my_real, ALLOCATABLE, DIMENSION(:) :: dira,dirb,dir1_crk,dir2_crk
      my_real, DIMENSION(:) ,POINTER     :: dir_a, dir_b,crkdir,crklen,dadv
      TARGET :: dira,dirb
      my_real :: dt1,dtinv,asrate,eps_m2,eps_k2
      my_real, dimension(nel)   :: epsd_pg
      my_real, dimension(mvsiz) :: fheat        ! variables for heat transfer
      my_real, dimension(mvsiz) :: ssp_eq        
!
C--- Variables for non-local
      INTEGER :: NDDL, K, INOD(4),NC1(MVSIZ), NC2(MVSIZ), NC3(MVSIZ), NC4(MVSIZ),
     .           IPOS(4), L_NLOC, INLOC
      my_real, DIMENSION(:,:), ALLOCATABLE :: var_reg
      my_real, DIMENSION(:), POINTER :: dnl
C-----
 
      TYPE(G_BUFEL_) ,POINTER :: GBUF
 
      TYPE(L_BUFEL_DIR_) ,POINTER :: LBUF_DIR
      TYPE (STACK_PLY) :: STACK
      INTEGER SDIR_A    ! Size of DIR_A
      INTEGER SDIR_B    ! Size of DIR_B
C=======================================================================
      gbuf => elbuf_str%GBUF
      idrape = elbuf_str%IDRAPE
C-----
      iun = 1
      ibid = 0
      bid = zero
      nlay = elbuf_str%NLAY
      sedrape = scdrape
      numel_drape = numelc_drape
      tempel(:) = zero
      fheat(: ) = zero
 
cc      NPT  = MAX(NLAY,NPTT) --> set to = IPARG(6) , keeping it original
C                                 to allow for NPT = 0 (global LAW_3)
      npg  = 1
      ir   = 1
      is   = 1
      ng   = 1
      ixel  = 0
      ixlay = 0
      irep = iparg(35)
      inloc = iparg(78)
      actifxfem = iparg(70)
C
      npttot  = 0
      DO ilay=1,nlay
        npttot = npttot + elbuf_str%BUFLY(ilay)%NPTT
      ENDDO
!-------------------------------------------      
! Array for non-local variable
      nddl = npttot
      ALLOCATE(var_reg(nel,nddl)) 
!-------------------------------------------
      IF (npt == 0) npttot = npt  !  compatibility with global integration
c--------------------------------------------
c     Front wave
c--------------------------------------------
      ifailwave = iparg(79)
      IF (ifailwave > 0) THEN
        fwave_el(:) = zero
        offly(:) = elbuf_str%BUFLY(1)%OFF(:)
        DO i=2,nlay
          DO j=1,nel
            offly(j) = max(offly(j), elbuf_str%BUFLY(i)%OFF(j))
          ENDDO
        ENDDO        
        dadv => gbuf%DMG
        CALL set_failwave_sh4n(failwave    ,fwave_el ,dadv     ,
     .       nel      ,ixc      ,itab      ,ngl      ,offly    )
      ENDIF
c-------------------------------------
      l_dira = elbuf_str%BUFLY(1)%LY_DIRA
      l_dirb = elbuf_str%BUFLY(1)%LY_DIRB
      igtyp= igeo(11,ixc(6,1))
      IF(idrape > 0 .AND. (igtyp == 51 .OR. igtyp == 52) ) THEN
        ALLOCATE(dira(npttot*nel*l_dira))
        ALLOCATE(dirb(npttot*nel*l_dirb))
        IF (l_dira == 0) THEN
            CONTINUE
        ELSEIF (irep == 0) THEN
           npttot = 0
           DO ilay=1,nlay
              nptt = elbuf_str%BUFLY(ilay)%NPTT
              DO it=1,nptt
                 j = npttot + it
                 lbuf_dir =>  elbuf_str%BUFLY(ilay)%LBUF_DIR(it)
                 j1 = 1+(j-1)*l_dira*nel
                 j2 = j*l_dira*nel
                 dira(j1:j2) = lbuf_dir%DIRA(1:nel*l_dira)
              ENDDO
              npttot = npttot + nptt
            ENDDO 
        ENDIF
        sdir_a=npttot*nel*l_dira
        sdir_b=npttot*nel*l_dirb
        dir_a => dira(1:npttot*nel*l_dira)
        dir_b => dirb(1:npttot*nel*l_dirb)
      ELSE ! idrape
        sdir_a=nlay*nel*l_dira
        sdir_b=nlay*nel*l_dirb
        ALLOCATE(dira(nlay*nel*l_dira))
        ALLOCATE(dirb(nlay*nel*l_dirb))
        dira=zero
        dirb=zero
        IF (l_dira == 0) THEN
          CONTINUE
        ELSEIF (irep == 0) THEN
           DO j=1,nlay
              j1 = 1+(j-1)*l_dira*nel
              j2 = j*l_dira*nel
              dira(j1:j2) = elbuf_str%BUFLY(j)%DIRA(1:nel*l_dira)
           ENDDO
         ENDIF
         sdir_a=nlay*nel*l_dira
         sdir_b=nlay*nel*l_dirb
         dir_a => dira(1:nlay*nel*l_dira)
         dir_b => dirb(1:nlay*nel*l_dirb)
      ENDIF ! IDRAPE    
c-------------------------------------
      nxlay = nlay
C---
      IF (ixfem > 0) THEN
        ALLOCATE(elcrkini(nxlaymax*nel))
        ALLOCATE(dir1_crk(nxlaymax*nel))
        ALLOCATE(dir2_crk(nxlaymax*nel))
        dir1_crk = zero
        dir2_crk = zero
        elcrkini = 0
        IF (nlevset > 0) THEN
          CALL precrklay(jft     ,jlt    ,nft     ,nxlay  ,elcrkini,
     .                   iel_crk,inod_crk,nodenr  ,crkedge,xedge4n )
        ENDIF
      ELSE
        ALLOCATE(elcrkini(0))
        ALLOCATE(dir1_crk(0))
        ALLOCATE(dir2_crk(0))
      ENDIF
C----------------------
C       X1G-> X1,X2G-> X2,X3G-> X3,X4G-> X4,
C       Y1G-> Y1,Y2G-> Y2,Y3G-> Y3,Y4G-> Y4,
C       Z1G-> Z1,Z2G-> Z2,Z3G-> Z3,Z4G-> Z4
        CALL ccoor3(jft  ,jlt  ,x    ,ixc  ,geo  ,gbuf%OFF,
     2              off  ,sigy ,pid ,v ,vr ,vl1,vl2,vl3,vl4,
     3              vrl1,vrl2,vrl3,vrl4,mat,dt1c,thke ,thk0 ,ngl,
     4              x1g  ,x2g  ,x3g  ,x4g  ,y1g  ,y2g    ,
     5              y3g  ,y4g  ,z1g  ,z2g  ,z3g  ,z4g    )
C
      icsen= igeo(3,pid(1))
      igtyp= igeo(11,pid(1))
!!      IF(IGTYP == 11 .AND. IDTMINS > 0)IGEO(98,PID(1)) = 1
       igmat= igeo(98,pid(1))
C
      IF (ishfram == 1) THEN
C       non convected frame (version 3) 
        CALL cevec3(elbuf_str ,dir_a,dir_b,
     1    jft ,jlt ,x1g ,x2g ,x3g ,x4g ,y1g ,y2g ,y3g ,y4g ,
     2    z1g ,z2g ,z3g ,z4g ,e1x ,e1y ,e1z ,e2x ,e2y ,e2z , 
     3    e3x ,e3y ,e3z ,irep,nlay,nel)
      ELSE
C       repere orthogonalise convecte
        CALL cnvec3(elbuf_str ,dir_a   ,dir_b   ,
     1       jft     ,jlt     ,irep    ,igtyp   ,nlay    ,
     2       x1g     ,x2g     ,x3g     ,x4g     ,y1g     ,y2g     ,       
     3       y3g     ,y4g     ,z1g     ,z2g     ,z3g     ,z4g     ,       
     4       e1x     ,e1y     ,e1z     ,e2x     ,e2y     ,e2z     ,       
     5       e3x     ,e3y     ,e3z     ,nel     )
      ENDIF
C
      IF (ismstr /= 3)THEN
       CALL cderi3(
     1   jft,       jlt,       gbuf%SMSTR,gbuf%OFF,
     2   sti,       stir,      area,      px1,
     3   px2,       py1,       py2,       x2,
     4   x3,        x4,        y2,        y3,
     5   y4,        z2,        x1g,       x2g,
     6   x3g,       x4g,       y1g,       y2g,
     7   y3g,       y4g,       z1g,       z2g,
     8   z3g,       z4g,       e1x,       e1y,
     9   e1z,       e2x,       e2y,       e2z,
     a   e3x,       e3y,       e3z,       vhx,
     b   vhy,       a_i,       ux1,       ux2,
     c   ux3,       ux4,       uy1,       uy2,
     d   uy3,       uy4,       nel,       ismstr)
      ELSE
       CALL cpxpy3(jft       ,jlt     ,pm      ,sti     ,stir    ,
     2             gbuf%SMSTR,px1     ,px2     ,py1     ,py2     ,
     3             ixc       ,area    ,x2      ,x3      ,x4      , 
     4             y2        ,y3      ,y4      ,z2      ,thk0    ,  
     5             mat       ,nel     )
      ENDIF
C
      CALL ccoef3(jft     ,jlt    ,pm     ,mat      ,geo     ,
     2            pid     ,off    ,area   ,sti      ,stir    ,
     3            shf     ,thk0   ,thk02  ,nu       ,
     4            g       ,ym     ,a11    ,a12      ,gbuf%THK,
     5            ssp     ,rho    ,h1     ,h2       ,h3      ,
     6            vol0    ,vol00  ,alpe   ,gs       ,mtn     ,
     7            ithk    ,ismstr  ,npttot ,kfts   ,
     8            srh1   ,srh2   ,srh3     ,igeo    ,
     9            a11r    ,isubstack      ,stack%PM )
C
      IF ((ismstr /= 3 .AND. idt1sh == 0) .OR.
     .     idtmin(3) /= 0 .OR. igtyp == 16 .OR. glob_therm%IDT_THERM == 1) THEN
C      CDTLEN must be always called for tissue property/law for time step correction
       CALL cdlen3(jft     ,jlt     ,pm      ,off     ,area,
     2             x2      ,x3      ,x4      ,y2      ,y3  ,
     3             y4      ,aldt    ,mat     ,geo     ,pid ,
     4             ihbe    )
      ENDIF
      CALL cdefo3(jft ,jlt,vl1,vl2,vl3,vl4,dt1c,px1 ,px2 ,py1,py2,area,
     2            exx ,eyy,exy,exz ,eyz ,x2  ,x3 ,x4 ,y2  ,y3,
     3            y4  ,z2 ,vx1,vx2 ,vx3 ,vx4 ,vy1,vy2,vy3 ,vy4 ,
     4            vz1 ,vz2,vz3,vz4 ,e1x ,e1y ,e1z,e2x,e2y ,e2z ,
     5            e3x ,e3y,e3z,ihbe)
      CALL ccurv3(jft ,jlt ,vrl1,vrl2,vrl3,vrl4 ,px1 
     1            ,px2 ,py1 ,py2 ,area,
     2            rx1 ,rx2 ,rx3 ,rx4 ,ry1 ,ry2 ,ry3 ,ry4 ,
     3            e1x ,e1y ,e1z ,e2x ,e2y ,e2z ,e3x ,e3y ,
     4            e3z ,kxx ,kyy ,kxy ,exz ,eyz ,
     5            wxy ,ismstr)
      IF (ismstr == 10 ) THEN
        CALL ccoort3(jft    ,jlt    ,x     ,ixc   ,gbuf%OFF   ,
     1               dr     ,px1t  ,px2t  ,py1t   ,py2t       ,
     2               e1x    ,e1y   ,e1z   ,e2x    ,e2y        ,
     3               e2z    ,e3x   ,e3y   ,e3z    ,areat      , 
     4               u13x   ,u24x  ,u13y  ,u24y   ,gbuf%SMSTR ,
     5               nel    )
C  
        CALL cdefot3(jft ,jlt  ,px1t  ,px2t  ,py1t   ,
     2              py2t ,u13x ,u24x  ,u13y  ,u24y   ,
     3              f_def )
      END IF !(ISMSTR == 10 ) THEN
      CALL cstra3(jft    ,jlt       ,gbuf%STRA,shf    ,area , 
     2            exx    ,eyy       ,exy      ,exz    ,eyz  ,   
     3            kxx    ,kyy       ,kxy      ,dt1c   ,tani , 
     4            gbuf%FOR,gbuf%MOM ,ismstr ,mtn  ,
     6            ihbe   ,nft       ,istrain  ,ux1    ,ux2  ,
     7            ux3    ,ux4       ,uy1      ,uy2    ,uy3  ,
     8            uy4    ,px1       ,px2      ,py1    ,py2  ,
     9            wxy    ,gbuf%STRW ,f_def    ,nel    )
c-------------------------------------------
      ! transfer nodal temperature to gauss point from /heat/mat
      IF (jthe /= 0) THEN      
        CALL tempcg(numnod,nel   ,ixc  ,temp ,tempel )
      ENDIF
c-------------------------------------------
c    COMPUTE Regularized non local variable in Gauss point
c-------------------------------------------     
      IF (inloc > 0) THEN
        l_nloc = nloc_dmg%L_NLOC
        dnl  => nloc_dmg%DNL(1:l_nloc) ! DNL = non local variable increment
        DO i=jft,jlt
          nc1(i)  = ixc(2,i)
          nc2(i)  = ixc(3,i)
          nc3(i)  = ixc(4,i)
          nc4(i)  = ixc(5,i)
        ENDDO
        DO k = 1,nddl
#include "vectorize.inc"
          DO i=jft,jlt
            inod(1) = nloc_dmg%IDXI(nc1(i))
            inod(2) = nloc_dmg%IDXI(nc2(i))
            inod(3) = nloc_dmg%IDXI(nc3(i))
            inod(4) = nloc_dmg%IDXI(nc4(i)) 
            ipos(1) = nloc_dmg%POSI(inod(1)) 
            ipos(2) = nloc_dmg%POSI(inod(2)) 
            ipos(3) = nloc_dmg%POSI(inod(3)) 
            ipos(4) = nloc_dmg%POSI(inod(4))
            var_reg(i,k) = fourth*(dnl(ipos(1)+k-1) + dnl(ipos(2)+k-1) 
     .                          + dnl(ipos(3)+k-1) + dnl(ipos(4)+k-1))
          ENDDO
        ENDDO
      ENDIF     
!-------------------------------------------------------------------------------
!       global element strain rate (shell energy equivalent)
!-----------------------------------------------------------
!       e = 1/t integ[1/2 e (eps_m + k z)^2 dz ]
!       e = 1/2 e eps_eq^2
!       eps_eq = sqrt[ eps_m^2 + 1/12 k^2t^2 ]
!-------------------------------------------------------------------------------
      dt1 = dt1c(1)
      dtinv = dt1 / max(dt1**2,em20)  ! inverse of dt
      asrate = one   ! to be changed for default value
#include "vectorize.inc"
      do i = 1,nel
        eps_k2 = (kxx(i)**2+kyy(i)**2+kxx(i)*kyy(i)+fourth*kxy(i)**2)
     .         * one_over_9*gbuf%thk(i)**2
        eps_m2 = four_over_3*(exx(i)**2+eyy(i)**2+exx(i)*eyy(i) + fourth*exy(i)**2)
        epsd_pg(i) = sqrt(eps_k2 + eps_m2)*dtinv
      end do
      gbuf%epsd(1:nel) = asrate * epsd_pg(1:nel) + (one - asrate) * gbuf%epsd(1:nel)      
!-------------------------------------------------------------------------------
      IF ((imon_mat==1).AND. itask == 0)CALL startime(timers,35)
C-----------------------------
        CALL cmain3(timers,
     1  elbuf_str   ,jft         ,jlt         ,nft        ,iparg      ,
     2  nel         ,mtn         ,ipla        ,ithk       ,group_param,
     3  pm          ,geo         ,npf         ,tf         ,bufmat     ,
     4  ssp         ,rho         ,viscmx      ,dt1c       ,sigy       ,
     5  area        ,exx         ,eyy         ,exy        ,exz        ,
     6  eyz         ,kxx         ,kyy         ,kxy        ,nu         ,
     7  off         ,thk0        ,mat         ,pid        ,mat_elem   ,
     8  gbuf%FOR    ,gbuf%MOM    ,gbuf%STRA   ,failwave   ,fwave_el   ,
     9  gbuf%THK    ,gbuf%EINT   ,iofc        ,
     a  g           ,a11         ,a12         ,vol0       ,indxof     ,
     b  ngl         ,zcfac       ,shf         ,gs         ,epsd_pg    ,
     c  kfts        ,ihbe        ,alpe        ,
     d  dir_a       ,dir_b       ,igeo        ,
     e  ipm         ,ifailure    ,npg         ,fheat      ,
     f  tempel      ,die         ,jthe        ,iexpan     ,gbuf%TEMP  ,
     g  ibid        ,bid         ,
     h  bid         ,bid         ,bid         ,bid        ,bid        ,
     i  bid         ,bid         ,bid         ,e1x        ,e1y        ,
     j  e1z         ,e2x         ,e2y         ,e2z        ,e3x        ,
     k  e3y         ,e3z         ,ng          ,table      ,ixfem      ,
     l  bid         ,sensors     ,bid         ,elcrkini   ,
     m  dir1_crk    ,dir2_crk    ,aldt        ,glob_therm%IDT_THERM ,glob_therm%THEACCFACT,
     n  ismstr      ,ir          ,is          ,nlay       ,npt        ,
     o  ixlay       ,ixel        ,isubstack   ,stack      ,
     p  f_def       ,itask       ,drape_sh4n  ,var_reg    ,nloc_dmg   ,
     r  indx_drape  ,thke        ,sedrape     ,numel_drape  ,dt        ,
     q  ncycle      ,snpc        ,stf         ,nxlaymax    ,idel7nok   , 
     s  userl_avail ,maxfunc     ,npttot      ,sbufmat     ,sdir_a     ,
     t  sdir_b      ,gbuf%FOR_G  ,ssp_eq ,
     x    ipart     ,lipart1     ,ipartc     )
C-----------------------------
      IF ((imon_mat==1).AND. itask == 0)CALL stoptime(timers,35)
C-----------------------------
C     FORCES ANTI-SABLIER
C-----------------------------
      DO i=jft,jlt
        viscmx(i) = sqrt(one + viscmx(i)*viscmx(i)) - viscmx(i)
      ENDDO
      IF (nodadt /= 0 .AND. mtn ==58 ) 
     .            CALL cssp2a11(pm,mat(jft),ssp ,a11  ,jlt   ,1 )
c---------------------------------------------------
      IF (iabs(npttot) == 1) THEN
        CALL mhvis3(jft    ,jlt    ,pm     ,gbuf%THK,gbuf%HOURG,
     2              off    ,px1    ,px2    ,py1     ,py2       ,
     3              ixc    ,dt1c   ,ssp    ,rho     ,sti       ,
     4              eani   ,geo    ,pid    ,stir    ,mat       ,
     5              thk0   ,viscmx ,alpe   ,ipartc  ,partsav   ,
     6              ihbe   ,nft    ,ismstr  ,rx1       ,
     7              rx2    ,rx3    ,rx4    ,ry1     ,ry2       ,
     8              ry3    ,ry4    ,vx1    ,vx2     ,vx3       ,
     9              vx4    ,vy1    ,vy2    ,vy3     ,vy4       ,
     a              vz1    ,vz2    ,vz3    ,vz4     ,b11       ,
     b              b12    ,b13    ,b14    ,b21     ,b22       ,
     c              b23    ,b24    ,area   ,ym      ,nu        ,
     d              vhx    ,vhy    ,h11    ,h12     ,h13       ,
     e              h14    ,h21    ,h22    ,h23     ,h24       ,
     f              h31    ,h32    ,h33    ,h34     ,h1        ,
     g              h2     ,igeo   ,nel    ,mtn     ,a11      )
      ELSEIF(ihbe == 2)THEN
        CALL chsti3(jft    ,jlt    ,gbuf%THK,gbuf%HOURG,off    ,px1 ,
     2              px2    ,py1    ,py2     ,sigy      ,ixc    ,dt1c,
     3              ssp    ,rho    ,sti     ,z2        ,eani   ,stir,
     4              shf    ,thk0   ,thk02   ,viscmx    ,g      ,a11 ,
     5   h1  ,h2  ,h3 ,ym   ,nu  , alpe   , 
     6   vhx ,vhy ,vx1 ,vx2 ,vx3 ,vx4 ,vy1 ,
     7   vy2 ,vy3 ,vy4,vz1  ,vz2  ,vz3 ,vz4 ,area   ,
     8   h11 ,h12 ,h13 ,h21 ,h22 ,h23  ,h31  ,h32  ,h33  ,
     9   b11 ,b12 ,b13 ,b14 ,b21 ,b22  ,b23  ,b24  ,
     a   rx1 ,rx2 ,rx3 ,rx4 ,ry1 ,ry2  ,ry3  ,ry4, 
     b   ipartc,partsav,
     c   ihbe   ,nft     ,ismstr,srh3,igtyp  ,
     d   igmat  ,a11r  ,nel)
      ELSE
        CALL chvis3(
     2   jft ,jlt ,gbuf%THK,gbuf%HOURG,off ,px1 ,px2 ,py1 ,py2 ,
     3   ixc ,dt1c,ssp,rho  ,sti  ,vx1 ,vx2 ,vx3 ,vx4 ,vy1 ,
     4   vy2 ,vy3 ,vy4,vz1  ,vz2  ,vz3 ,vz4 ,area,thk0,vhx ,
     5   vhy ,shf ,z2 ,eani ,stir,viscmx,g  ,a11 ,
     6   h1  ,h2  ,h3 ,ym   ,nu  ,thk02,alpe,h11 ,
     7   h12 ,h13 ,h21 ,h22 ,h23 ,h31  ,h32 ,h33 ,
     8   b11 ,b12 ,b13 ,b14 ,b21 ,b22  ,b23 ,b24 ,
     9   rx1 ,rx2 ,rx3 ,rx4 ,ry1 ,ry2  ,ry3 ,ry4 , 
     a   ipartc,partsav,
     b   ihbe   ,nft     ,ismstr,kfts   ,
     c    srh1, srh2, srh3   ,igtyp ,
     d   igmat   ,a11r  ,nel)
      ENDIF
C--------------------------
C     BALANCES BY MATERIAL
C--------------------------
c      IFLAG=MOD(NCYCLE,NCPRI)
      IF (ipri>0) THEN
         CALL cbilan(
     1   jft,        jlt,        pm,         v,
     2   ixc,        gbuf%THK,   gbuf%EINT,  partsav,
     3   area,       mat,        ipartc,     x,
     4   vr,         vol0,       vol00,      thk0,
     5   thk02,      1,          off,        nft,
     6   gresav,     grth,       igrth,      vl1,
     7   vl2,        vl3,        vl4,        vrl1,
     8   vrl2,       vrl3,       vrl4,       x1g,
     9   x2g,        x3g,        x4g,        y1g,
     a   y2g,        y3g,        y4g,        z1g,
     b   z2g,        z3g,        z4g,        ibid,
     c   iexpan,     gbuf%EINTTH,itask,      gbuf%VOL,
     d   actifxfem,  igre,       sensors,    nel,
     e   gbuf%G_WPLA,gbuf%WPLA )
       ENDIF
C--------------------------
C     TIME STEP
C--------------------------
      IF (ismstr /= 3.AND.(nodadt == 0.OR.idtmin(3) /= 0))THEN
        CALL cdt3(jft     ,jlt    ,ym        ,off      ,dt2t    ,
     2            neltst  ,ityptst,sti       ,stir     ,gbuf%OFF,
     3            dtcsub  ,ndt    ,dt2c      ,ixc      ,ssp     ,
     4            viscmx  ,px1    ,px2       ,py1      ,py2     ,
     5            vol0    ,vol00  ,rho       ,aldt     ,alpe    ,
     6            indxof  ,ngl    ,ismstr    ,iofc     ,msc     ,
     7            dmelc   ,jsms   ,gbuf%G_DT ,gbuf%DT)
      ENDIF
C----------------------------
C     FORCES INTERNES + ASSEMBLE
C----------------------------
      CALL cfint3(jft    ,jlt    ,gbuf%FOR,gbuf%MOM,thk0   ,thk02  ,
     2            px1    ,px2    ,py1     ,py2     ,area   ,z2     ,
     3            f11    ,f12    ,f13     ,f14     ,f21    ,f22    ,
     4            f23    ,f24    ,f31     ,f32     ,f33    ,f34    ,
     5            h11    ,h12    ,h13     ,h21     ,h22    ,h23    ,
     6            h31    ,h32    ,h33     ,b11     ,b12    ,b13    ,
     7            b14    ,b21    ,b22     ,b23     ,b24    ,nel    ,
     8            m11    ,m12    ,m13     ,m14     ,m21    ,m22    , 
     9            m23    ,m24    ,m31     ,m32     ,m33    ,m34    ,
     a            e1x    ,e1y    ,e1z     ,e2x     ,e2y    ,e2z    ,
     b            e3x    ,e3y    ,e3z     ,ihbe    ,npttot ,fzero  )
C-------------------------
c     Shell thermal
C--------------------------
C
       IF (jthe /= 0) THEN   ! /heat/mat
         IF (mat_elem%MAT_PARAM(mat(1))%HEAT_FLAG == 1) THEN
           CALL thermc(jft   ,jlt   ,pm     ,mat     ,thk0  ,ixc     ,
     .                 px1   ,px2   ,py1    ,py2     ,area  ,dt1c    ,
     .                 temp  ,tempel,fheat  , them   ,glob_therm%THEACCFACT) 
         ELSE
           CALL thermc(jft   ,jlt   ,pm     ,mat     ,thk0  ,ixc     ,
     .                 px1   ,px2   ,py1    ,py2     ,area  ,dt1c    ,
     .                 temp  ,tempel,die   , them    ,glob_therm%THEACCFACT) 
         END IF
       ENDIF
C--------------------------
C     THERMAL TIME STEP
C--------------------------
      imat = mat(1)
      IF (jthe /= 0 .AND. glob_therm%IDT_THERM == 1) THEN
           call dttherm(nel     ,pm(1,imat),npropm,glob_therm,mat_elem%mat_param(imat),
     .                  jtur    ,tempel  ,vol0    ,rho       ,
     .                  aldt    ,off     ,conde   ,gbuf%re   ,gbuf%rk )
      ENDIF 
c-------------------------
c     Virtual internal forces of regularized non local ddl 
c--------------------------
      IF (inloc > 0) THEN
        ! Warning : Y24 = PX1
        CALL cfint_reg(
     1   nloc_dmg,        var_reg,         gbuf%THK,        nel,
     2   off,             area,            nc1,             nc2,
     3   nc3,             nc4,             px1,             py1,
     4   px2,             py2,             elbuf_str%NLOC(1,1),
     5   ixc(1,jft),      nddl,            itask,           dt2t,
     6   aldt,            gbuf%THK_I,      gbuf%AREA,       nft)
      ENDIF
C--------------------------
      IF (iparit == 3) THEN
       CALL cupdt3f(jft     ,jlt  ,f      ,m    ,nvc      , 
     2              gbuf%OFF,off  ,sti    ,stir ,stifn    ,
     3              stifr   ,ixc  ,pm     ,area ,gbuf%THK ,
     4              f11     ,f12  ,f13    ,f14  ,f21      ,
     5              f22     ,f23  ,f24    ,f31  ,f32      ,
     6              f33     ,f34  ,m11    ,m12  ,m13      ,
     7              m14     ,m21  ,m22    ,m23  ,m24      ,
     8              m31     ,m32  ,m33    ,m34  ,gbuf%EINT,
     9              partsav ,mat  ,ipartc ,glob_therm%NODADT_THERM)
      ELSEIF (iparit == 0) THEN
       CALL cupdt3(jft     ,jlt   ,f     ,m    ,nvc      , 
     2             gbuf%OFF,off   ,sti   ,stir ,stifn    ,
     3             stifr   ,ixc   ,pm    ,area ,gbuf%THK ,
     4             f11     ,f12   ,f13   ,f14  ,f21      ,
     5             f22     ,f23   ,f24   ,f31  ,f32      ,
     6             f33     ,f34   ,m11   ,m12  ,m13      ,
     7             m14     ,m21   ,m22   ,m23  ,m24      ,
     8             m31     ,m32   ,m33   ,m34  ,gbuf%EINT,
     9             partsav ,mat   ,ipartc,jthe ,them     ,
     a             fthe    ,condn ,conde ,glob_therm%NODADT_THERM)
      ELSE
        CALL cupdt3p(jft      ,jlt     ,gbuf%OFF,off   ,sti    ,
     2               stir     ,fsky    ,fsky    ,iadc  ,ixc    ,
     3               f11      ,f12     ,f13     ,f14   ,f21    ,
     4               f22      ,f23     ,f24     ,f31   ,f32    ,
     5               f33      ,f34     ,m11     ,m12   ,m13    ,
     6               m14      ,m21     ,m22     ,m23   ,m24    ,
     7               m31      ,m32     ,m33     ,m34   ,
     8               gbuf%EINT,partsav ,mat     ,ipartc,pm     ,
     9               area     ,gbuf%THK,jthe    ,them  ,fthesky,
     a               condnsky ,conde   ,glob_therm%NODADT_THERM)
      ENDIF
c
      IF (icsen > 0)
     .  CALL csens3(jft    ,jlt    ,pid    ,igeo   ,epsd_pg)
C-------------------------
c     Shell cracking
C--------------------------
      IF (ixfem > 0) THEN
        DO ilay=1,nxlay                                                         
c
          ! crack length calculation for advancing crack                                            
          crklen => elbuf_str%BUFLY(ilay)%DMG(1:nel)
          CALL  crklen4n_adv(
     .           nel       ,nft       ,ilay      ,nlay      ,ixc       ,
     .           crklen    ,elcrkini  ,iel_crk   ,dir1_crk  ,dir2_crk  ,     
     .           nodedge   ,crkedge   ,xedge4n   ,ngl       ,x2        ,
     .           x3        ,x4        ,y2        ,y3        ,y4        ,
     .           aldt      )
c
          CALL crklayer4n_adv(
     .         xfem_str(ig,1:nxel) ,nel      ,nft      ,ixc      ,elcutc   ,    
     .         ilay     ,nxlay    ,iel_crk  ,inod_crk ,                
     .         iadc_crk ,nodenr   ,elcrkini ,dir1_crk ,dir2_crk ,                 
     .         nodedge  ,crknodiad,knod2elc ,crkedge  ,a_i      ,
     .         x2       ,x3       ,x4       ,y2       ,y3       ,
     .         y4       ,xedge4n  ,ngl      )                                                         
c
          CALL crklayer4n_ini(
     .         xfem_str(ig,1:nxel) ,nel      ,nft      ,ixc      ,elcutc   ,
     .         ilay     ,nxlay    ,iel_crk  ,inod_crk ,
     .         iadc_crk ,nodenr   ,elcrkini ,dir1_crk ,dir2_crk ,
     .         nodedge  ,crknodiad,knod2elc ,crkedge  ,a_i      ,
     .         x2       ,x3       ,x4       ,y2       ,y3       ,
     .         y4       ,xedge4n  ,ngl      )                                                         
        ENDDO                                                                   
C
        CALL crkoffc(elbuf_str,xfem_str(ig,1:nxel)  ,
     .               jft      ,jlt       ,nft    ,ir      ,is        ,
     .               nxlay    ,iel_crk   ,crkedge,xedge4n )
      END IF                                                             
c--------------------------------------------
c     Front wave
c--------------------------------------------
      IF (ifailwave > 0) THEN
        crkdir => elbuf_str%BUFLY(1)%CRKDIR
c        
        CALL set_failwave_nod4(failwave   ,fwave_el ,ngl      ,
     .       nel      ,ixc      ,itab     ,crkdir   ,dir_a    ,
     .       l_dira   ,x2       ,x3       ,x4       ,y2       ,
     .       y3       ,y4       )
      ENDIF
C-----------
      IF (ALLOCATED(dir2_crk)) DEALLOCATE(dir2_crk)                                                          
      IF (ALLOCATED(dir1_crk)) DEALLOCATE(dir1_crk)                                                          
      IF (ALLOCATED(elcrkini)) DEALLOCATE(elcrkini)                                                          
      IF (ALLOCATED(dirb))     DEALLOCATE(dirb)                                                          
      IF (ALLOCATED(dira))     DEALLOCATE(dira)           
      IF (ALLOCATED(var_reg))  DEALLOCATE(var_reg)  
C-----------
      RETURN