sms_mass_scale_2.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/.
 
!||====================================================================
!||    sms_mass_scale_2   ../engine/source/ams/sms_mass_scale_2.F
!||--- called by ------------------------------------------------------
!||    resol              ../engine/source/engine/resol.F
!||--- calls      -----------------------------------------------------
!||    ancmsg             ../engine/source/output/message/message.F
!||    arret              ../engine/source/system/arret.F
!||    my_barrier         ../engine/source/system/machine.F
!||    ngr2usr            ../engine/source/input/freform.F
!||    sms_admesh_0       ../engine/source/ams/sms_admesh.F
!||    sms_bcs            ../engine/source/ams/sms_bcs.F
!||    sms_bcscyc         ../engine/source/ams/sms_bcscyc.F
!||    sms_cjoint_0       ../engine/source/ams/sms_cjoint.F
!||    sms_cjoint_2       ../engine/source/ams/sms_cjoint.F
!||    sms_fixvel         ../engine/source/ams/sms_fixvel.F
!||    sms_gravit         ../engine/source/ams/sms_gravit.F
!||    sms_mav_lt         ../engine/source/ams/sms_pcg.F
!||    sms_pcg            ../engine/source/ams/sms_pcg.F
!||    sms_rbe3t1         ../engine/source/ams/sms_rbe3.F
!||    sms_rbe_cnds       ../engine/source/ams/sms_rbe2.F
!||    sms_thbcs          ../engine/source/ams/sms_thbcs.F
!||    spmd_exch_a_rb6    ../engine/source/mpi/kinematic_conditions/spmd_exch_a_rb6.F
!||    spmd_list_sms      ../engine/source/mpi/ams/spmd_sms.F
!||    spmd_mij_sms       ../engine/source/mpi/ams/spmd_sms.F
!||    startime           ../engine/source/system/timer_mod.F90
!||    stoptime           ../engine/source/system/timer_mod.F90
!||--- uses       -----------------------------------------------------
!||    ams_work_mod       ../engine/source/modules/ams_work_mod.F90
!||    element_mod        ../common_source/modules/elements/element_mod.F90
!||    groupdef_mod       ../common_source/modules/groupdef_mod.F
!||    intstamp_mod       ../engine/share/modules/intstamp_mod.F
!||    message_mod        ../engine/share/message_module/message_mod.F
!||    my_alloc_mod       ../common_source/tools/memory/my_alloc.F90
!||    python_funct_mod   ../common_source/modules/python_mod.f90
!||    sensor_mod         ../common_source/modules/sensor_mod.F90
!||    timer_mod          ../engine/source/system/timer_mod.F90
!||====================================================================
      SUBROUTINE sms_mass_scale_2(TIMERS,PYTHON,
     1        ITASK    ,NODFT    ,NODLT    ,NODII_SMS ,INDX2_SMS ,
     2        NODXI_SMS,MS       ,MS0      ,A         ,ICODT     ,
     3        ICODR    ,ISKEW    ,SKEW     ,JAD_SMS   ,JDI_SMS   ,
     4        LT_SMS   ,X_SMS    ,P_SMS    ,Z_SMS     ,Y_SMS     ,
     5        PREC_SMS ,INDX1_SMS,DIAG_SMS ,IAD_ELEM  ,FR_ELEM   ,
     6        WEIGHT   ,NPBY      ,LPBY      ,
     7        TAGSLV_RBY_SMS,LAD_SMS ,KAD_SMS ,JRB_SMS,IBFV      ,
     8        VEL      ,NPC      ,TF       ,V         ,X         ,
     9        D        ,SENSOR_TAB,NSENSOR ,IFRAME    ,XFRAME    ,
     A        JADI_SMS ,JDII_SMS ,LTI_SMS  ,FR_SMS    ,FR_RMS    ,
     B        ISKYI_SMS,MSKYI_SMS,RES_SMS  ,IGRV      ,AGRV      ,
     C        LGRAV    ,ILINK    ,RLINK    ,FR_RL     ,FRL6      ,
     D        NNLINK   ,LNLINK   ,FR_LL    ,FNL6      ,TAG_LNK_SMS,
     E        ITAB     ,FSAV     ,LJOINT   ,IADCJ     ,FR_CJ      ,
     F        AM       ,VR       ,IN       ,FRL       ,FNL        ,
     G        NPRW     ,LPRW     ,RWBUF    ,RWSAV     ,FOPT       ,
     H        FR_WALL  ,NRWL_SMS ,INTSTAMP ,KINET     ,IXC        ,
     I        IXTG     ,SH4TREE  ,SH3TREE  ,CPTREAC   ,NODREAC    ,
     J        FTHREAC  ,FRWL6    ,DIM      ,TAGSLV_RBY,DAMPR      ,
     K        DAMP     ,IGRNOD   ,DR       ,RBY       ,
     L        TAGMSR_RBY_SMS,JSM_SMS,IRBE2  ,LRBE2      ,
     N        IAD_RBE2 ,FR_RBE2M ,NMRBE2   ,R2SIZE    ,IRBE3      ,
     O        LRBE3      ,FRBE3    ,IAD_RBE3M ,FR_RBE3M ,FR_RBE3MP,
     P        RRBE3      ,RRBE3_PON,PREC_SMS3 ,DIAG_SMS3,IAD_RBY  ,
     Q        FR_RBY6    ,RBY6     ,R3SIZE ,BETATE      ,IBCSCYC  ,
     R        LBCSCYC    ,MSKYI_FI_SMS, LIST_SMS,LIST_RMS,CJWORK  ,
     S        FREA       ,IRWL_WORK,VFI,sz_mw6,MW6,WFEXT,ams_work)
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE timer_mod
      USE intstamp_mod
      USE message_mod
      USE groupdef_mod
      USE sensor_mod
      USE ams_work_mod
      USE my_alloc_mod
      use python_funct_mod, only : python_
      use element_mod , only : nixc,nixtg
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "com04_c.inc"
#include      "com06_c.inc"
#include      "com08_c.inc"
#include      "param_c.inc"
#include      "parit_c.inc"
#include      "remesh_c.inc"
#include      "scr03_c.inc"
#include      "sms_c.inc"
#include      "tabsiz_c.inc"
#include      "task_c.inc"
#include      "timeri_c.inc"
#include      "units_c.inc"
#include      "warn_c.inc"
#include      "stati_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      TYPE(timer_), INTENT(INOUT) :: TIMERS
      TYPE(python_), INTENT(INOUT) :: PYTHON
      INTEGER  ITASK, NODFT,NSENSOR,NODLT, NODII_SMS(*), INDX2_SMS(*),
     .         NODXI_SMS(*), ICODT(*), ICODR(*),
     .         ISKEW(*), JAD_SMS(*), JDI_SMS(*), INDX1_SMS(*),
     .         IAD_ELEM(2,NSPMD+1) ,FR_ELEM(*), WEIGHT(*),
     .         NPBY(NNPBY,*), LPBY(*), TAGSLV_RBY_SMS(*), TAGSLV_RBY(*),
     .         LAD_SMS(*), KAD_SMS(*), JRB_SMS(*),
     .         NPC(*), IBFV(NIFV,*),IFRAME(LISKN,*),
     .         JADI_SMS(*), JDII_SMS(*),
     .         FR_RMS(NSPMD+1), FR_SMS(NSPMD+1), ISKYI_SMS(LSKYI_SMS,*),
     .         IGRV(*),CPTREAC,NODREAC(*),
     .         ILINK(*),RLINK(*), FR_RL(NSPMD+2,*), NNLINK(10,*),
     .         LNLINK(*), FR_LL(NSPMD+2,*), TAG_LNK_SMS(*), ITAB(*),
     .         LJOINT(*), FR_CJ(*), IADCJ(*),
     .         NPRW(*), LPRW(*), FR_WALL(*), NRWL_SMS(*),
     .         KK, MAIN, KINET(*),
     .         IXC(NIXC,*), IXTG(NIXTG,*),
     .         SH4TREE(KSH4TREE,*), SH3TREE(KSH3TREE,*), DIM,
     .         TAGMSR_RBY_SMS(*), JSM_SMS(*),
     .         IRBE2(NRBE2L,*), LRBE2(*), IAD_RBE2(*),
     .         FR_RBE2M(*), NMRBE2, R2SIZE, IRBE3(NRBE3L,*), LRBE3(*),
     .         IAD_RBE3M(*) ,FR_RBE3M(*) ,FR_RBE3MP(*),
     .         FR_RBY6(*),IAD_RBY(*),R3SIZE,IBCSCYC(*),LBCSCYC(*)
      my_real
     .    MS(*), MS0(*), A(3,*), DIAG_SMS(*),
     .    SKEW(LSKEW,*), LT_SMS(*),
     .    x_sms(3,*), p_sms(3,*), y_sms(3,*), z_sms(3,*), prec_sms(*),
     .    v(3,*), x(3,*), d(3,*), tf(*), vel(lfxvelr,*),
     .    xframe(nxframe,*), lti_sms(*), mskyi_sms(*),
     .    res_sms(3,*), agrv(*),lgrav(*),
     .    fsav(nthvki,*), am(3,*), vr(3,*), in(*), frl(*), fnl(*),
     .    rwbuf(*), rwsav(*), fopt(*),fthreac(6,*),
     .    dampr(nrdamp,*), damp(dim,*), dr(3,*), rby(nrby,*), 
     .    frbe3(*), rrbe3(*), prec_sms3(3,*), diag_sms3(3,*),betate
      DOUBLE PRECISION FRL6(*), FNL6(*), FRWL6(*), RRBE3_PON(*)
      DOUBLE PRECISION RBY6(8,6,NRBYKIN)
      my_real,dimension(fr_rms(nspmd+1)),intent(inout) :: MSKYI_FI_SMS
      integer,dimension(fr_sms(nspmd+1)),intent(inout) :: LIST_SMS
      integer,dimension(fr_rms(nspmd+1)),intent(inout) :: LIST_RMS
      my_real, DIMENSION(18,NJOINT),intent(inout):: CJWORK
      my_real, DIMENSION(3,NUMNOD),intent(inout):: FREA
      integer, dimension(SLPRW),intent(inout):: IRWL_WORK
      my_real, DIMENSION(3,FR_RMS(NSPMD+1)+FR_SMS(NSPMD+1) ), intent(inout):: VFI
      integer, intent(in) :: sz_mw6
      DOUBLE PRECISION,dimension(6,sz_mw6),intent(inout) :: MW6
      DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
 
      TYPE(INTSTAMP_DATA) INTSTAMP(*)
      TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
