C------------------------------------------------------
C
      PROGRAM VARIC
C
C**** PROGRAMME PRINCIPAL.
C**** TEST DE L'ELEMENT FINI DE DEGRE 2 A 16 PARAMETRES. 
C**** VARIANTE I-C.
C
C**** INITIALISATIONS.
      DIMENSION X(4), Y(4), Z(4)
C
      WRITE(UNIT=*, FMT=120)
      WRITE(UNIT=*, FMT=*)
     1 '---*** THESE P. 150, TESTS 8.4-1A ET 8.4-1B ***---'
      WRITE(UNIT=*, FMT=*) 'CONTINUITE DES FONCTIONS DE REFERENCE'
      WRITE(UNIT=*, FMT=120)
      CALL CONT
      DO 30 NEL=1,2
C****    CHOIX D'UN TETRAEDRE, AU HASARD.
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*) ' *** NOUVEL ELEMENT : X, Y, Z'
         DO 20 I=1,4
            CALL RANDOM(YFL) 
            X(I)=4.*YFL-2. 
            CALL RANDOM(YFL) 
            Y(I)=4.*YFL-2. 
            CALL RANDOM(YFL) 
            Z(I)=4.*YFL-2.
   20    CONTINUE
         WRITE(UNIT=*, FMT=110)  (X(I), Y(I), Z(I), I=1,4)
C****    INITIALISATIONS RELATIVES AU TETRAEDRE CHOISI.
         CALL ELEM0 (X(1),Y(1),Z(1), X(2),Y(2),Z(2), 
     1               X(3),Y(3),Z(3), X(4),Y(4),Z(4))
         CALL ELEM1
         CALL ELEM2
         CALL ELEM3
C****    VERIFICATIONS DE L'ELEMENT.
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 150 (ET P. 119-120), TEST 8.4-2 ***---'
         WRITE(UNIT=*, FMT=*) 'BASE DE LAGRANGE (HERMITE)'
         WRITE(UNIT=*, FMT=120)
         CALL LAGR
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 150 (ET P. 120), TEST 8.4-3 ***---'
         WRITE(UNIT=*, FMT=*) 'POLYNOMES DE DEGRE DEUX'
         WRITE(UNIT=*, FMT=120)
         CALL POL2
         WRITE(UNIT=*, FMT=120)
         WRITE(UNIT=*, FMT=*)
     1   '---*** THESE P. 150 (ET P. 120), TEST 8.4-4 ***---'
         WRITE(UNIT=*, FMT=*) 'DERIVEE NORMALE DE DEGRE .LE. 1'
         WRITE(UNIT=*, FMT=120)
         CALL NORLIN
   30 CONTINUE
      WRITE(UNIT=*, FMT=120)
      WRITE(UNIT=*, FMT=*)
     1 '---*** THESE P. 150 (ET P. 121), TEST 8.4-5 ***---'
      WRITE(UNIT=*, FMT=*)
     1 'MOSAIQUE DE DEUX TETRAEDRES (RACCORDEMENT)'
      WRITE(UNIT=*, FMT=120)
      CALL MOSAIQ
      WRITE(UNIT=*, FMT=120)
      WRITE(UNIT=*, FMT=*)
     1 '---*** THESE P. 150 (ET P. 121), TEST 8.4-6 ***---'
      WRITE(UNIT=*, FMT=*) 'PERMUTATION DES SOMMETS'
      WRITE(UNIT=*, FMT=120)
      CALL PERM (X, Y, Z) 
      STOP
  110 FORMAT (3F18.12) 
  120 FORMAT (/)
      END
C--------------------------------------------------------------------
C
      SUBROUTINE CONT
C
C**** NOUS VERIFIONS SI LES FONCTIONS D'INTERPOLATION DE L'ELEMENT DE 
C**** REFERENCE, AINSI QUE LEURS DERIVEES PARTIELLES, SONT CONTINUES. 
      DATA EPS/0.5E-06/
      D=3.*EPS
      M=8
      CALL TABUL(0.3, 0.3, 0.3, -1., 1.,0., D, M) 
      CALL TABUL(0.3, 0.3, 0.3, -1., 0., 1., D, M) 
      CALL TABUL(0.6, 0.15, 0.15, 0., -1., 1., D, M) 
      CALL TABUL(0.15, 0.6, 0.15, -1.,0. ,1., D, M) 
      CALL TABUL(0.15, 0.15, 0.6, -1., 1., 0., D, M)
      CALL TABUL(0.2, 0.3, 0.5, -0.05, 0.05, 0.25, 4.*D, M) 
      CALL TABUL(0.0, 0.5, 0.5, -0.25, 0.25, 0.25, 4.*D, M) 
      RETURN
      END
C--------------------------------------------------
C
      SUBROUTINE TABUL (XT,YT,ZT, TX,TY,TZ, D, ME)
C
C**** TABULATION DES FONCTIONS WT11, RT10 ET RT13, AINSI QUE DE LEURS 
C**** DERIVEES PARTIELLES, LE LONG DU SEGMENT
C****                ( (XT,YT,ZT) + F*(TX,TY,TZ) ),     -D.LE.F.LE.D 
C**** EN ME+1 POINTS EQUIDISTANTS.
      DIMENSION UT(44), T(4,3)
      X0=XT 
      Y0=YT 
      Z0=ZT 
      T1=TX
      T2=TY
      T3=TZ
      WRITE(UNIT=*, FMT=100) X0,Y0,Z0,  T1,T2,T3
      M=ME
      MP1=M+1
      HHALF= D/FLOAT(M)
C**** F(I) = (2*I-M)*D/M,   I=0(1)M. 
      DO 10 IP1=1,MP1
         I=IP1-1
         FI=HHALF*FLOAT(2*I-M)
         XI=X0+FI*T1 
         YI=Y0+FI*T2 
         ZI=Z0+FI*T3
C***     J =0.
         CALL UT44(44, XI, YI, ZI, UT)
         T(1,1)=UT(17) 
         T(1,2)=UT(29) 
         T(1,3)=UT(35)
C****    J = 1.
         CALL DUT44(44, 1., 0., 0., XI, YI, ZI, UT)
         T(2,1)=UT(17) 
         T(2,2)=UT(29) 
         T(2,3)=UT(35) 
C***     J = 2.
         CALL DUT44(44, 0., 1., 0., XI, YI, ZI, UT)
         T(3,1)=UT(17) 
         T(3,2)=UT(29) 
         T(3,3)=UT(35)
C****    J = 3.
         CALL DUT44(44, 0., 0., 1., XI, YI, ZI, UT)
         T(4,1)=UT(17) 
         T(4,2)=UT(29) 
         T(4,3)=UT(35)
         WRITE(UNIT=*, FMT=110) ( (T(JP1,K), JP1=1,4), K=1,3 )
   10 CONTINUE
      WRITE(UNIT=*, FMT=120)
      RETURN
  100 FORMAT(11H **** TABUL   , 6X, 3F16.4, 16X, 3F16.4  /)
  110 FORMAT( 3(5X,4F12.8) ) 
  120 FORMAT(/)
      END
C----------------------------------------------------
C**** TEST.
      SUBROUTINE LAGR
