czforc3_crk.F File Reference

#include "implicit_f.inc"
#include "param_c.inc"
#include "mvsiz_p.inc"
#include "com_xfem1.inc"
#include "scr14_c.inc"
#include "parit_c.inc"
#include "timeri_c.inc"
#include "comlock.inc"
#include "com04_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine

czforc3_crk (timers, xfem_str, jft, jlt, nft, ityptst, ipari, mtn, ipri, ithk, neltst, istrain, ipla, tt, dt1, dt2t, pm, geo, partsav, ixc, group_param, bufmat, tf, npf, iadc, failwave, x, d, dr, v, vr, f, m, stifn, stifr, fsky, tani, offset, eani, indxof, ipartc, thke, nvc, iofc, ihbe, f11, f12, f13, f14, f21, f22, f23, f24, f31, f32, f33, f34, m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, kfts, fzero, ismstr, mat_elem, 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, crksky, sensors, ixel, isubstack, uxint_mean, uyint_mean, uzint_mean, nlevxf, nodedge, crkedge, stack, drape_sh4n, nloc_dmg, indx_drape, igre, dt, ncycle, snpc, stf, glob_therm, idel7nok, userl_avail, maxfunc, sbufmat, ipart, lipart1)

Function/Subroutine Documentation

czforc3_crk()

subroutine czforc3_crk	(	type(timer_), intent(inout)	timers,
		type (elbuf_struct_), target	xfem_str,
		integer	jft,
		integer	jlt,
		integer	nft,
		integer	ityptst,
		integer	ipari,
		integer	mtn,
		integer	ipri,
		integer	ithk,
		integer	neltst,
		integer	istrain,
		integer	ipla,
			tt,
			dt1,
			dt2t,
			pm,
			geo,
			partsav,
		integer, dimension(nixc,*)	ixc,
		type (group_param_)	group_param,
		target	bufmat,
			tf,
		integer, dimension(*)	npf,
		integer, dimension(4,*)	iadc,
		type (failwave_str_)	failwave,
			x,
			d,
			dr,
			v,
			vr,
			f,
			m,
			stifn,
			stifr,
			fsky,
			tani,
		integer	offset,
			eani,
		integer, dimension(mvsiz)	indxof,
		integer, dimension(*)	ipartc,
			thke,
		integer	nvc,
		integer	iofc,
		integer	ihbe,
			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	kfts,
			fzero,
		integer	ismstr,
		type (mat_elem_), intent(inout)	mat_elem,
		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 (xfem_sky_), dimension(*)	crksky,
		type (sensors_), intent(inout)	sensors,
		integer	ixel,
		integer	isubstack,
			uxint_mean,
			uyint_mean,
			uzint_mean,
		integer	nlevxf,
		integer, dimension(2,*)	nodedge,
		type (xfem_edge_), dimension(*)	crkedge,
		type (stack_ply)	stack,
		type (drape_), dimension(numelc_drape)	drape_sh4n,
		type (nlocal_str_)	nloc_dmg,
		integer, dimension(scdrape)	indx_drape,
		integer, intent(in)	igre,
		type (dt_), intent(in)	dt,
		integer, intent(in)	ncycle,
		integer, intent(in)	snpc,
		integer, intent(in)	stf,
		type (glob_therm_), intent(inout)	glob_therm,
		integer, intent(inout)	idel7nok,
		integer, intent(in)	userl_avail,
		integer, intent(in)	maxfunc,
		integer, intent(in)	sbufmat,
		integer, dimension(lipart1, npart), intent(in)	ipart,
		integer, intent(in)	lipart1 )

Definition at line 65 of file czforc3_crk.F.

