C------------------------------------------------      
      PROGRAM ATEST1
C
C***************************************************************
C     ADDENDA, VARIANTE II, TESTS NUMERIQUES.
C     VALEURS DES FONCTIONS DE BASE DE L'ELEMENT GENERIQUE ELARGI
C     ET DE SES DERIVEES.
C     L'OUTPUT EST A COMPARER AVEC CELUI DU PROGRAMME MATHEMATICA.
C***************************************************************
C
C     PAGE 106, SUBROUTINE UT28, BLOC 1, FONCTIONS NUMEROS 1 A 16.
C
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	      NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	      D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION UT(28), DUT(28), U(28), DU(28)
   10 FORMAT(4F15.9)
   15 FORMAT(/)
      WRITE(UNIT=*, FMT=15)
      WRITE(UNIT=*, FMT=*) 'PSEUDO-REFERENCE, VALEURS'
      WRITE(UNIT=*, FMT=15)
C
      CALL UT28(16, 0.1, 0.2, 0.3, UT)
      IND=0
      DO 20, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) UT(IND),UT(IND+1),UT(IND+2),UT(IND+3)
   20 CONTINUE
C
      CALL UT28(16, 0.3, 0.1, 0.2, UT)
      IND=0
      DO 30, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) UT(IND),UT(IND+1),UT(IND+2),UT(IND+3)
   30 CONTINUE
C
      CALL UT28(16, 0.2, 0.3, 0.1, UT)
      IND=0
      DO 40, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) UT(IND),UT(IND+1),UT(IND+2),UT(IND+3)
   40 CONTINUE
C
C     PAGE 106, SUBROUTINE UT28, BLOC 1, FONCTIONS NUMEROS 17 a 28.
C
      MUX(1)=0.5
      MUY(1)=1./3.
      MUZ(1)=1./6.
      MUX(2)=0.5
      MUY(2)=1./3.
      MUZ(2)=1./6.
      MUX(3)=0.5
      MUY(3)=1./3.
      MUZ(3)=1./6.
      MUX(4)=0.5
      MUY(4)=1./3.
      MUZ(4)=1./6.
      WRITE(UNIT=*, FMT=15)
C
      CALL UT28(28, 0.1, 0.2, 0.3, UT)
      DO 120, I=1,4
         IND=3*I+14
         WRITE(UNIT=*, FMT=10) UT(IND),UT(IND+1),UT(IND+2)
  120 CONTINUE
      WRITE(UNIT=*, FMT=15)
C
      CALL UT28(28, 0.3, 0.1, 0.2, UT)
      DO 130, I=1,4
         IND=3*I+14
         WRITE(UNIT=*, FMT=10) UT(IND),UT(IND+1),UT(IND+2)
  130 CONTINUE
      WRITE(UNIT=*, FMT=15)
C
      CALL UT28(28, 0.2, 0.3, 0.1, UT)
      DO 140, I=1,4
         IND=3*I+14
         WRITE(UNIT=*, FMT=10) UT(IND),UT(IND+1),UT(IND+2)
  140 CONTINUE
C
C     PAGE 109, SUBROUTINE DUT28, BLOC 1, FONCTIONS NUMEROS 1 A 16.
C
      WRITE(UNIT=*, FMT=15)
      WRITE(UNIT=*, FMT=*) 'PSEUDO-REFERENCE, DERIVEES'
      WRITE(UNIT=*, FMT=15)
      CALL DUT28(16, 1., 1., 1., 0.1, 0.2, 0.3, DUT)
      DO 220, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  220 CONTINUE
      CALL DUT28(16, 1., 1., 1., 0.3, 0.1, 0.2, DUT)
      DO 230, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  230 CONTINUE
      CALL DUT28(16, 1., 1., 1., 0.2, 0.3, 0.1, DUT)
      DO 240, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  240 CONTINUE
      WRITE(UNIT=*, FMT=15)
C
      CALL DUT28(16, 1., -1., 1., 0.1, 0.2, 0.3, DUT)
      DO 250, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  250 CONTINUE
      CALL DUT28(16, 1., -1., 1., 0.3, 0.1, 0.2, DUT)
      DO 260, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  260 CONTINUE
      CALL DUT28(16, 1., -1., 1., 0.2, 0.3, 0.1, DUT)
      DO 270, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  270 CONTINUE
      WRITE(UNIT=*, FMT=15)
C
      CALL DUT28(16, 1., 1., -1., 0.1, 0.2, 0.3, DUT)
      DO 280, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  280 CONTINUE
      CALL DUT28(16, 1., 1., -1., 0.3, 0.1, 0.2, DUT)
      DO 290, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  290 CONTINUE
      CALL DUT28(16, 1., 1., -1., 0.2, 0.3, 0.1, DUT)
      DO 300, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2),DUT(IND+3)
  300 CONTINUE
C
C     PAGE 109, SUBROUTINE DUT28, BLOC 2, FONCTIONS NUNEROS 17 A 28.
C
      TX=0.25
      TY=0.35
      TZ=0.4
      WRITE(UNIT=*, FMT=15)
C
      CALL DUT28(28, TX,TY,TZ, 0.1, 0.2, 0.3, DUT)
      DO 320, I=1,4
         IND=3*I+14
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2)
  320 CONTINUE
      CALL DUT28(28, TX,TY,TZ, 0.3, 0.1, 0.2, DUT)
      DO 330, I=1,4
         IND=3*I+14
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2)
  330 CONTINUE
      CALL DUT28(28, TX,TY,TZ, 0.2, 0.3, 0.1, DUT)
      DO 340, I=1,4
         IND=3*I+14
         WRITE(UNIT=*, FMT=10) DUT(IND),DUT(IND+1),DUT(IND+2)
  340 CONTINUE
      WRITE(UNIT=*, FMT=15)
C
      CALL ELEM0(0., 0., 0., 5., 2., 0., 0., 4., 1., 1., 0., 3.)
C
      WRITE(UNIT=*, FMT=15)
      WRITE(UNIT=*, FMT=*) 'GENERIQUE ELARGI, VALEURS'
      WRITE(UNIT=*, FMT=15)
      CALL U28(28, 1.4, 1.6, 1.5, U)
      DO 420, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) U(IND),U(IND+1),U(IND+2), U(IND+3) 
  420 CONTINUE
      WRITE(UNIT=*, FMT=15)
      DO 430, I=1,4
         IND=3*I+14
         WRITE(UNIT=*, FMT=10) U(IND),U(IND+1),U(IND+2)
  430 CONTINUE     
