partel_para.F

Go to the documentation of this file.

!                   *******************
                    PROGRAM partel_para
!                   *******************
!
!***********************************************************************
! PARALLEL   V6P2                                   20/02/2012
!***********************************************************************
!
!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  JAJ; MATTHIEU GONZALES DE LINARES
!+        27/01/2003
!+        SIXTH VERSION
!+    CORRECTED A WRONG DIMENSION OF THE ARRAY ALLVAR
!
!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.
!
!history   C. MOULINEC, P. VEZOLLE, O. BOITEAU
!+
!+      NEXT VERSION
!+      SOME CHANGES
!
!history   C. DENIS (EDF-SINETICS)
!+        31/05/2010
!+      V6P2
!+      FURTHER DEVELOPPED IN ORDER TO DECREASE THE AMOOUT OF MEMORY
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/02/2012
!+        V6P2
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET
!+        27/03/2014
!+    Arguments added to the call elebd.
!+    Note that all the history before this is wrong...
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_parallel
      USE declarations_special
      USE bief, ONLY : bief_ncsize => ncsize, nbmaxnshare, ipid, front2
      USE interface_parallel, ONLY : p_time, p_max_array
      IMPLICIT NONE
!
      ! MAX PARTITION NUMBER [00000..99999]
      INTEGER, PARAMETER :: maxnproc = 100000
      ! MAX ADDED SUFFIX LENGTH
      INTEGER, PARAMETER :: maxaddch = 10
      ! MAX NUMBER OF VARIABLES
      INTEGER, PARAMETER :: maxvar = 100
      ! MAXVAR*32 FOR ALLVAR
      INTEGER, PARAMETER :: maxallvarlength = 3200
!
      INTEGER pmethod
      INTEGER nvar, nplan, nptfr, nptir, nptfrmax
      INTEGER nelem, npoin, ndp, nelem2, npoin2, ndum
      INTEGER ib(10)
!
      INTEGER, ALLOCATABLE :: ikles(:), ikles_p(:)
      INTEGER, ALLOCATABLE :: ikles3d(:),ikles3d_p(:,:,:)
      INTEGER, ALLOCATABLE :: irand(:), irand_p(:)
      INTEGER, ALLOCATABLE :: lihbor(:), liubor(:), livbor(:)
      INTEGER, ALLOCATABLE :: litbor(:)
      INTEGER ::  npoin_p, nelem_p , nptfr_p,nptir_p
      INTEGER, ALLOCATABLE :: nbor(:), nbor_p(:)
      INTEGER, ALLOCATABLE :: numliq(:)
      INTEGER, ALLOCATABLE :: knolg(:), knogl(:),check(:)
      INTEGER, ALLOCATABLE :: elelg(:)
      INTEGER, ALLOCATABLE :: cut(:), cut_p(:), sort(:)
      INTEGER, ALLOCATABLE :: part_p(:,:)
!
      REAL, ALLOCATABLE    :: f(:,:), f_p(:,:)
      REAL, ALLOCATABLE    :: hbor(:)
      REAL, ALLOCATABLE    :: ubor(:), vbor(:), aubor(:)
      REAL, ALLOCATABLE    :: tbor(:), atbor(:), btbor(:)
!
      REAL times
!
      INTEGER :: ninp, ncli, nout, nclm
      INTEGER time(3), date(3)
!
      CHARACTER(LEN=80)  :: title
      CHARACTER(LEN=32)  :: vari, variable(maxvar)
      CHARACTER(LEN=PATH_LEN)  :: nameinp, namecli, nameout, nameclm
      CHARACTER(LEN=15)  :: fmt4
!
      INTEGER  max_n_neigh
      INTEGER i, j, k, l , m, err, iso
      INTEGER istop, istart, iseg, iloop
      INTEGER i_len, i_s, i_sp, i_lencli, i_leninp
!
      REAL xseg, yseg
      LOGICAL is, timecount
!
!   METISOLOGY
!
      INTEGER nparts
      INTEGER, ALLOCATABLE :: epart(:), npart(:)
