i15tott1.F File Reference

#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com04_c.inc"
#include "units_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	i15tott1 (noint, ndeb, ntc, x, ksurf, igrsurf, bufsf, ktc, ksi, nold, xp1, xp2, xp3, xtk, ytk, ztk, ntx, nty, ntz, penet, depth, xi, yi, zi, nxi, nyi, nzi, ansmx, hold, iactiv, itab)

Function/Subroutine Documentation

i15tott1()

subroutine i15tott1	(	integer	noint,
		integer	ndeb,
		integer	ntc,
			x,
		integer	ksurf,
		type (surf_), dimension(nsurf)	igrsurf,
			bufsf,
		integer, dimension(*)	ktc,
		integer, dimension(4,*)	ksi,
			nold,
			xp1,
			xp2,
			xp3,
			xtk,
			ytk,
			ztk,
			ntx,
			nty,
			ntz,
			penet,
			depth,
			xi,
			yi,
			zi,
			nxi,
			nyi,
			nzi,
			ansmx,
			hold,
		integer, dimension(4,*)	iactiv,
		integer, dimension(*)	itab )

Definition at line 30 of file i15tott1.F.

C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE groupdef_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   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
#include      "units_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER KSURF ,KSI(4,*) ,
     .        NOINT ,NDEB  ,NTC, IACTIV(4,*), KTC(*), ITAB(*)
C     REAL
      my_real 
     .  x(3,*)    ,bufsf(*) ,nold(3,*) ,
     .  xtk(*) ,ytk(*) ,ztk(*) ,ntx(*) ,nty(*) ,ntz(*) ,
     .  penet(*) ,depth(*) ,xi(*) ,yi(*)  ,zi(*)  ,
     .  nxi(*)   ,nyi(*) ,nzi(*) ,xp1(3,mvsiz), xp2(3,mvsiz),
     .  xp3(3,mvsiz), ansmx, hold(3,*)
      TYPE (SURF_)   , DIMENSION(NSURF)   :: IGRSURF
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER ADRBUF, I, IL, NLS, IDG,
     .        IN1, IN2, IN3,
     .        INSIDE1,INSIDE2
      my_real
     .   a, b, c, an, bn, cn, rot(9), dgr, expn,
     .   xg, yg, zg,
     .   ntn,
     .   n1, n2, n3, n, n11,n12,n13,nr1,n21,n22,n23,nr2,
     .   xkn1, ykn1, zkn1, sgnxkn, sgnykn, sgnzkn, xkn, ykn, zkn, eh,
     .   lambda1, lambda2, alp, bet,
     .   xh , yh , zh , mu, xh1, yh1, zh1, mu1, xh2, yh2, zh2, mu2,
     .   dx1, dy1, dz1, dx2, dy2, dz2, dx3, dy3, dz3,
     .   side1, side2, oldnx, oldny, oldnz, out, ps, nr,
     .   lambda, em, ep
      my_real
     .   x1(mvsiz),  y1(mvsiz),  z1(mvsiz),
     .   x2(mvsiz),  y2(mvsiz),  z2(mvsiz),
     .   x3(mvsiz),  y3(mvsiz),  z3(mvsiz),
     .   x12(mvsiz), y12(mvsiz), z12(mvsiz),
     .   x13(mvsiz), y13(mvsiz), z13(mvsiz),
     .   x23(mvsiz), y23(mvsiz), z23(mvsiz),
     .   xm(mvsiz) , ym(mvsiz) , zm(mvsiz) , 
     .   xtk2(mvsiz) , ytk2(mvsiz) , ztk2(mvsiz) , 
     .   d(mvsiz)
C-----------------------------------------------
      adrbuf=igrsurf(ksurf)%IAD_BUFR
      a =bufsf(adrbuf+1)
      b =bufsf(adrbuf+2)
      c =bufsf(adrbuf+3)
      dgr =bufsf(adrbuf+36)
      idg =nint(dgr)
      dgr =idg
      dgr =max(two,dgr)
      expn=dgr/(dgr-one)
      an=a**dgr
      bn=b**dgr
      cn=c**dgr
      an=one/max(em20,an)
      bn=one/max(em20,bn)
      cn=one/max(em20,cn)      
      DO i=1,9
       rot(i)=bufsf(adrbuf+7+i-1)
      END DO