C**** NOUS VERIFIONS SI LA BASE D'HERMITE POSSEDE LES PARAMETRES VOULUS.
C
      COMMON /TETRAS/ X(4),Y(4),Z(4), TETR(1513), IPERM(3)
      DIMENSION E(16) 
      WRITE(UNIT=*, FMT=*) ' *** LAGR : ERR REL SUR HERM ET DHERM'
      DO 10 I=1,4
C****    J = 0.
         IND=4*I-3
         CALL HERM(X(I), Y(I), Z(I), E)
         E(IND)=E(IND)-1.0
         WRITE(UNIT=*, FMT=110)   (E(K), K=1,16)
C****    J = 1.
         IND=IND+1
         CALL DHERM(1., 0., 0., X(I), Y(I), Z(I), E)
         E(IND)=E(IND)-1.0
         WRITE(UNIT=*, FMT=110)  (E(K), K=1,16)
C****	 J = 2. 
         IND=IND+1
         CALL DHERM(0., 1., 0., X(I), Y(I), Z(I), E)
         E(IND)=E(IND)-1.0
         WRITE(UNIT=*, FMT=110) (E(K), K=1,16)
C****	 J = 3. 
         IND=IND+1
         CALL DHERM(0., 0., 1., X(I), Y(I), Z(I), E)
         E(IND)=E(IND)-1.0
         WRITE(UNIT=*, FMT=110)  (E(K), K=1,16)
         WRITE(UNIT=*, FMT=120)
   10 CONTINUE
      RETURN
  110 FORMAT(16E10.1)
  120 FORMAT(/) 
      END
C----------------------------------------------------------------
C**** TEST.
      SUBROUTINE POL2
C**** NOUS VERIFIONS SI TOUS LES POLYNOMES DE DEGRE DEUX APPARTIENNENT 
C**** A L'ESPACE ELEMENT FINI.
C
      COMMON /TETRAS/ X(4),Y(4),Z(4), TETR(1513), IPERM(3)
      DIMENSION C(44), S(44), Q(10), E(20) 
      DATA Q/10*0.0/, H/1.0/
      WRITE(UNIT=*, FMT=*) ' *** POL2 : ERR SUR FINT ET DFINT'
      DO 60 IPI=1,10
C****    CHOIX D'UN POLYNOME DE DEGRE DEUX. 
         Q(IPI)=1.0
C****    CALCUL DES PARAMETRES DE CE POLYNOME.
         DO 20 I=1,4
            IND=4*I-3
            C(IND  )=P(Q,X(I),Y(I),Z(I))
            C(IND+1)=0.5*( P(Q,X(I)+H,Y(I),Z(I))
     1                    -P(Q,X(I)-H,Y(I),Z(I)) )/H
            C(IND+2)=0.5*( P(Q,X(I),Y(I)+H,Z(I))
     1                    -P(Q,X(I),Y(I)-H,Z(I)) )/H
            C(IND+3)=0.5*( P(Q,X(I),Y(I),Z(I)+H)
     1                    -P(Q,X(I),Y(I),Z(I)-H) )/H
   20    CONTINUE
C****    COORDONNEES DE LA FONCTION D'INTERPOLATION CORRESPONDANTE. 
         CALL FE2P16(C, S)
C****    COMPARAISON EN CERTAINS POINTS.
   40       CALL RANDOM(B2)
            CALL RANDOM(B3)
            CALL RANDOM(B4)
            B1=1.-B2-B3-B4
            IF (B1.LT.0.) GO TO 40
         PX=B1*X(1)+B2*X(2)+B3*X(3)+B4*X(4) 
         PY=B1*Y(1)+B2*Y(2)+B3*Y(3)+B4*Y(4) 
         PZ=B1*Z(1)+B2*Z(2)+B3*Z(3)+B4*Z(4)
         IND=2*IPI-1 
         VALP=P(Q,PX,PY,PZ)
         E(IND)=FINT(PX,PY,PZ,S) - VALP
         IND=IND+1
         CALL RANDOM(B1)
         B1=4.*B1-2.
         CALL RANDOM(B2)
         B2=4.*B2-2.
         CALL RANDOM(B3)
         B3=4.*B3-2.
         VALP=0.5*( P(Q,PX+B1*H,PY+B2*H,PZ+B3*H)
     1             -P(Q,PX-B1*H,PY-B2*H,PZ-B3*H) )/H
         E(IND)=DFINT(B1,B2,B3,PX,PY,PZ,S) - VALP
   60 CONTINUE
      WRITE(UNIT=*, FMT=110)  (E(K), K=1,20)
      RETURN
  110 FORMAT(20E10.1)
      END
C----------------------------------------------------------------
C**** TEST.
      FUNCTION P(Q,PX,PY,PZ)
C**** VALEUR DU POLYNOME DE DEGRE 2, DE COEFFICIENTS Q(1), ..., Q(10) 
C**** AU POINT (PX,PY,PZ).
C
      DIMENSION Q(10)
      X=PX 
      Y=PY 
      Z=PZ
      P= Q(1) + Q(2)*X + Q(3)*Y + Q(4)*Z
     1   + Q(5)*X*X + Q(6)*Y*Y + Q(7)*Z*Z
     2   + Q(8)*Y*Z + Q(9)*X*Z + Q(10)*X*Y 
      RETURN
      END
C----------------------------------------------------
C**** TEST.
      SUBROUTINE NORLIN
C**** ON VERIFIE SI LES DERIVEES NORMALES SONT LINEAIRES
C**** DANS LES FACES DU TETRAEDRE.
C 
      REAL MIX,MIY,MIZ
      COMMON /TETRAS/ X(4),Y(4),Z(4), TETR(1513),IPERM(3)
      DIMENSION C(44),S(44),E(16)
      WRITE(UNIT=*, FMT=*) ' *** NORLIN : ERR REL SUR DFINT'
C**** CHOIX DES PARAMETRES DE L'ELEMENT FINI. 
      DO 10 I=1,16
         CALL RANDOM(YFL)
         C(I)=4.*YFL-2.
   10 CONTINUE
C**** CALCUL DES COORDONNEES DE LA FONCTION D'INTERPOLATION.
      CALL FE2P16(C,S)
C**** VERIFICATIONS.
      DO 40 I=1,4
         J=MOD(I,4)+1
         K=MOD(J,4)+1
         L=MOD(K,4)+1 
         MIX=(Y(K)-Y(J))*(Z(L)-Z(J))-(Z(K)-Z(J))*(Y(L)-Y(J)) 
         MIY=(Z(K)-Z(J))*(X(L)-X(J))-(X(K)-X(J))*(Z(L)-Z(J)) 
         MIZ=(X(K)-X(J))*(Y(L)-Y(J))-(Y(K)-Y(J))*(X(L)-X(J)) 
         DJ=D1NOD(MIX,MIY,MIZ,J,C)
         DK=D1NOD(MIX,MIY,MIZ,K,C) 
         DL=D1NOD(MIX,MIY,MIZ,L,C) 
         DO 30 NU=1,4
   20       CALL RANDOM(YFL)
            B2=0.667*YFL
            CALL RANDOM(YFL)
            B3=0.667*YFL
            B1=1.-B2-B3
            IF (B1.LT.0.) GO TO 20
C****       INTERPOLATION LINEAIRE DE LA DERIVEE NORMALE DANS LA FACE. 
            DERLIN=B1*DJ+B2*DK+B3*DL
