sforc3.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    sforc3                 ../engine/source/elements/solid/solide/sforc3.f
!||--- called by ------------------------------------------------------
!||    alemain                ../engine/source/ale/alemain.F
!||    forint                 ../engine/source/elements/forint.F
!||--- calls      -----------------------------------------------------
!||    aleflow                ../engine/source/ale/porous/aleflow.F
!||    alefvm_epsdot          ../engine/source/ale/alefvm/alefvm_epsdot.F
!||    alefvm_freset          ../engine/source/ale/alefvm/alefvm_freset.F
!||    alefvm_gravity         ../engine/source/ale/alefvm/alefvm_gravity.F
!||    alefvm_gravity_int22   ../engine/source/ale/alefvm/alefvm_gravity_int22.f
!||    alefvm_stress          ../engine/source/ale/alefvm/alefvm_stress.F
!||    alefvm_stress_int22    ../engine/source/ale/alefvm/alefvm_stress_int22.F
!||    alefvm_tfext           ../engine/source/ale/alefvm/alefvm_tfext.F
!||    amass3                 ../engine/source/elements/solid/solide/amass3.F
!||    amass3f                ../engine/source/elements/solid/solide/amass3f.F
!||    amass3p                ../engine/source/elements/solid/solide/amass3p.F
!||    amass3pf               ../engine/source/elements/solid/solide/amass3pf.F
!||    amomt3                 ../engine/source/elements/solid/solide/amomt3.F
!||    boltst                 ../engine/source/elements/solid/solide/boltst.F
!||    check_off_ale          ../engine/source/elements/solid/solide/check_off_ale.F
!||    deltax22               ../engine/source/interfaces/int22/deltax22.F
!||    emomt3b                ../engine/source/elements/solid/solide/emomt3b.F
!||    epxle3                 ../engine/source/elements/solid/solide/epxle3.F
!||    fcumu3                 ../engine/source/elements/solid/solide/fcumu3.F
!||    fcumu3p                ../engine/source/elements/solid/solide/fcumu3p.F
!||    fdefo3                 ../engine/source/elements/solid/solide/fdefo3.F
!||    fderi3                 ../engine/source/elements/solid/solide/fderi3.F
!||    fdlen3                 ../engine/source/elements/solid/solide/fdlen3.F
!||    fe_close               ../engine/source/elements/solid/solide/fe_close.F
!||    ffint3                 ../engine/source/elements/solid/solide/ffint3.F
!||    fhvis3                 ../engine/source/elements/solid/solide/fhvis3.F
!||    frho3                  ../engine/source/elements/solid/solide/frho3.F
!||    mmain                  ../engine/source/materials/mat_share/mmain.F90
!||    mod_close              ../engine/source/elements/solid/solide/mod_close.F
!||    s11defo3               ../engine/source/elements/solid/solide/s11defo3.F
!||    s11fx3                 ../engine/source/elements/solid/solide/s11fx3.F
!||    s8for_distor           ../engine/source/elements/solid/solide/s8for_distor.F
!||    s8get_x3               ../engine/source/elements/solid/solide/s8get_x3.F
!||    s8sav12                ../engine/source/elements/solid/solide/s8sav12.F
!||    s8sav3                 ../engine/source/elements/solid/solide/s8sav3.F
!||    s8savsr3               ../engine/source/elements/solid/solide/s8savsr3.F
!||    s8upd11t12             ../engine/source/elements/solid/solide/s8upd11t12.F
!||    s_hg5                  ../engine/source/elements/solid/solide/s_hg5.F
!||    sbilan                 ../engine/source/elements/solid/solide/sbilan.F
!||    scoor3                 ../engine/source/elements/solid/solide/scoor3.F
!||    scoor3_fvm             ../engine/source/ale/alefvm/scoor3_fvm.F
!||    scoor_cp2sp            ../engine/source/elements/solid/solidez/scoor_cp2sp.F
!||    scumu3                 ../engine/source/elements/solid/solide/scumu3.F
!||    scumu3p                ../engine/source/elements/solid/solide/scumu3p.F
!||    sdefo3                 ../engine/source/elements/solid/solide/sdefo3.F
!||    sdefot3                ../engine/source/elements/solid/solide/sdefot3.F
!||    sderi3                 ../engine/source/elements/solid/solide/sderi3.F
!||    sderit3                ../engine/source/elements/solid/solide/sderit3.F
!||    sderito3               ../engine/source/elements/solid/solide/sderito3.F
!||    sdistor_ini            ../engine/source/elements/solid/solide/sdistror_ini.F90
!||    sdlen3                 ../engine/source/elements/solid/solide/sdlen3.F
!||    sdlenmax               ../engine/source/elements/solid/solide/sdlenmax.F
!||    sfillopt               ../engine/source/elements/solid/solide/sfillopt.F
!||    sfint3                 ../engine/source/elements/solid/solide/sfint3.F
!||    sfint3b                ../engine/source/elements/solid/solide/sfint3b.F
!||    sfint_reg              ../engine/source/elements/solid/solide/sfint_reg.F
!||    sgcoor3                ../engine/source/elements/solid/solide/sgcoor3.F
!||    sgeodel3               ../engine/source/elements/solid/solide/sgeodel3.F
!||    shvis3                 ../engine/source/elements/solid/solide/shvis3.F
!||    smalla3                ../engine/source/elements/solid/solide/smalla3.F
!||    smallb3                ../engine/source/elements/solid/solide/smallb3.F
!||    smallgeo3              ../engine/source/elements/solid/solide/smallgeo3.F
!||    sordef12               ../engine/source/elements/solid/solidez/sordef12.F
!||    sordeft12              ../engine/source/elements/solid/solidez/sordeft12.F
!||    sordeft3               ../engine/source/elements/solid/solidez/sordeft3.F
!||    srbilan                ../engine/source/elements/solid/solide/srbilan.F
!||    srcoor12               ../engine/source/elements/solid/solide/srcoor12.F
!||    srcoor3                ../engine/source/elements/solid/solide/srcoor3.F
!||    srho3                  ../engine/source/elements/solid/solide/srho3.F
!||    srmalla11              ../engine/source/elements/solid/solide/srmalla11.F
!||    srmalla3               ../engine/source/elements/solid/solide/srmall3.F
!||    srota3                 ../engine/source/elements/solid/solide/srota3.F
!||    sroto12_sig            ../engine/source/elements/solid/solidez/sroto12_sig.F
!||    srrota3                ../engine/source/elements/solid/solide/srrota3.F
!||    srrotadp               ../engine/source/elements/solid/solide/srrotadp.F
!||    sstra3                 ../engine/source/elements/solid/solide/sstra3.F
!||    startime               ../engine/source/system/timer_mod.F90
!||    stherm                 ../engine/source/elements/solid/solide/stherm.F
!||    stoptime               ../engine/source/system/timer_mod.F90
!||    voln22                 ../engine/source/interfaces/int22/voln22.F
!||--- uses       -----------------------------------------------------
!||    ale_connectivity_mod   ../common_source/modules/ale/ale_connectivity_mod.F
!||    ale_mod                ../common_source/modules/ale/ale_mod.F
!||    dt_mod                 ../engine/source/modules/dt_mod.F
!||    elbufdef_mod           ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    element_mod            ../common_source/modules/elements/element_mod.F90
!||    glob_therm_mod         ../common_source/modules/mat_elem/glob_therm_mod.F90
!||    mat_elem_mod           ../common_source/modules/mat_elem/mat_elem_mod.F90
!||    mmain_mod              ../engine/source/materials/mat_share/mmain.F90
!||    nlocal_reg_mod         ../common_source/modules/nlocal_reg_mod.F
!||    output_mod             ../common_source/modules/output/output_mod.F90
!||    sdistor_ini_mod        ../engine/source/elements/solid/solide/sdistror_ini.F90
!||    sensor_mod             ../common_source/modules/sensor_mod.F90
!||    table_mod              ../engine/share/modules/table_mod.F
!||    timer_mod              ../engine/source/system/timer_mod.F90
!||====================================================================
      SUBROUTINE sforc3(TIMERS, OUTPUT,
     1   ELBUF_TAB,  NG,         PM,         GEO,
     2   IXS,        X,          NV46,       A,
     3   V,          MS,         W,          FLUX,
     4   FLU1,       VEUL,       FV,         ALE_CONNECT,
     5   IPARG,      TF,         NPF,        BUFMAT,
     6   PARTSAV,    ITAB,       DT2T,       NELTST,
     7   ITYPTST,    STIFN,      FSKY,       IADS,
     8   OFFSET,     EANI,       IPARTS,     F11,
     9   F21,        F31,        F12,        F22,
     A   F32,        F13,        F23,        F33,
     B   F14,        F24,        F34,        F15,
     C   F25,        F35,        F16,        F26,
     D   F36,        F17,        F27,        F37,
     E   F18,        F28,        F38,        NEL,
     F   FSKYM,      MSNF,       ISKY,       FSKYI,
     G   NVC,        IPM,        IGEO,       AR,
     H   VR,         IN,         FR_WAVE,    DR,
     I   BUFVOIS,    ITASK,      QMV,        ISTRAIN,
     J   TEMP,       FTHE,       FTHESKY,    IEXPAN,
     K   GRESAV,     GRTH,       IGRTH,      MSSA,
     L   DMELS,      TABLE,      PHI1,       PHI2,
     M   VF,         AF,         DF,         WF,
     N   FFSKY,      AFGLOB,     MSF,        IPARG1,
     O   XDP,        POR,        ICONTACT,   IFOAM,
     P   VOLN,       CONDN,      CONDNSKY,   AGRAV,
     Q   IGRV,       LGRAV,      SENSORS,    SKEW,
     R   NALE,       D,          IOUTPRT,    NLOC_DMG,
     S   MAT_ELEM,   H3D_STRAIN, DT,         IDEL7NOK,NSVOIS,
     T   SZ_BUFVOIS, SNPC,STF,SBUFMAT,SVIS,IDTMINS,IRESP,
     U   IDEL7NG,    MAXFUNC,    USERL_AVAIL, GLOB_THERM, 
     V   impl_s,     idyna,      WFEXT)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE timer_mod
      USE output_mod, ONLY : output_
      USE mmain_mod
      USE table_mod
      USE mat_elem_mod         
      USE nlocal_reg_mod
      USE ale_connectivity_mod
      USE sensor_mod
      USE ale_mod
      USE dt_mod
      USE elbufdef_mod
      USE sdistor_ini_mod, ONLY : sdistor_ini
      use glob_therm_mod
      use element_mod , only : nixs
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      "vect01_c.inc"
#include      "com01_c.inc"
#include      "com04_c.inc"
#include      "com08_c.inc"
#include      "scr03_c.inc"
#include      "scr06_c.inc"
#include      "parit_c.inc"
#include      "param_c.inc"
#include      "timeri_c.inc"
#include      "inter22.inc"
#include      "scr18_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
      TYPE(output_), INTENT(INOUT) :: OUTPUT
      INTEGER, INTENT(IN) :: SNPC
      INTEGER, INTENT(IN) :: SZ_BUFVOIS
      INTEGER, INTENT(IN) :: NSVOIS
      INTEGER, INTENT(IN) :: STF
      INTEGER, INTENT(IN) :: SBUFMAT
      INTEGER, INTENT(INOUT) :: IDEL7NOK
      INTEGER, INTENT(IN)    :: IDTMINS
      INTEGER ,INTENT(IN) :: IRESP
      INTEGER ,INTENT(IN) :: IDEL7NG
      INTEGER ,INTENT(IN) :: MAXFUNC
      INTEGER, INTENT(IN) :: IMPL_S
      INTEGER, INTENT(IN) :: IDYNA
      INTEGER, INTENT(IN) :: USERL_AVAIL
      INTEGER, DIMENSION(*), INTENT(IN) :: IPARG,IPARTS,NPF,NALE
      INTEGER,INTENT(IN) :: NG,NEL
      INTEGER IXS(NIXS,NUMELS), IADS(8,*),IPM(NPROPMI,NUMMAT),
     .        IGEO(NPROPGI,*),GRTH(*),IGRTH(*),IPARG1(*),IFOAM(*),
     .        ICONTACT(*),ISKY(*),ITAB(*)
      INTEGER NELTST,ITYPTST,OFFSET,NVC,ITASK,
     .        ISTRAIN, IEXPAN, NV46,IOUTPRT,H3D_STRAIN
 
      DOUBLE PRECISION :: XDP(3,*)
 
      my_real, DIMENSION(MVSIZ,6), INTENT(INOUT) :: SVIS
      my_real :: DT2T
      my_real
     .   PM(NPROPM,NUMMAT), GEO(NPROPG,NUMGEO), X(*), A(*), V(3,*), MS(*), W(*),
     .   FLUX(6,*), FSKYM(*),
     .   FLU1(*), VEUL(*), FV(*), TF(*), BUFMAT(*),
     .   partsav(*),stifn(*), fsky(*),eani(*),
     .   f11(mvsiz),f21(mvsiz),f31(mvsiz),
     .   f12(mvsiz),f22(mvsiz),f32(mvsiz),
     .   f13(mvsiz),f23(mvsiz),f33(mvsiz),
     .   f14(mvsiz),f24(mvsiz),f34(mvsiz),
     .   f15(mvsiz),f25(mvsiz),f35(mvsiz),
     .   f16(mvsiz),f26(mvsiz),f36(mvsiz),
     .   f17(mvsiz),f27(mvsiz),f37(mvsiz),
     .   f18(mvsiz),f28(mvsiz),f38(mvsiz),
     .   temp(*), fthe(*),fthesky(*),gresav(*), mssa(*), dmels(*), 
     .   phi1(*),phi2(*),vf(3,*) ,d(3,*),
     .   af(3,*), ffsky(3,*),df(3,*),wf(3,*),msf(*),afglob(3,*),
     .   por(*),voln(mvsiz), condn(*),condnsky(*)
 
      my_real :: bufvois(6,*) ! Neighbour element numbers (ALE/CFD)
      my_real :: msnf(*),ar(*),vr(3,*),in(*),fr_wave(*),dr(3,*),qmv(12,*)
      TYPE (TTABLE) TABLE(*)
      TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
      TYPE (NLOCAL_STR_)  , TARGET :: NLOC_DMG
      TYPE (MAT_ELEM_) ,INTENT(INOUT) :: MAT_ELEM
 
      my_real :: p(mvsiz),fskyi(lskyi,nfskyi)
 
      INTEGER :: IGRV(*), LGRAV(*)
      my_real :: AGRAV(*),SKEW(LSKEW,*)
      TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
      TYPE (SENSORS_) ,INTENT(INOUT) :: SENSORS
      TYPE (DT_),   INTENT(INOUT) :: DT
      type (glob_therm_) ,intent(inout)   :: glob_therm
 
      DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,NF1,IBID,IPTR,IPTS,IPTT,ILAY,ITET,IFVM_SKIP
      INTEGER NBF1,NBF2,NBF3, IALEFVM_MAT
      INTEGER MXT(MVSIZ),NGL(MVSIZ),NGEO(MVSIZ),INOD(8)
      my_real VD2(MVSIZ) , DVOL(MVSIZ),DELTAX(MVSIZ),
     .        VIS(MVSIZ) , QVIS(MVSIZ), CXX(MVSIZ) ,
     .        S1(MVSIZ)  , S2(MVSIZ)  , S3(MVSIZ)  ,
     .        S4(MVSIZ)  , S5(MVSIZ)  , S6(MVSIZ)  ,
     .        B1(MVSIZ)  , B2(MVSIZ)  , B3(MVSIZ)  ,
     .        B4(MVSIZ)  , B5(MVSIZ)  , B6(MVSIZ)  ,
     .        DXX(MVSIZ) , DYY(MVSIZ) , DZZ(MVSIZ) ,
     .        D4(MVSIZ)  , D5(MVSIZ)  , D6(MVSIZ)  ,
     .        jac1(mvsiz), jac2(mvsiz), jac3(mvsiz),
     .        jac4(mvsiz), jac5(mvsiz), jac6(mvsiz),
     .        vdx(mvsiz) , vdy(mvsiz) , vdz(mvsiz),ssp_eq(mvsiz),
     .        aire(mvsiz),fqvis(mvsiz), fcxx(mvsiz)
      my_real, dimension(mvsiz) :: fheat