C
C     PAGE 112, subroutine DU28
C
      T1=0.25
      T2=0.35
      T3=0.4
C
      CALL ELEM0(0., 0., 0., 5., 2., 0.,
     1           0., 4., 1., 1., 0., 3.)
C
      WRITE(UNIT=*, FMT=15)
      WRITE(UNIT=*, FMT=*) 'GENERIQUE ELARGI, DERIVEES'
      WRITE(UNIT=*, FMT=15)
      CALL DU28(28, T1,T2,T3, 1.4, 1.6, 1.5, DU)
      DO 520, I=1,4
         IND=4*I-3
         WRITE(UNIT=*, FMT=10) DU(IND),DU(IND+1),DU(IND+2), DU(IND+3) 
  520 CONTINUE
      WRITE(UNIT=*, FMT=15)
      DO 530, I=1,4
         IND=3*I+14
         WRITE(UNIT=*, FMT=10) DU(IND),DU(IND+1),DU(IND+2)
  530 CONTINUE     
      STOP
      END
C*************************************************************
C*************************************************************
C*************************************************************
C****
C****
C****    E L E M E N T   F I N I    T E T R A E D R I Q U E 
C****
C****    DE CLASSE C1, DE DEGRE DEUX,  A SEIZE PARAMETRES. 
C****
C****
C****    VARIANTE II.                        MARCEL DELEZE
C****                                        LE CHENE 21
C****    FEVRIER 1978                        1633 MARSENS
C****    OCTOBRE 2008                        SWITZERLAND 
C****
C*************************************************************
C*************************************************************
C------------------------------------------------------------- 
C
      FUNCTION PT1 (J, WT,XT,YT,ZT)
C 
C-------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—REFERENCE.
C**** VALEUR DU POLYNOME DE DEGRE TROIS PT(1,J) AU POINT (XT,YT,ZT).
C**** J PEUT PRENDRE LES VALEURS 0 OU 1. WT=l.—XT—YT—ZT.
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (J) 10,10,11
   10 PT1=W*(2.-W-2.*(X*X+Y*Y+Z*Z+X*Y+Y*Z+Z*X))
      RETURN
   11 PT1=X*W*(W+0.5*(Y+Z))
      RETURN
      END
C---------------------------------------------------------------- 
C
      SUBROUTINE REFCT (WT,XT,YT,ZT, R0,R1,R2,R3)
C
C----------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—REFERENCE. 
C**** VALEURS DES FONCTIONS RATIONNELLES R0, Rl, R2 ET R3
C**** AU POINT (XT,YT,ZT).   WT=1.—XT—YT—ZT.
C**** NOUS SUPPOSONS QUE LE POINT (XT,YT,ZT) NE SOIT PAS SITUE SUR
C**** UNE ARETE DU TETRAEDRE.
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      R0=2.*X*Y*Y*Z*Z/((X+Y)*(X+Z))
      R =2.*X*X*Y*Y*Z*Z*W/((W+X)*(W+Y)*(W+Z)*(1.-Y)*(1.-Z)) 
      R1=R*(-1. - Y*Z*(1.+X)/((1.-Y)*(1.-Z)))
      R2=R*Z*(2.*X+Y)/(1.-Z)
      R3=R*Y*(2.*X+Z)/(1.-Y)
      RETURN
      END
C-------------------------------------------------------------------
C
      FUNCTION WT11 (MUX,MUY,MUZ, WT,XT,YT,ZT)
C
C-------------------------------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO-REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE WT(1,1,MU)
C**** AU POINT (XT,YT,ZT) DU TETRAEDRE DE REFERENCE.  WT=1.-XT-YT-ZT. 
C**** LA DIRECTION NORMALE DEPLACEE (MUX, MUY, MUZ) DEPEND DU TETRAEDRE 
C**** GENERIQUE.	NOUS SUPPOSONS QUE	MUX+MUY+MUZ = 1.
      REAL MUX,MUY,MUZ
      DATA EPS /0.5E-06/
      W=WT
      X=XT
      Y=YT
      Z=ZT
      IF (X.GT.EPS .AND. Y.GT.EPS .AND. Z.GT.EPS) GO TO 30
C**** DANS LES FACES 2, 3 ET 4, ON A 
      WT11=0.0
      RETURN
   30 CONTINUE
      IF (W.GT.EPS) GO TO 40
C**** A L'INTERIEUR DE LA FACE 1, ON A 
      WT11=2.*X*Y*Y*Z*Z/((X+Y)*(X+Z) ) 
      RETURN
   40 CONTINUE
C**** A L'INTERIEUR DU TETRAEDRE, ON A
      CALL REFCT(W,X,Y,Z,R0,R1,R2,R3) 
      WT11=R0 + MUX*R1 + MUY*R2 + MUZ*R3 
      RETURN
      END
c-----------------------------------------------------------------------------
C
      SUBROUTINE UT28 (NFCT, XT,YT,ZT, UT) 
C
C---------------------------–-----------------–----–-------------------------- 
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—REFERENCE.
C**** VALEUR DE UT(1), ..., UT(NFCT) AU POINT (XT,YT,ZT) DU TETRAEDRE DE 
C**** REFERENCE.
C**** NFCT PEUT PRENDRE LES VALEURS 16 OU 28.
C**** LES DIRECTIONS NORMALES DEPLACEES MUX(I), MUY(I), MUZ(I), I=1(1)4, 
C**** DEPENDENT DU TETRAEDRE GENERIQUE.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3) 
      DIMENSION UT(28)
      NFCS=NFCT
      XL=XT
      YL=YT
      ZL=ZT
      WL=1.-XL-YL-ZL
      DO 90 I=1,4
C****	 BLOC 1.	FONCTIONS NUMEROS 1 A 16. 
         IND=4*I-3
         UT(IND  )=PT1(0,WL,XL,YL,ZL) 
         UT(IND+1)=PT1(1,WL,XL,YL,ZL) 
         UT(IND+2)=PT1(1,WL,YL,ZL,XL) 
         UT(IND+3)=PT1(1,WL,ZL,XL,YL) 
         IF (NFCS.LE.16) GO TO 80