C****       DERIVEE DE LA FONCTION D'INTERPOLATION. 
            PX=B1*X(J)+B2*X(K)+B3*X(L) 
            PY=B1*Y(J)+B2*Y(K)+B3*Y(L) 
            PZ=B1*Z(J)+B2*Z(K)+B3*Z(L) 
            DER=DFINT(MIX,MIY,MIZ,PX,PY,PZ,S)
C****       COMPARAISON.
            IND=4*(I-1)+NU
            E(IND)=(DER-DERLIN)/DERLIN
   30    CONTINUE
   40 CONTINUE
      WRITE(UNIT=*, FMT=50)  (E(I), I=1,16)
   50 FORMAT(16E10.1) 
      RETURN
      END
C------------------------------------------------------------
C
      SUBROUTINE MOSAIQ
C
      DIMENSION PX(5),PY(5),PZ(5), GX(5),GY(5),GZ(5), S(44),C(44), E(20) 
C**** CHOIX DU PLAN DE RACCORDEMENT.
      DO 10 M=2,4
         CALL RANDOM (YFL) 
         PX(M)=4.*YFL - 2. 
         CALL RANDOM (YFL) 
         PY(M)=4.*YFL - 2. 
         CALL RANDOM (YFL) 
         PZ(M)=4.*YFL - 2.
   10 CONTINUE
C**** VECTEUR NORMAL AU PLAN DE RACCORDEMENT.
      RX=(PY(3)-PY(2))*(PZ(4)-PZ(2)) - (PZ(3)-PZ(2))*(PY(4)-PY(2))
      RY=(PZ(3)-PZ(2))*(PX(4)-PX(2)) - (PX(3)-PX(2))*(PZ(4)-PZ(2))
      RZ=(PX(3)-PX(2))*(PY(4)-PY(2)) - (PY(3)-PY(2))*(PX(4)-PX(2))
C**** CHOIX DE POINTS DANS LE PLAN DE RACCORDEMENT.
      DO 30 M=1,5
   20    CALL RANDOM (B2) 
         CALL RANDOM (B3) 
         B4=1.-B2-B3
         IF (B4.LE.0.) GO TO 20
         GX(M)=B2*PX(2) + B3*PX(3) + B4*PX(4)
         GY(M)=B2*PY(2) + B3*PY(3) + B4*PY(4)
         GZ(M)=B2*PZ(2) + B3*PZ(3) + B4*PZ(4)
   30 CONTINUE
C**** CHOIX DES PARAMETRES ATTACHES AU PLAN DE RACCORDEMENT.
      DO 40 M=5,16
         CALL RANDOM (YFL) 
         C(M)=4.*YFL - 2. 
   40 CONTINUE
C**** TEST DES DEUX TETRAEDRES M=1 ET M=5.
      DO 200 M=1,5,4
   50    CALL RANDOM (YFL) 
         PX(M)=4.*YFL - 2. 
         CALL RANDOM (YFL) 
         PY(M)=4.*YFL - 2. 
         CALL RANDOM (YFL) 
         PZ(M)=4.*YFL - 2.
         SC= RX*(PX(M) - (PX(2)+PX(3)+PX(4))/3. )
     1      +RY*(PY(M) - (PY(2)+PY(3)+PY(4))/3. )
     2      +RZ*(PZ(M) - (PZ(2)+PZ(3)+PZ(4))/3. )
         IF (M.EQ.1  .AND. SC.LT.0.1) GO TO 50
         IF (M.EQ.5  .AND. SC.GT.-0.1) GO TO 50
         WRITE(UNIT=*, FMT=*) ' *** MOSAIQUE : FINT, DFINT'
         CALL ELEM0 (PX(M),PY(M),PZ(M), PX(2),PY(2),PZ(2), 
     1               PX(3),PY(3),PZ(3), PX(4),PY(4),PZ(4)) 
         CALL ELEM1
         CALL ELEM2
         DO 70 N=1,4
            CALL RANDOM (YFL)
            C(N)=4.*YFL - 2.
   70    CONTINUE
         CALL FE2P16 (C, S)
         DO 80 N=1,5
            E(N   )= FINT (GX(N),GY(N),GZ(N), S)
            E(N+ 5)=DFINT (1., 0., 0., GX(N),GY(N),GZ(N), S) 
            E(N+10)=DFINT (0., 1., 0., GX(N),GY(N),GZ(N), S)
            E(N+15)=DFINT (0., 0., 1., GX(N),GY(N),GZ(N), S) 
   80    CONTINUE
         WRITE(UNIT=*, FMT=350)  (E(N), N=1,20)
  200 CONTINUE
      RETURN
  350 FORMAT (20F20.12)
      END
C---------------------------------------------------
C
      SUBROUTINE PERM (AX, AY, AZ)
C
      DIMENSION AX(4),AY(4),AZ(4), S(44), F0(4),FX(4),FY(4),FZ(4), 
     1          GX(15),GY(15),GZ(15), E(30), C(44) 
C**** CHOIX DES PARAMETRES.
      DO 10 M=1,4
         CALL RANDOM (YFL)
         F0(M)=2.*YFL - 1.
         CALL RANDOM (YFL)
         FX(M)=2.*YFL - 1.
         CALL RANDOM (YFL)
         FY(M)=2.*YFL - 1.
         CALL RANDOM (YFL)
         FZ(M)=2.*YFL - 1.
   10 CONTINUE
C**** CHOIX DE POINTS DANS LES FACES.
C**** FACE 1.
      DO 30 M=1,3
   20    CALL RANDOM (B2)
         CALL RANDOM (B3)
         B4=1.-B2-B3
         IF (B4.LE.0.) GO TO 20
         GX(M)=B2*AX(2) + B3*AX(3) + B4*AX(4) 
         GY(M)=B2*AY(2) + B3*AY(3) + B4*AY(4) 
         GZ(M)=B2*AZ(2) + B3*AZ(3) + B4*AZ(4)
   30 CONTINUE
C**** FACE 2.
      DO 50 M=4,6
   40    CALL RANDOM (B3)
         CALL RANDOM (B4)
         B1=1.-B3-B4
         IF (B1.LE.0.) GO TO 40
         GX(M)=B3*AX(3) + B4*AX(4) + B1*AX(1) 
         GY(M)=B3*AY(3) + B4*AY(4) + B1*AY(1) 
         GZ(M)=B3*AZ(3) + B4*AZ(4) + B1*AZ(1)
   50 CONTINUE
C**** FACE 3.
      DO 70 M=7,9
   60    CALL RANDOM (B4)
         CALL RANDOM (B1)
         B2=1.-B4-B1
         IF (B2.LE.0.) GO TO 60
         GX(M)=B4*AX(4) + B1*AX(1) + B2*AX(2) 
         GY(M)=B4*AY(4) + B1*AY(1) + B2*AY(2) 
         GZ(M)=B4*AZ(4) + B1*AZ(1) + B2*AZ(2)
   70 CONTINUE
C**** FACE 4.
      DO 90 M=10,12
   80    CALL RANDOM (B1)
         CALL RANDOM (B2)
         B3=1.-B1-B2
         IF (B3.LE.0.) GO TO 80
         GX(M)=B1*AX(1) + B2*AX(2) + B3*AX(3) 
         GY(M)=B1*AY(1) + B2*AY(2) + B3*AY(3) 
         GZ(M)=B1*AZ(1) + B2*AZ(2) + B3*AZ(3)
   90 CONTINUE
