partel.F

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE partel
!                       *****************
     &     (nameinp, namecli, nparts, pmethod, fformat,
     &      namesec, namezfi, nameseu)
!
!
!***********************************************************************
! PARTEL
!***********************************************************************
!
!BRIEF    PREPROCESSING STEP BEFORE A PARALLEL COMPUTATION
!
!HISTORY   R. KOPMANN (BAW)
!+
!+
!+         FIRST  VERSION JANUARY-MARCH 2000
!
!HISTORY   JAJ
!+      12/12/2000
!+      SECOND VERSION PINXIT
!+     PARTITIONING OF GEOMETRY AND 2D RESULT FILES POSSIBLE

!HISTORY   JAJ
!+      22/02/2002
!+      THIRD VERSION
!+     ERRORS IN BC VALUES IN DECOMPOSED BC FILES REMOVED
!+     ERRONEOUS TREATMENT OF ISLANDS DEBUGGED
!
!HISTORY   J-M HERVOUET ; JAJ
!+      17/04/2002
!+     FOURTH VERSION
!+     PARTITIONING FOR 3D RESULT FILES DONE BY JMH
!+     INCLUDING BOTH PARTITIONING METHODS AND BEAUTIFYING BY JAJ
!
!HISTORY  J-M HERVOUET
!+     21/01/2003
!+     FIFTH VERSION
!+     CORRECTED A WRONG DIMENSION OF THE ARRAY CUT, AN ERROR
!+     OCCURING BY A LARGER NUMBER OF PROCESSORS
!
!HISTORY  J-M HERVOUET
!+       12/03/2003
!+      SEVENTH VERSION
!+      ALGORITHM CHANGED : A SEGMENT IS IN A SUBDOMAIN IF IT BELONGS
!+      TO AN ELEMENT IN THE SUBDOMAIN NOT IF THE 2 POINTS OF THE
!+      SEGMENT BELONG TO THE SUBDOMAIN.
!+       SPECIFIC ELEBD INCLUDED, ALL REFERENCE TO MPI OR BIEF REMOVED
!
!HISTORY  J-M HERVOUET
!+        01/09/2003
!+      EIGHTH VERSION
!+      UBOR AND VBOR INVERTED LINE 613 WHEN READING THE CLI FILE.
!+    OTHER MODIFICATIONS PERFORMED
!
!HISTORY  C. DENIS J-M HERVOUET (SINETICS & LNHE)
!+        22/06/2012
!+        V6P2
!+      DOUBLE PRECISION SERAFIN NOW POSSIBLE.
!
!HISTORY  Y. AUDOUIN (STFC & LNHE)
!+        25/06/2012
!+        V6P2
!+      INTERFACE FOR LATEST RELEASE OF METIS (>= VERSION 5)
!
!HISTORY  J-M HERVOUET (EDF LAB, LNHE)
!+        27/03/2014
!+        V7P0
!+    ARGUMENTS ADDED TO THE CALL ELEBD.
!
!HISTORY Y AUDOUIN (LNHE)
!+       25/05/2015
!+       V7P0
!+       MODIFICATION TO COMPLY WITH THE HERMES MODULE
!
!HISTORY C. COULET (ARTELIA)
!+       01/09/2016
!+       V7P2
!+       MODIFICATION TO ADD THE WEIR FILE MANAGEMENT
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| NAMEGEO        |<--| NAME OF THE GEOMETRY FILE
!| NAMECLI        |<--| NAME OF THE BOUNDARY CONDITIONS FILE
!| NPARTS         |<--| NUMBER OF PARTITIONS
!| PMETHOD        |<--| 1: FOR METIS 2: FOR SCOTCH
!| FFORMAT        |<--| FORMAT OF THE GEOMETRY FILE
!| NAMESEC        |<--| NAME OF THE SECTION FILE ' ' IF THERE ARE NONE
!| NAMEZFI        |<--| NAME OF THE FRICTION ZONE FILE ' ' IF THERE ARE NONE
!| NAMESEU        |<--| NAME OF THE WEIR FILE ' ' IF THERE ARE NONE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_partel
      USE declarations_special
      USE bief, ONLY : nbmaxnshare, nptir,
     &                 read_mesh_info, ncsize
      USE interface_hermes
      USE mod_hash_table
      USE mod_numbering_open_boundaries
      USE mod_init_dataval
      USE mod_compute_boundary_and_interface
      USE mod_write_solutions
      USE mod_handle_sections
      USE mod_handle_friction_zones
      USE mod_handle_weirs
!
      IMPLICIT NONE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      CHARACTER(LEN=PATH_LEN), INTENT(IN) :: NAMEINP
      CHARACTER(LEN=PATH_LEN), INTENT(IN) :: NAMECLI
      INTEGER, INTENT(IN) :: NPARTS
      INTEGER, INTENT(IN) :: PMETHOD
      CHARACTER(LEN=8), INTENT(INOUT) :: FFORMAT
      CHARACTER(LEN=PATH_LEN), INTENT(IN) :: NAMESEC
      CHARACTER(LEN=PATH_LEN), INTENT(IN) :: NAMEZFI
      CHARACTER(LEN=PATH_LEN), INTENT(IN) :: NAMESEU
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!
      INTEGER NVAR, NPLAN, NPTFR
      INTEGER NELEM, NPOIN, NDP, NELEM2, NPOIN2, NDP_BND
      INTEGER DIM_MESH
      !USED IN SOME LOOP TO AVOID THE MULTIPLE CALL OF HASH_TABLE_GET
      INTEGER TEMP
