C SEPA2     SOURCE    PV        13/11/29    21:15:06     7872           
C****************************************************************************
C****************************************************************************
C***************SEPAR...trouve separation************************************
C****************************************************************************
C****************************************************************************

C    SEPAR   trouve la separation a  partir du domaine defini  par
C    MASQ=.TRUE. et du noeud appele PIVOT, renvoie DIMSEP,le nombre de 
C    noeuds contenant dans la separation, MASQ=.FALSE. pour les noeuds 
C    appartenant a la separation, renumerote celle-ci dans IADJ,

        SUBROUTINE SEPA2(IADJ,JADJC,PIVOTA,MASQUE,DIMSEP,N,IPOS,
     &                   NODES,IPOSMAX,nrelong,noelon,noel2,            
     >                   londim,nbthr)

      IMPLICIT INTEGER(I-N)
      common/cnumop/fcout2c(64),fcoutc(64),iadjc,jadjcc,nrelongc,
     > noelonc,isensc(64),dimlonc(64),masqc,nodesc,iposc,nbtotc,
     > lfrontc(64),londimc,lpivc(64),ipointc,noel2c,
     > boolc(64),ldim2c(64),nbthrc,pivotc(50,64),npointc(64),
     > npoint2c(64),imaxc,lfronc(64)
        external noepr2i
        integer dimlonc,boolc
        SEGMENT IADJ(0)
        SEGMENT JADJC(0)
        integer pivot(50),pivot2(50),pivota,pivots,pivotc,pivotb
        SEGMENT MASQUE
                LOGICAL MASQ(0)
        ENDSEGMENT
        SEGMENT IPOS(0)
        INTEGER NODES
        INTEGER BOOL
        INTEGER DIMSEP,N
        REAL*8 FCOUT,FCOUT2,xbtot,xbfr2,xdeseq2,fcout2c,fcoutc
        real*8 xran
        INTEGER LONG,LONG2,L,DIMLON
C       LONG2,LONG correspond au pseudo-diametre.
C       L  determine les noeuds appartenant a la separation.

        SEGMENT NRELONG(NODES)
C       NRELONG contient pour chaque noeud sa profondeur.

        SEGMENT NOELON(NODES)
        SEGMENT NOEL2(NODES)
        SEGMENT LONDIM(NODES)
C       NOELON contient les noeuds de profondeur LONG.
C       DIMLON= dimension de NOELON.
        
C       Initialisation des segments  de travail.
*       write (6,*) ' entree dans sepa2 nbthr ',nbthr
*  noepr2 travaille differamment suivant que nbtot est nul ou pas
        nbtot=0
        lpiv=1
        pivot(lpiv)=pivota
        MDOMN=IPOS(PIVOTA+NODES)
        CALL NOEPR2(IADJ(1),JADJC(1),PIVOT,NRELONG(1),NOELON(1),
     >  isens,DIMLON,MASQ(1),NODES,IPOS(1),NBTOT,lfront,lfron,
     >  londim(1),fcout,lpiv)
**      nbtot=0
**      pivot(lpiv)=noelon(1)
**      pivotb=pivot(lpiv)
**      do ix=1,dimlon
**        ii=noelon(ix)
**        nrelong(ii)=0
**      enddo
**      CALL NOEPR2(IADJ(1),JADJC(1),PIVOT,NRELONG(1),NOELON(1),
**   >  isens,DIMLON,MASQ(1),NODES,IPOS(1),NBTOT,lfront,lfron,
**   >  londim(1),fcout,lpiv)
        do ix=1,dimlon
          ii=noelon(ix)
          nrelong(ii)=0
        enddo
        if (fcout.lt.2) goto 11
        if (nbtot.le.13) goto 11
        ldim2=londim(2)
        do j=1,dimlon
          noel2(j)=noelon(j)
        enddo