C
      my_real :: sti(mvsiz), wxx(mvsiz) , wyy(mvsiz) , wzz(mvsiz)
      my_real :: muvoid(mvsiz)  ! used in case of SPH 
      my_real :: sigy(mvsiz),et(mvsiz),r3_free(mvsiz)
C
C Connectivities
      INTEGER NC1(MVSIZ), NC2(MVSIZ), NC3(MVSIZ), NC4(MVSIZ), 
     .        nc5(mvsiz), nc6(mvsiz), nc7(mvsiz), nc8(mvsiz)
C
      DOUBLE PRECISION 
     .   XD1(MVSIZ), XD2(MVSIZ), XD3(MVSIZ), XD4(MVSIZ),
     .   xd5(mvsiz), xd6(mvsiz), xd7(mvsiz), xd8(mvsiz),
     .   yd1(mvsiz), yd2(mvsiz), yd3(mvsiz), yd4(mvsiz),
     .   yd5(mvsiz), yd6(mvsiz), yd7(mvsiz), yd8(mvsiz),
     .   zd1(mvsiz), zd2(mvsiz), zd3(mvsiz), zd4(mvsiz),
     .   zd5(mvsiz), zd6(mvsiz), zd7(mvsiz), zd8(mvsiz),
     .   x0(mvsiz,8),y0(mvsiz,8),z0(mvsiz,8),voldp(mvsiz)
C
      my_real
     .   off(mvsiz) , rhoo(mvsiz),
     .   x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
     .   x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
     .   y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
     .   y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz),
     .   z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
     .   z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz),
     .  vx1(mvsiz),vx2(mvsiz),vx3(mvsiz),vx4(mvsiz),
     .  vx5(mvsiz),vx6(mvsiz),vx7(mvsiz),vx8(mvsiz),
     .  vy1(mvsiz),vy2(mvsiz),vy3(mvsiz),vy4(mvsiz),
     .  vy5(mvsiz),vy6(mvsiz),vy7(mvsiz),vy8(mvsiz),
     .  vz1(mvsiz),vz2(mvsiz),vz3(mvsiz),vz4(mvsiz),
     .  vz5(mvsiz),vz6(mvsiz),vz7(mvsiz),vz8(mvsiz),
     .  px1(mvsiz),px2(mvsiz),px3(mvsiz),px4(mvsiz),
     .  px5(mvsiz),px6(mvsiz),px7(mvsiz),px8(mvsiz),
     .  py1(mvsiz),py2(mvsiz),py3(mvsiz),py4(mvsiz),
     .  py5(mvsiz),py6(mvsiz),py7(mvsiz),py8(mvsiz),
     .  pz1(mvsiz),pz2(mvsiz),pz3(mvsiz),pz4(mvsiz),
     .  pz5(mvsiz),pz6(mvsiz),pz7(mvsiz),pz8(mvsiz),
     .  px1h1(mvsiz),px2h1(mvsiz),px3h1(mvsiz),px4h1(mvsiz),
     .  px1h2(mvsiz),px2h2(mvsiz),px3h2(mvsiz),px4h2(mvsiz),
     .  px1h3(mvsiz),px2h3(mvsiz),px3h3(mvsiz),px4h3(mvsiz),
     .  vdx1(mvsiz),vdx2(mvsiz),vdx3(mvsiz),vdx4(mvsiz),
     .  vdx5(mvsiz),vdx6(mvsiz),vdx7(mvsiz),vdx8(mvsiz),
     .  vdy1(mvsiz),vdy2(mvsiz),vdy3(mvsiz),vdy4(mvsiz),
     .  vdy5(mvsiz),vdy6(mvsiz),vdy7(mvsiz),vdy8(mvsiz),
     .  vdz1(mvsiz),vdz2(mvsiz),vdz3(mvsiz),vdz4(mvsiz),
     .  vdz5(mvsiz),vdz6(mvsiz),vdz7(mvsiz),vdz8(mvsiz),
     .  vgxa(mvsiz),vgya(mvsiz),vgza(mvsiz), vga2(mvsiz),
     .  dxy(mvsiz),dyx(mvsiz),dyz(mvsiz),
     .  dzy(mvsiz),dzx(mvsiz),dxz(mvsiz),
     .  r11(mvsiz),r12(mvsiz),r13(mvsiz),                        
     .  r21(mvsiz),r22(mvsiz),r23(mvsiz),                        
     .  r31(mvsiz),r32(mvsiz),r33(mvsiz),hh(mvsiz),gama(mvsiz,6),
     .  vx0(mvsiz,8),vy0(mvsiz,8),vz0(mvsiz,8),
     .  mfxx(mvsiz),mfxy(mvsiz),mfyx(mvsiz),
     .  mfyy(mvsiz),mfyz(mvsiz),mfzy(mvsiz),
     .  mfzz(mvsiz),mfzx(mvsiz),mfxz(mvsiz),
     .  dfe(mvsiz,3),sigf(mvsiz,6),offg0(mvsiz),
     .  xgxa(mvsiz),xgya(mvsiz),xgza(mvsiz),
     .  xgxya(mvsiz),xgyza(mvsiz),xgzxa(mvsiz),
     .  xgxa2(mvsiz),xgya2(mvsiz),xgza2(mvsiz)
C
      my_real
     . vxr(mvsiz),vyr(mvsiz),vzr(mvsiz),fr_wav(mvsiz),
     . tempel(mvsiz),them(mvsiz,8),die(mvsiz)
C
      my_real
     .  mx1(mvsiz),my1(mvsiz),mz1(mvsiz),
     .  mx2(mvsiz),my2(mvsiz),mz2(mvsiz),
     .  mx3(mvsiz),my3(mvsiz),mz3(mvsiz),
     .  mx4(mvsiz),my4(mvsiz),mz4(mvsiz),
     .  mx5(mvsiz),my5(mvsiz),mz5(mvsiz),
     .  mx6(mvsiz),my6(mvsiz),mz6(mvsiz),
     .  mx7(mvsiz),my7(mvsiz),mz7(mvsiz),
     .  mx8(mvsiz),my8(mvsiz),mz8(mvsiz)
C
      my_real
     .   n1x(mvsiz), n2x(mvsiz), n3x(mvsiz),
     .   n1y(mvsiz), n2y(mvsiz), n3y(mvsiz),
     .   n1z(mvsiz), n2z(mvsiz), n3z(mvsiz),
     .   n4x(mvsiz), n5x(mvsiz), n6x(mvsiz),
     .   n4y(mvsiz), n5y(mvsiz), n6y(mvsiz),
     .   n4z(mvsiz), n5z(mvsiz), n6z(mvsiz)
C
      my_real
     .  fvx1(mvsiz),fvx2(mvsiz),fvx3(mvsiz),fvx4(mvsiz),
     .  fvx5(mvsiz),fvx6(mvsiz),fvx7(mvsiz),fvx8(mvsiz),
     .  fvy1(mvsiz),fvy2(mvsiz),fvy3(mvsiz),fvy4(mvsiz),
     .  fvy5(mvsiz),fvy6(mvsiz),fvy7(mvsiz),fvy8(mvsiz),
     .  fvz1(mvsiz),fvz2(mvsiz),fvz3(mvsiz),fvz4(mvsiz),
     .  fvz5(mvsiz),fvz6(mvsiz),fvz7(mvsiz),fvz8(mvsiz),
     .  fvdx1(mvsiz),fvdx2(mvsiz),fvdx3(mvsiz),fvdx4(mvsiz),
     .  fvdx5(mvsiz),fvdx6(mvsiz),fvdx7(mvsiz),fvdx8(mvsiz),
     .  fvdy1(mvsiz),fvdy2(mvsiz),fvdy3(mvsiz),fvdy4(mvsiz),
     .  fvdy5(mvsiz),fvdy6(mvsiz),fvdy7(mvsiz),fvdy8(mvsiz),
     .  fvdz1(mvsiz),fvdz2(mvsiz),fvdz3(mvsiz),fvdz4(mvsiz),
     .  fvdz5(mvsiz),fvdz6(mvsiz),fvdz7(mvsiz),fvdz8(mvsiz),
     .  fdxx(mvsiz) ,fdxy(mvsiz)  ,fdxz(mvsiz)  ,fdyx(mvsiz),
     .  fdyy(mvsiz) ,fdyz(mvsiz)  ,fdzx(mvsiz)  ,fdzy(mvsiz),
     .  fdzz(mvsiz) ,fd4(mvsiz)   ,fd5(mvsiz)   ,fd6(mvsiz) ,
     .  fwxx(mvsiz) ,fwyy(mvsiz) ,fwzz(mvsiz)  ,fvxr(mvsiz),
     .  fvyr(mvsiz) ,fvzr(mvsiz) ,fvoln(mvsiz),
     .  fvd2(mvsiz)