C****	 BLOC 2.	FONCTIONS NUMEROS 17 A 28. 
         IND=3*I+14
         UT(IND  )=WT11(MUX(I),MUY(I),MUZ(I),WL,XL,YL,ZL)
         UT(IND+1)=WT11(MUY(I),MUZ(I),MUX(I),WL,YL,ZL,XL) 
         UT(IND+2)=WT11(MUZ(I),MUX(I),MUY(I),WL,ZL,XL,YL)
   80    CH=WL 
         WL=XL
         XL=YL 
         YL=ZL 
         ZL=CH
   90 CONTINUE
      RETURN
      END
C---------------------------------------------------------------
C
      FUNCTION DPT1 (J, TX,TY,TZ, WT,XT,YT,ZT)
C
c------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—REFERENCE.
C**** DERIVEE DU POLYNOME DE DEGRE TROIS PT(1,J) DANS LA DIRECTION
C**** (TX,TY,TZ) AU POINT (XT,YT,ZT).
C**** J PEUT PRENDRE LES VALEURS 0 OU 1.  WT=1.-XT-YT-ZT.
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      T1=TX 
      T2=TY 
      T3=TZ 
      T0=-T1-T2-T3
      IF (J) 10,10,11
   10 P= 2.-W - 2.*(X*X+Y*Y+Z*Z+X*Y+Y*Z+Z*X)
      DP =-T0 - 2.*( T1*(2.*X+Y+Z)+T2*(2.*Y+Z+X)+T3*(2.*Z+X+Y) )
      GO TO 20
   11 P=X*(W+0.5*(Y+Z))
      DP =T0*X + T1*(W+0.5*(Y+Z)) + (T2+T3)*0.5*X
   20 DPT1=T0*P + W*DP
      RETURN
      END
C---------------------------------------------------------------------
C
      SUBROUTINE DREFCT (TX,TY,TZ, WT,XT,YT,ZT, DR0,DR1,DR2,DR3)
C
C----------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO—REFERENCE.
C**** DERIVEES DES FONCTIONS RATIONNELLES R0, R1, R2 ET R3
C**** DANS LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT). WT=1.-XT-YT-ZT. 
C**** NOUS SUPPOSONS QUE LE POINT (XT,YT,ZT) NE SOIT PAS SITUE SUR
C**** UNE ARETE DU TETRAEDRE.
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      T1=TX
      T2=TY 
      T3=TZ 
      T0=-T1-T2-T3
C**** DERIVEE DE R0.
      P =2.*X*Y*Y*Z*Z
      DP =2.*Y*Z*(T1*Y*Z+2.*X*(T2*Z+T3*Y))
      Q =(X+Y)*(X+Z)
      DR0=(DP -  P*( (T1+T2)/(X+Y)+(T1+T3)/(X+Z) ))/Q
C**** DERIVEE DE R.
      P  =2.*X*X*Y*Y*Z*Z*W
      DP =2.*X*Y*Z*( T0*X*Y*Z+2.*W*(T1*Y*Z+T2*X*Z+T3*X*Y) )
      Q  =(W+X)*(W+Y)*(W+Z)*(1.-Y)*(1.-Z)
      R  =P/Q
      DR =(DP  -   P*( (T0+T1)/(W+X)+(T0+T2)/(W+Y)+(T0+T3)/(W+Z)
     1	               -T2/(1.-Y)-T3/(1.-Z)  ))/Q
C**** DERIVEE DE R1.
      P	 = R*Y*Z*(1.+X)
      DP =DR*Y*Z*(1.+X) + R*( T1*Y*Z+(1.+X)*(T2*Z+T3*Y) )
      Q	 =(1.-Y)*(1.-Z)
      DR1=-DR  -  (DP - P*( -T2/(1.-Y)-T3/(1.-Z)  ))/Q
C**** DERIVEE DE R2.
      P  =R*Z*(2.*X+Y)
      DP =DR*Z*(2.*X+Y) + R*( (T1*2.+T2)*Z+T3*(2.*X+Y) )
      Q  =1.-Z
      DR2=(DP - P*(-T3/Q))/Q
C**** DERIVEE DE R3.
      P  =R*Y*(2.*X+Z)
      DP =DR*Y*(2.*X+Z) + R*( Y*(2.*T1+T3)+T2*(2.*X+Z) )
      Q	 =1.-Y
      DR3=(DP - P*(-T2/Q))/Q
      RETURN
      END
C------------------------------------------------------------------
C
      FUNCTION DWT11 (TX,TY,TZ, MUX,MUY,MUZ, WT,XT,YT,ZT)
C
C------------------------------------------------------------------	 
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO-REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE WT(1,1,MU) DANS LA DIRECTION 
C**** (TX,TY,TZ) AU POINT (XT,YT,ZT) DU TETRAEDRE DE REFERENCE.
C**** WT=1.-XT-YT-ZT.
C**** LA DIRECTION NORMALE DEPLACEE (MUX, MUY, MUZ) DEPEND DU TETRAEDRE 
C**** GENERIQUE. NOUS SUPPOSONS QUE	MUX+MUY+MUZ = 1.
      REAL MUX,MUY,MUZ
      DATA  EPS /0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (Y.GT.EPS .AND. W+Y.GT.EPS .AND.
     1    Z.GT.EPS .AND. W+Z.GT.EPS) GO TO 30
C**** DANS LES FACES 3 ET 4, ON A
      DWT11=0.0
      RETURN
   30 CONTINUE
      T1=TX 
      T2=TY 
      T3=TZ
      IF (X.GT.EPS .AND. W+X.GT.EPS) GO TO 40
C**** DANS LA FACE 2, ON A
      DWT11=2.*T1*Y*Z
      RETURN
   40 CONTINUE
C**** A L'INTERIEUR DU TETRAEDRE ET A L'INTERIEUR DE LA FACE 1, ON A 
      CALL DREFCT(T1,T2,T3,W,X,Y,Z,DR0,DR1,DR2,DR3)
      DWT11=DR0 + MUX*DR1 + MUY*DR2 + MUZ*DR3
      RETURN
      END
C----------------------------------------------------------------------
C
      SUBROUTINE DUT28 (NFCT, TX,TY,TZ, XT,YT,ZT, DUT)
C
C-----------------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE PSEUDO-REFERENCE.
C**** DERIVEE DE UT(1), ..., UT(NFCT) DANS LA DIRECTION (TX,TY,TZ) 
C**** AU POINT (XT,YT,ZT) DU TETRAEDRE DE REFERENCE.
C**** NFCT PEUT PRENDRE LES VALEURS 16 OU 28.
C**** LES DIRECTIONS NORMALES DEPLACEES MUX(I), MUY(I), MUZ(I),  I=1(1)4,
C**** DEPENDENT DU TETRAEDRE GENERIQUE.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	      NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	      D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3) 
      DIMENSION DUT(28)
      NFCS=NFCT 
      XL=XT
      YL=YT
      ZL=ZT
      WL=1.-XL-YL-ZL
      T1=TX
      T2=TY
      T3=TZ
      T0=-T1-T2-T3
      DO 90 I=1,4
