vectos.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE vectos
!                   *****************
!
     &(svec,vec,op,formul,
     & xmul,f,g,h,u,v,w,sf,sg,sh,su,sv,sw,
     & t,lego,
     & xel,yel,xpt,ypt,zpt,surfac,lgseg,ikle,iklbor,nbor,
     & xnor,ynor,npt,nelem,neleb,nelmax,nelebx,
     & ielm1,lv,msk,maskel,mesh,dim1t,nelbor,nulone,asspar)
!
!***********************************************************************
! BIEF   V7P1
!***********************************************************************
!
!brief    COMPUTES VECTORS.
!+
!+            THE VECTOR IS IDENTIFIED BY THE FORMULATION IN
!+                THE CHARACTER STRING 'FORMUL'.
!code
!+  MEANING OF IELM1
!+
!+  TYPE OF ELEMENT      NUMBER OF POINTS
!+
!+  10 : P0 TRIANGLE            1
!+  11 : P1 TRIANGLE            3
!+  12 : QUASI-BUBBLE TRIANGLE  4
!+  13 : P2 TRIANGLE            6
!+
!+  20 : Q0 QUADRILATERAL       1
!+  21 : Q1 QUADRILATERAL       4
!+
!+  40 : TELEMAC-3D P0 PRISMS   1
!+  41 : TELEMAC-3D P1 PRISMS   6
!
!history  JM HERVOUET (LNHE)
!+        16/07/07
!+        V6P0
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET (LNHE)
!+        06/12/2011
!+        V6P2
!+   Call of VC13TT modified.
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        07/01/2013
!+        V6P3
!+   Arguments XPT, YPT and ZPT added, various XEL, YEL and ZEL changed
!+   into XPT, etc. in the calls to 3D vectors.
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        09/05/2014
!+        V7P0
!+   Adaptation to assembly with I8 integers.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        12/05/2014
!+        V7P0
!+   Discontinuous elements better treated: new types 15, 16 and 17 for
!+   discontinuous linear, quasi-bubble, and quadratic, rather than
!+   using component DIMDISC=11, 12 or 13.
!
!history  A. JOLY (EDF LAB, LNHE)
!+        27/08/2015
!+        V7P1
!+   Adding vc05aa.
!
!history  R.NHEILI (Univerte de Perpignan, DALI)
!+        24/02/2016
!+        V7P3
!+        ADD MODASS 3
!
!history  A. BOURGOIN (EDF LAB, LNHE)
!+        01/08/2016
!+        V7P
!+   Adding vc17aa et vc18aa for spalart allmalas.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DIM1T          |-->| FIRST DIMENSION OF T (NELMAX OR NELEBX)
!| F              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| FORMUL         |-->| STRING WITH THE FORMULA DESCRIBING THE VECTOR
!| G              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| H              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| IELM1          |-->| TYPE OF ELEMENT
!| IKLBOR         |-->| CONNECTIVITY TABLE OF BOUNDARY ELEMENTS.
!| IKLE           |-->| CONNECTIVITY TABLE.
!| LEGO           |-->| IF YES : THE VECTOR WILL BE ASSEMBLED
!| LGSEG          |-->| LENGTH OF BOUNDARY SEGMENTS
!| LV             |-->| VECTOR LENGTH OF THE MACHINE
!| MASKEL         |-->| MASKING OF ELEMENTS
!|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
!| MESH           |-->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS
!| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
!| NELBOR         |-->| ADJACENT ELEMENT NUMBER
!| NELEB          |-->| NUMBER OF BOUNDARY ELEMENTS
!| NELEBX         |-->| MAXIMUM NUMBER OF BOUNDARY ELEMENTS
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELEM_AS       |-->| NUMBER OF ELEMENTS FOR THE CALL TO ASSVEC
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NELMAX_AS      |-->| MAXIMUM NUMBER OF ELEMENTS FOR THE CALL TO ASSVEC
!| NPT            |-->| NUMBER OF POINTS OF VECTOR.
!| NULONE         |-->| LOCAL NUMBERING OF BOUNDARY ELEMENT IN ADJACENT
!|                |   | ELEMENT.
!| OP             |-->| OPERATION TO BE DONE (SEE ABOVE)
!| SF             |-->| BIEF_OBJ STRUCTURE OF F
!| SG             |-->| BIEF_OBJ STRUCTURE OF G
!| SH             |-->| BIEF_OBJ STRUCTURE OF H
!| SU             |-->| BIEF_OBJ STRUCTURE OF U
!| SV             |-->| BIEF_OBJ STRUCTURE OF V
!| SW             |-->| BIEF_OBJ STRUCTURE OF W
!| SURFAC         |-->| AREA OF TRIANGLES
!| T              |-->| WORK ARRAY WITH THE NON ASSEMBLED VECTOR
!| U              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| V              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| W              |-->| FUNCTION USED IN THE VECTOR FORMULA
!| VEC            |<->| RESULTING VECTOR
!| XEL            |-->| ABSCISSAE OF POINTS IN THE MESH, PER ELEMENT
!| YEL            |-->| ORDINATES OF POINTS IN THE MESH, PER ELEMENT
!| XMUL           |-->| MULTIPLICATION COEFFICIENT
!| XNOR           |-->| X-COMPONENT OF NORMAL VECTOR
!| YNOR           |-->| Y-COMPONENT OF NORMAL VECTOR
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_vectos => vectos
      USE declarations_telemac, ONLY : modass
      USE declarations_special
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NELMAX,NPT,NELEM,IELM1,LV,NELEB,NELEBX
      INTEGER, INTENT(IN) :: DIM1T
      INTEGER, INTENT(IN) :: IKLE(nelmax,*),NBOR(*),IKLBOR(nelebx,*)
      INTEGER, INTENT(IN) :: NELBOR(nelebx),NULONE(nelebx,*)