C
      my_real
     .   ff11(mvsiz),ff21(mvsiz),ff31(mvsiz),
     .   ff12(mvsiz),ff22(mvsiz),ff32(mvsiz),
     .   ff13(mvsiz),ff23(mvsiz),ff33(mvsiz),
     .   ff14(mvsiz),ff24(mvsiz),ff34(mvsiz),
     .   ff15(mvsiz),ff25(mvsiz),ff35(mvsiz),
     .   ff16(mvsiz),ff26(mvsiz),ff36(mvsiz),
     .   ff17(mvsiz),ff27(mvsiz),ff37(mvsiz),
     .   ff18(mvsiz),ff28(mvsiz),ff38(mvsiz),
     .   fvdx(mvsiz),fvdy(mvsiz),fvdz(mvsiz),
     .   frho0(mvsiz),alpha(mvsiz),
     .   fpx1h1(mvsiz),fpx1h2(mvsiz), fpx1h3(mvsiz), 
     .   fpx2h1(mvsiz),fpx2h2(mvsiz), fpx2h3(mvsiz),
     .   fpx3h1(mvsiz),fpx3h2(mvsiz), fpx3h3(mvsiz), 
     .   fpx4h1(mvsiz),fpx4h2(mvsiz), fpx4h3(mvsiz), 
     .   fn1x(mvsiz), fn2x(mvsiz), fn3x(mvsiz),
     .   fn1y(mvsiz), fn2y(mvsiz), fn3y(mvsiz),
     .   fn1z(mvsiz), fn2z(mvsiz), fn3z(mvsiz),
     .   fn4x(mvsiz), fn5x(mvsiz), fn6x(mvsiz),
     .   fn4y(mvsiz), fn5y(mvsiz), fn6y(mvsiz),
     .   fn4z(mvsiz), fn5z(mvsiz), fn6z(mvsiz),   
     .   fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz),
     .   fx5(mvsiz), fx6(mvsiz), fx7(mvsiz), fx8(mvsiz),
     .   fy1(mvsiz), fy2(mvsiz), fy3(mvsiz), fy4(mvsiz),
     .   fy5(mvsiz), fy6(mvsiz), fy7(mvsiz), fy8(mvsiz),
     .   fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz),
     .   fz5(mvsiz), fz6(mvsiz), fz7(mvsiz), fz8(mvsiz)
C
      my_real
     .  fpx1(mvsiz),fpx2(mvsiz),fpx3(mvsiz),fpx4(mvsiz),
     .  fpy1(mvsiz),fpy2(mvsiz),fpy3(mvsiz),fpy4(mvsiz),
     .  fpz1(mvsiz),fpz2(mvsiz),fpz3(mvsiz),fpz4(mvsiz),
     .  fdeltax(mvsiz) ,flux0(mvsiz,6),flu01(mvsiz)
      my_real, DIMENSION(:), ALLOCATABLE :: var_reg
      my_real :: conde(mvsiz),amu(mvsiz),divde(mvsiz),l_max(mvsiz),
     .           sti_c(mvsiz),ll(mvsiz),fld(mvsiz),
     .           cns2,fqmax
      my_real, DIMENSION(:), POINTER :: eint
      my_real :: elem_mass
      my_real :: sum_eps(9),sum_m,sum_vol
      INTEGER :: NUM_ELEM_ALE
C
      TYPE(G_BUFEL_) ,POINTER :: GBUF
      TYPE(L_BUFEL_) ,POINTER :: LBUF     
      TYPE(BUF_MAT_) ,POINTER :: MBUF  
C
      INTEGER IALEFVM_FLG, IFVM22, IMAT, INLOC, L_NLOC, J, IPOS(8),PID,
     .         NN_DEL,ISCTL,ISTAB(MVSIZ)
C
      INTEGER IBOLTP,NBPRELD,ISM12_11  ! IBOLTP Flag Bolt Preloading
      my_real, DIMENSION(:), POINTER :: BPRELD
      my_real, DIMENSION(:), POINTER :: DNL
      INTEGER SZ_IX
C-----------------------------------------------
C   S o u r c e  L i n e s
C=======================================================================
      GBUF => elbuf_tab(ng)%GBUF
      lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
      mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
      ism12_11 = elbuf_tab(ng)%BUFLY(1)%L_SIGL
      inloc = iparg(78)
      ALLOCATE(var_reg(nel))
      tempel(:) = zero
      fheat(:)  = zero
C
      iboltp = iparg1(72) !IPARG(72,NG)
      nbpreld = gbuf%G_BPRELD
      bpreld =>gbuf%BPRELD(1:nbpreld*nel)
C
      sz_ix=numelq+numels+nsvois
C
      ibid  = 0
      nf1=nft+1
C
      IF (jcvt == 0) THEN
C Gather nodal variables (coordinates, velocities) updated lagrangian formulation
       ialefvm_mat = ipm(251,ixs(1,nf1))
       IF(ialefvm_mat<=1)THEN
C FEM velocity
         CALL scoor3(x,ixs(1,nf1),v,w,gbuf%GAMA,gama,
     .     x1, x2, x3, x4, x5, x6, x7, x8,
     .     y1, y2, y3, y4, y5, y6, y7, y8,
     .     z1, z2, z3, z4, z5, z6, z7, z8,
     .     vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8,
     .     vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8,
     .     vz1, vz2, vz3, vz4, vz5, vz6, vz7, vz8,
     .     vdx1, vdx2, vdx3, vdx4, vdx5, vdx6, vdx7, vdx8,
     .     vdy1, vdy2, vdy3, vdy4, vdy5, vdy6, vdy7, vdy8,
     .     vdz1, vdz2, vdz3, vdz4, vdz5, vdz6, vdz7, vdz8,
     .     vdx,vdy,vdz,vd2,vis,gbuf%OFF,off,gbuf%SMSTR,gbuf%RHO,
     .     rhoo,nc1,nc2,nc3,nc4,nc5,nc6,nc7,nc8,ngl,mxt,ngeo,
     .     vr,vxr,vyr,vzr,fr_wave,fr_wav,
     .     xd1, xd2, xd3, xd4, xd5, xd6, xd7, xd8,
     .     yd1, yd2, yd3, yd4, yd5, yd6, yd7, yd8,
     .     zd1, zd2, zd3, zd4, zd5, zd6, zd7, zd8,
     .     xdp,iparg, ng , nel)
       ELSE
C FVM velocity
         CALL scoor3_fvm(x,ixs(1,nf1),v,w,gama,
     .     x1, x2, x3, x4, x5, x6, x7, x8,
     .     y1, y2, y3, y4, y5, y6, y7, y8,
     .     z1, z2, z3, z4, z5, z6, z7, z8,
     .     vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8,
     .     vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8,
     .     vz1, vz2, vz3, vz4, vz5, vz6, vz7, vz8,
     .     vdx1, vdx2, vdx3, vdx4, vdx5, vdx6, vdx7, vdx8,
     .     vdy1, vdy2, vdy3, vdy4, vdy5, vdy6, vdy7, vdy8,
     .     vdz1, vdz2, vdz3, vdz4, vdz5, vdz6, vdz7, vdz8,
     .     vdx,vdy,vdz,vd2,vis,gbuf%OFF,off,gbuf%RHO,
     .     rhoo,nc1,nc2,nc3,nc4,nc5,nc6,nc7,nc8,ngl,mxt,ngeo,
     .     xd1, xd2 , xd3, xd4, xd5, xd6, xd7, xd8,
     .     yd1, yd2 , yd3, yd4, yd5, yd6, yd7, yd8,
     .     zd1, zd2 , zd3, zd4, zd5, zd6, zd7, zd8,
     .     xdp, iparg, ng , nel, gbuf%MOM,gbuf%TAG22)
       ENDIF
      ELSE
C Gather nodal variables Belytschko co-rotational formulation
       CALL srcoor3(x,ixs(1,nf1),v,w,gbuf%GAMA,gama,
     .   x1, x2, x3, x4, x5, x6, x7, x8,
     .   y1, y2, y3, y4, y5, y6, y7, y8,
     .   z1, z2, z3, z4, z5, z6, z7, z8,
     .   vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8,
     .   vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8,
     .   vz1, vz2, vz3, vz4, vz5, vz6, vz7, vz8,
     .   vd2,vis,gbuf%OFF,off,gbuf%SMSTR,gbuf%RHO,rhoo,
     .   r11, r12, r13, r21, r22, r23, r31, r32, r33,
     .   nc1,nc2,nc3,nc4,nc5,nc6,nc7,nc8,ngl,mxt,ngeo,
     .   ioutprt, vgxa, vgya, vgza, vga2,
     .   xd1, xd2, xd3, xd4, xd5, xd6, xd7, xd8,
     .   yd1, yd2, yd3, yd4, yd5, yd6, yd7, yd8,
     .   zd1, zd2, zd3, zd4, zd5, zd6, zd7, zd8,         
     .   xdp, x0 , y0 , z0 , nel, xgxa, xgya, xgza,
     .   xgxa2,xgya2  ,xgza2 ,xgxya ,xgyza ,xgzxa,iparg1,
     .   gbuf%GAMA_R) 
      ENDIF  ! JCVT == 0
!
      nn_del = 0
      pid  = ngeo(1)
      imat = mxt(1)  
!
      IF (geo(190,pid)+geo(191,pid)+geo(192,pid)+geo(192,pid)>zero)
     .        nn_del=8
      IF (nn_del ==0 .AND. dt%IDEL_BRICK>0) nn_del=8
C      
      isctl = igeo(97,pid)
C
      IF ((ismstr >= 10.AND.ismstr <= 12).AND.jlag > 0) THEN
C Gather initial coordinates total strain formulation
        CALL sgcoor3(
     1   tt,        8,         x,         ixs(1,nf1),
     2   x0,        y0,        z0,        vx0,
     3   vy0,       vz0,       gbuf%SMSTR,d,
     4   gbuf%OFF,  offg0,     nel,       xdp,
     5   mtn,       ismstr)
C
       IF (ismstr == 11) THEN
        IF (jcvt /= 0) THEN
          CALL srrotadp(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     x0(1,1), x0(1,2), x0(1,3),
     4   x0(1,4), x0(1,5), x0(1,6), x0(1,7),
     5   x0(1,8), y0(1,1), y0(1,2), y0(1,3),
     6   y0(1,4), y0(1,5), y0(1,6), y0(1,7),
     7   y0(1,8), z0(1,1), z0(1,2), z0(1,3),
     8   z0(1,4), z0(1,5), z0(1,6), z0(1,7),
     9   z0(1,8), nel)
          CALL srrota3(
     1   r11,     r12,     r13,     r21,
     2   r22,     r23,     r31,     r32,
     3   r33,     vx0(1,1),vx0(1,2),vx0(1,3),
     4   vx0(1,4),vx0(1,5),vx0(1,6),vx0(1,7),
     5   vx0(1,8),vy0(1,1),vy0(1,2),vy0(1,3),
     6   vy0(1,4),vy0(1,5),vy0(1,6),vy0(1,7),
     7   vy0(1,8),vz0(1,1),vz0(1,2),vz0(1,3),
     8   vz0(1,4),vz0(1,5),vz0(1,6),vz0(1,7),
     9   vz0(1,8),nel)
        ENDIF
C Compute shape functions partial derivatives Pij total strain formulation
        CALL sderit3(off,voln,ngl,
     .   x0(1,1), x0(1,2), x0(1,3), x0(1,4), 
     .   x0(1,5), x0(1,6), x0(1,7), x0(1,8),
     .   y0(1,1), y0(1,2), y0(1,3), y0(1,4), 
     .   y0(1,5), y0(1,6), y0(1,7), y0(1,8),
     .   z0(1,1), z0(1,2), z0(1,3), z0(1,4), 
     .   z0(1,5), z0(1,6), z0(1,7), z0(1,8),
     .   px1, px2, px3, px4,
     .   py1, py2, py3, py4,
     .   pz1, pz2, pz3, pz4,     
     .   px1h1, px1h2, px1h3,
     .   px2h1, px2h2, px2h3,
     .   px3h1, px3h2, px3h3,
     .   px4h1, px4h2, px4h3,
     .   jac1,jac2,jac3,
     .   jac4,jac5,jac6,
     .   nel,jhbe)
       ELSE
        IF (idtmin(1)==3.AND.ismstr == 12.AND.jcvt == 0.AND.ism12_11==0) THEN
