s8zfint_reg.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/.
!||====================================================================
!||    s8zfint_reg      ../engine/source/elements/solid/solide8z/s8zfint_reg.F
!||--- called by ------------------------------------------------------
!||    s8eforc3         ../engine/source/elements/solid/solide8e/s8eforc3.F
!||    s8sforc3         ../engine/source/elements/solid/solide8s/s8sforc3.F
!||    s8zforc3         ../engine/source/elements/solid/solide8z/s8zforc3.f
!||--- uses       -----------------------------------------------------
!||    nlocal_reg_mod   ../common_source/modules/nlocal_reg_mod.F
!||====================================================================
      SUBROUTINE s8zfint_reg(
     1   NLOC_DMG,VAR_REG, NEL,     OFFG,
     2   VOL,     NC1,     NC2,     NC3,
     3   NC4,     NC5,     NC6,     NC7,
     4   NC8,     PX1,     PX2,     PX3,
     5   PX4,     PX5,     PX6,     PX7,
     6   PX8,     PY1,     PY2,     PY3,
     7   PY4,     PY5,     PY6,     PY7,
     8   PY8,     PZ1,     PZ2,     PZ3,
     9   PZ4,     PZ5,     PZ6,     PZ7,
     A   PZ8,     IMAT,    H,       WI,
     B   IP,      ITASK,   DT2T,    VOL0,
     C   NFT)
C-----------------------------------------------
C   M o d u l e s
C----------------------------------------------- 
      USE nlocal_reg_mod
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      "parit_c.inc"
#include      "scr02_c.inc"
#include      "scr18_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NFT
      INTEGER ::  NEL,IMAT,IP,ITASK
      INTEGER, DIMENSION(NEL) :: NC1,NC2,NC3,NC4,NC5,NC6,NC7,NC8
      my_real, DIMENSION(NEL), INTENT(IN) :: 
     .   OFFG
      my_real, INTENT(INOUT)                     ::
     .   DT2T
      my_real, DIMENSION(NEL), INTENT(IN) :: 
     .   VAR_REG,PX1,PX2,PX3,PX4,PX5,PX6,PX7,PX8,
     .   PY1,PY2,PY3,PY4,PY5,PY6,PY7,PY8,PZ1,PZ2,
     .   PZ3,PZ4,PZ5,PZ6,PZ7,PZ8,VOL,H(8),VOL0
      TYPE(NLOCAL_STR_), TARGET :: NLOC_DMG 
      my_real
     .   wi
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,II,K,NNOD,N1,N2,N3,N4,N5,N6,N7,N8,L_NLOC
      my_real 
     . dx, dy, dz, l2,xi,ntvar,a,
     . b1,b2,b3,b4,b5,b6,b7,b8,
     . a1,a2,a3,a4,a5,a6,a7,a8,c1,c2,c3,c4,c5,c6,c7,c8,
     . zeta,sspnl,dtnl,le_max,maxstif,minmasscal
      my_real, DIMENSION(NEL) :: 
     . f1,f2,f3,f4,f5,f6,f7,f8,lc
      my_real, DIMENSION(:) ,ALLOCATABLE   :: 
     . btb11,btb12,btb13,btb14,btb15,btb16,btb17,btb18,
     . btb22,btb23,btb24,btb25,btb26,btb27,btb28,btb33,
     . btb34,btb35,btb36,btb37,btb38,btb44,btb45,btb46,
     . btb47,btb48,btb55,btb56,btb57,btb58,btb66,btb67,
     . btb68,btb77,btb78,btb88,sti1,sti2,sti3,sti4,sti5,
     . sti6,sti7,sti8
      INTEGER, DIMENSION(:), ALLOCATABLE   ::
     . POS1,POS2,POS3,POS4,POS5,POS6,POS7,POS8
      my_real, POINTER, DIMENSION(:) :: 
     . vnl,fnl,unl,stifnl,mass,mass0,vnl0
      ! Coefficient for non-local stability to take into account damping
      my_real, PARAMETER :: cdamp = 0.7d0
C=======================================================================
      NULLIFY(mass)
      NULLIFY(mass0)
      l2     = nloc_dmg%LEN(imat)**2
      xi     = nloc_dmg%DAMP(imat)
      nnod   = nloc_dmg%NNOD
      l_nloc = nloc_dmg%L_NLOC
      zeta   = nloc_dmg%DENS(imat)
      sspnl  = nloc_dmg%SSPNL(imat)
      le_max = nloc_dmg%LE_MAX(imat) ! Maximal length of convergence
      lc(1:nel) = zero
      vnl  => nloc_dmg%VNL(1:l_nloc)
      vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
      unl  => nloc_dmg%UNL(1:l_nloc)
      ALLOCATE(btb11(nel),btb12(nel),btb13(nel),btb14(nel),btb15(nel),
     . btb16(nel),btb17(nel),btb18(nel),btb22(nel),btb23(nel),btb24(nel),
     . btb25(nel),btb26(nel),btb27(nel),btb28(nel),btb33(nel),btb34(nel),
     . btb35(nel),btb36(nel),btb37(nel),btb38(nel),btb44(nel),btb45(nel),
     . btb46(nel),btb47(nel),btb48(nel),btb55(nel),btb56(nel),btb57(nel),
     . btb58(nel),btb66(nel),btb67(nel),btb68(nel),btb77(nel),btb78(nel),
     . btb88(nel),pos1(nel),pos2(nel),pos3(nel),pos4(nel),pos5(nel),
     . pos6(nel),pos7(nel),pos8(nel))
      ! For nodal timestep
      IF (nodadt > 0) THEN
        ! Non-local nodal stifness
        ALLOCATE(sti1(nel),sti2(nel),sti3(nel),sti4(nel),sti5(nel),sti6(nel),
     .   sti7(nel),sti8(nel))
        mass =>  nloc_dmg%MASS(1:l_nloc)
        ! non-local mass
        mass =>  nloc_dmg%MASS(1:l_nloc)
        ! Initial non-local mass
        mass0 => nloc_dmg%MASS0(1:l_nloc)
      ENDIF
c
      !-----------------------------------------------------------------------
      ! Computation of the element volume and the BtB matrix product
      !-----------------------------------------------------------------------
      ! Loop over elements