!
      DOUBLE PRECISION, INTENT(IN)    :: SURFAC(nelmax),LGSEG(neleb)
      DOUBLE PRECISION, INTENT(IN)    :: XEL(*),YEL(*)
      DOUBLE PRECISION, INTENT(IN)    :: XPT(*),YPT(*),ZPT(*)
      DOUBLE PRECISION, INTENT(IN)    :: XNOR(*),YNOR(*)
      DOUBLE PRECISION, INTENT(INOUT) :: T(dim1t,*),VEC(*)
      DOUBLE PRECISION, INTENT(IN)    :: XMUL,MASKEL(*)
!
!     STRUCTURES OF FUNCTIONS F, G, H, U, V, W AND REAL DATA
!
      DOUBLE PRECISION, INTENT(IN) :: F(*),G(*),H(*),U(*),V(*),W(*)
!
      TYPE(bief_obj), INTENT(IN)     :: SF,SG,SH,SU,SV,SW
      TYPE(bief_obj), INTENT(INOUT)  :: SVEC
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
      LOGICAL, INTENT(IN) :: MSK,LEGO,ASSPAR
!
      CHARACTER(LEN=16), INTENT(IN) :: FORMUL
      CHARACTER(LEN=1), INTENT(IN) ::  OP
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER ICOORD,IDP,NDP,I
      LOGICAL INIT,SPECAD
      INTEGER NELEM2
!
      INTEGER(KIND=K8), POINTER :: PWI8(:)
!
      DOUBLE PRECISION QT
      INTEGER NPLAN
      IF(mesh%DIM1.EQ.2) THEN
        nplan=1
      ELSE
        nplan=bief_nbpts(41,mesh)/bief_nbpts(11,mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     CHECKS THE TYPE OF VECTOR
!
!=======================================================================
!     MASS MATRIX VECTOR X VECTOR F EQUALS 1
!=======================================================================
!
      IF(formul(1:16).EQ.'MASBAS          '.OR.
     &   formul(1:16).EQ.'MASBAS2         '     ) THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc00aa(xmul,surfac,nelem,nelmax,t(1,1),t(1,2),t(1,3))
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
        ELSEIF(ielm1.EQ.12) THEN
!
          CALL vc00bb(xmul,surfac,nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3),t(1,4))
!
!
!-----------------------------------------------------------------------
!
!       ELEMENT P2 TRIANGLE
!
        ELSEIF(ielm1.EQ.13) THEN
!
          CALL vc00cc(xmul,surfac,nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3),t(1,4),t(1,5),t(1,6))
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC00OO(XMUL,SURFAC,NELEM,NELMAX,T(1,1),T(1,2))
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 PRISM
!
        ELSEIF(ielm1.EQ.41) THEN
!
          CALL vc00pp(xmul,zpt,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6),formul)
