i17lagm.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/.
!||====================================================================
!||    i17lagm       ../engine/source/interfaces/int17/i17lagm.f
!||--- called by ------------------------------------------------------
!||    i17main       ../engine/source/interfaces/int17/i17main.F
!||--- calls      -----------------------------------------------------
!||    i17lll        ../engine/source/interfaces/int17/i17lagm.f
!||    my_barrier    ../engine/source/system/machine.f
!||--- uses       -----------------------------------------------------
!||    element_mod   ../common_source/modules/elements/element_mod.f90
!||====================================================================
      SUBROUTINE i17lagm(X    ,V     ,LLL     ,JLL   ,SLL   ,
     2                  XLL   ,CANDN ,CANDE   ,I_STOK,IXS   ,
     3                  IXS16 ,IADLL ,EMINX   ,NELES ,NELEM ,
     4                  NC    ,N_MUL_MX,ITASK ,A     ,ITIED ,
     5                  NINT  ,NKMAX ,EMINXS  ,COMNTAG)
      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      "task_c.inc"
#include      "com04_c.inc"
#include      "com08_c.inc"
      COMMON /i17globi/ie_min,ies_min
      COMMON /i17globr/vit_min
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NC,I_STOK,N_MUL_MX,ITASK,ITIED,NINT,NKMAX ,
     .  LLL(*),JLL(*),SLL(*),CANDN(*),CANDE(*),COMNTAG(*),
     .  IXS(NIXS,*),IXS16(8,*),IADLL(*),NELES(*)  ,NELEM(*) 
C     REAL
      my_real
     .  x(3,*),v(3,*),xll(*),
     .  eminx(6,*),eminxs(6,*),a(3,*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,K,IK,IE,IS,IC,NK,III(MVSIZ,17),LLT,NFT,LE,FIRST,LAST,
     .        I16,LES,IES,IIIS(MVSIZ,16),NC_SAV,IEL(MVSIZ),IESL(MVSIZ),
     .        IE_MIN,IES_MIN
      my_real
     .   xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17),
     .   xxs(mvsiz,16),yys(mvsiz,16),zzs(mvsiz,16),
     .   aa,xmin,ymin,zmin,xmax,ymax,zmax,dist,vit_min
C-----------------------------------------------
C
C
C       | M | Lt| | a    | M ao
C       |---+---| |    = |
C       | L | 0 | | la   | bo
C
C        [M] a + [L]t la = [M] ao
C        [L] a = bo
C
C         a = -[M]-1[L]t la + ao
C        [L][M]-1[L]t la  = [L] ao - bo
C
C on pose:
C        [H] = [L][M]-1[L]t
C        b  = [L] ao - bo
C
C        [H] la = b
C
C        a = ao - [M]-1[L]t la
C-----------------------------------------------
C
C        la              : LAMBDA(NC)
C        ao              : A(NUMNOD)
C        L               : XLL(NK,NC)
C        M               : MAS(NUMNOD)
C        [L][M]-1[L]t la : HLA(NC)
C        [L] ao - b      : B(NC)
C        [M]-1[L]t la    : LTLA(NUMNOD)
C
C        NC : number of contacts
C        NK: Number of node for contact (8+1.16+1.8+8.16+16)
C
C        IC : contact number (1,NC)
C        IK : local node number for a contact (1,NK)
C        I  : global node number (1,NUMNOD)
C
C        IADLL(NC)        : IAD = IADLL(IC)
C        LLL(NC*(17,51))  : I  = LLL(IAD+1,2...IADNEXT-1)
C-----------------------------------------------
C  evaluation of b:
C
C         Vs = Somme(Ni Vi)
C         Vs_ + dt As = Somme(Ni Vi_) + Somme(dt Ni Ai)
C         Somme(dt Ni Ai) - dt As =  Vs_ -Somme(Ni Vi_)
C         [L] = dt {N1,N2,..,N15,-1}
C         bo = [L] a = -[L]/dt v_
C         b  = [L] ao - bo
C         b  = [L] ao + [L]/dt v_ = [L] (v_ + ao dt)/dt
C-----------------------------------------------
C                 b  = [L] vo+/dt   +   vout
C-----------------------------------------------
C-----------------------------------------------------------------------
C     loop over contact candidates
C-----------------------------------------------------------------------
      nc_sav = nc
      vit_min = zero
      ie_min  = 99999999
      ies_min = 99999999
      CALL my_barrier
      first = 1 + i_stok * itask / nthread
      last = i_stok*(itask+1) / nthread
      llt = 0
      nft=llt+1
      DO ic=first,last
        le  = cande(ic)
        les = candn(ic)
        ie  = nelem(le)
        ies = neles(les)
C-----------------------------------------------------------------------
C       test if inside the box
C-----------------------------------------------------------------------
        IF(le >0.AND.les>0.AND.
     .           eminxs(4,les)>eminx(1,le).AND.
     .           eminxs(5,les)>eminx(2,le).AND.
     .           eminxs(6,les)>eminx(3,le).AND.
     .           eminxs(1,les)<eminx(4,le).AND.
     .           eminxs(2,les)<eminx(5,le).AND.
     .           eminxs(3,les)<eminx(6,le))THEN
c
c      print *, "in la boite",XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
c
          llt = llt+1
            iel(llt)=ie
            iesl(llt)=ies
          DO k=1,4
C III 1,2,3,4,5,6,7,8, 9,10,11,12,13,14,15,16
            iii(llt,k)   =ixs(k+1,ie)
            iii(llt,k+4) =ixs(k+5,ie)
            iii(llt,k+8) =ixs16(k,ie-numels8-numels10-numels20)
            iii(llt,k+12)=ixs16(k+4,ie-numels8-numels10-numels20)
C IIIS 1,2,3,4,9,10,11,12, 5,6,7,8,13,14,15,16
            iiis(llt,k)  =ixs(k+1,ies)
            iiis(llt,k+8) =ixs(k+5,ies)
            iiis(llt,k+4)=ixs16(k,ies-numels8-numels10-numels20)
            iiis(llt,k+12)=ixs16(k+4,ies-numels8-numels10-numels20)
          ENDDO
          DO k=1,16
            i = iii(llt,k)
c            XX(LLT,K)=X(1,I)+DT2*(V(1,I)+DT12*A(1,I))
c            YY(LLT,K)=X(2,I)+DT2*(V(2,I)+DT12*A(2,I))
c            ZZ(LLT,K)=X(3,I)+DT2*(V(3,I)+DT12*A(3,I))
c            XX(LLT,K)=X(1,I)
c            YY(LLT,K)=X(2,I)
c            ZZ(LLT,K)=X(3,I)
            xx(llt,k)=x(1,i)+half*dt2*(v(1,i)+dt12*a(1,i))
            yy(llt,k)=x(2,i)+half*dt2*(v(2,i)+dt12*a(2,i))
            zz(llt,k)=x(3,i)+half*dt2*(v(3,i)+dt12*a(3,i))
            i = iiis(llt,k)
            xxs(llt,k)=x(1,i)+half*dt2*(v(1,i)+dt12*a(1,i))
            yys(llt,k)=x(2,i)+half*dt2*(v(2,i)+dt12*a(2,i))
            zzs(llt,k)=x(3,i)+half*dt2*(v(3,i)+dt12*a(3,i))
          ENDDO
c
c      print *, "XX(1,1),XX(1,9)",XX(1,1),XX(1,9)
c
C-----------------------------------------------------------------------
C     calculation of [L] by mvsiz packet
C-----------------------------------------------------------------------
          IF(llt==mvsiz-1)THEN
            CALL i17lll(
     1       llt ,lll       ,jll       ,sll       ,xll     ,v       ,
     2       xx  ,yy        ,zz        ,iii       ,nc      ,iadll   ,
     3       n_mul_mx ,a    ,x         ,itied     ,nint    ,nkmax   ,
     4       xxs ,yys       ,zzs       ,iiis      ,nc_sav  ,vit_min ,
     5       iel  ,iesl     ,ie_min    ,ies_min   ,itask   ,comntag )
            nft=llt+1
            llt = 0
          ENDIF
        ELSE
c debug
          k=0
        ENDIF
      ENDDO
C-----------------------------------------------------------------------
C     calculation of [L] for last packet
C-----------------------------------------------------------------------
      IF(llt/=0) CALL i17lll(
     1       llt ,lll       ,jll       ,sll       ,xll     ,v       ,
     2       xx  ,yy        ,zz        ,iii       ,nc      ,iadll   ,
     3       n_mul_mx ,a    ,x         ,itied     ,nint    ,nkmax   ,
     4       xxs ,yys       ,zzs       ,iiis      ,nc_sav  ,vit_min ,
     5       iel  ,iesl     ,ie_min    ,ies_min   ,itask   ,comntag )
C
C-----------------------------------------------
      CALL my_barrier
      RETURN
      END