# include "vectorize.inc"
      DO i=1,nel
c
        ! Number of the element nodes
        n1 = nloc_dmg%IDXI(nc1(i))
        n2 = nloc_dmg%IDXI(nc2(i))
        n3 = nloc_dmg%IDXI(nc3(i))
        n4 = nloc_dmg%IDXI(nc4(i))
        n5 = nloc_dmg%IDXI(nc5(i))
        n6 = nloc_dmg%IDXI(nc6(i))
        n7 = nloc_dmg%IDXI(nc7(i))
        n8 = nloc_dmg%IDXI(nc8(i))
c        
        ! Recovering the position of the non-local d.o.f.s
        pos1(i) = nloc_dmg%POSI(n1)
        pos2(i) = nloc_dmg%POSI(n2)
        pos3(i) = nloc_dmg%POSI(n3)
        pos4(i) = nloc_dmg%POSI(n4) 
        pos5(i) = nloc_dmg%POSI(n5)
        pos6(i) = nloc_dmg%POSI(n6)
        pos7(i) = nloc_dmg%POSI(n7)
        pos8(i) = nloc_dmg%POSI(n8)
c        
        ! Computation of the product BtxB 
        btb11(i) = px1(i)**2 + py1(i)**2 + pz1(i)**2
        btb12(i) = px1(i)*px2(i) + py1(i)*py2(i) + pz1(i)*pz2(i)
        btb13(i) = px1(i)*px3(i) + py1(i)*py3(i) + pz1(i)*pz3(i)
        btb14(i) = px1(i)*px4(i) + py1(i)*py4(i) + pz1(i)*pz4(i)
        btb15(i) = px1(i)*px5(i) + py1(i)*py5(i) + pz1(i)*pz5(i)
        btb16(i) = px1(i)*px6(i) + py1(i)*py6(i) + pz1(i)*pz6(i)
        btb17(i) = px1(i)*px7(i) + py1(i)*py7(i) + pz1(i)*pz7(i)
        btb18(i) = px1(i)*px8(i) + py1(i)*py8(i) + pz1(i)*pz8(i)
        btb22(i) = px2(i)**2 + py2(i)**2 + pz2(i)**2
        btb23(i) = px2(i)*px3(i) + py2(i)*py3(i) + pz2(i)*pz3(i)
        btb24(i) = px2(i)*px4(i) + py2(i)*py4(i) + pz2(i)*pz4(i)
        btb25(i) = px2(i)*px5(i) + py2(i)*py5(i) + pz2(i)*pz5(i)
        btb26(i) = px2(i)*px6(i) + py2(i)*py6(i) + pz2(i)*pz6(i)
        btb27(i) = px2(i)*px7(i) + py2(i)*py7(i) + pz2(i)*pz7(i)
        btb28(i) = px2(i)*px8(i) + py2(i)*py8(i) + pz2(i)*pz8(i)
        btb33(i) = px3(i)**2 + py3(i)**2 + pz3(i)**2
        btb34(i) = px3(i)*px4(i) + py3(i)*py4(i) + pz3(i)*pz4(i)
        btb35(i) = px3(i)*px5(i) + py3(i)*py5(i) + pz3(i)*pz5(i)
        btb36(i) = px3(i)*px6(i) + py3(i)*py6(i) + pz3(i)*pz6(i)
        btb37(i) = px3(i)*px7(i) + py3(i)*py7(i) + pz3(i)*pz7(i)
        btb38(i) = px3(i)*px8(i) + py3(i)*py8(i) + pz3(i)*pz8(i)
        btb44(i) = px4(i)**2 + py4(i)**2 + pz4(i)**2
        btb45(i) = px4(i)*px5(i) + py4(i)*py5(i) + pz4(i)*pz5(i)
        btb46(i) = px4(i)*px6(i) + py4(i)*py6(i) + pz4(i)*pz6(i)
        btb47(i) = px4(i)*px7(i) + py4(i)*py7(i) + pz4(i)*pz7(i)
        btb48(i) = px4(i)*px8(i) + py4(i)*py8(i) + pz4(i)*pz8(i)
        btb55(i) = px5(i)**2 + py5(i)**2 + pz5(i)**2
        btb56(i) = px5(i)*px6(i) + py5(i)*py6(i) + pz5(i)*pz6(i)
        btb57(i) = px5(i)*px7(i) + py5(i)*py7(i) + pz5(i)*pz7(i)
        btb58(i) = px5(i)*px8(i) + py5(i)*py8(i) + pz5(i)*pz8(i)
        btb66(i) = px6(i)**2 + py6(i)**2 + pz6(i)**2
        btb67(i) = px6(i)*px7(i) + py6(i)*py7(i) + pz6(i)*pz7(i)
        btb68(i) = px6(i)*px8(i) + py6(i)*py8(i) + pz6(i)*pz8(i)
        btb77(i) = px7(i)**2 + py7(i)**2 + pz7(i)**2
        btb78(i) = px7(i)*px8(i) + py7(i)*py8(i) + pz7(i)*pz8(i)
        btb88(i) = px8(i)**2 + py8(i)**2 + pz8(i)**2
c     
      ENDDO        
c
      !-----------------------------------------------------------------------
      ! Computation of non-local forces
      !-----------------------------------------------------------------------
      ! Loop over elements
# include "vectorize.inc"
      DO i=1,nel
c
        ! If the element is not broken, normal computation
        IF (offg(i)/=zero) THEN  
          ! Computation of LEN**2*BTB*Unl product and NTN*UNL, NTN*VNL
          a1 = vol(i) * (h(1)*h(1)*unl(pos1(i)) + h(1)*h(2)*unl(pos2(i)) + h(1)*h(3)*unl(pos3(i))
     .       + h(1)*h(4)*unl(pos4(i)) + h(1)*h(5)*unl(pos5(i)) + h(1)*h(6)*unl(pos6(i))
     .       + h(1)*h(7)*unl(pos7(i)) + h(1)*h(8)*unl(pos8(i)))
