almesh.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE almesh
!                   *****************
!
     &(mesh,nom,ielm,spheri,cfg,fformat,nfic,equa,refine,nplan,npmax,
     & nptfrx,nelmax,projection,lati0,longi0,convergence,rlevel)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    ALLOCATES A BIEF_MESH MESH STRUCTURE.
!
!history  J-M HERVOUET (LNHE)
!+        05/02/2010
!+        V6P0
!+   EDGE-BASED STUCTURES ALWAYS ALLOCATED
!
!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)
!+        28/11/2011
!+        V6P2
!+   Calls to eleb3d and eleb3dt changed. Calls of READGEO3 and CPIKLE2
!+   and CPIKLE3 swapped (now KNOLG used in CPIKLE3).
!
!history  J-M HERVOUET (LNHE)
!+        20/07/2012
!+        V6P2
!+   Finding the original number of nodes in parallel, and completing
!+   KNOLG for upper planes in 3D.
!
!history  J-M HERVOUET (LNHE)
!+        07/02/2013
!+        V6P3
!+   Argument NELMAX removed in call to SEGBOR.
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        11/03/2013
!+        V6P3
!+   Dimension of LIMVOI now set to (11,2).
!history  R.NHEILI (Univerte de Perpignan, DALI)
!+        24/02/2016
!+        V7P3
!     ALLOCATE BUFER FOR ERRORS
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        15/05/2015
!+        V7P1
!+   Call to a new subroutine that checks the mesh (coordinates...).
!
!history Y AUDOUIN (LNHE)
!+       25/05/2015
!+       V7P0
!+       Modification to comply with the hermes module
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        30/07/2015
!+        V7P1
!+   Introducing MESH%IFAC, to replace MESH%FAC.
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        08/09/2017
!+        V7P3
!+   New optional argument RLEVEL added. Parts of the subroutine not
!+   executed when RLEVEL > 0 and NPOIN, NPTFR, NELEM retrieved in
!+   structure MESH.
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        23/09/2017
!+        V7P3
!+   In the convergence procedure, exact formulas are now used for
!+   dimensioning the arrays. Moreover IKLES and IPOBO, and a number
!+   of other arrays need not be overdimensioned.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CONVERGENCE    |-->| IF YES, WE ARE IN A REFINEMENT PROCEDURE
!| EQUA           |-->| NAME IN 20 CHARACTERS TO ENABLE DIFFERENT
!|                |   | OPTIONS. OPTIONS ARE:
!|                |   | "SAINT-VENANT EF"
!|                |   | "SAINT-VENANT VF"
!|                |   | "BOUSSINESQ"
!| IELM           |-->| ELEMENT TYPE WITH THE LARGET NUMBER OF DEGREES
!|                |   | OF FREEDOM THAT WILL BE USED
!| LATI0          |-->| LATITUDE OF ORIGIN POINT
!| LONGI0         |-->| LONGITUDE OF ORIGIN POINT
!| MESH           |-->| MESH STRUCTURE TO BE ALLOCATED
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS IN THE MESH
!| NFIC           |-->| LOGICAL UNIT WHERE TO READ THE MESH
!| NOM            |-->| NAME OF THE MESH
!| NPLAN          |-->| NUMBER OF PLANES (OPTIONAL,3D MESHES OF PRISMS)
!| NPMAX          |-->| MAXIMUM NUMBER OF POINTS IN THE MESH
!| NPTFRX         |-->| MAXIMUM NUMBER OF BOUNDARY NODES
!| PROJECTION     |<->| SPATIAL PROJECTION TYPE
!| REFINE         |-->| NUMBER OF REFINEMENT LEVELS FOR CONVERGENCE
!| RLEVEL         |-->| REFINEMENT LEVEL FOR CONVERGENCE
!| SPHERI         |-->| LOGICAL, IF YES : SPHERICAL COORDINATES
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_almesh => almesh
      USE declarations_telemac
      USE interface_hermes
!
      USE declarations_special
      USE interface_parallel, ONLY : p_max
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(bief_mesh)  , INTENT(INOUT)           :: MESH
      INTEGER          , INTENT(IN)              :: IELM
      INTEGER          , INTENT(IN)              :: NFIC
      CHARACTER(LEN=8) , INTENT(IN)              :: FFORMAT
      LOGICAL          , INTENT(IN)              :: SPHERI
      CHARACTER(LEN=6) , INTENT(IN)              :: NOM
      CHARACTER(LEN=20), INTENT(IN)              :: EQUA
      INTEGER          , INTENT(INOUT)           :: CFG(2)
      INTEGER          , INTENT(IN)              :: REFINE
      INTEGER          , INTENT(IN),    OPTIONAL :: NPLAN
      INTEGER          , INTENT(IN),    OPTIONAL :: NPMAX
      INTEGER          , INTENT(IN),    OPTIONAL :: NPTFRX
      INTEGER          , INTENT(IN),    OPTIONAL :: NELMAX,RLEVEL
      INTEGER          , INTENT(IN),    OPTIONAL :: PROJECTION
      DOUBLE PRECISION , INTENT(IN),    OPTIONAL :: LATI0,LONGI0
      LOGICAL          , INTENT(IN),    OPTIONAL :: CONVERGENCE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER D,IELM0,IELM1,STOCFG,IELB0,IELB1,NSEG,NNPMAX
      INTEGER NNPTFRX,NNELEB,ERR,NNELMAX,NNPLAN,NPOIN,NPTFR,NELEM
      INTEGER MXPTVS,MXELVS,NDP,IELEM,NSEGBOR,IPLAN,NPOIN_GLOB
      INTEGER NVAR,TYP_ELEM,TYP_BND_ELEM,IERR,REFLEVEL
      LOGICAL CONV