C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE timer_mod
      USE table_mod
      USE crackxfem_mod
      USE stack_mod
      USE failwave_mod
      USE mat_elem_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      "param_c.inc"
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C O M M O N   B L O C K S
C-----------------------------------------------
#include      "com_xfem1.inc"
#include      "scr14_c.inc"
#include      "parit_c.inc"
#include      "timeri_c.inc"
#include      "comlock.inc"
#include      "com04_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(INOUT) :: IDEL7NOK    ! Element deletion flag for IDEL flag in contact interfaces
      INTEGER,INTENT(IN)    :: SBUFMAT     ! Size of BUFMAT
      INTEGER,INTENT(IN)    :: STF         ! Size of TF
      INTEGER,INTENT(IN)    :: SNPC        ! Size of NPF
      INTEGER JFT ,JLT,NFT,IPARI,MTN,IPRI,ITHK,NELTST,
     .        ITYPTST,ISTRAIN,IPLA,OFFSET,NVC,JSMS,IOFC,IHBE,
     .        KFTS,ISMSTR,IFAILURE,JTHE,IXEL,ISUBSTACK,NLEVXF
      INTEGER IXC(NIXC,*),IADC(4,*),IPARTC(*),NPF(*),IGEO(NPROPGI,*),
     .        IPM(NPROPMI,*),INDXOF(MVSIZ),ITASK,IEXPAN,GRTH(*),IGRTH(*),
     .        IPARG(*),IXFEM,INOD_CRK(*),IEL_CRK(*),IADC_CRK(4,*),
     .        ELCUTC(2,*),NODEDGE(2,*),INDX_DRAPE(SCDRAPE)
      INTEGER, INTENT(IN) :: IGRE, NCYCLE
      integer, intent(in) :: LIPART1      
      INTEGER, DIMENSION(LIPART1, NPART), INTENT(IN)     :: IPART
     
C     REAL OR REAL*8
      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),
     .   tf(*),pm(npropm,*),geo(npropg,*),partsav(*),
     .   bufmat(*),x(3,*),d(*),dr(*),v(3,*),vr(3,*),
     .   f(3,*),m(3,*),stifn(*), stifr(*),fsky(8,*),tani(6,*),
     .   eani(*),thke(*),fzero(3,4,*),temp(*),fthe(*),
     .   fthesky(*),gresav(*),msc(*), dmelc(*),
     .   uxint_mean(nlevxf,mvsiz),
     .   uyint_mean(nlevxf,mvsiz),uzint_mean(nlevxf,mvsiz)
      my_real    
     .   tt,dt1,dt2t
      TARGET :: bufmat
      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
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
C   L O C A L   V A R I A B L E S
C--------------------------------
      LOGICAL PLAT(MVSIZ)
      INTEGER 
     .   I,NEL,NNOD,NPG,IBID,ISROT,IXLAY,NXLAY,NLAYER,NPTT,
     .   IR,IS,IT,NG,L_DIRA,L_DIRB,J1,J2,IXFEM2,IPTHK,
     .   IORTH,ILEV,IGMAT,IGTYP,IREP,IMAT,
     .   ACTIFXFEM, SEDRAPE,NUMEL_DRAPE
      INTEGER  MAT(MVSIZ),PID(MVSIZ),NGL(MVSIZ),FWAVE(MVSIZ)
      parameter(nnod = 4)
      my_real 
     .  ll(mvsiz),vl13(mvsiz,3),vl24(mvsiz,3),vlhi(mvsiz,3),
     .  rlxyz(mvsiz,8),x13(mvsiz),x24(mvsiz),y13(mvsiz),y24(mvsiz),
     .  mx13(mvsiz),mx23(mvsiz),mx34(mvsiz),my13(mvsiz),my23(mvsiz),
     .  my34(mvsiz),z1(mvsiz),z2(mvsiz)
      my_real 
     .   vq(mvsiz,9),vdef(mvsiz,8),off(mvsiz),dhg(mvsiz,6),
     .   area(mvsiz),vqn(mvsiz,12),l13(mvsiz),l24(mvsiz),
     .   zcfac(mvsiz,2),a_i(mvsiz),thk02(mvsiz)
      my_real 
     .   exx(mvsiz),eyy(mvsiz),exy(mvsiz),exz(mvsiz),eyz(mvsiz),
     .   kxx(mvsiz),kyy(mvsiz),kxy(mvsiz),sigy(mvsiz), 
     .   dt1c(mvsiz),ssp(mvsiz),viscmx(mvsiz),rho(mvsiz) ,
     .   nu(mvsiz),g(mvsiz),a11(mvsiz),a12(mvsiz),vol0(mvsiz),
     .   thk0(mvsiz),sti(mvsiz),stir(mvsiz),shf(mvsiz) ,
     .   vf(mvsiz,12),vm(mvsiz,8),gs(mvsiz),fac1(mvsiz),
     .   alpe(mvsiz),ym(mvsiz),facn(mvsiz,2),lxyz(mvsiz,8),
     .   dd(mvsiz,6),db(mvsiz,12),amu(mvsiz),gsr(mvsiz),
     .   a11sr(mvsiz),a12sr(mvsiz),nusr(mvsiz),shfsr(mvsiz),
     .   px2(mvsiz),py1(mvsiz),py2(mvsiz),die(mvsiz),tempel(mvsiz),
     .   them(mvsiz,4),r11(mvsiz),r12(mvsiz),r13(mvsiz),
     .   r21(mvsiz),r22(mvsiz),r23(mvsiz),r31(mvsiz),
     .   r32(mvsiz),r33(mvsiz),rlz(mvsiz,4),vhgzk(mvsiz,5),
     .   vhgze(mvsiz,5),vrlz(mvsiz),bm0rz(mvsiz,4,4),bmkrz(mvsiz,4,4),
     .   bmerz(mvsiz,4,4),vmz(mvsiz,4),krz(mvsiz),diz(mvsiz,3)
      my_real 
     .   x1g(mvsiz),x2g(mvsiz),x3g(mvsiz),x4g(mvsiz),
     .   y1g(mvsiz),y2g(mvsiz),y3g(mvsiz),y4g(mvsiz),
     .   z1g(mvsiz),z2g(mvsiz),z3g(mvsiz),z4g(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),
     .   vx1(mvsiz),vx2(mvsiz),vx3(mvsiz),vx4(mvsiz),
     .   vy1(mvsiz),vy2(mvsiz),vy3(mvsiz),vy4(mvsiz),
     .   vz1(mvsiz),vz2(mvsiz),vz3(mvsiz),vz4(mvsiz),
     .   vrx1(mvsiz),vrx2(mvsiz),vrx3(mvsiz),vrx4(mvsiz),
     .   vry1(mvsiz),vry2(mvsiz),vry3(mvsiz),vry4(mvsiz),
     .   vrz1(mvsiz),vrz2(mvsiz),vrz3(mvsiz),vrz4(mvsiz),
     .   ux1(mvsiz),ux2(mvsiz),ux3(mvsiz),ux4(mvsiz),
     .   uy1(mvsiz),uy2(mvsiz),uy3(mvsiz),uy4(mvsiz),a11r(mvsiz),
     .   thke0(mvsiz),xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),
     .   yl2(mvsiz),yl3(mvsiz),yl4(mvsiz),fac58(mvsiz,2)
      my_real  bid,thkr
      !
      my_real, DIMENSION(MVSIZ) :: zoffset