c
          IF (nodadt == 0) THEN 
            a1 = a1 + vol(i) * xi * (h(1)*h(1)*vnl(pos1(i)) + h(1)*h(2)*vnl(pos2(i)) + h(1)*h(3)*vnl(pos3(i))
     .         + h(1)*h(4)*vnl(pos4(i)) + h(1)*h(5)*vnl(pos5(i)) + h(1)*h(6)*vnl(pos6(i))
     .         + h(1)*h(7)*vnl(pos7(i)) + h(1)*h(8)*vnl(pos8(i)))
          ELSE
            minmasscal = min(sqrt(mass(pos1(i))/mass0(pos1(i))),
     .                       sqrt(mass(pos2(i))/mass0(pos2(i))),
     .                       sqrt(mass(pos3(i))/mass0(pos3(i))),
     .                       sqrt(mass(pos4(i))/mass0(pos4(i))),
     .                       sqrt(mass(pos5(i))/mass0(pos5(i))),
     .                       sqrt(mass(pos6(i))/mass0(pos6(i))),
     .                       sqrt(mass(pos7(i))/mass0(pos7(i))),
     .                       sqrt(mass(pos8(i))/mass0(pos8(i))))
            a1 = a1 + vol(i) * xi * (h(1)*h(1)*minmasscal*vnl(pos1(i)) + 
     .                               h(1)*h(2)*minmasscal*vnl(pos2(i)) + 
     .                               h(1)*h(3)*minmasscal*vnl(pos3(i)) + 
     .                               h(1)*h(4)*minmasscal*vnl(pos4(i)) + 
     .                               h(1)*h(5)*minmasscal*vnl(pos5(i)) + 
     .                               h(1)*h(6)*minmasscal*vnl(pos6(i)) + 
     .                               h(1)*h(7)*minmasscal*vnl(pos7(i)) + 
     .                               h(1)*h(8)*minmasscal*vnl(pos8(i)))
          ENDIF
c        
          b1 = l2 * vol(i) * ( btb11(i)*unl(pos1(i)) + btb12(i)*unl(pos2(i)) 
     .       + btb13(i)*unl(pos3(i)) + btb14(i)*unl(pos4(i)) + btb15(i)*unl(pos5(i))
     .       + btb16(i)*unl(pos6(i)) + btb17(i)*unl(pos7(i)) + btb18(i)*unl(pos8(i)))
c     
          c1 = vol(i) * h(1) * var_reg(i)     
c     
          a2 = vol(i) * (h(2)*h(1)*unl(pos1(i)) + h(2)*h(2)*unl(pos2(i)) + h(2)*h(3)*unl(pos3(i))
     .       + h(2)*h(4)*unl(pos4(i)) + h(2)*h(5)*unl(pos5(i)) + h(2)*h(6)*unl(pos6(i))
     .       + h(2)*h(7)*unl(pos7(i)) + h(2)*h(8)*unl(pos8(i)))
c
          IF (nodadt == 0) THEN 
            a2 = a2 + vol(i) * xi * (h(2)*h(1)*vnl(pos1(i)) + h(2)*h(2)*vnl(pos2(i)) + h(2)*h(3)*vnl(pos3(i))
     .         + h(2)*h(4)*vnl(pos4(i)) + h(2)*h(5)*vnl(pos5(i)) + h(2)*h(6)*vnl(pos6(i))
     .         + h(2)*h(7)*vnl(pos7(i)) + h(2)*h(8)*vnl(pos8(i)))
          ELSE
            a2 = a2 + vol(i) * xi * (h(2)*h(1)*minmasscal*vnl(pos1(i)) + 
     .                               h(2)*h(2)*minmasscal*vnl(pos2(i)) + 
     .                               h(2)*h(3)*minmasscal*vnl(pos3(i)) + 
     .                               h(2)*h(4)*minmasscal*vnl(pos4(i)) + 
     .                               h(2)*h(5)*minmasscal*vnl(pos5(i)) + 
     .                               h(2)*h(6)*minmasscal*vnl(pos6(i)) + 
     .                               h(2)*h(7)*minmasscal*vnl(pos7(i)) + 
     .                               h(2)*h(8)*minmasscal*vnl(pos8(i)))
          ENDIF
c        
          b2 = l2 * vol(i) * ( btb12(i)*unl(pos1(i)) + btb22(i)*unl(pos2(i)) 
     .       + btb23(i)*unl(pos3(i)) + btb24(i)*unl(pos4(i)) + btb25(i)*unl(pos5(i))
     .       + btb26(i)*unl(pos6(i)) + btb27(i)*unl(pos7(i)) + btb28(i)*unl(pos8(i)))
c     
          c2 = vol(i) * h(2) * var_reg(i) 
c     
          a3 = vol(i) * (h(3)*h(1)*unl(pos1(i)) + h(3)*h(2)*unl(pos2(i)) + h(3)*h(3)*unl(pos3(i))
     .       + h(3)*h(4)*unl(pos4(i)) + h(3)*h(5)*unl(pos5(i)) + h(3)*h(6)*unl(pos6(i))
     .       + h(3)*h(7)*unl(pos7(i)) + h(3)*h(8)*unl(pos8(i)))
c
          IF (nodadt == 0) THEN 
            a3 = a3 + vol(i) * xi * (h(3)*h(1)*vnl(pos1(i)) + h(3)*h(2)*vnl(pos2(i)) + h(3)*h(3)*vnl(pos3(i))
     .         + h(3)*h(4)*vnl(pos4(i)) + h(3)*h(5)*vnl(pos5(i)) + h(3)*h(6)*vnl(pos6(i))
     .         + h(3)*h(7)*vnl(pos7(i)) + h(3)*h(8)*vnl(pos8(i))) 
          ELSE
            a3 = a3 + vol(i) * xi * (h(3)*h(1)*minmasscal*vnl(pos1(i)) + 
     .                               h(3)*h(2)*minmasscal*vnl(pos2(i)) + 
     .                               h(3)*h(3)*minmasscal*vnl(pos3(i)) + 
     .                               h(3)*h(4)*minmasscal*vnl(pos4(i)) + 
     .                               h(3)*h(5)*minmasscal*vnl(pos5(i)) + 
     .                               h(3)*h(6)*minmasscal*vnl(pos6(i)) + 
     .                               h(3)*h(7)*minmasscal*vnl(pos7(i)) + 
     .                               h(3)*h(8)*minmasscal*vnl(pos8(i)))
          ENDIF
