Sources Esope | Cast3M

bbcalc
C BBCALC    SOURCE    OF166741  23/04/25    21:15:05     11608          
 
      SUBROUTINE BBCALC(MELE,NBNN,MFR,IDIM,XE,
     1                  NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
     2                  NSTRS,LRE,IFOUR,NHRM,A,BB,
     3                  SHPTOT,SHP,BGENE,XDPGE,YDPGE,P)
 
*******************************************************************
* projection de dN/dx dans l espace de dimension inferieure
c=======================================================================
c   calcul des composantes de la matrice b-barre-dlatation
c           - en 2d: elements ict3,icq4,ict6
c           - en 3d: elements ....
c   calcul des coefficients de modification de la matrice b-barre-dilatation
c           - en 2d: elements icq8
c           - en 3d: elements ....
c                    coefficients {a}={a1,a2,a3}
c=======================================================================
c   input
c     mele = numero de l'element dans nomtp
c     nbnn = nombre de noeuds
c     mfr  = formulation de l'element
c     idim = dimension espace
c     xe = coordonnees de l'element
c     nbpgau = nombre de points de gauss pour la rigidite
c     poigau = poids des points de gauss
c     qsigau,etagau,dzegau = coordonnees des points de gauss
c     nstrs = nombre de composantes de contraintes
c     lre = nombre de colonnes de la matrice b
c     ifour = ifour dans ccoptio
c     nhrm = numero du mode de fourier
c     xdpge,ydpge = coordonnees du point autour duquel se fait
c                   le mouvement de la section
c
c   output
c      tableau des coefficients de modification: a(4,*)
c      tableau des composantes bbarre-dilatation: bb(3,*)
c=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC CCREEL
      PARAMETER (XZer=0.D0, XUn=1.D0, X1s2=0.5D0, X1s4 = 0.25D0,
     &           X2s3 = 0.666666666666666666666666666666666666666667D0,
     &           X1s3 = 0.333333333333333333333333333333333333333333D0,
     &           X1s6 = 0.166666666666666666666666666666666666666667D0)
 
      DIMENSION XE(3,*), POIGAU(*),QSIGAU(*),ETAGAU(*),DZEGAU(*),
     &          SHP(6,*),SHPTOT(6,NBNN,*),
     &          A(4,*),BB(3,*),BGENE(NSTRS,*),P(4,*)
 
      DIMENSION XNUM1(8),XNUM2(8),XNUM3(8)
 
*      write(*,*) 'Entree dans le sousprogramme BBCALC'
*      write(*,*) 'Element :',MELE
*      write(*,*) 'Nombre de points Gauss : ', nbpgau
 
      IF (MFR.NE.31) GOTO 666
C=======================================================================
C= Elements incompressibles implementes :
C=======================================================================
C= NOM  :   ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15,
C= MELE :    69 ,  70 ,  71 ,  72 ,  73 ,  74 ,  75 ,  76 ,  77 ,  78 ,
C= NOM  :   ICY5, IC13
C= MELE :    273,  274
C=======================================================================
 
C= 0 = Quelques initialisations ========================================
      IFR = IFOUR+4
      IF (IFR.LE.4 .OR. IFR.GE.7) THEN
        KINC = 2
      ELSE
        KINC = 3
      ENDIF
 
      CALL ZERO(BB,3,60)
      IF (MELE.EQ.72 .OR. MELE.GE.76) THEN
        CALL ZERO(A,4,60)
        CALL ZERO(P,4,4)
      ENDIF
 
C= 1 = Branchement selon l'Element Fini BBAR ============================
      IF (MELE.EQ.273 .OR. MELE.EQ.274) GOTO 100
      GOTO ( 10, 20, 30, 40,100,100,100,100,100,100), (MELE-68)
      GOTO 999
 
  3   FORMAT (4(A,I1,A,I1,A,F10.4,2X))
  4   FORMAT (3(A,I1,A,I1,A,F10.4,2X))
 
C=======================================================================
C========== Elements MASSIFS INCOMPRESSIBLES BIDIMENSIONNELS ===========
C=======================================================================
 
*-----------------------------------------------------------------------
*----------------------------- Element ICT3 ----------------------------
*-----------------------------------------------------------------------
 10   CONTINUE
*      write(6,*) 'Element ICT3',ifour
*      write(6,*) 'Ecriture de la matrice xe[ ]'
*      write(6,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,3)
      GOTO (999,101,101,103),IFR
      GOTO 999
 
*--------contraintes planes ou deformations planes (ifour=-2,-1)
 101  CONTINUE
*--------donnee des composantes de bb-dil
      r_z = - XE(1,1)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,2)*XE(2,3)
     &      + XE(1,1)*XE(2,2) + XE(1,3)*XE(2,1) - XE(1,3)*XE(2,2)
*      AIRE = X1s2 * r_z
      r_z = XUn / r_z
      BB(1,1) = (XE(2,2)-XE(2,3)) * r_z
      BB(2,2) = (XE(1,3)-XE(1,2)) * r_z
      BB(1,3) = (XE(2,3)-XE(2,1)) * r_z
      BB(2,4) = (XE(1,1)-XE(1,3)) * r_z
      BB(1,5) = (XE(2,1)-XE(2,2)) * r_z
      BB(2,6) = (XE(1,2)-XE(1,1)) * r_z
*      write(6,*) 'calcul de aire pour bbdil'
*      write(6,*) 'aire = ',AIRE
*      write(6,*) 'ecriture de la matrice bb'
*      write(6,3) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
      GOTO 666
*--------axisymetrique (ifour=0)
 103  CONTINUE
*--------donnee des composantes de bb-dil
      r_z = - XE(1,1)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,2)*XE(2,3)
     &      + XE(1,1)*XE(2,2) + XE(1,3)*XE(2,1) - XE(1,3)*XE(2,2)
*      AIRE = X1s6 * r_z * (XE(1,1)+XE(1,2)+XE(1,3))
      r_z = XUn / r_z
      BB(1,1) = r_z * (XE(2,2)-XE(2,3))
      BB(2,2) = r_z * (XE(1,3)-XE(1,2))
      BB(1,3) = r_z * (XE(2,3)-XE(2,1))
      BB(2,4) = r_z * (XE(1,1)-XE(1,3))
      BB(1,5) = r_z * (XE(2,1)-XE(2,2))
      BB(2,6) = r_z * (XE(1,2)-XE(1,1))