!
!   FOR CALLING FRONT2
!
      ! MAX NUMBER OF BOUNDARIES
      INTEGER, PARAMETER :: maxfro = 3000
      INTEGER nfrliq, nfrsol, debliq(maxfro), finliq(maxfro)
      INTEGER debsol(maxfro), finsol(maxfro)
      INTEGER, ALLOCATABLE :: dejavu(:), kp1bor(:,:)
      INTEGER, ALLOCATABLE :: nachb(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: x_d(:), y_d(:)
!
!   FOR CALLING BIEF MESH SUBROUTINES (TO BE OPTIMISED SOON):
!
      INTEGER, ALLOCATABLE :: ifabor(:,:), ifanum(:,:), nelbor(:)
      INTEGER, ALLOCATABLE :: nulone(:,:)
      INTEGER, ALLOCATABLE :: ikle(:,:), iklbor(:,:), isegf(:)
      INTEGER, ALLOCATABLE :: it1(:), it2(:), it3(:)
!
!   TIME MEASURING
!
      INTEGER  tdebp, tfinp, temps, parsec
!     EXTENS
      CHARACTER(LEN=11) :: extens
!     TENTATIVELY CALLING ONLY FOR SERAFIN ? VERIFY
      CHARACTER(LEN=8),PARAMETER  :: fformat='SERAFIN '
      EXTERNAL extens
!
!----------------------------------------------------------------------
!
!   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
!
      DOUBLE PRECISION :: starttime,endtime,startiotime,startctime
      INTEGER, ALLOCATABLE :: gelegl(:),gelegl1(:)
!
      ! THE HALO ELEMENTS NEIGHBOURHOOD DESCRIPTION FOR A HALO CELL
      INTEGER, ALLOCATABLE :: ifapar(:,:)
!
      ! THE NUMBER OF HALO CELLS PRO PARTITION
      INTEGER :: nhalo
!
      ! WORK VARIABLES
      INTEGER ifaloc(3)
      INTEGER ndp_2d,nb_inter_glob,compt
      INTEGER ef
      INTEGER, ALLOCATABLE :: nbre_ef(:),nbre_ef_loc(:),ef_i(:),
     &     tab_tmp(:),ef_ii(:),global_inter_node_reord(:)
      INTEGER, ALLOCATABLE :: part_p_tmp1(:),part_p_tmp2(:)
      LOGICAL halo
      INTEGER noeud,nbre_noeud_interne
      INTEGER nbre_ef_i,ier,nbre_nptir
      LOGICAL interface
!
!   FOR SCOTCH
!
      CHARACTER(LEN=MAXLENTMPDIR) :: path
      INTEGER :: ncar, ncsize
      INTEGER :: id_input
      ndp_2d=3
      CALL p_init(path, ncar, ipid, ncsize)
!
!----------------------------------------------------------------------
!
      timecount = .true.
      IF (parsec==0) timecount = .false.  ! COUNT_RATE == 0 : NO CLOCK
!
      lng=2    ! ENGLISH PLEASE
      lu=6     ! FORTRAN STANDARD OUPUT CHANNEL
!
!----------------------------------------------------------------------
!   INTRODUCE YOURSELF
!
      IF (ipid .EQ. 0) THEN
!
        WRITE(lu,*) ' '
        WRITE(lu,*) '+----------------------------------------------+'
        WRITE(lu,*) ' '
        WRITE(lu,*) '  PARTEL //'
        WRITE(lu,*) '  PARALLEL VERSION DEVELOPPED BY'
        WRITE(lu,*) '  CHRISTOPHE DENIS (SINETICS)'
        WRITE(lu,*) ' '
        WRITE(lu,*) '  PARTEL'
        WRITE(lu,*) '  BUNDESANSTALT FUER WASSERBAU, KARLSRUHE'
        WRITE(lu,*) ' '
        WRITE(lu,*) '  METIS 5.0.2'
        WRITE(lu,*) '  REGENTS OF THE UNIVERSITY OF MINNESOTA'
        WRITE(lu,*) ' '
        WRITE(lu,*) '+----------------------------------------------+'
!
      END IF
!
!----------------------------------------------------------------------
!   NAMES OF THE INPUT FILES:
!
      INQUIRE (file='PARTEL.PAR',exist=is)
      IF( .NOT.is ) THEN
        WRITE (lu,'('' FILE DOES NOT EXIST: '',A30)') 'PARTEL.PAR'
        CALL plante(1)
        stop
      END IF
!
      CALL get_free_id(id_input)
      OPEN(unit=id_input,file='partel.par')
      READ(id_input,*) nameinp
      READ(id_input,*) namecli
      READ(id_input,*) nparts
      READ(id_input,*) pmethod
      ! SI PMETHOD = 3 or 4 ON NE  DECOUPE PAS AVEC METIS ON LIT UN FICHIER
      ! RESULT_SEQ_METIS CREE PAR PARTEL_PRELIM (OU AUTRE !!)
      ! SINON ON DECOUPE AVEC METIS
      CLOSE(id_input)
!
      IF (pmethod .EQ. 3) THEN
        WRITE(lu,*) ' '
        WRITE(lu,*) 'PARTITIONING USING PARMETIS'
        WRITE(lu,*) ' '
      ELSEIF (pmethod .EQ. 4) THEN
        WRITE(lu,*) ' '
        WRITE(lu,*) 'PARTITIONING USING PTSCOTCH'
        WRITE(lu,*) ' '
      ENDIF
!
      INQUIRE (file=nameinp,exist=is)
      IF( .NOT.is ) THEN
        WRITE (lu,'('' FILE DOES NOT EXIST: '',A30)') nameinp
        CALL plante(1)
        stop
      END IF
!
      INQUIRE (file=namecli,exist=is)
      IF (.NOT.is) THEN
        WRITE (lu,'('' FILE DOES NOT EXIST: '',A30)') namecli
        CALL plante(1)
        stop
      END IF
!
      i_s  = len(nameinp)
      i_sp = i_s + 1
      DO i=1,i_s
        IF (nameinp(i_sp-i:i_sp-i) .NE. ' ') EXIT
      ENDDO
      i_len=i_sp - i
      i_leninp = i_len
!
      IF (i_leninp > path_len) THEN
        WRITE(lu,*) ' '
        WRITE(lu,*) 'ATTENTION:'
        WRITE(lu,*) 'THE NAME OF THE INPUT FILE:'
        WRITE(lu,*) nameinp
        WRITE(lu,*) 'IS LONGER THAN ',path_len,' CHARACTERS'
        WRITE(lu,*) 'WHICH IS THE LONGEST APPLICABLE NAME FOR TELEMAC'
        WRITE(lu,*) 'INPUT AND OUTPUT FILES. STOPPED. '
        CALL plante(1)
        stop
      ENDIF
!
      IF (ipid .EQ. 0)  THEN
        starttime=p_time()
      END IF
      CALL get_free_id(ninp)
      OPEN(ninp,file=nameinp,status='OLD',form='UNFORMATTED')
      rewind ninp

      READ (ninp) title
      READ (ninp) i, j
      nvar = i + j
      istart = 42
      DO i=1,nvar
        READ(ninp) vari
        variable(i) = vari
        DO j=1,32
          IF(vari(j:j).EQ.' ') vari(j:j) = '-'
        END DO
        istop = istart+20
        IF (istop.GT.maxallvarlength) THEN
          WRITE(lu,*) 'VARIABLE NAMES TOO LONG FOR STRING ALLVAR'
          WRITE(lu,*) 'STOPPED.'
          CALL plante(1)
          stop
        ENDIF
        istart=istop+1
      ENDDO
!
!   READ THE REST OF THE SELAFIN FILE
!   10 INTEGERS, THE FIRST IS THE NUMBER OF RECORDS (TIMESTEPS)
!
      READ (ninp) (ib(i), i=1,10)
      IF (ib(8).NE.0.OR.ib(9).NE.0) THEN
        WRITE(lu,*) 'THIS IS A PARTIAL OUTPUT FILE'
        WRITE(lu,*) 'MAYBE MEET GRETEL BEFORE...'
      ENDIF
      nplan = ib(7)
!
      IF (ib(10).EQ.1) THEN
        READ(ninp) date(1), date(2), date(3), time(1), time(2), time(3)
      ENDIF
!
      READ (ninp) nelem,npoin,ndp,ndum
      IF (nplan.GT.1) THEN
!
        npoin2 = npoin/nplan
        nelem2 = nelem/(nplan-1)
        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
      ELSE
!
        npoin2 = npoin
        nelem2 = nelem
      ENDIF
!
      IF (ndp.EQ.3) THEN
        WRITE(lu,*) 'THE INPUT FILE ASSUMED TO BE 2D SELAFIN'
      ELSEIF (ndp.EQ.6) THEN
        WRITE(lu,*) 'THE INPUT FILE ASSUMED TO BE 3D SELAFIN'
      ELSE
!       WRITE(LU,*) 'THE ELEMENTS ARE NEITHER TRIANGLES NOR PRISMS!'
!       WRITE(LU,*) 'NDP = ',NDP
        CALL plante(1)
        stop
      ENDIF
!
!   NOW LET US ALLOCATE
!
      ALLOCATE (ikles(nelem2*3),stat=err)
      CALL check_allocate(err, 'IKLES')
      IF(nplan.GT.1) THEN
        ALLOCATE (ikles3d(nelem*ndp),stat=err)
        CALL check_allocate(err, 'IKLES3D')
      ENDIF
      ALLOCATE (irand(npoin),stat=err)
      CALL check_allocate(err, 'IRAND')
      ! NVAR+2 : FIRST TWO FUNCTIONS ARE X AND Y
      ! NPOIN IS 3D HERE IN 3D
      ALLOCATE (f(npoin,nvar+2),stat=err)
      CALL check_allocate(err, 'IRAND')
!
!   CONNECTIVITY TABLE:
!
      IF(nplan.EQ.0) THEN
        READ(ninp) ((ikles((k-1)*ndp+j),j=1,ndp),k=1,nelem)
      ELSE
        READ(ninp) ((ikles3d((k-1)*ndp+j),j=1,ndp),k=1,nelem)
!       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
!
      READ(ninp) (irand(j),j=1,npoin)
!
!
      READ(ninp) (f(j,1),j=1,npoin)
      READ(ninp) (f(j,2),j=1,npoin)
!
      iloop = 0
!      DO
!
!   READ THE TIME STEP
!
        READ(ninp) times
!         WRITE(NINPFORMAT,*) TIMES
        iloop = iloop + 1
!
!
! READ THE TIME VARIABLES; NO 1 AND 2 ARE X,Y
!
        DO k=3,nvar+2
!            WRITE(LU,*) 'NOW READING VARIABLE',K-2
          READ(ninp) (f(j,k), j=1,npoin)
!            WRITE(NINPFORMAT,*) (F(J,K), J=1,NPOIN)
!            WRITE(LU,*) 'READING VARIABLE',K-2,' SUCCESSFUL'
        END DO
!      END DO
      CLOSE (ninp)
      CALL p_sync()
      IF (ipid .EQ. 0) THEN
        startctime=p_time()
        WRITE(lu,*) 'TEMPS LECTURE :',startctime-starttime
      END IF

!
!----------------------------------------------------------------------
! READ THE BOUNDARY CONDITIONS FILE
!
!   BUT ALLOCATE FIRST
!
      nptfrmax = npoin2   ! BETTER IDEA ?
!
      ALLOCATE (lihbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'LIHBOR')
      ALLOCATE (liubor(nptfrmax),stat=err)
      CALL check_allocate(err, 'LIUBOR')
      ALLOCATE (livbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'LIVBOR')
      ALLOCATE (hbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'HBOR')
      ALLOCATE (ubor(nptfrmax),stat=err)
      CALL check_allocate(err, 'UBOR')
      ALLOCATE (vbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'VBOR')
      ALLOCATE (aubor(nptfrmax),stat=err)
      CALL check_allocate(err, 'AUBOR')
      ALLOCATE (tbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'TBOR')
      ALLOCATE (atbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'ATBOR')
      ALLOCATE (btbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'BTBOR')
      ALLOCATE (litbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'LITBOR')
      ALLOCATE (nbor(nptfrmax),stat=err)
      CALL check_allocate(err, 'NBOR')
      ALLOCATE (numliq(nptfrmax),stat=err)
      CALL check_allocate(err, 'NUMLIQ')
      ALLOCATE (check(nptfrmax),stat=err)
      CALL check_allocate(err, 'CHECK')
!
!   CORE NAME LENGTH
!
      i_s  = len(namecli)
      i_sp = i_s + 1
      DO i=1,i_s
        IF (namecli(i_sp-i:i_sp-i) .NE. ' ') EXIT
      ENDDO
      i_len=i_sp - i
      i_lencli = i_len
!
      IF (i_leninp > path_len) THEN
        WRITE(lu,*) ' '
        WRITE(lu,*) 'ATTENTION:'
        WRITE(lu,*) 'THE NAME OF THE BOUNDARY CONDITIONS FILE:'
        WRITE(lu,*) namecli
        WRITE(lu,*) 'IS LONGER THAN ',path_len,' CHARACTERS'
        WRITE(lu,*) 'WHICH IS THE LONGEST APPLICABLE NAME FOR TELEMAC '
        WRITE(lu,*) 'INPUT AND OUTPUT FILES. STOPPED. '
        CALL plante(1)
        stop
      ENDIF
!
      CALL get_free_id(ncli)
      OPEN(ncli,file=namecli,status='OLD',form='FORMATTED')
      rewind ncli
!
!   READING BOUNDARY FILE AND COUNTING BOUNDARY POINTS
!
      k=1
 900  CONTINUE
      READ(ncli,*,end=901,err=901) lihbor(k),liubor(k),
     &                             livbor(k),
     &             hbor(k),ubor(k),vbor(k),aubor(k),litbor(k),
     &             tbor(k),atbor(k),btbor(k),nbor(k),check(k)
!
!   NOW CHECK IS THE BOUNDARY NODE COLOUR
!     IF(CHECK(K).NE.K) THEN
!        WRITE(LU,*) 'ERROR IN BOUNDARY CONDITIONS FILE AT LINE ',K
!        CALL PLANTE(1)
!        STOP
!     ENDIF
      k=k+1
      GOTO 900
 901  CONTINUE
      nptfr = k-1
!     WRITE (LU,*) ' '
!      WRITE (LU,*) 'NUMBER OF BOUNDARY NODES IN 2D MESH: ',NPTFR
!      WRITE (LU,*) ' '
      CLOSE(ncli)
!
!----------------------------------------------------------------------
!   NUMBERING OF OPEN BOUNDARIES
!   NUMBERING OF LIQUID BOUNDARY, IF 0 = SOLID
!   OPN: NUMBER OF OPEN BOUNDARY
!   IN ORDER TO DO IT IN THE SAME WAY AS TELEMAC DOES,
!   IT IS BEST TO CALL FRONT2 HERE
!
!   FOR CALLING BIEF MESH SUBROUTINES
!   CAN BE OPTIMISED / USES A LOT OF MEMORY
!   THE ONLY REASON IS TO OBTAIN KP1BOR AND NUMLIQ
!
      ALLOCATE (dejavu(nptfr),stat=err)
      CALL check_allocate(err, 'DEJAVU')
      ALLOCATE (kp1bor(nptfr,2),stat=err)
      CALL check_allocate(err, 'KP1BOR')
!JAJ----------V ////
!     CHANGED NELEM TO NELEM2, NDP TO 3 HUH!
!     CAUSING ERRORS WHEN 3D RESTART/REFERENCE FILES ARE PARTITIONED
!     AND BC FILE IS WRITTEN AGAIN (WHAT FOR, ACTUALLY???)
!     CAUSE: CALLING VOISIN WITH NELEM2 BUT IFABOR(NELEM=NELEM3,NDP=6)
      ALLOCATE (ifabor(nelem2,3),stat=err)
      CALL check_allocate(err, 'IFABOR')
      ALLOCATE (ifanum(nelem2,3),stat=err)
      CALL check_allocate(err, 'IFANUM')
      ALLOCATE (iklbor(nptfr,2),stat=err)
      CALL check_allocate(err, 'IKLBOR')
      ALLOCATE (nelbor(nptfr),stat=err)
      CALL check_allocate(err, 'NELBOR')
      ALLOCATE (nulone(nptfr,2),stat=err)
      CALL check_allocate(err, 'NULONE')
      ALLOCATE (isegf(nptfr),stat=err)
      CALL check_allocate(err, 'ISEGF')
      ALLOCATE (ikle(nelem2,3),stat=err)
      CALL check_allocate(err, 'IKLE')
      ALLOCATE (it1(npoin),stat=err)
      CALL check_allocate(err, 'IT1')
      ALLOCATE (it2(npoin),stat=err)
      CALL check_allocate(err, 'IT2')
      ALLOCATE (it3(npoin),stat=err)
      CALL check_allocate(err, 'IT3')
      nptir = 1
      ALLOCATE (nachb(nbmaxnshare,nptir),stat=err)
      CALL check_allocate(err, 'NACHB')
      ALLOCATE (x_d(npoin2),stat=err)
      CALL check_allocate(err, 'X_D')
      ALLOCATE (y_d(npoin2),stat=err)
      CALL check_allocate(err, 'Y_D')
!
!   TRANSFORM IKLES--> IKLE FOR 2D ROUTINES  (AN OLD TELEMAC DISEASE)
!
      DO i = 1,3
        DO j  = 1,nelem2
          ikle(j,i) = ikles((j-1)*3+i)
        ENDDO
      ENDDO
!
!     Forcong voisn elebd and front2 to run on serila mode
      bief_ncsize = 1
      CALL voisin(ifabor, nelem2, nelem2, 11, ikle, nelem2,
     &                   npoin2, nachb, nbor, nptfr, it1, it2)
!
!      WRITE(LU,'(/,'' CALLING ELEBD'')')
!
      CALL elebd (nelbor, nulone, kp1bor, ifabor, nbor(1:nptfr), ikle,
     &                   nelem2, iklbor, nelem2, nelem2,
     &                   npoin2, nptfr, 11, lihbor, 2,
     &                   isegf, it1, it2, it3 ,
     &                   nptfr ,nptfr)
!                        NELEBX,NELEB (EQUAL NPTFR HERE)
!
!      WRITE(LU,'(/,'' BOUNDARY TYPE NUMBERING USING FRONT2'')')
!
      nfrliq = 0
      nfrsol = 0
      debsol(:) = 0
      finsol(:) = 0
      debliq(:) = 0
      finliq(:) = 0
      DO i=1,npoin2
        x_d(i) = dble(f(i,1))
        y_d(i) = dble(f(i,2))
      ENDDO
      CALL front2 (nfrliq,
     &             lihbor,liubor,x_d,y_d,
     &             nbor,kp1bor(1:nptfr,1),dejavu,npoin2,nptfr,
     &             2,.false.,numliq,maxfro)
!
      DEALLOCATE (dejavu)
      DEALLOCATE (nachb)
      DEALLOCATE (x_d)
      DEALLOCATE (y_d)
!JAJ //// IFABOR APPLIED LATER FOR FINDING HALO CELL NEIGHBOURHOODS
!!!!      DEALLOCATE (IFABOR)
      DEALLOCATE (ifanum)
      DEALLOCATE (iklbor)
!     DEALLOCATE (NELBOR)
      DEALLOCATE (nulone)
      DEALLOCATE (isegf)
      DEALLOCATE (ikle)
      DEALLOCATE (it1)
      DEALLOCATE (it2)
      DEALLOCATE (it3)
!
!======================================================================
!   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
!======================================================================
!
!     NEW METIS INTERFACE (>= VERSION 5) :
!
!     EPTR, EIND: THESE ARRAYS SPECIFIES CONTAINS
!     IKLES IN CSR FORMAT
!
      ALLOCATE (epart(nelem2),stat=err)
      CALL check_allocate(err, 'EPART')
      ALLOCATE (npart(npoin2),stat=err)
      CALL check_allocate(err, 'NPART')
!
      ! If metis was run previously just read the results
      IF (pmethod.NE.3 .AND. pmethod.NE.4) THEN
        OPEN(id_input,file='RESULT_SEQ_METIS')
        DO i=1,nelem2
          READ(id_input,*) epart(i)
        END DO
        CLOSE(id_input)
      ELSE
!
!       Choosing partitionning method
!       1 : METIS IF PARMETIS NOT INSTALLED PARMETIS OTHERWISE
!       2 : PT SCOTCH
!
        WRITE(lu,*) ' THE MESH PARTITIONING STEP  STARTS'
        IF (timecount) THEN
          CALL system_clock (count=temps, count_rate=parsec)
          tdebp = temps
        ENDIF
        CALL partitioner_para(pmethod, nelem2, npoin2, ndp, nparts,
     &                           ikles, epart, npart)
        WRITE(lu,*) ' THE MESH PARTITIONING STEP HAS FINISHED'
        IF (timecount) THEN
          CALL system_clock (count=temps, count_rate=parsec)
          tfinp = temps
          WRITE(lu,*) ' RUNTIME OF PARTITIONER ',
     &              (1.0*(tfinp-tdebp))/(1.0*parsec),' SECONDS'
        ENDIF
      ENDIF
!

!======================================================================
!   STEP 3 : ALLOCATE THE GLOBAL  ARRAYS NOT DEPENDING OF THE PARTITION
!
!======================================================================
!
!      KNOGL(I) =>  GLOBAL LABEL OF THE LOCAL POINT I
      ALLOCATE (knogl(npoin2),stat=err)
      CALL check_allocate(err, 'KNOGL')
      knogl(:)=0

!   NBRE_EF(I) => NUMBER OF FINITE ELEMENT CONTAINING I
!   I IS A GLOBAL LABEL
      ALLOCATE (nbre_ef(npoin2),stat=err)
      CALL check_allocate(err, 'NBRE_EF')
!
      IF(nplan.EQ.0) THEN
        ALLOCATE (f_p(npoin2,nvar+2),stat=err)
      ELSE
        ALLOCATE (f_p(npoin2*nplan,nvar+2),stat=err)
      ENDIF
      CALL check_allocate(err, 'F_P')
!
      ALLOCATE (gelegl(nelem2),stat=err)
      CALL check_allocate(err, 'GELEGL')
      ALLOCATE (gelegl1(nelem2),stat=err)
      CALL check_allocate(err, 'GELEGL1')
!
      ALLOCATE (sort(npoin2),stat=err)
      CALL check_allocate(err, 'CUT_P')
!
      ALLOCATE (cut(npoin2),stat=err)
      CALL check_allocate(err, 'CUT_P')
      sort=0
!
      ALLOCATE(nbor_p(npoin2),stat=err)
      ALLOCATE(tab_tmp( nbmaxnshare),stat=err)
!
!======================================================================
!   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
      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
!
!     THE PROCESS MPI IPID  WORK ON THE SUBMESH IPID+1
      i=ipid+1
!
!   LOOP OVER THE FINITE ELEMENTD 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
        IF (epart(ef) .EQ. i) THEN
          nelem_p=nelem_p+1
          DO k=1,ndp_2d
            noeud=ikles((ef-1)*3+k)
            IF (knogl(noeud) .EQ. 0) THEN
              npoin_p=npoin_p+1
              knogl(noeud)=npoin_p
            END IF
          END DO
        END IF
      END DO
!
!======================================================================
!     STEP 4 : ALLOCATION OF LOCAL ARRAYS NEEDED BY MPI PROCESSUS IPID
!              WORKING ON SUBMESH IPID+1
!======================================================================
!
!   ELEGL(E) => GLOBAL LABEL OF THE FINITE ELEMENT E
!   E IS THE LOCAL LABEL ON SUBMESH I
      ALLOCATE (elelg(nelem_p),stat=err)
      CALL check_allocate(err, 'ELELG')
      elelg(:)=0
!   KNOLG(I) => GLOBAL LABEL OF THE POINT I
!   I IS THE LOCAL LABEL ON SUBDOMAIN I
      IF(nplan.EQ.0) THEN
        ALLOCATE (knolg(npoin_p),stat=err)
      ELSE
        ALLOCATE (knolg(npoin_p*nplan),stat=err)
      ENDIF
      CALL check_allocate(err, 'KNOLG')
      knolg(:)=0
!   NBRE_EF_LOC(I) : NUMBER OF FINITE ELEMENTS CONTAINING THE POINT I
!                    ON SUBMESH I
!   I IS THE LOCAL LABEL ON SUBMESH I
      ALLOCATE (nbre_ef_loc(npoin_p),stat=err)
      CALL check_allocate(err, 'NBRE_EF_LOC')
!
!   EF_I(E) IS THE GLOBAL LABEL OF THE INTERFACE FINITE ELEMENT NUMBER E
      ALLOCATE (ef_i(nelem_p),stat=err)
      CALL check_allocate(err, 'EF_I')
!   EF_II(E) IS THE LOCAL LABEL OF THE INTERFACE FINITE ELEMENT NUMBER E
      ALLOCATE (ef_ii(nelem_p),stat=err)
      CALL check_allocate(err, 'EF_II')
!
!======================================================================
!     STEP 5 : INITIALISATION  OF LOCAL ARRAYS
!                  (GELELG AND ELELG, NBRE_EF_LOC)
!
!======================================================================
!
      nelem_p=0
      gelegl(:)=0
      DO ef=1,nelem2
        IF (epart(ef) .EQ. i) THEN
          nelem_p=nelem_p+1
          elelg(nelem_p)=ef
          gelegl(ef)=nelem_p
        END IF
      END DO
      CALL p_isum_array(gelegl,gelegl1,nelem2,ier)
      DEALLOCATE(gelegl)
      DEALLOCATE(npart)
      nptfr_p=0
      DO j=1,npoin_p
        nbre_ef_loc(j)=0
      END DO
!
!======================================================================
!     STEP 5 : COMPUTE THE NUMBER OF BOUNDARY AND INTERFACE POINTS
!              INITIALISATION OF NBRE_EF_LOC AND F_P
!
!======================================================================
!
      npoin_p=0
      nbre_noeud_interne=0
!
      DO j=1,nelem_p
        ef=elelg(j)
        DO k=1,3
          noeud=ikles((ef-1)*3+k)
          nbre_ef_loc(knogl(noeud))=nbre_ef_loc(knogl(noeud))+1
          IF (nbre_ef_loc(knogl(noeud)) .EQ. 1) THEN
!   THE POINT NOEUD IS ENCOUNTERED FOR THE FIRST TIME
            npoin_p=npoin_p+1
!   IS NOEUD A BOUNDARY POINT ?
            IF (irand(noeud) .NE. 0) THEN
              nptfr_p= nptfr_p+1
            END IF
!   MODIFICATION OF   KNOGL ET F_P
            knolg(npoin_p)=noeud
            DO l=1,nvar+2
              f_p(npoin_p,l)=f(noeud,l)
            END DO
          END IF
!
!   NOEUD IS A INTERNAL POINT IF ALL FINITE ELEMENTS
!   CONTAINING IT BELONGS TO THE SAME SUBMESH
          IF (nbre_ef_loc(knogl(noeud)) .EQ. nbre_ef(noeud)) THEN
            nbre_noeud_interne=nbre_noeud_interne+1
          END IF
        END DO
      END DO
!
!======================================================================
!     STEP 5 : FIND THE INTERFACES POINTS AND THE INTERFACE FINITE
!              ELEMENTS OF SUBMESH I
!======================================================================
!
!   COMPUTE THE NUMBER OF INTERFACES POINTS TO ALLOCATE CORRECLY  CUT_P
      nbre_nptir=npoin_p-nbre_noeud_interne
      ALLOCATE (cut_p(nbre_nptir),stat=err)
      CALL check_allocate(err, 'CUT_P')
      nbre_ef_i=0 ! THE NUMBER OF INTERFACE FINITE ELEMENT IS SET TO 0
      nptir_p=0 !THE NUMBER OF INTERFACE POINTS IS SET TO 0
!   LOOP ON THE FINITE ELEMENT BELONGING TO SUBMESH I
      DO j=1,nelem_p
        interface=.false.
        ef=elelg(j)         ! GLOBAL LABEL OF THE FINITE ELEMENT
        DO k=1,ndp_2d       ! LOOP ON THE NODES CONTAINED IN EF
          noeud=ikles((ef-1)*3+k)
          IF (abs(nbre_ef_loc(knogl(noeud))) .NE. nbre_ef(noeud))
     &         THEN
            ! THE GLOBAL NUMBER OF FINITE ELEMENTS CONTAINING NOEUD
            ! IS GREATER THAN THE LOCAL NUMBER OF FINITE ELEMENTS
            ! CONTAINING NOEUD => NOEUD IS AN INTERFACE POINT
            !    => EF IS AN INTERFACE FINITE ELEMENT
            interface=.true.
          END IF
          IF (nbre_ef_loc(knogl(noeud)) .NE.  nbre_ef(noeud).AND.
     &         nbre_ef_loc(knogl(noeud)) .GT. 0) THEN
            ! IT IS THE FIRST TIME THAT NOEUD IS DETECTED AS INTERFACE POINT
            ! AS AFTER  NBRE_EF_LOC(KNOGL(NOEUD))= -1*NBRE_EF_LOC(KNOGL(NOEUD))
            interface=.true.
            nptir_p=nptir_p+1
            cut_p(nptir_p)=noeud
            sort(noeud)=1
            nbre_ef_loc(knogl(noeud))=
     &            -1*nbre_ef_loc(knogl(noeud))
          END IF
        END DO
        IF (INTERFACE .EQV. .true.) THEN
          nbre_ef_i=nbre_ef_i+1
          ef_i(nbre_ef_i)=ef
          ef_ii(nbre_ef_i)=j
        END IF
      END DO
      WRITE(22,*) 'NPTIR_P',nptir_p
      nb_inter_glob=0
      ! IF (IPID. EQ. 0) WRITE(89,*) SORT
      ! CUT=0
      CALL p_max_array(sort,cut,npoin2,ier)
!
      DO j=1,npoin2
        IF (cut(j) .NE. 0) THEN
          nb_inter_glob=nb_inter_glob+1
        END IF
      END DO
      ALLOCATE(global_inter_node_reord(nb_inter_glob),stat=err)
      CALL check_allocate(err, 'GLOBAL_INTER_REORD')
      ALLOCATE(part_p(nb_inter_glob,nbmaxnshare+1),stat=err)
      CALL check_allocate(err, 'PART_P')
      ALLOCATE(part_p_tmp1(nparts+1),stat=err)
      CALL check_allocate(err, 'PART_P_TMP1')
      ALLOCATE(part_p_tmp2(nparts+1),stat=err)
      CALL check_allocate(err, 'PART_P_TMP2')
      nb_inter_glob=0
      DO j=1,npoin2
        IF (cut(j) .NE. 0) THEN
          nb_inter_glob=nb_inter_glob+1
          global_inter_node_reord( nb_inter_glob)=j
!            IF (NBRE_EF_LOC(KNOGL(NOEUD)) .LT. 0) THEN
!               PART_P(NB_INTER_GLOB,I)=I
!            END IF
        END IF
      END DO
      part_p=0
      DO j=1,nb_inter_glob
        part_p_tmp1=0
        part_p_tmp2=0
        compt=0
        noeud=global_inter_node_reord(j)
        IF (sort(noeud) .NE. 0) THEN
          part_p_tmp1(i)=i
          part_p_tmp1(nparts+1)=part_p_tmp1(nparts+1)+1
        END IF
        CALL p_isum_array(part_p_tmp1,part_p_tmp2,nparts+1,
     &        ier)
        DO k=1,nparts
          IF (part_p_tmp2(k) .NE. 0 .AND. k .NE. i) THEN
            compt=compt+1
            part_p(j,compt)=k
            part_p(j,nbmaxnshare+1)=part_p_tmp2(nparts+1)
          END IF
        END DO
      END DO
      DEALLOCATE(part_p_tmp1)
      DEALLOCATE(part_p_tmp2)
      DEALLOCATE(cut)
      max_n_neigh=maxval(part_p(:,nbmaxnshare+1))
!      WRITE(LU,*) 'MAX', MAX_N_NEIGH
      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!'
      END IF
      IF (max_n_neigh > 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 < nbmaxnshare-1) max_n_neigh = nbmaxnshare-1
!
!======================================================================
!     STEP 6 : COMPUTE THE HALO
!
!======================================================================
!
      nhalo=0
!   LOOP ON THE INTERFACE FINITE ELEMENT
      DO j=1,nbre_ef_i
        ef=ef_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 .EQV. .true.) THEN
          nhalo=nhalo+1
        END IF
      END DO
      ALLOCATE (ifapar(7,nhalo),stat=err)
      CALL check_allocate(err, 'IFAPAR')
      ifapar(:,:)=0