C****    BLOC 1.	FONCTIONS NUMEROS 1 A 16.
         IND=4*I-3
         DUT(IND  )=DPT1(0,T1,T2,T3,WL,XL,YL,ZL) 
         DUT(IND+1)=DPT1(1,T1,T2,T3,WL,XL,YL,ZL) 
         DUT(IND+2)=DPT1(1,T2,T3,T1,WL,YL,ZL,XL) 
         DUT(IND+3)=DPT1(1,T3,T1,T2,WL,ZL,XL,YL)
         IF (NFCS.LE.16) GO TO 80
C****    BLOC 2.	FONCTIONS NUMEROS 17 A 28.
         IND=3*I+14
         DUT(IND  )=DWT11(T1,T2,T3,MUX(I),MUY(I),MUZ(I),WL,XL,YL,ZL)
         DUT(IND+1)=DWT11(T2,T3,T1,MUY(I),MUZ(I),MUX(I),WL,YL,ZL,XL)
         DUT(IND+2)=DWT11(T3,T1,T2,MUZ(I),MUX(I),MUY(I),WL,ZL,XL,YL)
   80    CH=WL 
         WL=XL 
         XL=YL 
         YL=ZL 
         ZL=CH 
         CH=T0 
         T0=T1 
         T1=T2 
         T2=T3 
         T3=CH
   90 CONTINUE
      RETURN
      END
C-----------------------------------------------------------------------------
C
      SUBROUTINE ELEM0 (X1,Y1,Z1, X2,Y2,Z2, X3,Y3,Z3, X4,Y4,Z4)
C
C------------------------------------------------------------------------------
C**** ELEMENT GENERIQUE ELARGI.
C**** INITIALISATIONS RELATIVES AU TETRAEDRE GENERIQUE.
C**** ENTREES	: X1,Y1,Z1, X2,Y2,Z2, X3,Y3,Z3, X4,Y4,Z4,  SOMMETS DU
C****                    TETRAEDRE GENERIQUE.
C**** SORTIES : TETRAS BLOC COMMUN PAR LEQUEL SE FONT LES SORTIES.
C****       LL      PARTIE LINEAIRE DE L'APPLICATION AFFINE L,
C****               FACTORISEE SELON LA METHODE DE GAUSS.
C****       IPERM   VECTEUR CONTENANT LA NOUVELLE NUMEROTATION DES
C****	            LIGNES DE LL, CONSECUTIVE AU CHOIX DES PIVOTS.
C****	    NX(I,J), NY(I,J), NZ(I,J), J=1,2,3,  I=1(1)4, VECTEURS 
C****               PARALLELES A LA I-EME FACE ET NORMAUX AUX
C****	            ARETES DU TETRAEDRE.
C****       MUX(I), MUY(I), MUZ(I), I=1(1)4, DIRECTIONS NORMALES 
C****               DEPLACEES. CHAQUE VECTEUR A ETE MULTIPLIE PAR 
C****	            UN SCALAIRE TEL QUE LE PRODUIT SATISFASSE
C****               MUX(I) + MUY(I) + MUZ(I) = 1,	I=1(1)4.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	     D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION GEN(3),REF(3)
      X(1)=X1 
      Y(1)=Y1 
      Z(1)=Z1 
      X(2)=X2 
      Y(2)=Y2 
      Z(2)=Z2 
      X(3)=X3 
      Y(3)=Y3 
      Z(3)=Z3 
      X(4)=X4 
      Y(4)=Y4 
      Z(4)=Z4
C**** PARTIE LINEAIRE DE L'APPLICATION AFFINE L.
      LL(1,1)=X(2)-X(1) 
      LL(2,1)=Y(2)-Y(1) 
      LL(3,1)=Z(2)-Z(1)
      LL(1,2)=X(3)-X(1)
      LL(2,2)=Y(3)-Y(1) 
      LL(3,2)=Z(3)-Z(1) 
      LL(1,3)=X(4)-X(1) 
      LL(2,3)=Y(4)-Y(1) 
      LL(3,3)=Z(4)-Z(1)
      CALL DECOMP(3,3,LL,GEN,IPERM,DET)
C**** VOL = ABS(DET)/6. EST LE VOLUME DU TETRAEDRE.
C**** SI LE TETRAEDRE EST DEGENERE (VOL=0.), LA SUBROUTINE DECOMP 
C**** IMPRIME UN MESSAGE D'ERREUR.
      DO 60 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1 
         L=MOD(K,4)+1
C****	 VECTEUR NORMAL A LA I-EME FACE, PLACE DANS GEN. 
         GEN(1)=(Y(K)-Y(J))*(Z(L)-Z(J))-(Z(K)-Z(J))*(Y(L)-Y(J)) 
         GEN(2)=(Z(K)-Z(J))*(X(L)-X(J))-(X(K)-X(J))*(Z(L)-Z(J)) 
         GEN(3)=(X(K)-X(J))*(Y(L)-Y(J))-(Y(K)-Y(J))*(X(L)-X(J))
C****    VECTEURS PARALLELES A LA I-EME FACE ET NORMAUX AUX ARETES. 
         NX(I,1)=GEN(2)*(Z(L)-Z(K))-GEN(3)*(Y(L)-Y(K)) 
         NY(I,1)=GEN(3)*(X(L)-X(K))-GEN(1)*(Z(L)-Z(K)) 
         NZ(I,1)=GEN(1)*(Y(L)-Y(K))-GEN(2)*(X(L)-X(K)) 
         S=SQRT(NX(I,1)**2+NY(I,1)**2+NZ(I,1)**2)
         NX(I,1)=NX(I,1)/S 
         NY(I,1)=NY(I,1)/S
         NZ(I,1)=NZ(I,1)/S 
         NX(I,2)=GEN(2)*(Z(J)-Z(L))-GEN(3)*(Y(J)-Y(L)) 
         NY(I,2)=GEN(3)*(X(J)-X(L))-GEN(1)*(Z(J)-Z(L))
         NZ(I,2)=GEN(1)*(Y(J)-Y(L))-GEN(2)*(X(J)-X(L)) 
         S=SQRT(NX(I,2)**2+NY(I,2)**2+NZ(I,2)**2)
         NX(I,2)=NX(I,2)/S
         NY(I,2)=NY(I,2)/S
         NZ(I,2)=NZ(I,2)/S 
         NX(I,3)=GEN(2)*(Z(K)-Z(J))-GEN(3)*(Y(K)-Y(J)) 
         NY(I,3)=GEN(3)*(X(K)-X(J))-GEN(1)*(Z(K)-Z(J)) 
         NZ(I,3)=GEN(1)*(Y(K)-Y(J))-GEN(2)*(X(K)-X(J)) 
         S=SQRT(NX(I,3)**2+NY(I,3)**2+NZ(I,3)**2)
         NX(I,3)=NX(I,3)/S
         NY(I,3)=NY(I,3)/S
         NZ(I,3)=NZ(I,3)/S