!||====================================================================
!||    i17lll    ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17lagm   ../engine/source/interfaces/int17/i17lagm.F
!||--- calls      -----------------------------------------------------
!||    i17lll4   ../engine/source/interfaces/int17/i17lagm.F
!||    i17mini   ../engine/source/interfaces/int17/i17lagm.F
!||    i17rst    ../engine/source/interfaces/int17/i17lagm.F
!||    i17vit4   ../engine/source/interfaces/int17/i17lagm.F
!||====================================================================
      SUBROUTINE i17lll(LLT  ,LLL    ,JLL   ,SLL    ,XLL   ,V      ,
     2                  XX   ,YY     ,ZZ    ,III    ,NC    ,IADLL  ,
     3                  N_MUL_MX,A   ,X     ,ITIED  ,NINT  ,NKMAX  ,
     4                  XXS  ,YYS    ,ZZS   ,IIIS   ,NC_SAV,VIT_MIN,
     5                  IE   ,IES    ,IE_MIN,IES_MIN,ITASK ,COMNTAG)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
Cbbdebug +1
#include      "comlock.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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT,NC,N_MUL_MX,ITIED,NINT ,NKMAX ,NC_SAV ,IE_MIN,IES_MIN 
      INTEGER LLL(*),JLL(*),SLL(*),ITASK,COMNTAG(*),
     .        III(MVSIZ,17),IADLL(*) ,IIIS(MVSIZ,16),IE(*) ,IES(*)
C     REAL
      my_real
     .  XLL(*),V(3,*),A(3,*)
      my_real
     .   XX(MVSIZ,17),YY(MVSIZ,17),ZZ(MVSIZ,17),X(3,*),
     .   xxs(mvsiz,16) ,yys(mvsiz,16) ,zzs(mvsiz,16) ,vit_min       
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,IK,NK,I1,I2,I3,I4,IAD,ICON,
     +        ICONT(MVSIZ,4)
      my_real
     .   vx,vy,vz,vn,aa,vv,pene
      my_real
     .   r_cm(mvsiz),t_cm(mvsiz),s_cm(mvsiz),si_s(mvsiz,8),
     .   ri_s(mvsiz,8),ti_s(mvsiz,8),
     .   nx(mvsiz), ny(mvsiz), nz(mvsiz),
     .   ni_m(mvsiz,17),ni_s(mvsiz,8) ,
     .   r_1s(mvsiz) ,r_2s(mvsiz) ,t_1s(mvsiz) ,t_2s(mvsiz),
     .   r_1m(mvsiz) ,r_2m(mvsiz) ,t_1m(mvsiz) ,t_2m(mvsiz),
     .   r_3m(mvsiz) ,r_4m(mvsiz) ,t_3m(mvsiz) ,t_4m(mvsiz),
     .   r_cs(mvsiz) ,s_cs(mvsiz) ,t_cs(mvsiz) ,vit(mvsiz),
     .   r_s(mvsiz) ,t_s(mvsiz) 
C-----------------------------------------------
C      calculation de r,s,t
C-----------------------------------------------
c
c      print *, "XX(1,1),XX(1,9)",XX(1,1),XX(1,9)
c
C-----------------------------------------------
C      central node face 1 2 3 4
C-----------------------------------------------
      DO i=1,llt
        xx(i,17) = half *(xxs(i,5) +xxs(i,6) +xxs(i,7) +xxs(i,8))
     .           - fourth*(xxs(i,1) +xxs(i,2) +xxs(i,3) +xxs(i,4))
        yy(i,17) = half *(yys(i,5) +yys(i,6) +yys(i,7) +yys(i,8))
     .           - fourth*(yys(i,1) +yys(i,2) +yys(i,3) +yys(i,4))
        zz(i,17) = half *(zzs(i,5) +zzs(i,6) +zzs(i,7) +zzs(i,8))
     .           - fourth*(zzs(i,1) +zzs(i,2) +zzs(i,3) +zzs(i,4))     
      ENDDO
      CALL i17rst(llt   ,ri_s(1,1),si_s(1,1),ti_s(1,1),ni_m  ,
     2            xx    ,yy       ,zz       )
C-----------------------------------------------
C      Central node face 5 6 7 8
C-----------------------------------------------
      DO i=1,llt
        xx(i,17) = half *(xxs(i,13)+xxs(i,14)+xxs(i,15)+xxs(i,16))
     .           - fourth*(xxs(i,9) +xxs(i,10)+xxs(i,11)+xxs(i,12))
        yy(i,17) = half *(yys(i,13)+yys(i,14)+yys(i,15)+yys(i,16))
     .           - fourth*(yys(i,9) +yys(i,10)+yys(i,11)+yys(i,12))
        zz(i,17) = half *(zzs(i,13)+zzs(i,14)+zzs(i,15)+zzs(i,16))
     .           - fourth*(zzs(i,9) +zzs(i,10)+zzs(i,11)+zzs(i,12))     
      ENDDO
      CALL i17rst(llt   ,ri_s(1,2),si_s(1,2),ti_s(1,2),ni_m  ,
     2            xx    ,yy       ,zz       )
C-----------------------------------------------
C      choix face 1 2 3 4 face 5 6 7 8 cote second
C-----------------------------------------------
      DO i=1,llt
        IF(abs(si_s(i,1))<=abs(si_s(i,2)))THEN
C         face 1 2 3 4
          s_cs(i) = -one
        ELSE
C         face 5 6 7 8
          s_cs(i) =  one
          iiis(i,1) = iiis(i,9)
          iiis(i,2) = iiis(i,10)
          iiis(i,3) = iiis(i,11)
          iiis(i,4) = iiis(i,12)
          iiis(i,5) = iiis(i,13)
          iiis(i,6) = iiis(i,14)
          iiis(i,7) = iiis(i,15)
          iiis(i,8) = iiis(i,16)
C
          xxs(i,1) = xxs(i,9)
          xxs(i,2) = xxs(i,10)
          xxs(i,3) = xxs(i,11)
          xxs(i,4) = xxs(i,12)
          xxs(i,5) = xxs(i,13)
          xxs(i,6) = xxs(i,14)
          xxs(i,7) = xxs(i,15)
          xxs(i,8) = xxs(i,16)
C
          yys(i,1) = yys(i,9)
          yys(i,2) = yys(i,10)
          yys(i,3) = yys(i,11)
          yys(i,4) = yys(i,12)
          yys(i,5) = yys(i,13)
          yys(i,6) = yys(i,14)
          yys(i,7) = yys(i,15)
          yys(i,8) = yys(i,16)
C
          zzs(i,1) = zzs(i,9)
          zzs(i,2) = zzs(i,10)
          zzs(i,3) = zzs(i,11)
          zzs(i,4) = zzs(i,12)
          zzs(i,5) = zzs(i,13)
          zzs(i,6) = zzs(i,14)
          zzs(i,7) = zzs(i,15)
          zzs(i,8) = zzs(i,16)
C
        ENDIF
      ENDDO
C-----------------------------------------------
C      calculation de SI_S=s(relatif element main) 
c      the 8 nodes of the face of the second element
C-----------------------------------------------
      DO j=1,8
        DO i=1,llt
          xx(i,17) = xxs(i,j)
          yy(i,17) = yys(i,j)
          zz(i,17) = zzs(i,j)
        ENDDO
        CALL i17rst(llt   ,ri_s(1,j),si_s(1,j),ti_s(1,j),ni_m  ,
     2              xx    ,yy    ,zz    )
      ENDDO