c        
          b3 = l2 * vol(i) * ( btb13(i)*unl(pos1(i)) + btb23(i)*unl(pos2(i)) 
     .       + btb33(i)*unl(pos3(i)) + btb34(i)*unl(pos4(i)) + btb35(i)*unl(pos5(i))
     .       + btb36(i)*unl(pos6(i)) + btb37(i)*unl(pos7(i)) + btb38(i)*unl(pos8(i)))
c     
          c3 = vol(i) * h(3) * var_reg(i) 
c     
          a4 = vol(i) * (h(4)*h(1)*unl(pos1(i)) + h(4)*h(2)*unl(pos2(i)) + h(4)*h(3)*unl(pos3(i))
     .       + h(4)*h(4)*unl(pos4(i)) + h(4)*h(5)*unl(pos5(i)) + h(4)*h(6)*unl(pos6(i))
     .       + h(4)*h(7)*unl(pos7(i)) + h(4)*h(8)*unl(pos8(i)))
c
          IF (nodadt == 0) THEN 
            a4 = a4 + vol(i) * xi * (h(4)*h(1)*vnl(pos1(i)) + h(4)*h(2)*vnl(pos2(i)) + h(4)*h(3)*vnl(pos3(i))
     .         + h(4)*h(4)*vnl(pos4(i)) + h(4)*h(5)*vnl(pos5(i)) + h(4)*h(6)*vnl(pos6(i))
     .         + h(4)*h(7)*vnl(pos7(i)) + h(4)*h(8)*vnl(pos8(i)))  
          ELSE
            a4 = a4 + vol(i) * xi * (h(4)*h(1)*minmasscal*vnl(pos1(i)) + 
     .                               h(4)*h(2)*minmasscal*vnl(pos2(i)) + 
     .                               h(4)*h(3)*minmasscal*vnl(pos3(i)) + 
     .                               h(4)*h(4)*minmasscal*vnl(pos4(i)) + 
     .                               h(4)*h(5)*minmasscal*vnl(pos5(i)) + 
     .                               h(4)*h(6)*minmasscal*vnl(pos6(i)) + 
     .                               h(4)*h(7)*minmasscal*vnl(pos7(i)) + 
     .                               h(4)*h(8)*minmasscal*vnl(pos8(i)))
          ENDIF
c        
          b4 = l2 * vol(i) * ( btb14(i)*unl(pos1(i)) + btb24(i)*unl(pos2(i)) 
     .       + btb34(i)*unl(pos3(i)) + btb44(i)*unl(pos4(i)) + btb45(i)*unl(pos5(i))
     .       + btb46(i)*unl(pos6(i)) + btb47(i)*unl(pos7(i)) + btb48(i)*unl(pos8(i)))
c     
          c4 = vol(i) * h(4) * var_reg(i)      
c     
          a5 = vol(i) * (h(5)*h(1)*unl(pos1(i)) + h(5)*h(2)*unl(pos2(i)) + h(5)*h(3)*unl(pos3(i))
     .       + h(5)*h(4)*unl(pos4(i)) + h(5)*h(5)*unl(pos5(i)) + h(5)*h(6)*unl(pos6(i))
     .       + h(5)*h(7)*unl(pos7(i)) + h(5)*h(8)*unl(pos8(i)))
c
          IF (nodadt == 0) THEN
            a5 = a5 + vol(i) * xi * (h(5)*h(1)*vnl(pos1(i)) + h(5)*h(2)*vnl(pos2(i)) + h(5)*h(3)*vnl(pos3(i))
     .         + h(5)*h(4)*vnl(pos4(i)) + h(5)*h(5)*vnl(pos5(i)) + h(5)*h(6)*vnl(pos6(i))
     .         + h(5)*h(7)*vnl(pos7(i)) + h(5)*h(8)*vnl(pos8(i)))  
          ELSE
            a5 = a5 + vol(i) * xi * (h(5)*h(1)*minmasscal*vnl(pos1(i)) + 
     .                               h(5)*h(2)*minmasscal*vnl(pos2(i)) + 
     .                               h(5)*h(3)*minmasscal*vnl(pos3(i)) + 
     .                               h(5)*h(4)*minmasscal*vnl(pos4(i)) + 
     .                               h(5)*h(5)*minmasscal*vnl(pos5(i)) + 
     .                               h(5)*h(6)*minmasscal*vnl(pos6(i)) + 
     .                               h(5)*h(7)*minmasscal*vnl(pos7(i)) + 
     .                               h(5)*h(8)*minmasscal*vnl(pos8(i)))
          ENDIF
c        
          b5 = l2 * vol(i) * ( btb15(i)*unl(pos1(i)) + btb25(i)*unl(pos2(i)) 
     .       + btb35(i)*unl(pos3(i)) + btb45(i)*unl(pos4(i)) + btb55(i)*unl(pos5(i))
     .       + btb56(i)*unl(pos6(i)) + btb57(i)*unl(pos7(i)) + btb58(i)*unl(pos8(i)))
c     
          c5 = vol(i) * h(5) * var_reg(i)      
c     
          a6 = vol(i) * (h(6)*h(1)*unl(pos1(i)) + h(6)*h(2)*unl(pos2(i)) + h(6)*h(3)*unl(pos3(i))
     .       + h(6)*h(4)*unl(pos4(i)) + h(6)*h(5)*unl(pos5(i)) + h(6)*h(6)*unl(pos6(i))
     .       + h(6)*h(7)*unl(pos7(i)) + h(6)*h(8)*unl(pos8(i)))