*  recherche aléatoire
        npoint=2
        do i=lpiv+1,npoint
         pivot(i)=0
        enddo
        lpiv=npoint
        bool=2
        ipoint=0
        imax=0
 499    continue
        do 500 itent=1,2**30
         if (bool.le.0 ) goto 510
         bool=bool-1
         ipoint=ipoint+1
         if (ipoint.gt.npoint) then 
           imax=1       
           ipoint=1
         endif
         npoint2=npoint
         nbthrl=min(nbthr,nbtot/16+1)
         pivots=pivot(ipoint)
         do 521 ith=1,nbthrl
          iadjc=iadj
          jadjcc=jadjc
          nrelongc=nrelong
          noelonc=noelon
          isensc(ith)=isens
          dimlonc(ith)=dimlon
          masqc=masque
          nodesc=nodes
          iposc=ipos
          nbtotc=nbtot
          lfrontc(ith)=lfront
          lfronc(ith)=lfron
          londimc=londim
          fcout2c(ith)=fcout2
          lpivc(ith)=lpiv
          ipointc=ipoint
          noel2c=noel2
          boolc(ith)=bool      
          fcoutc(ith)=fcout
          ldim2c(ith)=ldim2
          npointc(ith)=npoint
          npoint2c(ith)=npoint2                        
          nbthrc=nbthrl
          imaxc=imax
          do jj=1,lfron
           londim(jj+(ith-1)*nodes)=londim(jj)
          enddo
          do jj=1,npoint2
            pivotc(jj,ith)=pivot(jj)
          enddo
          if (ith.ne.nbthrl) then
*          write (6,*) ' appel threadid ',ith
           call threadid(ith,noepr2i)
          else
           call noepr2i(ith)
          endif
 521     continue
         do 523 ith=1,nbthrl-1
          call threadif(ith)
 523     continue
         do 522 ith=1,nbthrl
          if (fcoutc(ith).lt.fcout.or.(fcoutc(ith).eq.fcout.and.
     >       pivot(ipoint).ne.pivotc(ipoint,ith))) then
*          write (6,*) ' sepa2 ith fcoutc fcout ',ith,fcoutc(ith),
*    >       fcout,' pivotc ',   pivotc(ipoint,ith)
           if (fcoutc(ith).ne.fcout) bool=2
           isens=isensc(ith)
           fcout=fcoutc(ith)
           ldim2=ldim2c(ith)
           npoint=npointc(ith)
           npoint2=npoint2c(ith)
           pivot(ipoint)=pivotc(ipoint,ith)
           pivot(npoint2)=0           
          endif
 522     continue
         npoint=npoint2
         lpiv=npoint
 500    continue
 510    continue
*       write (6,*) ' nb pts maitres ',npoint
*  petit gradient local
        istop=0
        do 560 itent=1,2**30
        if (istop.eq.1) goto 570
        istop=1
        do 550 ipoint=  npoint,1,-1
         lpoint=pivot(ipoint)
         if (lpoint.le.0) goto 550
*  on ne le fait que si il n'y a pas trop de noeuds connectes
         if (iadj(lpoint+1)-iadj(lpoint).gt.dimlon/5) goto 550
         if (ipoint.ne.1) then
          pivot(ipoint)=0
          CALL NOEPR2(IADJ(1),JADJC(1),PIVOT,NRELONG(1),NOELON(1),
     >     isens,DIMLON,MASQ(1),NODES,IPOS(1),NBTOT,lfront,lfron,
     >     londim(1),fcout2,lpiv)
          if (fcout.lt.2) goto 11
          if (fcout2.lt.fcout) then
*     write (6,*) ' amelioration 2 fcout fcout2 ',fcout,fcout2
           istop=0
           fcout=fcout2
           lpoint=pivot(ipoint)
           goto 560
          endif
         endif
         do 580 kk=iadj(lpoint),iadj(lpoint+1)-1
          k=jadjc(kk)
          if (k.eq.0) goto 580
          if (IPOS(k+NODES).ne.mdomn) goto 580
          pivot(ipoint)=k
          CALL NOEPR2(IADJ(1),JADJC(1),PIVOT,NRELONG(1),NOELON(1),
     >     isens,DIMLON,MASQ(1),NODES,IPOS(1),NBTOT,lfront,lfron,
     >     londim(1),fcout2,lpiv)
          if (fcout.lt.2) goto 11
          if (fcout2.lt.fcout) then
*     write (6,*) ' amelioration 1 fcout fcout2 ',fcout,fcout2
           istop=0
           fcout=fcout2
           lpoint=pivot(ipoint)
*          goto 560
          else
           pivot(ipoint)=lpoint
          endif
 580    continue
 550    continue
*       if (istop.eq.0) goto 499
 560    continue
 570    continue


 11     continue