!
      nhalo=0
!   LOOP ON THE INTERFACE FINITE ELEMENT
      DO j=1,nbre_ef_i
        ef=ef_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 .EQV. .true.) THEN
          nhalo=nhalo+1
          ifapar(1,nhalo)=ef_ii(j)
          ifapar(2:4,nhalo)=ifaloc(:)
          ifapar(5:7,nhalo)=ifabor(ef_i(j),:)
!            DO K=1,3
!               NPART(IFALOC(K))=NPART(IFALOC(K))+1
!            END DO
        END IF
      END DO
!
!       WRITE(LU,*) 'SOUS DOMAINE ',I,'NBRE POINTS',NPOIN_P,
!     &        'NBRE NOEUD INTE',NBRE_NOEUD_INTERNE,'INTERFACE',
!!     &        NPTIR_P,'NBRE FRONT',NPTFR_P, 'HALO',NHALO,
!     &        'NBRE_EFRONT',NBRE_EF_I
!
      IF (.NOT. ALLOCATED(nbor_p)) THEN
        ALLOCATE(nbor_p(npoin2),stat=err)
        CALL check_allocate(err, 'NBOR_P')
      END IF
!     END DO

!     DO I=1,NPARTS
      CALL p_sync()
      IF (ipid .EQ. 0) THEN
        startiotime=p_time()
        WRITE(lu,*) 'TEMPS CALCUL :',startiotime-startctime
      END IF