C**** CHOIX DE POINTS A L'INTERIEUR.
      DO 110 M=13,15
  100    CALL RANDOM (B1) 
         CALL RANDOM (B2) 
         CALL RANDOM (B3) 
         B4=1.-B1-B2-B3
         IF (B4.LE.0.) GO TO 100
         GX(M)=B1*AX(1) + B2*AX(2) + B3*AX(3) + B4*AX(4) 
         GY(M)=B1*AY(1) + B2*AY(2) + B3*AY(3) + B4*AY(4) 
         GZ(M)=B1*AZ(1) + B2*AZ(2) + B3*AZ(3) + B4*AZ(4)
  110 CONTINUE
C**** TEST DES 24 PERMUTATIONS DES SOMMETS DU TETRAEDRE.
      DO 200 I=1,4
         DO 190 J=1,4
            IF (J.EQ.I) GO TO 190
            DO 180 K=1,4
               IF (K.EQ.I .OR. K.EQ.J) GO TO 180
               DO 170 L=1,4
                  IF (L.EQ.I .OR. L.EQ.J .OR. L.EQ.K) GO TO 170
                  WRITE(UNIT=*, FMT=*) ' *** PERM : ',I,J,K,L,
     1                                 ' FINT, DFINT'
                  CALL ELEM0   (AX(I),AY(I),AZ(I), AX(J),AY(J),AZ(J),
     1                          AX(K),AY(K),AZ(K), AX(L),AY(L),AZ(L)) 
                  CALL ELEM1
                  CALL ELEM2 
                  C(1)=F0(I) 
                  C(2)=FX(I)
                  C(3)=FY(I) 
                  C(4)=FZ(I) 
                  C(5)=F0(J) 
                  C(6)=FX(J) 
                  C(7)=FY(J) 
                  C(8)=FZ(J) 
                  C(9)=F0(K) 
                  C(10)=FX(K) 
                  C(11)=FY(K) 
                  C(12)=FZ(K) 
                  C(13)=F0(L) 
                  C(14)=FX(L) 
                  C(15)=FY(L) 
                  C(16)=FZ(L)
                  CALL FE2P16 (C, S)
                  DO 150 M=1,15
                     IND=2*M-1
                     E(IND  )= FINT(GX(M),GY(M),GZ(M), S) 
                     E(IND+1)=DFINT(1., -3., 2., GX(M),GY(M),GZ(M), S) 
  150             CONTINUE
                  WRITE(UNIT=*, FMT=350)   (E(M), M=1,30)
  170          CONTINUE
  180       CONTINUE
  190    CONTINUE
  200 CONTINUE
      RETURN
  350 FORMAT (30F20.12) 
      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****   VARIANTE I-C.
C****                                        MARCEL DELEZE 
C****                                        LE CHENE 21 
C****   FEVRIER 1978.                        1633 MARSENS
C****   NOVEMBRE 2008                        SWITZERLAND   
C****                           
C****	
C********************************************************
C********************************************************
C--------------------------------------------------------
C
       FUNCTION PT1 (J, WT,XT,YT,ZT)
C
C--------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE 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=1.-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
      FUNCTION WT11 (WT,XT,YT,ZT)
C 
C---------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DU POLYNOME PAR MORCEAUX WT(1,1) AU POINT (XT,YT,ZT). 
C**** WT EST INACTIF.
      X=XT
      Y=YT
      Z=ZT
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 20
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 30
C**** MORCEAU 1.
      WT11=X*( 2.*X*X/3.-X*Y-X*Z+2.*Y*Z )
      RETURN
C**** MORCEAU 2.
   20 WT11=Y*Y*( Z-Y/3. )
      RETURN
C**** MORCEAU 3.
   30 WT11=Z*Z*( Y-Z/3. )
      RETURN
      END
C-----------------------------------------------------------
C
      FUNCTION ST10 (WT,XT,YT,ZT)
C
C------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE RT(1,0), RESTREINTE AU
C**** MORCEAU 1, 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 DES CINQ ARETES DU MORCEAU 1 APPARTENANT AU BORD DU TETRAEDRE. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT
      Z=ZT 
      YMX=Y-X
      ZMX=Z-X
      IF (YMX.GT.EPS .OR. ZMX.GT.EPS) GO TO 10
C**** SUR LA DROITE X=Y=Z, ON A
      ST10=3.*X*X*W
      RETURN
   10 CONTINUE
C**** DANS LE RESTE DU MORCEAU 1, ON A
      VAR=-X*W*( 3./(1.-3.*X) + 2./(Y+ZMX) )
      VAR=-YMX*ZMX*( 6.*X-2.*W-3.+VAR )
      VAR=VAR*( ZMX/(Y*(W+YMX))+YMX/(Z*(W+ZMX)) )/(YMX+ZMX)
      VAR=W*( 2.*(YMX+ZMX)+VAR )/((1.-3.*X)*(Y+ZMX))
      ST10=X*X*W*(3.+VAR)
      RETURN
      END
C-----------------------------------------------------------------
C
      FUNCTION RT10 (WT,XT,YT,ZT)
C
C-----------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE PAR MORCEAUX RT(1,0) AU POINT 
C**** (XT,YT,ZT) DU TETRAEDRE DE REFERENCE.	WT=1.-XT-YT-ZT.
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (X.GT.EPS .AND. Y.GT.EPS .AND.
     1	  Z.GT.EPS .AND. W.GT.EPS ) GO TO 50
C**** SUR LE BORD DU TETRAEDRE, ON A
      RT10=0.0
      RETURN
   50 CONTINUE
C**** A L'INTERIEUR DU TETRAEDRE, ON A
      IF (X.GT.Y .AND. Z.GE.Y) GO TO 200
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 300
  100 CONTINUE
C**** MORCEAU 1.
      RT10=ST10(W,X,Y,Z)
      RETURN
  200 CONTINUE
C**** MORCEAU 2.
      CH=X
      X=Y
      Y=CH
      GO TO 100
  300 CONTINUE
C**** MORCEAU 3.
      CH=X
      X=Z
      Z=CH
      GO TO 100 
      END
C--------------------------------------------------------------
C
      FUNCTION ST13 (WT,XT,YT,ZT)
C
C--------------------------------------------------------------	
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE RT(1,3), RESTREINTE AU
C**** MORCEAU 1, 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 DES CINQ ARETES DU MORCEAU 1 APPARTENANT AU BORD DU TETRAEDRE. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      YMX=Y-X
      ZMX=Z-X
      IF (YMX.GT.EPS .OR. ZMX.GT.EPS) GO TO 10
C**** SUR LA DROITE X=Y=Z, ON A
      ST13=0.0
      RETURN
   10 CONTINUE
C**** DANS LE RESTE DU MORCEAU 1, ON A
      VAR=1./(1.-3.*X) - 1./(W+Y)
      VAR=-4.*X + 2.*Z + X*ZMX*VAR
      VAR=W*YMX*VAR /( Y*(YMX+ZMX)*(W+YMX)*(W+Y)*(1.-3.*X) )
      VAR=VAR + 2./( (W+X)*(W+ZMX) )
      ST13=2.*X*X*ZMX*ZMX*W*VAR
      RETURN
      END