*       write (6,*) ' sepa2 pivot qivot rivot sivot nbtot ',pivot,qivot,
*    >    rivot,sivot,nbtot
        nbtot=0
        CALL NOEPR2(IADJ(1),JADJC(1),PIVOT,NRELONG(1),NOELON(1),
     >   isens,DIMLON,MASQ(1),NODES,IPOS(1),NBTOT,lfront,lfron,
     >   londim(1),fcout,lpiv)

C       PIVOT correspond au noeud pseudo-peripherique.
C       LONG2 correspond au pseudo-diametre.
C

        DIMSEP=0
C       MDOMN est le numéro de mon domaine

C       pour l'instant,aucun noeud n'appartient a la separation.


        l=lfront
C       L correspond a la distance moyenne pour aller d'un bout a l'autre 
C       du domaine.

C       on arrete de separer si LONG < = 5.
C       on masque alors tous les noeuds repondant a cette condition.

        IF(Lfront.le.1.OR.NBTOT.LE.13) THEN
C       insertion d'une seule zone
          IPOSV=IPOS(MDOMN+1)
C       IPOSV est son rang
          DO I=0,IPOSMAX
           IF (IPOS(I+1).GT.IPOSV) IPOS(I+1)=IPOS(I+1)+1
          ENDDO
         IPOSMAX=IPOSMAX+1
          if (iposmax.gt.nodes-1) then
*          write (6,*) ' iposmax > modes '
           iposmax=nodes-1
          endif
        IPOS(IPOSMAX+1-0)=IPOSV+1
C         IPOSV est son rang
              DO 40 IX=1,DIMLON
              I=NOELON(IX)
                IF(MDOMN.EQ.IPOS(I+NODES).AND.NRELONG(I).NE.0) THEN
                        MASQ(I)=.FALSE.
                        IPOS(I+NODES)=IPOSMAX
                        DIMSEP=DIMSEP+1
                ENDIF
 40                CONTINUE
*        do 41 i=1,nodes
*         IF(MDOMN.EQ.IPOS(I+NODES).AND.NRELONG(I).EQ.0) THEN
*            write (6,*) ' zone non connexe 1 ',mdomn,i
*         endif
*41      continue
*         CALL CMG2(IADJ,JADJC,PIVOT,LONG,NRELONG,NOELON,
*    *               DIMLON,MASQUE,NODES,IPOS,DIMSEP)
                GOTO 50
        ENDIF

C       calcul des nouveaux indices de zones
        IPOSV=IPOS(MDOMN+1)
C       IPOSV est son rang
        DO I=0,IPOSMAX
         IF (IPOS(I+1).GT.IPOSV) IPOS(I+1)=IPOS(I+1)+3
        ENDDO
        IPOSMAX=IPOSMAX+3
          if (iposmax.gt.nodes-1) then
*          write (6,*) ' iposmax > modes '
           iposmax=nodes-1
          endif
        isens=1
        IPOS(IPOSMAX+1-2)=IPOSV+1
        IPOS(IPOSMAX+1-1)=IPOSV+1+isens
        IPOS(IPOSMAX+1-0)=IPOSV+4-isens

C       On calcule la position de chaque noeud I.
C       si I a une profondeur L, I appartient a la separation.

        iwr=1
        izcti1=2**30
        izcti2=2**30
        izcti3=2**30
        izcti4=2**30
        izcti5=2**30
        izcti6=2**30
        izcte1=2**30
        izcte2=2**30
        izcte3=2**30
        izcte4=2**30
        izcte5=2**30
        izcte6=2**30
        dimseq=dimsep
        DO 20 IX=1,DIMLON
         I=NOELON(IX)
*  L'intérieur
             IF(MDOMN.EQ.IPOS(I+NODES).AND.(NRELONG(I).GT.0)) THEN
                IF((NRELONG(I).LT.L)) THEN
