cforc3_crk.F File Reference

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

Go to the source code of this file.

Functions/Subroutines

subroutine

cforc3_crk (timers, xfem_str, jft, jlt, pm, ixc, x, f, m, v, vr, failwave, nvc, mtn, geo, tf, npf, bufmat, partsav, dt2t, neltst, ityptst, stifn, stifr, fsky, crksky, iadc, 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, kfts, ithk, iofc, ipla, nft, ismstr, fzero, igeo, ipm, ifailure, itask, jthe, temp, fthe, fthesky, iexpan, gresav, grth, igrth, msc, dmelc, jsms, table, iparg, ixfem, inod_crk, iel_crk, iadc_crk, elcutc, sensors, ixel, stack, isubstack, uxint_mean, uyint_mean, uzint_mean, nlevxf, nodedge, crkedge, drape_sh4n, ipri, nloc_dmg, indx_drape, igre, dt, snpc, stf, glob_therm, userl_avail, maxfunc, sbufmat, ipart)

Function/Subroutine Documentation

cforc3_crk()

subroutine cforc3_crk	(	type(timer_), intent(inout)	timers,
		type (elbuf_struct_), target	xfem_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,
		type (xfem_sky_), dimension(*)	crksky,
		integer, dimension(4,*)	iadc,
			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_)	group_param,
		type (mat_elem_), intent(inout)	mat_elem,
		integer	nel,
		integer	istrain,
		integer	ihbe,
		integer	kfts,
		integer	ithk,
		integer	iofc,
		integer	ipla,
		integer	nft,
		integer	ismstr,
			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(*)	igrth,
			msc,
			dmelc,
		integer	jsms,
		type(ttable), dimension(*)	table,
		integer, dimension(*)	iparg,
		integer	ixfem,
		integer, dimension(*)	inod_crk,
		integer, dimension(*)	iel_crk,
		integer, dimension(4,*)	iadc_crk,
		integer, dimension(2,*)	elcutc,
		type (sensors_), intent(inout)	sensors,
		integer	ixel,
		type (stack_ply)	stack,
		integer	isubstack,
			uxint_mean,
			uyint_mean,
			uzint_mean,
		integer	nlevxf,
		integer, dimension(2,*)	nodedge,
		type (xfem_edge_), dimension(*)	crkedge,
		type (drape_), dimension(numelc_drape)	drape_sh4n,
		integer	ipri,
		type (nlocal_str_)	nloc_dmg,
		integer, dimension(scdrape)	indx_drape,
		integer, intent(in)	igre,
		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 )

Definition at line 66 of file cforc3_crk.F.

C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE timer_mod
      USE table_mod
      USE crackxfem_mod
      USE mat_elem_mod
      USE stack_mod
      USE failwave_mod
      USE drape_mod  
      USE nlocal_reg_mod  
      USE sensor_mod
      USE elbufdef_mod
      USE dt_mod
      use glob_therm_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      "com_xfem1.inc"
#include      "scr02_c.inc"
#include      "scr17_c.inc"
#include      "scr18_c.inc"
#include      "param_c.inc"
#include      "parit_c.inc"
#include      "timeri_c.inc"
#include      "comlock.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 TF
      INTEGER,INTENT(IN)    :: SNPC        ! Size of NPC
      INTEGER, INTENT(IN) :: IGRE
      INTEGER JFT,JLT,NVC,MTN,NELTST,ITYPTST,OFFSET,
     .   NEL,ISTRAIN,IHBE,ITHK,IOFC,IPLA,NFT,
     .   ISMSTR,KFTS,IFAILURE,IEXPAN,IXEL,NLEVXF,IPRI
      INTEGER NPF(*),IXC(NIXC,*),IADC(4,*),INDXOF(MVSIZ),IPARTC(*),
     .   IGEO(NPROPGI,*),IPM(NPROPMI,*),ITASK,JTHE, JSMS,
     .   GRTH(*),IGRTH(*),IPARG(*),IXFEM,INOD_CRK(*),IEL_CRK(*),
     .   IADC_CRK(4,*),ELCUTC(2,*),ISUBSTACK,NODEDGE(2,*),INDX_DRAPE(SCDRAPE)