!
!-----------------------------------------------------------------------
!
!       ELEMENT TETRAHEDRON T1
!
        ELSEIF(ielm1.EQ.31.OR.ielm1.EQ.51) THEN
!
          CALL vc00tt(xmul,xpt,ypt,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3),t(1,4),formul,
     &                bief_nbpts(11,mesh),bief_nbpts(10,mesh),ielm1)
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE (SIDE FACES OF A MESH OF PRISMS
!                            SPLIT IN TETRAHEDRONS)
!
        ELSEIF(ielm1.EQ.61) THEN
!
          CALL vc00ft(xmul,xpt,ypt,zpt,
     &                iklbor(1,1),iklbor(1,2),iklbor(1,3),nbor,
     &                neleb,nelebx,t(1,1),t(1,2),t(1,3))
!
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 QUADRILATERAL
!
        ELSEIF(ielm1.EQ.71) THEN
!
!         FOR VERTICAL RECTANGULAR SIDES OF THE PRISMS
          CALL vc00ff(xmul,xel,yel,zpt,
     &                iklbor(1,1),iklbor(1,2),
     &                iklbor(1,3),iklbor(1,4),nbor,
     &                neleb,nelebx,t(1,1),t(1,2),t(1,3),t(1,4),
     &                nelbor,nulone,nelmax)
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     MASS MATRIX VECTOR X VECTOR F
!=======================================================================
!
      ELSEIF(formul(1:6).EQ.'MASVEC') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc01aa(xmul,sf,f,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3) )
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE (SIDE FACES OF A MESH OF PRISMS
!                            SPLIT IN TETRAHEDRONS)
!
        ELSEIF(ielm1.EQ.61.OR.ielm1.EQ.81) THEN
!
          IF(formul(7:7).NE.'2') THEN
!
            CALL vc01ft(xmul,sf,f,xpt,ypt,zpt,
     &                  iklbor(1,1),iklbor(1,2),
     &                  iklbor(1,3),nbor,
     &                  neleb,nelebx,t(1,1),t(1,2),t(1,3))
          ELSE
!
            CALL vc01ft2(xmul,sf,f,sg,g,xpt,ypt,zpt,
     &                   iklbor(1,1),iklbor(1,2),
     &                   iklbor(1,3),nbor,
     &                   neleb,nelebx,t(1,1),t(1,2),t(1,3))
          ENDIF
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
        ELSEIF(ielm1.EQ.12) THEN
!
          CALL vc01bb(xmul,sf,f,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3),t(1,4) )
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
        ELSEIF(ielm1.EQ.1) THEN
!
          CALL vc01oo(xmul,sf,f,lgseg,
     &                iklbor(1,1),iklbor(1,2),nbor,neleb,nelebx,
     &                t(1,1),t(1,2)  )
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 QUADRILATERAL (IN A 3D MESH OF PRISMS)
!
        ELSEIF(ielm1.EQ.71) THEN
!
!         FOR VERTICAL RECTANGULAR SIDES OF THE PRISMS
          CALL vc01ff(xmul,sf,f,xel,yel,zpt,
     &                iklbor(1,1),iklbor(1,2),
     &                iklbor(1,3),iklbor(1,4),nbor,
     &                neleb,nelebx,t(1,1),t(1,2),t(1,3),t(1,4),
     &                nelbor,nulone,nelmax)
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 PRISM
!
        ELSEIF(ielm1.EQ.41) THEN
!
          CALL vc01pp(xmul,sf,f,zpt,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6))
!
!-----------------------------------------------------------------------
!
!       ELEMENT T1 TETRAHEDRON
!
        ELSEIF(ielm1.EQ.31.OR.ielm1.EQ.51) THEN
!
          CALL vc01tt(xmul,sf,f,xpt,ypt,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3),t(1,4))
!
!-----------------------------------------------------------------------
!
!       ELEMENT T0 TETRAHEDRON
!
        ELSEIF(ielm1.EQ.30) THEN
!
          CALL vc01tt0(xmul,sf,f,xpt,ypt,zpt,
     &                ikle(:,1),ikle(:,2),ikle(:,3),ikle(:,4),
     &                nelem,nelmax,vec)
!
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!                 -    ---->   --->
!     VECTOR:     F  . GRAD(U).GRAD(PSI)
!
!     F DIFFUSION TENSOR WITH DIAGONAL COMPONENTS F, G, H
!
!     EQUIVALENT OF MATDIF * U
!
!=======================================================================
!
      ELSEIF(formul(1:6).EQ.'VECDIF') THEN