*-- Integration analytique des termes BB(3,1),BB(3,3),BB(3,5)
      BB(3,1) = XUn / (XE(1,1)+XE(1,2)+XE(1,3))
      BB(3,3) = BB(3,1)
      BB(3,5) = BB(3,1)
*      write(6,*) 'calcul de aire pour bbdil'
*      write(6,*) 'aire = ',AIRE
*      write(6,*) 'ecriture de la matrice bb'
*      write(6,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
      GOTO 666
 
*-----------------------------------------------------------------------
*----------------------------- Element ICQ4 ----------------------------
*-----------------------------------------------------------------------
 20   CONTINUE
*      write(6,*) 'Element ICQ4',ifour
*      write(6,*) 'Ecriture de la matrice xe[ ]'
*      write(6,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,4)
      GOTO (999,201,201,203),IFR
      GOTO 999
*--------contraintes planes ou deformations planes (ifour=-2,-1)
 201  CONTINUE
*--------donnee des composantes de bb-dil
      r_z = - XE(1,4)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,1)*XE(2,2)
     &      - XE(1,1)*XE(2,4) - XE(1,3)*XE(2,2) + XE(1,2)*XE(2,3)
     &      + XE(1,3)*XE(2,4) + XE(1,4)*XE(2,1)
*      AIRE = X1s2 * r_z
      r_z = XUn / r_z
      BB(1,1) = r_z * (XE(2,2)-XE(2,4))
      BB(2,2) = r_z * (XE(1,4)-XE(1,2))
      BB(1,3) = r_z * (XE(2,3)-XE(2,1))
      BB(2,4) = r_z * (XE(1,1)-XE(1,3))
      BB(1,5) = r_z * (XE(2,4)-XE(2,2))
      BB(2,6) = r_z * (XE(1,2)-XE(1,4))
      BB(1,7) = r_z * (XE(2,1)-XE(2,3))
      BB(2,8) = r_z * (XE(1,3)-XE(1,1))
*      write (*,*) 'calcul de aire pour bbdil'
*      write (*,*) 'aire = ',aire
*      write (*,*) 'ecriture de la matrice bb'
*      write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
      GOTO 666
*
*--------axisymetriques(ifour=0)
 203  CONTINUE
      r_z =  XE(1,1)*( XE(1,1)*XE(2,2)-XE(1,1)*XE(2,4)-XE(1,4)*XE(2,4)
     &                +XE(1,4)*XE(2,1)-XE(1,2)*XE(2,1)+XE(1,2)*XE(2,2))
     &     + XE(1,2)*(-XE(1,2)*XE(2,1)+XE(1,2)*XE(2,3)+XE(1,3)*XE(2,3)
     &                -XE(1,3)*XE(2,2))
     &     + XE(1,3)*( XE(1,3)*XE(2,4)+XE(1,4)*XE(2,4)-XE(1,4)*XE(2,3)
     &                -XE(1,3)*XE(2,2))
     &     + XE(1,4)*XE(1,4)*(XE(2,1)-XE(2,3))
*      AIRE = x1s6 * r_z
      r_z = xUn / r_z
      BB(1,1) = ( XE(2,2)*(XE(1,2)+XE(1,1))-XE(2,4)*(XE(1,1)+XE(1,4))
     &           +X1s2*( XE(1,3)*(XE(2,2)-XE(2,4))
     &                  +XE(1,4)*(XE(2,2)+XE(2,3))
     &                  -XE(1,2)*(XE(2,3)+XE(2,4))) ) * r_z
      BB(2,2) = ( XE(1,4)*(XE(1,1)+XE(1,4))
     &           -XE(1,2)*(XE(1,1)+XE(1,2)) ) * r_z
      BB(1,3) = ( XE(2,3)*(XE(1,2)+XE(1,3))-XE(2,1)*(XE(1,2)+XE(1,1))
     &           +X1s2*( XE(1,4)*(XE(2,3)-XE(2,1))
     &                  +XE(1,1)*(XE(2,3)+XE(2,4))
     &                  -XE(1,3)*(XE(2,4)+XE(2,1))) ) * r_z
      BB(2,4) = ( XE(1,1)*(XE(1,2)+XE(1,1))
     &           -XE(1,3)*(XE(1,3)+XE(1,2)) ) * r_z
      BB(1,5) = ( XE(2,4)*(XE(1,3)+XE(1,4))-XE(2,2)*(XE(1,3)+XE(1,2))
     &           +X1s2*( XE(1,1)*(XE(2,4)-XE(2,2))
     &                  +XE(1,2)*(XE(2,4)+XE(2,1))
     &                  -XE(1,4)*(XE(2,1)+XE(2,2))) ) * r_z
      BB(2,6) = ( XE(1,2)*(XE(1,2)+XE(1,3))
     &           -XE(1,4)*(XE(1,4)+XE(1,3)) ) * r_z
      BB(1,7) = ( XE(2,1)*(XE(1,4)+XE(1,1))-XE(2,3)*(XE(1,4)+XE(1,3))
     &           +X1s2*( XE(1,2)*(XE(2,1)-XE(2,3))
     &                  +XE(1,3)*(XE(2,1)+XE(2,2))
     &                  -XE(1,1)*(XE(2,2)+XE(2,3))) ) * r_z
      BB(2,8) = ( XE(1,3)*(XE(1,3)+XE(1,4))
     &           -XE(1,1)*(XE(1,1)+XE(1,4)) ) * r_z