C     REAL
      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(*), x(3,*), f(*), m(*), v(*), vr(*),
     .   geo(npropg,*),tf(*), bufmat(*),partsav(*),stifn(*),stifr(*),
     .   fsky(8,*),d(*),dr(*),tani(6,*),eani(*), thke(*),dt2t,
     .   fzero(3,4,*),temp(*), fthe(*),fthesky(*),gresav(*),
     .   msc(*), dmelc(*),uxint_mean(nlevxf,mvsiz),
     .   uyint_mean(nlevxf,mvsiz),uzint_mean(nlevxf,mvsiz)
      TYPE(TTABLE) TABLE(*)
      TYPE (ELBUF_STRUCT_), TARGET :: XFEM_STR
      TYPE (XFEM_EDGE_)   , DIMENSION(*) :: CRKEDGE
      TYPE (XFEM_SKY_)    , DIMENSION(*) :: CRKSKY
      TYPE (STACK_PLY) :: STACK
      TYPE (FAILWAVE_STR_) :: FAILWAVE 
      TYPE (GROUP_PARAM_)  :: GROUP_PARAM
      TYPE (DRAPE_) , DIMENSION(NUMELC_DRAPE)  ::  DRAPE_SH4N
      TYPE (MAT_ELEM_)   ,INTENT(INOUT) :: MAT_ELEM
      TYPE (NLOCAL_STR_) :: NLOC_DMG 
      TYPE (SENSORS_) ,INTENT(INOUT) :: SENSORS
      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 I,IGTYP,ICSEN,IXLAY,NXLAY,NLAYER,IUN,NPG,NPTT,IBID,
     .   IR,IS,IT,NG,L_DIRA,L_DIRB,J1,J2,ILEV,IGMAT,IPTHK,IREP,
     .   ACTIFXFEM,SEDRAPE,NUMEL_DRAPE
      INTEGER MAT(MVSIZ),PID(MVSIZ),NDT(MVSIZ),NGL(MVSIZ),FWAVE(MVSIZ)
C
      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),
     .   dtcsub(mvsiz),a_i(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),
     .   die(mvsiz),tempel(mvsiz),them(mvsiz,4),
     .   ux1(mvsiz),ux2(mvsiz),ux3(mvsiz),ux4(mvsiz),
     .   uy1(mvsiz),uy2(mvsiz),uy3(mvsiz),uy4(mvsiz),a11r(mvsiz),
     .   thke0(mvsiz)
      my_real 
     .   bid,thkr
C---
      INTEGER, ALLOCATABLE, DIMENSION(:) :: ELCRKINI
      my_real, 
     .    ALLOCATABLE, DIMENSION(:) :: dira,dirb,dir1_crk,dir2_crk
      my_real,
     .  DIMENSION(:) ,POINTER  ::  dir_a,dir_b
      TARGET :: dira,dirb
!     variables for heat transfer
      my_real, dimension(mvsiz) :: fheat
      my_real, dimension(mvsiz) :: ssp_eq
!
C---
      TYPE(BUF_LAY_) ,POINTER :: BUFLY
      TYPE(G_BUFEL_) ,POINTER :: GBUF
      TYPE(L_BUFEL_) ,POINTER :: LBUF
      my_real,
     .  DIMENSION(:) ,POINTER  ::  offg,thkg,strag,forg,momg,
     .  eintg,epsdg,tempg,eintth,hourgg
      DOUBLE PRECISION,
     .  DIMENSION(:) ,POINTER  ::  SMSTRG
      my_real,
     .  DIMENSION(:,:), ALLOCATABLE :: varnl
      INTEGER NPTTOT
      INTEGER SDIR_A    ! Size of DIR_A
      INTEGER SDIR_B    ! Size of DIR_B