!
!     TEMPORARY CONNECTIVITY TABLE
!
      INTEGER, ALLOCATABLE :: IKLES(:)
!
!     TEMPORARY TABLE TO NUMBER THE BOUNDARY NODES
!
      INTEGER, ALLOCATABLE :: IPOBO(:)
!
      INTEGER IELB0V,IELB1V,I
!
      CHARACTER(LEN=80) TITLE
!
      INTEGER PROJEC
      DOUBLE PRECISION LATI,LONGI
!
!     EARTH RADIUS
      DOUBLE PRECISION, PARAMETER :: R=6.37d6
!
!-----------------------------------------------------------------------
!
      IF(PRESENT(convergence)) THEN
        conv = convergence
      ELSE
        conv = .false.
      ENDIF
      IF(PRESENT(rlevel).AND.conv) THEN
        reflevel=rlevel
      ELSE
        reflevel=0
      ENDIF
      IF(PRESENT(projection)) THEN
        projec=projection
      ELSE
        projec=1
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(PRESENT(lati0)) THEN
        lati=lati0
      ELSE
        lati=0.d0
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(PRESENT(longi0)) THEN
        longi=longi0
      ELSE
        longi=0.d0
      ENDIF
!
!-----------------------------------------------------------------------
!
!     FIRST READS THE GEOMETRY FILE TO GET NPOIN,..
!
      mesh%NAME = nom
!
!     IN PARALLEL MODE, THIS IS WHERE NPTIR IS READ
      ALLOCATE(mesh%X_ORIG)
      ALLOCATE(mesh%Y_ORIG)
!
      CALL read_mesh_info(fformat,nfic,title,nvar,npoin,typ_elem,
     &                    nelem,nptfr,nptir,ndp,nnplan,
     &                    mesh%X_ORIG,mesh%Y_ORIG,
     &                    typ_bnd_elem,nneleb)
!
!     ALLOCATES THE TEMPORARY CONNECTIVITY TABLE
!
      ALLOCATE(ikles(nelem*ndp),stat=err)
      CALL check_allocate(err,'ALMESH:IKLES')
!
!     ALLOCATES THE TEMPORARY TABLE FOR THE BOUNDARY NODES
!
      ALLOCATE(ipobo(npoin),stat=err)
      CALL check_allocate(err,'ALMESH:IPOBO')
!
!     GETS CONNECTIVITY TABLE IKLE AND READS IPOBO IF NOT IN PARALLEL
!
      CALL read_mesh_conn(fformat,nfic,npoin,typ_elem,nelem,ndp,
     &                    typ_bnd_elem,nneleb,ikles,ipobo)
!
!     CALCULATES THE MAXIMUM NUMBER OF ELEMENTS AROUND A NODE MXELVS
!     AND THE MAXIMUM NUMBER OF SURROUNDING NODES, MXPTVS
!
      CALL mxptel(mxptvs,mxelvs,ikles,ielm,
     &            npoin,nelem,ndp,ipobo,.true.)
!
      DEALLOCATE(ipobo)
!
!-----------------------------------------------------------------------
!
      IF(ncsize.GT.1) THEN
!       IN PARALLEL MODE NSEGBOR (COMPUTED IN SEGBOR) IS NOT NPTFR
!       IT INCLUDES INTERNAL SEGMENTS
        CALL segbor(nsegbor,ikles,nelem,npoin)
      ELSE
        nsegbor=nptfr
      ENDIF
!
!-----------------------------------------------------------------------
!
!     INITIALISES COMMONS DIMS AND NODES
!
      IF(PRESENT(npmax)) THEN
        nnpmax = npmax
      ELSEIF(conv) THEN
!       EXACT FORMULA FOR ELEMENTS SUCCESSIVELY SPLIT INTO 4
!       UP TO LEVEL REFINE
        nnpmax = npoin+(nelem*(4**refine-1)+nsegbor*(2**refine-1))/2
      ELSE
        nnpmax = npoin
      ENDIF
      IF(PRESENT(nptfrx)) THEN
        nnptfrx = nptfrx
      ELSEIF(conv) THEN
!       EXACT FORMULA
        nnptfrx = nsegbor*2**refine
      ELSE
        nnptfrx = nptfr
      ENDIF
      IF(PRESENT(nelmax)) THEN
        nnelmax = nelmax
      ELSEIF(conv) THEN
!       EXACT FORMULA
        nnelmax = nelem*4**refine
      ELSE
        nnelmax = nelem
      ENDIF
      IF(PRESENT(nplan)) THEN
        nnplan = nplan
      ELSE
        nnplan = 1
      ENDIF
!
!-----------------------------------------------------------------------
!
!     IF IN A REFINEMENT PROCEDURE, THE DIMENSIONS OF THE INITIAL MESH
!     ARE NOT THE CURRENT DIMENSIONS, THEY ARE ALREADY IN THE MESH
!     STRUCTURE.
!
      IF(reflevel.GT.0) THEN
        npoin=mesh%NPOIN
        nptfr=mesh%NPTFR
        nelem=mesh%NELEM
        IF(ncsize.GT.1) THEN
          nsegbor=nsegbor*2**reflevel
        ELSE
          nsegbor=nptfr
        ENDIF
      ENDIF