*- Integration analytique des termes BB(3,1),BB(3,3),BB(3,5),BB(3,7)
      BB(3,1) = ( XE(1,1)*(XE(2,2)-XE(2,4))+XE(2,1)*(XE(1,4)-XE(1,2))
     &           +X1s2*( XE(1,2)*(XE(2,3)+XE(2,4))
     &                  +XE(1,3)*(XE(2,4)-XE(2,2))
     &                  -XE(1,4)*(XE(2,3)+XE(2,2))) ) * r_z
      BB(3,3) = ( XE(1,2)*(XE(2,3)-XE(2,1))+XE(2,2)*(XE(1,1)-XE(1,3))
     &           +X1s2*( XE(1,3)*(XE(2,1)+XE(2,4))
     &                  +XE(1,4)*(XE(2,1)-XE(2,3))
     &                  -XE(1,1)*(XE(2,3)+XE(2,4))) ) * r_z
      BB(3,5) = ( XE(1,3)*(XE(2,4)-XE(2,2))+XE(2,3)*(XE(1,2)-XE(1,4))
     &           +X1s2*( XE(1,4)*(XE(2,1)+XE(2,2))
     &                  +XE(1,1)*(XE(2,2)-XE(2,4))
     &                  -XE(1,2)*(XE(2,4)+XE(2,1))) ) * r_z
      BB(3,7) = ( XE(1,4)*(XE(2,1)-XE(2,3))+XE(2,4)*(XE(1,3)-XE(1,1))
     &           +X1s2*( XE(1,1)*(XE(2,2)+XE(2,3))
     &                  +XE(1,2)*(XE(2,3)-XE(2,1))
     &                  -XE(1,3)*(XE(2,1)+XE(2,2))) ) * r_z
*      write (*,*) 'calcul de aire pour bbdil'
*      write (*,*) 'aire = ',aire
*      write (*,*) 'ecriture de la matrice bb'
*      write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
      GOTO 666
 
*-----------------------------------------------------------------------
*----------------------------- Element ICT6 ----------------------------
*-----------------------------------------------------------------------
 30   CONTINUE
*      write(*,*) 'Element ICT6',ifour
*      write (*,*) 'ecriture de la matrice xe[ ]'
*      write (*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,6)
      GOTO (999,301,301,303),IFR
      GOTO 999
*
*--------contraintes planes ou deformations planes(ifour=-2,-1)
 301  CONTINUE
*--------donnee des composantes de bb-dil
      r_z =  (  XE(1,1)*(XE(2,2)-XE(2,6)) + XE(1,2)*(XE(2,3)-XE(2,1))
     &        + XE(1,3)*(XE(2,4)-XE(2,2)) + XE(1,4)*(XE(2,5)-XE(2,3))
     &        + XE(1,5)*(XE(2,6)-XE(2,4)) + XE(1,6)*(XE(2,1)-XE(2,5)) )
     &     + X1s4 * (  XE(1,1)*(XE(2,5)-XE(2,3))
     &               + XE(1,3)*(XE(2,1)-XE(2,5))
     &               + XE(1,5)*(XE(2,3)-XE(2,1)) )
      AIRE = X2s3 * r_z
      r_z = xUn / r_z
      BB(1, 1) = ( (XE(2,2)-XE(2,6)) + X1s4 * (XE(2,5)-XE(2,3)) ) * r_z
      BB(2, 2) = ( (XE(1,6)-XE(1,2)) + X1s4 * (XE(1,3)-XE(1,5)) ) * r_z
      BB(1, 3) = (XE(2,3)-XE(2,1)) * r_z
      BB(2, 4) = (XE(1,1)-XE(1,3)) * r_z
      BB(1, 5) = ( (XE(2,4)-XE(2,2)) + X1s4 * (XE(2,1)-XE(2,5)) ) * r_z
      BB(2, 6) = ( (XE(1,2)-XE(1,4)) + X1s4 * (XE(1,5)-XE(1,1)) ) * r_z
      BB(1, 7) = (XE(2,5)-XE(2,3)) * r_z
      BB(2, 8) = (XE(1,3)-XE(1,5)) * r_z
      BB(1, 9) = ( (XE(2,6)-XE(2,4)) + X1s4 * (XE(2,3)-XE(2,1)) ) * r_z
      BB(2,10) = ( (XE(1,4)-XE(1,6)) + X1s4 * (XE(1,1)-XE(1,3)) ) * r_z
      BB(1,11) = (XE(2,1)-XE(2,5)) * r_z
      BB(2,12) = (XE(1,5)-XE(1,1)) * r_z
 
*      write (*,*) 'calcul de aire pour bbdil'
*      write (*,*) 'aire = ',aire
*      write (*,*) 'ecriture de la matrice bb'
*      write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
      GOTO 666
*
*--------axisymetriques (ifour=0)
 303  CONTINUE
*--------donnee des composantes de bb-dil
*-----integration numerique des composantes bb-dil------
      AIRE = XZer
      DO i = 1,nbnn
        XNUM1(i) = XZer
        XNUM2(i) = XZer
        XNUM3(i) = XZer
      ENDDO
      DO IGAU = 1,NBPGAU
        DO i = 1,nbnn
          SHP(1,i) = SHPTOT(1,i,IGAU)
          SHP(2,i) = SHPTOT(2,i,IGAU)
          SHP(3,i) = SHPTOT(3,i,IGAU)
        ENDDO
        CALL DISTRR(XE,SHP,nbnn,RR)
        CALL DISTR2(XE,SHP,nbnn,RR2)
        CALL DISTR3(XE,SHP,nbnn,RR3)
        CALL DISTZ2(XE,SHP,nbnn,ZZ2)
        CALL DISTZ3(XE,SHP,nbnn,ZZ3)
        CALL JACOBI(XE,SHP,2,nbnn,DJAC)
        r_z  = POIGAU(IGAU) * RR
        r_zz2 = r_z * ZZ2
        r_zz3 = r_z * ZZ3
        r_rr2 = r_z * RR2
        r_rr3 = r_z * RR3
        r_z  = POIGAU(IGAU) * DJAC
        AIRE = AIRE + (r_z * RR)
        DO i = 1,nbnn
          XNUM1(i) = XNUM1(i) + (  SHPTOT(2,i,IGAU)*r_zz3
     &                           - SHPTOT(3,i,IGAU)*r_zz2 )
          XNUM2(i) = XNUM2(i) + (  SHPTOT(3,i,IGAU)*r_rr2
     &                           - SHPTOT(2,i,IGAU)*r_rr3 )
          XNUM3(i) = XNUM3(i) + (  SHPTOT(1,i,IGAU)*r_z )
        ENDDO
      ENDDO
      r_z = XUn / AIRE
      DO i = 1,nbnn
         k = 2*i
         BB(1,k-1) = XNUM1(i) * r_z
         BB(2,k)   = XNUM2(i) * r_z
         BB(3,k-1) = XNUM3(i) * r_z
      ENDDO