c
          IF (nodadt == 0) THEN
            a6 = a6 + vol(i) * xi * (h(6)*h(1)*vnl(pos1(i)) + h(6)*h(2)*vnl(pos2(i)) + h(6)*h(3)*vnl(pos3(i))
     .         + h(6)*h(4)*vnl(pos4(i)) + h(6)*h(5)*vnl(pos5(i)) + h(6)*h(6)*vnl(pos6(i))
     .         + h(6)*h(7)*vnl(pos7(i)) + h(6)*h(8)*vnl(pos8(i))) 
          ELSE
            a6 = a6 + vol(i) * xi * (h(6)*h(1)*minmasscal*vnl(pos1(i)) + 
     .                               h(6)*h(2)*minmasscal*vnl(pos2(i)) + 
     .                               h(6)*h(3)*minmasscal*vnl(pos3(i)) + 
     .                               h(6)*h(4)*minmasscal*vnl(pos4(i)) + 
     .                               h(6)*h(5)*minmasscal*vnl(pos5(i)) + 
     .                               h(6)*h(6)*minmasscal*vnl(pos6(i)) + 
     .                               h(6)*h(7)*minmasscal*vnl(pos7(i)) + 
     .                               h(6)*h(8)*minmasscal*vnl(pos8(i)))
          ENDIF
c        
          b6 = l2 * vol(i) * ( btb16(i)*unl(pos1(i)) + btb26(i)*unl(pos2(i)) 
     .       + btb36(i)*unl(pos3(i)) + btb46(i)*unl(pos4(i)) + btb56(i)*unl(pos5(i))
     .       + btb66(i)*unl(pos6(i)) + btb67(i)*unl(pos7(i)) + btb68(i)*unl(pos8(i)))
c     
          c6 = vol(i) * h(6) * var_reg(i)      
c     
          a7 = vol(i) * (h(7)*h(1)*unl(pos1(i)) + h(7)*h(2)*unl(pos2(i)) + h(7)*h(3)*unl(pos3(i))
     .       + h(7)*h(4)*unl(pos4(i)) + h(7)*h(5)*unl(pos5(i)) + h(7)*h(6)*unl(pos6(i))
     .       + h(7)*h(7)*unl(pos7(i)) + h(7)*h(8)*unl(pos8(i)))
c
          IF (nodadt == 0) THEN
            a7 = a7 + vol(i) * xi * (h(7)*h(1)*vnl(pos1(i)) + h(7)*h(2)*vnl(pos2(i)) + h(7)*h(3)*vnl(pos3(i))
     .         + h(7)*h(4)*vnl(pos4(i)) + h(7)*h(5)*vnl(pos5(i)) + h(7)*h(6)*vnl(pos6(i))
     .         + h(7)*h(7)*vnl(pos7(i)) + h(7)*h(8)*vnl(pos8(i)))
          ELSE
            a7 = a7 + vol(i) * xi * (h(7)*h(1)*minmasscal*vnl(pos1(i)) + 
     .                               h(7)*h(2)*minmasscal*vnl(pos2(i)) + 
     .                               h(7)*h(3)*minmasscal*vnl(pos3(i)) + 
     .                               h(7)*h(4)*minmasscal*vnl(pos4(i)) + 
     .                               h(7)*h(5)*minmasscal*vnl(pos5(i)) + 
     .                               h(7)*h(6)*minmasscal*vnl(pos6(i)) + 
     .                               h(7)*h(7)*minmasscal*vnl(pos7(i)) + 
     .                               h(7)*h(8)*minmasscal*vnl(pos8(i)))
          ENDIF  
c        
          b7 = l2 * vol(i) * ( btb17(i)*unl(pos1(i)) + btb27(i)*unl(pos2(i)) 
     .       + btb37(i)*unl(pos3(i)) + btb47(i)*unl(pos4(i)) + btb57(i)*unl(pos5(i))
     .       + btb67(i)*unl(pos6(i)) + btb77(i)*unl(pos7(i)) + btb78(i)*unl(pos8(i)))
c     
          c7 = vol(i) * h(7) * var_reg(i)      
c     
          a8 = vol(i) * (h(8)*h(1)*unl(pos1(i)) + h(8)*h(2)*unl(pos2(i)) + h(8)*h(3)*unl(pos3(i))
     .       + h(8)*h(4)*unl(pos4(i)) + h(8)*h(5)*unl(pos5(i)) + h(8)*h(6)*unl(pos6(i))
     .       + h(8)*h(7)*unl(pos7(i)) + h(8)*h(8)*unl(pos8(i)))
c
          IF (nodadt == 0) THEN
            a8 = a8 + vol(i) * xi * (h(8)*h(1)*vnl(pos1(i)) + h(8)*h(2)*vnl(pos2(i)) + h(8)*h(3)*vnl(pos3(i))
     .         + h(8)*h(4)*vnl(pos4(i)) + h(8)*h(5)*vnl(pos5(i)) + h(8)*h(6)*vnl(pos6(i))
     .         + h(8)*h(7)*vnl(pos7(i)) + h(8)*h(8)*vnl(pos8(i)))  
          ELSE
            a8 = a8 + vol(i) * xi * (h(8)*h(1)*minmasscal*vnl(pos1(i)) + 
     .                               h(8)*h(2)*minmasscal*vnl(pos2(i)) + 
     .                               h(8)*h(3)*minmasscal*vnl(pos3(i)) + 
     .                               h(8)*h(4)*minmasscal*vnl(pos4(i)) + 
     .                               h(8)*h(5)*minmasscal*vnl(pos5(i)) + 
     .                               h(8)*h(6)*minmasscal*vnl(pos6(i)) + 
     .                               h(8)*h(7)*minmasscal*vnl(pos7(i)) + 
     .                               h(8)*h(8)*minmasscal*vnl(pos8(i)))
          ENDIF
c        
          b8 = l2 * vol(i) * ( btb18(i)*unl(pos1(i)) + btb28(i)*unl(pos2(i)) 
     .       + btb38(i)*unl(pos3(i)) + btb48(i)*unl(pos4(i)) + btb58(i)*unl(pos5(i))
     .       + btb68(i)*unl(pos6(i)) + btb78(i)*unl(pos7(i)) + btb88(i)*unl(pos8(i)))
c     
          c8 = vol(i) * h(8) * var_reg(i)      