!
!     IN CALL ININDS, ALL VALUES ARE 2D VALUES
!     ONLY NNPLAN TELLS IF IT'S 2D OR 3D
!
!     3D TETRAHEDRONS MESH:
!     ADD THE OPTIONAL ARGUMENT NNELEB
!
      IF(reflevel.EQ.0) ALLOCATE(mesh%NDS(0:81,7))
      CALL bief_ininds(npoin,nptfr,nelem,nnpmax,nnptfrx,nnelmax,
     &                 nnplan,nsegbor,mesh%NDS,nneleb)
!
!     P0 AND P1 ELEMENTS
!
      ielm0  = 10*(ielm/10)
      ielm1  = ielm0 + 1
      d      = dimens(ielm0)
!
! BOUNDARY ELEMENTS (AT THE SURFACE AND AT THE BOTTOM FOR PRISMS)
!
      ielb0  = ielbor(ielm0,1)
      ielb1  = ielbor(ielm1,1)
!
! LATERAL BOUNDARY ELEMENTS (DIFFERENT FOR PRISMS)
!
      ielb0v = ielbor(ielm0,2)
      ielb1v = ielbor(ielm1,2)
!
!-----------------------------------------------------------------------
!
!  ALLOCATES THE ARRAYS OF REALS
!
!     COORDINATES BY ELEMENTS: XEL, YEL, ZEL
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%XEL)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%YEL)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ZEL)
!
      CALL bief_allvec(1,mesh%XEL,'XEL   ',
     &                 ielm0,bief_nbpel(ielm1,mesh),1,mesh)
      CALL bief_allvec(1,mesh%YEL,'YEL   ',
     &                 ielm0,bief_nbpel(ielm1,mesh),1,mesh)
!
      IF(d.GE.3) THEN
        CALL bief_allvec(1,mesh%ZEL,'ZEL   ',
     &                   ielm0,bief_nbpel(ielm1,mesh),1,mesh)
      ELSE
        CALL bief_allvec(1,mesh%ZEL,'ZEL   ',    0,           1,0,mesh)
      ENDIF
!
!     SURFACES OF THE ELEMENTS: SURFAC
!     JAJ CAN BE USED FOR ELEMENT VOLUMES...
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%SURFAC)
      CALL bief_allvec(1,mesh%SURFAC,'SURFAC',ielm0,1,1,mesh)
!
!     1/DET : SURDET ! NOT USED IN 3D, WHY?
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%SURDET)
      CALL bief_allvec(1,mesh%SURDET,'SURDET',ielm0,1,1,mesh)
!
!     LENGTHS OF THE SEGMENTS: LGSEG
!     CAN BE USED (IN THEORY) FOR LATERAL SURFACES IN 3D,
!     BUT THEN IELB0V INSTEAD OF IELB0! (2D CASE NOT AFFECTED)
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%LGSEG)
      CALL bief_allvec(1,mesh%LGSEG,'LGSEG ',ielb0v,1,1,mesh)
!
!     NORMALS TO THE SEGMENTS: XSGBOR, YSGBOR, ZSGBOR
! CAN BE (IN THEORY) USED FOR "NON-SIGMA" MESH FOR LATERAL NORMAL VECTORS
! PER LATERAL BOUNDARY ELEMENT, BUT THEN IELB0V INSTEAD OF IELB0!
! 2D CASE NOT AFFECTED
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%XSGBOR)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%YSGBOR)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ZSGBOR)
!     SEE NORMAB FOR MEANING OF 4 DIMENSIONS
      CALL bief_allvec(1,mesh%XSGBOR,'XSGBOR',ielb0v,4,1,mesh)
      CALL bief_allvec(1,mesh%YSGBOR,'YSGBOR',ielb0v,4,1,mesh)
      IF(d.GE.3) THEN
        CALL bief_allvec(1,mesh%ZSGBOR,'ZSGBOR',ielb0v,4,1,mesh)
      ELSE
        CALL bief_allvec(1,mesh%ZSGBOR,'ZSGBOR',     0,4,0,mesh)
      ENDIF
!
!     NORMALS AT THE NODES: XNEBOR, YNEBOR, ZNEBOR
!
! IN 3D THEY ARE NORMAL VECTORS AT THE BOTTOM
! SO THAT IELB1 REMAINS
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%XNEBOR)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%YNEBOR)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ZNEBOR)
!
      CALL bief_allvec(1,mesh%XNEBOR,'XNEBOR',ielb1,2,1,mesh)
      CALL bief_allvec(1,mesh%YNEBOR,'YNEBOR',ielb1,2,1,mesh)
!
      IF(d.GE.3) THEN !JAJ NOT USED, ACTUALLY
        CALL bief_allvec(1,mesh%ZNEBOR,'ZNEBOR',ielb1,2,1,mesh)
      ELSE
        CALL bief_allvec(1,mesh%ZNEBOR,'ZNEBOR',    0,2,0,mesh)
      ENDIF
!
!     COORDINATES BY POINTS: X, Y AND Z
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%X)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%Y)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%Z)
      CALL bief_allvec(1,mesh%X,'X     ',ielm1,1,1,mesh)
      CALL bief_allvec(1,mesh%Y,'Y     ',ielm1,1,1,mesh)
      IF(d.GE.3) THEN
        CALL bief_allvec(1,mesh%Z,'Z     ',ielm1,1,1,mesh)
      ELSE
        CALL bief_allvec(1,mesh%Z,'Z     ',    0,1,0,mesh)
      ENDIF