C------------------------------------------------------------------ 
C
      FUNCTION RT13 (WT,XT,YT,ZT)
C
C------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE LA FONCTION RATIONNELLE PAR MORCEAUX RT(1,3) AU 
C**** POINT (XT,YT,ZT).	       WT=1.-XY-YT-ZT.
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 300
      IF (X.GT.EPS .AND. Y.GT.EPS .AND.
     1	  Z.GT.EPS .AND. W.GT.EPS ) GO TO 50
C**** SUR LE BORD DU TETRAEDRE, ON A
      RT13=0.0
      RETURN
   50 CONTINUE
C**** A L'INTERIEUR DU TETRAEDRE, ON A 
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 200
  100 CONTINUE 
C**** MORCEAU 1. 
      RT13=ST13(W,X,Y,Z)
      RETURN
  200 CONTINUE
C**** MORCEAU 2.
      CH=X
      X=Y
      Y=CH
      GO TO 100
  300 CONTINUE
C**** MORCEAU 3. 
      RT13=0.0 
      RETURN
      END
C-----------------------------------------------------------------
C
      SUBROUTINE UT44 (NFCT, XT,YT,ZT, UT)
C
C-----------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** VALEUR DE UT(1), ... , UT(NFCT) AU POINT (XT,YT,ZT) DU TETRAEDRE 
C**** DE REFERENCE.
C**** NFCT PEUT PRENDRE LES VALEURS 16, 28, 32 OU 44.
      DIMENSION UT(44)
      N=NFCT
      X=XT
      Y=YT
      Z=ZT
      W=1.-X-Y-Z
      DO 20 I=1,4
C****    BLOC 1.    FONCTIONS NUMERO 1 A 16.
         IND=4*I - 3
         UT(IND  )=PT1(0,W,X,Y,Z) 
         UT(IND+1)=PT1(1,W,X,Y,Z) 
         UT(IND+2)=PT1(1,W,Y,Z,X) 
         UT(IND+3)=PT1(1,W,Z,X,Y) 
         IF (N.LE.16) GO TO 10
C****    BLOC 2.    FONCTIONS NUMERO 17 A 28.
         IND=3*I+14
         UT(IND  )=WT11(W,X,Y,Z) 
         UT(IND+1)=WT11(W,Y,Z,X) 
         UT(IND+2)=WT11(W,Z,X,Y) 
         IF (N.LE.28) GO TO 10
C****    BLOC 3.    FONCTIONS NUMERO 29 A 32.
         UT(I+28 )=RT10(W,X,Y,Z) 
         IF (N.LE.32) GO TO 10
C****    BLOC 4.    FONCTIONS NUMERO 33 A 44.
         IND=3*I+30
         UT(IND  )=RT13(W,Y,Z,X) 
         UT(IND+1)=RT13(W,Z,X,Y) 
         UT(IND+2)=RT13(W,X,Y,Z)
   10    CH=W
         W=X
         X=Y
         Y=Z
         Z=CH
   20 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 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
      FUNCTION DWT11 (TX,TY,TZ, WT,XT,YT,ZT)
C
C------------------------------------------------------------------	
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DU POLYNOME PAR MORCEAUX WT(1,1) DANS LA DIRECTION 
C**** (TX,TY,TZ) AU POINT (XT,YT,ZT).
C**** WT EST INACTIF.
      X=XT
      Y=YT
      Z=ZT
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 20
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 30
C**** MORCEAU 1.
      DWT11=2.*TX*(X*(X-Y-Z)+Y*Z) + X*(TY*(2.*Z-X)+TZ*(2.*Y-X)) 
      RETURN
C**** MORCEAU 2.
   20 DWT11=Y*(TY*(2.*Z-Y) + TZ*Y)
      RETURN
C**** MORCEAU 3.
   30 DWT11=Z*(TY*Z + TZ*(2.*Y-Z))
      RETURN
      END
C------------------------------------------------------------------ 
C
      FUNCTION DST10 (TX,TY,TZ, WT,XT,YT,ZT)
C 
C------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE RT(1,0), RESTREINTE AU 
C**** MORCEAU 1, DANS LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT). 
C**** WT=1.-XT-YT-ZT.
C**** NOUS SUPPOSONS QUE LE POINT (XT,YT,ZT) NE SOIT PAS SITUE SUR
C**** UNE DES CINQ ARETES DU MORCEAU 1 APPARTENANT AU BORD DU TETRAEDRE. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT 
      YMX=Y-X
      ZMX=Z-X
      T1=TX 
      T2=TY 
      T3=TZ
      T0=-T1-T2-T3
      IF (YMX.GT.EPS .OR. ZMX.GT.EPS) GO TO 10
C**** SUR LA DROITE X=Y=Z, LA VALEUR ET LA DERIVEE DE C SONT NULLES. 
      V1=0.0
      DV1=0.0 
      GO TO 20
   10 CONTINUE
C**** DANS LE RESTE DU MORCEAU 1, ON A
C**** VALEUR ET DERIVEE DE A (PLACEES DANS V1, DV1).
      P  =X*W
      DP =T1*W + T0*X
      V1 =-P*( 3./(1.-3.*X) + 2./(Y+ZMX) )
      DV1=-3.*(DP - P*( -3.*T1/(1.-3.*X) ))/(1.-3.*X)
     1    -2.*(DP - P*( (T2+T3-T1)/(Y+ZMX) ))/(Y+ZMX)
C**** VALEUR ET DERIVEE DE B (PLACEES DANS V2,DV2).
      P  =6.*X-2.*W-3.+V1
      V2 =-YMX*ZMX*P
      DV2=-((T2-T1)*ZMX+YMX*(T3-T1))*P - YMX*ZMX*(6.*T1-2.*T0+DV1)
C**** VALEUR ET DERIVEE DE C (PLACEES DANS V1, DV1).
      P =ZMX*V2
      DP =(T3-T1)*V2 + ZMX*DV2
      Q =(YMX+ZMX)*Y*(W+YMX)
      V1 =P/Q
      DV1=(DP - P*( (T2+T3-2.*T1)/(YMX+ZMX)+T2/Y+(T0+T2-T1)/(W+YMX) ))/Q 
      P =YMX*V2
      DP =(T2-T1)*V2 + YMX*DV2
      Q  =(YMX+ZMX)*Z*(W+ZMX)
      V1 =V1 + P/Q
      DV1=DV1 + (DP - P*( (T2+T3-2.*T1)/(YMX+ZMX)+T3/Z 
     1                   +(T0+T3-T1)/(W+ZMX) ))/Q 
   20 CONTINUE
C**** VALEUR ET DERIVEE DE D (PLACEES DANS V2, DV2). 
      P =W*( 2.*(YMX+ZMX)+V1 )
      DP =T0*( 2.*(YMX+ZMX)+V1)+ W*( 2.*(T2+T3-2.*T1)+DV1 ) 
      Q =(1.-3.*X)*(Y+ZMX)
      V2 =P/Q
      DV2=(DP - P*( -3.*T1/(1.-3.*X)+(T2+T3-T1)/(Y+ZMX) ))/Q 