!
!       ELEMENT SEGMENT P1
!
        IF(ielm1.EQ.41.AND.formul(8:8).EQ.'*') THEN
!
          CALL vc02pp_star(xmul,sf,sg,sh,su,f,g,h,u,xel,yel,zpt,surfac,
     &                     ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                     ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),
     &                     t(1,2),t(1,3),t(1,4),t(1,5),t(1,6),formul)
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR K GRAD(PSI) U.GRAD(F)
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'SUPG            ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc03aa(xmul,sf,sg,sh,su,sv,f,g,h,u,v,
     &                xel,yel,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3) )
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
        ELSEIF(ielm1.EQ.12) THEN
!
          CALL vc03bb(xmul,sf,sg,sh,su,sv,f,g,h,u,v,
     &                xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3),t(1,4) )
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC03OO(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NELEM,NELMAX,
!    &                T(1,1),T(1,2)  )
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR U GRAD(PSI)
!=======================================================================
!
      ELSEIF(formul(1:6).EQ.'VGRADP') THEN
!
        specad = .false.
        IF(formul(8:8).EQ.'2') specad = .true.
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc04aa(xmul,su,sv,u,v,xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3),specad)
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 PRISM
!
        ELSEIF(ielm1.EQ.41) THEN
!
          nelem2 = bief_nbpts(10,mesh)
          CALL vc04pp(xmul,su,sv,sw,u,v,w,f,h(1:nelem2*2),
     &                xel,yel,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6),specad,formul,
!                     THIS IS TO RETRIEVE THE NUMBER OF TRIANGLES
     &                nelem2)
!
!-----------------------------------------------------------------------
!
!       ELEMENT T1 TETRAHEDRON
!
        ELSEIF(ielm1.EQ.31.OR.ielm1.EQ.51) THEN
!
          CALL vc04tt(xmul,su,sv,sw,u,v,w,f,h,xpt,ypt,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),formul,specad,
!                     THIS IS TO RETRIEVE THE NUMBER OF 2D POINTS
!                     AND TRIANGLES
     &                bief_nbpts(11,mesh),bief_nbpts(10,mesh))
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR U.N    (N VECTOR NORMAL TO THE ELEMENT)
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'FLUBOR          ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE FOR VERTICAL SIDES OF THE PRISMS SPLIT
!       IN TETRAHEDRONS
!
        IF(ielm1.EQ.61) THEN
!
          CALL vc05ft(xmul,su,sv,u,v,xpt,ypt,zpt,
     &                iklbor(1,1),iklbor(1,2),iklbor(1,3),nbor,
     &                neleb,nelebx,t(1,1),t(1,2),t(1,3))
!
!       ELEMENT P1 QUADRILATERAL FOR VERTICAL SIDES OF THE PRISMS
!
        ELSEIF(ielm1.EQ.71) THEN
!
          CALL vc05ff(xmul,su,sv,u,v,xel,yel,zpt,
     &                iklbor(1,1),iklbor(1,2),iklbor(1,3),iklbor(1,4),
     &                nbor,neleb,nelebx,t(1,1),t(1,2),t(1,3),t(1,4),
     &                nelbor,nulone,nelmax)

! !
! !       ELEMENT P1 PRISM
! !
!         ELSEIF(IELM1.EQ.41) THEN
!
!           CALL VC05PP(XMUL,SF,F,ZPT,SURFAC,
!      &                IKLE(1,1),IKLE(1,2),IKLE(1,3),
!      &                IKLE(1,4),IKLE(1,5),IKLE(1,6),
!      &                NELEM,NELMAX,T(1,1),T(1,2),T(1,3),
!      &                T(1,4),T(1,5),T(1,6) )
!
!
!       ELEMENT P1 TRIANGLES
!
        ELSEIF(ielm1.EQ.11) THEN
!
          CALL vc05aa(xmul,sf,f,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3) )
!
!       ELEMENT LINEAR SEGMENT
!
        ELSEIF(ielm1.EQ.1) THEN
!
          CALL vc05oo(xmul,su,sv,u,v,xnor,ynor,lgseg,
     &                iklbor,nbor,neleb,nelebx,t(1,1),t(1,2))
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR U GRAD(F)
!=======================================================================
!
      ELSEIF(formul(1:13).EQ.'VGRADF       ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc08aa(xmul,sf,su,sv,f,u,v,xel,yel,ikle,
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3) , formul )
!
!-----------------------------------------------------------------------
!
!       ELEMENT P2 TRIANGLE
!
        ELSEIF(ielm1.EQ.13) THEN
