      subroutine fonmit(npo,npt,x,y,z,xnorm,t,vdtq2,vp2,f,fz,ind,gzz,tt)
c ............................................................
c
c   calcul de F = g_*,* aux tying points F= F+fz*z
c                       ou aux points d'integration

c   npo        : nombre de points
c   npt        : somme sur les points d'interpolation 
c   x,y,z      : coordonnes des noeuds de l'element courant
c   xnorm      : normale aux noeuds
c   t          : epaisseur de la coque (donee aux noeuds)
c   vdtq2(1,*) : derivee par rapport a r aux pt d'interpolation (tying)
c   vdtq2(2,*) : derivee par rapport a s ux pt d'interpolation
c   vp2        : valeur des polynomes de base aux pt d'interpolation
c
C  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     PROGRAMMEUR  :  Marina Vidrascu INRIA 2000
C  .................................................................. 
      double precision xnorm(3,npo),vdtq2(2,npo,npt),t(npo),tt(npt),
     *                 gzz(npt),vp2(npo,npt),F(3,3,npt),fz(3,3,npt),zero
      double precision x(npo),y(npo),z(npo)
      parameter (zero = 0.)



      CALL Dcopy(9*npt,zero,0,f,1)
      CALL Dcopy(9*npt,zero,0,fz,1)

      if (ind .eq. 1)then
          CALL Dcopy(npt,zero,0,gzz,1)
          CALL Dcopy(npt,zero,0,tt,1)
      endif
      do 1 l = 1,npt
      do 1 i = 1,npo

         f(1,1,l)  = f(1,1,l)  + vdtq2(1,i,l)* x(i)
         fz(1,1,l) = fz(1,1,l) + .5*t(i)*vdtq2(1,i,l)*xnorm(1,i)
   
         f(1,2,l)  = f(1,2,l)  + vdtq2(1,i,l)* y(i)
         fz(1,2,l) = fz(1,2,l) + .5*t(i)*vdtq2(1,i,l)*xnorm(2,i)

         f(1,3,l)  = f(1,3,l)  + vdtq2(1,i,l)* z(i)
         fz(1,3,l) = fz(1,3,l) +  .5*t(i)*vdtq2(1,i,l)*xnorm(3,i)

         f(2,1,l)  = f(2,1,l)  +  vdtq2(2,i,l)* x(i)
         fz(2,1,l) = fz(2,1,l) +  .5*t(i)*vdtq2(2,i,l)*xnorm(1,i)

         f(2,2,l)  = f(2,2,l)  +  vdtq2(2,i,l)* y(i)
         fz(2,2,l) = fz(2,2,l) + .5*t(i)*vdtq2(2,i,l)*xnorm(2,i)

         f(2,3,l)  = f(2,3,l)  +  vdtq2(2,i,l)* z(i)
         fz(2,3,l) = fz(2,3,l) +  .5*t(i)*vdtq2(2,i,l)*xnorm(3,i)

         f(3,1,l)  = f(3,1,l)  +  .5*t(i)*vp2(i,l)*xnorm(1,i)
         f(3,2,l)  = f(3,2,l)  +  .5*t(i)*vp2(i,l)*xnorm(2,i)
         f(3,3,l)  = f(3,3,l)  +  .5*t(i)*vp2(i,l)*xnorm(3,i)

         if (ind .eq. 1) then

             gzz(l)   = gzz(l) + .25*(t(i)**2)*vp2(i,l)
             tt(l)    = tt(l) + t(i)*vp2(i,l)

         endif

    1 continue

      end
      