!
!     COS AND SIN OF THE LATITUDE
!     WITH IELM (EXAMPLE : VELOCITY IN CORIOLIS)
!     COSLAT AND SINLAT ARE WORKING ARRAYS, TO WHICH
!     THE STRUCTURE OF X IS GIVEN TO BEGIN WITH.
!     THEY CAN BE EXTENDED TO ELEMENT IELM AT A LATER DATE.
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%COSLAT)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%SINLAT)
      IF(spheri) THEN
        CALL bief_allvec(1,mesh%COSLAT,'COSLAT',ielm,1,2,mesh)
        CALL bief_allvec(1,mesh%SINLAT,'SINLAT',ielm,1,2,mesh)
        CALL cpstvc(mesh%X,mesh%COSLAT)
        CALL cpstvc(mesh%X,mesh%SINLAT)
      ELSE
        CALL bief_allvec(1,mesh%COSLAT,'COSLAT',   0,1,0,mesh)
        CALL bief_allvec(1,mesh%SINLAT,'SINLAT',   0,1,0,mesh)
      ENDIF
!
!     DISTANCES TO BOUNDARIES : DISBOR
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%DISBOR)
      CALL bief_allvec(1,mesh%DISBOR,'DISBOR',ielb0,1,1,mesh)
!
!     WORKING MATRIX (INTERNAL TO BIEF), WITH CLASSICAL STORAGE
!
      stocfg = cfg(1)
      cfg(1) = 1
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%M)
      CALL bief_allmat(mesh%M,'M     ',ielm,ielm,cfg,'Q','Q',mesh)
      cfg(1) = stocfg
!
!     WORKING MATRIX BY SEGMENT
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%MSEG)
!     FROM 5.9 ON, ALWAYS DONE IN 2D (IELM=11,12,13 OR 14)
      IF(cfg(1).EQ.3.OR.10*(ielm/10).EQ.10) THEN
        CALL bief_allmat(mesh%MSEG,'MSEG  ',ielm,ielm,cfg,'Q','Q',mesh)
      ELSE
        CALL bief_allmat(mesh%MSEG,'MSEG  ',ielm,ielm,cfg,'0','0',mesh)
      ENDIF
!
!     WORKING ARRAY FOR A NOT ASSEMBLED VECTOR
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%W)
      CALL bief_allvec(1,mesh%W,'W     ',
     &                 ielm0,bief_nbpel(ielm,mesh),2,mesh)
!
!     WORKING ARRAY FOR A NORMAL VECTOR
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%T)
      CALL bief_allvec(1,mesh%T,'T     ',ielm,1,2,mesh)
!
!     VNOIN: ARRAY WITH NORMALS VNOIN FOR FINITE VOLUMES
!     CMI: COORDINATES OF THE MIDDLE OF THE SEGMENTS (KINETIC SCHEMES)
!     DTHAUT:
!     DPX,DPY: GRADIENTS OF THE BASE FUNCTIONS
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%VNOIN)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%CMI)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%AIRST)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%DTHAUT)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%DPX)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%DPY)
      IF(equa(1:15).EQ.'SAINT-VENANT VF') THEN
        nseg=bief_nbseg(ielm1,mesh)
        CALL bief_allvec(1,mesh%VNOIN ,'VNOIN ',3*nseg,1,0,mesh)
        CALL bief_allvec(1,mesh%CMI   ,'CMI   ',2*nseg,1,0,mesh)
        CALL bief_allvec(1,mesh%AIRST ,'AIRST ',2*nseg,1,0,mesh)
        CALL bief_allvec(1,mesh%DTHAUT,'DTHAUT',ielm1 ,1,2,mesh)
        CALL bief_allvec(1,mesh%DPX   ,'DPX   ',ielm0 ,3,2,mesh)
        CALL bief_allvec(1,mesh%DPY   ,'DPY   ',ielm0 ,3,2,mesh)
      ELSE
        CALL bief_allvec(1,mesh%VNOIN ,'VNOIN ',     0,1,0,mesh)
        CALL bief_allvec(1,mesh%CMI   ,'CMI   ',     0,1,0,mesh)
        CALL bief_allvec(1,mesh%AIRST ,'AIRST ',     0,1,0,mesh)
        CALL bief_allvec(1,mesh%DTHAUT,'DTHAUT',     0,1,0,mesh)
        CALL bief_allvec(1,mesh%DPX   ,'DPX   ',     0,1,0,mesh)
        CALL bief_allvec(1,mesh%DPY   ,'DPY   ',     0,1,0,mesh)
      ENDIF
!     FOR COORDINATES OF CENTER OF GRAVITY OF ELEMENTS NEIGHBORING EDGES
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%COORDG)
      nseg=bief_nbseg(ielm1,mesh)
      CALL bief_allvec(1,mesh%COORDG ,'COORDG',4*nseg,1,0,mesh)
!
!     FOR PARALLEL MODE
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%XSEG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%YSEG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%IFAC)
!     THEIR ALLVEC IS IN PARINI
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%BUF_SEND)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%BUF_RECV)
!
      IF (modass .EQ.3) THEN
        IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%BUF_SEND_ERR)
        IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%BUF_RECV_ERR)
      ENDIF
!
      IF(ncsize.GT.1) THEN
!       XSEG
        CALL bief_allvec(1,mesh%XSEG,'XSEG  ',ielbor(ielm1,2),1,2,mesh)
!       YSEG
        CALL bief_allvec(1,mesh%YSEG,'YSEG  ',ielbor(ielm1,2),1,2,mesh)