C---
      INTEGER, ALLOCATABLE, DIMENSION(:) :: ELCRKINI
      my_real, 
     .    ALLOCATABLE, DIMENSION(:) :: dira,dirb,dir1_crk,dir2_crk
      my_real,
     .  DIMENSION(:) ,POINTER  ::  dir_a,dir_b,uvar
      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 
     .   hm(mvsiz,6), hf(mvsiz,6), hc(mvsiz,2),hmfor(mvsiz,6)
      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=======================================================================
      nel  = jlt-jft+1
      bid  = zero
      ibid = 0
      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 (ismst r>= 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
      DO i=jft,jlt                
        fac1(i) = geo(17,pid(i))  
      ENDDO
C      
      isrot = iparg(41)
      igtyp = igeo(11,pid(1))
      igmat = igeo(98,pid(1))
      irep  = iparg(35)
      actifxfem = iparg(70)
      imat  = ixc(1,jft)
C-----------------------------------------
      gbuf => xfem_str%GBUF
C
      nxlay  = xfem_str%NLAY
      nlayer = 1    ! only one current xfem layer passes to the cmain3 & mulawc
C-----------------------------------------
      ixfem2 = 1               ! flag for CZCORC1 and BALANCE
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)
      uvar  => xfem_str%BUFLY(1)%MAT(ir,is,it)%VAR
c-------------------------------------------------
c     Loop over all (Xfem) layers
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
        ENDDO
C-----------------------------------------
C       CALCULS PRELIMINAIRES GLOBAL-->LOCAL
C-----------------------------------------
        ilev = nxel*(ixlay-1) + ixel
C---
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
          thkr  = geo(ipthk+ixlay,pid(1))
          DO i=jft,jlt
            thke0(i) = thke(i) * thkr
          ENDDO
        ENDIF
c
        IF (nxlay > 1) THEN
          lbuf   => xfem_str%BUFLY(ixlay)%LBUF(ir,is,it)
          bufly  => xfem_str%BUFLY(ixlay)
          offg   => lbuf%OFF
          smstrg => lbuf%SMSTR
          thkg   => lbuf%THK     ! layer thickness, not integration pt ! 
          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