!
          CALL vc08cc(xmul,sf,su,sv,f,u,v,xel,yel,
     &         ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &         ikle(1,5),ikle(1,6),nelem,nelmax,
     &         t(1,1),t(1,2),t(1,3),t(1,4),t(1,5),t(1,6),formul)
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
        ELSEIF(ielm1.EQ.12) THEN
!
          CALL vc08bb(xmul,sf,su,sv,f,u,v,xel,yel,
     &         ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &         nelem,nelmax,t(1,1),t(1,2),t(1,3),t(1,4),formul)
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC08OO(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NELEM,NELMAX,
!    &                T(1,1),T(1,2)  )
!
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 PRISM
!
        ELSEIF(ielm1.EQ.41) THEN
!
          CALL vc08pp(xmul,sf,su,sv,sw,f,u,v,w,xel,yel,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6))
!
!-----------------------------------------------------------------------
!
!       ELEMENT T1 TETRAHEDRON
!
        ELSEIF(ielm1.EQ.31.OR.ielm1.EQ.51) THEN
!
          CALL vc08tt(xmul,sf,su,sv,sw,f,u,v,w,xpt,ypt,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4))
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR U GRAD(F) 2
!=======================================================================
!
      ELSEIF(formul(1:13).EQ.'VGRADF2      ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 PRISM
!
        IF(ielm1.EQ.41) THEN
!
          CALL vc18pp(xmul,sf,su,sv,f,u,v,xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3))
!
!-----------------------------------------------------------------------
!
!       ELEMENT T1 TETRAHEDRON
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR Q GRAD(F)
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'QGRADF          ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc09aa(xmul,sf,sg,su,sv,f,g,u,v,xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3) )
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC09OO(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NELEM,NELMAX,
!    &                T(1,1),T(1,2)  )
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR F U.N    (N VECTOR NORMAL TO THE ELEMENT)
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'FLUBDF          ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
!       IF(IELM1.EQ.11) THEN
!
!         CALL VC10AA(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),IKLE(1,3),NELEM,NELMAX,
!    &                T(1,1)   ,T(1,2)   ,T(1,3))
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
        IF(ielm1.EQ.1) THEN
!
          CALL vc10oo(xmul,sf,su,sv,f,u,v,xnor,ynor,lgseg,
     &                iklbor,nbor,neleb,nelebx,t(1,1),t(1,2))
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR G GRADIENT(F)
!=======================================================================
!
      ELSEIF(formul(1:15).EQ.'GGRADF         ') THEN
!
!     CHARACTER 16 GIVES THE SELECTED COORDINATE
!
      IF(formul(16:16).EQ.'X') THEN
        icoord=1
      ELSEIF(formul(16:16).EQ.'Y') THEN
        icoord=2
      ELSEIF(formul(16:16).EQ.'Z') THEN
        icoord=3
      ENDIF
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
!         CHECKS IF G IS DISCONTINUOUS P1
!
          IF(sg%ELM.EQ.15) THEN
!
            CALL vc11aa2(xmul,sf,sg,sh,f,g,h,xel,yel,
     &                   ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                   t(1,1),t(1,2),t(1,3) , icoord )
!
!         CLASSICAL CASE: G IS P1
!
          ELSE
!
            CALL vc11aa(xmul,sf,sg,f,g,xel,yel,
     &                  ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                  t(1,1),t(1,2),t(1,3) , icoord )
!
          ENDIF
!
!-----------------------------------------------------------------------
!
        ELSEIF(ielm1.EQ.12) THEN
!
          CALL vc11bb(xmul,sf,sg,f,g,xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3),t(1,4) , icoord )
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC11OO(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,XNOR,YNOR,ZNOR,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NBOR,NELEM,NELMAX,
!    &                T(1,1),T(1,2)  )
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 PRISM
!
        ELSEIF(ielm1.EQ.41) THEN
!
          CALL vc11pp(xmul,sf,sg,f,g,
     &                xel,yel,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6),icoord)
!
!-----------------------------------------------------------------------
!
!       ELEMENT T1 TETRAHEDRON
!
        ELSEIF(ielm1.EQ.31.OR.ielm1.EQ.51) THEN