!       IFAC
        CALL bief_allvec(2,mesh%IFAC,'IFAC  ',ielm ,1,2,mesh)
      ELSE
        CALL bief_allvec(1,mesh%XSEG  ,'XSEG  ',0,1,0,mesh)
        CALL bief_allvec(1,mesh%YSEG  ,'YSEG  ',0,1,0,mesh)
        CALL bief_allvec(2,mesh%IFAC  ,'IFAC  ',0,1,0,mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     1) INTEGER VALUES (ALLOCATE BECAUSE THEY ARE POINTERS)
!
      IF(reflevel.EQ.0) ALLOCATE(mesh%NELEM)
      mesh%NELEM  = bief_nbpts(ielm0,mesh)
      IF(reflevel.EQ.0) ALLOCATE(mesh%NELMAX)
!     WILL GIVE BACK NNELMAX FEEDED TO BIEF_ININDS ABOVE...
      mesh%NELMAX = bief_nbmpts(ielm0,mesh)
!
!
! I DO USE MESH%NPTFR FOR THE NUMBER OF LATERAL BOUNDARY NODES
! IELBOR(IELM0,1) CHANGED TO IELBOR(IELM1,1)
! THE PROBLEM IS, THAT FOR 3D (IELM=41):
! BIEF_NBPTS(IELBOR(IELM0,1),MESH) IS THE NUMBER OF HORIZONTAL BOUNDARY ELEMENTS
! BIEF_NBPTS(IELBOR(IELM0,2),MESH) IS THE NUMBER OF VERTICAL BOUNDARY ELEMENTS
! BIEF_NBPTS(IELBOR(IELM1,1),MESH) IS THE NUMBER OF HORIZONTAL BOUNDARY NODES
! BIEF_NBPTS(IELBOR(IELM1,2),MESH) IS THE NUMBER OF VERTICAL BOUNDARY NODES
!
! FUNNY, BUT THE 2D CASE IS NOT AFFECTED, BECAUSE THE NUMBER OF BOUNDARY
! SEGMENTS IS EQUAL TO THE NUMBER OF BOUNDARY NODES.
!
      IF(reflevel.EQ.0) ALLOCATE(mesh%NPTFR)
      mesh%NPTFR  = bief_nbpts(ielbor(ielm1,2),mesh)
      IF(reflevel.EQ.0) ALLOCATE(mesh%NPTFRX)
      mesh%NPTFRX = bief_nbmpts(ielbor(ielm1,2),mesh)
!
! NUMBER OF LATERAL BOUNDARY ELEMENTS
!
      IF(reflevel.EQ.0) ALLOCATE(mesh%NELEB)
      IF(reflevel.EQ.0) ALLOCATE(mesh%NELEBX)
!
!     3D MESH
!
      IF(ielm.EQ.31) THEN
!       HERE NNELEB HAS BEEN READ IN THE GEOMETRY FILE
        mesh%NELEB   = nneleb
        mesh%NELEBX  = nneleb
      ELSEIF(ielm.EQ.11.OR.ielm.EQ.12.OR.ielm.EQ.13.OR.ielm.EQ.14) THEN
        mesh%NELEB   = nptfr
        mesh%NELEBX  = nnptfrx
      ELSEIF(ielm.EQ.41) THEN
        mesh%NELEB   = nptfr*(nnplan-1)
        mesh%NELEBX  = nnptfrx*(nnplan-1)
      ELSEIF(ielm.EQ.51) THEN
        mesh%NELEB   = 2*nptfr*(nnplan-1)
        mesh%NELEBX  = 2*nnptfrx*(nnplan-1)
      ELSE
        WRITE(lu,*) 'ALMESH, UNEXPECTED ELEMENT FOR NELEB:',ielm
        CALL plante(1)
        stop
      ENDIF
!
      IF(reflevel.EQ.0) ALLOCATE(mesh%DIM1)
      mesh%DIM1   = dimens(ielm0)
      IF(reflevel.EQ.0) ALLOCATE(mesh%TYPELM)
      mesh%TYPELM = ielm0
      IF(reflevel.EQ.0) ALLOCATE(mesh%TYPELMBND)
      mesh%TYPELMBND = typ_bnd_elem
      IF(reflevel.EQ.0) ALLOCATE(mesh%NPOIN)
      mesh%NPOIN  = bief_nbpts(ielm1,mesh)
      IF(reflevel.EQ.0) ALLOCATE(mesh%NPMAX)
      mesh%NPMAX  = bief_nbmpts(ielm1,mesh)
      IF(reflevel.EQ.0) ALLOCATE(mesh%MXPTVS)
      mesh%MXPTVS = mxptvs
      IF(reflevel.EQ.0) ALLOCATE(mesh%MXELVS)
      mesh%MXELVS = mxelvs
!     LV WILL BE RECOMPUTED LATER
      IF(reflevel.EQ.0) ALLOCATE(mesh%LV)
      mesh%LV     = 1
      IF(reflevel.EQ.0) ALLOCATE(mesh%NSEG)
      mesh%NSEG = bief_nbseg(ielm1,mesh)
!
!     2) ARRAYS OF INTEGERS
!
!     ALLOCATES IKLE AND KLEI (SAME SIZE, 2 INVERTED DIMENSIONS)
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%IKLE)
      CALL bief_allvec(2,mesh%IKLE,'IKLE  ',
     &                 ielm0,bief_nbpel(ielm,mesh),1,mesh)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%KLEI)
      CALL bief_allvec(2,mesh%KLEI,'KLEI  ',
     &                 ielm0,bief_nbpel(ielm,mesh),1,mesh)
!
!     IFABOR
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%IFABOR)
      CALL bief_allvec(2,mesh%IFABOR,'IFABOR',
     &                 ielm0,bief_nbfel(ielm,mesh),1,mesh)