!
!-----------------------------------------------------------------------
!  THE CORE NAMES FOR THE OUTPUT BC FILES ACCORDING TO THE NUMBER OF PARTS
!
      nameclm = namecli    ! CORE NAME LENGTH IS I_LENCLI
      nameout = nameinp    ! CORE NAME LENGTH IS I_LENINP
!
!----------------------------------------------------------------------
!
      nameclm(i_lencli+1:i_lencli+11) = extens(nparts-1,ipid)
!
      CALL get_free_id(nclm)
      OPEN(nclm,file=nameclm,
     &        status='UNKNOWN',form='FORMATTED')
      rewind(nclm)
!
!   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
!
!      NPTIR = 0
      j = 0
!
      DO k=1,nptfr
!
!   BOUNDARY NODES BELONGING TO THIS PARTITION
!
        IF ( knogl(nbor(k)) /= 0) THEN
          j = j + 1
!            NBOR_P(J) = NBOR(K)
          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 = f(iseg,1)
            yseg = 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 = f(-iseg,1)
            yseg = f(-iseg,2)
            iso = iso + 1
          ENDIF
!
!     WHEN BOTH NEIGHBOURS BOUNDARY NODES BELONG TO ANOTHER PARTITION
!
          IF (iso == 2) THEN
            iseg = -9999
            iso = 0
          ENDIF