*      write (*,*) 'calcul de aire pour bbdil'
*      write (*,*) 'aire = ',aire
*      write (*,*) 'ecriture de la matrice bb'
*      write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn)
      GOTO 666
 
*-----------------------------------------------------------------------
*----------------------------- Element ICQ8 ----------------------------
*-----------------------------------------------------------------------
 40   CONTINUE
*      write(*,*) 'Element ICQ8',ifour
*      write(*,*) 'Ecriture de la matrice xe[ ]'
*      write(*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,6)
      GOTO (999,401,401,403),IFR
      GOTO 999
 
*--------contraintes planes ou deformations planes (ifour=-2,-1)
 401  CONTINUE
*--------donnee de la matrice symetrique p(3,3) (1 fois par element)
      P(1,1) =  X2s3 * ( XE(1,1)*(XE(2,2)-XE(2,8))
     1                 + XE(1,2)*(XE(2,3)-XE(2,1))
     1         +XE(1,3)*(XE(2,4)-XE(2,2))+XE(1,4)*(XE(2,5)-XE(2,3))
     2         +XE(1,5)*(XE(2,6)-XE(2,4))+XE(1,6)*(XE(2,7)-XE(2,5))
     3         +XE(1,7)*(XE(2,8)-XE(2,6))+XE(1,8)*(XE(2,1)-XE(2,7)) )
     4        + X1s6 * ( (XE(1,1)-XE(1,5))*(XE(2,7)-XE(2,3))
     5                  +(XE(1,3)-XE(1,7))*(XE(2,1)-XE(2,5)) )
      P(1,2) = 5.D0/9.D0*(-XE(2,2)*(XE(1,1)+XE(1,3))
     1         +XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,6)*(-XE(2,5)-XE(2,7))
     1         +XE(2,6)*(XE(1,5)+XE(1,7)))
     2         +4.D0/9.D0*(XE(1,8)*(XE(2,7)+XE(2,2)-XE(2,6)-XE(2,1))
     2         +XE(2,8)*(-XE(1,7)-XE(1,2)+XE(1,6)+XE(1,1))
     2         +XE(1,4)*(XE(2,2)-XE(2,6)+XE(2,5)-XE(2,3))
     2         +XE(2,4)*(-XE(1,2)+XE(1,6)-XE(1,5)+XE(1,3)))
     3         +7.D0/45.D0*(XE(1,2)*(XE(2,5)+XE(2,7))
     3         -XE(2,2)*(XE(1,5)+XE(1,7))-XE(1,6)*(XE(2,1)+XE(2,3))
     3         +XE(2,6)*(XE(1,1)+XE(1,3)))
     4         +8.D0/15.D0*(XE(1,6)*XE(2,2)-XE(1,2)*XE(2,6))
     5         +7.D0/90.D0*(XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
     5         +XE(1,5)*XE(2,3)-XE(1,1)*XE(2,7))
     6         +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7)
     6         -XE(1,5)*XE(2,1)+XE(1,1)*XE(2,5))
      P(1,3) = 5.D0/9.D0*(-XE(1,8)*(XE(2,7)+XE(2,1))
     1         +XE(2,8)*(XE(1,7)+XE(1,1))
     1         +XE(1,4)*(XE(2,3)+XE(2,5))
     1         -XE(2,4)*(XE(1,3)+XE(1,5)))
     2        +4.D0/9.D0*((XE(1,8)-XE(1,4))*(XE(2,2)+XE(2,6))
     2         +(-XE(2,8)+XE(2,4))*(XE(1,2)+XE(1,6))
     2         +XE(1,2)*(XE(2,1)-XE(2,3))+XE(2,2)*(-XE(1,1)+XE(1,3))
     2         +XE(1,6)*(XE(2,7)-XE(2,5))+XE(2,6)*(-XE(1,7)+XE(1,5)))
     3         +7.D0/45.D0*(-XE(1,8)*(XE(2,3)+XE(2,5))
     3         +XE(2,8)*(XE(1,3)+XE(1,5))
     3         +XE(1,4)*(XE(2,1)+XE(2,7))
     3         -XE(2,4)*(XE(1,1)+XE(1,7)))
     4         +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7)
     4         +XE(1,5)*XE(2,1)-XE(1,1)*XE(2,5))
     5         +7.D0/90.D0*(XE(1,7)*XE(2,5)-XE(1,3)*XE(2,1)
     5         -XE(1,5)*XE(2,7)+XE(1,1)*XE(2,3))
     6         +8.D0/15.D0*(XE(1,8)*XE(2,4)-XE(1,4)*XE(2,8))
      P(1,4) = XZer
 
      P(2,1) = P(1,2)
      P(2,2) = 16.D0/45.D0*(XE(2,6)*(XE(1,5)-XE(1,7))
     1         +XE(1,6)*(XE(2,7)-XE(2,5))+XE(1,2)*(XE(2,3)-XE(2,1))
     1         +XE(2,2)*(XE(1,1)-XE(1,3)))
     2         +14.D0/45.D0*(XE(2,4)*(XE(1,3)-XE(1,5))
     2         +XE(1,4)*(XE(2,5)-XE(2,3))+XE(1,8)*(XE(2,1)-XE(2,7))
     2         +XE(2,8)*(XE(1,7)-XE(1,1)))
     3         +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6))
     3         +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4))
     3         +XE(1,6)*(XE(2,4)-XE(2,8)))
     4         +2.D0/45.D0*(XE(2,2)*(XE(1,5)-XE(1,7))
     4         +XE(1,2)*(XE(2,7)-XE(2,5)) +XE(1,6)*(XE(2,3)-XE(2,1))
     4         +XE(2,6)*(XE(1,1)-XE(1,3)))
     5         +4.D0/45.D0*(XE(2,4)*(XE(1,1)-XE(1,7))
     5         +XE(1,4)*(XE(2,7)-XE(2,1)) +XE(1,8)*(XE(2,3)-XE(2,5))
     5         +XE(2,8)*(XE(1,5)-XE(1,3)))
     6         +17.D0/90.D0*(XE(1,7)*XE(2,5)+XE(1,3)*XE(2,1)
     6         -XE(1,5)*XE(2,7)-XE(1,1)*XE(2,3))
     7         +1.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
     7         +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7))
      P(2,3) = 1.D0/3.D0*(XE(1,5)*(XE(2,6)-XE(2,4))
     1         +XE(1,8)*(XE(2,7)+XE(2,1))+XE(1,7)*(XE(2,6)-XE(2,8))
     1         -XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,1)*(XE(2,2)-XE(2,8))
     1         +XE(1,4)*(XE(2,5)+XE(2,3))-XE(1,6)*(XE(2,7)+XE(2,5))
     1         +XE(1,3)*(XE(2,2)-XE(2,4)))
     2         +4.D0/9.D0*(-XE(1,4)*(XE(2,2)+XE(2,6))
     2         -XE(1,8)*(XE(2,6)+XE(2,2))+XE(1,2)*(XE(2,8)+XE(2,4))
     2         +XE(1,6)*(XE(2,4)+XE(2,8)))
     3         +1.D0/9.D0*(XE(1,7)*(XE(2,2)-XE(2,4))
     3         +XE(1,8)*(XE(2,3)+XE(2,5))
     3         +XE(1,5)*(XE(2,2)-XE(2,8))+XE(1,3)*(XE(2,6)-XE(2,8))
     3         -XE(1,2)*(XE(2,5)+XE(2,7))+XE(1,1)*(XE(2,6)-XE(2,4))
     3         +XE(1,4)*(XE(2,1)+XE(2,7))-XE(1,6)*(XE(2,1)+XE(2,3)))
      P(2,4) = XZer
 
      P(3,1) = P(1,3)
      P(3,2) = P(2,3)
      P(3,3) =  14.D0/45.D0*(XE(1,6)*(XE(2,7)-XE(2,5))
     1         +XE(2,6)*(XE(1,5)-XE(1,7))+XE(1,2)*(XE(2,3)-XE(2,1))
     1         +XE(2,2)*(XE(1,1)-XE(1,3)))
     2         +16.D0/45.D0*(XE(1,8)*(XE(2,1)-XE(2,7))
     2         +XE(2,8)*(XE(1,7)-XE(1,1))+XE(1,4)*(XE(2,5)-XE(2,3))
     2         +XE(2,4)*(XE(1,3)-XE(1,5)))
     3         +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6))
     3         +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4))
     3         +XE(1,6)*(XE(2,4)-XE(2,8)))
     4         +2.D0/45.D0*(XE(2,4)*(XE(1,7)-XE(1,1))
     4         +XE(1,8)*(XE(2,5)-XE(2,3))+XE(2,8)*(-XE(1,5)+XE(1,3))
     4         +XE(1,4)*(XE(2,1)-XE(2,7)))
     5         +4.D0/45.D0*(XE(2,2)*(XE(1,7)-XE(1,5))
     5         +XE(1,2)*(XE(2,5)-XE(2,7))+XE(1,6)*(XE(2,1)-XE(2,3))
     5         +XE(2,6)*(XE(1,3)-XE(1,1)))
     6         +1.D0/90.D0*(-XE(1,5)*XE(2,7)+XE(1,3)*XE(2,1)
     6         +XE(1,7)*XE(2,5)-XE(1,1)*XE(2,3))
     7         +17.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5)
     7         +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7))
      P(3,4) = XZer
 
      P(4,1) = XZer
      P(4,2) = XZer
      P(4,3) = XZer
      P(4,4) = XZer
 
