      subroutine mitlin(nno,npt1,npt2,npt3,npt,irint,npi,poids,npiz,
     *                  zint,poiz,xint,yint,vdtq2,vt2,f,fz,inmitc,
     *                  e,enu,t,x,y,z, vdpq2,vp2,b,bz,bzz,xnorm,v1,v2,
     *                  fi,fiz,gzz,tt,ijt,ae,nmat)

c ............................................................
C  BUT :   CALCUL DE LA MATRICE ELEMENTAIRE DE RAIDEUR  pour 
C  ---     une coque mitc
C
c   calcul de la matrice d'elasticite elas au points d'integration
c   (polynome de degre 2 en z)
c   calcul de tB*elas*B ou B=
c ............................................................
c
c   nno           : nombre de noeuds de l'element
c   npi           : nombre de points d'integration (en r,s)
c   poids         : poids de la formule d'integration en r,s
c   npiz          : nombre de points d'integration en z
c   zint(npiz)    : coord des nzpi points d'integration en z sur [-1,1]
c   poiz          : poids de la formule d'integration en z
c   e,enu         : coeff de Young, Poisson  aux points d'integration
c   t             : epaisseur de la coque (donee aux noeuds)
c   x,y,z         : coordonnes des noeuds de l'element courant
c   vdpq2(1,*)    : derivee par rapport a r aux pt integration
c   vdpq2(2,*)    : derivee par rapport a s aux pt integration
c   vp2           : valeur des polynomes de base aux pt d'integration
c   b,bz,bzz      : B = b+ bz*z+ bzz*z**2
c   xnorm         : normale aux noeuds
c   ijt           : permutation composante / noeud
c   ae            : matrice elementaire triangulaire sup
C  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     PROGRAMMEUR  :  Marina Vidrascu INRIA 2002
C  ..................................................................
      double precision ae(nmat),vdpq2(2,nno,npi),vdtq2(2,nno,npt),     
     *                 xnorm(3,nno),ddot,
     *                 vp2(nno,npi),vt2(nno,npt),f(3,3,npt),
     *                 fz(3,3,npt),gzz(npi),fi(3,3,npi),fiz(3,3,npi),
     *                 zero ,grr, grrz, grrzz,
     *                 grs, grsz, grszz, gss,gssz,gsszz,b(5,5*nno),
     *                 bz(5,5*nno),bzz(5,5*nno),elas(9),
     *                 e,enu,t(nno),tt(npi),deltat,
     *                 s,poids(npi),poiz(npiz) ,un, 
     *                 girr,giss,girs,v1(3,nno),v2(3,nno)

      character*4 na
      integer ijt(nno*5)
      double precision x(nno),y(nno),z(nno),
     *          xint(npi),yint(npi),zint(npiz)
      external inmitc
      parameter ( zero = 0. , un=1.)

c
c     calcul de la matrice de rigidite
c     --------------------------------
c     Somme(omega) Somme(-1,1) {BT * C * B}
c
c
c     boucle sur les points d'integration (r,s)
c
      do 2 l = 1,npi
      
c
c     calcul de B au points d'integration (utiliser B au tying points)
c
          call bremit(nno,npt1,npt2,npt3,npt,xnorm,v1,v2,t,npi,xint,
     *            yint,vdtq2,vt2,f,fz,inmitc,b,bz,bzz,
     *            irint,fi(1,1,l),fiz(1,1,l),vdpq2,vp2,l)


c
c        calcul de grr gss grs comme polynomes de degre 2 en z
c        gr (lig 1 de fi) gs (lig 2 de fi)

         grr   = ddot(3,fi(1,1,l),3,fi(1,1,l),3) 
         grrz  = 2*ddot(3,fi(1,1,l),3,fiz(1,1,l),3) 
         grrzz = ddot(3,fiz(1,1,l),3,fiz(1,1,l),3) 

         gss   = ddot(3,fi(2,1,l),3,fi(2,1,l),3) 
         gssz  = 2*ddot(3,fi(2,1,l),3,fiz(2,1,l),3) 
         gsszz = ddot(3,fiz(2,1,l),3,fiz(2,1,l),3) 

         grs   = ddot(3,fi(1,1,l),3,fi(2,1,l),3) 
         grsz  = ddot(3,fi(1,1,l),3,fiz(2,1,l),3) + 
     *           ddot(3,fi(2,1,l),3,fiz(1,1,l),3)
         grszz = ddot(3,fiz(1,1,l),3,fiz(2,1,l),3) 


c
c        boucle sur les points d'integration z
c
         do 10 lz = 1,npiz
c
c           calcul de la matrice d'elasticite
c
            call melmit(grr,grrz,grrzz,grs,grsz,grszz,gss,gssz,gsszz,
     *                  gzz,zint,npiz,npi,enu,e,tt,l,lz,deltat,elas,
     *                  girr,giss,girs)
c
c           calcul de TB*elas*B
c
            do 3 j = 1,5*nno
            do 3 i = 1,j
c
c              permutation ds la matrice
c
               II = MIN(IJT(I),IJT(J))
               JJ = MAX(IJT(I),IJT(J))
               K1 = JJ*(JJ-1)/2+II
               s  = 0.
c              2 blocs 1 3*3 et un 3*2
               do 4 m = 1,3
               do 4 k = 1,3
                  if (k .le. m) then
                     ie = (m-1)*m/2+k
                  else
                     ie = (k-1)*k/2+m
                  endif
c                 s = s+ somme (k=1,3) B(m,i)*elas(m,k)*B(k,j)
                  s = s + elas(ie)*
     *             (b(m,i)+zint(lz)*bz(m,i)+zint(lz)**2*bzz(m,i))*
     *             (b(k,j)+zint(lz)*bz(k,j)+zint(lz)**2*bzz(k,j))
 4             continue
               do 5 m = 4,5
               do 5 k = 4,5
                  if (k .le. m) then
                     ie = 6+(m-4)*(m-3)/2+k-3
                  else
                     ie = 6+(k-4)*(k-3)/2+m-3
                  endif
c                 s = s+ somme (k=1,3) B(m,i)*elas(m,k)*B(k,j)
                  s = s + elas(ie)*b(m,i)*b(k,j)
 5             continue
               ae(k1) = ae(k1)+s*deltat*poids(l)*poiz(lz) 
 3          continue
 10    continue

 2     continue
       end