!
          WRITE (nclm,4000)
     &         lihbor(k), liubor(k), livbor(k),
     &         hbor(k), ubor(k), vbor(k),
     &          aubor(k), litbor(k), tbor(k), atbor(k), btbor(k),
     &          nbor(k),check(k), iseg, xseg, yseg, numliq(k)

 4000       FORMAT (1x,i2,1x,2(i1,1x),3(f24.12,1x),1x,
     &           f24.12,3x,i1,1x,3(f24.12,1x),1i15,1x,1i15,
     &           1x,i15,1x,2(f27.15,1x),i15)
        ENDIF
!
      END DO

      fmt4='(I15)'
      WRITE (nclm,*) nptir_p
!      WRITE(LU,*) '######### ', NPTIR_P
      IF (max_n_neigh < nbmaxnshare-1) max_n_neigh = nbmaxnshare-1
      fmt4='(   (I15,1X))'
      WRITE (fmt4(2:4),'(I3)') max_n_neigh+1
!
      DO j=1, nb_inter_glob
        noeud=global_inter_node_reord(j)
        IF (sort(noeud) .NE. 0) THEN

          WRITE(nclm,fmt=fmt4) global_inter_node_reord(j),
     &        (part_p(j,k)-1, k=1,max_n_neigh)
        END IF
      END DO