C=======================================================================
      ibid = 0
      iun  = 1
      bid  = zero
      npg  = 0
      ng   = 1
      ir   = 1
      is   = 1
      it   = 1
      sedrape = scdrape
      numel_drape = numelc_drape
      npttot = 1
      ALLOCATE(varnl(nel,1))
      varnl = zero
C -----will put it in starter     
      IF (ismstr>=10) ismstr=4
C
      DO i=jft,jlt
        mat(i) = ixc(1,i)
        pid(i) = ixc(6,i)
        ngl(i) = ixc(7,i)
      ENDDO
C
      icsen = igeo(3,pid(1))
      igtyp = igeo(11,pid(1))
      igmat = igeo(98,pid(1))
      irep  = iparg(35)
      actifxfem = iparg(70)
C-----------------------------------------
      gbuf => xfem_str%GBUF
      nxlay = xfem_str%NLAY
      nlayer = 1    ! only one current xfem layer passes to the cmain3 & mulawc
C-----------------------------------------
      ALLOCATE(elcrkini(nxlay*nel))
      ALLOCATE(dir1_crk(nxlay*nel))
      ALLOCATE(dir2_crk(nxlay*nel))
      elcrkini = 0
      dir1_crk = zero
      dir2_crk = zero
C
      l_dira = xfem_str%BUFLY(1)%LY_DIRA
      l_dirb = xfem_str%BUFLY(1)%LY_DIRB
C
      ALLOCATE (dira(nxlay*nel*l_dira))
      ALLOCATE (dirb(nxlay*nel*l_dirb))
      dira=zero
      dirb=zero
      sdir_a=nxlay*nel*l_dira
      sdir_b=nxlay*nel*l_dirb
      dir_a => dira(1:nxlay*nel*l_dira)
      dir_b => dirb(1:nxlay*nel*l_dirb)
C
      DO ixlay=1,nxlay
        nptt = xfem_str%BUFLY(ixlay)%NPTT
        IF (l_dira == 0) THEN
          dira = zero
        ELSEIF (irep == 0) THEN
          j1 = 1+(ixlay-1)*l_dira*nel
          j2 = ixlay*l_dira*nel
          dira(j1:j2) = xfem_str%BUFLY(ixlay)%DIRA(1:nel*l_dira)
        ENDIF
C
        DO i=jft,jlt
          sigy(i)    = ep30
          zcfac(i,1) = one
          zcfac(i,2) = one
          alpe(i)    = one
          tempel(i)  = zero
        ENDDO
C---
        ilev = nxel*(ixlay-1) + ixel
C---
        IF (igtyp == 1 .or. igtyp == 9) THEN
          DO i=jft,jlt
            thke0(i) = thke(i)
          ENDDO
        ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
          ipthk = 1 + nxlay
          thkr  = stack%GEO(ipthk+ixlay,isubstack)
          DO i=jft,jlt
            thke0(i) = thke(i) * thkr     ! layer thickness (real)
          ENDDO
        ELSE    !  IGTYP == 11...
          ipthk = 300
          DO i=jft,jlt
            thke0(i) = thke(i) * geo(ipthk+ixlay,pid(i))
          ENDDO
        ENDIF
c
        IF (nxlay > 1) THEN
          lbuf => xfem_str%BUFLY(ixlay)%LBUF(ir,is,1)
          bufly  => xfem_str%BUFLY(ixlay)
          offg   => lbuf%OFF
          smstrg => lbuf%SMSTR
          thkg   => lbuf%THK
          strag  => lbuf%STRA
          forg   => lbuf%FOR
          momg   => lbuf%MOM
          eintg  => lbuf%EINT
          epsdg  => lbuf%EPSD
          tempg  => lbuf%TEMP
          eintth => lbuf%EINTTH
          hourgg => bufly%HOURG
        ELSEIF (nxlay == 1) THEN
          offg   => gbuf%OFF
          smstrg => gbuf%SMSTR
          thkg   => gbuf%THK
          strag  => gbuf%STRA
          forg   => gbuf%FOR
          momg   => gbuf%MOM
          eintg  => gbuf%EINT
          epsdg  => gbuf%EPSD
          tempg  => gbuf%TEMP
          eintth => gbuf%EINTTH
          hourgg => gbuf%HOURG
        ENDIF