C Compute Rij local reference frame for ISMSTR=12 (OFFG>1)
          CALL srcoor12(
     1   gbuf%OFF,x,       xdp,     nc1,
     2   nc2,     nc3,     nc4,     nc5,
     3   nc6,     nc7,     nc8,     r11,
     4   r12,     r13,     r21,     r22,
     5   r23,     r31,     r32,     r33,
     6   nel,     jhbe)
        END IF 
C Compute Pij in global system using inverse Jacobian matrix JAC_I for ISMSTR=10 
        CALL sderito3(off,voln,
     .    px1, px2, px3, px4,
     .    py1, py2, py3, py4,
     .    pz1, pz2, pz3, pz4,
     .    gbuf%JAC_I,nel)
       END IF
C Compute strain rates MFij total strain formulation
       CALL sdefot3(
     1   px1,     px2,     px3,     px4,
     2   py1,     py2,     py3,     py4,
     3   pz1,     pz2,     pz3,     pz4,
     4   vx0(1,1),vx0(1,2),vx0(1,3),vx0(1,4),
     5   vx0(1,5),vx0(1,6),vx0(1,7),vx0(1,8),
     6   vy0(1,1),vy0(1,2),vy0(1,3),vy0(1,4),
     7   vy0(1,5),vy0(1,6),vy0(1,7),vy0(1,8),
     8   vz0(1,1),vz0(1,2),vz0(1,3),vz0(1,4),
     9   vz0(1,5),vz0(1,6),vz0(1,7),vz0(1,8),
     a   mfxx,    mfxy,    mfxz,    mfyx,
     b   mfyy,    mfyz,    mfzx,    mfzy,
     c   mfzz,    nel)
C
       IF (jcvt /= 0 .AND.(ismstr == 10 .OR.ismstr == 12)) THEN
          CALL sordeft3(lft,llt,mfxx, mfxy, mfxz, mfyx, mfyy, mfyz,
     .         mfzx, mfzy, mfzz,
     .         r11, r21, r31, r12, r22, r32, r13, r23, r33)
       ELSEIF (idtmin(1)==3.AND.ismstr == 12.AND.jcvt==0.AND.ism12_11==0) THEN
          CALL sordeft12(lft,llt,mfxx, mfxy, mfxz,
     .         mfyx, mfyy, mfyz,
     .         mfzx, mfzy, mfzz,
     .         r11, r21, r31, r12, r22, r32, r13, r23, r33,offg0)
       END IF
      ENDIF ! End total strain formulation
C
      IF(jale+jlag /= 0)THEN
C Lagrangian or Arbitrary Euler Lagrange simulation
       IF(jclose /= 0) THEN
C ALE option /ALE/CLOSE
         CALL mod_close(
     1   geo,     ngeo,    x1,      x2,
     2   x3,      x4,      x5,      x6,
     3   x7,      x8,      y1,      y2,
     4   y3,      y4,      y5,      y6,
     5   y7,      y8,      z1,      z2,
     6   z3,      z4,      z5,      z6,
     7   z7,      z8,      hh,      xd1,
     8   xd2,     xd3,     xd4,     xd5,
     9   xd6,     xd7,     xd8,     yd1,
     a   yd2,     yd3,     yd4,     yd5,
     b   yd6,     yd7,     yd8,     zd1,
     c   zd2,     zd3,     zd4,     zd5,
     d   zd6,     zd7,     zd8,     nel)
       ENDIF
C
       IF (ismstr==11) THEN   
        CALL scoor_cp2sp(
     1   x0,      y0,      z0,      x1,
     2   x2,      x3,      x4,      x5,
     3   x6,      x7,      x8,      y1,
     4   y2,      y3,      y4,      y5,
     5   y6,      y7,      y8,      z1,
     6   z2,      z3,      z4,      z5,
     7   z6,      z7,      z8,      nel)
       ELSE
C Compute shape functions partial derivatives Pij
        CALL sderi3(
     1   off,       voln,      ngl,       xd1,
     2   xd2,       xd3,       xd4,       xd5,
     3   xd6,       xd7,       xd8,       yd1,
     4   yd2,       yd3,       yd4,       yd5,
     5   yd6,       yd7,       yd8,       zd1,
     6   zd2,       zd3,       zd4,       zd5,
     7   zd6,       zd7,       zd8,       px1,
     8   px2,       px3,       px4,       py1,
     9   py2,       py3,       py4,       pz1,
     a   pz2,       pz3,       pz4,       px1h1,
     b   px1h2,     px1h3,     px2h1,     px2h2,
     c   px2h3,     px3h1,     px3h2,     px3h3,
     d   px4h1,     px4h2,     px4h3,     jac1,
     e   jac2,      jac3,      jac4,      jac5,
     f   jac6,      gbuf%SMSTR,gbuf%OFF,  nel,
     g   voldp,     jhbe,      ismstr,    jlag)
       END IF
C Compute characteristic length
        CALL sdlen3(
     1   voln,    deltax,  x1,      x2,
     2   x3,      x4,      x5,      x6,
     3   x7,      x8,      y1,      y2,
     4   y3,      y4,      y5,      y6,
     5   y7,      y8,      z1,      z2,
     6   z3,      z4,      z5,      z6,
     7   z7,      z8,      n1x,     n2x,
     8   n3x,     n4x,     n5x,     n6x,
     9   n1y,     n2y,     n3y,     n4y,
     a   n5y,     n6y,     n1z,     n2z,
     b   n3z,     n4z,     n5z,     n6z,
     c   nel,     mtn,     jale,    jeul)
C
      ELSEIF(jeul /= 0)THEN
C Eulerian simulation
       CALL epxle3(
     1   gbuf%VOL,veul,    x1,      x2,
     2   x3,      x4,      x5,      x6,
     3   x7,      x8,      y1,      y2,
     4   y3,      y4,      y5,      y6,
     5   y7,      y8,      z1,      z2,
     6   z3,      z4,      z5,      z6,
     7   z7,      z8,      px1,     px2,
     8   px3,     px4,     py1,     py2,
     9   py3,     py4,     pz1,     pz2,
     a   pz3,     pz4,     px5,     px6,
     b   px7,     px8,     py5,     py6,
     c   py7,     py8,     pz5,     pz6,
     d   pz7,     pz8,     px1h1,   px1h2,
     e   px1h3,   px2h1,   px2h2,   px2h3,
     f   px3h1,   px3h2,   px3h3,   px4h1,
     g   px4h2,   px4h3,   voln,    deltax,
     h   n1x,     n2x,     n3x,     n4x,
     i   n5x,     n6x,     n1y,     n2y,
     j   n3y,     n4y,     n5y,     n6y,
     k   n1z,     n2z,     n3z,     n4z,
     l   n5z,     n6z,     nel,     nft,
     m   jhbe)
C
C SuperCell characteristic length for interface type 22   
       CALL deltax22(
     1   gbuf%VOL,  veul,      voln,      deltax,
     2   gbuf%TAG22,ixs,       lft,       llt,
     3   jale,      jeul)
      ENDIF
C
      IF(int22>0)THEN
C SDERI3 is updating GBUF%VOL with uncut brick volume 
C One must have polyhedra maincell volume for interface type 22
        IF(i22_aleul == 1)THEN
          CALL voln22(
     1   voln,      gbuf%TAG22,ixs,       lft,
     2   llt,       jale,      jeul)
        ENDIF
      ENDIF
C
      IF(jale+jeul > 0 .AND. iparg1(64)==1)THEN
C Compute strain rates Dij (silent boundary element)
        CALL s11defo3(
     1   pm,         v,          veul,       x,
     2   ixs,        ale_connect,dxx,        dxy,
     3   dxz,        dyx,        dyy,        dyz,
     4   dzx,        dzy,        dzz,        d4,
     5   d5,         d6,         wxx,        wyy,
     6   wzz,        bufmat,     nel,        nft,
     7   ipm,        nummat,     numels+nsvois      )
      ELSE
C Compute strain rates Dij
        CALL sdefo3(output%DATA%FANI_CELL,
     1   px1,     px2,     px3,     px4,
     2   py1,     py2,     py3,     py4,
     3   pz1,     pz2,     pz3,     pz4,
     4   px5,     px6,     px7,     px8,
     5   py5,     py6,     py7,     py8,
     6   pz5,     pz6,     pz7,     pz8,
     7   vx1,     vx2,     vx3,     vx4,
     8   vx5,     vx6,     vx7,     vx8,
     9   vy1,     vy2,     vy3,     vy4,
     a   vy5,     vy6,     vy7,     vy8,
     b   vz1,     vz2,     vz3,     vz4,
     c   vz5,     vz6,     vz7,     vz8,
     d   dxx,     dxy,     dxz,     dyx,
     e   dyy,     dyz,     dzx,     dzy,
     f   dzz,     d4,      d5,      d6,
     g   wxx,     wyy,     wzz,     vxr,
     h   vyr,     vzr,     nel,
     i   nft,     ismstr,  jeul,    jhbe,
     j   jcvt,    isrot)
C
        IF (idtmin(1)==3.AND.ismstr == 12.AND.ism12_11==0.AND.jcvt==0) THEN
          CALL sordef12(lft,llt,dxx, dyy, dzz,
     .         d4, d5, d6,
     .         r11, r21, r31, r12, r22, r32, r13, r23, r33,gbuf%OFF)
        ENDIF    
      ENDIF 
C
      CALL s11fx3(
     1   pm,         flux(1,nf1),ale_connect,ixs,
     2   ipm,        bufmat,     nel,        nft,
     3   jale,       jeul  ,     nummat,     numels+nsvois)
C
      divde(1:nel) = dt1*(dxx(1:nel)+ dyy(1:nel)+ dzz(1:nel))  
C Compute new density
      CALL srho3(
     1   pm,         gbuf%VOL,   gbuf%RHO,   gbuf%EINT,
     2   divde,      flux(1,nf1),flu1(nf1),  voln,
     3   dvol,       ngl,        mxt,        off,
     4   iparg1(64), gbuf%TAG22, voldp,      lbuf%VOL0DP,
     5   amu,        gbuf%OFF,   nel,        mtn,
     6   jale,       ismstr,     jeul,       jlag)
C
      IF (ismstr == 12.AND.idtmin(1)==3.AND.ism12_11==0.AND.jcvt==0) THEN
         CALL sroto12_sig(lft,llt,lbuf%SIG,nel,
     .           r11, r21, r31, r12, r22, r32, r13, r23, r33,offg0)
      ENDIF
C
      IF (jcvt == 0) THEN
C Jaumann objective stress tensor
        CALL srota3(
     1   gbuf%SIG,s1,      s2,      s3,
     2   s4,      s5,      s6,      wxx,
     3   wyy,     wzz,     nel,     mtn,
     4   ismstr)
C Small Strain formulation
        CALL smalla3(
     1   gbuf%SMSTR,gbuf%OFF,  off,       wxx,
     2   wyy,       wzz,       nel,       ismstr,
     3   jlag)
      ELSE  
C Stress tensor co-rotational formulation
        CALL srmalla3(gbuf%SIG,s1,s2,s3,s4,s5,s6,
     .                gbuf%OFF,off,nel)
C Small strain formulation
        CALL srmalla11(
     1   gbuf%SMSTR,gbuf%OFF,  wxx,       wyy,
     2   wzz,       r11,       r12,       r13,
     3   r21,       r22,       r23,       r31,
     4   r32,       r33,       nel,       ismstr)
      ENDIF
C
      IF (ismstr <= 3.OR.(ismstr==4.AND.jlag>0)) THEN
C Update the reference configuration (possible future change to small strain option)
C Total strain option doesn't change the reference configuration
       IF (jcvt /= 0.AND.isorth/=0 ) THEN
         CALL s8savsr3(
     1   gbuf%OFF,  gbuf%SMSTR,x0,        y0,
     2   z0,        nel)
       ELSE
         CALL s8sav3(
     1   gbuf%OFF,  gbuf%SMSTR,xd1,       xd2,
     2   xd3,       xd4,       xd5,       xd6,
     3   xd7,       xd8,       yd1,       yd2,
     4   yd3,       yd4,       yd5,       yd6,
     5   yd7,       yd8,       zd1,       zd2,
     6   zd3,       zd4,       zd5,       zd6,
     7   zd7,       zd8,       nel)
       ENDIF
      END IF