!
      DO j=1,nhalo
        DO m=0,2
          IF (ifapar(2+m,j)>0) THEN
            ifapar(5+m,j)=gelegl1(ifapar(5+m,j))
          END IF
        END DO
      END DO
!      IFAPAR(2:4,:)=IFAPAR(2:4,:)-1
      DO j=1,nhalo
        DO m=0,2
          IF (ifapar(2+m,j)>0) THEN
            ifapar(2+m,j)=ifapar(2+m,j)-1
          END IF
        END DO
      END DO
!
      WRITE(nclm,'(I15)') nhalo
      DO k=1,nhalo

        WRITE (nclm,'(7(I15,1X))') ifapar(:,k)
      END DO
      CLOSE(nclm)
!
      ALLOCATE (irand_p(npoin2),stat=err)
      CALL check_allocate(err, 'IRAND_P')
!
      CALL check_allocate(err, 'F_P')
      ALLOCATE(ikles_p(nelem2*3),stat=err)
      IF(nplan.GT.1) THEN
        ALLOCATE(ikles3d_p(6,nelem2,nplan-1),stat=err)
      ENDIF
      CALL check_allocate(err, 'IKLES3D_P')
!
!***************************************************************
!     WRITING GEOMETRY FILES FOR ALL PARTS/PROCESSORS
!
      nameout(i_leninp+1:i_leninp+11) = extens(nparts-1,ipid)