C-----------------------------------------------
      TYPE (GROUP_)  , DIMENSION(NGRNOD) :: IGRNOD
      TYPE (AMS_WORK_) , INTENT(INOUT) :: AMS_WORK
C-----------------------------------------------
C   L o c a l  V a r i a b l e s
C-----------------------------------------------
      INTEGER I, N, ISP, IT, IX, IERROR
      INTEGER ICOUNT, J, K, L, NSN, IMOV, ITYP, ILAGM, IFLAG,
     .        N2, N3, N4, N5, N6, N7, ND, IGR, ISK,
     .        M, IAD, MSR, KAD, KI, KJ, JI, NSR,
     .        LOC_PROC, P, NN, LENR, SIZE, NRBDIM
      INTEGER NODFT1_SMS, NODLT1_SMS
      INTEGER NODFT2_SMS, NODLT2_SMS,NGR2USR
      my_real
     .   vx,vy,vz, mvx, mvy, mvz,
     .   vxj, vyj, vzj, mas,wfextt, errtet, dw, dt15, dt25, rbid,
     .   omega, betasdt, dampt, factb, d_tstart, d_tstop, da, adt,
     .   p1, p2, p3, uomega, domega
C-----
      INTEGER, DIMENSION(:), ALLOCATABLE :: IMV
      my_real
     .       , DIMENSION(:), ALLOCATABLE :: mv
      my_real, DIMENSION(:,:), ALLOCATABLE :: mvskw
      my_real, DIMENSION(:,:), ALLOCATABLE :: vskw
      my_real, DIMENSION(:,:), ALLOCATABLE :: rskw
      my_real, DIMENSION(:,:), ALLOCATABLE :: dampskw
      double precision
     .       , DIMENSION(:,:), ALLOCATABLE :: mv6
      EXTERNAL ngr2usr