*      write (*,*) 'ecriture de la matrice p[ ]'
*      write (*,4) (('p(',i,',',j,')=',p(i,j),i=1,3),j=1,3)
 
*-----------------donnee des vecteurs {t}={t1,t2,t3}--------------------
      A(1,1) = 2.D0/3.D0*(XE(2,2)-XE(2,8))
     1        +1.D0/6.D0*(XE(2,7)-XE(2,3))
      A(2,1) = 1.D0/9.D0*(4.D0*XE(2,8)-5.D0*XE(2,2))
     1        +7.D0/90.D0*(-XE(2,7)+2.D0*XE(2,6))
     1        +1.D0/30.D0*XE(2,5)
      A(3,1) = 1.D0/9.D0*(-4.D0*XE(2,2)+5.D0*XE(2,8))
     1        +7.D0/90.D0*(XE(2,3)-2.D0*XE(2,4))
     1        -1.D0/30.D0*XE(2,5)
 
      A(1,2) = 2.D0/3.D0*(XE(1,8)-XE(1,2))
     1        +1.D0/6.D0*(XE(1,3)-XE(1,7))
      A(2,2) = 1.D0/9.D0*(-4.D0*XE(1,8)+5.D0*XE(1,2))
     1        +7.D0/90.D0*(XE(1,7)-2.D0*XE(1,6))
     1        -1.D0/30.D0*XE(1,5)
      A(3,2) = 1.D0/9.D0*(4.D0*XE(1,2)-5.D0*XE(1,8))
     1        +7.D0/90.D0*(-XE(1,3)+2.D0*XE(1,4))
     1        +1.D0/30.D0*XE(1,5)
 
      A(1,3) = 2.D0/3.D0*(-XE(2,1)+XE(2,3))
      A(2,3) = 5.D0/9.D0*(XE(2,1)+XE(2,3))
     1        +7.D0/45.D0*(XE(2,5)+XE(2,7))
     1        -4.D0/9.D0*(XE(2,4)+XE(2,8))
     1        -8.D0/15.D0*XE(2,6)
      A(3,3) = 4.D0/9.D0*(XE(2,1)-XE(2,3)+XE(2,4)-XE(2,8))
 
      A(1,4) = 2.D0/3.D0*(-XE(1,3)+XE(1,1))
      A(2,4) = -5.D0/9.D0*(XE(1,1)+XE(1,3))
     1        -7.D0/45.D0*(XE(1,5)+XE(1,7))
     1        +4.D0/9.D0*(XE(1,4)+XE(1,8))
     1        +8.D0/15.D0*XE(1,6)
      A(3,4) = 4.D0/9.D0*(-XE(1,1)+XE(1,3)-XE(1,4)+XE(1,8))
 
      A(1,5) = 2.D0/3.D0*(XE(2,4)-XE(2,2))
     1         +1.D0/6.D0*(-XE(2,5)+XE(2,1))
      A(2,5) = 1.D0/9.D0*(4.D0*XE(2,4)-5.D0*XE(2,2))
     1        +7.D0/90.D0*(-XE(2,5)+2.D0*XE(2,6))
     1        +1.D0/30.D0*XE(2,7)
      A(3,5) = 1.D0/9.D0*(4.D0*XE(2,2)-5.D0*XE(2,4))
     1        +7.D0/90.D0*(-XE(2,1)+2.D0*XE(2,8))
     1        +1.D0/30.D0*XE(2,7)
 
      A(1,6) = 2.D0/3.D0*(-XE(1,4)+XE(1,2))
     1         +1.D0/6.D0*(XE(1,5)-XE(1,1))
      A(2,6) = 1.D0/9.D0*(-4.D0*XE(1,4)+5.D0*XE(1,2))
     1        +7.D0/90.D0*(XE(1,5)-2.D0*XE(1,6))
     1        -1.D0/30.D0*XE(1,7)
      A(3,6) = 1.D0/9.D0*(-4.D0*XE(1,2)+5.D0*XE(1,4))
     1        +7.D0/90.D0*(XE(1,1)-2.D0*XE(1,8))
     1        -1.D0/30.D0*XE(1,7)
 
      A(1,7) = 2.D0/3.D0*(-XE(2,3)+XE(2,5))
      A(2,7) = 4.D0/9.D0*(XE(2,5)-XE(2,3)+XE(2,2)-XE(2,6))
      A(3,7) = 5.D0/9.D0*(XE(2,3)+XE(2,5))
     1        +7.D0/45.D0*(XE(2,1)+XE(2,7))
     1        -4.D0/9.D0*(XE(2,2)+XE(2,6))
     1        -8.D0/15.D0*XE(2,8)
 
      A(1,8) = 2.D0/3.D0*(XE(1,3)-XE(1,5))
      A(2,8) = 4.D0/9.D0*(-XE(1,5)+XE(1,3)-XE(1,2)+XE(1,6))
      A(3,8) = -5.D0/9.D0*(XE(1,3)+XE(1,5))
     1        -7.D0/45.D0*(XE(1,1)+XE(1,7))
     1        +4.D0/9.D0*(XE(1,2)+XE(1,6))
     1        +8.D0/15.D0*XE(1,8)
 
      A(1,9) = 2.D0/3.D0*(-XE(2,4)+XE(2,6))
     1         +1.D0/6.D0*(XE(2,3)-XE(2,7))
      A(2,9) = 1.D0/9.D0*(-4.D0*XE(2,4)+5.D0*XE(2,6))
     1        +7.D0/90.D0*(XE(2,3)-2.D0*XE(2,2))
     1        -1.D0/30.D0*XE(2,1)
      A(3,9) = 1.D0/9.D0*(4.D0*XE(2,6)-5.D0*XE(2,4))
     1        +7.D0/90.D0*(-XE(2,7)+2.D0*XE(2,8))
     1        +1.D0/30.D0*XE(2,1)
 
      A(1,10) = 2.D0/3.D0*(XE(1,4)-XE(1,6))
     1         +1.D0/6.D0*(-XE(1,3)+XE(1,7))
      A(2,10) = 1.D0/9.D0*(4.D0*XE(1,4)-5.D0*XE(1,6))
     1        +7.D0/90.D0*(-XE(1,3)+2.D0*XE(1,2))
     1        +1.D0/30.D0*XE(1,1)
      A(3,10) = 1.D0/9.D0*(-4.D0*XE(1,6)+5.D0*XE(1,4))
     1        +7.D0/90.D0*(XE(1,7)-2.D0*XE(1,8))
     1        -1.D0/30.D0*XE(1,1)
 
      A(1,11) = 2.D0/3.D0*(-XE(2,5)+XE(2,7))
      A(2,11) = -5.D0/9.D0*(XE(2,7)+XE(2,5))
     1        -7.D0/45.D0*(XE(2,1)+XE(2,3))
     1        +4.D0/9.D0*(XE(2,4)+XE(2,8))
     1        +8.D0/15.D0*XE(2,2)
      A(3,11) = 4.D0/9.D0*(XE(2,4)-XE(2,5)-XE(2,8)+XE(2,7))
 
      A(1,12) = 2.D0/3.D0*(XE(1,5)-XE(1,7))
      A(2,12) = 5.D0/9.D0*(XE(1,7)+XE(1,5))
     1        +7.D0/45.D0*(XE(1,1)+XE(1,3))
     1        -4.D0/9.D0*(XE(1,4)+XE(1,8))
     1        -8.D0/15.D0*XE(1,2)
      A(3,12) = 4.D0/9.D0*(-XE(1,4)+XE(1,5)+XE(1,8)-XE(1,7))
 
      A(1,13) = 2.D0/3.D0*(-XE(2,6)+XE(2,8))
     1         +1.D0/6.D0*(XE(2,5)-XE(2,1))
      A(2,13) = 1.D0/9.D0*(-4.D0*XE(2,8)+5.D0*XE(2,6))
     1        +7.D0/90.D0*(XE(2,1)-2.D0*XE(2,2))
     1        -1.D0/30.D0*XE(2,3)
      A(3,13) = 1.D0/9.D0*(-4.D0*XE(2,6)+5.D0*XE(2,8))
     1        +7.D0/90.D0*(XE(2,5)-2.D0*XE(2,4))
     1        -1.D0/30.D0*XE(2,3)
 
      A(1,14) = 2.D0/3.D0*(XE(1,6)-XE(1,8))
     1         +1.D0/6.D0*(-XE(1,5)+XE(1,1))
      A(2,14) = 1.D0/9.D0*(4.D0*XE(1,8)-5.D0*XE(1,6))
     1        +7.D0/90.D0*(-XE(1,1)+2.D0*XE(1,2))
     1        +1.D0/30.D0*XE(1,3)
      A(3,14) = 1.D0/9.D0*(4.D0*XE(1,6)-5.D0*XE(1,8))
     1        +7.D0/90.D0*(-XE(1,5)+2.D0*XE(1,4))
     1        +1.D0/30.D0*XE(1,3)
 
      A(1,15) = 2.D0/3.D0*(-XE(2,7)+XE(2,1))
      A(2,15) = 4.D0/9.D0*(-XE(2,6)+XE(2,7)+XE(2,2)-XE(2,1))
      A(3,15) = -5.D0/9.D0*(XE(2,1)+XE(2,7))
     1        -7.D0/45.D0*(XE(2,3)+XE(2,5))
     1        +4.D0/9.D0*(XE(2,2)+XE(2,6))
     1        +8.D0/15.D0*XE(2,4)
 
      A(1,16) = 2.D0/3.D0*(XE(1,7)-XE(1,1))
      A(2,16) = 4.D0/9.D0*(XE(1,6)-XE(1,7)-XE(1,2)+XE(1,1))
      A(3,16) = 5.D0/9.D0*(XE(1,1)+XE(1,7))
     1        +7.D0/45.D0*(XE(1,3)+XE(1,5))
     1        -4.D0/9.D0*(XE(1,2)+XE(1,6))
     1        -8.D0/15.D0*XE(1,4)
 