C
        CALL ccoor3_crk(jft     ,jlt  ,nft   ,ixc    ,ngl     ,
     2                  pid     ,mat  ,ilev  ,iel_crk,iadc_crk,
     3                  offg    ,off  ,sigy  ,thke0  ,thk0    ,
     4                  dt1c    ,vl1  ,vl2   ,vl3    ,vl4     ,
     5                  vrl1    ,vrl2 ,vrl3  ,vrl4   ,x1g     ,
     6                  x2g     ,x3g  ,x4g   ,y1g    ,y2g     ,
     7                  y3g     ,y4g  ,z1g   ,z2g    ,z3g     ,
     8                  z4g     )
c
        CALL cnvec3(xfem_str,dir_a   ,dir_b   ,                     
     1     jft     ,jlt      ,irep    ,igtyp   ,nxlay   ,           
     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     )                             
C
        IF (ismstr /= 3) THEN
          CALL cderi3(
     1   jft,     jlt,     smstrg,  offg,
     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                smstrg    ,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      ,thkg    ,
     5              ssp     ,rho    ,h1     ,h2       ,h3      ,
     6              vol0    ,vol00  ,alpe   ,gs       ,mtn     ,
     7              ithk    ,ismstr  ,nptt   ,kfts   ,
     8              srh1   ,srh2   ,srh3     ,igeo    ,
     9              a11r    ,isubstack      ,stack%PM)
C
        IF ((nodadt == 0 .AND. ismstr /= 3 .AND. idt1sh == 0) .OR.
     .       idtmin(3) /= 0 .OR. igtyp == 16)
c         CDLEN3 must be always called for tissue property/law for time step correction
     .    CALL cdlen3(jft     ,jlt     ,pm      ,off     ,area,
     .                x2      ,x3      ,x4      ,y2      ,y3  ,
     .                y4      ,aldt    ,mat     ,geo     ,pid ,
     .                ihbe    )
        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              bid ,ismstr   )
        CALL cstra3(jft    ,jlt      ,strag ,shf    ,area, 
     2              exx    ,eyy      ,exy      ,exz    ,eyz ,   
     3              kxx    ,kyy      ,kxy      ,dt1c   ,tani, 
     4              forg   ,momg     ,ismstr   ,mtn ,
     6              ihbe   ,nft      ,istrain  ,ux1    ,ux2 ,
     7              ux3    ,ux4      ,uy1      ,uy2    ,uy3 ,
     8              uy4    ,px1      ,px2      ,py1    ,py2 ,
     9              bid    ,bid      ,bid      ,nel    )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
        DO i = jlt,jft
          tempel(i) = zero
        ENDDO
C
        IF (jthe > 0) CALL tempcg(numnod, nel  ,ixc ,temp ,tempel)
C-----------------------------
        IF ((itask==0).AND.(imon_mat==1)) CALL startime(timers,35)