C**** DERIVEE DE RT(1,0).
      DST10=X*( (T1*2.*W+X*T0)*(3.+V2) + X*W*DV2 ) 
      RETURN
      END
C------------------------------------------------------------------
C
      FUNCTION DRT10 (TX,TY,TZ, WT,XT,YT,ZT)
C 
C------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE PAR MORCEAUX RT(1,0) DANS 
C**** LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT).  WT=1.-XT-YT-ZT.
      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 50
C**** DANS LES FACES 2, 3 ET 4, ON A
      DRT10=0.0
      RETURN
   50 CONTINUE
      T1=TX 
      T2=TY 
      T3=TZ
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 200
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 300
  100 CONTINUE
C**** MORCEAU 1.
      IF (W.GT.EPS) GO TO 150
C**** DANS L'INTERSECTION DU MORCEAU 1 ET DE LA FACE 1, ON A
      DRT10= (-T1-T2-T3)*3.*X*X
      RETURN 
  150 CONTINUE
C**** DANS L'INTERSECTION DU MORCEAU 1 ET DE L'INTERIEUR DU TETRAEDRE, 
C**** ON A
      DRT10=DST10(T1,T2,T3,W,X,Y,Z)
      RETURN 
  200 CONTINUE
C**** MORCEAU 2.
      CH=X
      X=Y
      Y=CH 
      CH=T1 
      T1=T2 
      T2=CH
      GO TO 100
  300 CONTINUE
C**** MORCEAU 3.
      CH=X
      X=Z
      Z=CH 
      CH=T1 
      T1=T3 
      T3=CH
      GO TO 100
      END
C------------------------------------------------------------- 
C
      FUNCTION DST13 (TX,TY,TZ, WT,XT,YT,ZT)
C
C-------------------------------------------------------------	
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE RT(1,3), RESTREINTE AU
C**** MORCEAU 1, DANS LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT).
C**** WT=1.-XT-YT-ZT.
C**** NOUS SUPPOSONS QUE LE POINT (XT,YT,ZT) NE SOIT PAS SITUE SUR
C**** UNE DES CINQ ARETES DU MORCEAU 1 APPARTENANT AU BORD DU TETRAEDRE. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT
      Z=ZT
      YMX=Y-X 
      ZMX=Z-X
      IF (YMX.GT.EPS .OR. ZMX.GT.EPS) GO TO 10
C**** SUR LA DROITE X=Y=Z, ON A
      DST13=0.0
      RETURN
   10 CONTINUE
C**** DANS LE RESTE DU MORCEAU 1, ON A
      T1=TX 
      T2=TY 
      T3=TZ 
      T0=-T1-T2-T3
C**** VALEUR ET DERIVEE DE E.
      P =X*ZMX
      DP=T1*ZMX + X*(T3-T1)
      Q =1.-3.*X
      E =P/Q
      DE=(DP - E*(-3.*T1))/Q
      Q =W+Y
      E =-4.*X+2.*Z + E - P/Q
      DE=-4.*T1+2.*T3 + DE - ( DP - P*(T0+T2)/Q )/Q
C**** VALEUR ET DERIVEE DE F.
      P =W*YMX*E
      DP=T0*YMX*E + W*(T2-T1)*E + W*YMX*DE
      Q =Y*(YMX+ZMX)*(W+YMX)*(W+Y)*(1.-3.*X)
      F =P/Q
      DF=(DP - P*( T2/Y+(T2+T3-2.*T1)/(YMX+ZMX)+(T0+T2-T1)/(W+YMX)
     1          +(T0+T2)/(W+Y)-3.*T1/(1.-3.*X) ))/Q
      Q =(W+X)*(W+ZMX)
      F =F + 2./Q
      DF=DF - 2.*( (T0+T1)/(W+X)+(T0+T3-T1)/(W+ZMX) )/Q
C**** DERIVEE DE DST13.
      DST13=2.*X*ZMX*( 2.*(T1*ZMX+X*(T3-T1))*W*F + X*ZMX*(T0*F+W*DF) ) 
      RETURN
      END
C------------------------------------------------------------------
C
      FUNCTION DRT13 (TX,TY,TZ, WT,XT,YT,ZT)
C 
C-------------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE LA FONCTION RATIONNELLE PAR MORCEAUX RT(1,3) DANS 
C**** LA DIRECTION (TX,TY,TZ) AU POINT (XT,YT,ZT).   WT=1.-XT-YT-ZT. 
      DATA EPS/0.5E-06/
      W=WT 
      X=XT 
      Y=YT 
      Z=ZT
      IF (X.GE.Z .AND. Y.GE.Z) GO TO 300
      IF (X.GT.EPS .AND. Y.GT.EPS .AND. Z.GT.EPS) GO TO 50
C**** DANS LES FACES 2, 3 ET 4, ON A
      DRT13=0.0
      RETURN
   50 CONTINUE
      T1=TX 
      T2=TY 
      T3=TZ
      IF (X.GE.Y .AND. Z.GE.Y) GO TO 200
  100 CONTINUE
C**** MORCEAU 1.
      IF (W.GT.EPS) GO TO 150
C**** DANS L'INTERSECTION DU MORCEAU 1 ET DE LA FACE 1, ON A 
      DRT13=4.*(-T1-T2-T3)*X*(Z-X)
      RETURN
  150 CONTINUE
C**** DANS L'INTERSECTION DU MORCEAU 1 ET DE L'INTERIEUR DU TETRAEDRE, 
C**** ON A
      DRT13=DST13(T1,T2,T3,W,X,Y,Z)
      RETURN
  200 CONTINUE
C**** MORCEAU 2.
      CH=X
      X=Y
      Y=CH
      CH=T1
      T1=T2
      T2=CH
      GO TO 100
  300 CONTINUE
C**** MORCEAU 3. 
      DRT13=0.0 
      RETURN
      END
C-------------------------------------------------------------
C
      SUBROUTINE DUT44 (NFCT, TX,TY,TZ, XT,YT,ZT, DUT)
C 
C-------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT DE REFERENCE.
C**** DERIVEE DE UT(1), ..., UT(NFCT) DANS LA DIRECTION (TX,TY,TZ) 
C**** AU POINT (XT,YT,ZT).
C**** NFCT PEUT PRENDRE LES VALEURS 16, 28, 32 OU 44. 
      DIMENSION DUT(44)
      N=NFCT
      T1=TX
      T2=TY
      T3=TZ
      T0=-T1-T2-T3
      X=XT
      Y=YT
      Z=ZT
      W=1.-X-Y-Z
      DO 20 I=1,4
C****    BLOC 1.   FONCTIONS NUMERO 1 A 16. 
         IND=4*I-3
         DUT(IND  )=DPT1(0,T1,T2,T3,W,X,Y,Z) 
         DUT(IND+1)=DPT1(1,T1,T2,T3,W,X,Y,Z) 
         DUT(IND+2)=DPT1(1,T2,T3,T1,W,Y,Z,X) 
         DUT(IND+3)=DPT1(1,T3,T1,T2,W,Z,X,Y) 
         IF (N.LE.16) GO TO 10
C****    BLOC 2.   FONCTIONS NUMERO 17 A 28. 
         IND=3*I+14
         DUT(IND  )=DWT11(T1,T2,T3,W,X,Y,Z) 
         DUT(IND+1)=DWT11(T2,T3,T1,W,Y,Z,X) 
         DUT(IND+2)=DWT11(T3,T1,T2,W,Z,X,Y) 
         IF (N.LE.28) GO TO 10