!
          CALL vc11tt(xmul,sf,sg,f,g,
     &                xpt,ypt,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),icoord)
!
        ELSEIF(ielm1.EQ.30) THEN
!
! COMMENT JP RENAUD 21/09/2005
! BEWARE: THIS CALL CREATES A P0 VECTOR. THEREFORE
! THERE IS NO NEED TO "ASSEMBLE" IT AFTERWARDS. NOTE
! THAT VEC ITSELF (AND NOT T) IS GIVEN TO VC11TT0 AND
! THAT LEGO IS ALSO SET TO .FALSE. TO AVOID
! CALLING ASSVEC AT THE END OF THE SUBROUTINE.
!
          CALL vc11tt0(xmul,sf,sg,f,g,
     &                 xpt,ypt,zpt,
     &                 ikle(1:nelem,1),ikle(1:nelem,2),
     &                 ikle(1:nelem,3),ikle(1:nelem,4),
     &                 nelem,mesh%NPOIN,
     &                 vec,icoord)
!
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR GRADIENT(F)
!=======================================================================
!
      ELSEIF(formul(1:5).EQ.'GRADF') THEN
!
!     CHARACTER 16 GIVES THE SELECTED COORDINATE
!
      IF(formul(16:16).EQ.'X') THEN
        icoord=1
      ELSEIF(formul(16:16).EQ.'Y') THEN
        icoord=2
      ELSEIF(formul(16:16).EQ.'Z') THEN
        icoord=3
      ENDIF
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc13aa(xmul,sf,f,xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,1),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3),icoord)
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
        ELSEIF(ielm1.EQ.12) THEN
!
          CALL vc13bb(xmul,sf,f,xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),nelem,
     &                nelmax,t(1,1),t(1,2),t(1,3),t(1,4),icoord)
!
!-----------------------------------------------------------------------
!
!       ELEMENT P2 TRIANGLE
!
        ELSEIF(ielm1.EQ.13) THEN
!
          CALL vc13cc(xmul,sf,f,xel,yel,
     &                 ikle(1,1),ikle(1,2),ikle(1,3),
     &                 ikle(1,4),ikle(1,5),ikle(1,6),
     &                 nelem,nelmax,
     &                 t(1,1),t(1,2),t(1,3),
     &                 t(1,4),t(1,5),t(1,6),icoord)
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC13OO(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,XNOR,YNOR,ZNOR,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NBOR,NELEM,NELMAX,
!    &                T(1,1),T(1,2)  )
!
!
!
!-----------------------------------------------------------------------
!
!       ELEMENT P2 SEGMENT
!
!       ELSEIF(IELM1.EQ.2) THEN
!
!         CALL VC13OC(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,XNOR,YNOR,ZNOR,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NBOR,NELEM,NELMAX,
!    &                T(1,1),T(1,2),T(1,3)  )
!
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 PRISM
!
        ELSEIF(ielm1.EQ.41) THEN
!
          IF(formul(1:15).EQ.'GRADF(X,Y)     ') THEN
!         SIMPLIFIED FORMULATION FOR EFFICIENCY AND ACCURACY
          CALL vc13pp2(xmul,sf,f,xel,yel,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6),icoord)
          ELSEIF( formul(8:15).EQ.'        ') THEN
!                 FORMUL(6:7) IS LEFT FOR OPTIONS
          CALL vc13pp(xmul,sf,f,xel,yel,zpt,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6),icoord,formul)
!
          ELSE
            WRITE(lu,1001) formul
            CALL plante(1)
            stop
          ENDIF
!
!-----------------------------------------------------------------------
!
!       ELEMENT T1 TETRAHEDRON
!
        ELSEIF(ielm1.EQ.31.OR.ielm1.EQ.51) THEN
!
          IF(formul(1:5).EQ.'GRADF'.AND.
     &       formul(8:15).EQ.'        ') THEN
!           FORMUL(6:7) IS LEFT FOR OPTIONS
            CALL vc13tt(xmul,sf,f,xpt,ypt,zpt,
     &                  ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                  nelem,nelmax,t(1,1),t(1,2),
     &                  t(1,3),t(1,4),icoord,formul)