C-----------------------------
        CALL cmain3(timers,
     1    xfem_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    forg      ,momg      ,strag     ,failwave  ,fwave     ,
     9    thkg      ,eintg     ,iofc      ,
     a    g         ,a11       ,a12       ,vol0      ,indxof    ,
     b    ngl       ,zcfac     ,shf       ,gs        ,epsdg     ,
     c    kfts      ,ihbe      ,alpe      ,
     d    dir_a     ,dir_b     ,igeo      ,
     e    ipm       ,ifailure  ,npg       ,fheat     ,
     f    tempel    ,die       ,jthe      ,iexpan    ,tempg     ,
     g    ibid      ,bid       ,
     h    bid       ,bid       ,bid       ,bid       ,bid       ,
     i    bid       ,bid       ,bid       ,e1x       ,e1y       ,
     j    e1z       ,e2x       ,e2y       ,e2z       ,e3x       ,
     k    e3y       ,e3z       ,ibid      ,table     ,ixfem     ,
     l    bid       ,sensors   ,bid       ,elcrkini  ,
     m    dir1_crk  ,dir2_crk  ,aldt      ,glob_therm%IDT_THERM ,glob_therm%THEACCFACT,
     n    ismstr    ,ir        ,is        ,nlayer    ,nptt      ,
     o    ixlay     ,ixel      ,isubstack ,stack     ,
     p    bid       ,itask     ,drape_sh4n  ,varnl     ,nloc_dmg ,
     r    indx_drape ,thke, sedrape  ,numel_drape,dt   ,
     q    ncycle     ,snpc      ,  stf,nxlaymax ,
     s    idel7nok   ,userl_avail ,maxfunc  ,npttot    ,sbufmat,
     t    sdir_a     ,sdir_b   ,gbuf%FOR_G,ssp_eq,
     x    ipart      ,lipart1  ,ipartc   )
C-----------------------------
        IF ((itask==0).AND.(imon_mat==1)) 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 (nptt == 1) THEN
          CALL mhvis3(jft    ,jlt    ,pm     ,thkg    ,hourgg    ,
     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(
     1            jft    ,jlt    ,thkg   ,hourgg ,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 ,thkg,hourgg,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,      thkg,     eintg,    partsav,
     3   area,     mat,      ipartc,   x,
     4   vr,       vol0,     vol00,    thk0,
     5   thk02,    iun,      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,      ixfem,
     c   iexpan,   eintth,   itask,    gbuf%VOL,
     d   actifxfem,igre,sensors,nel,gbuf%G_WPLA,
     e   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   ,offg    ,        
     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    ,forg    ,momg    ,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    ,nptt    ,fzero  )
C-------------------------
c     Shell thermal
C--------------------------
        IF (jthe /= 0) THEN
          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       --- to be added --- (see cforc3.F)
C--------------------------
C
C a faire:
cc      IF(ISUB /= 0)THEN
cc        CALL CSUBF3(JFT  ,JLT  ,NDT ,AREAS,AREA,
cc     4            F11    ,F12  ,F13 ,F14  ,F21 ,
cc     5            F22    ,F23  ,F24 ,F31  ,F32 ,
cc     6            F33    ,F34  ,M11 ,M12  ,M13 ,
cc     7            M14    ,M21  ,M22 ,M23  ,M24 ,
cc     8            M31    ,M32  ,M33 ,M34  )
cc      ENDIF
C
        IF (icsen > 0)
     .    CALL csens3(jft    ,jlt    ,pid    ,igeo   ,epsdg)
C--------------------------
C     ASSEMBLE
C--------------------------
        IF (iparit == 1)
     .    CALL cupdt3_crk(
     .         jft  ,jlt  ,nft  ,ixc   ,off     ,iadc   ,
     .         f11  ,f21  ,f31  ,f12   ,f22     ,f32    ,
     .         f13  ,f23  ,f33  ,f14   ,f24     ,f34    ,
     .         m11  ,m21  ,m31  ,m12   ,m22     ,m32    ,
     .         m13  ,m23  ,m33  ,m14   ,m24     ,m34    ,
     .         sti  ,stir ,fsky ,elcutc,iadc_crk,iel_crk,
     .         ilev ,inod_crk,offg,eintg,partsav,ipartc,
     .         ixlay,crksky)
C-------------------------
      ENDDO  !  DO IXLAY=1,NXLAY
C-------------------------
      IF (ALLOCATED(dira))     DEALLOCATE(dira)
      IF (ALLOCATED(dirb))     DEALLOCATE(dirb)
      IF (ALLOCATED(elcrkini)) DEALLOCATE(elcrkini)
      IF (ALLOCATED(dir1_crk)) DEALLOCATE(dir1_crk)
      IF (ALLOCATED(dir2_crk)) DEALLOCATE(dir2_crk)
      IF (ALLOCATED(varnl))    DEALLOCATE(varnl)
C------------------------- 
      RETURN