!
      INTEGER, ALLOCATABLE :: IKLES(:), IKLES_P(:)
      INTEGER, ALLOCATABLE :: IKLES3D(:),IKLES3D_P(:,:,:)
      INTEGER, ALLOCATABLE :: IRAND(:)
      INTEGER, ALLOCATABLE :: LIHBOR(:),LIUBOR(:)
      INTEGER, ALLOCATABLE :: LIVBOR(:),LITBOR(:),COLOR(:)
      DOUBLE PRECISION, ALLOCATABLE :: HBOR(:),UBOR(:),VBOR(:)
      DOUBLE PRECISION, ALLOCATABLE :: CHBORD(:)
      DOUBLE PRECISION, ALLOCATABLE :: TBOR(:),ATBOR(:),BTBOR(:)
      INTEGER, ALLOCATABLE :: NBOR(:),IKLE_BND(:)
      INTEGER, ALLOCATABLE :: NPOIN_P(:), NELEM_P(:), NPTFR_P(:)
      INTEGER, ALLOCATABLE :: NPTIR_P(:)
      INTEGER, ALLOCATABLE :: NUMLIQ(:)
      INTEGER, ALLOCATABLE :: KNOLG(:,:)
      INTEGER, ALLOCATABLE :: ELELG(:,:)
      INTEGER, ALLOCATABLE :: CUT(:), SORT(:)
      TYPE(hash_table) :: CUT_P, KNOGL
      INTEGER, ALLOCATABLE :: PART_P(:,:)
!
!     FOR DOUBLE PRECISION SERAFIN FORMAT
!
      DOUBLE PRECISION, ALLOCATABLE    :: F(:,:)
!
      DOUBLE PRECISION TIMES
!
      INTEGER :: NINP=10
      INTEGER :: NOUT=17, nclm=18
      INTEGER TIME(3),DATE(3), DATE_TMP(6)
!
      CHARACTER(LEN=80)  :: TITLE
      CHARACTER(LEN=32), ALLOCATABLE  :: VARIABLE(:)
      CHARACTER(LEN=PATH_LEN) :: NAMECLM
      CHARACTER(LEN=12)  :: FMT4
!
      INTEGER MAX_NELEM_P
      INTEGER  MAX_NPOIN_P,MAX_N_NEIGH
      INTEGER I, J, K, L, M, ISO, IDUM
      INTEGER ISEG, NTIMESTEP
!
      REAL XSEG, YSEG
      LOGICAL TIMECOUNT
!
! METISOLOGY
!
      INTEGER, ALLOCATABLE :: EPART(:), NPART(:)
!
! FOR CALLING FRONT2
!
      INTEGER, ALLOCATABLE :: KP1BOR(:,:)
!
! FOR CALLING BIEF MESH SUBROUTINES (TO BE OPTIMISED SOON):
!
      INTEGER, ALLOCATABLE :: IFABOR(:,:), NELBOR(:)
      INTEGER, ALLOCATABLE :: IKLE(:,:)
!
! TIME MEASURING
!
      INTEGER  TDEB, TFIN, TDEBP, TFINP, TEMPS, PARSEC
!     HERMES TEMPORARY ARRAYS
      DOUBLE PRECISION,ALLOCATABLE :: DATAVAL(:,:,:)
      INTEGER :: IERR
      CHARACTER(LEN=16), ALLOCATABLE :: VAR_NAME(:), VAR_UNIT(:)
      INTEGER NELEBD
      INTEGER TYP_ELEM, TYP_BND_ELEM
!
!----------------------------------------------------------------------
!
!JAJ NEW FOR PARALLEL CHARACTERISTICS ////
! HALO ELEMENTS: THESE ADJACENT TO THE INTERFACE EDGES HAVING
! NEIGHBOURS BEHIND A BOUNDARY
!
      ! THE ELEMENTAL GLOBAL->LOCAL NUMBERING TRANSLATION TABLE
      ! THIS IS ELEGL SAVED FROM ALL PARTITIONS FOR FURTHER USE
      !INTEGER, ALLOCATABLE :: GELEGL(:,:)
      TYPE(hash_table) :: GELEGL
!
      ! THE HALO ELEMENTS NEIGHBOURHOOD DESCRIPTION FOR A HALO CELL
      INTEGER, ALLOCATABLE :: IFAPAR(:,:)
!
      ! THE NUMBER OF HALO CELLS PRO PARTITION
      INTEGER, ALLOCATABLE :: NHALO(:)
!
!     WORK VARIABLES
!
      INTEGER IFALOC(3)
      INTEGER NDP_2D
      INTEGER EF
      INTEGER, ALLOCATABLE :: NBRE_EF(:),EF_I(:,:)
      INTEGER, ALLOCATABLE :: TAB_TMP(:),EF_II(:,:)
      LOGICAL HALO
      INTEGER NOEUD
      INTEGER, ALLOCATABLE :: NBRE_EF_I(:)
!
! #### FOR SECTIONS
!
      LOGICAL :: WITH_SECTIONS
!
! #### FOR ZONES
!
      LOGICAL :: WITH_ZONES
!
!     CD: FLAG FOR SERAFIN FORMAT (.TRUE. IF DOUBLE PRECISION)
!
!     PARTEL EXTENS
      CHARACTER(LEN=11) :: EXTENS
      EXTERNAL extens
!
! #### FOR WEIRS
!
      LOGICAL :: WITH_WEIRS
!
!     FOR PARTEL CONCAT
      INTEGER :: NCLM_IDX, OFFSET_BEGIN, OFFSET_END, WRITTEN_LINES
      CHARACTER(LEN=PATH_LEN) :: NAMECLM_IDX