C-----------------------------------------------
      CALL my_alloc(mvskw,3,numnod)
      CALL my_alloc(vskw,3,numnod)
      CALL my_alloc(rskw,3,numnod)
      CALL my_alloc(dampskw,3,numnod)
C-----------------------------------------------
      frea(1:3,nodft:nodlt)=zero
C
      IF(iparit/=0)THEN
        IF(debug(9)==0)THEN
          ALLOCATE(imv(2*nisky_sms+fr_rms(nspmd+1)),
     .           mv(3*(2*nisky_sms+fr_rms(nspmd+1))),
     .           mv6(6,3*(2*nisky_sms+fr_rms(nspmd+1))),stat=ierror)
        ELSE
          ALLOCATE(imv(numnod+nnz_sms+2*nisky_sms+fr_rms(nspmd+1)),
     .           mv(3*(numnod+nnz_sms+2*nisky_sms+fr_rms(nspmd+1))),
     .           mv6(6,3*(numnod+nnz_sms+2*nisky_sms+fr_rms(nspmd+1))),
     .           stat=ierror)
        END IF
        IF(ierror/=0) THEN
          WRITE(istdo,*)
          CALL ancmsg(msgid=19,anmode=aninfo,
     .                c1='(/DT/.../AMS)')
          CALL arret(2)
        ENDIF
      END IF
C
      IF(nspmd > 1)THEN
        IF(itask==0)THEN
          CALL spmd_list_sms(iskyi_sms,fr_sms,fr_rms,list_sms,list_rms,
     .                       npby     ,tagslv_rby_sms)
        END IF
C
        CALL my_barrier
C
      END IF
C
C----
C
      IF(nspmd > 1)THEN
C
        CALL my_barrier()
C
        IF(itask==0) THEN   ! communication on the first thread
          CALL spmd_mij_sms(
     1           iskyi_sms,fr_sms,fr_rms,list_rms,mskyi_sms,
     2           mskyi_fi_sms)
        END IF
      END IF
C----
C
      CALL sms_gravit(igrv  ,agrv   ,npc   ,tf   ,a     ,
     2                v     ,x      ,skew  ,ms   ,sensor_tab,
     3                weight,lgrav ,itask,tagslv_rby_sms,nsensor,wfext, python)
C
      CALL my_barrier