C****	 DIRECTION NORMALE DEPLACEE RELATIVE A LA I-EME FACE.
         CALL SOLVE (3,3,LL,GEN,REF,IPERM)
         GO TO (10,20,30,40), I
   10    MUX(1)=REF(1) 
         MUY(1)=REF(2) 
         MUZ(1)=REF(3) 
         GO TO 50
   20    MUX(2)=REF(2) 
         MUY(2)=REF(3) 
         MUZ(2)=-REF(1)-REF(2)-REF(3)
         GO TO 50
   30	 MUX(3)=REF(3) 
         MUY(3)=-REF(1)-REF(2)-REF(3)
         MUZ(3)=REF(1) 
         GO TO 50
   40	 MUX(4)=-REF(1)-REF(2)-REF(3)
         MUY(4)=REF(1) 
         MUZ(4)=REF(2)
   50	 S=MUX(I)+MUY(I)+MUZ(I)
         MUX(I)=MUX(I)/S 
         MUY(I)=MUY(I)/S 
         MUZ(I)=MUZ(I)/S
   60 CONTINUE
      RETURN
      END
C---------------------------------------------------------------
C
      SUBROUTINE U28 (NFCT, PX,PY,PZ, U)
C 
C---------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** VALEUR DE U(1), ..., U(NFCT) AU POINT (PX,PY,PZ) DU TETRAEDRE. 
C**** NFCT PEUT PRENDRE LES VALEURS 16 OU 28.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	     NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	     D(4,3,3),D5(12),Q1(12,16),T(28,16), LL(3,3),IPERM(3) 
      DIMENSION U(28),GEN(3),REF(3)
      DATA	NERR/0/,	EPS/-1.E-03/
      GEN(1)=PX-X(1) 
      GEN(2)=PY-Y(1) 
      GEN(3)=PZ-Z(1)
      CALL SOLVE(3,3,LL,GEN,REF,IPERM)
      XT=REF(1)
      YT=REF(2)
      ZT=REF(3)
      IF (XT.GT.EPS  .AND.  YT.GT.EPS  .AND.  ZT.GT.EPS  .AND.
     1	  1.-XT-YT-ZT.GT.EPS) GO TO 20
C**** LE POINT (PX, PY, PZ) EST A L'EXTERIEUR DU TETRAEDRE GENERIQUE. 
      IF (NERR.LT.10) WRITE(UNIT=*, FMT=*)
     1                ' *** ERROR MESSAGE CALLED FROM U28 :',
     2                ' THE POINT IS OUTSIDE THE ELEMENT.'
      NERR=NERR+1
   20 CONTINUE
C**** VALEURS DES FONCTIONS DE BASE.
      CALL UT28(NFCT,XT,YT,ZT,U)
      RETURN
      END
C------------------------------------------------------------------- 
C
      SUBROUTINE DU28 (NFCT, VX,VY,VZ, PX,PY,PZ, DU)
C 
C-------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** DERIVEE DE U(1), .., U(NFCT) DANS LA DIRECTION (VX, VY, VZ) 
C**** AU POINT (PX,PY,PZ) DU TETRAEDRE GENERIQUE.
C**** NFCT PEUT PRENDRE LES VALEURS 16 OU 28.
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	     NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	     D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION DU(28),GEN(3),REF(3)
      DATA    NERR/0/,   EPS/-1.E-03/
      GEN(1)=VX 
      GEN(2)=VY 
      GEN(3)=VZ
      CALL SOLVE(3,3,LL,GEN,REF,IPERM)
      TX=REF(1) 
      TY=REF(2) 
      TZ=REF(3) 
      GEN(1)=PX-X(1)
      GEN(2)=PY-Y(1)
      GEN(3)=PZ-Z(1)
      CALL SOLVE(3,3,LL,GEN,REF,IPERM)
      XT=REF(1) 
      YT=REF(2) 
      ZT=REF(3)
      IF (XT.GT.EPS .AND. YT.GT.EPS .AND. ZT.GT.EPS .AND.
     1	  1.-XT-YT-ZT.GT.EPS) GO TO 20
C**** LE POINT (PX, PY, PZ) EST A L'EXTERIEUR DU TETRAEDRE GENERIQUE.
      IF (NERR.LT.10) WRITE(UNIT=*, FMT=*)
     1                ' *** ERROR MESSAGE CALLED FROM DU28 :',
     2                ' THE POINT IS OUTSIDE THE ELEMENT.'
      NERR=NERR+1 
   20 CONTINUE
C**** DERIVEES DES FONCTIONS DE BASE.
      CALL DUT28(NFCT,TX,TY,TZ,XT,YT,ZT,DU)
      RETURN
      END
C-------------------------------------------------------------------------- 
C
      SUBROUTINE ELEM1
C
c--------------------------------------------------------------------------	
C**** MATRICE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** INITIALISATION DE L'INVERSE DES BLOCS DIAGONAUX DE LA MATRICE Q. 
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	     NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2	     D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DO 10 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1 
         L=MOD(K,4)+1
C****    INITIALISATION DE D1, D2, D3 OU D4.
         D(I,1,1)=X(J)-X(I)
         D(I,1,2)=Y(J)-Y(I)
         D(I,1,3)=Z(J)-Z(I)
         D(I,2,1)=X(K)-X(I)
         D(I,2,2)=Y(K)-Y(I)
         D(I,2,3)=Z(K)-Z(I) 
         D(I,3,1)=X(L)-X(I)
         D(I,3,2)=Y(L)-Y(I)
         D(I,3,3)=Z(L)-Z(I)