C       !Coords + Calculation + lives of ghost elements instead of STD
        CALL ccoor3z_crk(jft  ,jlt  ,nft  ,iel_crk,iadc_crk,
     2                   vl1  ,vl2  ,vl3  ,vl4    ,vrl1    ,
     3                   vrl2 ,vrl3 ,vrl4 ,x1g    ,x2g     ,
     4                   x3g  ,x4g  ,y1g  ,y2g    ,y3g     ,
     5                   y4g  ,z1g  ,z2g  ,z3g    ,z4g     ,
     6                   vx1  ,vx2  ,vx3  ,vx4    ,vy1     ,
     7                   vy2  ,vy3  ,vy4  ,vz1    ,vz2     ,
     8                   vz3  ,vz4  ,vrx1 ,vrx2   ,vrx3    ,
     9                   vrx4 ,vry1 ,vry2 ,vry3   ,vry4    ,
     a                   vrz1 ,vrz2 ,vrz3 ,vrz4   ,ilev    ,
     b                   offg )
C
        CALL czcorc1(numnod ,numelc ,xfem_str,
     1               jft    ,jlt    ,x       ,v       ,vr     ,
     2               ixc    ,pm     ,plat    ,area    ,
     3               a_i    ,vl13   ,vl24    ,vlhi    ,rlxyz  ,
     4               vqn    ,vq     ,ll      ,l13     ,l24    ,
     5               x13    ,x24    ,y13     ,y24     ,mx13   ,
     6               mx23   ,mx34   ,my13    ,my23    ,my34   ,
     7               z1     ,lxyz   ,dd      ,db      ,smstrg  ,
     9               irep   ,nptt   ,nxlay   ,ismstr  ,
     a               dir_a  ,dir_b  ,offg    ,rlxyz   ,lxyz   ,
     b               facn   ,py1    ,px2     ,py2     ,r11    ,
     c               r12    ,r13    ,r21     ,r22     ,r23    ,
     d               r31    ,r32    ,r33     ,rlz     ,isrot  ,
     e               ixfem2 ,vx1    ,vx2     ,vx3     ,vx4    ,
     f               vy1    ,vy2    ,vy3     ,vy4     ,vz1    ,
     g               vz2    ,vz3    ,vz4     ,vrx1    ,vrx2   ,
     h               vrx3   ,vrx4   ,vry1    ,vry2    ,vry3   ,
     i               vry4   ,vrz1   ,vrz2    ,vrz3    ,vrz4   ,
     j               x1g    ,x2g    ,x3g     ,x4g     ,y1g    ,
     k               y2g    ,y3g    ,y4g     ,z1g     ,z2g    ,
     l               z3g    ,z4g    ,thke0   ,diz     ,ux1    ,
     m               ux2    ,ux3    ,ux4     ,uy1     ,uy2    ,
     n               uy3    ,uy4    ,xl2     ,xl3     ,xl4    ,
     o               yl2    ,yl3    ,yl4     ,vl1     ,vl2    ,
     p               vl3    ,vl4    ,nel     ,z2      )
C
        CALL cncoef3b(jft    ,jlt     ,pm      ,mat    ,geo    ,
     2                pid    ,area    ,shf     ,thk0   ,
     3                thk02  ,nu      ,g       ,ym     ,
     4                a11    ,a12     ,thkg    ,thke0  ,ssp    ,
     5                rho    ,vol0    ,gs      ,mtn    ,ithk   ,
     6                nptt   ,dt1c    ,dt1     ,ihbe   ,amu    ,
     7                gsr    ,a11sr   ,a12sr   ,nusr   ,shfsr  ,
     8                krz    ,igeo    ,a11r  ,isubstack, stack%PM, 
     9                mat_elem%MAT_PARAM(imat)%UPARAM ,dira    ,dirb    ,
     a                uvar   ,fac58  ,nel     ,zoffset)
C
        CALL cncoefort(jft    ,jlt    ,pm       ,mat     ,geo    ,  
     1                 pid    ,mtn    ,nptt     ,hm      ,hf     ,  
     2                 hc     ,hmfor  ,iorth    ,dir_a   ,igeo   ,  
     3                 isubstack,stack,xfem_str ,nxlay   ,thkg   ,
     4                 drape_sh4n ,nft  ,nel    ,indx_drape , thke,
     5                 sedrape,numel_drape , mat_elem)           