C-----------------------------------------------
C      calculation of the minimum distance point
C      dSI_S/dr = 0 dSI_S/dt = 0     (r,t de l'element second) 
C-----------------------------------------------
      CALL i17mini(llt   ,r_cs  ,s_cs  ,t_cs  ,ri_s  ,si_s   ,
     2             ti_s  ,ni_s  ,xxs   ,yys   ,zzs   ,xx     ,
     3             yy    ,zz    ,r_cm  ,s_cm  ,t_cm  ,nx     ,
     4             ny    ,nz    ,r_1s  ,r_2s  ,t_1s  ,t_2s   ,
     5             r_1m  ,r_2m  ,r_3m  ,r_4m  ,t_1m  ,t_2m   ,
     6             t_3m  ,t_4m  ,icont )
C-----------------------------------------------
C      Choice Facing 1 2 3 4 Face 5 6 7 8 main side
C-----------------------------------------------
      DO i=1,llt
        IF(s_cm(i)<=0.)THEN
C         face 1 2 3 4
          iii(i,5) = iii(i,9)
          iii(i,6) = iii(i,10)
          iii(i,7) = iii(i,11)
          iii(i,8) = iii(i,12)
        ELSE
C         face 5 6 7 8
          iii(i,1) = iii(i,5)
          iii(i,2) = iii(i,6)
          iii(i,3) = iii(i,7)
          iii(i,4) = iii(i,8)
          iii(i,5) = iii(i,13)
          iii(i,6) = iii(i,14)
          iii(i,7) = iii(i,15)
          iii(i,8) = iii(i,16)
        ENDIF
      ENDDO
C-----------------------------------------------------------------------
C     calculation of penetration velocities
C-----------------------------------------------------------------------
      DO i=1,llt
        vit(i) = em20
      ENDDO
      CALL i17vit4(llt       ,nint   ,v      ,a      ,iii    ,iiis   ,
     2             ni_m      ,ni_s   ,nx     ,ny     ,nz     ,vit    ,
     3             icont(1,1),r_1m   ,t_1m   ,r_1s   ,t_cs   ,s_cs   ,
     4             r_cm      ,t_cm   ,r_s    ,t_s    ,icont(1,1))
      CALL i17vit4(llt       ,nint   ,v      ,a      ,iii    ,iiis   ,
     2             ni_m      ,ni_s   ,nx     ,ny     ,nz     ,vit    ,
     3             icont(1,2),r_2m   ,t_2m   ,r_2s   ,t_cs   ,s_cs   ,
     4             r_cm      ,t_cm   ,r_s    ,t_s    ,icont(1,1))
      CALL i17vit4(llt       ,nint   ,v      ,a      ,iii    ,iiis   ,
     2             ni_m      ,ni_s   ,nx     ,ny     ,nz     ,vit    ,
     3             icont(1,3),r_3m   ,t_3m   ,r_cs   ,t_1s   ,s_cs   ,
     4             r_cm      ,t_cm   ,r_s    ,t_s    ,icont(1,1))
      CALL i17vit4(llt       ,nint   ,v      ,a      ,iii    ,iiis   ,
     2             ni_m      ,ni_s   ,nx     ,ny     ,nz     ,vit    ,
     3             icont(1,4),r_4m   ,t_4m   ,r_cs   ,t_2s   ,s_cs   ,
     4             r_cm      ,t_cm   ,r_s    ,t_s    ,icont(1,1))
c tmp !!!!!!!!!!!!!!!!!!!!!!!!!!! single contact
c      VV = 0.
c      DO I=1,LLT
c        VV = MIN(VV,VIT(I))
c      ENDDO
c      DO I=1,LLT
c        IF(VV==VIT(I))THEN
c          VV=1.
c        ELSE
c          ICONT(I,1)=0
c        ENDIF
c      ENDDO
c tmp !!!!!!!!!!!!!!!!!!!!!!!!!!!
c      DO I=1,LLT
c        PENE = 1. - ABS(S_CM(I))
c        IF(ICONT(I,1)==1.AND.PENE>PENE_MAX)THEN
c          PENE_MAX = PENE
c          NC = NC_SAV
c          J  = I
c        ENDIF
c        ICONT(I,1) = 0
c      ENDDO
c      ICONT(J,1) = 1
#include "lockon.inc"
      j=0
      DO i=1,llt
        IF(vit(i)<vit_min.OR.
     .    (vit(i)==vit_min.AND.
     .          (ie(i)<ie_min.OR.
     .          (ie(i)==ie_min.AND.ies(i)<ies_min))))THEN
          vit_min = vit(i)
          ie_min  = ie(i)
          ies_min = ies(i)
          nc = nc_sav
          j  = i
          icon=icont(i,1)
        ENDIF
        icont(i,1) = 0
      ENDDO
      IF(j/=0)THEN
        icont(j,1) = icon
cc  print *,ITASK,':',VIT(J),IE(J),IES(J),NC,J,ICON
      ENDIF
C-----------------------------------------------------------------------
C     contact calculation
C-----------------------------------------------------------------------
      CALL i17lll4(llt       ,lll    ,jll    ,sll    ,xll    ,n_mul_mx,
     2             itied     ,nint   ,nkmax  ,nc     ,v      ,a       ,
     3             iadll     ,iii    ,iiis   ,ni_m   ,ni_s   ,
     4             nx        ,ny     ,nz     ,vit    ,comntag,
     5             icont(1,1),r_cm   ,t_cm   ,r_s    ,t_s    )
c      IF(J/=0)THEN
c        do i=1,nc*16
c  print *,ITASK,':',i,':',LLL(i),JLL(i),XLL(i)
c  enddo
c      ENDIF
C
#include "lockoff.inc"
      RETURN
      END
!||====================================================================
!||    i17vit4   ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17lll    ../engine/source/interfaces/int17/i17lagm.F
!||--- calls      -----------------------------------------------------
!||    i17ni     ../engine/source/interfaces/int17/i17lagm.F
!||====================================================================
      SUBROUTINE i17vit4(LLT,NINT ,V      ,A      ,III    ,IIIS   ,
     2             NI_M   ,NI_S   ,NX     ,NY     ,NZ     ,VIT    ,
     3             ICONT  ,RM     ,TM     ,RS     ,TS     ,SM     ,
     4             R_M    ,T_M    ,R_S    ,T_S    ,ICONTN )
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT,NINT   
      INTEGER III(MVSIZ,17),IIIS(MVSIZ,16),
     +        ICONT(MVSIZ),ICONTN(MVSIZ)
C     REAL
      my_real
     .  V(3,*),A(3,*),VIT(*)
      my_real
     .   RM(MVSIZ) ,RS(MVSIZ) ,TM(MVSIZ) ,TS(MVSIZ) ,SM(MVSIZ),
     .   R_M(MVSIZ) ,R_S(MVSIZ) ,T_M(MVSIZ) ,T_S(MVSIZ) ,
     .   NI_M(MVSIZ,*) ,NI_S(MVSIZ,*),NX(MVSIZ) ,NY(MVSIZ) ,NZ(MVSIZ)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,IK,NK,I1,I2,I3,I4
      my_real
     .   vx,vy,vz,vn,aa
      my_real
     .   ni_ml(mvsiz,8) ,ni_sl(mvsiz,8)
C-----------------------------------------------------------------------
C     calculation of Ni on face 1 2 3 4 or 5 6 7 8 of the main (s=+-1)
C-----------------------------------------------------------------------
      CALL i17ni(llt,rm ,tm ,ni_ml )
C-----------------------------------------------------------------------
C     calculation of Ni on face 1 2 3 4 or 5 6 7 8 of the second (s=+-1)
C-----------------------------------------------------------------------
      CALL i17ni(llt,rs ,ts ,ni_sl )
C-----------------------------------------------
C      calculation of [L]
C-----------------------------------------------
      DO i=1,llt
C-----------------------------------------------
C       Test if contact
C-----------------------------------------------
        IF(icont(i)/=0)THEN
C
          vx = zero
          vy = zero
          vz = zero
          DO ik=1,8
c            VX = VX - (V(1,III(I,IK))+DT12*A(1,III(I,IK)))*NI_M(I,IK)
c            VY = VY - (V(2,III(I,IK))+DT12*A(2,III(I,IK)))*NI_M(I,IK)
c            VZ = VZ - (V(3,III(I,IK))+DT12*A(3,III(I,IK)))*NI_M(I,IK)
c            VX = VX + (V(1,IIIS(I,IK))+DT12*A(1,IIIS(I,IK)))*NI_S(I,IK)
c            VY = VY + (V(2,IIIS(I,IK))+DT12*A(2,IIIS(I,IK)))*NI_S(I,IK)
c            VZ = VZ + (V(3,IIIS(I,IK))+DT12*A(3,IIIS(I,IK)))*NI_S(I,IK)
            vx = vx - (v(1,iii(i,ik)))*ni_ml(i,ik)
            vy = vy - (v(2,iii(i,ik)))*ni_ml(i,ik)
            vz = vz - (v(3,iii(i,ik)))*ni_ml(i,ik)
            vx = vx + (v(1,iiis(i,ik)))*ni_sl(i,ik)
            vy = vy + (v(2,iiis(i,ik)))*ni_sl(i,ik)
            vz = vz + (v(3,iiis(i,ik)))*ni_sl(i,ik)
          ENDDO
c
          vn = nx(i)*vx + ny(i)*vy + nz(i)*vz
C-----------------------------------------------
C         Test if incoming velocity in S
C-----------------------------------------------
          IF(sm(i)*vn<=vit(i))THEN
c
c      print *, "velocity entrante",vn
c      print *, "s = ",S(I)
c
            vit(i)=sm(i)*vn
            r_s(i)=rs(i)
            r_m(i)=rm(i)
            t_s(i)=ts(i)
            t_m(i)=tm(i)
            DO ik=1,8
              ni_m(i,ik)=ni_ml(i,ik)
              ni_s(i,ik)=ni_sl(i,ik)
            ENDDO
            icontn(i)=1
          ENDIF
        ENDIF
      ENDDO
c
c      print *, "r,s,t",r(1),s(1),t(1)
C
      RETURN
      END
C
!||====================================================================
!||    i17lll4       ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17lll        ../engine/source/interfaces/int17/i17lagm.F
!||--- calls      -----------------------------------------------------
!||    ancmsg        ../engine/source/output/message/message.F
!||    arret         ../engine/source/system/arret.F
!||--- uses       -----------------------------------------------------
!||    message_mod   ../engine/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE i17lll4(LLT,LLL  ,JLL    ,SLL    ,XLL    ,N_MUL_MX,
     2             ITIED  ,NINT   ,NKMAX  ,NC     ,V      ,A      ,
     3             IADLL  ,III    ,IIIS   ,NI_M   ,NI_S   ,
     4             NX     ,NY     ,NZ     ,VIT    ,COMNTAG,
     5             ICONT  ,RM     ,TM     ,RS     ,TS     )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
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   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT,NC,N_MUL_MX,ITIED,NINT ,NKMAX   
      INTEGER LLL(*),JLL(*),SLL(*),COMNTAG(*),
     .        iii(mvsiz,17),iadll(*) ,iiis(mvsiz,16),
     +        icont(mvsiz)
C     REAL
      my_real
     .  xll(*),v(3,*),a(3,*),vit(*)
      my_real
     .   rm(mvsiz) ,rs(mvsiz) ,tm(mvsiz) ,ts(mvsiz) ,
     .   ni_m(mvsiz,*) ,ni_s(mvsiz,*),nx(mvsiz) ,ny(mvsiz) ,nz(mvsiz)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,IK,NK,I1,I2,I3,I4,IAD,NN
      my_real
     .   vx,vy,vz,vn,aa
C-----------------------------------------------
      IF(itied==0)THEN
       DO i=1,llt
C-----------------------------------------------
C       Test if contact
C-----------------------------------------------
        IF(icont(i)/=0)THEN
C
           nk = 16
            nc=nc+1
            IF(nc>n_mul_mx)THEN
              CALL ancmsg(msgid=89,anmode=aninfo)
              CALL arret(2)
            ENDIF
            iadll(nc+1)=iadll(nc) + 48
            IF(iadll(nc+1)-1>nkmax)THEN
              CALL ancmsg(msgid=89,anmode=aninfo)
              CALL arret(2)
            ENDIF
            iad = iadll(nc) - 1
            DO ik=1,8
                lll(iad+ik)    = iii(i,ik)
                jll(iad+ik)    = 1
                sll(iad+ik)    = 0
                xll(iad+ik)    = nx(i)*ni_m(i,ik)
                lll(iad+ik+16) = iii(i,ik)
                jll(iad+ik+16) = 2
                sll(iad+ik+16) = 0
                xll(iad+ik+16) = ny(i)*ni_m(i,ik)
                lll(iad+ik+32) = iii(i,ik)
                jll(iad+ik+32) = 3
                sll(iad+ik+32) = 0
                xll(iad+ik+32) = nz(i)*ni_m(i,ik)
C
                lll(iad+ik+8)  = iiis(i,ik)
                jll(iad+ik+8)  = 1
                sll(iad+ik+8)  = nint
                xll(iad+ik+8)  = -nx(i)*ni_s(i,ik)
                lll(iad+ik+24) = iiis(i,ik)
                jll(iad+ik+24) = 2
                sll(iad+ik+24) = nint
                xll(iad+ik+24) = -ny(i)*ni_s(i,ik)
                lll(iad+ik+40) = iiis(i,ik)
                jll(iad+ik+40) = 3
                sll(iad+ik+40) = nint
                xll(iad+ik+40) = -nz(i)*ni_s(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
c#include "lockoff.inc"
        ENDIF
       ENDDO
      ELSEIF(itied==1)THEN
C-----------------------------------------------
C      ITIED = 1
C-----------------------------------------------
       DO i=1,llt
C-----------------------------------------------
C       Test if contact
C-----------------------------------------------
        IF(icont(i)/=0)THEN
C
           nk = 16
            IF(nc+3>n_mul_mx)THEN
              CALL ancmsg(msgid=89,anmode=aninfo)
              CALL arret(2)
            ENDIF
            IF(iadll(nc+1)-1+16*3>nkmax)THEN
              CALL ancmsg(msgid=89,anmode=aninfo)
              CALL arret(2)
            ENDIF
C
            nc=nc+1
            iadll(nc+1)=iadll(nc) + 16
            iad = iadll(nc) - 1
            DO ik=1,8
                lll(iad+ik)   = iii(i,ik)
                jll(iad+ik)   = 1
                sll(iad+ik)   = 0
                xll(iad+ik)   = ni_m(i,ik)
                lll(iad+ik+8) = iiis(i,ik)
                jll(iad+ik+8) = 1
                sll(iad+ik+8) = nint
                xll(iad+ik+8) = -ni_s(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
C
            nc=nc+1
            iadll(nc+1)=iadll(nc) + 16
            iad = iadll(nc) - 1
            DO ik=1,8
                lll(iad+ik)   = iii(i,ik)
                jll(iad+ik)   = 2
                sll(iad+ik)   = 0
                xll(iad+ik)   = ni_m(i,ik)
                lll(iad+ik+8) = iiis(i,ik)
                jll(iad+ik+8) = 2
                sll(iad+ik+8) = nint
                xll(iad+ik+8) = -ni_s(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
            ENDDO
C
            nc=nc+1
            iadll(nc+1)=iadll(nc) + 16
            iad = iadll(nc) - 1
            DO ik=1,8
                lll(iad+ik)   = iii(i,ik)
                jll(iad+ik)   = 3
                sll(iad+ik)   = 0
                xll(iad+ik)   = ni_m(i,ik)
                lll(iad+ik+8) = iiis(i,ik)
                jll(iad+ik+8) = 3
                sll(iad+ik+8) = nint
                xll(iad+ik+8) = -ni_s(i,ik)
                nn = lll(iad+ik)
                comntag(nn) = comntag(nn) + 1
           ENDDO
c#include "lockoff.inc"
        ENDIF
       ENDDO
      ENDIF
c
c      print *, "r,s,t",r(1),s(1),t(1)
C
      RETURN
      END
C
!||====================================================================
!||    i17rst        ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17lll        ../engine/source/interfaces/int17/i17lagm.F
!||    i17lll_pena   ../engine/source/interfaces/int17/i17for3.F
!||--- calls      -----------------------------------------------------
!||    i16deri       ../engine/source/interfaces/int16/i16lagm.F
!||    i16ni         ../engine/source/interfaces/int16/i16lagm.F
!||    i16rstn       ../engine/source/interfaces/int16/i16lagm.F
!||====================================================================
      SUBROUTINE i17rst(LLT,R  ,S     ,T     ,NI    ,
     2            XX    ,YY    ,ZZ    )
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT
C     REAL
      my_real
     .   xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17)
      my_real
     .   r(mvsiz),s(mvsiz),t(mvsiz),ni(mvsiz,17) 
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,IK,NK,ITER,NITERMAX,JTER,NJTERMAX,CONV
      my_real
     .   VX,VY,VZ,VN
      my_real
     .   DRDX(MVSIZ),DRDY(MVSIZ),DRDZ(MVSIZ),
     .   DSDX(MVSIZ),DSDY(MVSIZ),DSDZ(MVSIZ),     
     .   DTDX(MVSIZ),DTDY(MVSIZ),DTDZ(MVSIZ),
     .   dxdr(mvsiz),dydr(mvsiz),dzdr(mvsiz),
     .   dxdt(mvsiz),dydt(mvsiz),dzdt(mvsiz),
     .   rr(mvsiz),ss(mvsiz),tt(mvsiz)     
C-----------------------------------------------
C
C      r=s=t=0
C
C  +---> iter 
C  |                  
C  |     Ni(r,s,t) =  
C  |     dNi/dr    =  
C  |     ...         _ 
C  |                \
C  |     dx/dr    = /_ (xi * dNi/dr)
C  |     ...        
C  | 
C  |            [dx/dr dy/dr dz/dr]            
C  |      [J] = |dx/ds dy/ds dz/ds|                
C  |            [dx/dt dy/dt dz/dt]               
C  | 
C  | +--> jter                                 
C  | |                  _ 
C  | |                 \
C  | |      x(r,s,t) = /_ (xi * Ni(r,s,t))
C  | |      ...        
C  | |
C  | |      |r|   |r|      -1  |xs-x(r,s,t)|
C  | |      {s} = {s} + [J]    {ys-y(r,s,t)}
C  | |      |t|   |t|          |zs-z(r,s,t)|
C  | |                   
C  | |      Ni(r,s,t) =  
C  +-+---
C-----------------------------------------------
       nitermax = 3
       njtermax = 3
       conv = 0
C
       DO i=1,llt
         rr(i) = 0.
         ss(i) = 0.
         tt(i) = 0.
       ENDDO
C-----------------------------------------------
C      calculation de r,s,t   et   Ni(r,s,t)
C-----------------------------------------------
       DO iter=1,nitermax
c
c      print *, "iter",iter
c
C-----------------------------------------------
C          calculation de Ni(r,s,t); [J]; [J]-1
C-----------------------------------------------
           CALL i16deri(llt,rr ,ss    ,tt    ,ni    ,
     2            drdx  ,drdy  ,drdz  ,dsdx  ,dsdy  ,dsdz  ,
     3            dtdx  ,dtdy  ,dtdz  ,dxdr  ,dydr  ,dzdr  ,
     4            dxdt  ,dydt  ,dzdt  ,xx    ,yy    ,zz    )
C
           DO jter=1,njtermax
c
c      print *, "jter",jter
c
C-----------------------------------------------
C            calculation de r,s,t new
C-----------------------------------------------
             CALL i16rstn(llt,rr,ss   ,tt    ,ni    ,conv  ,
     2            drdx  ,drdy  ,drdz  ,dsdx  ,dsdy  ,dsdz  ,
     3            dtdx  ,dtdy  ,dtdz  ,xx    ,yy    ,zz    ,
     4            r     ,s     ,t     )
c
c      print *, "r,s,t",r(1),s(1),t(1)
c      print *, "rr,ss,tt",rr(1),ss(1),tt(1)
c
C-----------------------------------------------
C            calculation de Ni(-1<r<1 , -1<s<1 , -1<t<1)
C-----------------------------------------------
             CALL i16ni(llt,rr ,ss ,tt ,ni  )
C
           ENDDO
      ENDDO
C
      RETURN
      END
C
!||====================================================================
!||    i17mini     ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17lll      ../engine/source/interfaces/int17/i17lagm.F
!||--- calls      -----------------------------------------------------
!||    i17abc      ../engine/source/interfaces/int17/i17lagm.F
!||    i17borne    ../engine/source/interfaces/int17/i17lagm.F
!||    i17norm     ../engine/source/interfaces/int17/i17lagm.F
!||    i17racine   ../engine/source/interfaces/int17/i17lagm.F
!||====================================================================
      SUBROUTINE i17mini(LLT   ,R_CS  ,S_CS  ,T_CS  ,RI_S  ,SI_S  ,
     2                   TI_S  ,NI_S  ,XXS   ,YYS   ,ZZS   ,XX    ,
     3                   YY    ,ZZ    ,R_CM  ,S_CM  ,T_CM  ,NX    ,
     4                   NY    ,NZ    ,R_1S  ,R_2S  ,T_1S  ,T_2S  ,
     5                   R_1M  ,R_2M  ,R_3M  ,R_4M  ,T_1M  ,T_2M  ,
     6                   T_3M  ,T_4M  ,ICONT )
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT,
     +        ICONT(MVSIZ,4)
C     REAL
      my_real
     + si_s(mvsiz,*),ni_s(mvsiz,*),ri_s,ti_s ,
     + xxs(mvsiz,*) ,yys(mvsiz,*) ,zzs(mvsiz,*) ,
     + xx(mvsiz,*)  ,yy(mvsiz,*)  ,zz(mvsiz,*) ,
     + r_cm(mvsiz) ,s_cm(mvsiz) ,t_cm(mvsiz),
     + r_cs(mvsiz) ,s_cs(mvsiz) ,t_cs(mvsiz),
     + r_1s(mvsiz) ,r_2s(mvsiz) ,t_1s(mvsiz) ,t_2s(mvsiz),
     + r_1m(mvsiz) ,r_2m(mvsiz) ,t_1m(mvsiz) ,t_2m(mvsiz),
     + r_3m(mvsiz) ,r_4m(mvsiz) ,t_3m(mvsiz) ,t_4m(mvsiz),
     + nx(mvsiz) ,ny(mvsiz) ,nz(mvsiz)   
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,ITER,NITERMAX
      my_real
     +        A1(MVSIZ),A2(MVSIZ),A3(MVSIZ),A4(MVSIZ),A5(MVSIZ),
     +        B1(MVSIZ),B2(MVSIZ),B3(MVSIZ),B4(MVSIZ),B5(MVSIZ),
     +        C1(MVSIZ),C2(MVSIZ),C3(MVSIZ),
     +        f1,f2,f3,f4,f5,f6,f7,f8,
     +        cc1,cc2,cc3,dd1,dd2,dd3,dd,d,
     +        a0,ab,ba,a4r,b4t,a5t,b5r,eps
C-----------------------------------------------
C ro = r ri    to = t ti
C
C i=1,4  
C    ri=+-1  ti=+-1
C    Ni = 1/4 (1+ro)(1+to)(ro+to-1)
C    Ni = 1/4 (r^2+r^2to+t^2+roto+rot^2-1)
C    Ni = 1/4 (r^2(1+to)+t^2(1+ro)+roto-1)
C    dNi/dr = ri/4 (1+to)(2ro+to)
C    dNi/dt = ti/4 (1+ro)(2to+ro)
C    dNi/dr = 1/4 (2 r + ri ti t + 2ti rt + ri t^2)
C    dNi/dt = 1/4 (2 t + ri ti r + 2ri rt + ti r^2)
C
C i=6;8
C    ri=0    ti=+-1
C    Ni = 1/2 (1-r^2)(1+to)
C    dNi/dr = -r (1+to)
C    dNi/dt = ti/2 (1-r^2)
C    dNi/dr = 1/4 (-4 r - 4ti r t)
C    dNi/dt = 1/4 (2ti - 2ti r^2 )
C
C
C i=5;7
C    ri=+-1  ti=0
C    Ni = 1/2 (1-t^2)(1+ro)
C    dNi/dr = ri/2 (1-t^2)
C    dNi/dt = -t (1+ro)
C    dNi/dr = 1/4 (2ri  - 2ri t^2 )
C    dNi/dt = 1/4 (-4 t - 4ri rt  )
C-----------------------------------------------
C  df/dr = Somme( fi dNi/dr )
C  df/dt = Somme( fi dNi/dt )
C
C  df/dr = A1 + A2 r + A3 t + A4 rt + A5 t^2
C  df/dt = B1 + B2 r + B3 t + B4 rt + B5 r^2
C 
C  A1 = (                      -f5          +   f7        )/2
C  A2 = ( f1 + f2 + f3 + f4         - 2 f6         - 2 f8)/2
C  A3 = ( f1 - f2 + f3 - f4                                )/4
C  A4 = (-f1 + f2 + f3 - f4         - 2 f6         + 2 f8)/2
C  A5 = (-f1 - f2 + f3 + f4 + 2 f5          - 2 f7        )/4
C
C  B1 = (                               f6         -   f8)/2
C  B2 = ( f1 - f2 + f3 - f4                                )/4
C  B3 = ( f1 + f2 + f3 + f4 - 2 f5          - 2 f7        )/2
C  B4 = (-f1 - f2 + f3 + f4 + 2 f5          - 2 f7        )/2
C  B5 = (-f1 + f2 + f3 - f4         - 2 f6         + 2 f8)/4
C
C  df/dr = A1 + A2 r + A3 t + A4 rt + A5 t^2 = 0
C  df/dt = B1 + B2 r + B3 t + B4 rt + B5 r^2 = 0
C  r = -(A1 + A3 t + A5 t^2 ) / (A2 + A4 t)
C  t = -(B1 + B2 r + B5 r^2 ) / (B3 + B4 r)
c
c
c                                r
C                                ^
c                                |
C                          .     |7    
C            4 O-----------------O-----------------O 3        
C              |           .     |                 |          
C              |           .     |                 |          
C              |                 |      C          |
C              |           .    r+------+df/dt=0   |  
C              |           .     |      |df/dr=0   |    
C              |                 |      |          |
C              |                 |      |          |   
C              |                 |      |          |6       
C            8 0                 +------+----------0----> t  
C              |                        t          |       
C              |                                   |      
C              |                                   |        
C              |                                   |      
C              |                                   |    
C              |                                   |    
C              |                                   |    
C              |'                                  |    
C              O-----------------O-----------------O    
C             1                  5                  2
C
C-----------------------------------------------
C
c
c
      nitermax = 5
C
      DO i=1,llt
        d = si_s(i,1)*si_s(i,1)+si_s(i,2)*si_s(i,2)
     +    + si_s(i,3)*si_s(i,3)+si_s(i,4)*si_s(i,4)
     +    + si_s(i,5)*si_s(i,5)+si_s(i,6)*si_s(i,6)
     +    + si_s(i,7)*si_s(i,7)+si_s(i,8)*si_s(i,8)
        d = 1./max(em20,sqrt(d))
        f1 = d * si_s(i,1)
        f2 = d * si_s(i,2)
        f3 = d * si_s(i,3)
        f4 = d * si_s(i,4)
        f5 = d * si_s(i,5)
        f6 = d * si_s(i,6)
        f7 = d * si_s(i,7)
        f8 = d * si_s(i,8)
        a0 = ( f1 + f2 + f3 + f4 )*half
        ab = f5 - f7
        ba = f6 - f8
C
        a1(i) = -ab * half
        a2(i) = a0 - f6 - f8
        a3(i) = ( f1 - f2 + f3 - f4 )*fourth
        a4(i) = (-f1 + f2 + f3 - f4 )*half  - ba
        a5(i) = (-f1 - f2 + f3 + f4 )*fourth - a1(i)
C
        b1(i) = ba*half
        b2(i) = a3(i)
        b3(i) = a0 - f5  - f7
        b4(i) = (-f1 - f2 + f3 + f4 )*half  + ab
        b5(i) = (-f1 + f2 + f3 - f4 )*fourth - b1(i)
c
        r_cs(i) = zero
        t_cs(i) = zero
      ENDDO
c------------------------------------------------
c   Newton ITER: lineari_sation in r and t
c------------------------------------------------
C  fr = df/dr = A1 + A2 r + A3 t + A4 t r + A5 t t = 0
C  ft = df/dt = B1 + B2 r + B3 t + B4 r t + B5 r r = 0 
c
C  fr = fr_ + dfr/dr dr + dfr/dt dt
C  ft = ft_ + dft/dr dr + dft/dt dt
c
C  fr = fr_ + dfr/dr (r-r_) + dfr/dt (t-t_) = 0
C  ft = ft_ + dft/dr (r-r_) + dft/dt (t-t_) = 0
c
C  dfr/dr = A2 + A4 t_ = C2
C  dfr/dt = A3 + A4 r_ + 2 A5 t_ = C3
C  dft/dr = B2 + B4 t_ + 2 B5 r_ = D2
C  dft/dt = B3 + B4 r_ = D3
c
c  C1 = A1 - A4 r_ t_ - A5 t_^2
c  C2 = A2 + A4 t_
c  C3 = A3 + A4 r_ + 2 A5 t_
c
c  D1 = B1 - B4 r_ t_ - B5 r_^2
c  D2 = B2 + B4 t_ + 2 B5 r_
c  D3 = B3 + B4 r_ 
c
C  fr = C1 + C2 r + C3 t = 0
C  ft = D1 + D2 r + D3 t = 0
c
C  r = (C3 D1 - D3 C1) / (D3 C2 - C3 D2) 
C  t = (D2 C1 - C2 D1) / (D3 C2 - C3 D2)
c------------------------------------------------
      DO iter=1,nitermax
        DO i=1,llt
          a4r = a4(i) * r_cs(i)
          a5t = a5(i) * t_cs(i)
          b4t = b4(i) * t_cs(i)
          b5r = b5(i) * r_cs(i)
          cc1  = a1(i) -(a4r + a5t) * t_cs(i)
          cc2  = a2(i) + a4(i) * t_cs(i)
          cc3  = a3(i) + a4r + a5t + a5t
          dd1  = b1(i) -(b4t + b5r )* r_cs(i)
          dd2  = b2(i) + b4t + b5r + b5r
          dd3  = b3(i) + b4(i) * r_cs(i)
          d    = dd3 * cc2 - cc3 * dd2 
          IF(abs(d)<em20)THEN
            cc2 = cc2 + em10
            dd3 = dd3 + em10
            d   = dd3 * cc2 - cc3 * dd2 
          ENDIF
          r_cs(i) = (cc3 * dd1 - dd3 * cc1) / d
          t_cs(i) = (dd2 * cc1 - cc2 * dd1) / d
          r_cs(i) = max(-one,min(one,r_cs(i)))
          t_cs(i) = max(-one,min(one,t_cs(i)))
c          dfdr = A1 + A2 * r + A3 * t + A4 * t * r + A5 * t * t 
c          dfdt = B1 + B2 * r + B3 * t + B4 * r * t + B5 * r * r 
c          print *,'it=',iter,' r=',r,' t=',t,' fr=',dfdr,' ft=',dfdt
        ENDDO
      ENDDO
C-----------------------------------------------------------------------
C     calculation of 4 points of the iso S_M = +- 1
C 
C          second.           main
C     1: <R_1S , T_CS>    <R_1M , T_1M>
C     2: <R_2S , T_CS>    <R_2M , T_2M>
C     3: <R_CS , T_1S>    <R_3M , T_3M>
C     4: <R_CS , T_2S>    <R_4M , T_4M> 
C
C     We limit R and you +- 1 on the second element.::
C                -1 < R_1S < 1  et -1 < R_2S < 1 
C                -1 < T_1S < 1  et -1 < T_2S < 1 
c
C     We limit R and you +- 1 on the main element:
C                -1 < R_M < 1   et -1 < T_M < 1 
C-----------------------------------------------------------------------
      CALL  i17abc(llt ,si_s,r_cs,t_cs,
     +             b1  ,b2  ,b3  ,c1  ,c2  ,c3  )
C
      DO i=1,llt
        s_cm(i) = b1(i) + b2(i)*r_cs(i) + b3(i)*r_cs(i)*r_cs(i)
      ENDDO
c
      DO i=1,llt
          IF(s_cm(i)>zero)THEN
            b1(i) = b1(i) - one
            c1(i) = c1(i) - one
          ELSE
            b1(i) = b1(i) + one
            c1(i) = c1(i) + one
          ENDIF
      ENDDO
c
C  f = B1 + B2 r + B3 r^2
c
      CALL i17racine(llt,b3,b2,b1,r_1s,r_2s)
 
C  f = C1 + C2 t + C3 t^2
C
      CALL i17racine(llt,c3,c2,c1,t_1s,t_2s)
c-----------------------------------------------------------------------
c      bound r at +-1 on the main
c-----------------------------------------------------------------------
C
      CALL  i17abc(llt ,ri_s,r_cs,t_cs,
     +             b1  ,b2  ,b3  ,c1  ,c2  ,c3  )
C
      DO i=1,llt
        r_cm(i) = b1(i) + (b2(i) + b3(i)*r_cs(i))*r_cs(i)
        icont(i,1) = 1
        icont(i,2) = 1
        icont(i,3) = 1
        icont(i,4) = 1
      ENDDO
C
c      CALL I17BORNE(LLT ,R_CS,B3  ,B2  ,B1  ,ICONT(1,1),R_CM)
c      CALL I17BORNE(LLT ,T_CS,C3  ,C2  ,C1  ,ICONT(1,1),R_CM)
      CALL i17borne(llt ,r_1s,b3  ,b2  ,b1  ,icont(1,1),r_1m)
      CALL i17borne(llt ,r_2s,b3  ,b2  ,b1  ,icont(1,2),r_2m)
      CALL i17borne(llt ,t_1s,c3  ,c2  ,c1  ,icont(1,3),r_3m)
      CALL i17borne(llt ,t_2s,c3  ,c2  ,c1  ,icont(1,4),r_4m)
c-----------------------------------------------------------------------
c  2: bound t at +-1 on the main
c-----------------------------------------------------------------------
      CALL  i17abc(llt ,ti_s,r_cs,t_cs,
     +             b1  ,b2  ,b3  ,c1  ,c2  ,c3  )
C
      DO i=1,llt
        t_cm(i) = b1(i) + (b2(i) + b3(i)*r_cs(i))*r_cs(i)
      ENDDO
C
c      CALL I17BORNE(LLT ,R_CS,B3  ,B2  ,B1  ,ICONT(1,1),T_CM)
c      CALL I17BORNE(LLT ,T_CS,C3  ,C2  ,C1  ,ICONT(1,1),T_CM)
      CALL i17borne(llt ,r_1s,b3  ,b2  ,b1  ,icont(1,1),t_1m)
      CALL i17borne(llt ,r_2s,b3  ,b2  ,b1  ,icont(1,2),t_2m)
      CALL i17borne(llt ,t_1s,c3  ,c2  ,c1  ,icont(1,3),t_3m)
      CALL i17borne(llt ,t_2s,c3  ,c2  ,c1  ,icont(1,4),t_4m)
c-----------------------------------------------------------------------
c  Test if contact
c-----------------------------------------------------------------------
      CALL  i17abc(llt ,si_s,r_cs,t_cs,
     +             b1  ,b2  ,b3  ,c1  ,c2  ,c3  )
      eps = em3
      DO i=1,llt
        d = b1(i) + (b2(i) + b3(i)*r_1s(i)) * r_1s(i)
        IF(d<-one-eps.OR.d>one+eps) icont(i,1)=0
        d = b1(i) + (b2(i) + b3(i)*r_2s(i)) * r_2s(i)
        IF(d<-one-eps.OR.d>one+eps) icont(i,2)=0
        d = c1(i) + (c2(i) + c3(i)*t_1s(i)) * t_1s(i)
        IF(d<-one-eps.OR.d>one+eps) icont(i,3)=0
        d = c1(i) + (c2(i) + c3(i)*t_2s(i)) * t_2s(i)
        IF(d<-one-eps.OR.d>one+eps) icont(i,4)=0
c        R_CS(I) = ICONT(I,1)*R_1S(I) + ICONT(I,2)*R_2S(I) +  
c     +            ICONT(I,3)*R_CS(I) + ICONT(I,4)*R_CS(I)
c        T_CS(I) = ICONT(I,1)*T_CS(I) + ICONT(I,2)*T_CS(I) +  
c     +            ICONT(I,3)*T_1S(I) + ICONT(I,4)*T_2S(I)
cc        R_CM(I) = ICONT(I,1)*R_1M(I) + ICONT(I,2)*R_2M(I) +  
cc     +            ICONT(I,3)*R_3M(I) + ICONT(I,4)*R_4M(I)
cc        T_CM(I) = ICONT(I,1)*T_1M(I) + ICONT(I,2)*T_2M(I) +  
cc     +            ICONT(I,3)*T_3S(I) + ICONT(I,4)*T_4S(I)
c        ICONT(I,1) = ICONT(I,1) + ICONT(I,2) +  
c     +               ICONT(I,3) + ICONT(I,4)
c        IF(ICONT(I,1)/=0)THEN
c          R_CS(I) = R_CS(I) / ICONT(I,1) 
c          T_CS(I) = T_CS(I) / ICONT(I,1) 
c          ICONT(I,1) = 1 
c        ENDIF
      ENDDO
C
      CALL  i17abc(llt ,si_s,r_cs,t_cs,
     +             b1  ,b2  ,b3  ,c1  ,c2  ,c3  )
      DO i=1,llt
        s_cm(i) = b1(i) + (b2(i) + b3(i)*r_cs(i))*r_cs(i)
      ENDDO
C
      CALL  i17abc(llt ,ri_s,r_cs,t_cs,
     +             b1  ,b2  ,b3  ,c1  ,c2  ,c3  )
      DO i=1,llt
        r_cm(i) = b1(i) + (b2(i) + b3(i)*r_cs(i))*r_cs(i)
      ENDDO
C
      CALL  i17abc(llt ,ti_s,r_cs,t_cs,
     +             b1  ,b2  ,b3  ,c1  ,c2  ,c3  )
      DO i=1,llt
        t_cm(i) = b1(i) + (b2(i) + b3(i)*r_cs(i))*r_cs(i)
      ENDDO
C-----------------------------------------------
C          calculation de Ni(r,s,t); dNi/dr; dNi/dt; normale second(-n)
C-----------------------------------------------
      CALL i17norm(llt  ,r_cs  ,s_cs  ,t_cs  ,
     2             nx   ,ny    ,nz    ,xxs   ,yys   ,zzs   )
c-----------------------------------------------------------------------
c  
c-----------------------------------------------------------------------
      RETURN
      END
!||====================================================================
!||    i17ni         ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17vit4       ../engine/source/interfaces/int17/i17lagm.F
!||    i17vit_pena   ../engine/source/interfaces/int17/i17for3.F
!||====================================================================
      SUBROUTINE i17ni(LLT,R ,T ,NI    )
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT
      my_real
     .   r(mvsiz),t(mvsiz),ni(mvsiz,*)  
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  U_M_R,U_P_R,U_M_S,U_P_S,U_M_T,U_P_T,
     .  ums_umt,ums_upt,ups_umt,ups_upt,
     .  umr_ums,umr_ups,upr_ums,upr_ups,
     .  umt_umr,umt_upr,upt_umr,upt_upr,
     .  a,r05,s05,t05
C-----------------------------------------------------------------------
C     calculation of Ni
C-----------------------------------------------------------------------
      DO i=1,llt
C
C
        r05 = half*r(i)
        t05 = half*t(i)
C
        u_m_r = half - r05
        u_p_r = half + r05
C
        u_m_t = half - t05
        u_p_t = half + t05
C
        ni(i,1) = u_m_t * u_m_r * (-r(i)-t(i)-one)
        ni(i,2) = u_p_t * u_m_r * (-r(i)+t(i)-one)
        ni(i,3) = u_p_t * u_p_r * ( r(i)+t(i)-one)
        ni(i,4) = u_m_t * u_p_r * ( r(i)-t(i)-one)
C------------------------------------
        a = (one-r(i)*r(i))
        ni(i,6) = a * u_p_t
        ni(i,8) = a * u_m_t
C------------------------------------
        a = (one-t(i)*t(i))
        ni(i,5) = a * u_m_r
        ni(i,7) = a * u_p_r
C------------------------------------
      ENDDO
C-----------------------------------------------
      RETURN
      END
!||====================================================================
!||    i17racine       ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17borne        ../engine/source/interfaces/int17/i17lagm.F
!||    i17borne_pena   ../engine/source/interfaces/int17/i17for3.F
!||    i17mini         ../engine/source/interfaces/int17/i17lagm.F
!||====================================================================
      SUBROUTINE i17racine(LLT,A,B,C,R1,R2)
C-----------------------------------------------
C
C    calculation des racines r1,r2 bornees a +- 1
C
C    de  a x^2 + b x + c = 0
C 
C    The routine returns -1,+1 if there are no roots
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT
      my_real
     .   c(mvsiz),b(mvsiz),a(mvsiz),r1(mvsiz),r2(mvsiz) 
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  D      
c
C
      DO I=1,llt
        d =  b(i)*b(i) - four*a(i)*c(i) 
        IF(abs(a(i))<em20)THEN
          IF(abs(b(i))<em20)THEN
            r1(i) = -one
            r2(i) =  one
          ELSE
            r1(i) = -c(i) / b(i)
            r2(i) = b(i) / abs(b(i))
          ENDIF
        ELSEIF(d<zero)THEN
            r1(i) = -one
            r2(i) =  one      
        ELSE
            d = sqrt( d )
            r1(i) = (-b(i) - d) / (two*a(i))
            r2(i) = (-b(i) + d) / (two*a(i))
        ENDIF
      ENDDO
c-----------------------------------------------------------------------
c      bound r (or t) at +-1 on the second
c-----------------------------------------------------------------------
      DO i=1,llt
          r1(i) = max(-one,min(one,r1(i)))
          r2(i) = max(-one,min(one,r2(i)))
      ENDDO
c
      RETURN
      END
!||====================================================================
!||    i17borne    ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17mini     ../engine/source/interfaces/int17/i17lagm.f
!||--- calls      -----------------------------------------------------
!||    i17racine   ../engine/source/interfaces/int17/i17lagm.f
!||====================================================================
      SUBROUTINE i17borne(LLT ,R_S ,A   ,B  ,C  ,ICONT,RS)
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT,ICONT(MVSIZ)
      my_real
     +        r_s(mvsiz),c(mvsiz),b(mvsiz),a(mvsiz),rs(mvsiz)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     +   CC(MVSIZ),R1(MVSIZ),R2(MVSIZ),RS1,RS2
c-----------------------------------------------------------------------
c  1: bound r (or s on the 2nd call) at +-1 on the main
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c  1.1: bound <r1,t> and <r,t1> at +-1 on the main
c-----------------------------------------------------------------------
        DO i=1,llt
          rs(i) = c(i) + ( b(i) + a(i)*r_s(i) )*r_s(i)
          IF(rs(i)>one)THEN
            cc(i) = c(i) - one
          ELSEIF(rs(i)<-one)THEN
            cc(i) = c(i) + one
          ELSE 
            cc(i) = one
          ENDIF
        ENDDO
c
C  f = C + B r + A r^2
c
        CALL i17racine(llt,a,b,cc,r1,r2)
c
        DO i=1,llt
          IF(rs(i)>one.OR.rs(i)<-one)THEN
            rs1 = c(i) + ( b(i) + a(i)*r1(i) )*r1(i)
            rs2 = c(i) + ( b(i) + a(i)*r2(i) )*r2(i)
            IF(rs1>=-one.AND.rs1<=one.AND.
     +         rs2>=-one.AND.rs2<=one)THEN
              IF(abs(rs(i)-r1(i))<abs(rs(i)-r2(i)))THEN
                r_s(i) = r1(i)
                rs(i)  = rs1
              ELSE
                r_s(i) = r2(i)
                rs(i)  = rs2
              ENDIF
            ELSEIF(rs1>=-one.AND.rs1<=one)THEN
                r_s(i) = r1(i)
                rs(i)  = rs1
            ELSEIF(rs2>=-one.AND.rs2<=one)THEN
                r_s(i) = r2(i)
                rs(i)  = rs2
            ELSE
C No main/second cover
                icont(i)=0
            ENDIF
          ENDIF
        ENDDO
c
      RETURN
      END
!||====================================================================
!||    i17abc         ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17mini        ../engine/source/interfaces/int17/i17lagm.F
!||    i17mini_pena   ../engine/source/interfaces/int17/i17for3.F
!||====================================================================
      SUBROUTINE i17abc(LLT ,F   ,R   ,T   ,
     +                  B1  ,B2  ,B3  ,C1  ,C2  ,C3  )
C-----------------------------------------------
C i=1,4  
C    ri=+-1  ti=+-1
C    Ni = 1/4 (1+ro)(1+to)(ro+to-1)
C    Ni = 1/4 (r^2 + ti r^2 t + t^2 + ri ti r t + ri r t^2 - 1)
C
C i=6;8
C    ri=0    ti=+-1
C    Ni = 1/2 ( 1 + ti t - r^2 - ti t r^2)
C
C
C i=5;7
C    ri=+-1  ti=0
C    Ni = 1/2 (1 + ri r - t^2 - ri r t^2)
C-----------------------------------------------
C  f = Somme( fi Ni )
C
C  f = A1 + A2 r + A3 t + A4 rt + A5 r^2 + A6 t^2 + A7 tr^2 + A8 rt^2
C  f = (A1 + A3 t + A6 t^2) + (A2 + A4 t + A8 t^2)r + (A5 + A7 t)r^2
C  f = B1 + B2 r + B3 r^2
C  f = (A1 + A2 r + A5 r^2) + (A3 + A4 r + A7 r^2)t + (A6 + A8 r)t^2
C  f = C1 + C2 t + C3 t^2
C 
C  A1 = (-f1 - f2 - f3 - f4 + 2 f5  + 2 f6  + 2 f7 + 2 f8)/4
C  A2 = (                   -   f5          +   f7       )/2
C  A3 = (                           +   f6         -   f8)/2
C  A4 = (+f1 - f2 + f3 - f4                              )/4
C  A5 = (+f1 + f2 + f3 + f4         - 2 f6         - 2 f8)/4                     
C  A6 = (+f1 + f2 + f3 + f4 - 2 f5         - 2 f7        )/4   
C  A7 = (-f1 + f2 + f3 - f4         - 2 f6         + 2 f8)/4   
C  A8 = (-f1 - f2 + f3 + f4 + 2 f5         - 2 f7        )/4  
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT
      my_real
     +        F(MVSIZ,*),R(MVSIZ),T(MVSIZ),
     +        B1(MVSIZ),B2(MVSIZ),B3(MVSIZ),
     +        C1(MVSIZ),C2(MVSIZ),C3(MVSIZ)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     +        A1,A2,A3,A4,A5,A6,A7,A8,R2,T2,FF5,FF6,FF7,FF8
C
      DO I=1,llt
        ff5 = f(i,5) + f(i,5)
        ff6 = f(i,6) + f(i,6)
        ff7 = f(i,7) + f(i,7)
        ff8 = f(i,8) + f(i,8)
c
c        A1 = (-F(I,1)-F(I,2)-F(I,3)-F(I,4)+FF5+FF6+FF7+FF8)*0.25 ...
c        B1(I) = A1 + ( A3 + A6 * T(I) ) *T(I) ...
c
        a1 = (-f(i,1)-f(i,2)-f(i,3)-f(i,4)+ff5+ff6+ff7+ff8)
        a2 = (                            -ff5    +ff7    )
        a3 = (                                +ff6    -ff8)
        a4 = (+f(i,1)-f(i,2)+f(i,3)-f(i,4)                )
        a5 = (+f(i,1)+f(i,2)+f(i,3)+f(i,4)    -ff6    -ff8)
        a6 = (+f(i,1)+f(i,2)+f(i,3)+f(i,4)-ff5    -ff7    ) 
        a7 = (-f(i,1)+f(i,2)+f(i,3)-f(i,4)    -ff6    +ff8) 
        a8 = (-f(i,1)-f(i,2)+f(i,3)+f(i,4)+ff5    -ff7    )
c
c
        b1(i) = ( a1 + ( a3 + a6 * t(i) ) *t(i) )*fourth
        b2(i) = ( a2 + ( a4 + a8 * t(i) ) *t(i) )*fourth
        b3(i) = ( a5 + a7 * t(i)                )*fourth
c
        c1(i) = ( a1 + ( a2 + a5 * r(i) ) *r(i) )*fourth
        c2(i) = ( a3 + ( a4 + a7 * r(i) ) *r(i) )*fourth
        c3(i) = ( a6 + a8 * r(i)                )*fourth
      ENDDO
c
      RETURN
      END
!||====================================================================
!||    i17norm        ../engine/source/interfaces/int17/i17lagm.F
!||--- called by ------------------------------------------------------
!||    i17mini        ../engine/source/interfaces/int17/i17lagm.F
!||    i17mini_pena   ../engine/source/interfaces/int17/i17for3.F
!||====================================================================
      SUBROUTINE i17norm(LLT ,RR  ,SS  ,TT    ,
     2                   NX  ,NY  ,NZ  ,XX    ,YY    ,ZZ    )
C-----------------------------------------------
C          calculation de Ni(r,t); dNi/dr; dNi/dt; normale
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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER LLT
C     REAL
      my_real
     .   xx(mvsiz,*) ,yy(mvsiz,*),zz(mvsiz,*),
     .   nx(mvsiz) ,ny(mvsiz) ,nz(mvsiz),
     .   rr(mvsiz) ,ss(mvsiz) ,tt(mvsiz)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,N
      my_real
     .   DXDR(MVSIZ), DYDR(MVSIZ), DZDR(MVSIZ),
     .   DXDT(MVSIZ), DYDT(MVSIZ), DZDT(MVSIZ),
     .   DNIDR(8),DNIDS(8),DNIDT(8)
      my_real
     .  U_M_R,U_P_R,U_M_S,U_P_S,U_M_T,U_P_T,
     .  UMS_UMT,UMS_UPT,UPS_UMT,UPS_UPT,
     .  UMR_UMS,UMR_UPS,UPR_UMS,UPR_UPS,
     .  UMT_UMR,UMT_UPR,UPT_UMR,UPT_UPR,
     .  A,R05,S05,T05,AA
C-----------------------------------------------------------------------
C     Ni; dNi/dr; dNi/ds; dNi/dt  s=+-1
C-----------------------------------------------------------------------
      DO i=1,llt
        r05 = half*rr(i)
        s05 = half*ss(i)
        t05 = half*tt(i)
C
        u_m_r = half - r05
        u_p_r = half + r05
C
        u_m_s = half - s05
        u_p_s = half + s05
C
        u_m_t = half - t05
        u_p_t = half + t05
C
        ums_umt = u_m_s * u_m_t
        ums_upt = u_m_s * u_p_t
        ups_umt = u_p_s * u_m_t
        ups_upt = u_p_s * u_p_t
C
        umr_ums = u_m_r * u_m_s
        umr_ups = u_m_r * u_p_s
        upr_ums = u_p_r * u_m_s
        upr_ups = u_p_r * u_p_s
C
        umt_umr = u_m_t * u_m_r
        umt_upr = u_m_t * u_p_r
        upt_umr = u_p_t * u_m_r
        upt_upr = u_p_t * u_p_r
C
        a = -t05 - rr(i)
        dnidr(1) = (-ums_umt - ups_umt) * a
        a =  t05 - rr(i)
        dnidr(2) = (-ums_upt - ups_upt) * a
        a =  t05 + rr(i)
        dnidr(3) =  (ums_upt + ups_upt) * a
        a = -t05 + rr(i)
        dnidr(4) =  (ums_umt + ups_umt) * a
C
        a = -r05 - tt(i)
        dnidt(1) = (-umr_ums - umr_ups) * a
        a = -r05 + tt(i)
        dnidt(2) =  (umr_ums + umr_ups) * a
        a = +r05 + tt(i)
        dnidt(3) =  (upr_ums + upr_ups) * a
        a = +r05 - tt(i)
        dnidt(4) = (-upr_ums - upr_ups) * a
C------------------------------------
       a = half*(one-rr(i)*rr(i))
        dnidt(6) =  a * (u_m_s + u_p_s)
C    
        a = -two*rr(i)
        dnidr(6) = a * (ums_upt + ups_upt)
        dnidr(8) = a * (ums_umt + ups_umt)
C------------------------------------
        a = half*(one-tt(i)*tt(i))
        dnidr(5)  = -a * (u_m_s + u_p_s)
C
        a = -two*tt(i)
        dnidt(5) = a * (umr_ums + umr_ups)
        dnidt(7) = a * (upr_ums + upr_ups)
C-----------------------------------------------------------------------
C     DX/DR;DX/DT
C-----------------------------------------------------------------------
        dxdr(i) = dnidr(1)*xx(i,1) + dnidr(2)*xx(i,2) + dnidr(3)*xx(i,3)
     +          + dnidr(4)*xx(i,4)   
     +    + dnidr(5)*(xx(i,5) - xx(i,7)) + dnidr(6)*xx(i,6)
     +    + dnidr(8)*xx(i,8) 
C  
        dxdt(i) = dnidt(1)*xx(i,1) + dnidt(2)*xx(i,2) + dnidt(3)*xx(i,3)
     +          + dnidt(4)*xx(i,4)   
     +    + dnidt(5)*xx(i,5)   + dnidt(6)*(xx(i,6) - xx(i,8))
     +    + dnidt(7)*xx(i,7)
C-----------------------------------------------------------------------
C     dy/dr; dy/dt
C-----------------------------------------------------------------------
        dydr(i) = dnidr(1)*yy(i,1) + dnidr(2)*yy(i,2) + dnidr(3)*yy(i,3)
     +          + dnidr(4)*yy(i,4)   
     +    + dnidr(5)*(yy(i,5) - yy(i,7))   + dnidr(6)*yy(i,6)
     +    + dnidr(8)*yy(i,8) 
C  
        dydt(i) = dnidt(1)*yy(i,1) + dnidt(2)*yy(i,2) + dnidt(3)*yy(i,3)
     +          + dnidt(4)*yy(i,4)   
     +    + dnidt(5)*yy(i,5)   + dnidt(6)*(yy(i,6) - yy(i,8))
     +    + dnidt(7)*yy(i,7)
C-----------------------------------------------------------------------
C     dz/dr; dz/dt
C-----------------------------------------------------------------------
        dzdr(i) = dnidr(1)*zz(i,1) + dnidr(2)*zz(i,2) + dnidr(3)*zz(i,3)
     +          + dnidr(4)*zz(i,4)   
     +    + dnidr(5)*(zz(i,5) - zz(i,7))   + dnidr(6)*zz(i,6)
     +    + dnidr(8)*zz(i,8) 
C  
        dzdt(i) = dnidt(1)*zz(i,1) + dnidt(2)*zz(i,2) + dnidt(3)*zz(i,3)
     +          + dnidt(4)*zz(i,4)   
     +    + dnidt(5)*zz(i,5)   + dnidt(6)*(zz(i,6) - zz(i,8))
     +    + dnidt(7)*zz(i,7)
C-----------------------------------------------------------------------
C     calculation of the normal -n
C-----------------------------------------------------------------------
        nx(i) = -dydt(i)*dzdr(i) + dzdt(i)*dydr(i)
        ny(i) = -dzdt(i)*dxdr(i) + dxdt(i)*dzdr(i)
        nz(i) = -dxdt(i)*dydr(i) + dydt(i)*dxdr(i)
C
        aa = one/sqrt(nx(i)*nx(i)+ny(i)*ny(i)+nz(i)*nz(i))
        nx(i) = nx(i)*aa
        ny(i) = ny(i)*aa
        nz(i) = nz(i)*aa
      ENDDO
c
      RETURN
      END