!
      CALL get_free_id(nout)
      OPEN(nout,file=nameout,form='UNFORMATTED'
     &     ,status='UNKNOWN')
      rewind(nout)
!
!   TITLE, THE NUMBER OF VARIABLES
      WRITE(nout) title
      WRITE(nout) nvar,0
!
      DO k=1,nvar
        WRITE(nout) variable(k)
      END DO
!     WRITE(*,*) 'HERE '
! 10 INTEGERS...
! 1.  IS THE NUMBER OF RECORDINGS IN FILES
! 8.  IS THE NUMBER OF BOUNDARY POINTS (NPTFR_P)
! 9.  IS THE NUMBER OF INTERFACE POINTS (NPTIR_P)
! 10. IS 0 WHEN NO DATE PASSED; 1 IF A DATE/TIME RECORD FOLLOWS
!
!      IB(7) = NPLAN   (ALREADY DONE)
      ib(8) = nptfr_p
      ib(9) = nptir_p
!      WRITE(*,*) 'SD',I,'NPTFR_P(I)',NPTFR_P, 'NPTIR_P(I)',
!     &       NPTIR_P,'NHALO(I)',NHALO,NPOIN_P
      WRITE(nout) (ib(k), k=1,10)
      IF (ib(10).EQ.1) THEN
        WRITE(nout) date(1), date(2), date(3),
     &               time(1), time(2), time(3)
      ENDIF
      IF(nplan.LE.1) THEN
        WRITE(nout) nelem_p, npoin_p, ndp, ndum
      ELSE
        WRITE(nout) nelem_p*(nplan-1),
     &               npoin_p*nplan, ndp, ndum
      ENDIF
      DO j=1,nelem_p
        ef=elelg(j)
        DO k=1,3
          ikles_p((j-1)*3+k) = knogl(ikles((ef-1)*3+k))
        END DO
      END DO
      IF(nplan > 1) THEN
        DO k = 1,nplan-1
          DO j = 1,nelem_p
            ikles3d_p(1,j,k) = ikles_p(1+(j-1)*3) + (k-1)*npoin_p
            ikles3d_p(2,j,k) = ikles_p(2+(j-1)*3) + (k-1)*npoin_p
            ikles3d_p(3,j,k) = ikles_p(3+(j-1)*3) + (k-1)*npoin_p
            ikles3d_p(4,j,k) = ikles_p(1+(j-1)*3) +  k   *npoin_p
            ikles3d_p(5,j,k) = ikles_p(2+(j-1)*3) +  k   *npoin_p
            ikles3d_p(6,j,k) = ikles_p(3+(j-1)*3) +  k   *npoin_p
          ENDDO
        ENDDO
      ENDIF