C
      nodft1_sms=1+itask*nindx1_sms/nthread
      nodlt1_sms=(itask+1)*nindx1_sms/nthread
C
      nodft2_sms=1+itask*nindx2_sms/nthread
      nodlt2_sms=(itask+1)*nindx2_sms/nthread
C----
C
      DO n=nodft,nodlt
 
        a(1,n)=a(1,n)+res_sms(1,n)
        a(2,n)=a(2,n)+res_sms(2,n)
        a(3,n)=a(3,n)+res_sms(3,n)
 
        res_sms(1,n)=zero
        res_sms(2,n)=zero
        res_sms(3,n)=zero
 
      END DO
C
      CALL my_barrier
C
C--------------------------------------------
C     RAILEIGH DAMPING
C--------------------------------------------
      IF(ndamp/=0.OR.istat==1.OR.istat==3)THEN
C
        DO n=nodft,nodlt
         IF(nodxi_sms(n)==0)THEN
           z_sms(1,n)=ms(n)*v(1,n)
           z_sms(2,n)=ms(n)*v(2,n)
           z_sms(3,n)=ms(n)*v(3,n)
         ELSE
           x_sms(1,n)=v(1,n)
           x_sms(2,n)=v(2,n)
           x_sms(3,n)=v(3,n)
         END IF
        ENDDO
C-----------------------------------
        IF(nrbody/=0)THEN
C
          CALL my_barrier()
C 
          DO n=nodft1_sms,nodlt1_sms
           i=indx1_sms(n)
           m=tagslv_rby_sms(i)
           IF(m /= 0)THEN
             msr=npby(1,m)
             x_sms(1,i)=x_sms(1,msr)
             x_sms(2,i)=x_sms(2,msr)
             x_sms(3,i)=x_sms(3,msr)
           END IF
          END DO
C 
        END IF
C
        CALL my_barrier
C
C       Z_SMS temporarily uses for [m] V
        CALL sms_mav_lt(timers,
     1           nodft   ,nodlt  ,numnod ,jad_sms  ,jdi_sms  ,
     2           itask   ,diag_sms,lt_sms,x_sms ,z_sms ,
     3           nodft1_sms,nodlt1_sms,indx1_sms,nodxi_sms,iad_elem ,
     4           fr_elem   ,weight    ,jadi_sms ,jdii_sms ,lti_sms  ,
     5           iskyi_sms ,mskyi_sms ,fr_sms   ,fr_rms   ,list_sms ,
     6           list_rms  ,mskyi_fi_sms ,vfi   ,imv      ,mv      ,
     7           mv6       ,mw6       ,nodft2_sms,nodlt2_sms,indx2_sms,
     8           nodii_sms )
C
        CALL my_barrier
C
C-----------------------------------
C     remontee Yi => Ym
C-----------------------------------
        IF(nrbody/=0)THEN
C
!$OMP DO SCHEDULE(DYNAMIC,1)
          DO m =1,nrbody
            DO k = 1, 6
              rby6(1,k,m) = zero
              rby6(2,k,m) = zero
              rby6(3,k,m) = zero
            END DO
C
            msr=npby(1,m)
            IF(msr < 0) cycle
C
            IF(tagmsr_rby_sms(msr) /= 0) THEN
              rby6(1,1,m)=z_sms(1,msr)*weight(msr)
              rby6(2,1,m)=z_sms(2,msr)*weight(msr)
              rby6(3,1,m)=z_sms(3,msr)*weight(msr)
            END IF
C
          END DO
!$OMP END DO
 
!$OMP SINGLE
          DO n=1,nindx1_sms
           i=indx1_sms(n)
           m=tagslv_rby_sms(i)
           IF(m /= 0)THEN
             IF(weight(i) /= 0)THEN
               rby6(1,1,m)=rby6(1,1,m)+z_sms(1,i)
               rby6(2,1,m)=rby6(2,1,m)+z_sms(2,i)
               rby6(3,1,m)=rby6(3,1,m)+z_sms(3,i)
             END IF
             z_sms(1,i)=zero
             z_sms(2,i)=zero
             z_sms(3,i)=zero
           END IF
          END DO
!$OMP END SINGLE
 
         IF (nspmd > 1) THEN
!$OMP SINGLE
           nrbdim=3
           CALL spmd_exch_a_rb6(
     1       nrbdim,iad_rby,fr_rby6,iad_rby(nspmd+1),rby6)
!$OMP END SINGLE
         END IF
 
!$OMP DO SCHEDULE(DYNAMIC,1)
          DO m =1,nrbody
            msr=npby(1,m)
            IF(msr < 0) cycle
            IF(tagmsr_rby_sms(msr) /= 0) THEN
              z_sms(1,msr)=rby6(1,1,m)
              z_sms(2,msr)=rby6(2,1,m)
              z_sms(3,msr)=rby6(3,1,m)
            END IF
          END DO
!$OMP END DO
C
        END IF
C
        CALL my_barrier