C----------------------------------
C     STRAIN RATE CALCULATION
C----------------------------------
        CALL czdef(jft     ,jlt  ,area ,a_i  ,vl13   ,vl24  ,
     2             vlhi    ,rlxyz,vdef ,dhg  ,x13    ,
     3             x24     ,y13  ,y24  ,mx13 ,mx23   ,mx34  ,
     4             my13    ,my23 ,my34 ,z1   ,dt1    ,off   ,
     5             offg    ,rlxyz)
        IF (isrot > 0) THEN
          CALL czdefrz(jft  ,jlt  ,area ,a_i ,rlz   ,
     1                 vdef ,vhgzk,vhgze,x13 ,x24   ,
     2                 y13  ,y24  ,mx13 ,mx23,mx34  ,
     3                 my13 ,my23 ,my34 ,z1  ,dhg   ,
     4                 bm0rz,bmkrz,bmerz,vl13,vl24  ,
     5                 vrlz )
        ENDIF
C----------------------------------
C     STRAIN CALCULATION
C----------------------------------
        CALL czstra3(jft    ,jlt    ,nft, vdef,strag    ,
     2               exx    ,eyy    ,exy, exz ,eyz      ,  
     3               kxx    ,kyy    ,kxy, dt1c,tani     ,
     4               iepsdot,istrain,ux1 ,ux2 ,ux3      ,
     5               ux4    ,uy1    ,uy2 ,uy3 ,uy4      ,
     6               y24    ,px2    ,py1 ,py2 ,area     ,
     7               ismstr ,mtn    ,bid ,bid ,bid      ,
     8               nel    )
C-----------------
C     CONTRAINTES
C-----------------
        DO i = jft,jlt
          tempel(i) = zero
        ENDDO
        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       ,r11       ,r12     ,
     j       r13       ,r21       ,r22       ,r23       ,r31     ,
     k       r32       ,r33       ,ibid      ,table     ,ixfem   ,
     l       bid       ,sensors   ,bid       ,elcrkini,
     m       dir1_crk  ,dir2_crk  ,ll        ,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  ,idel7nok  ,
     s       userl_avail ,maxfunc ,npttot      ,sbufmat   ,sdir_a    ,
     t       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     BALANCES BY MATERIAL
C--------------------------
        IF (ipri == 1)
     1    CALL cbilan(
     1   jft,      jlt,      pm,       v,
     2   ixc,      thkg,     eintg,    partsav,
     3   area,     mat,      ipartc,   bid,
     4   bid,      bid,      bid,      bid,
     5   bid,      ibid,     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,      ixfem2,
     c   iexpan,   eintth,   itask,    gbuf%VOL,
     d   actifxfem,igre,sensors,nel,gbuf%G_WPLA,
     e   gbuf%WPLA )