C-------------------------------
C     Passage au repere local :
C-------------------------------
      DO 75 nls=ndeb+1,min(ndeb+mvsiz,ntc)
      il =ktc(nls)
      i  =nls-ndeb
      in1=ksi(1,il)
      in2=ksi(2,il)
      in3=ksi(3,il)
C     Passage au repere local.
      xg=x(1,in1)-bufsf(adrbuf+4)
      yg=x(2,in1)-bufsf(adrbuf+5)
      zg=x(3,in1)-bufsf(adrbuf+6)
      xp1(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
      xp1(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
      xp1(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
      xg=x(1,in2)-bufsf(adrbuf+4)
      yg=x(2,in2)-bufsf(adrbuf+5)
      zg=x(3,in2)-bufsf(adrbuf+6)
      xp2(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
      xp2(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
      xp2(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
      xg=x(1,in3)-bufsf(adrbuf+4)
      yg=x(2,in3)-bufsf(adrbuf+5)
      zg=x(3,in3)-bufsf(adrbuf+6)
      xp3(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
      xp3(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
      xp3(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
 75   CONTINUE
C-------------------------------
      DO 100 nls=ndeb+1,min(ndeb+mvsiz,ntc)
      il =ktc(nls)
      i  =nls-ndeb
C-----
      x1(i)=xp1(1,i)
      y1(i)=xp1(2,i)
      z1(i)=xp1(3,i)
      x2(i)=xp2(1,i)
      y2(i)=xp2(2,i)
      z2(i)=xp2(3,i)
      x3(i)=xp3(1,i)
      y3(i)=xp3(2,i)
      z3(i)=xp3(3,i)
      x12(i)=x2(i)-x1(i)
      y12(i)=y2(i)-y1(i)
      z12(i)=z2(i)-z1(i)
      x13(i)=x3(i)-x1(i)
      y13(i)=y3(i)-y1(i)
      z13(i)=z3(i)-z1(i)
      n1=y12(i)*z13(i)-z12(i)*y13(i)
      n2=z12(i)*x13(i)-x12(i)*z13(i)
      n3=x12(i)*y13(i)-y12(i)*x13(i)
      ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
      ntx(i)=ntn*n1
      nty(i)=ntn*n2
      ntz(i)=ntn*n3
      d(i)    =-ntx(i)*x1(i)-nty(i)*y1(i)-ntz(i)*z1(i)
C-----
      x23(i)=x3(i)-x2(i)
      y23(i)=y3(i)-y2(i)
      z23(i)=z3(i)-z2(i)
 100  CONTINUE
C-------------------------------
C     POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
C-------------------------------
      DO 125 nls=ndeb+1,min(ndeb+mvsiz,ntc)
      il =ktc(nls)
      i  =nls-ndeb
C-----
      eh  =(abs(ntx(i)/(dgr*an))**expn)*an
     .    +(abs(nty(i)/(dgr*bn))**expn)*bn
     .    +(abs(ntz(i)/(dgr*cn))**expn)*cn
C     X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
C     LAMBDA1=EXP(LOG(1./EH)/EXPN)
      lambda1=(max(em20,eh))**(-one/expn)
      xh1   =abs(lambda1*ntx(i)/(dgr*an))**(one/(dgr-one))
      IF (ntx(i)<zero) xh1=-xh1
      yh1   =abs(lambda1*nty(i)/(dgr*bn))**(one/(dgr-one))
      IF (nty(i)<zero) yh1=-yh1
      zh1   =abs(lambda1*ntz(i)/(dgr*cn))**(one/(dgr-one))
      IF (ntz(i)<zero) zh1=-zh1
      mu1   =-ntx(i)*xh1-nty(i)*yh1-ntz(i)*zh1-d(i)
C     LAMBDA2=-LAMBDA1
      xh2   =-xh1
      yh2   =-yh1
      zh2   =-zh1
C      MU2   =-NTX(I)*XH2-NTY(I)*YH2-NTZ(I)*ZH2-D(I)
      mu2=-mu1-two*d(i)
      xtk(i) =xh1+mu1*ntx(i)
      ytk(i) =yh1+mu1*nty(i)
      ztk(i) =zh1+mu1*ntz(i)
      xtk2(i)=xh2+mu2*ntx(i)
      ytk2(i)=yh2+mu2*nty(i)
      ztk2(i)=zh2+mu2*ntz(i)
C-------------------------------
 125  CONTINUE
C-------------------------------
C     RAMENER LE PT PK SUR LE TRIANGLE.
C-------------------------------
      DO 150 nls=ndeb+1,min(ndeb+mvsiz,ntc)
      il =ktc(nls)
      i  =nls-ndeb
C-----
      dx1=xtk(i)-x1(i)
      dy1=ytk(i)-y1(i)
      dz1=ztk(i)-z1(i)
      dx2=xtk(i)-x2(i)
      dy2=ytk(i)-y2(i)
      dz2=ztk(i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(i)
     .     +(dz1*dx2-dz2*dx1)*nty(i)
     .     +(dx1*dy2-dy1*dx2)*ntz(i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(i)=alp*x1(i)+bet*x2(i)
        ytk(i)=alp*y1(i)+bet*y2(i)
        ztk(i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk(i)-x3(i)
      dy3=ytk(i)-y3(i)
      dz3=ztk(i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(i)
     .     +(dz2*dx3-dz3*dx2)*nty(i)
     .     +(dx2*dy3-dy2*dx3)*ntz(i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(i)=alp*x2(i)+bet*x3(i)
        ytk(i)=alp*y2(i)+bet*y3(i)
        ztk(i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(i)
     .     +(dz3*dx1-dz1*dx3)*nty(i)
     .     +(dx3*dy1-dy3*dx1)*ntz(i)
      IF (out<0.) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk(i)=alp*x3(i)+bet*x1(i)
        ytk(i)=alp*y3(i)+bet*y1(i)
        ztk(i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-----
      dx1=xtk2(i)-x1(i)
      dy1=ytk2(i)-y1(i)
      dz1=ztk2(i)-z1(i)
      dx2=xtk2(i)-x2(i)
      dy2=ytk2(i)-y2(i)
      dz2=ztk2(i)-z2(i)
      out = (dy1*dz2-dy2*dz1)*ntx(i)
     .      +(dz1*dx2-dz2*dx1)*nty(i)
     .      +(dx1*dy2-dy1*dx2)*ntz(i)
      IF (out<zero) THEN
C       PROJECTION SUR 1,2
        ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
        nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk2(i)=alp*x1(i)+bet*x2(i)
        ytk2(i)=alp*y1(i)+bet*y2(i)
        ztk2(i)=alp*z1(i)+bet*z2(i)
      ENDIF 
      dx3=xtk2(i)-x3(i)
      dy3=ytk2(i)-y3(i)
      dz3=ztk2(i)-z3(i)
      out = (dy2*dz3-dy3*dz2)*ntx(i)
     .      +(dz2*dx3-dz3*dx2)*nty(i)
     .      +(dx2*dy3-dy2*dx3)*ntz(i)
      IF (out<zero) THEN
C       PROJECTION SUR 2,3
        ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
        nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk2(i)=alp*x2(i)+bet*x3(i)
        ytk2(i)=alp*y2(i)+bet*y3(i)
        ztk2(i)=alp*z2(i)+bet*z3(i)
      ENDIF 
      out = (dy3*dz1-dy1*dz3)*ntx(i)
     .      +(dz3*dx1-dz1*dx3)*nty(i)
     .      +(dx3*dy1-dy3*dx1)*ntz(i)
      IF (out<0.) THEN
C       PROJECTION SUR 3,1
        ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
        nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
        bet=ps/max(em20,nr)
        bet=max(bet,zero)
        bet=min(bet,one)
        alp=1.-bet
        xtk2(i)=alp*x3(i)+bet*x1(i)
        ytk2(i)=alp*y3(i)+bet*y1(i)
        ztk2(i)=alp*z3(i)+bet*z1(i)
      ENDIF 
C-------------------------------
 150  CONTINUE
C-------------------------------
C     PROJECTION DE PK SUR L ET PENETRATION.
C-------------------------------
      DO 175 nls=ndeb+1,min(ndeb+mvsiz,ntc)
      il =ktc(nls)
      i  =nls-ndeb
      IF (iactiv(1,il)==-1) GOTO 175
C-----
      xh =xtk(i)
      yh =ytk(i)
      zh =ztk(i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=zero) sgnykn=one      
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one
      n11 =sgnxkn*xkn1*an
      n12 =sgnykn*ykn1*bn
      n13 =sgnzkn*zkn1*cn
      nr1 =n11*n11+n12*n12+n13*n13
      nr1 =sqrt(nr1)
      em =n11*xtk(i)+n12*ytk(i)+n13*ztk(i)
      IF (em<=one) THEN
        lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr1)
        inside1=1
      ELSE
        inside1=0
      ENDIF
C-----
      xh =xtk2(i)
      yh =ytk2(i)
      zh =ztk2(i)
      xkn1 =xh**(idg-1)
      sgnxkn=-one
      IF (xkn1*xh>=zero) sgnxkn=one
      ykn1 =yh**(idg-1)
      sgnykn=-one
      IF (ykn1*yh>=zero) sgnykn=one
      zkn1 =zh**(idg-1)
      sgnzkn=-one
      IF (zkn1*zh>=zero) sgnzkn=one
      n21 =sgnxkn*xkn1*an
      n22 =sgnykn*ykn1*bn
      n23 =sgnzkn*zkn1*cn
      nr2 =n21*n21+n22*n22+n23*n23
      nr2 =sqrt(nr2)
      em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
      IF (em<=one) THEN
        lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
     .         / max(exp((dgr-one)*log(em20)/dgr),nr2)
        inside2=1
      ELSE
        inside2=0
      ENDIF
C-----
      IF (inside1==0.AND.inside2==0) THEN
        iactiv(1,il)=0
      ELSE
C-----
      IF (iactiv(1,il)==0) THEN
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
          xm(i)=xtk(i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(i)-lambda1*n13/max(em20,nr1)
         ELSE
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(i)=xtk2(i)
          ytk(i)=ytk2(i)
          ztk(i)=ztk2(i)
         ENDIF
        ELSEIF(inside1/=0) THEN
          xm(i)=xtk(i)-lambda1*n11/max(em20,nr1)
          ym(i)=ytk(i)-lambda1*n12/max(em20,nr1)
          zm(i)=ztk(i)-lambda1*n13/max(em20,nr1)
        ELSEIF(inside2/=0) THEN
          xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
          ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
          zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
          xtk(i)=xtk2(i)
          ytk(i)=ytk2(i)
          ztk(i)=ztk2(i)
        ENDIF
C-----
      ELSE
        xh=hold(1,4*(il-1)+1)
        yh=hold(2,4*(il-1)+1)
        zh=hold(3,4*(il-1)+1)
        n1=nold(1,4*(il-1)+1)
        n2=nold(2,4*(il-1)+1)
        n3=nold(3,4*(il-1)+1)
        lambda1=(xh-xtk(i))*n1
     .         +(yh-ytk(i))*n2
     .         +(zh-ztk(i))*n3
        lambda2=(xh-xtk2(i))*n1
     .         +(yh-ytk2(i))*n2
     .         +(zh-ztk2(i))*n3
        IF (inside1/=0.AND.inside2/=0) THEN
         IF (abs(lambda1)>=abs(lambda2)) THEN
           xm(i)=xtk(i)+lambda1*n1
           ym(i)=ytk(i)+lambda1*n2
           zm(i)=ztk(i)+lambda1*n3
          ELSE
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(i)=xtk2(i)
           ytk(i)=ytk2(i)
           ztk(i)=ztk2(i)
          ENDIF
        ELSEIF(inside1/=0) THEN
           xm(i)=xtk(i)+lambda1*n1
           ym(i)=ytk(i)+lambda1*n2
           zm(i)=ztk(i)+lambda1*n3
        ELSEIF(inside2/=0) THEN
           xm(i)=xtk2(i)+lambda2*n1
           ym(i)=ytk2(i)+lambda2*n2
           zm(i)=ztk2(i)+lambda2*n3
           xtk(i)=xtk2(i)
           ytk(i)=ytk2(i)
           ztk(i)=ztk2(i)
        ENDIF
C       one more iteration.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one        
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        em=n1*xm(i)+n2*ym(i)+n3*zm(i)
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
        n1 =sgnxkn*xkn1*an
        n2 =sgnykn*ykn1*bn
        n3 =sgnzkn*zkn1*cn
        nr =n1*n1+n2*n2+n3*n3
        nr =one/max(em20,sqrt(nr))
        n1 =n1*nr
        n2 =n2*nr
        n3 =n3*nr
        lambda1=(xm(i)-xtk(i))*n1
     .         +(ym(i)-ytk(i))*n2
     .         +(zm(i)-ztk(i))*n3
        xm(i)=xtk(i)+lambda1*n1
        ym(i)=ytk(i)+lambda1*n2
        zm(i)=ztk(i)+lambda1*n3
      ENDIF 
      iactiv(1,il)=iactiv(1,il)+1
      ENDIF
 175  CONTINUE
C-----
#include "vectorize.inc"
      DO 200 nls=ndeb+1,min(ndeb+mvsiz,ntc)
      il =ktc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(1,il)<=0) GOTO 200
C-----
C      projection radiale de M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       em=n1*xm(i)+n2*ym(i)+n3*zm(i)
C------
        xm(i)=xm(i)/max(em20,em**(one/dgr))
        ym(i)=ym(i)/max(em20,em**(one/dgr))
        zm(i)=zm(i)/max(em20,em**(one/dgr))        
C-----
C      normale a l'ellipsoide en M.
        xkn1 =xm(i)**(idg-1)
        sgnxkn=-one
        IF (xkn1*xm(i)>=zero) sgnxkn=one        
        ykn1 =ym(i)**(idg-1)
        sgnykn=-one
        IF (ykn1*ym(i)>=zero) sgnykn=one        
        zkn1 =zm(i)**(idg-1)
        sgnzkn=-one
        IF (zkn1*zm(i)>=zero) sgnzkn=one
       n1 =sgnxkn*xkn1*an
       n2 =sgnykn*ykn1*bn
       n3 =sgnzkn*zkn1*cn
       nr =n1*n1+n2*n2+n3*n3
       nr =one/max(em20,sqrt(nr))
C-----
       nxi(i)=n1*nr
       nyi(i)=n2*nr
       nzi(i)=n3*nr
       xi(i)=xm(i)
       yi(i)=ym(i)
       zi(i)=zm(i)
       depth(i)=xm(i)*nxi(i)+ym(i)*nyi(i)+zm(i)*nzi(i)
C------
      penet(i)=(xtk(i)-xm(i))*nxi(i)
     .        +(ytk(i)-ym(i))*nyi(i)
     .        +(ztk(i)-zm(i))*nzi(i)
      penet(i)=-penet(i)
      penet(i)=max(zero,penet(i))
      IF (depth(i)-penet(i)<em10*depth(i)) THEN
         penet(i)    =zero
         iactiv(1,il)=-2
      ENDIF
      ansmx=max(ansmx,penet(i))
 200  CONTINUE
C-------------------------------
C     MESSAGES DESACTIVATION
C-------------------------------
      DO 210 nls=ndeb+1,min(ndeb+mvsiz,ntc)
      il =ktc(nls)
      i  =nls-ndeb
C-----
      IF (iactiv(1,il)==-2) THEN
         in1=itab(ksi(1,il))
         in2=itab(ksi(2,il))
         in3=itab(ksi(3,il))
#include "lockon.inc"
         WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(istdo,'(A,A,A,3I8)')' ELEMENT/SEGMENT',
     .                   ' DE-ACTIVATED FROM INTERFACE;',                 
     .                   ' ELEMENT/SEGMENT NODES:',
     .                   in1,in2,in3
         WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
         WRITE(iout,'(A,A,A,3I8)')' ELEMENT/SEGMENT',
     .                   ' DE-ACTIVATED FROM INTERFACE;',                 
     .                   ' ELEMENT/SEGMENT NODES:',
     .                   in1,in2,in3
#include "lockoff.inc"
         iactiv(1,il)=-1
      ENDIF
 210  CONTINUE
C------------------------------------------------------------
      RETURN