C
      IF (jthe < 0) THEN 
C Element Temperature
        DO i=1,nel
          tempel(i) = one_over_8 *(  temp(nc1(i)) + temp(nc2(i))  
     .                             + temp(nc3(i)) + temp(nc4(i)) 
     .                             + temp(nc5(i)) + temp(nc6(i)) 
     .                             + temp(nc7(i)) + temp(nc8(i)))       
        ENDDO
      ENDIF
C
c===========================================================================
c     Computation of the non-local variable at Gauss point from nodal d.o.fs
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=1,nel
          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))
          inod(5) = nloc_dmg%IDXI(nc5(i))
          inod(6) = nloc_dmg%IDXI(nc6(i))
          inod(7) = nloc_dmg%IDXI(nc7(i))
          inod(8) = nloc_dmg%IDXI(nc8(i))
          DO j = 1, 8
            ipos(j) = nloc_dmg%POSI(inod(j))
          ENDDO
          var_reg(i) = dnl(ipos(1)) + dnl(ipos(2)) + dnl(ipos(3)) + dnl(ipos(4)) 
     .               + dnl(ipos(5)) + dnl(ipos(6)) + dnl(ipos(7)) + dnl(ipos(8))
          var_reg(i) = var_reg(i)*one_over_8
        ENDDO
      ENDIF
c===========================================================================
C
      IF (jale == 3 .AND. jlag > 0) THEN
C Law 77 internal gas flow of open cell solved with ALE formulation
C AFGLOB ----> ZERO is not activated
        CALL aleflow (
     .     pm     ,mxt   ,x      ,
     .     mbuf   ,nel   ,
     .     ixs    ,ale_connect   ,vf     ,afglob,
     .     wf     ,dfe   ,frho0  ,alpha  ,                  
     .     nc1    ,nc2   ,nc3    ,nc4    ,nc5    ,nc6   ,              
     .     nc7    ,nc8   ,fx1    ,fx2    ,fx3    ,fx4   ,              
     .     fx5    ,fx6   ,fx7    ,fx8    ,fy1    ,fy2   ,              
     .     fy3    ,fy4   ,fy5    ,fy6    ,fy7    ,fy8   ,              
     .     fz1    ,fz2   ,fz3    ,fz4    ,fz5    ,fz6   ,              
     .     fz7    ,fz8   ,phi1   ,phi2   ,flux0  ,flu01 ,   
     .     fvx1   ,fvx2  ,fvx3   ,fvx4   ,fvx5   ,fvx6  ,              
     .     fvx7   ,fvx8  ,fvy1   ,fvy2   ,fvy3   ,fvy4  ,              
     .     fvy5   ,fvy6  ,fvy7   ,fvy8   ,fvz1   ,fvz2  ,              
     .     fvz3   ,fvz4  ,fvz5   ,fvz6   ,fvz7   ,fvz8  ,              
     .     fvdx1  ,fvdx2 ,fvdx3  ,fvdx4  ,fvdx5  ,fvdx6 ,              
     .     fvdx7  ,fvdx8 ,fvdy1  ,fvdy2  ,fvdy3  ,fvdy4 ,              
     .     fvdy5  ,fvdy6 ,fvdy7  ,fvdy8  ,fvdz1  ,fvdz2 ,              
     .     fvdz3  ,fvdz4 ,fvdz5  ,fvdz6  ,fvdz7  ,fvdz8 ,              
     .     fvdx   ,fvdy  ,fvdz   ,fvd2   ,
     .     vx1    ,vx2   ,vx3    ,vx4    ,vx5    ,vx6   ,              
     .     vx7    ,vx8   ,vy1    ,vy2    ,vy3    ,vy4   ,              
     .     vy5    ,vy6   ,vy7    ,vy8    ,vz1    ,vz2   ,              
     .     vz3    ,vz4   ,vz5    ,vz6    ,vz7    ,vz8   ,
     .     por   ,icontact ,ifoam)    
C
        CALL fderi3(
     1   off,     fvoln,   ngl,     fx1,
     2   fx2,     fx3,     fx4,     fx5,
     3   fx6,     fx7,     fx8,     fy1,
     4   fy2,     fy3,     fy4,     fy5,
     5   fy6,     fy7,     fy8,     fz1,
     6   fz2,     fz3,     fz4,     fz5,
     7   fz6,     fz7,     fz8,     fpx1,
     8   fpx2,    fpx3,    fpx4,    fpy1,
     9   fpy2,    fpy3,    fpy4,    fpz1,
     a   fpz2,    fpz3,    fpz4,    fpx1h1,
     b   fpx1h2,  fpx1h3,  fpx2h1,  fpx2h2,
     c   fpx2h3,  fpx3h1,  fpx3h2,  fpx3h3,
     d   fpx4h1,  fpx4h2,  fpx4h3,  nel,
     e   jhbe)
C
        CALL fdlen3(
     1   fvoln,   fdeltax, fx1,     fx2,
     2   fx3,     fx4,     fx5,     fx6,
     3   fx7,     fx8,     fy1,     fy2,
     4   fy3,     fy4,     fy5,     fy6,
     5   fy7,     fy8,     fz1,     fz2,
     6   fz3,     fz4,     fz5,     fz6,
     7   fz7,     fz8,     fn1x,    fn2x,
     8   fn3x,    fn4x,    fn5x,    fn6x,
     9   fn1y,    fn2y,    fn3y,    fn4y,
     a   fn5y,    fn6y,    fn1z,    fn2z,
     b   fn3z,    fn4z,    fn5z,    fn6z,
     c   nel)
C
C Strain flow    
        CALL fdefo3(
     1   fpx1,    fpx2,    fpx3,    fpx4,
     2   fpy1,    fpy2,    fpy3,    fpy4,
     3   fpz1,    fpz2,    fpz3,    fpz4,
     4   fvx1,    fvx2,    fvx3,    fvx4,
     5   fvx5,    fvx6,    fvx7,    fvx8,
     6   fvy1,    fvy2,    fvy3,    fvy4,
     7   fvy5,    fvy6,    fvy7,    fvy8,
     8   fvz1,    fvz2,    fvz3,    fvz4,
     9   fvz5,    fvz6,    fvz7,    fvz8,
     a   fdxx,    fdxy,    fdxz,    fdyx,
     b   fdyy,    fdyz,    fdzx,    fdzy,
     c   fdzz,    fd4,     fd5,     fd6,
     d   fwxx,    fwyy,    fwzz,    fvxr,
     e   fvyr,    fvzr,    fqvis,   nel,
     f   jhbe,    jcvt,    isrot)
C
        DO i=1,nel                                   
           fvoln(i) = alpha(i)*fvoln(i)                
           dfe(i,1) = dfe(i,1)*fvoln(i)                
           dfe(i,2) = dfe(i,2)*fvoln(i)                
           dfe(i,3) = dfe(i,3)*fvoln(i)                                                                      
        ENDDO 
C
        CALL frho3(mbuf%VAR,frho0,fdxx, fdyy, fdzz,    
     .             flux0,flu01,fvoln,ngl,mxt,off ,nel) 
      ENDIF    ! JALE == 3 .AND. JLAG > 0
C
      ialefvm_mat = ipm(251,ixs(1,nf1))
      IF (ialefvm_mat  > 1) THEN
C ALEFVM strain rates
          CALL alefvm_epsdot(
     1   ixs,
     2   voln,    n1x,     n2x,     n3x,
     3   n4x,     n5x,     n6x,     n1y,
     4   n2y,     n3y,     n4y,     n5y,
     5   n6y,     n1z,     n2z,     n3z,
     6   n4z,     n5z,     n6z,     dxx,
     7   dyy,     dzz,     x,       ipm,
     8   nel)
      ENDIF
C
      IF ((imon_mat==1).AND. itask == 0) CALL startime(timers,35)
C
      ilay = 1                                                
      iptr = 1                                                
      ipts = 1                                                
      iptt = 1   
C
      IF(iboltp /= 0) THEN
        IF(jcvt==0)
     .          CALL srota3(
     1   bpreld(3*nel+1),b1,             b2,             b3,
     2   b4,             b5,             b6,             wxx,
     3   wyy,            wzz,            nel,            mtn,
     4   ismstr)
 
        CALL boltst(iptr  ,bpreld,lbuf%SIG       ,tt    ,nel   ,        
     .              npt   ,sensors%NSENSOR,sensors%SENSOR_TAB ,
     .            iparg1(67),iparg1(68))
      END IF
C
C Stress computation according to material laws 
C
      CALL mmain(timers, output,
     1   elbuf_tab,   ng,          pm,          geo,
     2   ale_connect, ixs,         iparg,
     3   v,           tf,          npf,         bufmat,
     4   sti,         x,           dt2t,        neltst,
     5   ityptst,     offset,      nel,         w,
     6   off,         ngeo,        mxt,         ngl,
     7   voln,        vd2,         dvol,        deltax,
     8   vis,         qvis,        cxx,         s1,
     9   s2,          s3,          s4,          s5,
     a   s6,          dxx,         dyy,         dzz,
     b   d4,          d5,          d6,          wxx,
     c   wyy,         wzz,         jac1,        jac2,
     d   jac3,        jac4,        jac5,        jac6,
     e   vdx,         vdy,         vdz,         muvoid,
     f   ssp_eq,      aire,        sigy,        et,
     g   bufvois,     lbuf%PLA,    r3_free,     amu,
     h   mfxx,        mfxy,        mfxz,        mfyx,
     i   mfyy,        mfyz,        mfzx,        mfzy,
     j   mfzz,        ipm,         gama,        fr_wav,
     k   dxy,         dyx,         dyz,         dzy,
     l   dzx,         dxz,         istrain,     tempel,
     m   die,         iexpan,      ilay,        mssa,
     n   dmels,       iptr,        ipts,        iptt,
     o   table,       fvd2,        fdeltax,     fcxx,
     p   fqvis,       iparg1,      igeo,        conde,
     q   itask,       nloc_dmg,    var_reg,     mat_elem,
     r   h3d_strain,  jplasol,     jsph,        sz_bufvois,
     t   snpc,        stf,         sbufmat,     glob_therm,
     .   svis,        sz_ix,       iresp,       
     .   n2d,         th_strain,   ngroup,      tt,
     .   dt1,         ntable,      numelq,      nummat,
     .   numgeo,      numnod,      numels,            
     .   idel7nok,    idtmin,      maxfunc,
     .   imon_mat,    userl_avail, impl_s,
     .   idyna,       dt         , fheat     ,sensors)
C
      IF (ismstr == 12.AND.idtmin(1)==3) THEN
        IF (ism12_11==0.AND.jcvt==0) THEN
          CALL sroto12_sig(lft,llt,lbuf%SIG,nel,
     .           r11, r12, r13, r21, r22, r23, r31, r32, r33,offg0)
C
          IF (istrain == 1) THEN 
            CALL sordef12(lft,llt,dxx, dxy, dxz,
     .         d4, d5, d6,
     .         r11, r12, r13, r21, r22, r23, r31, r32, r33,offg0)
          ENDIF
        END IF 
      ENDIF ! ISMSTR == 12.AND.IDTMIN(1)==3
C
      IF (istrain == 1) THEN 
C Compute Strains for output
        CALL sstra3(
     1   dxx,      dyy,      dzz,      d4,
     2   d5,       d6,       lbuf%STRA,wxx,
     3   wyy,      wzz,      off,      nel,
     4   jcvt)
      ENDIF
      IF ((imon_mat==1).AND. itask == 0) CALL stoptime(timers,35)
C
      IF (nn_del >0) THEN
        CALL sdlenmax(l_max, 
     1   x1,      x2,      x3,      x4,      
     2   x5,      x6,      x7,      x8,      
     3   y1,      y2,      y3,      y4,      
     4   y5,      y6,      y7,      y8,      
     5   z1,      z2,      z3,      z4,      
     6   z5,      z6,      z7,      z8,      
     7   nel)
       CALL sgeodel3(ngl,gbuf%OFF,voln,deltax,gbuf%VOL,geo(1,ngeo(1)),l_max,dt,nel,idel7nok )
      END IF