C****    INITIALISATION DE D5.
         IND=3*I-2
         D5(IND  )=NX(I,1)*(2.*X(J)-X(K)-X(L))
     1            +NY(I,1)*(2.*Y(J)-Y(K)-Y(L))
     2            +NZ(I,1)*(2.*Z(J)-Z(K)-Z(L))
         D5(IND+1)=NX(I,2)*(2.*X(K)-X(L)-X(J))
     1            +NY(I,2)*(2.*Y(K)-Y(L)-Y(J))
     2            +NZ(I,2)*(2.*Z(K)-Z(L)-Z(J))
         D5(IND+2)=NX(I,3)*(2.*X(L)-X(J)-X(K))
     1            +NY(I,3)*(2.*Y(L)-Y(J)-Y(K))
     2            +NZ(I,3)*(2.*Z(L)-Z(J)-Z(K))
   10 CONTINUE
      RETURN
      END
C---------------------------------------------------------------
C
      SUBROUTINE ELEM2
C 
C---------------------------------------------------------------
C**** MATRICE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** INITIALISATION DU BLOC Q1 DE Q. 
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ 
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION DU(28)
      DO 40 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1 
         L=MOD(K,4)+1 
         BX=0.5*(X(K)+X(L))
         BY=0.5*(Y(K)+Y(L))
         BZ=0.5*(Z(K)+Z(L))
         CALL DU28(16,NX(I,1),NY(I,1),NZ(I,1),BX,BY,BZ,DU)
         IND=3*I-2
         DO 10 M=1,16
            Q1(IND,M)=DU(M) 
   10    CONTINUE
         BX=0.5*(X(J)+X(L)) 
         BY=0.5*(Y(J)+Y(L)) 
         BZ=0.5*(Z(J)+Z(L))
         CALL DU28(16,NX(I,2),NY(I,2),NZ(I,2),BX,BY,BZ,DU)
         IND=IND+1
         DO 20 M=1,16
            Q1(IND,M)=DU(M) 
   20    CONTINUE
         BX=0.5*(X(K)+X(J)) 
         BY=0.5*(Y(K)+Y(J)) 
         BZ=0.5*(Z(K)+Z(J))
         CALL DU28(16,NX(I,3),NY(I,3),NZ(I,3),BX,BY,BZ,DU)
         IND=IND+1
         DO 30 M=1,16
            Q1(IND,M)=DU(M) 
   30    CONTINUE
   40 CONTINUE
      RETURN
      END
C-------------------------------------------------------------------	
C
      SUBROUTINE FE2P28 (C, S)
C 
C-------------------------------------------------------------------
C**** ELEMENT GENERIQUE ELARGI.
C**** ENTREES : TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT
C****                  INITIALISE A L'AIDE DES SUBROUTINES ELEM0,
C****                  ELEM1 ET ELEM2.
C****           C(1), ..., C(28) VALEURS DES PARAMETRES DE L'ELEMENT.
C**** SORTIES : S(1), ..., S(28) COORDONNEES DE LA FONCTION D'INTER-
C****                  POLATION DANS LA BASE D'INTERPOLATION DE
C****                  L'ELEMENT GENERIQUE ELARGI.
C**** METHODE : RESOLUTION DU SYSTEME LINEAIRE Q*S = C. 
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION C(28),S(28)
C**** CALCUL DE S(1), ..., S(16).
      DO 10 I=1,4
         K=4*I-3
         S(K)=C(K) 
         S(K+1)=D(I,1,1)*C(K+1)+D(I,1,2)*C(K+2)+D(I,1,3)*C(K+3) 
         S(K+2)=D(I,2,1)*C(K+1)+D(I,2,2)*C(K+2)+D(I,2,3)*C(K+3) 
         S(K+3)=D(I,3,1)*C(K+1)+D(I,3,2)*C(K+2)+D(I,3,3)*C(K+3)
   10 CONTINUE
C**** CALCUL DE S(17), ..., S(28).
      DO 30 I=1,12
         SI=0.0
         DO 20 J=1,16
            SI=SI+Q1(I,J)*S(J)
   20    CONTINUE
         K=I+16
         S(K)=D5(I)*(C(K)-SI)
   30 CONTINUE
      RETURN
      END
C------------------------------------------------------------------------ 
C
      FUNCTION D1NOD (VX,VY,VZ, M, C)
C 
C------------------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX,
C**** A 16 PARAMETRES.
C**** DONNEES :   M      PEUT PRENDRE LES VALEURS 1, 2, 3 OU 4.
C****             C(4*M-2), C(4*M-1) ET C(4*M) DERIVEES PARTIELLES DX,
C****                    DY ET DZ AU SOMMET A(M).
C****             (VX, VY, VZ) UN VECTEUR.
C**** SORTIE   :  D1NOD  DERIVEE DANS LA DIRECTION (VX, VY, VZ)
C****                    AU SOMMET A(M).
      DIMENSION C(28)
      IND=4*M-2
      D1NOD=VX*C(IND)+VY*C(IND+1)+VZ*C(IND+2)
      RETURN
      END
C----------------------------------------------------------------
C
      SUBROUTINE FE2P16 (C, S)
C
C----------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX, 
C**** A 16 PARAMETRES.
C**** RESOLUTION DU PROBLEME D'INTERPOLATION.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT
C****                    INITIALISE A L'AIDE DES SUBROUTINES ELEM0,
C****                    ELEM1 ET ELEM2.
C****             C(1), ..., C(16) VALEURS DES PARAMETRES DE L'ELEMENT.
C**** SORTIES  :  S(1), ..., S(28) COORDONNEES DE LA FONCTION D'INTER-
C****                    POLATION DANS LA BASE D'INTERPOLATION DE
C****                    L'ELEMENT GENERIQUE ELARGI.
C**** LES MEMOIRES C(17), ..., C(28) SONT DES PLACES DE TRAVAIL. 
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION C(28),S(28)
C**** CALCUL DES PARAMETRES C(17), ..., C(28).
      DO 10 I=1,4
         J=MOD(I,4)+1
         K=MOD(J,4)+1
         L=MOD(K,4)+1
         IND=3*I+14
         C(IND  )=0.5*(D1NOD(NX(I,1),NY(I,1),NZ(I,1),K,C)
     1                +D1NOD(NX(I,1),NY(I,1),NZ(I,1),L,C))
         C(IND+1)=0.5*(D1NOD(NX(I,2),NY(I,2),NZ(I,2),L,C)
     1                +D1NOD(NX(I,2),NY(I,2),NZ(I,2),J,C))
         C(IND+2)=0.5*(D1NOD(NX(I,3),NY(I,3),NZ(I,3),J,C)
     1                +D1NOD(NX(I,3),NY(I,3),NZ(I,3),K,C))
   10 CONTINUE