C
        CALL cndt3(
     1             jft    ,jlt    ,off    , dt2t   ,amu   ,
     2             neltst ,ityptst,sti    , stir   ,offg  ,
     3             ssp    ,viscmx ,rho    , vol0   ,thk0  ,thk02,
     4             a11    ,ll     ,alpe   , ngl    ,ismstr,
     5             iofc   ,nnod   ,area   , g      ,shf   ,
     6             msc    ,dmelc  ,jsms   , bid    ,igtyp ,
     7             igmat  ,a11r   ,gbuf%G_DT, gbuf%DT,mtn   ,
     8             pm     ,mat(jft),nel     ,zoffset,ssp_eq)
        CALL czfintce(jft  ,jlt  ,thk0 ,thk02,a_i    ,x13    ,
     2                x24  ,y13  ,y24  ,z1   ,mx23   ,mx13   ,
     3                mx34 ,my13 ,my23 ,my34 ,forg   ,momg   ,
     4                vf   ,vm   ,nel  )
        IF (isrot > 0) THEN
          CALL czfintcrz(jft      ,jlt  ,thk0 ,vol0 ,area  ,x13   ,
     2                   x24      ,y13  ,y24  ,z1   ,mx23  ,mx13  ,
     3                   mx34     ,my13 ,my23 ,my34 ,forg  ,hourgg,
     4                   vf       ,vmz  ,bm0rz,krz  ,vrlz  ,dt1c  ,
     5                   eintg    ,off  ,nel  )
          CALL czfintnrz_or(
     1         jft     ,jlt     ,thk0   ,thk02    ,a_i  ,dhg       ,
     2         x13     ,x24     ,y13    ,y24      ,z1   ,mx23      ,
     3         mx13    ,mx34    ,my13   ,my23     ,my34 ,hourgg    ,
     4         forg    ,momg    ,vf     ,vm       ,zcfac,a11       ,
     5         a12     ,g       ,gs     ,sigy     ,off  ,fac1      , 
     6         rho     ,area    ,dt1    ,eintg    ,amu  ,vlhi      ,
     7         nptt    ,ipartc  ,partsav,kfts     ,gsr       ,
     8         a11sr   ,a12sr   ,nusr   ,shfsr    ,bmkrz,bmerz     ,
     9         vhgzk   ,vhgze   ,krz    ,vmz      ,iorth,hm        ,
     a         hf      ,hc      ,hmfor  ,mtn      ,nel  )
        ELSE
         IF (iorth == 0) THEN
          CALL czfintn1(jft    ,jlt    ,thk0   ,thk02,a_i  ,dhg   ,
     2                  x13    ,x24    ,y13    ,y24  ,z1   ,mx23  ,
     3                  mx13   ,mx34   ,my13   ,my23 ,my34 ,hourgg,
     4                  forg   ,momg   ,vf     ,vm   ,zcfac,a11   ,
     5                  a12    ,g      ,shf    ,sigy ,off  ,fac1  ,
     6                  rho    ,area   ,dt1    ,eintg,amu  ,vlhi  ,
     7                  nptt   ,ipartc ,partsav,kfts ,gsr  ,nel   ,
     8                  a11sr  ,a12sr  ,nusr   ,shfsr,mtn  ,fac58 )
         ELSE
          CALL czfintn_or(jft  ,jlt  ,thk0 ,thk02,a_i  ,dhg     ,
     2                    x13  ,x24  ,y13  ,y24  ,z1   ,mx23    ,
     3                    mx13 ,mx34 ,my13 ,my23 ,my34 ,hourgg  ,
     4                    forg ,momg ,vf   ,vm  ,zcfac  ,a11    ,
     5                    a12  ,g    ,gs   ,sigy ,off  ,fac1    , 
     6                    rho  ,area ,  dt1,eintg,amu  ,vlhi    ,
     7                    nptt ,ipartc,partsav,kfts  ,gsr  ,
     8                    a11sr,a12sr ,nusr,shfsr ,iorth  ,hm   ,
     9                    hf   ,hc    ,hmfor,mtn  ,nel)
         ENDIF  !  IF (IORTH == 0)
        ENDIF !  (ISROT > 0) THEN
C
        IF (nptt == 1)
     1    CALL czfintnm(jft  ,jlt   ,thk0   ,a_i  ,dhg  ,
     2                  x13  ,x24   ,y13    ,y24  ,vf   ,
     3                  g    ,rho   ,area   ,amu  ,dt1  ,
     4                  off  ,ipartc,partsav,kfts )
        CALL czproj1(
     1               jft    ,jlt    ,vqn    ,vq     ,vf   ,
     2               vm     ,plat   ,
     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    ,fzero,
     8               z1     ,lxyz   ,dd     ,db     ,lxyz ,
     9               isrot  ,diz    ,vmz    )
C-------------------------
c     Shell thermal
C--------------------------
C
        IF (jthe /= 0) THEN
          IF (mat_elem%MAT_PARAM(mat(1))%HEAT_FLAG == 1) THEN
            CALL thermc(jft    ,jlt     ,pm      ,mat  ,thk0 ,ixc  ,
     .                  y24    ,px2     ,py1     ,py2  ,area ,dt1c ,
     .                  temp   ,tempel  ,fheat   ,them ,glob_therm%THEACCFACT)
          ELSE
            CALL thermc(jft    ,jlt     ,pm      ,mat  ,thk0 ,ixc  ,
     .                  y24    ,px2     ,py1     ,py2  ,area ,dt1c ,
     .                  temp   ,tempel  ,die     ,them ,glob_therm%THEACCFACT)
          ENDIF 
        ENDIF 
C
C--------------------------
C     FORCE ASSEMBLY in the sky addresses of phantom nodes
C--------------------------
        IF (iparit == 1)
     .    CALL cupdtn3_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,facn ,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