c
c       Fint = Vol * (L*L * BT*B*Unl + NT*N*Unl + Damp*NT*N*Vnl - NT*Vreg  )    
c       L    = longueur stocke ds la structure 
c       Unl  = cumul nodal def
c       vnl  = vitesse non locale
c          
          f1(i) = a1 + b1 - c1
          f2(i) = a2 + b2 - c2
          f3(i) = a3 + b3 - c3
          f4(i) = a4 + b4 - c4
          f5(i) = a5 + b5 - c5
          f6(i) = a6 + b6 - c6
          f7(i) = a7 + b7 - c7
          f8(i) = a8 + b8 - c8
c
          ! Computing nodal equivalent stiffness
          IF (nodadt > 0) THEN 
            sti1(i) = (abs(l2*btb11(i)  + h(1)*h(1)) + abs(l2*btb12(i)  + h(1)*h(2)) + abs(l2*btb13(i)  + h(1)*h(3)) +
     .                 abs(l2*btb14(i)  + h(1)*h(4)) + abs(l2*btb15(i)  + h(1)*h(5)) + abs(l2*btb16(i)  + h(1)*h(6)) +
     .                 abs(l2*btb17(i)  + h(1)*h(7)) + abs(l2*btb18(i)  + h(1)*h(8)))*vol(i)
            sti2(i) = (abs(l2*btb12(i)  + h(2)*h(1)) + abs(l2*btb22(i)  + h(2)*h(2)) + abs(l2*btb23(i)  + h(2)*h(3)) +
     .                 abs(l2*btb24(i)  + h(2)*h(4)) + abs(l2*btb25(i)  + h(2)*h(5)) + abs(l2*btb26(i)  + h(2)*h(6)) +
     .                 abs(l2*btb27(i)  + h(2)*h(7)) + abs(l2*btb28(i)  + h(2)*h(8)))*vol(i)
            sti3(i) = (abs(l2*btb13(i)  + h(3)*h(1)) + abs(l2*btb23(i)  + h(3)*h(2)) + abs(l2*btb33(i)  + h(3)*h(3)) +
     .                 abs(l2*btb34(i)  + h(3)*h(4)) + abs(l2*btb35(i)  + h(3)*h(5)) + abs(l2*btb36(i)  + h(3)*h(6)) +
     .                 abs(l2*btb37(i)  + h(3)*h(7)) + abs(l2*btb38(i)  + h(3)*h(8)))*vol(i)
            sti4(i) = (abs(l2*btb14(i)  + h(4)*h(1)) + abs(l2*btb24(i)  + h(4)*h(2)) + abs(l2*btb34(i)  + h(4)*h(3)) +
     .                 abs(l2*btb44(i)  + h(4)*h(4)) + abs(l2*btb45(i)  + h(4)*h(5)) + abs(l2*btb46(i)  + h(4)*h(6)) +
     .                 abs(l2*btb47(i)  + h(4)*h(7)) + abs(l2*btb48(i)  + h(4)*h(8)))*vol(i)
            sti5(i) = (abs(l2*btb15(i)  + h(5)*h(1)) + abs(l2*btb25(i)  + h(5)*h(2)) + abs(l2*btb35(i)  + h(5)*h(3)) +
     .                 abs(l2*btb45(i)  + h(5)*h(4)) + abs(l2*btb55(i)  + h(5)*h(5)) + abs(l2*btb56(i)  + h(5)*h(6)) +
     .                 abs(l2*btb57(i)  + h(5)*h(7)) + abs(l2*btb58(i)  + h(5)*h(8)))*vol(i)
            sti6(i) = (abs(l2*btb16(i)  + h(6)*h(1)) + abs(l2*btb26(i)  + h(6)*h(2)) + abs(l2*btb36(i)  + h(6)*h(3)) +
     .                 abs(l2*btb46(i)  + h(6)*h(4)) + abs(l2*btb56(i)  + h(6)*h(5)) + abs(l2*btb66(i)  + h(6)*h(6)) +
     .                 abs(l2*btb67(i)  + h(6)*h(7)) + abs(l2*btb68(i)  + h(6)*h(8)))*vol(i)
            sti7(i) = (abs(l2*btb17(i)  + h(7)*h(1)) + abs(l2*btb27(i)  + h(7)*h(2)) + abs(l2*btb37(i)  + h(7)*h(3)) +
     .                 abs(l2*btb47(i)  + h(7)*h(4)) + abs(l2*btb57(i)  + h(7)*h(5)) + abs(l2*btb67(i)  + h(7)*h(6)) +
     .                 abs(l2*btb77(i)  + h(7)*h(7)) + abs(l2*btb78(i)  + h(7)*h(8)))*vol(i)
            sti8(i) = (abs(l2*btb18(i)  + h(8)*h(1)) + abs(l2*btb28(i)  + h(8)*h(2)) + abs(l2*btb38(i)  + h(8)*h(3)) +
     .                 abs(l2*btb48(i)  + h(8)*h(4)) + abs(l2*btb58(i)  + h(8)*h(5)) + abs(l2*btb68(i)  + h(8)*h(6)) +
     .                 abs(l2*btb78(i)  + h(8)*h(7)) + abs(l2*btb88(i)  + h(8)*h(8)))*vol(i)
          ENDIF
c
        ! If the element is broken
        ELSE
c
          ! Initial element characteristic length
          lc(i) = ((wi/eight)*vol0(i))**third 