*      write (*,*) 'ecriture de la matrice t[ ]'
*      write (*,4) (('t(',i,',',j,')=',a(i,j),i=1,3),j=1,16)
*--------------resolution matricielle de [p]{a}={t}---------------------
*      write (*,*) 'Appel a GAUSSJ'
      N = 3
      NP = 4
      M = 16
      MP = 60
      CALL GAUSSJ(P,N,NP,A,M,MP)
      GOTO 666
 
*--------axisymetriques (ifour=0)
 403  CONTINUE
*--------donnee de la matrice symetrique p(3,3) (1 fois par element)
      CALL ZERO(P,4,4)
      DO IGAU = 1,NBPGAU
        DO I = 1,nbnn
          SHP(1,I) = SHPTOT(1,I,IGAU)
          SHP(2,I) = SHPTOT(2,I,IGAU)
          SHP(3,I) = SHPTOT(3,I,IGAU)
        ENDDO
        CALL DISTRR(XE,SHP,nbnn,RR)
        CALL JACOBI(XE,SHP,2,nbnn,DJAC)
        r_z = POIGAU(IGAU)*RR*DJAC
        P(1,1) = P(1,1) + r_z
        P(2,1) = P(2,1) + (r_z * QSIGAU(IGAU))
        P(3,1) = P(3,1) + (r_z * ETAGAU(IGAU))
        P(2,2) = P(2,2) + (r_z * QSIGAU(IGAU) * QSIGAU(IGAU))
        P(3,2) = P(3,2) + (r_z * QSIGAU(IGAU) * ETAGAU(IGAU))
        P(3,3) = P(3,3) + (r_z * ETAGAU(IGAU) * ETAGAU(IGAU))
      ENDDO
      P(1,2) = P(2,1)
      P(1,3) = P(3,1)
      P(2,3) = P(3,2)