!
!     NELBOR
!
! NELBOR & NULONE
! IT IS NOW CHANGED TO VERTICAL BOUNDARY ELEMENT...
! 2D NOT AFFECTED, 3D USAGE NELBO3(NPTFR,NETAGE) - NO. OF LAT. BD. ELEMENTS
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NELBOR)
      CALL bief_allvec(2,mesh%NELBOR,'NELBOR',ielbor(ielm0,2),1,1,mesh)
!
!     NULONE
!
!
! EXTRAORDINARILY STRANGE GEOMETRICALLY
! IN 2D NUMBER OF BOUNDARY NODES IS EQUAL TO THE NUMBER OF BOUNDARY
! ELEMENTS... IN 3D IT IS NOT THE CASE!
! NULONE 3D IS USED INTERNALLY AS: NULONE(NPTFR,NETAGE,4)
! "ASSOCIATES THE LOCAL BOUNDARY NUMBERING TO LOCAL 3D NUMBERING"
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NULONE)
!
      CALL bief_allvec(2,mesh%NULONE,'NULONE',
     &                 ielbor(ielm0,2),
     &                 bief_nbpel(ielbor(ielm,2),mesh),1,mesh)
!
!!! NOTE : FOR THE TETRAHEDRONS, THIS IS NO LONGER THE CASE. WE READ THE
!          BOUNDARY CONNECTIVITY TABLE BEFORE INITIALISING THE NUMBER OF
!          BOUNDARY NODES AND ELEMENTS. THIS IS WELL DEFINED NOW.
!
! IN 2D IT IS CALL BIEF_ALLVEC(2, MESH%NULONE, 'NULONE', 0,  2, 1,MESH)
! WHICH, IN 2D ONLY IS EQUIVALENT TO
!             CALL BIEF_ALLVEC(2, MESH%NULONE, 'NULONE', 1,  2, 1,MESH)
! IN 3D IT IS CALL BIEF_ALLVEC(2, MESH%NULONE, 'NULONE', 20, 4, 1,MESH)
!
!
!     KP1BOR
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%KP1BOR)
      CALL bief_allvec(2,mesh%KP1BOR,'KP1BOR', ielbor(ielm,1),2,1,mesh)
!
!     NBOR: GLOBAL NUMBERS OF THE BOUNDARY NODES
!     ALLOCATES NBOR, IT WILL BE READ IN LECLIM
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NBOR)
      CALL bief_allvec(2,mesh%NBOR,'NBOR  ',ielbor(ielm,2),1,1,mesh)
!
!     IKLBOR: IKLE FOR THE SEGMENTS OR BOUNDARY SIDES
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%IKLBOR)
!
      IF(ielm.EQ.41.OR.ielm.EQ.51) THEN
!       LATERAL BOUNDARY SIDES
!       SEE ININDS FOR ELEMENTS 20 AND 21
!       CALL BIEF_ALLVEC(2,MESH%IKLBOR,'IKLBOR',
!    &                   IELBOR(IELM0,2),
!    &                   BIEF_NBPEL(IELBOR(IELM1,2),MESH),1,MESH)
        CALL bief_allvec(2,mesh%IKLBOR,'IKLBOR',
     &                   mesh%NELEBX,
     &                   bief_nbpel(ielbor(ielm1,2),mesh),0,mesh)
      ELSEIF(ielm.EQ.11.OR.ielm.EQ.12.OR.ielm.EQ.13
     &   .OR.ielm.EQ.31) THEN
        CALL bief_allvec(2,mesh%IKLBOR,'IKLBOR',
     &                   ielbor(ielm0,1),
     &                   bief_nbpel(ielbor(ielm ,2),mesh),1,mesh)
      ELSE
        WRITE(lu,*) 'ALMESH : UNKNOWN ELEMENT FOR IKLBOR:',ielm
        CALL plante(1)
        stop
      ENDIF
!
!     IFANUM: NUMBER OF THE SIDE IN ADJACENT ELEMENT
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%IFANUM)
      IF(cfg(1).EQ.2) THEN
        CALL bief_allvec(2,mesh%IFANUM,'IFANUM',
     &                   3*bief_nbmpts(ielm0,mesh)
     &                   + bief_nbmpts(01   ,mesh) ,1,0 ,mesh)
      ELSEIF(cfg(1).NE.1.AND.cfg(1).NE.3) THEN
        WRITE(lu,99) cfg(1)
99      FORMAT(1x,'ALMESH : UNKNOWN STORAGE:',1i6)
        CALL plante(1)
        stop
      ELSE
        CALL bief_allvec(2,mesh%IFANUM,'IFANUM', 0 , 1,0,mesh )
      ENDIF
!
!     IKLEM1: INVERSE CONNECTIVITY TABLE FOR FRONTAL PRODUCT
!     LIMVOI: LIMITING NUMBER OF A GIVEN NUMBER OF NEIGHBOURS
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%IKLEM1)
      CALL first_all_biefobj(mesh%LIMVOI)
      IF(cfg(2).EQ.2) THEN
!       CALL BIEF_ALLVEC(2,MESH%IKLEM1,'IKLEM1',
!                        NBMPTS(IELM1)*MXPTVS,4,0,MESH)
!       FOR OPTASS=3: SYM AND NOT SYM ARE DIFFERENT
        CALL bief_allvec(2,mesh%IKLEM1,'IKLEM1',
     &                   bief_nbmpts(ielm1,mesh)*mxptvs,8,0,mesh)