C****    BLOC 3.   FONCTIONS NUMERO 29 A 32.
         DUT(I+28 )=DRT10(T1,T2,T3,W,X,Y,Z) 
         IF (N.LE.32) GO TO 10
C****    BLOC 4.   FONCTIONS NUMERO 33 A 44. 
         IND=3*I+30
         DUT(IND  )=DRT13(T2,T3,T1,W,Y,Z,X) 
         DUT(IND+1)=DRT13(T3,T1,T2,W,Z,X,Y) 
         DUT(IND+2)=DRT13(T1,T2,T3,W,X,Y,Z)
   10    CH=W
         W=X 
         X=Y 
         Y=Z
         Z=CH
         CH=T0 
         T0=T1 
         T1=T2 
         T2=T3 
         T3=CH
   20 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
C****                    DU 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****             MX(I), MY(I), MZ(I), VECTEUR NORMAL A LA I-EME FACE DU
C****                    TETRAEDRE, I=1(1)4.
C****             NX(I,J), NY(I,J), NZ(I,J), J=1,2,3, VECTEURS PARAL-
C****                    LELES A LA I-EME FACE ET NORMAUX AUX ARETES DU
C****                    TETRAEDRE, I=1(1)4.
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION AUX(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)
C**** FACTORISATION DE LL.
      CALL DECOMP(3,3,LL,AUX,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 10 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1 
         L=MOD(K,4)+1
C****    INITIALISATION DU VECTEUR (MX(I), MY(I), MZ(I)). 
         MX(I)=(Y(K)-Y(J))*(Z(L)-Z(J)) - (Z(K)-Z(J))*(Y(L)-Y(J)) 
         MY(I)=(Z(K)-Z(J))*(X(L)-X(J)) - (X(K)-X(J))*(Z(L)-Z(J)) 
         MZ(I)=(X(K)-X(J))*(Y(L)-Y(J)) - (Y(K)-Y(J))*(X(L)-X(J)) 
         S=SQRT(MX(I)**2 + MY(I)**2 + MZ(I)**2)
         MX(I)=MX(I)/S
         MY(I)=MY(I)/S
         MZ(I)=MZ(I)/S
C****    VECTEURS   (NX(I,J), NY(I,J), NZ(I,J)),  J=1,2,3. 
         NX(I,1)=MY(I)*(Z(L)-Z(K)) - MZ(I)*(Y(L)-Y(K)) 
         NY(I,1)=MZ(I)*(X(L)-X(K)) - MX(I)*(Z(L)-Z(K)) 
         NZ(I,1)=MX(I)*(Y(L)-Y(K)) - MY(I)*(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)=MY(I)*(Z(J)-Z(L)) - MZ(I)*(Y(J)-Y(L)) 
         NY(I,2)=MZ(I)*(X(J)-X(L)) - MX(I)*(Z(J)-Z(L)) 
         NZ(I,2)=MX(I)*(Y(J)-Y(L)) - MY(I)*(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)=MY(I)*(Z(K)-Z(J)) - MZ(I)*(Y(K)-Y(J)) 
         NY(I,3)=MZ(I)*(X(K)-X(J)) - MX(I)*(Z(K)-Z(J)) 
         NZ(I,3)=MX(I)*(Y(K)-Y(J)) - MY(I)*(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
   10 CONTINUE
      RETURN
      END
C----------------------------------------------------------------
C
      SUBROUTINE U44 (NFCT, PX,PY,PZ, U)
C 
C----------------------------------------------------------------
C**** BASE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****                    INITIALISE A L'AIDE DE LA SUBROUTINE ELEM0.
C****             NFCT   PEUT PRENDRE LES VALEURS 16, 28, 32 OU 44.
C****             (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  U(1),	..., U(NFCT) VALEURS DES NFCT PREMIERES FONC-
C****                    TIONS DE BASE AU POINT (PX, PY, PZ). 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION U(44), 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 U44 :',
     2                ' THE POINT IS OUTSIDE THE ELEMENT.'
      NERR=NERR+1
   20 CONTINUE
C**** VALEURS DES FONCTIONS DE BASE.
      CALL UT44(NFCT,XT,YT,ZT,U)
      RETURN
      END
C------------------------------------------------------------
C
      SUBROUTINE DU44 (NFCT, VX,VY,VZ, PX,PY,PZ, DU) 
C
C------------------------------------------------------------ 	
C**** BASE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI.
C**** ENTREES  :  TETRAS BLOC COMMUN QUI DOIT AVOIR ETE PREALABLEMENT 
C****                    INITIALISE A L'AIDE DE LA SUBROUTINE ELEM0.
C****             NFCT   PEUT PRENDRE LES VALEURS 16, 28, 32 OU 44.
C****	          (VX, VY, VZ) UN VECTEUR.
C****             (PX, PY, PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  DU(1), ..., DU(NFCT) DERIVEES DES NFCT PREMIERES
C****                    FONCTIONS DE BASE DANS LA DIRECTION (VX,VY,VZ)
C****                    AU POINT (PX, PY, PZ). 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION DU(44), 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 DU44 :',
     2                ' THE POINT IS OUTSIDE THE ELEMENT.'
      NERR=NERR+1
   20 CONTINUE
C**** DERIVEES DES FONCTIONS DE BASE.
      CALL DUT44(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, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,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 ET 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))
C****    INITIALISATION DE D6.
         D6(I)=MX(I)*(3.*X(I)-X(J)-X(K)-X(L))
     1        +MY(I)*(3.*Y(I)-Y(J)-Y(K)-Y(L))
     2        +MZ(I)*(3.*Z(I)-Z(J)-Z(K)-Z(L))
C****	 INITIALISATION DE D7.
         IND=3*I-2
         D7(IND  )=D6(I)
         D7(IND+1)=D6(I) 
         D7(IND+2)=D6(I) 
   10 CONTINUE
      RETURN
      END
C-------------------------------------------------------------
C
      SUBROUTINE ELEM2
C
C-------------------------------------------------------------	
C**** MATRICE D'INTERPOLATION DE L'ELEMENT GENERIQUE ELARGI. 
C**** INITIALISATION DES BLOCS SOUS-DIAGONAUX DE LA MATRICE Q. 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1	     NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION DU(44)
      DO 80 I=1,4
         J=MOD(I,4)+1
         K=MOD(J,4)+1
         L=MOD(K,4)+1
C****    INITIALISATION DU BLOC Q1.
         PX=0.5*(X(K)+X(L)) 
         PY=0.5*(Y(K)+Y(L)) 
         PZ=0.5*(Z(K)+Z(L))
         CALL DU44(16,NX(I,1),NY(I,1),NZ(I,1),PX,PY,PZ,DU)
         IND=3*I-2
         DO 10 M=1,16
            Q1(IND,M)=DU(M)
   10    CONTINUE
         PX=0.5*(X(J)+X(L)) 
         PY=0.5*(Y(J)+Y(L)) 
         PZ=0.5*(Z(J)+Z(L))
         CALL DU44(16,NX(I,2),NY(I,2),NZ(I,2),PX,PY,PZ,DU)
         IND=IND+1
         DO 20 M=1,16
            Q1(IND,M)=DU(M)
   20    CONTINUE
         PX=0.5*(X(K)+X(J)) 
         PY=0.5*(Y(K)+Y(J)) 
         PZ=0.5*(Z(K)+Z(J))
         CALL DU44(16,NX(I,3),NY(I,3),NZ(I,3),PX,PY,PZ,DU)
         IND=IND+1
         DO 30 M=1,16
            Q1(IND,M)=DU(M)
   30    CONTINUE
C****    INITIALISATION DU BLOC Q2. 
         PX=(X(J)+X(K)+X(L))/3. 
         PY=(Y(J)+Y(K)+Y(L))/3. 
         PZ=(Z(J)+Z(K)+Z(L))/3.
         CALL DU44(28,MX(I),MY(I),MZ(I),PX,PY,PZ,DU)
         DO 40 M=1,28
            Q2(I,M)=DU(M)
   40    CONTINUE
C****    INITIALISATION DU BLOC Q3. 
         PX=(4.*X(J)+X(K)+X(L))/6. 
         PY=(4.*Y(J)+Y(K)+Y(L))/6. 
         PZ=(4.*Z(J)+Z(K)+Z(L))/6.
         CALL DU44(32,MX(I),MY(I),MZ(I),PX,PY,PZ,DU)
         IND=3*I-2
         DO 50 M=1,32
            Q3(IND,M)=DU(M)
   50    CONTINUE
         PX=(X(J)+4.*X(K)+X(L))/6. 
         PY=(Y(J)+4.*Y(K)+Y(L))/6.
         PZ=(Z(J)+4.*Z(K)+Z(L))/6.
         CALL DU44(32,MX(I),MY(I),MZ(I),PX,PY,PZ,DU) 
         IND=IND+1
         DO 60 M=1,32
            Q3(IND,M)=DU(M) 
   60    CONTINUE
         PX=(X(J)+X(K)+4.*X(L))/6.
         PY=(Y(J)+Y(K)+4.*Y(L))/6. 
         PZ=(Z(J)+Z(K)+4.*Z(L))/6.
         CALL DU44(32,MX(I),MY(I),MZ(I),PX,PY,PZ,DU) 
         IND=IND+1
         DO 70 M=1,32
            Q3(IND,M)=DU(M) 
   70    CONTINUE
   80 CONTINUE
      RETURN
      END
C--------------------------------------------------------------
C
      SUBROUTINE FE2P44 (C, S)
C
C--------------------------------------------------------------
C**** ELEMENT GENERIQUE ELARGI.
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(44) VALEURS DES PARAMETRES DE L'ELEMENT.
C**** SORTIE   :  S(1), ..., S(44) 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, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION C(44), S(44)
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
C**** CALCUL DE S(29), ..., S(32).
      DO 50 I=1,4
      SI=0.0
         DO 40 J=1,28 
            SI=SI+Q2(I,J)*S(J)
   40    CONTINUE
         K=28+I
         S(K)=D6(I)*(C(K)-SI)
   50 CONTINUE
C**** CALCUL DE S(33),	..., S(44).
      DO 70 I=1,12
         SI=0.0
         DO 60 J=1,32 
            SI=SI+Q3(I,J)*S(J)
   60    CONTINUE
         K=I+32 
         S(K)=D7(I)*(C(K)-SI) 
   70 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(44)
      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**** SORTIE   :  S(1), ..., S(44) 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(44) SONT DES PLACES DE TRAVAIL. 
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION C(44), S(44)
      DO 10 I=1,4
         J=MOD(I,4)+1 
         K=MOD(J,4)+1
         L=MOD(K,4)+1
C****    CALCUL DES PARAMETRES C(17) A C(28).
         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) ) 
C****    CALCUL DES PARAMETRES C(29) A C(32). 
         DJ=D1NOD(MX(I),MY(I),MZ(I),J,C)
         DK=D1NOD(MX(I),MY(I),MZ(I),K,C)
         DL=D1NOD(MX(I),MY(I),MZ(I),L,C)
         C(I+28)=(DJ+DK+DL)/3.
C****    CALCUL DES PARAMETRES C(33) A C(44).
         IND=3*I+30
         C(IND  )=(4.*DJ+DK+DL)/6.
         C(IND+1)=(DJ+4.*DK+DL)/6.
         C(IND+2)=(DJ+DK+4.*DL)/6.
   10 CONTINUE
C**** CALCUL DES COORDONNEES S(1), ... , S(44).
      CALL FE2P44(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(44) 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(44), S(44)
      CALL U44(44,PX,PY,PZ,U)
      V=0.0
      DO 10 I=1,44
         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(44) 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(44), S(44)
      CALL DU44(44,VX,VY,VZ,PX,PY,PZ,DU)
      DV=0.0
      DO 10 I=1,44
         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, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3)
      DIMENSION C(44), S(44) 
      DATA C/44*0.0/
      DO 20 I=1,16
C****    CALCUL DES COORDONNEES T(1,I), ..., T(44,I) DU I-EME ELEMENT
C****	 DE LA BASE D'HERMITE PAR RAPPORT A LA BASE D'INTERPOLATION
C****	 DE L'ELEMENT GENERIQUE ELARGI.
         C(I)=1.0
         CALL FE2P16(C,S)
         DO 10 J=1,44
            T(J,I)=S(J)
   10    CONTINUE 
         C(I)=0.0 
   20 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**** SORTIE   :  V(1), ..., V(16) VALEURS DES 16 FONCTIONS DE LA BASE
C****                    D'HERMITE AU POINT (PX, PY, PZ).
      REAL LL, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3) 
      DIMENSION V(16), U(44)
      CALL U44(44,PX,PY,PZ, U)
      DO 20 I=1,16
         SI=0.0
         DO 10 J=1,44
            SI=SI + T(J,I)*U(J)
   10    CONTINUE
         V(I)=SI
   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**** 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****             (VX, VY, VZ) UN VECTEUR.