c
          IF (nodadt > 0) THEN           
            ! Non-local absorbing forces
            f1(i) = sqrt(mass(pos1(i))/mass0(pos1(i)))*h(1)*zeta*sspnl*half*
     .                                (vnl(pos1(i))+vnl0(pos1(i)))*(three/four)*(lc(i)**2)
            f2(i) = sqrt(mass(pos2(i))/mass0(pos2(i)))*h(2)*zeta*sspnl*half*
     .                                (vnl(pos2(i))+vnl0(pos2(i)))*(three/four)*(lc(i)**2)
            f3(i) = sqrt(mass(pos3(i))/mass0(pos3(i)))*h(3)*zeta*sspnl*half*
     .                                (vnl(pos3(i))+vnl0(pos3(i)))*(three/four)*(lc(i)**2)
            f4(i) = sqrt(mass(pos4(i))/mass0(pos4(i)))*h(4)*zeta*sspnl*half*
     .                                (vnl(pos4(i))+vnl0(pos4(i)))*(three/four)*(lc(i)**2)
            f5(i) = sqrt(mass(pos5(i))/mass0(pos5(i)))*h(5)*zeta*sspnl*half*
     .                                (vnl(pos5(i))+vnl0(pos5(i)))*(three/four)*(lc(i)**2)
            f6(i) = sqrt(mass(pos6(i))/mass0(pos6(i)))*h(6)*zeta*sspnl*half*
     .                                (vnl(pos6(i))+vnl0(pos6(i)))*(three/four)*(lc(i)**2)
            f7(i) = sqrt(mass(pos7(i))/mass0(pos7(i)))*h(7)*zeta*sspnl*half*
     .                                (vnl(pos7(i))+vnl0(pos7(i)))*(three/four)*(lc(i)**2)
            f8(i) = sqrt(mass(pos8(i))/mass0(pos8(i)))*h(8)*zeta*sspnl*half*
     .                                (vnl(pos8(i))+vnl0(pos8(i)))*(three/four)*(lc(i)**2)
            ! Computing nodal equivalent stiffness
            sti1(i) = em20
            sti2(i) = em20
            sti3(i) = em20
            sti4(i) = em20
            sti5(i) = em20
            sti6(i) = em20
            sti7(i) = em20
            sti8(i) = em20
          ELSE
            ! Non-local absorbing forces
            f1(i) = h(1)*zeta*sspnl*half*(vnl(pos1(i))+vnl0(pos1(i)))*(three/four)*(lc(i)**2)
            f2(i) = h(2)*zeta*sspnl*half*(vnl(pos2(i))+vnl0(pos2(i)))*(three/four)*(lc(i)**2)
            f3(i) = h(3)*zeta*sspnl*half*(vnl(pos3(i))+vnl0(pos3(i)))*(three/four)*(lc(i)**2)
            f4(i) = h(4)*zeta*sspnl*half*(vnl(pos4(i))+vnl0(pos4(i)))*(three/four)*(lc(i)**2)
            f5(i) = h(5)*zeta*sspnl*half*(vnl(pos5(i))+vnl0(pos5(i)))*(three/four)*(lc(i)**2)
            f6(i) = h(6)*zeta*sspnl*half*(vnl(pos6(i))+vnl0(pos6(i)))*(three/four)*(lc(i)**2)
            f7(i) = h(7)*zeta*sspnl*half*(vnl(pos7(i))+vnl0(pos7(i)))*(three/four)*(lc(i)**2)
            f8(i) = h(8)*zeta*sspnl*half*(vnl(pos8(i))+vnl0(pos8(i)))*(three/four)*(lc(i)**2)
          ENDIF
        ENDIF
      ENDDO
c-----------------------------------------------------------------------
c     Assemblage       
c-----------------------------------------------------------------------
      IF (iparit == 0) THEN 
c
        fnl => nloc_dmg%FNL(1:l_nloc,itask+1)
        IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1) ! Non-local equivalent nodal stiffness
        DO i=1,nel
          ! Assembling the forces in the classic way 
          fnl(pos1(i)) = fnl(pos1(i)) - f1(i)
          fnl(pos2(i)) = fnl(pos2(i)) - f2(i)
          fnl(pos3(i)) = fnl(pos3(i)) - f3(i)
          fnl(pos4(i)) = fnl(pos4(i)) - f4(i)
          fnl(pos5(i)) = fnl(pos5(i)) - f5(i)
          fnl(pos6(i)) = fnl(pos6(i)) - f6(i)
          fnl(pos7(i)) = fnl(pos7(i)) - f7(i)
          fnl(pos8(i)) = fnl(pos8(i)) - f8(i)
          IF (nodadt > 0) THEN
            ! Spectral radius of stiffness matrix
            maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i),sti5(i),sti6(i),sti7(i),sti8(i))
            ! Computing nodal stiffness
            stifnl(pos1(i)) = stifnl(pos1(i)) + maxstif
            stifnl(pos2(i)) = stifnl(pos2(i)) + maxstif
            stifnl(pos3(i)) = stifnl(pos3(i)) + maxstif
            stifnl(pos4(i)) = stifnl(pos4(i)) + maxstif
            stifnl(pos5(i)) = stifnl(pos5(i)) + maxstif
            stifnl(pos6(i)) = stifnl(pos6(i)) + maxstif
            stifnl(pos7(i)) = stifnl(pos7(i)) + maxstif
            stifnl(pos8(i)) = stifnl(pos8(i)) + maxstif
          ENDIF
        ENDDO
c
      ! If PARITH/ON
      ELSE
c
        DO i=1,nel
          ii  = i + nft
c
          ! Spectral radius of stiffness matrix
          IF (nodadt > 0) THEN
            maxstif = max(sti1(i),sti2(i),sti3(i),sti4(i),sti5(i),sti6(i),sti7(i),sti8(i))
          ENDIF