!
      IF (nplan.LE.1) THEN
        WRITE(nout)
     &     ((ikles_p((j-1)*3+k),k=1,3),j=1,nelem_p)
      ELSE
        WRITE(nout)
     &     (((ikles3d_p(l,j,k),l=1,6),j=1,nelem_p),k=1,nplan-1)
      ENDIF
!
!   INSTEAD OF IRAND, KNOLG IS WRITTEN !!!
!   I.E. THE TABLE PROCESSOR-LOCAL -> PROCESSOR-GLOBAL NODE NUMBERS
!
      IF (nplan.EQ.0) THEN
        WRITE(nout) (knolg(j), j=1,npoin_p)
      ELSE
!   BEYOND NPOIN_P(I) : DUMMY VALUES IN KNOLG, NEVER USED
        WRITE(nout) (knolg(j), j=1,npoin_p*nplan)
      ENDIF
!
! COMPLETING PLANES OTHER THAN 1 FOR FUNCTIONS F_P
      IF (nplan.GT.1) THEN
        DO l=2,nplan
          DO j=1,npoin_p
            DO k=1,nvar+2
              f_p(j+(l-1)*npoin_p,k) = f(knolg(j)+(l-1)*npoin2,k)
            END DO
          ENDDO
        ENDDO
      ENDIF
!
!   NODE COORDINATES X AND Y
!
      IF (nplan.EQ.0) THEN
        WRITE(nout) (f_p(j,1),j=1,npoin_p)
        WRITE(nout) (f_p(j,2),j=1,npoin_p)
      ELSE
        WRITE(nout) (f_p(j,1),j=1,npoin_p*nplan)
        WRITE(nout) (f_p(j,2),j=1,npoin_p*nplan)
      ENDIF
!
!   TIME STAMP (SECONDS)
      WRITE(nout) times
!   NOW THE TIME-DEPENDENT VARIABLES
!
      DO k=3,nvar+2
        IF(nplan.EQ.0) THEN
          WRITE(nout) (f_p(j,k),j=1,npoin_p)
        ELSE
          WRITE(nout) (f_p(j,k),j=1,npoin_p*nplan)
        ENDIF
      END DO
      CLOSE (nout)
      CALL p_sync()
!
!======================================================================
!   WRITING EPART AND NPART
!
! //// JAJ: LA FINITA COMMEDIA FOR PARALLEL CHARACTERISTICS, BYE!
!----------------------------------------------------------------------
!
!   NOTE BY J-M HERVOUET : DEALLOCATE CAUSES ERRORS ON HP
!   (POSSIBLE REMAINING BUG ?)
!   NOTE BY JAJ: DEALLOCATE(HP) ,^)
!
!     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(IRAND_P)
!     DEALLOCATE(F)
!     DEALLOCATE(F_P)
!
!     DEALLOCATE(KNOLG)
!     DEALLOCATE(KNOGL)
!     DEALLOCATE(ELELG)
!
!     DEALLOCATE(LIHBOR)
!     DEALLOCATE(LIUBOR)
!     DEALLOCATE(LIVBOR)
!     DEALLOCATE(HBOR)
!     DEALLOCATE(UBOR)
!     DEALLOCATE(VBOR)
!     DEALLOCATE(LITBOR)
!     DEALLOCATE(AUBOR)
!     DEALLOCATE(TBOR)
!     DEALLOCATE(ATBOR)
!     DEALLOCATE(BTBOR)
!     DEALLOCATE(NBOR)
!     DEALLOCATE(NUMLIQ)
!     DEALLOCATE(KP1BOR)
!     DEALLOCATE(CHECK)
!
!     DEALLOCATE(NBOR_P)
!     DEALLOCATE(CUT)
!     DEALLOCATE(CUT_P)
!     DEALLOCATE(PART_P)
!     DEALLOCATE(SORT)
!
!----------------------------------------------------------------------
!
      CALL p_sync()
      IF (ipid .EQ. 0) THEN
        ENDTIME=P_TIME()
        WRITE(lu,*) 'TEMPS ECRITURE :',endtime-startiotime
        WRITE(lu,*) 'TEMPS TOTAL :',endtime-starttime
        WRITE(lu,*) '+---- PARTEL: NORMAL TERMINATION ----+'
      END IF
!
      CALL p_exit()
      stop 0
!
      END PROGRAM partel_para