!
          ELSE
            WRITE(lu,1001) formul
            CALL plante(1)
            stop
          ENDIF
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     TURBULENT PRODUCTION VECTOR
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'PRODF           ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc14aa(xmul,su,sv,u,v,xel,yel,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3))
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
!       ELSEIF(IELM1.EQ.12) THEN
!
!         CALL VC14BB(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),IKLE(1,3),IKLE(1,4),NELEM,
!    &                NELMAX,T(1,1),T(1,2),T(1,3),T(1,4))
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC14OO(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,XNOR,YNOR,ZNOR,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NBOR,NELEM,NELMAX,
!    &                T(1,1),T(1,2))
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     SPALART ALLMARAS PRODUCTION VECTOR
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'PRSAF           ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
          CALL vc17aa(xmul,su,sv,u,v,xel,yel,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3))
!
!
!       ELEMENT P1 PRISM
!
        ELSEIF(ielm1.EQ.41) THEN
          CALL vc17pp(xmul,surfac,su,sv,sw,u,v,w,xel,yel,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6))

!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     SPALART ALLMARAS DESTRUCTION VECTOR
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'TRSAF           ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
          CALL vc21aa(xmul,su,sv,u,v,xel,yel,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3))
        ELSEIF(ielm1.EQ.41) THEN
          CALL vc21pp(xmul,sf,f,xel,yel,zpt,surfac,ikle(1,1),ikle(1,2),
     &                ikle(1,3),ikle(1,4),ikle(1,5),ikle(1,6),nelem,
     &                nelmax,t(1,1),t(1,2),t(1,3),t(1,4),t(1,5),t(1,6))

!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          IF (lng.EQ.1) WRITE(lu,1000) formul
          IF (lng.EQ.2) WRITE(lu,1001) formul
          IF (lng.EQ.1) WRITE(lu,2000) ielm1
          IF (lng.EQ.2) WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF

      ELSEIF(formul(1:16).EQ.'STRAIN          ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
          CALL vc20aa(xmul,su,sv,u,v,xel,yel,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3))
        ELSEIF(ielm1.EQ.41) THEN
          CALL vc20pp(xmul,surfac,su,sv,sw,u,v,w,xel,yel,zpt,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                ikle(1,5),ikle(1,6),nelem,nelmax,t(1,1),t(1,2),
     &                t(1,3),t(1,4),t(1,5),t(1,6))

!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR DIV(HU)
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'DIVQ            ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc15aa(xmul,sf,su,sv,f,u,v,xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),ikle(1,4),
     &                nelem,nelmax,t(1,1),t(1,2),t(1,3))
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
!       ELSEIF(IELM1.EQ.12) THEN
!
!         CALL VC15BB(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &         XEL,YEL,ZEL,SURFAC,
!    &         IKLE(1,1),IKLE(1,2),IKLE(1,3),IKLE(1,4),
!    &         NELEM,NELMAX,T(1,1),T(1,2),T(1,3),T(1,4),FORMUL)
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC15OO(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NELEM,NELMAX,
!    &                T(1,1),T(1,2)  )
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     VECTOR K.GRAD(PSI) DIV(U)
!=======================================================================
!
      ELSEIF(formul(1:16).EQ.'SUPGDIVU        ') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc16aa(xmul,sf,sg,su,sv,f,g,u,v,xel,yel,surfac,
     &                ikle(1,1),ikle(1,2),ikle(1,3),
     &                nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3) )
!
!-----------------------------------------------------------------------
!
!       ELEMENT QUASI-BUBBLE TRIANGLE
!
!       ELSEIF(IELM1.EQ.12) THEN
!
!         CALL VC16BB(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &         XEL,YEL,ZEL,SURFAC,
!    &         IKLE(1,1),IKLE(1,2),IKLE(1,3),IKLE(1,4),
!    &         NELEM,NELMAX,T(1,1),T(1,2),T(1,3),T(1,4),FORMUL)
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 SEGMENT
!
!       ELSEIF(IELM1.EQ.1) THEN
!
!         CALL VC16OO(XMUL,SF,SG,SH,SU,SV,SW,F,G,H,U,V,W,
!    &                XEL,YEL,ZEL,SURFAC,
!    &                IKLE(1,1),IKLE(1,2),NELEM,NELMAX,
!    &                T(1,1),T(1,2)  )
!
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!
!=======================================================================
!     VECTOR H U.GRAD(PSI)
!=======================================================================
!
      ELSEIF(formul(1:7).EQ.'HUGRADP') THEN
!
!-----------------------------------------------------------------------
!
!       ELEMENT P1 TRIANGLE
!
        IF(ielm1.EQ.11) THEN