C Small Strain
      CALL smallb3(
     1   gbuf%OFF,off,     nel,     ismstr)
      itet = 0
      CALL smallgeo3(
     1   ngl,     gbuf%OFF,voln,    deltax,
     2   gbuf%VOL,itet,    nel,     ismstr,dt )
      IF (ismstr == 12.AND.idtmin(1)==3) THEN
        CALL s8sav12(
     1   gbuf%OFF,  offg0,     gbuf%SMSTR,x,
     2   xdp,       nc1,       nc2,       nc3,
     3   nc4,       nc5,       nc6,       nc7,
     4   nc8,       nel,       jcvt)
        IF (ism12_11>0 .AND. isorth == 0) THEN
          CALL s8upd11t12(gbuf%OFF,offg0  ,x  ,xdp  ,
     .     nc1   ,nc2   ,nc3   ,nc4   ,nc5   ,nc6   ,nc7   ,nc8   ,
     .     gbuf%JAC_I,gbuf%SIG,lbuf%SIGL  ,nel   )
        END IF
      ENDIF
C
      IF (jlag+jale+jeul == 0) THEN
C Energy balance thermal material
        IF (ioutprt>0)THEN
          IF (iparg1(64) == 1) THEN                                    
            eint => elbuf_tab(ng)%GBUF%EINS(1:nel)                   
          ELSE                                                   
            eint => elbuf_tab(ng)%GBUF%EINT(1:nel)                   
          ENDIF                                                  
C
          IF (jcvt == 0) THEN  
            ifvm22      = 0
            ialefvm_flg = ipm(251,imat)
            IF(ialefvm_flg>=2)ifvm22 = 1
            CALL sbilan(partsav,eint,gbuf%RHO,gbuf%RK,gbuf%VOL,  
     .                  vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8,
     .                  vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8,
     .                  vz1, vz2, vz3, vz4, vz5, vz6, vz7, vz8,
     .                  voln,iparts,gresav,grth,igrth, gbuf%OFF,
     .                  iexpan,gbuf%EINTTH,gbuf%FILL,gbuf%MOM,nel,ifvm22,
     .                  x1, x2, x3, x4, x5, x6, x7, x8,
     .                  y1, y2, y3, y4, y5, y6, y7, y8,
     .                  z1, z2, z3, z4, z5, z6, z7, z8,
     .                  itask,iparg1,sensors,gbuf%G_WPLA,gbuf%WPLA)
          ELSE
            CALL srbilan(partsav,eint,gbuf%RHO,gbuf%RK,gbuf%VOL, 
     .                   vgxa, vgya, vgza, vga2, voln,iparts,
     .                   gresav,grth,igrth,gbuf%OFF,iexpan,gbuf%EINTTH,
     .                   gbuf%FILL,xgxa, xgya, xgza,
     .                   xgxa2,xgya2,xgza2,xgxya,xgyza,xgzxa,itask,iparg1,sensors,
     .                   nel,gbuf%G_WPLA,gbuf%WPLA)
          ENDIF
        ENDIF
        RETURN
      ENDIF
C
C Mass update
      IF (jale == 3 .AND. jlag == 1 .AND. ale%GLOBAL%INCOMP == 0) THEN  ! law 77 
       IF(iparit == 0)THEN         
          CALL amass3f(
     1   msf,     mbuf%VAR,fvoln,   nc1,
     2   nc2,     nc3,     nc4,     nc5,
     3   nc6,     nc7,     nc8,     msnf,
     4   nvc,     off,     nel)
        ELSE
          CALL amass3pf(
     1   fskym,   mbuf%VAR,fvoln,   iads,
     2   off,     nel,     nft)
        ENDIF
      ELSEIF (jale+jeul > 0 .AND. ale%GLOBAL%INCOMP == 0)THEN
        IF(iparit == 0)THEN
          CALL amass3(
     1   ms,                gbuf%RHO,          veul(lveul*nft+44),gbuf%TAG22,
     2   voln,              nc1,               nc2,               nc3,
     3   nc4,               nc5,               nc6,               nc7,
     4   nc8,               msnf,              nvc,               off,
     5   ixs,               nel,               jeul)
        ELSE
          CALL amass3p(
     1   fskym,             gbuf%RHO,          veul(lveul*nft+44),gbuf%TAG22,
     2   voln,              iads,              off,               ixs,
     3   nel,               nft,               jeul)
        ENDIF
      ENDIF
C
      ialefvm_mat = ipm(251,ixs(1,nf1))
      ifvm_skip=0
      IF(ialefvm_mat > 1)ifvm_skip=1
 
      !---------------------------------------------
      ! /ALE/GRID/FLOW-TRACKING - BUFFER UPDATED
      !---------------------------------------------
      ! this grid formulation needs an averaged tensor (E=GRAD U)
      IF(ale%GRID%NWALE == 7 .AND. jale == 1)THEN
        sum_eps(1:9) = zero
        sum_m = zero
        sum_vol = zero
        num_elem_ale = nel
        DO i=1,nel
          elem_mass = gbuf%RHO(i)*gbuf%VOL(i)
          sum_eps(1) = sum_eps(1) + elem_mass*dxx(i)
          sum_eps(2) = sum_eps(2) + elem_mass*dyy(i)
          sum_eps(3) = sum_eps(3) + elem_mass*dzz(i)
          sum_eps(4) = sum_eps(4) + elem_mass*dxy(i)
          sum_eps(5) = sum_eps(5) + elem_mass*dxz(i)
          sum_eps(6) = sum_eps(6) + elem_mass*dyz(i)
          sum_eps(7) = sum_eps(7) + elem_mass*dyx(i)
          sum_eps(8) = sum_eps(8) + elem_mass*dzx(i)
          sum_eps(9) = sum_eps(9) + elem_mass*dzy(i)
          sum_m = sum_m + elem_mass
          sum_vol = sum_vol + gbuf%VOL(i)
        ENDDO
#include "lockon.inc"
      ale%GRID%flow_tracking_data%EP(1) = ale%GRID%flow_tracking_data%EP(1) + sum_eps(1)
      ale%GRID%flow_tracking_data%EP(2) = ale%GRID%flow_tracking_data%EP(2) + sum_eps(2)
      ale%GRID%flow_tracking_data%EP(3) = ale%GRID%flow_tracking_data%EP(3) + sum_eps(3)
      ale%GRID%flow_tracking_data%EP(4) = ale%GRID%flow_tracking_data%EP(4) + sum_eps(4)
      ale%GRID%flow_tracking_data%EP(5) = ale%GRID%flow_tracking_data%EP(5) + sum_eps(5)
      ale%GRID%flow_tracking_data%EP(6) = ale%GRID%flow_tracking_data%EP(6) + sum_eps(6)
      ale%GRID%flow_tracking_data%EP(7) = ale%GRID%flow_tracking_data%EP(7) + sum_eps(7)
      ale%GRID%flow_tracking_data%EP(8) = ale%GRID%flow_tracking_data%EP(8) + sum_eps(8)
      ale%GRID%flow_tracking_data%EP(9) = ale%GRID%flow_tracking_data%EP(9) + sum_eps(9)
      ale%GRID%flow_tracking_data%SUM_M = ale%GRID%flow_tracking_data%SUM_M + sum_m
      ale%GRID%flow_tracking_data%SUM_VOL = ale%GRID%flow_tracking_data%SUM_VOL + sum_vol
      ale%GRID%flow_tracking_data%NUM_ELEM_ALE = ale%GRID%flow_tracking_data%NUM_ELEM_ALE + num_elem_ale
 
#include "lockoff.inc"
      endif!(ALE%GRID%NWALE == 7)
 
C
C Anti hourglass forces
      IF (iint==3) THEN
         CALL s_hg5(
     1   pm, gbuf%RHO,     off,     vx1,
     2   vx2,     vx3,     vx4,     vx5,
     3   vx6,     vx7,     vx8,     vy1,
     4   vy2,     vy3,     vy4,     vy5,
     5   vy6,     vy7,     vy8,     vz1,
     6   vz2,     vz3,     vz4,     vz5,
     7   vz6,     vz7,     vz8,     f11,
     8   f21,     f31,     f12,     f22,
     9   f32,     f13,     f23,     f33,
     a   f14,     f24,     f34,     f15,
     b   f25,     f35,     f16,     f26,
     c   f36,     f17,     f27,     f37,
     d   f18,     f28,     f38,     px1h1,
     e   px1h2,   px1h3,   px2h1,   px2h2,
     f   px2h3,   px3h1,   px3h2,   px3h3,
     g   px4h1,   px4h2,   px4h3,   voln,
     i   gbuf%HOURG,partsav,iparts, nft,
     j   mtn,     jlag,    eani   , ngeo ,
     k   ivector ,dt1 ,    mxt,     geo,
     l   nel,     isctl,   sti,   cxx  )
      ELSE
        CALL shvis3(
     1   pm,       gbuf%RHO, off,      vx1,
     2   vx2,      vx3,      vx4,      vx5,
     3   vx6,      vx7,      vx8,      vy1,
     4   vy2,      vy3,      vy4,      vy5,
     5   vy6,      vy7,      vy8,      vz1,
     6   vz2,      vz3,      vz4,      vz5,
     7   vz6,      vz7,      vz8,      f11,
     8   f21,      f31,      f12,      f22,
     9   f32,      f13,      f23,      f33,
     a   f14,      f24,      f34,      f15,
     b   f25,      f35,      f16,      f26,
     c   f36,      f17,      f27,      f37,
     d   f18,      f28,      f38,      px1h1,
     e   px1h2,    px1h3,    px2h1,    px2h2,
     f   px2h3,    px3h1,    px3h2,    px3h3,
     g   px4h1,    px4h2,    px4h3,    voln,
     h   mxt,      cxx,      vis,      vd2,
     i   deltax,   eani,     ngeo,     geo,
     j   partsav,  iparts,   gbuf%OFF, gbuf%VOL,
     k   iparg1,   ifvm_skip,nel,      nft,
     l   mtn,      ismstr,   jlag,     jhbe)
      END IF !(IINT==3) THEN
C
      IF (jale == 3) THEN   ! law 77
        CALL fhvis3(
     1   pm,      mbuf%VAR,off,     fvx1,
     2   fvx2,    fvx3,    fvx4,    fvx5,
     3   fvx6,    fvx7,    fvx8,    fvy1,
     4   fvy2,    fvy3,    fvy4,    fvy5,
     5   fvy6,    fvy7,    fvy8,    fvz1,
     6   fvz2,    fvz3,    fvz4,    fvz5,
     7   fvz6,    fvz7,    fvz8,    ff11,
     8   ff21,    ff31,    ff12,    ff22,
     9   ff32,    ff13,    ff23,    ff33,
     a   ff14,    ff24,    ff34,    ff15,
     b   ff25,    ff35,    ff16,    ff26,
     c   ff36,    ff17,    ff27,    ff37,
     d   ff18,    ff28,    ff38,    fpx1h1,
     e   fpx1h2,  fpx1h3,  fpx2h1,  fpx2h2,
     f   fpx2h3,  fpx3h1,  fpx3h2,  fpx3h3,
     g   fpx4h1,  fpx4h2,  fpx4h3,  fvoln,
     h   mxt,     fcxx,    vis,     fvd2,
     i   fdeltax, eani,    ngeo,    geo,
     j   partsav, iparts,  nel,     nft,
     k   jlag,    jhbe)
C
      ENDIF 
C
      IF (jclose /= 0) 
     .  CALL fe_close(
     1   hh,      gbuf%RHO,voln,    vx1,
     2   vx2,     vx3,     vx4,     vx5,
     3   vx6,     vx7,     vx8,     vy1,
     4   vy2,     vy3,     vy4,     vy5,
     5   vy6,     vy7,     vy8,     vz1,
     6   vz2,     vz3,     vz4,     vz5,
     7   vz6,     vz7,     vz8,     f11,
     8   f21,     f31,     f12,     f22,
     9   f32,     f13,     f23,     f33,
     a   f14,     f24,     f34,     f15,
     b   f25,     f35,     f16,     f26,
     c   f36,     f17,     f27,     f37,
     d   f18,     f28,     f38,     nel)
C
      ialefvm_mat = ipm(251,ixs(1,nf1))
      IF(ialefvm_mat <= 1)THEN
       IF (jale == 1 .OR. (jeul == 1.AND.integ8 == 0)) THEN