!       11 IS HERE THE MAXIMUM MXPTVS PROGRAMMED IN OPASS
        CALL bief_allvec(2,mesh%LIMVOI,'LIMVOI',11,2,0,mesh)
      ELSEIF(cfg(2).NE.1) THEN
        WRITE(lu,97) cfg(2)
97      FORMAT(1x,'ALMESH: UNKNOWN MATRIX-VECTOR PRODUCT:',1i6)
        CALL plante(1)
        stop
      ELSE
        CALL bief_allvec(2,mesh%IKLEM1,'IKLEM1',0,4,0,mesh)
        CALL bief_allvec(2,mesh%LIMVOI,'LIMVOI',0,2,0,mesh)
      ENDIF
!
!     INTEGER ARRAYS FOR SEGMENT-BASED STORAGE
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%GLOSEG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ELTSEG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ORISEG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%GLOSEGBOR)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ELTSEGBOR)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ORISEGBOR)
!
      CALL bief_allvec(2,mesh%GLOSEG,'GLOSEG',
     &                 bief_nbseg(ielm,mesh),2,0,mesh)
      CALL bief_allvec(2,mesh%ELTSEG,'ELTSEG',
     &                 bief_nbmpts(ielm0,mesh),
     &                 bief_nbsegel(ielm,mesh),0,mesh)
      CALL bief_allvec(2,mesh%ORISEG,'ORISEG',
     &                 bief_nbmpts(ielm0,mesh),
     &                 bief_nbsegel(ielm,mesh),0,mesh)
      CALL bief_allvec(2,mesh%GLOSEGBOR,'GSGBOR',
     &                 bief_nbseg(ielb1,mesh),2,0,mesh)
      CALL bief_allvec(2,mesh%ELTSEGBOR,'ESGBOR',
     &                 bief_nbmpts(ielb0,mesh),
     &                 bief_nbsegel(ielb1,mesh),0,mesh)
      CALL bief_allvec(2,mesh%ORISEGBOR,'OSGBOR',
     &                 bief_nbmpts(ielb0,mesh),
     &                 bief_nbsegel(ielb1,mesh),0,mesh)
!
!     INTEGER ARRAY FOR THE METHOD OF CHARACTERISTICS
!     THE STARTING ELEMENT (0 IF NOT IN THIS SUBDOMAIN)
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ELTCAR)
      CALL bief_allvec(2,mesh%ELTCAR,'ELTCAR',ielm,1,1,mesh)
!
!     INTEGER ARRAYS FOR PARALLEL MODE
!
!     KNOLG
!     NACHB
!     ISEG
!     INDPU
!     NHP
!     NHM
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%KNOLG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NACHB)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%ISEG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%INDPU)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NHP)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NHM)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%IFAPAR)
      IF (reflevel.EQ.0) ALLOCATE(mesh%NB_NEIGHB)
      IF (reflevel.EQ.0) ALLOCATE(mesh%NB_NEIGHB_SEG)
!     THERE ALLVEC IS IN PARINI
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NB_NEIGHB_PT)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%LIST_SEND)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NH_COM)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NB_NEIGHB_PT_SEG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%LIST_SEND_SEG)
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NH_COM_SEG)
!
!     DATA STRUCTURE IN PARALLEL
!
      IF(ncsize.GT.1) THEN
!
        CALL bief_allvec(2,mesh%KNOLG,'KNOLG ',
     &                   ielm1            ,1,1,mesh)
        CALL bief_allvec(2,mesh%NACHB,'NACHB ',
     &                   nbmaxnshare*nptir,1,0,mesh)
        CALL bief_allvec(2,mesh%ISEG ,'ISEG  ',
     &                   ielbor(ielm1,1)  ,1,1,mesh)
        CALL bief_allvec(2,mesh%INDPU,'INDPU ',
     &                   ielm1              ,1,1,mesh)
        CALL bief_allvec(2,mesh%NHP  ,'NHP   ',
     &                   nbmaxdshare*2*nptir,1,0,mesh)
        CALL bief_allvec(2,mesh%NHM  ,'NHM   ',
     &                   nbmaxdshare*2*nptir,1,0,mesh)
        CALL bief_allvec(2,mesh%IFAPAR,'IFAPAR',10,6,1,mesh)
!
        DO i=1,6*bief_nbpts(10,mesh)
          mesh%IFAPAR%I(i)=0
        ENDDO
!
      ELSE
!
        CALL bief_allvec(2,mesh%KNOLG ,'KNOLG ',0,1,0,mesh)
        CALL bief_allvec(2,mesh%NACHB ,'NACHB ',0,1,0,mesh)
        CALL bief_allvec(2,mesh%ISEG  ,'ISEG  ',0,1,0,mesh)
        CALL bief_allvec(2,mesh%INDPU ,'INDPU ',0,1,0,mesh)
        CALL bief_allvec(2,mesh%NHP   ,'NHP   ',0,1,0,mesh)
        CALL bief_allvec(2,mesh%NHM   ,'NHM   ',0,1,0,mesh)
        CALL bief_allvec(2,mesh%IFAPAR,'IFAPAR',0,1,0,mesh)
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!  FINITE VOLUMES
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%NUBO)
      IF(equa(1:15).EQ.'SAINT-VENANT VF') THEN
        CALL bief_allvec(2,mesh%NUBO,'NUBO  ',2*nseg,1,0,mesh)
      ELSE
        CALL bief_allvec(2,mesh%NUBO,'NUBO  ',     0,1,0,mesh)
      ENDIF