C****             (PX, PY,PZ) UN POINT DU TETRAEDRE GENERIQUE.
C**** SORTIE   :  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, MX,MY,MZ, NX,NY,NZ
      COMMON /TETRAS/ X(4),Y(4),Z(4),
     1       NX(4,3),NY(4,3),NZ(4,3), MX(4),MY(4),MZ(4),
     2       D(4,3,3),D5(12),D6(4),D7(12), Q1(12,16),Q2(4,28),Q3(12,32),
     3       T(44,16), LL(3,3),IPERM(3) 
      DIMENSION DV(16), DU(44)
      CALL DU44(44,VX,VY,VZ,PX,PY,PZ,DU) 
      DO 20 I=1,16
         SI=0.0
         DO 10 J=1,44
            SI=SI + T(J,I)*DU(J)
   10    CONTINUE 
         DV(I)=SI 
   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
C------------------------------------------
      SUBROUTINE RANDOM(ALEAT)
C
C**** GENERATEUR DE NOMBRES ALEATOIRES UNIFORMEMENT DISTRIBUES
C**** ENTRE 0.0 ET 1.0
      DATA IA,IX,IC/6133,4409,2741/
      DATA M,RM/81799,81799./
      IX=MOD(IA*IX+IC,M)
      X=IX
      ALEAT=X/RM
      RETURN
      END