*  recherche de izcti
               do kk=iadj(i),iadj(i+1)-1
                 k=jadjc(kk)
               if (k.ne.0) then
         if (nrelong(k).eq.0) then
         if (izcti1.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti1.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=izcti4
                   izcti4=izcti3
                   izcti3=izcti2
                   izcti2=izcti1
                   izcti1=ipos(ipos(k+nodes)+1)
         elseif (izcti2.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti2.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=izcti4
                   izcti4=izcti3
                   izcti3=izcti2
                   izcti2=ipos(ipos(k+nodes)+1)
         ELSEIF (IZCTI3.EQ.IPOS(IPOS(K+NODES)+1)) THEN
         elseif (izcti3.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=izcti4
                   izcti4=izcti3
                   izcti3=ipos(ipos(k+nodes)+1)
         elseif (izcti4.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti4.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=izcti4
                   izcti4=ipos(ipos(k+nodes)+1)
         elseif (izcti5.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti5.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=ipos(ipos(k+nodes)+1)
         elseif (izcti6.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti6.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=ipos(ipos(k+nodes)+1)
         endif
         endif
               endif
               enddo
                        IPOS(I+NODES)=IPOSMAX-1
                        IPOS(I+2*NODES)=IPOS(I+2*NODES)+1
*                 write (6,*) ' zone ',i,IPOS(I+NODES)
                        GOTO 20
                ENDIF

                IF(NRELONG(I).EQ.L) THEN
*  pour que le noeud soit bien dans la separation, il faut qu'il soit
*  connecter a un noeud de profondeur plus élevée (prob des noeuds milieux)
                    DO 100 J=1,IADJ(I+1)-IADJ(I)
                       K=JADJC(IADJ(I)+J-1)
C  K:voisin de I
*                      if (k.eq.0) write (6,*) ' k peut valoir 0 '
                       IF(K.EQ.0) GOTO 100
                       IF (NRELONG(K).EQ.0) GOTO 100
                       IF(MDOMN.NE.IPOS(k+NODES).and.
     >                    ipos(k+nodes).lt.iposmax-2) then
                          write (6,*) ' sepa2 mauvais domaine ',mdomn,
     >                    ipos(k+nodes),nrelong(k)
                          GOTO 100
                       endif
                       IF(NRELONG(K).GT.L) GOTO 110
 100                CONTINUE
*  recherche de izcti
               do kk=iadj(i),iadj(i+1)-1
                 k=jadjc(kk)
               if (k.ne.0) then
         if (nrelong(k).eq.0) then
         if (izcti1.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti1.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=izcti4
                   izcti4=izcti3
                   izcti3=izcti2
                   izcti2=izcti1
                   izcti1=ipos(ipos(k+nodes)+1)
         elseif (izcti2.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti2.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=izcti4
                   izcti4=izcti3
                   izcti3=izcti2
                   izcti2=ipos(ipos(k+nodes)+1)
         elseif (izcti3.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti3.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=izcti4
                   izcti4=izcti3
                   izcti3=ipos(ipos(k+nodes)+1)
         elseif (izcti4.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti4.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=izcti4
                   izcti4=ipos(ipos(k+nodes)+1)
         elseif (izcti5.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti5.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=izcti5
                   izcti5=ipos(ipos(k+nodes)+1)
         elseif (izcti6.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcti6.gt.ipos(ipos(k+nodes)+1)) then
                   izcti6=ipos(ipos(k+nodes)+1)
         endif
         endif
               endif
               enddo
                    IPOS(I+NODES)=IPOSMAX-1
                        IPOS(I+2*NODES)=IPOS(I+2*NODES)+1
*                 write (6,*) ' zone ',i,IPOS(I+NODES)
                    GOTO 20
                ENDIF
*
*  L'extérieur
                IF(NRELONG(I).GT.L) THEN
*  recherche de izcte
               do kk=iadj(i),iadj(i+1)-1
                 k=jadjc(kk)
               if (k.ne.0) then
         if (nrelong(k).eq.0) then
         if (izcte1.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcte1.gt.ipos(ipos(k+nodes)+1)) then
                   izcte6=izcte5
                   izcte5=izcte4
                   izcte4=izcte3
                   izcte3=izcte2
                   izcte2=izcte1
                   izcte1=ipos(ipos(k+nodes)+1)
         elseif (izcte2.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcte2.gt.ipos(ipos(k+nodes)+1)) then
                   izcte6=izcte5
                   izcte5=izcte4
                   izcte4=izcte3
                   izcte3=izcte2
                   izcte2=ipos(ipos(k+nodes)+1)
         elseif (izcte3.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcte3.gt.ipos(ipos(k+nodes)+1)) then
                   izcte6=izcte5
                   izcte5=izcte4
                   izcte4=izcte3
                   izcte3=ipos(ipos(k+nodes)+1)
         elseif (izcte4.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcte4.gt.ipos(ipos(k+nodes)+1)) then
                   izcte6=izcte5
                   izcte5=izcte4
                   izcte4=ipos(ipos(k+nodes)+1)
         elseif (izcte5.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcte5.gt.ipos(ipos(k+nodes)+1)) then
                   izcte6=izcte5
                   izcte5=ipos(ipos(k+nodes)+1)
         elseif (izcte6.eq.ipos(ipos(k+nodes)+1)) then
         elseif (izcte6.gt.ipos(ipos(k+nodes)+1)) then
                   izcte6=ipos(ipos(k+nodes)+1)
         endif
         endif
               endif
               enddo

                        IPOS(I+NODES)=IPOSMAX
                        IPOS(I+2*NODES)=IPOS(I+2*NODES)+1
*                 write (6,*) ' zone ',i,IPOS(I+NODES)
                        GOTO 20
                ENDIF
                GOTO 21
*  la frontiere
 110                    CONTINUE
*                 write (6,*) ' zone ',i,IPOS(I+NODES)
                        MASQ(I)=.FALSE.
                        IPOS(I+NODES)=IPOSMAX-2
                        DIMSEP=DIMSEP+1
                        GOTO 20
            ENDIF
* mars
21      CONTINUE

20      CONTINUE
*        do 23 i=1,nodes
*          IF(MDOMN.EQ.IPOS(I+NODES).AND.NRELONG(I).EQ.0) THEN
*          write (6,*) ' zone non connexe 2 ',mdomn,i
*          endif
* 23     continue


*       write (6,*) ' taille de frontière ',dimsep-dimseq
*       write (6,*) ' rang des zones ',izcti,izcte
        ipos1=max(IPOS(IPOSMAX+1-1),IPOS(IPOSMAX+1+0))
        ipos2=min(IPOS(IPOSMAX+1-1),IPOS(IPOSMAX+1+0))
        if (izcti1.gt.izcte1) then
          IPOS(IPOSMAX+1-1)=ipos1
          IPOS(IPOSMAX+1+0)=ipos2
        elseif (izcti1.lt.izcte1) then
          IPOS(IPOSMAX+1-1)=ipos2
          IPOS(IPOSMAX+1+0)=ipos1
        elseif (izcti2.gt.izcte2) then
          IPOS(IPOSMAX+1-1)=ipos1
          IPOS(IPOSMAX+1+0)=ipos2
        elseif (izcti2.lt.izcte2) then
          IPOS(IPOSMAX+1-1)=ipos2
          IPOS(IPOSMAX+1+0)=ipos1
        elseif (izcti3.gt.izcte3) then
          IPOS(IPOSMAX+1-1)=ipos1
          IPOS(IPOSMAX+1+0)=ipos2
        elseif (izcti3.lt.izcte3) then
          IPOS(IPOSMAX+1-1)=ipos2
          IPOS(IPOSMAX+1+0)=ipos1
        elseif (izcti4.gt.izcte4) then
          IPOS(IPOSMAX+1-1)=ipos1
          IPOS(IPOSMAX+1+0)=ipos2
        elseif (izcti4.lt.izcte4) then
          IPOS(IPOSMAX+1-1)=ipos2
          IPOS(IPOSMAX+1+0)=ipos1
        elseif (izcti5.gt.izcte5) then
          IPOS(IPOSMAX+1-1)=ipos1
          IPOS(IPOSMAX+1+0)=ipos2
        elseif (izcti5.lt.izcte5) then
          IPOS(IPOSMAX+1-1)=ipos2
          IPOS(IPOSMAX+1+0)=ipos1
        elseif (izcti6.gt.izcte6) then
          IPOS(IPOSMAX+1-1)=ipos1
          IPOS(IPOSMAX+1+0)=ipos2
        elseif (izcti6.lt.izcte6) then
          IPOS(IPOSMAX+1-1)=ipos2
          IPOS(IPOSMAX+1+0)=ipos1
        else 
*         write (6,*) ' ordonnancement impossible',
*    > izcti1,izcti2,izcti3,izcti4,izcti5,izcti6,
*    > izcte1,izcte2,izcte3,izcte4,izcte5,izcte6
        endif

50      continue
         do ix=1,dimlon
           ii=noelon(ix)
           nrelong(ii)=0
         enddo

        RETURN
        END