!
!----------------------------------------------------------------------
!
      ndp_2d=3
      ! SET NCSIZE TO 1 TO USE VOISIN AND READ_MESH_INFO IN SERIAL MODE
      ncsize=1
!
      WRITE(lu,*) '+---- PARTEL: BEGINNING -------------+'
      CALL system_clock (count=temps, count_rate=parsec)
      timecount = .true.
      IF (parsec==0) timecount = .false.  ! COUNT_RATE == 0 : NO CLOCK
      IF (timecount) tdeb = temps
!
      WRITE(lu,*)'FICHIER:',nameinp
      CALL open_mesh(fformat, nameinp, ninp, 'READ     ', ierr)
      CALL check_call(ierr, 'PARTEL:OPENMESH:INP')
!
      CALL open_bnd(fformat,namecli,ninp,'READ     ',ierr)
      CALL check_call(ierr,'PARTEL:OPEN_BND:NCLI')
!
!----------------------------------------------------------------------
!
! START READING THE GEOMETRY OR RESULT FILE
!
      CALL read_mesh_info(fformat,ninp,title,nvar,npoin,typ_elem,nelem,
     &                    nptfr,nptir,ndp,nplan,x_orig,y_orig,
     &                    typ_bnd_elem,nelebd)
!
      ALLOCATE(nbre_ef_i(nparts), stat=ierr)
      CALL check_allocate(ierr, 'NBRE_EF_I')

      ! GET THE VARIABLES NAMES
      ALLOCATE(var_name(nvar),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:VAR_NAME')
      ALLOCATE(var_unit(nvar),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:VAR_UNIT')
      ALLOCATE(variable(nvar),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:VAR_UNIT')

      CALL get_data_var_list(fformat,ninp,nvar,var_name,var_unit,ierr)
      CALL check_call(ierr,'PARTEL:GET_DATA_VAR_LIST:INP')
      DO i=1,nvar
        variable(i)(1:16) = var_name(i)
        variable(i)(17:32) = var_unit(i)
      ENDDO
      DEALLOCATE(var_name)
      DEALLOCATE(var_unit)
!
! READ THE REST OF THE SELAFIN FILE
! 10 INTEGERS, THE FIRST IS THE NUMBER OF RECORDS (TIMESTEPS)
!
      CALL get_mesh_date(fformat,ninp,date_tmp,ierr)
      CALL check_call(ierr,'PARTEL:GET_MESH_DATE:INP')
      DO i=1,3
        date(i) = date_tmp(i)
        time(i) = date_tmp(i+3)
      ENDDO
!
      IF(nplan.GT.1) THEN
        WRITE(lu,*) ' '
        WRITE(lu,*) '3D MESH DETECTED'
        npoin2 = npoin/nplan
        nelem2 = nelem/(nplan-1)
        WRITE(lu,*) 'NDP NODES PER ELEMENT:             ',ndp
        WRITE(lu,*) 'NPLAN NUMBER OF MESH LEVELS:       ',nplan
        WRITE(lu,*) 'NPOIN2 NUMBER OF 2D MESH NODES:    ',npoin2
        WRITE(lu,*) 'NPOIN NUMBER OF 3D MESH NODES:     ',npoin
        WRITE(lu,*) 'NELEM2 NUMBER OF 2D MESH ELEMENTS: ',nelem2
        WRITE(lu,*) 'NELEM NUMBER OF 3D MESH ELEMENTS:  ',nelem
        IF (mod(npoin,nplan).NE.0) THEN
          WRITE (lu,*) 'BUT NPOIN2 /= NPOIN3/NPLAN!'
          CALL plante(1)
          stop
        ENDIF
        IF (mod(nelem,(nplan-1)).NE.0) THEN
          WRITE (lu,*) 'BUT NELEM2 /= NELEM3/NPLAN!'
          CALL plante(1)
          stop
        ENDIF
        WRITE(lu,*) ' '
        WRITE(lu,*) 'THE INPUT FILE ASSUMED TO BE 3D'
        dim_mesh = 3
      ELSE
        WRITE(lu,*) ' '
        WRITE(lu,*) 'ONE-LEVEL MESH.'
        WRITE(lu,*) 'NDP NODES PER ELEMENT:         ',ndp
        WRITE(lu,*) 'ELEMENT TYPE :                 ',typ_elem
        WRITE(lu,*) 'NPOIN NUMBER OF MESH NODES:    ',npoin
        WRITE(lu,*) 'NELEM NUMBER OF MESH ELEMENTS: ',nelem
        WRITE(lu,*) ' '
        npoin2 = npoin
        nelem2 = nelem
        WRITE(lu,*) 'THE INPUT FILE ASSUMED TO BE 2D'
        dim_mesh = 2
      ENDIF
!
! NOW LET US ALLOCATE
!
      ALLOCATE (ikles(nelem2*3),stat=ierr)
      CALL check_allocate(ierr, 'IKLES')
      IF(nplan.GT.1) THEN
        ALLOCATE (ikles3d(nelem*ndp),stat=ierr)
        CALL check_allocate(ierr, 'IKLES3D')
      ENDIF
      ALLOCATE (irand(npoin),stat=ierr)
      CALL check_allocate(ierr, 'IRAND')
!
!     SIZE 3: FIRST TWO FUNCTIONS ARE X AND Y, 3 IS ALL OTHER
!             VARIABLES (THEY WILL BE COPIED AND WRITTEN
!             ONE AFTER THE OTHER...)
!     NPOIN IS 3D HERE IN 3D
!
      ALLOCATE (f(npoin,2),stat=ierr)
      CALL check_allocate(ierr, 'F')
!
! CONNECTIVITY TABLE:
!
      IF(nplan.LE.1) THEN
        CALL get_mesh_connectivity(fformat,ninp,typ_elem,ikles,
     &                             nelem,ndp,ierr)
        CALL check_call(ierr,'PARTEL:GET_MESH_CONNECTIVITY:2D')
      ELSE
        CALL get_mesh_connectivity(fformat,ninp,typ_elem,ikles3d,
     &                             nelem,ndp,ierr)
        CALL check_call(ierr,'PARTEL:GET_MESH_CONNECTIVITY:3D')
!       BUILDING IKLES
        DO j=1,3
          DO k=1,nelem2
            ikles((k-1)*3+j)=ikles3d((k-1)*6+j)
          ENDDO
        ENDDO
      ENDIF
!
! BOUNDARY NODES INDICATIONS
!
      CALL get_bnd_ipobo(fformat,ninp,npoin,nptfr,typ_bnd_elem,
     &                   irand,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_IPOBO:NINP')
!
! IRAND IS NOT ALWAYS CORRECT AND MAY LEAD TO ERRORS
! THE BO0UNDARY FILE IS USED INSTEAD
!
! X-, Y-COORDINATES
!
      CALL get_mesh_coord(fformat,ninp,1,2,npoin,f(:,1),ierr)
      CALL check_call(ierr,'PARTEL:GET_MESH_COORD:X')
      CALL get_mesh_coord(fformat,ninp,2,2,npoin,f(:,2),ierr)
      CALL check_call(ierr,'PARTEL:GET_MESH_COORD:Y')
!
      CALL get_data_ntimestep(fformat,ninp,ntimestep,ierr)
      CALL check_call(ierr,'PARTEL:GET_DATA_NTIMESTEP')
!
      WRITE(lu,*) 'THERE ARE ',ntimestep,' TIME-DEPENDENT RECORDINGS'
!
!----------------------------------------------------------------------
!
!     READ THE BOUNDARY CONDITIONS FILE
!
      !
      CALL get_nodes_per_element(typ_bnd_elem,ndp_bnd)
      ! GET THE NUMBER OF BOUNDARY POINTS AND ELEMENTS
      CALL get_bnd_nelem(fformat,ninp,typ_bnd_elem,nelebd,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_NELEBD:NCLI')
      CALL get_bnd_npoin(fformat,ninp,typ_bnd_elem,nptfr,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_NPOIN:NCLI')
      !
      ALLOCATE(ikle_bnd(nelebd*ndp_bnd),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:IKLE_BND')
      ALLOCATE(nbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:NBOR')
      ! ALLOCATING ARRAY FOR THE BOUNDARY CONDITIONS
      ALLOCATE(lihbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:LIHBOR')
      ALLOCATE(liubor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:LIUBOR')
      ALLOCATE(livbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:LIVBOR')
      ALLOCATE(hbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:HBOR')
      ALLOCATE(ubor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:UBOR')
      ALLOCATE(vbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:VBOR')
      ALLOCATE(chbord(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:CHBORD')
      ALLOCATE(litbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:LITBOR')
      ALLOCATE(tbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:TBOR')
      ALLOCATE(atbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:ATBOR')
      ALLOCATE(btbor(nptfr),stat=ierr)
      CALL check_allocate(ierr,'PARTEL:BTBOR')
      ALLOCATE (color(nptfr),stat=ierr)
      CALL check_allocate(ierr, 'COLOR')
      ! Get the connectivity table for the boundary elements
      CALL get_bnd_connectivity(fformat,ninp,typ_bnd_elem,nelebd,
     &                    ndp_bnd,ikle_bnd,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_CONNECTIVITY:NCLI')
      ! FILL NBOR
      CALL get_bnd_numbering(fformat,ninp,typ_bnd_elem,nptfr,nbor,ierr)
      ! GET THE VALUE OF EACH BOUNDARY
      CALL get_bnd_value(fformat,ninp,typ_bnd_elem,nelebd,
     &                   lihbor,liubor,
     &                   livbor,hbor,ubor,vbor,chbord,.true.,
     &                   litbor,tbor,atbor,btbor,nptfr,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_VALUE:NCLI')
      CALL get_bnd_color(fformat,ninp,typ_bnd_elem,nelebd,
     &                   color,ierr)
      CALL check_call(ierr,'PARTEL:GET_BND_COLOR:NCLI')
!
      CALL numbering_open_boundaries(nameinp, ikle, ikles, kp1bor,
     &                   numliq, dim_mesh, npoin2, nptfr, npoin, nelem2,
     &                   nelbor, liubor, lihbor, nbor, ifabor, f,.true.)
!
!======================================================================
! PARTITIONING
!
!
!
!======================================================================
! STEP 2 : PARTITIONING THE MESH
!
! OTHER PARTITIONING METHODS SHOULD BE USED (SCOTCH FOR EXAMPLE)
!     ALL PROCESSORS PERFORM THIS TASK TO AVOID COMMUNICATION
!     THE USE OF PARMETIS OR PTSCOTCH COULD BE USED FOR LARGER MESHES
!     IF THERE WILL BE SOME MEMORY ALLOCATION PROBLEM
!======================================================================
!
      ALLOCATE(epart(nelem2),stat=ierr)
      CALL check_allocate(ierr, 'EPART')
      ALLOCATE(npart(npoin2),stat=ierr)
      CALL check_allocate(ierr, 'NPART')
!
!     PARTITIONNING METHOD
!     1 : METIS
!     2 : SCOTCH
      WRITE(lu,*) ' THE MESH PARTITIONING STEP STARTS'
      IF(timecount) THEN
!       CALL SYSTEM_CLOCK (COUNT=TEMPS, COUNT_RATE=PARSEC)
        tdebp = temps
      ENDIF
      CALL partitioner(pmethod, nelem2, npoin2, 3, nparts,
     &                 ikles, epart, npart)
      IF (timecount) THEN
!       CALL SYSTEM_CLOCK (COUNT=TEMPS, COUNT_RATE=PARSEC)
        tfinp = temps
        WRITE(lu,*) ' RUNTIME OF METIS ',
     &            (1.0*(tfinp-tdebp))/(1.0*parsec),' SECONDS'
      ENDIF
      WRITE(lu,*) ' THE MESH PARTITIONING STEP HAS FINISHED'
!
!
!----------------------------------------------------------------------
!
!======================================================================
! STEP 3 : ALLOCATE THE GLOBAL  ARRAYS NOT DEPENDING OF THE PARTITION
!
!======================================================================
!
!     KNOGL(I) =>  GLOBAL LABEL OF THE LOCAL POINT I
!
      CALL hash_table_create(knogl, 2**20)
!
!     NBRE_EF(I) => NUMBER OF FINITE ELEMENT CONTAINING I
!     I IS A GLOBAL LABEL
      ALLOCATE (nbre_ef(npoin2),stat=ierr)
      CALL check_allocate(ierr, 'NBRE_EF')
!
      ALLOCATE (part_p(npoin2,0:nbmaxnshare),stat=ierr)
      CALL check_allocate(ierr, 'PART_P')
      part_p(:,:)=0

      CALL hash_table_create(cut_p, 2**20)

      CALL hash_table_create(gelegl, nelem2)

      ALLOCATE (sort(npoin2),stat=ierr)
      CALL check_allocate(ierr, 'SORT')

      ALLOCATE (cut(npoin2),stat=ierr)
      CALL check_allocate(ierr, 'CUT')

      ALLOCATE (nelem_p(nparts),stat=ierr)
      CALL check_allocate(ierr, 'NELEM_P')

      ALLOCATE (npoin_p(nparts),stat=ierr)
      CALL check_allocate(ierr, 'NPOIN_P')

      ALLOCATE (nptfr_p(nparts),stat=ierr)
      CALL check_allocate(ierr, 'NPTFR_P')

      ALLOCATE (nptir_p(nparts),stat=ierr)
      CALL check_allocate(ierr, 'NPTIR_P')

      ALLOCATE (nhalo(nparts),stat=ierr)
      CALL check_allocate(ierr, 'NHALO')

      ALLOCATE(tab_tmp( nbmaxnshare),stat=ierr)
      CALL check_allocate(ierr, 'TAB_TMP')

      ALLOCATE(ifapar(7,nbmaxhalo),stat=ierr)
      CALL check_allocate(ierr, 'IFAPAR')
      ifapar(:,:)=0
!
!======================================================================
! STEP 4 : COMPUTE THE NUMBER OF FINITE ELEMENTS AND POINTS
!     BELONGING TO SUBMESH I
!
!======================================================================
!

!     FIRSTLY, ALL MPI PROCESSES  WORK ON THE WHOLE MESH
!     ----------------------------------------------
!
!     LOOP OVER THE FINITE ELEMENT OF THE MESH
!     TO COMPUTE THE NUMBER OF FINITE ELEMENTS CONTAINING EACH POINT NOEUD
      IF (code(1:3) == 'ART') THEN
        DO ef=1,nelem2
          DO k=1,ndp_2d
            noeud=ikles((ef-1)*3+k)
            IF (irand(noeud) .NE. 0) THEN
              epart(ef)=1
            END IF
          END DO
        END DO
      END IF

      nbre_ef(:)=0
      DO ef=1,nelem2
        DO k=1,ndp_2d
          noeud=ikles((ef-1)*3+k)
          nbre_ef(noeud)=nbre_ef(noeud)+1
        END DO
      END DO
!
!     LOOP OVER THE FINITE ELEMENT OF THE MESH TO COMPUTE
!     THE NUMBER OF THE FINITE ELEMENT AND POINTS BELONGING
!     TO SUBMESH I
!
      nelem_p=0
      npoin_p=0
      DO ef=1,nelem2
        i = epart(ef)
        IF(i >= 1 .AND. i <= nparts) THEN
          nelem_p(i)=nelem_p(i)+1
          DO k=1,ndp_2d
            noeud=ikles((ef-1)*3+k)
            IF (hash_table_get(knogl,noeud,i) .EQ. 0) THEN
              npoin_p(i)=npoin_p(i)+1
              CALL hash_table_insert(knogl,noeud,i,npoin_p(i))
            END IF
          END DO
        END IF
      END DO
!
!======================================================================
!     STEP 4 : ALLOCATION OF LOCAL ARRAYS NEEDED BY MPI PROCESSUS ID
!              WORKING ON SUBMESH ID+1
!======================================================================
!
      max_nelem_p=maxval(nelem_p)
      max_npoin_p=maxval(npoin_p)
!
!     ELEGL(E) => GLOBAL LABEL OF THE FINITE ELEMENT E
!     E IS THE LOCAL LABEL ON SUBMESH I
      ALLOCATE (elelg(max_nelem_p,nparts),stat=ierr)
      CALL check_allocate(ierr, 'ELELG')
      elelg(:,:)=0
!     KNOLG(I) => GLOBAL LABEL OF THE POINT I
!     I IS THE LOCAL LABEL ON SUBDOMAIN I
      IF(nplan.LE.0) THEN
        ALLOCATE(knolg(max_npoin_p,nparts),stat=ierr)
      ELSE
        ALLOCATE(knolg(max_npoin_p*nplan,nparts),stat=ierr)
      ENDIF
      CALL check_allocate(ierr, 'KNOLG')
      knolg(:,:)=0
!
!     EF_I(E) IS THE GLOBAL LABEL OF THE INTERFACE FINITE ELEMENT NUMBER E
      ALLOCATE (ef_i(nparts, max_nelem_p),stat=ierr)
      CALL check_allocate(ierr, 'EF_I')
!     EF_II(E) IS THE LOCAL LABEL OF THE INTERFACE FINITE ELEMENT NUMBER E
      ALLOCATE (ef_ii(nparts, max_nelem_p),stat=ierr)
      CALL check_allocate(ierr, 'EF_II')
!
!======================================================================
!     STEP 5 : INITIALISATION  OF LOCAL ARRAYS
!              (GELELG AND ELELG)
!======================================================================
!
      nelem_p=0
      DO ef=1,nelem2
        i=epart(ef)
        IF(i >= 1 .AND. i <= nparts) THEN
          nelem_p(i)=nelem_p(i)+1
          elelg(nelem_p(i),i)=ef
          CALL hash_table_insert(gelegl,ef,i,nelem_p(i))
        ENDIF
      ENDDO
!
      CALL compute_boundary_and_interface(nparts,ndp_2d,npoin_p,
     &                                    nptfr_p, elelg,
     &                                    nelem_p,ikles,knogl,cut_p,ef_i
     &                                    ,ef_ii,nptir_p,nbre_ef_i,knolg
     &                                    ,irand,part_p,nbre_ef)

!
      max_n_neigh=maxval(part_p(:,0))
      IF ( max_n_neigh > nbmaxnshare-1 ) THEN
        WRITE(lu,*) 'SERIOUS WARNING: '
        WRITE(lu,*)
     &       'AN INTERFACE NODE BELONGS TO ',
     &       'MORE THAN NBMAXNSHARE-1 SUBDOMAINS'
        WRITE(lu,*) 'TELEMAC MAY PROTEST!'
      ENDIF
      IF(max_n_neigh.GT.maxnproc) THEN
        WRITE (lu,*) 'THERE IS A NODE WHICH BELONGS TO MORE THAN ',
     &       maxnproc,' PROCESSORS, HOW COME?'
        CALL plante(1)
        stop
      ENDIF
      IF (max_n_neigh.LT.nbmaxnshare-1) max_n_neigh = nbmaxnshare-1
!
      IF(partel_concat)THEN
        IF (code.NE.'  ') THEN
          nameclm = code//'PAR'//'-CONCAT'
          nameclm_idx = code//'PAR'//'-INDEX'
        ELSE
          nameclm = namecli(1:3)//'PAR'//'-CONCAT'
          nameclm_idx = namecli(1:3)//'PAR'//'-INDEX'
        ENDIF
        CALL get_free_id(nclm)
        OPEN(nclm,file=nameclm,status='UNKNOWN',form='FORMATTED')
        rewind(nclm)
        CALL get_free_id(nclm_idx)
        OPEN(nclm_idx,file=nameclm_idx,action='WRITE')
        offset_begin=1
        offset_end=1
      ENDIF
      DO i=1,nparts
        written_lines=1
!
!       FIRST LOOP TO COMPUTE THE NUMBER OF HALO TO ALLOCATE IFAPAR
!
!       FILLING IFAPAR
!
        nhalo(i)=0
        DO j=1,nbre_ef_i(i)  ! ON PARCOURT JUSTE LES ELEMENTS FINIS INTERFACES POUR
                          ! DETERMINER DES HALO
          ef=ef_i(i, j)
          halo=.false.
          ifaloc(:)=ifabor(ef,:)
          WHERE (ifaloc .GT. 0)
            ifaloc=epart(ifaloc)
          END WHERE
          halo=any(ifaloc .GT. 0 .AND. ifaloc .NE. i)
          IF(halo) THEN
            nhalo(i)=nhalo(i)+1
            IF(nhalo(i) > nbmaxhalo) THEN
              WRITE(lu,*)  'ERROR : NBMAXHALO TOO SMALL'
              CALL plante(1)
              stop
            ENDIF
            ifapar(1,nhalo(i))=ef_ii(i, j)
            ifapar(2:4,nhalo(i))=ifaloc(:)
            ifapar(5:7,nhalo(i))=ifabor(ef_i(i, j),:)
          ENDIF
        ENDDO

!
!-----------------------------------------------------------------------
! THE CORE NAMES FOR THE OUTPUT BC FILES ACCORDING TO THE NUMBER OF PARTS
!
        ! WORKING ON THE PARALLEL INFORMATIONS
        IF(dim_mesh.NE.3) THEN
        IF (code.NE.'  ') THEN
          nameclm = code//'PAR'//extens(nparts-1,i-1)
        ELSE
          nameclm = namecli(1:3)//'PAR'//extens(nparts-1,i-1)
        ENDIF
        IF(.NOT. partel_concat)THEN
          CALL get_free_id(nclm)
          OPEN(nclm,file=nameclm,status='UNKNOWN',form='FORMATTED')
          rewind(nclm)
        ENDIF
        WRITE(nclm,*) nptfr_p(i)
!
! FILE OPENED, NOW WORK ON BOUNDARIES
! -----------------------------------
!
! WHEN THE BOUNDARY NODE BELONGS TO THIS SUBDOMAIN IT WILL BE TAKEN
! J IS THE RUNNING BOUNDARY NODE NUMBER
!
        j = 0
        l = 0
!
        DO k=1,npoin_p(i)
          IF(irand(knolg(k,i)).NE.0) THEN
            l = l + 1
          ENDIF
        ENDDO
        DO k=1,nptfr
!
!         BOUNDARY NODES BELONGING TO THIS PARTITION
!
          temp=hash_table_get(knogl,nbor(k),i)
          IF(temp.NE.0) THEN
            j = j + 1
            iseg = 0
            xseg = 0.0
            yseg = 0.0
!
!           IF THE ORIGINAL (GLOBAL) BOUNDARY LEADS FURTHER INTO
!           ANOTHER PARTITION THEN ISEG IS SET NOT EQUAL TO ZERO
!           THE NEXT NODE ALONG THE GLOBAL BOUNDARY HAS IPTFR = M
!           (BUT CHECK THE CASE THE CIRCLE CLOSES)
!
            m = kp1bor(k,1)
!
!           NBOR_P CANNOT BE USED, IT IS NOT FULLY FILLED WITH DATA
!
            iso = 0
!           CHECKING IF THE ADJACENT ELEMENT IS NOT IN THE
!           SUB-DOMAIN
            IF (epart(nelbor(k)).NE.i) THEN
!
              iseg = nbor(m)
              xseg = REAL(f(iseg,1))
              yseg = REAL(f(iseg,2))
              iso = iso + 1
            ENDIF
!
            m = kp1bor(k,2)
!
!           SAME AS ABOVE, BUT PREVIOUS SEGMENT ,THUS M, NOT K
            IF (epart(nelbor(m)).NE.i) THEN
              iseg = -nbor(m)
              xseg = REAL(f(-iseg,1))
              yseg = REAL(f(-iseg,2))
              iso = iso + 1
            ENDIF
!
!           WHEN BOTH NEIGHBOURS BOUNDARY NODES BELONG TO ANOTHER PARTITION
!
            IF (iso == 2) THEN
              iseg = -9999
              iso = 0
              WRITE(lu,*) 'ISOLATED BOUNDARY POINT', nbor(k), temp
            ENDIF
!
!           WRITE BOUNDARY PARALLEL INFORMATIOSN
!           CONCERNING THE NODE WHICH HAS BEEN RESEARCHED
!
            WRITE (nclm,*) color(k), hash_table_get(knogl,nbor(k),i),
     &                     iseg, xseg, yseg, numliq(k)
            written_lines=written_lines+1
          ENDIF
!
        END DO
!
!       CHECKING THAT THERE IS NOT AN ERROR IN THE BOUNDARIES
        IF(j.NE.nptfr_p(i)) THEN
          WRITE(lu,*) 'ERROR IN BOUNDARIES J=',j
          WRITE(lu,*) 'WHILE NPTFR_P(I)=',nptfr_p(i)
          WRITE(lu,*) 'WHILE L=',l
          CALL plante(1)
          stop
        ENDIF


!
        fmt4='(I7)'
        WRITE (nclm,*) nptir_p(i)
        written_lines=written_lines+1
        IF (max_n_neigh < nbmaxnshare-1) max_n_neigh = nbmaxnshare-1
        fmt4='(   (I7,1X))'
        WRITE (fmt4(2:4),'(I3)') max_n_neigh+1
!
!       SORTING NODE NUMBERS TO SORT(J) SO THAT CUT_P(SORT(J)) IS ORDERED
!       CUT IS OVERWRITTEN NOW
!
        DO j=1,nptir_p(i)
          cut(j)=hash_table_get(cut_p, j, i)
        ENDDO
!
!       IF A NODE HAS BEEN ALREADY FOUND AS MIN, CUT(NODE) GETS 0
!
        DO j=1,nptir_p(i)
          idum = npoin2+1  ! LARGEST POSSIBLE NODE NUMBER + 1
          k=0
 401      CONTINUE
          k = k + 1
          temp=hash_table_get(cut_p,k,i)
          IF ( cut(k) /= 0 .AND. temp < idum ) THEN
            sort(j) = k
            idum = temp
          ENDIF
          IF ( k < nptir_p(i) ) THEN
            GOTO 401
          ELSE
            cut(sort(j)) = 0
          ENDIF
        ENDDO
!
        DO j=1,nptir_p(i)
          tab_tmp=0
          l=0
          DO k=1,max_n_neigh
            temp=hash_table_get(cut_p,sort(j),i)
            IF(part_p(temp,k) .NE. i .AND.part_p(temp,k) .NE. 0) THEN
              l=l+1
              tab_tmp(l)=part_p(temp, k)
            ENDIF
          ENDDO
          WRITE(nclm,fmt=fmt4)
     &      hash_table_get(knogl,temp,i),(tab_tmp(k)-1, k=1,max_n_neigh)
        ENDDO
        written_lines=written_lines+nptir_p(i)
!                                !
        DO j=1,nhalo(i)
          DO m=0,2
            IF (ifapar(2+m,j)>0) THEN
              ifapar(5+m,j)=hash_table_get(gelegl,ifapar(5+m,j),
     &             ifapar(2+m,j))
            END IF
          ENDDO
        ENDDO
        DO j=1,nhalo(i)
          DO m=0,2
            IF (ifapar(2+m,j)>0) THEN
              ifapar(2+m,j)=ifapar(2+m,j)-1
            END IF
          ENDDO
        ENDDO
!
        WRITE(nclm,'(I9)') nhalo(i)
        DO k=1,nhalo(i)
          WRITE (nclm,'(7(I9,1X))') ifapar(:,k)
        END DO
        written_lines=written_lines+nhalo(i)+1
!
        IF(partel_concat)THEN
          offset_end=offset_end+written_lines
          WRITE(nclm_idx,*)offset_begin
          offset_begin=offset_end
        ELSE
          CLOSE(nclm)
        ENDIF
        !
        ENDIF
      END DO
      IF(partel_concat)THEN
        CLOSE(nclm)
        CLOSE(nclm_idx)
      ENDIF
!
      DEALLOCATE(ifapar)
      DEALLOCATE(part_p)
      DEALLOCATE(numliq)
      DEALLOCATE(tab_tmp)
      CALL hash_table_destroy(gelegl)
      DEALLOCATE(cut)
      CALL hash_table_destroy(cut_p)
      DEALLOCATE(sort)
!
      ALLOCATE(ikles_p(max_nelem_p*3),stat=ierr)
      IF(nplan.GT.1) THEN
        ALLOCATE(ikles3d_p(6,max_nelem_p,nplan-1),stat=ierr)
      ENDIF
      CALL check_allocate(ierr, 'IKLES3D_P')
!
      CALL init_dataval(dataval, npoin, nvar, ntimestep, fformat,
     &                  ninp, variable)
!
      CALL write_solutions(fformat,nbor, knogl,ndp_bnd,nelebd,ninp,
     &                     nplan,ikles,times,knolg, date,nelem_p,
     &                     npoin_p,nptfr_p,ubor,hbor,chbord,tbor,vbor,
     &                     btbor,atbor,liubor,lihbor,litbor,livbor,
     &                     dataval,ikle_bnd,elelg,typ_elem,time,title,
     &                     nvar,ntimestep,nptfr,npoin2,nparts,npoin,
     &                     nout,ndp,f,variable,typ_bnd_elem,
     &                     nptir_p,nameinp,namecli,color)

      DEALLOCATE(nbor)
      DEALLOCATE(lihbor)
      DEALLOCATE(liubor)
      DEALLOCATE(livbor)
      DEALLOCATE(hbor)
      DEALLOCATE(ubor)
      DEALLOCATE(vbor)
      DEALLOCATE(chbord)
      DEALLOCATE(litbor)
      DEALLOCATE(tbor)
      DEALLOCATE(atbor)
      DEALLOCATE(btbor)
      DEALLOCATE(dataval)
      DEALLOCATE(color)
!   -------------------------------------------------------------------
! //// JAJ: LA FINITA COMMEDIA FOR PARALLEL CHARACTERISTICS, BYE!
!----------------------------------------------------------------------
! !JAJ #### DEAL WITH SECTIONS
!
      with_sections = namesec(1:1) .NE. ' '
      IF (nplan.NE.0) with_sections=.false.
      IF (with_sections) THEN ! PRESENTLY, FOR TELEMAC2D, EV. SISYPHE
!
        WRITE(lu,*) 'DEALING WITH SECTIONS WITH FILE ',trim(namesec)
        CALL handle_sections(namesec, nparts, nelem, ndp, ikle, npoin,
     &                       f, knogl)
        WRITE(lu,*) 'FINISHED DEALING WITH SECTIONS'
!
      ENDIF ! NPLAN==0

! !YA #### DEAL WITH ZONES
!
!----------------------------------------------------------------------
!
      with_zones = namezfi(1:1) .NE. ' '
      IF (nplan.NE.0) with_zones = .false.
      IF(with_zones) THEN
!
        WRITE(lu,*) 'DEALING WITH ZONES WITH FILE ',trim(namezfi)
        CALL handle_friction_zones(namezfi, nparts, npoin, npoin_p,
     &                             max_npoin_p, knolg)
        WRITE(lu,*) 'FINISHED DEALING WITH ZONES'
!
      ENDIF ! WITH_ZONES
!
! !CCT #### DEAL WITH WEIRS
!
!----------------------------------------------------------------------
!
      with_weirs = nameseu(1:1) .NE. ' '
      IF (nplan.NE.0) with_weirs = .false.
      IF(with_weirs) THEN
!
        WRITE(lu,*) 'DEALING WITH WEIRS WITH FILE ',trim(nameseu)
        CALL handle_weirs(nameseu, nparts, knogl, npoin2)
        WRITE(lu,*) 'FINISHED DEALING WITH WEIRS'
!
      ENDIF ! WITH_WEIRS
!
!----------------------------------------------------------------------
!
      DEALLOCATE (ikle) ! #### MOVED FROM FAR ABOVE
      DEALLOCATE(npart)
      DEALLOCATE(epart)
      DEALLOCATE(npoin_p)
      DEALLOCATE(nelem_p)
      DEALLOCATE(nptfr_p)
      DEALLOCATE(nptir_p)
!
      DEALLOCATE(ikles)
      IF(nplan.GT.1) THEN
        DEALLOCATE(ikles3d)
        DEALLOCATE(ikles3d_p)
      ENDIF
      DEALLOCATE(ikles_p)
      DEALLOCATE(irand)
      DEALLOCATE(f)
!
      DEALLOCATE(knolg)
      CALL hash_table_destroy(knogl)
      DEALLOCATE(elelg)
      DEALLOCATE(kp1bor)
      DEALLOCATE(variable)
!
!----------------------------------------------------------------------
!
      IF (timecount) THEN
        CALL system_clock (count=temps, count_rate=parsec)
        tfin = temps
        WRITE(lu,*) 'OVERALL TIMING: ',
     &    (1.0*(tfin-tdeb))/(1.0*parsec),' SECONDS'
        WRITE(lu,*) ' '
      ENDIF
!
      WRITE(lu,*) '+---- PARTEL: NORMAL TERMINATION ----+'
      WRITE(lu,*) ' '
!
!----------------------------------------------------------------------
!
      END SUBROUTINE partel