!
      IF (reflevel.EQ.0) CALL first_all_biefobj(mesh%JMI)
      IF(equa(1:15).EQ.'SAINT-VENANT VF') THEN
        CALL bief_allvec(2,mesh%JMI,'JMI   ',nseg,1,0,mesh)
      ELSE
        CALL bief_allvec(2,mesh%JMI,'JMI   ',   0,1,0,mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     PART SKIPPED WHEN STARTING REFINEMENT
!
      IF(reflevel.EQ.0) THEN
!
!-----------------------------------------------------------------------
!
!     FILLS ARRAYS IKLE, X AND Y (AND Z)
!
!     READS TO KNOLG ARRAY LOCAL TO GLOBAL NUMBERING IN PARALLEL
      IF(ncsize.GT.1) THEN
        CALL get_mesh_l2g_numbering(fformat,nfic,mesh%KNOLG%I,
     &                              npoin,ierr)
        CALL check_call(ierr,'ALMESH:GET_MESH_L2G_NUMBERING')
      ENDIF
!     PRISMS CUT INTO TETRAHEDRONS
!
      IF(ielm.EQ.51) THEN
!
!       NOTE : NO Z HERE, AS IELM.EQ.41, SEE NOTE BELOW
        CALL read_mesh_coord(fformat,nfic,mesh%X%R,mesh%Y%R,npoin,
     &                projec,lati,longi)
!
        CALL cpikle3(mesh%IKLE%I,ikles,nelem,nnelmax,npoin,nnplan,
     &               mesh%KNOLG%I)
!
!     PRISMS
!
      ELSEIF(ielm.EQ.41) THEN
!
!       NOTE : WITH PRISMS Z IS COMPUTED WITH ZF AND H, OR
!              READ IN THE PREVIOUS COMPUTATION FILE, HENCE NO Z HERE
        CALL read_mesh_coord(fformat,nfic,mesh%X%R,mesh%Y%R,npoin,
     &                       projec,lati,longi)
!
        CALL cpikle2(mesh%IKLE%I,mesh%KLEI%I,ikles,
     &               nelem,nnelmax,npoin,nnplan)
!
!     TRIANGLES OR TETRAHEDRONS
!
      ELSEIF(ielm.EQ.11.OR.ielm.EQ.12.OR.ielm.EQ.13.OR.ielm.EQ.14
     &   .OR.ielm.EQ.31) THEN
!
!       IKLES(NDP,NELEM) COPIED INTO IKLE(NELMAX,NDP) AND KLEI(NDP,NELMAX)
        DO i = 1,ndp
          DO ielem  = 1,nelem
            mesh%IKLE%I((i-1)*nnelmax+ielem) = ikles((ielem-1)*ndp+i)
            mesh%KLEI%I((ielem-1)*ndp+i)     = ikles((ielem-1)*ndp+i)
          ENDDO
        ENDDO
        IF(ielm.EQ.11.OR.ielm.EQ.12.OR.ielm.EQ.13.OR.ielm.EQ.14) THEN
          CALL read_mesh_coord(fformat,nfic,mesh%X%R,mesh%Y%R,npoin,
     &                         projec,lati,longi)
        ELSEIF(ielm.EQ.31) THEN
!         TETRAHEDRONS: READS THE Z COORDINATE AFTER X AND Y
          CALL read_mesh_coord(fformat,nfic,mesh%X%R,mesh%Y%R,npoin,
     &                         projec,lati,longi,z=mesh%Z%R)
        ENDIF
!
      ELSE
!
! OTHER ELEMENT TYPES
!
        WRITE(lu,*) 'ALMESH: UNKNOWN ELEMENT:',ielm
        CALL plante(1)
        stop
!
      ENDIF
!
!     NOW WE HAVE KNOLG IN 2D
!     FINDING THE NUMBER OF POINTS IN THE ORIGINAL MESH
!     = MAXIMUM OF ALL KNOLG OF ALL SUB-DOMAINS
!
      npoin_glob=0
      IF(ncsize.GT.1) THEN
        DO i=1,npoin
          npoin_glob=max(npoin_glob,mesh%KNOLG%I(i))
        ENDDO
        npoin_glob=p_max(npoin_glob)
      ENDIF
!
!     COMPLEMENTS: ARRAYS X, Y FOR PRISMS AND TETRAHEDRONS
!                  KNOLG FOR PRISMS AND TETRAHEDRONS
!
      IF(ielm.EQ.41.OR.ielm.EQ.51) THEN
        DO iplan = 2,nnplan
          CALL ov_2( 'X=Y     ' , mesh%X%R,iplan, mesh%X%R,1,
     &                            mesh%X%R,1, 0.d0, nnpmax,npoin)
          CALL ov_2( 'X=Y     ' , mesh%Y%R,iplan, mesh%Y%R,1,
     &                            mesh%Y%R,1, 0.d0, nnpmax,npoin)
        ENDDO
        IF(ncsize.GT.1) THEN
          DO iplan = 2,nnplan
            DO i=1,npoin
              mesh%KNOLG%I(i+(iplan-1)*npoin)=
     &        mesh%KNOLG%I(i)+(iplan-1)*npoin_glob
            ENDDO
          ENDDO
        ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
!
!     END OF PART SKIPPED WHEN STARTING REFINEMENT
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! DEALLOCATES TEMPORARY ARRAYS
!
      DEALLOCATE(ikles)
!
!-----------------------------------------------------------------------
!
      WRITE(lu,*) 'MESH: ',nom,' ALLOCATED'
!
!-----------------------------------------------------------------------
!
      RETURN
      END