C
C-----------------------------------
        IF(itask==0)THEN
         IF (imon>0) CALL startime(timers,5)
         dw = zero
         DO nd=1,ndamp
          igr   = nint(dampr(2,nd))
          isk   = nint(dampr(15,nd))
          factb = dampr(16,nd)
          dampt  = min(dt1,dt2)*factb
          d_tstart = dampr(17,nd)
          d_tstop  = dampr(18,nd)
          IF (tt>=d_tstart .AND. tt<=d_tstop) THEN
            IF(isk<=1)THEN
C----- Damping on rotation dof and only -----
             IF (dampr(19,nd)>0) cycle
             dampa = dampr(3,nd)
             dampb = dampr(4,nd)
             betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
             omega  = one/ (one + half * dampa * dt1)
             DO n=1,igrnod(igr)%NENTITY
               i=igrnod(igr)%ENTITY(n)
               IF(tagslv_rby(i)/=0) cycle
               da=a(1,i)-dampa*z_sms(1,i)-betasdt *(a(1,i) - damp(1,i))
               da = da * omega - a(1,i)
               damp(1,i) = a(1,i)
               a(1,i)    = a(1,i) + da
C              DW =DW+DA*(V(1,I)+HALF*ACC(1,I)*DT1)*DT12*WEIGHT(I)
C              2nd order error
               dw =dw+da*v(1,i)*dt12*weight(i)
             END DO
             dampa = dampr(5,nd)
             dampb = dampr(6,nd)
             betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
             omega  = one/ (one + half * dampa * dt1)
             DO n=1,igrnod(igr)%NENTITY
               i=igrnod(igr)%ENTITY(n)
               IF(tagslv_rby(i)/=0) cycle
               da=a(2,i)-dampa*z_sms(2,i)-betasdt *(a(2,i) - damp(2,i))
               da = da * omega - a(2,i)
               damp(2,i) = a(2,i)
               a(2,i)    = a(2,i) + da
C              2nd order error
               dw =dw+da*v(2,i)*dt12*weight(i)
             END DO
             dampa = dampr(7,nd)
             dampb = dampr(8,nd)
             betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
             omega  = one/ (one + half * dampa * dt1)
             DO n=1,igrnod(igr)%NENTITY
               i=igrnod(igr)%ENTITY(n)
               IF(tagslv_rby(i)/=0) cycle
               da=a(3,i)-dampa*z_sms(3,i)-betasdt *(a(3,i) - damp(3,i))
               da = da * omega - a(3,i)
               damp(3,i) = a(3,i)
               a(3,i)    = a(3,i) + da
C              2nd order error
               dw =dw+da*v(3,i)*dt12*weight(i)
             END DO
            ELSE
#include "vectorize.inc"
             DO n=1,igrnod(igr)%NENTITY
               i=igrnod(igr)%ENTITY(n)
               IF(tagslv_rby(i)/=0) cycle
               mvskw(1,i)= skew(1,isk)*z_sms(1,i)
     .                    +skew(2,isk)*z_sms(2,i)
     .                    +skew(3,isk)*z_sms(3,i)
               mvskw(2,i)= skew(4,isk)*z_sms(1,i)
     .                    +skew(5,isk)*z_sms(2,i)
     .                    +skew(6,isk)*z_sms(3,i)
               mvskw(3,i)= skew(7,isk)*z_sms(1,i)
     .                    +skew(8,isk)*z_sms(2,i)
     .                    +skew(9,isk)*z_sms(3,i)
               vskw(1,i)= skew(1,isk)*v(1,i)
     .                   +skew(2,isk)*v(2,i)
     .                   +skew(3,isk)*v(3,i)
               vskw(2,i)= skew(4,isk)*v(1,i)
     .                   +skew(5,isk)*v(2,i)
     .                   +skew(6,isk)*v(3,i)
               vskw(3,i)= skew(7,isk)*v(1,i)
     .                   +skew(8,isk)*v(2,i)
     .                   +skew(9,isk)*v(3,i)
               rskw(1,i)= skew(1,isk)*a(1,i)
     .                   +skew(2,isk)*a(2,i)
     .                   +skew(3,isk)*a(3,i)
               rskw(2,i)= skew(4,isk)*a(1,i)
     .                   +skew(5,isk)*a(2,i)
     .                   +skew(6,isk)*a(3,i)
               rskw(3,i)= skew(7,isk)*a(1,i)
     .                   +skew(8,isk)*a(2,i)
     .                   +skew(9,isk)*a(3,i)
               dampskw(1,i)= skew(1,isk)*damp(1,i)
     .                      +skew(2,isk)*damp(2,i)
     .                      +skew(3,isk)*damp(3,i)
               dampskw(2,i)= skew(4,isk)*damp(1,i)
     .                      +skew(5,isk)*damp(2,i)
     .                      +skew(6,isk)*damp(3,i)
               dampskw(3,i)= skew(7,isk)*damp(1,i)
     .                      +skew(8,isk)*damp(2,i)
     .                      +skew(9,isk)*damp(3,i)
             END DO
             dampa = dampr(3,nd)
             dampb = dampr(4,nd)
             betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
             omega  = one/ (one + half * dampa * dt1)