C ALE & Euler momentum transport forces 
        CALL amomt3(
     1   pm,        gbuf%RHO,  voln,      x1,
     2   x2,        x3,        x4,        x5,
     3   x6,        x7,        x8,        y1,
     4   y2,        y3,        y4,        y5,
     5   y6,        y7,        y8,        z1,
     6   z2,        z3,        z4,        z5,
     7   z6,        z7,        z8,        vx1,
     8   vx2,       vx3,       vx4,       vx5,
     9   vx6,       vx7,       vx8,       vy1,
     a   vy2,       vy3,       vy4,       vy5,
     b   vy6,       vy7,       vy8,       vz1,
     c   vz2,       vz3,       vz4,       vz5,
     d   vz6,       vz7,       vz8,       f11,
     e   f21,       f31,       f12,       f22,
     f   f32,       f13,       f23,       f33,
     g   f14,       f24,       f34,       f15,
     h   f25,       f35,       f16,       f26,
     i   f36,       f17,       f27,       f37,
     j   f18,       f28,       f38,       px1,
     k   px2,       px3,       px4,       py1,
     l   py2,       py3,       py4,       pz1,
     m   pz2,       pz3,       pz4,       dxx,
     n   dxy,       dxz,       dyx,       dyy,
     o   dyz,       dzx,       dzy,       dzz,
     p   vdx1,      vdx2,      vdx3,      vdx4,
     q   vdx5,      vdx6,      vdx7,      vdx8,
     r   vdy1,      vdy2,      vdy3,      vdy4,
     s   vdy5,      vdy6,      vdy7,      vdy8,
     t   vdz1,      vdz2,      vdz3,      vdz4,
     u   vdz5,      vdz6,      vdz7,      vdz8,
     v   vdx,       vdy,       vdz,       deltax,
     w   vis,       mxt,       qmv,       bufmat,
     x   iparg1,    ixs,       gbuf%TAG22,nc1,
     y   nc2,       nc3,       nc4,       nc5,
     z   nc6,       nc7,       nc8,       nale,
     1   nel,       nft,       mtn,       ipm)
C
       ELSEIF (jeul == 1 .AND. integ8 == 1) THEN
C SUPG and TG UPWIND not included 
        CALL emomt3b(
     1   pm,                gbuf%RHO,          veul(lveul*nft+44),f11,
     2   f21,               f31,               f12,               f22,
     3   f32,               f13,               f23,               f33,
     4   f14,               f24,               f34,               f15,
     5   f25,               f35,               f16,               f26,
     6   f36,               f17,               f27,               f37,
     7   f18,               f28,               f38,               px1,
     8   px2,               px3,               px4,               py1,
     9   py2,               py3,               py4,               pz1,
     a   pz2,               pz3,               pz4,               px5,
     b   px6,               px7,               px8,               py5,
     c   py6,               py7,               py8,               pz5,
     d   pz6,               pz7,               pz8,               dxx,
     e   dxy,               dxz,               dyx,               dyy,
     f   dyz,               dzx,               dzy,               dzz,
     g   vdx,               vdy,               vdz,               mxt,
     h   qmv,               bufmat,            vx1,               vx2,
     i   vx3,               vx4,               vx5,               vx6,
     j   vx7,               vx8,               vy1,               vy2,
     k   vy3,               vy4,               vy5,               vy6,
     l   vy7,               vy8,               vz1,               vz2,
     m   vz3,               vz4,               vz5,               vz6,
     n   vz7,               vz8,               iparg1,            nel,
     o   mtn,               ipm)
       ENDIF
C
        IF(jeul+jale/=0) CALL check_off_ale(f11,f21,f31,f12,f22,
     1                           f32,f13,f23,f33,f14,
     2                           f24,f34,f15,f25,f35,
     3                           f16,f26,f36,f17,f27,
     4                           f37,f18,f28,f38,gbuf%OFF,
     5                           lft,llt,nel)
C
       IF (jale  == 3 )THEN       
         CALL amomt3(
     1   pm,        mbuf%VAR,  fvoln,     fx1,
     2   fx2,       fx3,       fx4,       fx5,
     3   fx6,       fx7,       fx8,       fy1,
     4   fy2,       fy3,       fy4,       fy5,
     5   fy6,       fy7,       fy8,       fz1,
     6   fz2,       fz3,       fz4,       fz5,
     7   fz6,       fz7,       fz8,       fvx1,
     8   fvx2,      fvx3,      fvx4,      fvx5,
     9   fvx6,      fvx7,      fvx8,      fvy1,
     a   fvy2,      fvy3,      fvy4,      fvy5,
     b   fvy6,      fvy7,      fvy8,      fvz1,
     c   fvz2,      fvz3,      fvz4,      fvz5,
     d   fvz6,      fvz7,      fvz8,      ff11,
     e   ff21,      ff31,      ff12,      ff22,
     f   ff32,      ff13,      ff23,      ff33,
     g   ff14,      ff24,      ff34,      ff15,
     h   ff25,      ff35,      ff16,      ff26,
     i   ff36,      ff17,      ff27,      ff37,
     j   ff18,      ff28,      ff38,      fpx1,
     k   fpx2,      fpx3,      fpx4,      fpy1,
     l   fpy2,      fpy3,      fpy4,      fpz1,
     m   fpz2,      fpz3,      fpz4,      fdxx,
     n   fdxy,      fdxz,      fdyx,      fdyy,
     o   fdyz,      fdzx,      fdzy,      fdzz,
     p   fvdx1,     fvdx2,     fvdx3,     fvdx4,
     q   fvdx5,     fvdx6,     fvdx7,     fvdx8,
     r   fvdy1,     fvdy2,     fvdy3,     fvdy4,
     s   fvdy5,     fvdy6,     fvdy7,     fvdy8,
     t   fvdz1,     fvdz2,     fvdz3,     fvdz4,
     u   fvdz5,     fvdz6,     fvdz7,     fvdz8,
     v   fvdx,      fvdy,      fvdz,      fdeltax,
     w   vis,       mxt,       qmv,       bufmat,
     x   iparg1,    ixs,       gbuf%TAG22,nc1,
     y   nc2,       nc3,       nc4,       nc5,
     z   nc6,       nc7,       nc8,       nale,
     1   nel,       nft,       mtn,       ipm)
       ENDIF     
      ENDIF !  IALEFVM_MAT <= 1
C
      IF (ioutprt>0)THEN
C Energy, momentum, mass balance by part in case of print-out
           eint => elbuf_tab(ng)%GBUF%EINT(1:nel)
           IF (jcvt == 0) THEN
            ifvm22      = 0
            ialefvm_flg = ipm(251,imat)
            IF(ialefvm_flg>=2)ifvm22 = 1
            CALL sbilan(partsav,eint,gbuf%RHO,gbuf%RK,gbuf%VOL,
     .                  vx1, vx2, vx3, vx4, vx5, vx6, vx7, vx8,
     .                  vy1, vy2, vy3, vy4, vy5, vy6, vy7, vy8,
     .                  vz1, vz2, vz3, vz4, vz5, vz6, vz7, vz8,
     .                  voln,iparts,gresav,grth,igrth,gbuf%OFF,
     .                  iexpan,gbuf%EINTTH,gbuf%FILL,gbuf%MOM,nel,ifvm22,
     .                  x1, x2, x3, x4, x5, x6, x7, x8,
     .                  y1, y2, y3, y4, y5, y6, y7, y8,
     .                  z1, z2, z3, z4, z5, z6, z7, z8,
     .                  itask,iparg1,sensors,gbuf%G_WPLA,gbuf%WPLA)
           ELSE
            CALL srbilan(partsav,eint    ,gbuf%RHO,gbuf%RK,gbuf%VOL ,
     .                   vgxa   ,vgya    ,vgza    ,vga2    ,voln     ,
     .                   iparts ,gresav  ,grth    ,igrth   ,gbuf%OFF ,
     .                   iexpan ,gbuf%EINTTH,gbuf%FILL,xgxa, xgya, xgza,
     .                   xgxa2,xgya2,xgza2,xgxya,xgyza,xgzxa,itask,iparg1,sensors,
     .                   nel,gbuf%G_WPLA,gbuf%WPLA)
           ENDIF
      ENDIF
C
C Compute nodal internal forces
      IF (isrot == 0)THEN
        CALL sfint3(
     1   gbuf%SIG,px1,     px2,     px3,
     2   px4,     py1,     py2,     py3,
     3   py4,     pz1,     pz2,     pz3,
     4   pz4,     px5,     px6,     px7,
     5   px8,     py5,     py6,     py7,
     6   py8,     pz5,     pz6,     pz7,
     7   pz8,     f11,     f21,     f31,
     8   f12,     f22,     f32,     f13,
     9   f23,     f33,     f14,     f24,
     a   f34,     f15,     f25,     f35,
     b   f16,     f26,     f36,     f17,
     c   f27,     f37,     f18,     f28,
     d   f38,     voln,    qvis,    n1x,
     e   n2x,     n3x,     n4x,     n5x,
     f   n6x,     n1y,     n2y,     n3y,
     g   n4y,     n5y,     n6y,     n1z,
     h   n2z,     n3z,     n4z,     n5z,
     i   n6z,     dfe,     p,       ixs,
     j   nel,     nft,     jale,    jeul,
     k   svis)
      ELSE
       iisrot = 1
        CALL sfint3b(
     1   gbuf%SIG, lbuf%SIGD,x,        nc1,
     2   nc2,      nc3,      nc4,      nc5,
     3   nc6,      nc7,      nc8,      px1,
     4   px2,      px3,      px4,      py1,
     5   py2,      py3,      py4,      pz1,
     6   pz2,      pz3,      pz4,      px5,
     7   px6,      px7,      px8,      py5,
     8   py6,      py7,      py8,      pz5,
     9   pz6,      pz7,      pz8,      f11,
     a   f21,      f31,      f12,      f22,
     b   f32,      f13,      f23,      f33,
     c   f14,      f24,      f34,      f15,
     d   f25,      f35,      f16,      f26,
     e   f36,      f17,      f27,      f37,
     f   f18,      f28,      f38,      voln,
     g   qvis,     dr,       gbuf%RHO, cxx,
     h   mx1,      my1,      mz1,      mx2,
     i   my2,      mz2,      mx3,      my3,
     j   mz3,      mx4,      my4,      mz4,
     k   mx5,      my5,      mz5,      mx6,
     l   my6,      mz6,      mx7,      my7,
     m   mz7,      mx8,      my8,      mz8,
     n   ixs,      nel,      ismstr,   svis)
      ENDIF
C
      IF (jale  == 3) THEN
C Flow internal forces
          nbf1 = nel*4
          nbf2 = nel + nbf1
          nbf3 = nel + nbf2
          DO i=1,nel
             sigf(i,1) = mbuf%VAR(nbf1 + i)                              
             sigf(i,2) = mbuf%VAR(nbf2 + i)                              
             sigf(i,3) = mbuf%VAR(nbf3 + i)                              
             sigf(i,4) = zero                                            
             sigf(i,5) = zero                                            
             sigf(i,6) = zero                                            
          ENDDO 
          CALL ffint3(
     1   sigf,    dfe,     fpx1,    fpx2,
     2   fpx3,    fpx4,    fpy1,    fpy2,
     3   fpy3,    fpy4,    fpz1,    fpz2,
     4   fpz3,    fpz4,    ff11,    ff21,
     5   ff31,    ff12,    ff22,    ff32,
     6   ff13,    ff23,    ff33,    ff14,
     7   ff24,    ff34,    ff15,    ff25,
     8   ff35,    ff16,    ff26,    ff36,
     9   ff17,    ff27,    ff37,    ff18,
     a   ff28,    ff38,    fvoln,   fqvis,
     b   fn1x,    fn2x,    fn3x,    fn4x,
     c   fn5x,    fn6x,    fn1y,    fn2y,
     d   fn3y,    fn4y,    fn5y,    fn6y,
     e   fn1z,    fn2z,    fn3z,    fn4z,
     f   fn5z,    fn6z,    llt)
      ENDIF
C
      ialefvm_mat = ipm(251,ixs(1,nf1))
      IF (ialefvm_mat > 1) THEN