C**** CALCUL DES COORDONNEES S(1), ..., S(28).
      CALL FE2P28(C,S)
      RETURN
      END
C----------------------------------------------------------------- 
C
      FUNCTION FINT (PX,PY,PZ, S)
C 
C-----------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX,
C**** A 16 PARAMETRES.
C**** DONNEES : TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****                  INITIALISE A L'AIDE DE LA SUBROUTINE ELEM0.
C****	        S(1), ..., S(28), COORDONNEES DE LA FONCTION D'INTER-
C****	               POLATION DANS LA BASE D'INTERPOLATION DE
C****                  L'ELEMENT GENERIQUE ELARGI.
C****           (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE  :  FINT    VALEUR DE LA FONCTION D'INTERPOLATION DE
C****                    COORDONNEES S AU POINT (PX, PY, PZ).
      DIMENSION U(28),S(28)
      CALL U28(28,PX,PY,PZ,U)
      V=0.0
      DO 10 I=1,28
         V=V+S(I)*U(I)
   10 CONTINUE 
      FINT=V 
      RETURN 
      END
C-------------------------------------------------------------------	
C
      FUNCTION DFINT (VX,VY,VZ, PX,PY,PZ, S)
C
C-------------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX,
C**** A 16 PARAMETRES.
C**** DONNEES : TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****                  INITIALISE A L'AIDE DE LA SUBROUTINE ELEM0.
C****	        S(1),... ,S(28), COORDONNEES DE LA FONCTION D'INTER-
C****                  POLATION DANS LA BASE D'INTERPOLATION DE
C****                  L'ELEMENT GENERIQUE ELARGI.
C****           (VX, VY, VZ) UN VECTEUR.
C****           (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE  : DFINT DERIVEE DE LA FONCTION D'INTERPOLATION DE
C****                 COORDONNEES S DANS LA DIRECTION (VX, VY, VZ)
C****                 AU POINT (PX,PY,PZ).
      DIMENSION DU(28),S(28)
      CALL DU28(28,VX,VY,VZ,PX,PY,PZ,DU)
      DV=0.0
      DO 10 I=1,28
         DV=DV+S(I)*DU(I)
   10 CONTINUE 
      DFINT=DV 
      RETURN
      END
C--------------------------------------------------------------------
C
      SUBROUTINE ELEM3
C
C--------------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX, 
C**** A 16 PARAMETRES,
C**** BASE D'HERMITE.	INITIALISATION DE LA MATRICE T. 
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION C(28), S(28)
      DATA C/28*0.0/
      DO 30 I=1,16
C****    CALCUL DES COEFFICIENTS T(1,I), ... , T(28,I) DU I-EME ELEMENT
C****    DE LA BASE D'HERMITE PAR RAPPORT A LA BASE D'INTERPOLATION DE
C****    L'ELEMENT GENERIQUE ELARGI.
         C(I)=1.0
         CALL FE2P16(C,S)
         DO 20 J=1,28
            T(J,I)=S(J)
   20    CONTINUE
         C(I)=0.0 
   30 CONTINUE
      RETURN
      END
C------------------------------------------------------------------
C
      SUBROUTINE HERM (PX,PY,PZ, V)
C
C------------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX, 
C**** A 16 PARAMETRES.
C**** BASE D'HERMITE.
C**** ENTREES : TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT
C****                  INITIALISE A L'AIDE DES SUBROUTINES ELEM0,
C****                  ELEM1,ELEM2 ET ELEM3.
C****	        (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIES : V(1), ..., V(16) VALEURS DES 16 FONCTIONS DE LA BASE
C****	               D'HERMITE AU POINT (PX, PY, PZ).
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION V(16), U(28)
      CALL U28(28,PX,PY,PZ,U)
      DO 20 I=1,16
         SUM=0.0
         DO 10 J=1,28
            SUM=SUM + T(J,I)*U(J)
   10    CONTINUE 
         V(I)=SUM 
   20 CONTINUE
      RETURN
      END
C--------------------------------------------------------------------
C
      SUBROUTINE DHERM (VX,VY,VZ, PX,PY,PZ, DV)
C
C--------------------------------------------------------------------
C**** ELEMENT FINI GENERIQUE, DE CLASSE C1, DE DEGRE DEUX, 
C**** A 16 PARAMETRES.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT
C****                    INITIALISE A L'AIDE DES SUBROUTINES ELEM0,
C****                    ELEM1, ELEM2 ET ELEM3.
C****             (VX, VY, VZ) UN VECTEUR.
C****             (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIES  :  DV(1), ..., DV(16) DERIVEES DES 16 FONCTIONS DE LA
C****                    BASE D'HERMITE DANS LA DIRECTION (VX, VY, VZ)
C****                    AU POINT (PX, PY, PZ). 
      REAL LL, MUX,MUY,MUZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MUX(4),MUY(4),MUZ(4),
     2       D(4,3,3),D5(12),Q1(12,16), T(28,16), LL(3,3),IPERM(3)
      DIMENSION DV(16), DU(28)
      CALL DU28(28,VX,VY,VZ,PX,PY,PZ,DU)
      DO 20 I=1,16
         SUM=0.0
         DO 10 J=1,28
            SUM=SUM + T(J,I)*DU(J)
   10    CONTINUE 
         DV(I)=SUM
   20 CONTINUE 
      RETURN 
      END
c-------------------------------------------------------------    
C
      SUBROUTINE DECOMP (NX,NN,A,SCALES,IPS,D) 
C
C     THIS ROUTINE DOES AN LU DECOMPOSITION ON A SQUARE MATRIX.
C     IT USES ROW SCALING AND ROW PIVOTING TO MAKE THE PROCESS
C     AS ACCURATE AS POSSIBLE.
C     THE PARAMETERS OF THE CALLING SEQUENCE ARE
C     N       THE ORDER OF THE MATRIX TO BE DECOMPOSED
C     NX      THE NUMBER OF LINES OF THE MATRIX A
C	      IN THE CALLING PROGRAM
C     A       THE MATRIX TO BE DECOMPOSED
C	      THE DECOMPOSED MATRIX. THE MAIN DIAGONAL AND EVERYTHING
C	      ABOVE IT MAKE UP THE UPPER TRIANGULAR MATRIX. THE LOWER
C	      TRIANGULAR MATRIX IS MADE UP OF EVERYTHING BELOW THE MAIN
C	      DIAGONAL. THE MAIN DIAGONAL FOR THE LOWER TRIANGULAR
C	      MATRIX IS ASSUMED TO BE ALL ONES.
C     D       THE DETERMINANT OF THE COEFFICIENT MATRIX
C     KERR    ERROR PRINT CONTROL AND RETURN ERROR FLAG
C             IF KERR SET NEGATIVE BY USER NO MESSAGE PRINTED
C             IF SET TO ZERO, ERROR MESSAGE WILL BE PRINTED.
C             ON RETURN, FOLLOWING CODE -
C             =1-MATRIX WITH ZERO ROW
C             =2-SINGULAR MATRIX
C             =3- SOLUTION OBTAINED
C     IPS    VECTOR CONTAINING PERMUTATED ROW NUMBERS OF A 
C
      DIMENSION A(NX,NX), SCALES(NX), IPS(NX)
      N=NN
      KERR=0
C
C     FIND THE SOLUTION IN THE CASE N=1
C
      IF(N.GT.1) GO TO 20 
      D=A(1,1)
      IPS(1)=1
      IF(D.EQ.0.) GO TO 10
      KERR=3 
      RETURN
   10 IF(KERR.EQ.0) CALL SING(1)
      KERR=1 
      RETURN
C
C	INITIALIZE
C
   20 ODD=1.
C
C	FIND THE SCALING FACTOR FOR EACH ROW AND STORE IN SCALES
C
      DO 40 I=1,N 
      IPS(I)=I 
      RMAX=0.
      DO 30 J=1,N
         RMAX=AMAX1(RMAX,ABS(A(I,J))) 
   30 CONTINUE
      IF(RMAX.EQ.0.) GO TO 110 
      SCALES(I)=1./RMAX
   40 CONTINUE
C
C     MAIN LOOP. GAUSSIAN ELIMINATION IS CARRIED OUT. LOOP IS GONE
C     THROUGH ONCE FOR EACH ELEMENT ON THE MAIN DIAGONAL . K IS THE
C     DIAGONAL ELEMENT NUMBER. LAST DIAGONAL ELEMENT IS AN EXCEPTION. 
C
      NM1=N-1
      DO 90 K=1,NM1 
         BIG=0.
         DO 50 I=K,N 
            IP=IPS(I) 
            SIZE=ABS(A(IP,K))*SCALES(IP)
            IF(BIG.GT.SIZE) GO TO 50
            BIG=SIZE
            IPIV=I
   50    CONTINUE
         IF(BIG.EQ.0.) GO TO 120
         IF(IPIV.EQ.K) GO TO 60
         J=IPS(K)
         IPS(K)=IPS(IPIV)
         IPS(IPIV)=J
         ODD=-ODD
   60    CONTINUE 
         KP=IPS(K) 
         PIVOT=A(KP,K)
         KP1=K+1
         DO 80 I=KP1,N
            IP=IPS(I)
            EM=A(IP,K)/PIVOT
            A(IP,K)=EM
            IF(EM.EQ.0.) GO TO 80
            DO 70 J=KP1,N
               A(IP,J)=A(IP,J)-EM*A(KP,J) 
   70       CONTINUE
   80    CONTINUE
   90 CONTINUE
      KP=IPS(N)
      IF(A(KP,N).EQ.0.) GO TO 120
C
C	CALCULATE THE DETERMINANT
C
      D=1.
      DO 100 I=1,N
         IP=IPS(I)
         D=D*A(IP,I)
  100 CONTINUE 
      D=D*ODD 
      KERR=3 
      RETURN
C
C	ERROR EXIST
C
  110 IF(KERR.EQ.0) CALL SING(1) 
      D=0.
      KERR=1
      RETURN
  120 IF(KERR.EQ.0) CALL SING(2) 
      D=0.
      KERR=2 
      RETURN 
      END
C
      SUBROUTINE SING(IWHY)
C
      GO TO (10,30),IWHY
   10 WRITE(UNIT=*, FMT=*) ' *** MATRIX WITH ZERO ROW IN DECOMPOSE'
      RETURN
   30 WRITE(UNIT=*, FMT=*) ' *** SINGULAR MATRIX IN DECOMPOSE.'
      WRITE(UNIT=*, FMT=*) '     ZERO DIVIDE IN SOLVE.'
      RETURN
      END
C
      SUBROUTINE SOLVE (NX,NN,A,B,X,IPS) 
C
C     THIS ROUTINE DOES FORWARD ELIMINATION AND BACK SUBSTITUTION
C     TO SOLVE THE EQUATION LU*X=B
C     THE PARAMETERS OF THE CALLING SEQUENCE ARE
C     N       THE ORDER OF LU
C     NX      THE NUMBER OF LINES OF THE MATRIX LU
C	      IN THE CALLING PROGRAM
C     A       THE DECOMPOSED MATRIX, L AND U ARE PACKED
C	      INTO ONE MATRIX BY DECOMP
C     B       THE RIGHT HAND SIDE OF THE EQUATION
C     IPS     VECTOR CONTAINING PERMUTATED ROW NUMBERS OF A
C	      THIS VECTOR IS BY DECOMP
C
      DIMENSION A(NX,NX), B(NX), X(NX), IPS(NX)
      N=NN
C
C     FIND THE SOLUTION IN THE CASE N=1
C
      IF(N.GT.1) GO TO 10
      X(1)=B(1)/A(1,1)
      RETURN
C
C     FORWARD ELIMINATION USING L AND STORING THE RESULT IN X
C
   10 NP1=N+1
      IP=IPS(1)
      X(1)=B(IP) 
      DO 30 I=2,N 
         IP=IPS(I)
         IM1=I-1
         SUM=0.
         DO 20 J=1,IM1
            SUM=SUM+A(IP,J)*X(J) 
   20    CONTINUE
         X(I)=B(IP)-SUM
   30 CONTINUE
C
C	BACK SUBSTITUTION USING U AND X
C
      IP=IPS(N) 
      X(N)=X(N)/A(IP,N)
      DO 50 IBACK=2,N
         I=NP1-IBACK
         IP=IPS(I)
         IP1=I+1 
         SUM=0.
         DO 40 J=IP1,N
               SUM=SUM+A(IP,J)*X(J)
   40    CONTINUE 
         X(I)=(X(I)-SUM)/A(IP,I)
   50 CONTINUE
      RETURN
      END