#include "vectorize.inc"
             DO n=1,igrnod(igr)%NENTITY
              i=igrnod(igr)%ENTITY(n)
              IF(tagslv_rby(i)/=0) cycle
              da = rskw(1,i) - dampa*mvskw(1,i)
     .                      - betasdt *(rskw(1,i) - dampskw(1,i))
              da = da * omega - rskw(1,i)
              dampskw(1,i) = rskw(1,i)
              rskw(1,i)    = rskw(1,i) + da
C             2nd order error
              dw =dw+da*vskw(1,i)*dt12*weight(i)
             ENDDO
             dampa = dampr(5,nd)
             dampb = dampr(6,nd)
             betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
             omega  = one/ (one + half * dampa * dt1)
#include "vectorize.inc"
             DO n=1,igrnod(igr)%NENTITY
              i=igrnod(igr)%ENTITY(n)
              IF(tagslv_rby(i)/=0) cycle
              da = rskw(2,i) - dampa*mvskw(2,i)
     .                      - betasdt *(rskw(2,i) - dampskw(2,i))
              da = da * omega - rskw(2,i)
              dampskw(2,i) = rskw(2,i)
              rskw(2,i)    = rskw(2,i) + da
C             2nd order error
              dw =dw+da*vskw(2,i)*dt12*weight(i)
             ENDDO
             dampa = dampr(7,nd)
             dampb = dampr(8,nd)
             betasdt= -min(dampb,dampt)*dt1/max(dt1*dt1,em30)
             omega  = one/ (one + half * dampa * dt1)
#include "vectorize.inc"
             DO n=1,igrnod(igr)%NENTITY
              i=igrnod(igr)%ENTITY(n)
              IF(tagslv_rby(i)/=0) cycle
              da = rskw(3,i) - dampa*mvskw(3,i)
     .                      - betasdt *(rskw(3,i) - dampskw(3,i))
              da = da * omega - rskw(3,i)
              dampskw(3,i) = rskw(3,i)
              rskw(3,i)    = rskw(3,i) + da
C             2nd order error
              dw =dw+da*vskw(3,i)*dt12*weight(i)
             ENDDO
#include "vectorize.inc"
             DO n=1,igrnod(igr)%NENTITY
               i=igrnod(igr)%ENTITY(n)
               IF(tagslv_rby(i)/=0) cycle
               a(1,i)= skew(1,isk)*rskw(1,i)
     .                +skew(4,isk)*rskw(2,i)
     .                +skew(7,isk)*rskw(3,i)
               a(2,i)= skew(2,isk)*rskw(1,i)
     .                +skew(5,isk)*rskw(2,i)
     .                +skew(8,isk)*rskw(3,i)
               a(3,i)= skew(3,isk)*rskw(1,i)
     .                +skew(6,isk)*rskw(2,i)
     .                +skew(9,isk)*rskw(3,i)
               damp(1,i)= skew(1,isk)*dampskw(1,i)
     .                   +skew(4,isk)*dampskw(2,i)
     .                   +skew(7,isk)*dampskw(3,i)
               damp(2,i)= skew(2,isk)*dampskw(1,i)
     .                   +skew(5,isk)*dampskw(2,i)
     .                   +skew(8,isk)*dampskw(3,i)
               damp(3,i)= skew(3,isk)*dampskw(1,i)
     .                   +skew(6,isk)*dampskw(2,i)
     .                   +skew(9,isk)*dampskw(3,i)
             END DO
            END IF
          END IF
         END DO
#include "lockon.inc"
         wfext = wfext + dw
#include "lockoff.inc"
         IF (imon>0) CALL stoptime(timers,5)
        END IF
C
        CALL my_barrier
C-----------------------------------
       IF (istat==1.OR.istat==3) THEN
!$OMP SINGLE
         omega  = betate * dt12
         uomega = one - omega
         domega = two*betate
         dw = zero
        IF(istatg==0)THEN
          DO j= 1,3
             DO i=1,numnod
               IF(tagslv_rby(i)/=0) cycle
               da = a(j,i)
               a(j,i)  = uomega*a(j,i) -domega*z_sms(j,i)  
               da = a(j,i) -da
               dw =dw+da*v(j,i)*dt12*weight(i)
             END DO
          END DO
          ELSE
         IF(istatg<0)THEN
          istatg=ngr2usr(-istatg,igrnod,ngrnod)
         ENDIF
         DO j= 1,3
#include "vectorize.inc"
             DO n=1,igrnod(istatg)%NENTITY
               i=igrnod(istatg)%ENTITY(n)
               IF(tagslv_rby(i)/=0) cycle
               da = a(j,i)
               a(j,i)  = uomega*a(j,i) -domega*z_sms(j,i)  
               da = a(j,i) -da
               dw =dw+da*v(j,i)*dt12*weight(i)
             END DO
         END DO
        END IF !(ISTATG==0)THEN
