C QUADRI    SOURCE    CHAT      05/01/13    02:40:46     5004           
      SUBROUTINE QUADRI(ICM,IDIM,XE,XEC,XX,SHPC,SGA,QSI,IG)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)  
      DIMENSION A(2,8),AJ(2,2),X(2),XE(3,*),XEC(3,2),XX(3,2),
     & BJ(2,2),F(2),SHPC(6,20,*),SGA(3,5,5),QSI(3,2),DT(2),SIB(2)
      DATA EPSI/1.D-7/
C --- IG: NUMERO DU POINT DE GAUSS
      DO 10 I=1,2
      A(I,1)=-XE(I,1)-XE(I,3)-XE(I,5)-XE(I,7)
     &  +2.D0*(XE(I,2)+XE(I,4)+XE(I,6)+XE(I,8))-4.D0*XX(I,IG)
      A(I,2)= 2.D0*( XE(I,4)-XE(I,8))
      A(I,3)= 2.D0*(-XE(I,2)+XE(I,6))
      A(I,4)= XE(I,1)-2.D0*XE(I,2)+XE(I,3)+XE(I,5)-2.D0*XE(I,6)+XE(I,7)
      A(I,5)= XE(I,1)-XE(I,3)+XE(I,5)-XE(I,7)
      A(I,6)= XE(I,1)+XE(I,3)-2.D0*XE(I,4)+XE(I,5)+XE(I,7)-2.D0*XE(I,8)
      A(I,7)=-XE(I,1)+2.D0*XE(I,2)-XE(I,3)+XE(I,5)-2.D0*XE(I,6)+XE(I,7)
  10  A(I,8)=-XE(I,1)+XE(I,3)-2.D0*XE(I,4)+XE(I,5)-XE(I,7)+2.D0*XE(I,8)
C --- Initialisation  
      DO 15 I=1,2
  15  X(I)=0.D0
C
      ITER=0
      ICRI=0
 100  CONTINUE
C --- LES DERIVEES DES FONCTIONS DE TRANSFORMATION
      QP=X(1)*X(2)
      Q1=X(1)*X(1)
      Q2=X(2)*X(2)
      DO 20 I=1,IDIM
      AJ(I,1)=A(I,2)+2.D0*A(I,4)*X(1)+A(I,5)*X(2)
     &        +2.D0*A(I,7)*QP+A(I,8)*Q2
  20  AJ(I,2)=A(I,3)+A(I,5)*X(1)+2.D0*A(I,6)*X(2)
     &        +A(I,7)*Q1+2.D0*A(I,8)*QP
C --- Calcul de l'inverse de la matrice AJ et son determinant :BJ=inverse(AJ)
      CALL INVER2(AJ,BJ,DET,ICRI)
      DO 40 I=1,IDIM
  40  F(I)=A(I,1)+A(I,2)*X(1)+A(I,3)*X(2)+A(I,4)*Q1+A(I,5)*QP
     &    +A(I,6)*Q2+A(I,7)*Q1*X(2)+A(I,8)*Q2*X(1)
      DD=SQRT(F(1)*F(1)+F(2)*F(2))
      IF(DD.GT.EPSI) THEN
      ITER=ITER+1
      DO 50 I=1,IDIM
  50  X(I)=X(I)-BJ(I,1)*F(1)-BJ(I,2)*F(2)
      GOTO 100
      END IF
      DO 60 I=1,IDIM
  60  QSI(I,IG)=X(I)
      C=0.25D0
      D=0.5D0
      YK1=1.D0+X(1)
      YK2=1.D0-X(1)
      YN1=1.D0+X(2)
      YN2=1.D0-X(2)
      YKN=X(1)+X(2)
      YNK=X(1)-X(2)
C --- Calcul des fonctions d'interpolation leur derivees
C --- Fonctions d'interpolation Ni
      SHPC(1,1,IG)=-C*YK2*YN2*(YK1+X(2))
      SHPC(1,2,IG)= D*YK1*YK2*YN2
      SHPC(1,3,IG)=-C*YK1*YN2*(YK2+X(2))
      SHPC(1,4,IG)= D*YK1*YN1*YN2
      SHPC(1,5,IG)=-C*YK1*YN1*(YK2-X(2))
      SHPC(1,6,IG)= D*YK1*YK2*YN1
      SHPC(1,7,IG)=-C*YK2*YN1*(YK1-X(2))
      SHPC(1,8,IG)= D*YK2*YN1*YN2
      IF(ICM.EQ.0) GOTO 1000
C --- Derivee de Ni par rapport a QSI
      SHPC(2,1,IG)= C*YN2*(X(1)+YKN)
      SHPC(2,2,IG)=-YN2*X(1)
      SHPC(2,3,IG)= C*YN2*(X(1)+YNK)
      SHPC(2,4,IG)= D*YN1*YN2
      SHPC(2,5,IG)= C*YN1*(X(1)+YKN)
      SHPC(2,6,IG)=-YN1*X(1)
      SHPC(2,7,IG)= C*YN1*(X(1)+YNK)
      SHPC(2,8,IG)=-SHPC(2,4,IG)
C --- Derivee de Ni par rapport a ETA
      SHPC(3,1,IG)= C*YK2*(YKN+X(2))
      SHPC(3,2,IG)=-D*YK1*YK2
      SHPC(3,3,IG)=-C*YK1*(YNK-X(2))
      SHPC(3,4,IG)=-YK1*X(2)
      SHPC(3,5,IG)= C*YK1*(YKN+X(2))
      SHPC(3,6,IG)=-SHPC(3,2,IG)
      SHPC(3,7,IG)=-C*YK2*(YNK-X(2))
      SHPC(3,8,IG)=-YK2*X(2)
C
      DJ=0.D0
      DO 70 I=1,IDIM
      DT(I)=0.50D0*(XEC(I,2)-XEC(I,1))
      SIB(I)=DT(I)
      DJ=DJ+DT(I)*DT(I)
      DO 70 J=1,IDIM
  70  BJ(I,J)=0.D0
      DO 90 I=1,IDIM
      DO 90 J=1,IDIM
      DO 80 K=1,8
  80  BJ(I,J)=BJ(I,J)+XE(I,K)*SHPC(J+1,K,IG)
  90  CONTINUE

C --- LES DRIVEES
      CALL LINEQ1(BJ,DT,X)
      DJ=SQRT(DJ)
      DO 120 I=1,IDIM
 120  SIB(I)=SIB(I)/DJ
1000  CONTINUE
      DO 2000 IS=1,IDIM
C --- Coordonnees du point d'integration dans le systeme reel (x,y,z)
      SGA(IS,1,IG)=XX(IS,IG)
C --- Coordonnees du point d'integration dans le systeme de reference 
C     (qsi,eta,dze)
      SGA(IS,2,IG)=QSI(IS,IG)
C --- Derivee de (qsi,eta) par rapport a miu
      SGA(IS,3,IG)=X(IS)
C --- Cosinus du cable avec les axes
      SGA(IS,4,IG)=SIB(IS)
C --- Jacobien
      SGA(IS,5,IG)=DJ
 2000 CONTINUE
      END