*      write (*,*) 'ecriture de la matrice p[ ]'
*      write (*,4) (('p(',i,',',j,')=',p(i,j),i=1,3),j=1,3)
*-----------calcul des vecteurs {t}={t1,t2,t3}--integration numerique---
      CALL ZERO(A,4,24)
      DO 4031 IGAU = 1,NBPGAU
        DO I = 1,nbnn
          SHP(1,I) = SHPTOT(1,I,IGAU)
          SHP(2,I) = SHPTOT(2,I,IGAU)
          SHP(3,I) = SHPTOT(3,I,IGAU)
        ENDDO
        CALL DISTRR(XE,SHP,nbnn,RR)
        CALL DISTR2(XE,SHP,nbnn,RR2)
        CALL DISTR3(XE,SHP,nbnn,RR3)
        CALL DISTZ2(XE,SHP,nbnn,ZZ2)
        CALL DISTZ3(XE,SHP,nbnn,ZZ3)
        CALL JACOBI(XE,SHP,2,nbnn,DJAC)
        K = 0
        DO I = 1,nbnn
          K = K+1
          r_z =  POIGAU(IGAU)*RR
     &          *(SHPTOT(2,I,IGAU)*ZZ3-SHPTOT(3,I,IGAU)*ZZ2)
          A(1,K) = A(1,K) + r_z
          A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
          A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
          K = K+1
          r_z = POIGAU(IGAU)*RR
     &          *(-SHPTOT(2,I,IGAU)*RR3+SHPTOT(3,I,IGAU)*RR2)
          A(1,K) = A(1,K) + r_z
          A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
          A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
          K = K+1
          r_z = POIGAU(IGAU)*SHPTOT(1,I,IGAU)*DJAC
          A(1,K) = A(1,K) + r_z
          A(2,K) = A(2,K) + r_z * QSIGAU(IGAU)
          A(3,K) = A(3,K) + r_z * ETAGAU(IGAU)
        ENDDO
 4031 CONTINUE
*      write(*,*) 'Ecriture de la matrice t[ ]'
*      write(*,4) (('t(',i,',',j,')=',a(i,j),i=1,3),j=1,24)
*------------resolution matricielle de [p]{a}={t}-----------------
      N = 3
      NP = 4
      M = 24
      MP = 60
*      write (*,*) 'appel a GAUSSJ'
      CALL GAUSSJ(P,N,NP,A,M,MP)
*      write(*,*) 'Ecriture de la matrice a[ ]'
*      write(*,4) (('a(',i,',',j,')=',a(i,j),i=1,3),j=1,24)
      GOTO 666
 
C=======================================================================
C========== Elements MASSIFS INCOMPRESSIBLES TRIDIMENSIONNELS ==========
C=======================================================================
 100  CONTINUE
      IF (ifour.NE.2) GOTO 999
 
      DIM3   = 1.D0
      FACAR  = 1.D0
      FACSCA = 1.D0
      ESTEL  = XZero
 
C= Boucle sur les points d'integration
      DO iGau=1,NBPGAU
 
* BGENE au point IGAU
        CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
     &              MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,DIM3,
     &              XE,SHPTOT,SHP,BGENE,DJAC,XDPGE,YDPGE)
 