#include "lockon.inc"
         wfext = wfext + dw
#include "lockoff.inc"
!$OMP END SINGLE
C
      CALL my_barrier
       END IF !(ISTAT==3) THEN
C
      END IF
 
C-----------------------------------
C RBE2
C-----------------------------------
      IF (nrbe2>0.OR.r2size>0) THEN      
        IF(itask==0)THEN
          CALL sms_rbe_cnds(
     1     irbe2 ,lrbe2 ,x      ,a      ,am     ,
     1     ms    ,in    ,skew   ,weight ,iad_rbe2,
     2     fr_rbe2m,nmrbe2)
        END IF
C
        CALL my_barrier
C
      END IF
C-----------------------------------
C RBE3
C-----------------------------------
      IF (nrbe3>0)THEN
        IF(itask==0)THEN
          CALL sms_rbe3t1(
     1      irbe3 ,lrbe3  ,x        ,a       ,frbe3    ,
     2      skew  ,weight ,iad_rbe3m,fr_rbe3m,fr_rbe3mp,
     3      rrbe3 ,rrbe3_pon ,r3size)
        END IF
C
        CALL my_barrier
C
      END IF
C-----------------------------------
C     CONDITIONS AUX LIMITES
      CALL sms_thbcs(nodft1_sms,nodlt1_sms,indx1_sms,icodt  ,icodr    ,
     2               iskew  ,skew      ,a         ,am       ,fthreac  ,
     3               nodreac,cptreac)
C
      CALL sms_bcs(nodft1_sms,nodlt1_sms,indx1_sms,icodt  ,iskew ,
     2             skew      ,a         ,nodlt1_sms)
C
      IF(iroddl/=0)
     1  CALL sms_bcs(nodft1_sms,nodlt1_sms,indx1_sms,icodr  ,iskew ,
     2               skew      ,am        ,nodlt1_sms)
C
      IF (nbcscyc>0) CALL sms_bcscyc(ibcscyc,lbcscyc,skew,x,a)
C      
      CALL my_barrier
C
C-----------------------------------
C
C     PREC_SMS used to store the true diagonal (cf rbodies)
      prec_sms(nodft:nodlt)=diag_sms(nodft:nodlt)
C
      CALL my_barrier()
C 
      IF(nrbody/=0)THEN
C
!$OMP DO SCHEDULE(DYNAMIC,1)
        DO m =1,nrbody
          DO k = 1, 6
            rby6(1,k,m) = zero
          END DO
C
          msr=npby(1,m)
          IF(msr < 0) cycle
C
          IF(tagmsr_rby_sms(msr) /= 0) THEN
            rby6(1,1,m)=diag_sms(msr)*weight(msr)
          END IF
C
        END DO
!$OMP  END DO
 
!$OMP SINGLE
       DO n=1,nindx1_sms
        i=indx1_sms(n)
        m=tagslv_rby_sms(i)
        IF(m /= 0)THEN
          IF(weight(i) /= 0)THEN
            rby6(1,1,m)=rby6(1,1,m)+diag_sms(i)
          END IF
        END IF
       END DO
!$OMP END SINGLE
 
       IF (nspmd > 1) THEN
!$OMP SINGLE
         nrbdim=1
         CALL spmd_exch_a_rb6(
     1     nrbdim,iad_rby,fr_rby6,iad_rby(nspmd+1),rby6)
!$OMP END SINGLE
       END IF
 
!$OMP DO SCHEDULE(DYNAMIC,1)
        DO m =1,nrbody
C
          msr=npby(1,m)
C
          IF(msr < 0) cycle
C
          IF(tagmsr_rby_sms(msr) /= 0) THEN
            prec_sms(msr)=rby6(1,1,m)
          END IF
C
        END DO
!$OMP  END DO
C 
        DO n=nodft1_sms,nodlt1_sms
         i=indx1_sms(n)
         m=tagslv_rby_sms(i)
         IF(m /= 0)THEN
           msr=npby(1,m)
           prec_sms(i)=prec_sms(msr)
         END IF
        END DO
C
        CALL my_barrier()
C 
      END IF
C-----------------------------------
C
      IF(nfxvel > 0)THEN
        IF(itask==0)THEN
          it=0
          CALL sms_fixvel(ibfv   ,a       ,v        ,npc    ,tf     ,
     2                    vel    ,ms      ,x        ,skew   ,sensor_tab,
     3                    weight  ,d        ,iframe ,xframe ,nsensor ,
     4                    it     ,prec_sms,nodxi_sms,cptreac,nodreac,
     5                    fthreac,am      ,vr       ,dr     ,in     ,
     6                    rby    ,wfext)
        END IF
C
        CALL my_barrier
C
      END IF