!
          CALL vc19aa(xmul,sf,sg,sh,su,sv,f,g,h,u,v,
     &                xel,yel,
     &                ikle(1,1),ikle(1,2),ikle(1,3),nelem,nelmax,
     &                t(1,1),t(1,2),t(1,3),formul)
!
!-----------------------------------------------------------------------
!       OTHER
!-----------------------------------------------------------------------
!
!       ELSEIF
!
!-----------------------------------------------------------------------
!       ERROR ON THE ELEMENT TYPE
!-----------------------------------------------------------------------
!
        ELSE
!
          WRITE(lu,1001) formul
          WRITE(lu,2001) ielm1
          CALL plante(1)
          stop
!
        ENDIF
!
!=======================================================================
!     OTHER VECTOR
!=======================================================================
!
!
!=======================================================================
!     ERROR: TYPE OF VECTOR NOT CODED UP
!=======================================================================
!
      ELSE
        WRITE(lu,1001) formul
        CALL plante(1)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!  POSSIBLE ASSEMBLY OF THE VECTOR
!
      IF(lego) THEN
!
        IF(op(1:1).EQ.'=') THEN
          init = .true.
        ELSEIF(op(1:1).EQ.'+') THEN
          init = .false.
        ELSE
          WRITE(lu,3001) op
3001      FORMAT(1x,'VECTOS (BIEF) : OP NOT RECOGNISED:',a1)
          CALL plante(1)
          stop
        ENDIF
!
        ndp=bief_nbpel(ielm1,mesh)
!
        IF(dimens(ielm1).EQ.mesh%DIM1) THEN
          IF(modass.EQ.1) THEN
            CALL assvec(vec, ikle, npt ,nelem,nelmax,
     &                  t,init,lv,msk,maskel,ndp)
!
          ELSEIF(modass.EQ.3 ) THEN
            CALL assvec(vec, ikle, npt ,nelem,nelmax,
     &                  t,init,lv,msk,maskel,ndp,svec%E)
!
          ELSEIF(modass.EQ.2) THEN
            CALL double_to_integer(t,mesh%WI8,ndp*nelem,qt,mesh%MXELVS)
            DO i=1,npt
              mesh%TI8(i)=0
            ENDDO
            DO idp = 1 , ndp
              pwi8=>mesh%WI8(1+(idp-1)*nelmax:idp*nelmax)
              CALL assve1i8(mesh%TI8,ikle(1,idp),pwi8,nelem)
            ENDDO
            IF(asspar) THEN
              CALL parcom2i8(mesh%TI8,mesh%TI8,mesh%TI8,
     &                       npt,nplan,2,1,mesh)
            ENDIF
            IF(init) THEN
              CALL integer_to_double('=',mesh%TI8,vec,npt,qt)
            ELSE
              CALL integer_to_double('+',mesh%TI8,vec,npt,qt)
            ENDIF
          ELSE
            WRITE(lu,*) 'VECTOS: MODASS=',modass,' UNEXPECTED CASE'
            CALL plante(1)
            stop
          ENDIF
        ELSEIF(neleb.GT.0) THEN
          CALL assvec(vec, iklbor, npt ,neleb,nelebx,
     &                t,init,lv,msk,maskel,ndp)
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!     OPTIONAL ASSEMBLY IN PARALLEL (ASSPAR DONE IN VECTOR)
!-----------------------------------------------------------------------
!
      IF(asspar.AND.modass.EQ.1) CALL parcom(svec,2,mesh)
      IF(ncsize.GT.1) THEN
        IF(asspar .AND. modass.EQ.3) THEN
          CALL parcom_comp(svec,svec%E,2,mesh)
        ENDIF
      ENDIF
!
!     CompSum:Correct the VEC with accumuletad fperrors
!
      IF(asspar.AND. modass .EQ. 3) THEN
        DO i = 1 , mesh%NPOIN
          vec(i)= vec(i)+svec%E(i)
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
1000  FORMAT(1x,'VECTOS (BIEF) : VECTEUR NON PREVU : ',a16)
1001  FORMAT(1x,'VECTOS (BIEF) : VECTOR NOT IMPLEMENTED:',a16)
2000  FORMAT(1x,'                POUR IELM1 = ',1i6)
2001  FORMAT(1x,'                FOR IELM1 = ',1i6)
!
!-----------------------------------------------------------------------
!
      RETURN
      END