c
          k = nloc_dmg%IADS(1,ii)
          IF (ip == 1) THEN 
            nloc_dmg%FSKY(k,1) = -f1(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
          ELSE
            nloc_dmg%FSKY(k,1) = nloc_dmg%FSKY(k,1) - f1(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = nloc_dmg%STSKY(k,1) + maxstif
          ENDIF
c
          k = nloc_dmg%IADS(2,ii)
          IF (ip == 1) THEN 
            nloc_dmg%FSKY(k,1) = -f2(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
          ELSE
            nloc_dmg%FSKY(k,1) = nloc_dmg%FSKY(k,1) - f2(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = nloc_dmg%STSKY(k,1) + maxstif
          ENDIF
c
          k = nloc_dmg%IADS(3,ii)
          IF (ip == 1) THEN 
            nloc_dmg%FSKY(k,1) = -f3(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
          ELSE
            nloc_dmg%FSKY(k,1) = nloc_dmg%FSKY(k,1) - f3(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = nloc_dmg%STSKY(k,1) + maxstif
          ENDIF
c
          k = nloc_dmg%IADS(4,ii)
          IF (ip == 1) THEN 
            nloc_dmg%FSKY(k,1) = -f4(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
          ELSE
            nloc_dmg%FSKY(k,1) = nloc_dmg%FSKY(k,1) - f4(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = nloc_dmg%STSKY(k,1) + maxstif
          ENDIF
c
          k = nloc_dmg%IADS(5,ii)
          IF (ip == 1) THEN 
            nloc_dmg%FSKY(k,1) = -f5(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
          ELSE
            nloc_dmg%FSKY(k,1) = nloc_dmg%FSKY(k,1) - f5(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = nloc_dmg%STSKY(k,1) + maxstif
          ENDIF
c
          k = nloc_dmg%IADS(6,ii)
          IF (ip == 1) THEN 
            nloc_dmg%FSKY(k,1) = -f6(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
          ELSE
            nloc_dmg%FSKY(k,1) = nloc_dmg%FSKY(k,1) - f6(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = nloc_dmg%STSKY(k,1) + maxstif
          ENDIF
c
          k = nloc_dmg%IADS(7,ii)
          IF (ip == 1) THEN 
            nloc_dmg%FSKY(k,1) = -f7(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
          ELSE
            nloc_dmg%FSKY(k,1) = nloc_dmg%FSKY(k,1) - f7(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = nloc_dmg%STSKY(k,1) + maxstif
          ENDIF
c
          k = nloc_dmg%IADS(8,ii)
          IF (ip == 1) THEN 
            nloc_dmg%FSKY(k,1) = -f8(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = maxstif
          ELSE
            nloc_dmg%FSKY(k,1) = nloc_dmg%FSKY(k,1) - f8(i)
            IF (nodadt > 0) nloc_dmg%STSKY(k,1) = nloc_dmg%STSKY(k,1) + maxstif
          ENDIF
c
        ENDDO
      ENDIF      
c      
      !-----------------------------------------------------------------------
      ! Computing non-local timestep
      !-----------------------------------------------------------------------   
      IF (nodadt == 0) THEN
        DO i = 1,nel
          ! If the element is not broken, normal computation
          IF (offg(i)/=zero) THEN
            ! Non-local critical time-step in the plane
            dtnl = (two*(min(vol(i)**third,le_max))*sqrt(three*zeta))/
     .                  sqrt(twelve*l2 + (min(vol(i)**third,le_max))**2)
            ! Retaining the minimal value
            dt2t = min(dt2t,dtfac1(1)*cdamp*dtnl)
          ENDIF
        ENDDO
      ENDIF
c-----------
      ! Deallocation of tables
      IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
      IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
      IF (ALLOCATED(btb13)) DEALLOCATE(btb13)
      IF (ALLOCATED(btb14)) DEALLOCATE(btb14)
      IF (ALLOCATED(btb15)) DEALLOCATE(btb15)
      IF (ALLOCATED(btb16)) DEALLOCATE(btb16)
      IF (ALLOCATED(btb17)) DEALLOCATE(btb17)
      IF (ALLOCATED(btb18)) DEALLOCATE(btb18)
      IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
      IF (ALLOCATED(btb23)) DEALLOCATE(btb23)      
      IF (ALLOCATED(btb24)) DEALLOCATE(btb24)
      IF (ALLOCATED(btb25)) DEALLOCATE(btb25)
      IF (ALLOCATED(btb26)) DEALLOCATE(btb26)
      IF (ALLOCATED(btb27)) DEALLOCATE(btb27)
      IF (ALLOCATED(btb28)) DEALLOCATE(btb28)
      IF (ALLOCATED(btb33)) DEALLOCATE(btb33)      
      IF (ALLOCATED(btb34)) DEALLOCATE(btb34)
      IF (ALLOCATED(btb35)) DEALLOCATE(btb35)
      IF (ALLOCATED(btb36)) DEALLOCATE(btb36)
      IF (ALLOCATED(btb37)) DEALLOCATE(btb37)
      IF (ALLOCATED(btb38)) DEALLOCATE(btb38)
      IF (ALLOCATED(btb44)) DEALLOCATE(btb44)
      IF (ALLOCATED(btb45)) DEALLOCATE(btb45)
      IF (ALLOCATED(btb46)) DEALLOCATE(btb46)
      IF (ALLOCATED(btb47)) DEALLOCATE(btb47)
      IF (ALLOCATED(btb48)) DEALLOCATE(btb48)
      IF (ALLOCATED(btb55)) DEALLOCATE(btb55)
      IF (ALLOCATED(btb56)) DEALLOCATE(btb56)
      IF (ALLOCATED(btb57)) DEALLOCATE(btb57)
      IF (ALLOCATED(btb58)) DEALLOCATE(btb58)
      IF (ALLOCATED(btb66)) DEALLOCATE(btb66)
      IF (ALLOCATED(btb67)) DEALLOCATE(btb67)
      IF (ALLOCATED(btb68)) DEALLOCATE(btb68)
      IF (ALLOCATED(btb77)) DEALLOCATE(btb77)
      IF (ALLOCATED(btb78)) DEALLOCATE(btb78)
      IF (ALLOCATED(btb88)) DEALLOCATE(btb88)
      IF (ALLOCATED(pos1))  DEALLOCATE(pos1)
      IF (ALLOCATED(pos2))  DEALLOCATE(pos2) 
      IF (ALLOCATED(pos3))  DEALLOCATE(pos3)
      IF (ALLOCATED(pos4))  DEALLOCATE(pos4) 
      IF (ALLOCATED(pos5))  DEALLOCATE(pos5)
      IF (ALLOCATED(pos6))  DEALLOCATE(pos6) 
      IF (ALLOCATED(pos7))  DEALLOCATE(pos7)
      IF (ALLOCATED(pos8))  DEALLOCATE(pos8)
      IF (ALLOCATED(sti1))  DEALLOCATE(sti1)
      IF (ALLOCATED(sti2))  DEALLOCATE(sti2) 
      IF (ALLOCATED(sti3))  DEALLOCATE(sti3)
      IF (ALLOCATED(sti4))  DEALLOCATE(sti4) 
      IF (ALLOCATED(sti5))  DEALLOCATE(sti5)
      IF (ALLOCATED(sti6))  DEALLOCATE(sti6) 
      IF (ALLOCATED(sti7))  DEALLOCATE(sti7)
      IF (ALLOCATED(sti8))  DEALLOCATE(sti8)
c
      END