C-----------------------------------
      IF(njoint > 0)THEN
        CALL sms_cjoint_0(a    ,am    ,v ,vr,x    ,
     2                    fsav ,ljoint,ms,in,iadcj,
     3                    fr_cj,cjwork,tag_lnk_sms(nrlink+nlink+1),
     .                          prec_sms,itask)
C
        CALL my_barrier
C
      END IF
C-----------------------------------
      IF(nadmesh/=0)THEN
        IF(itask==0)THEN
          CALL sms_admesh_0(a, prec_sms, ixc, ixtg,sh4tree  ,
     .                        sh3tree  )
        END IF
C
        CALL my_barrier
C
      END IF
C-----------------------------------
      CALL sms_pcg(timers, nodft  ,nodlt    ,nnz_sms,jad_sms  ,
     2            jdi_sms ,diag_sms ,lt_sms ,a      ,isp      ,
     3            x_sms   ,p_sms    ,z_sms  ,y_sms  ,prec_sms ,
     4            nodft1_sms,nodlt1_sms,indx1_sms,icodt ,icodr ,
     5            iskew   ,skew     ,itask  ,nodxi_sms,iad_elem,
     6            fr_elem ,weight   ,ibfv   ,vel      ,npc     ,
     7            tf       ,v      ,x       ,d        ,sensor_tab,
     8            iframe ,xframe  ,jadi_sms ,jdii_sms ,nsensor   ,
     9            lti_sms  ,fr_sms ,fr_rms  ,list_sms ,list_rms,
     a            mskyi_fi_sms,vfi ,iskyi_sms,mskyi_sms,
     b            res_sms  ,ilink  ,rlink   ,fr_rl    ,frl6    ,
     c            nnlink   ,lnlink ,fr_ll   ,fnl6     ,ms      ,
     d            tag_lnk_sms,itab ,fsav    ,ljoint   ,iadcj   ,
     e            fr_cj    ,cjwork ,frl     ,fnl      ,nprw    ,
     f            lprw     ,rwbuf  ,rwsav   ,fopt     ,fr_wall ,
     g            irwl_work,nrwl_sms,frea   ,intstamp ,imv     ,
     h            mv       ,mv6     ,mw6    ,kinet    ,ixc     ,
     i            ixtg     ,sh4tree ,sh3tree,cptreac  ,nodreac ,
     j            fthreac  ,frwl6  ,am       ,vr      ,
     k            dr       ,in      ,rby    ,npby     ,lpby    ,
     l      tagmsr_rby_sms ,irbe2   ,lrbe2  ,iad_rbe2 ,fr_rbe2m,
     m            nmrbe2   ,r2size  ,irbe3  ,lrbe3    ,frbe3   ,
     n            iad_rbe3m ,fr_rbe3m ,fr_rbe3mp,rrbe3,rrbe3_pon,
     o            prec_sms3,diag_sms3,iad_rby  ,fr_rby6    ,rby6,
     p            tagslv_rby_sms,r3size,nodft2_sms,nodlt2_sms,indx2_sms,
     q            nodii_sms,ibcscyc ,lbcscyc ,wfext,ams_work )
C
      CALL my_barrier
C
c      DT15=HALF*DT1
c      DT25=HALF*DT2
c      TFEXTT=ERRTE_SMS
c      ERRTET =ZERO
c      DO N=NODFT1_SMS,NODLT1_SMS
c        I = INDX1_SMS(N)
c        VX = V(1,I)+DT05*A(1,I)
c        VY = V(2,I)+DT05*A(2,I)
c        VZ = V(3,I)+DT05*A(3,I)
c        MVX = RES_SMS(1,I)
c        MVY = RES_SMS(2,I)
c        MVZ = RES_SMS(3,I)
c        DW = (VX*MVX + VY*MVY + VZ*MVZ)*WEIGHT(I)
c        TFEXTT =TFEXTT + DT15*DW
c        ERRTET =ERRTET  + DT25*DW
c      END DO
c#include "lockon.inc"
c      TFEXT=TFEXT+TFEXTT
c      ERRTE_SMS=ERRTE_SMS+ERRTET
c#include "lockoff.inc"
c
c      CALL MY_BARRIER
c
C
      IF(njoint > 0)THEN
        CALL sms_cjoint_2(a    ,am    ,v ,vr,x    ,
     2                    ljoint,ms,in,iadcj,fr_cj,
     3                    cjwork,tag_lnk_sms(nrlink+nlink+1),itask)
C
        CALL my_barrier
C
      END IF
C-----------------------------------
C
      DO n=nodft1_sms,nodlt1_sms
        i=indx1_sms(n)
        a(1,i) = a(1,i)*ms(i)
        a(2,i) = a(2,i)*ms(i)
        a(3,i) = a(3,i)*ms(i)
      ENDDO
C
      CALL my_barrier
C
      IF(iparit/=0)THEN
        DEALLOCATE(imv, mv, mv6)
      END IF
      DEALLOCATE(mvskw)
      DEALLOCATE(vskw)
      DEALLOCATE(rskw)
      DEALLOCATE(dampskw)
C
      RETURN
      END