C= Calcul de la composante integree en ce point de Gauss
        DJAC=ABS(DJAC)*POIGAU(iGau)*FACAR*DIM3
        ESTEL=ESTEL+FACSCA*DJAC
 
*----- Element ICQ8 ICCU20(76) ICTE10(77) ICPR15(78) ICPY13(274)
        IF (MELE.EQ.72 .or. (MELE.ge.76.and.MELE.le.78) .or.
     &      MELE.EQ.274) then
* second membre T
          K  = 1
          jj = 0
          DO i = 1,NBNN
            jj = jj + 1
            A(1,jj) = A(1,jj) + BGENE(1,K)*DJAC
            A(2,jj) = A(2,jj) + BGENE(1,K)*qsigau(igau)*DJAC
            A(3,jj) = A(3,jj) + BGENE(1,K)*etagau(igau)*DJAC
            if (mele.ge.76)
     &        A(4,jj) = A(4,jj) + BGENE(1,K)*dzegau(igau)*DJAC
            jj = jj + 1
            A(1,jj) = A(1,jj) + BGENE(2,K+1)*DJAC
            A(2,jj) = A(2,jj) + BGENE(2,K+1)*qsigau(igau)*DJAC
            A(3,jj) = A(3,jj) + BGENE(2,K+1)*etagau(igau)*DJAC
            if (mele.ge.76) then
              A(4,jj) = A(4,jj) + BGENE(2,K+1)*dzegau(igau)*DJAC
              jj = jj + 1
              A(1,jj) = A(1,jj) + BGENE(3,K+2)*DJAC
              A(2,jj) = A(2,jj) + BGENE(3,K+2)*qsigau(igau)*DJAC
              A(3,jj) = A(3,jj) + BGENE(3,K+2)*etagau(igau)*DJAC
              A(4,jj) = A(4,jj) + BGENE(3,K+2)*dzegau(igau)*DJAC
            endif
            K = K + KINC
          ENDDO
* matrice M
          P(1,1) = P(1,1) + DJAC
          P(1,2) = P(1,2) + qsigau(igau)*DJAC
          P(1,3) = P(1,3) + etagau(igau)*DJAC
          P(2,2) = P(2,2) + qsigau(igau)*qsigau(igau)*DJAC
          P(2,3) = P(2,3) + qsigau(igau)*etagau(igau)*DJAC
          P(3,3) = P(3,3) + etagau(igau)*etagau(igau)*DJAC
          if (mele.ge.76) then
            P(1,4) = P(1,4) + dzegau(igau)*DJAC
            P(2,4) = P(2,4) + qsigau(igau)*dzegau(igau)*DJAC
            P(3,4) = P(3,4) + etagau(igau)*dzegau(igau)*DJAC
            P(4,4) = P(4,4) + dzegau(igau)*dzegau(igau)*DJAC
          endif
 
        else if(mele.le.71) then
 
          K = 1
          DO i = 1,NBNN
            bb(1,K) = bb(1,K) + BGENE(1,K)*DJAC
            bb(2,K+1) = bb(2,K+1) + BGENE(2,K+1)*DJAC
            bb(3,K) = bb(3,K) + BGENE(3,K)*DJAC
            K = K + KINC
          ENDDO
 
        else if ( (mele.ge.73.and.mele.le.75)
     &.or.mele.eq.273) then
 
          K = 1
          DO i = 1,NBNN
            bb(1,K) = bb(1,K) + BGENE(1,K)*DJAC
            bb(2,K+1) = bb(2,K+1) + BGENE(2,K+1)*DJAC
            bb(3,K+2) = bb(3,K+2) + BGENE(3,K+2)*DJAC
            K = K + KINC
          ENDDO
 
        endif
 
      ENDDO
 
      if (MELE.EQ.72) then
        P(1,4) = XZer
        P(2,1) = P(1,2)
        P(2,4) = XZer
        P(3,1) = P(1,3)
        P(3,2) = P(2,3)
        P(3,4) = XZer
 
        P(4,1) = XZer
        P(4,2) = XZer
        P(4,3) = XZer
        P(4,4) = XZer
 
      ELSE IF (mele.GE.76) THEN
        P(2,1) = P(1,2)
        P(3,1) = P(1,3)
        P(3,2) = P(2,3)
        P(4,1) = P(1,4)
        P(4,2) = P(2,4)
        P(4,3) = P(3,4)
      ENDIF
*
** synthese
*
      if (MELE.EQ.72) then
        N = 3
        NP = 4
        M = 16
        MP = 60
*        write (*,*) 'appel a GAUSSJ'
        CALL GAUSSJ(P,N,NP,A,M,MP)
*        write(*,*) 'Ecriture de la matrice a[ ]'
*        write(*,4) (('a(',i,',',j,')=',a(i,j),i=1,3),j=1,24)
 
              elseif (mele.ge.76.and.mele.ne.273) then
      N = 4
      NP = 4
      M = NBNN*3
      MP = 60
*      write (*,*) 'appel a GAUSSJ'
      CALL GAUSSJ(P,N,NP,A,M,MP)
 
              elseif(mele.le.71) then
 
           K = 1
           DO i = 1,NBNN
             bb(1,K) = bb(1,K)/ESTEL
             bb(2,K+1) = bb(2,K+1)/ESTEL
             bb(3,K) = bb(3,K)/ESTEL
             K = K + KINC
           ENDDO
 
              elseif((mele.ge.73.and.mele.le.75).or.
     &mele.eq.273) then
 
           K = 1
           DO i = 1,NBNN
             bb(1,K) = bb(1,K)/ESTEL
             bb(2,K+1) = bb(2,K+1)/ESTEL
             bb(3,K+2) = bb(3,K+2)/ESTEL
             K = K + KINC
           ENDDO
 
              endif
 
      GOTO 666
 
C=======================================================================
C======================= ERREUR : Cas non prevus =======================
C=======================================================================
 999  CONTINUE
      write(6,*) 'Erreur d aiguillage - cas non prevu',MELE,IFOUR
      CALL ERREUR(5)
 
 666  CONTINUE
*      write(6,*) 'Sortie de BBCALC'
 
c      RETURN
      END