C ALEFVM forces
C External forces work
          CALL alefvm_tfext(
     1   ixs,  nv46 ,   gbuf%SIG,  ale_connect,
     2   n1x,  n2x,     n3x,       n4x,         n5x,   n6x,
     4   n1y,  n2y,     n3y,       n4y,         n5y,   n6y,
     5   n1z,  n2z,     n3z,       n4z,         n5z,   n6z,
     7   x,    ipm,     nale,
     8   w,    p,       nel,        wfext)
C
C Initial forces
          CALL alefvm_freset(gbuf%TAG22)
C
C Gravity forces
          CALL alefvm_gravity(gbuf%VOL, ixs, gbuf%RHO)
          IF(int22/=0)THEN
            CALL alefvm_gravity_int22(gbuf%VOL, ixs, gbuf%RHO,  gbuf%TAG22)
          ENDIF
C
C Internal forces
          CALL alefvm_stress(
     1   ixs,       gbuf%SIG,  qvis,
     2   n1x,       n2x,       n3x,       n4x,
     3   n5x,       n6x,       n1y,       n2y,
     4   n3y,       n4y,       n5y,       n6y,
     5   n1z,       n2z,       n3z,       n4z,
     6   n5z,       n6z,       x,         ipm,
     7   gbuf%RHO,  gbuf%VOL,
     8   nel,       gbuf%MOM,  lbuf%SSP,  svis)
          IF(int22/=0)THEN
            CALL alefvm_stress_int22(
     1   ixs,       gbuf%SIG,  qvis,
     2   n1x,       n2x,       n3x,       n4x,      n5x,       n6x,
     3   n1y,       n2y,       n3y,       n4y,      n5y,       n6y,
     4   n1z,       n2z,       n3z,       n4z,      n5z,       n6z,
     5   ipm,       gbuf%RHO,  gbuf%VOL,  gbuf%TAG22,
     6   nel,       gbuf%MOM,  lbuf%SSP,  svis)
          ENDIF    
      ENDIF ! IALEFVM_MAT > 1 
C
c===========================================================================
c     Computation of the internal non-local forces
c===========================================================================
      IF (inloc > 0) THEN    
C Virtual internal forces of regularized non local dof
        CALL sfint_reg(
     1    nloc_dmg ,var_reg  ,nel      ,off      ,
     2    voln     ,nc1      ,nc2      ,nc3      ,
     3    nc4      ,nc5      ,nc6      ,nc7      ,
     4    nc8      ,px1      ,px2      ,px3      ,
     5    px4      ,py1      ,py2      ,py3      ,
     6    py4      ,pz1      ,pz2      ,pz3      ,
     7    pz4      ,mxt(lft) ,itask    ,dt2t     ,
     8    gbuf%VOL ,nft      )
      ENDIF 
c===========================================================================  
C
      IF (jthe < 0) THEN
C Finite element heat transfert forces
        IF (mat_elem%MAT_PARAM(imat)%HEAT_FLAG == 1) THEN
          CALL stherm(mat_elem%MAT_PARAM(imat),
     1         voln,    nc1,     nc2,     nc3,     nc4,     
     2         nc5,     nc6,     nc7,     nc8,     px1,
     3         px2,     px3,     px4,     py1,
     4         py2,     py3,     py4,     pz1,
     5         pz2,     pz3,     pz4,     dt1,
     6         temp,    tempel,  fheat,   them,
     7         gbuf%OFF,lbuf%OFF,partsav, iparts,
     8         gbuf%VOL,nel     ,numnod ,glob_therm%THEACCFACT)
        ELSE
          CALL stherm(mat_elem%MAT_PARAM(imat),
     1         voln,    nc1,     nc2,     nc3,     nc4,
     2         nc5,     nc6,     nc7,     nc8,     px1,
     3         px2,     px3,     px4,     py1,
     4         py2,     py3,     py4,     pz1,
     5         pz2,     pz3,     pz4,     dt1,
     6         temp,    tempel,  die,     them,
     7         gbuf%OFF,lbuf%OFF,partsav, iparts,
     8         gbuf%VOL,nel     ,numnod ,glob_therm%THEACCFACT)
        END IF
      ENDIF
C
      IF (jcvt /= 0) THEN
C Transform forces from convected frame to global frame
       CALL srrota3(
     1   r11,     r21,     r31,     r12,
     2   r22,     r32,     r13,     r23,
     3   r33,     f11,     f12,     f13,
     4   f14,     f15,     f16,     f17,
     5   f18,     f21,     f22,     f23,
     6   f24,     f25,     f26,     f27,
     7   f28,     f31,     f32,     f33,
     8   f34,     f35,     f36,     f37,
     9   f38,     nel)
        IF(jale == 3) THEN
         CALL srrota3(
     1   r11,     r21,     r31,     r12,
     2   r22,     r32,     r13,     r23,
     3   r33,     ff11,    ff12,    ff13,
     4   ff14,    ff15,    ff16,    ff17,
     5   ff18,    ff21,    ff22,    ff23,
     6   ff24,    ff25,    ff26,    ff27,
     7   ff28,    ff31,    ff32,    ff33,
     8   ff34,    ff35,    ff36,    ff37,
     9   ff38,    nel)
        ENDIF 
      ENDIF
C----------------------------
C     distortion control  
C----------------------------
      IF (isctl > 0) THEN
         CALL sdistor_ini(                                      
     1                nel      ,sti_c    ,npropm     ,nummat  ,       
     2                ismstr   ,mxt      ,istab      ,pm      ,       
     3                gbuf%SIG ,gbuf%RHO ,cxx        ,off     ,       
     4                gbuf%OFF ,ll       ,voln       ,fld     ,       
     5                cns2     ,fqmax    )
!      x,v sould be in global sys  
         IF (ismstr == 1.OR.ismstr >= 11.OR.jcvt > 0) 
     .     CALL s8get_x3(      
     .     nc1,      nc2,      nc3,      nc4,
     .     nc5,      nc6,      nc7,      nc8,
     .      x1,       x2,       x3,       x4,      
     .      x5,       x6,       x7,       x8,      
     .      y1,       y2,       y3,       y4,      
     .      y5,       y6,       y7,       y8,      
     .      z1,       z2,       z3,       z4,      
     .      z5,       z6,       z7,       z8,      
     .       x,      xdp,   numnod,     nel )
         IF (jcvt > 0) THEN
           DO i=1,nel
            vx1(i)=v(1,nc1(i))
            vy1(i)=v(2,nc1(i))
            vz1(i)=v(3,nc1(i))
            vx2(i)=v(1,nc2(i))
            vy2(i)=v(2,nc2(i))
            vz2(i)=v(3,nc2(i))
            vx3(i)=v(1,nc3(i))
            vy3(i)=v(2,nc3(i))
            vz3(i)=v(3,nc3(i))
            vx4(i)=v(1,nc4(i))
            vy4(i)=v(2,nc4(i))
            vz4(i)=v(3,nc4(i))
            vx5(i)=v(1,nc5(i))
            vy5(i)=v(2,nc5(i))
            vz5(i)=v(3,nc5(i))
            vx6(i)=v(1,nc6(i))
            vy6(i)=v(2,nc6(i))
            vz6(i)=v(3,nc6(i))
            vx7(i)=v(1,nc7(i))
            vy7(i)=v(2,nc7(i))
            vz7(i)=v(3,nc7(i))
            vx8(i)=v(1,nc8(i))
            vy8(i)=v(2,nc8(i))
            vz8(i)=v(3,nc8(i))
           ENDDO
         END IF !(JCVT > 0) THEN
         CALL s8for_distor(
     .   x1,      x2,      x3,      x4,      
     .   x5,      x6,      x7,      x8,      
     .   y1,      y2,      y3,      y4,      
     .   y5,      y6,      y7,      y8,      
     .   z1,      z2,      z3,      z4,      
     .   z5,      z6,      z7,      z8,      
     .   vx1,     vx2,     vx3,     vx4,
     .   vx5,     vx6,     vx7,     vx8,
     .   vy1,     vy2,     vy3,     vy4,
     .   vy5,     vy6,     vy7,     vy8,
     .   vz1,     vz2,     vz3,     vz4,
     .   vz5,     vz6,     vz7,     vz8,
     .   f11,     f12,     f13,     f14,
     .   f15,     f16,     f17,     f18,
     .   f21,     f22,     f23,     f24,
     .   f25,     f26,     f27,     f28,
     .   f31,     f32,     f33,     f34,
     .   f35,     f36,     f37,     f38,
     .   sti,   sti_c,   fld  ,    cns2,
     .   ll ,   istab,   fqmax,     nel,
     .   gbuf%EINT_DISTOR, dt1)
      ENDIF
C
      IF(nfilsol/=0) CALL sfillopt(
     1   gbuf%FILL,sti,      f11,      f21,
     2   f31,      f12,      f22,      f32,
     3   f13,      f23,      f33,      f14,
     4   f24,      f34,      f15,      f25,
     5   f35,      f16,      f26,      f36,
     6   f17,      f27,      f37,      f18,
     7   f28,      f38,      nel)
C
C Assemble internal forces
      IF(iparit == 0)THEN
          CALL scumu3(
     1   gbuf%OFF,a,       nc1,     nc2,
     2   nc3,     nc4,     nc5,     nc6,
     3   nc7,     nc8,     stifn,   sti,
     4   f11,     f21,     f31,     f12,
     5   f22,     f32,     f13,     f23,
     6   f33,     f14,     f24,     f34,
     7   f15,     f25,     f35,     f16,
     8   f26,     f36,     f17,     f27,
     9   f37,     f18,     f28,     f38,
     a   nvc,     ar,      fr_wave, fr_wav,
     b   mx1,     my1,     mz1,     mx2,
     c   my2,     mz2,     mx3,     my3,
     d   mz3,     mx4,     my4,     mz4,
     e   mx5,     my5,     mz5,     mx6,
     f   my6,     mz6,     mx7,     my7,
     g   mz7,     mx8,     my8,     mz8,
     h   them,    fthe,    condn,   conde,
     i   nel,     jthe,    isrot,   ipartsph,glob_therm%NODADT_THERM)
      ELSE
          CALL scumu3p(
     1   gbuf%OFF,sti,     fsky,    fsky,
     2   iads,    f11,     f21,     f31,
     3   f12,     f22,     f32,     f13,
     4   f23,     f33,     f14,     f24,
     5   f34,     f15,     f25,     f35,
     6   f16,     f26,     f36,     f17,
     7   f27,     f37,     f18,     f28,
     8   f38,     nc1,     nc2,     nc3,
     9   nc4,     nc5,     nc6,     nc7,
     a   nc8,     ar,      fr_wave, fr_wav,
     b   mx1,     my1,     mz1,     mx2,
     c   my2,     mz2,     mx3,     my3,
     d   mz3,     mx4,     my4,     mz4,
     e   mx5,     my5,     mz5,     mx6,
     f   my6,     mz6,     mx7,     my7,
     g   mz7,     mx8,     my8,     mz8,
     h   them,    fthesky, condnsky,conde,
     i   nel,     nft,     jthe,    isrot,
     j   ipartsph,glob_therm%NODADT_THERM)
      ENDIF
C
      IF(jale  == 3) THEN   
       IF(iparit == 0 ) THEN
          CALL fcumu3(
     1   gbuf%OFF,af,      nc1,     nc2,
     2   nc3,     nc4,     nc5,     nc6,
     3   nc7,     nc8,     stifn,   sti,
     4   ff11,    ff21,    ff31,    ff12,
     5   ff22,    ff32,    ff13,    ff23,
     6   ff33,    ff14,    ff24,    ff34,
     7   ff15,    ff25,    ff35,    ff16,
     8   ff26,    ff36,    ff17,    ff27,
     9   ff37,    ff18,    ff28,    ff38,
     a   nvc,     nel)
       
       ELSE  
          CALL fcumu3p(
     1   gbuf%OFF,sti,     ffsky,   iads,
     2   ff11,    ff21,    ff31,    ff12,
     3   ff22,    ff32,    ff13,    ff23,
     4   ff33,    ff14,    ff24,    ff34,
     5   ff15,    ff25,    ff35,    ff16,
     6   ff26,    ff36,    ff17,    ff27,
     7   ff37,    ff18,    ff28,    ff38,
     8   nc1,     nc2,     nc3,     nc4,
     9   nc5,     nc6,     nc7,     nc8,
     a   nel,     nft)
       ENDIF   
      ENDIF
C
      IF (ALLOCATED(var_reg)) DEALLOCATE(var_reg)
      RETURN
      END