Initial commit from Polytechnique Montreal

1 files changed, 619 insertions, 0 deletions

diff --git a/Dragon/src/NXTLCU.f b/Dragon/src/NXTLCU.f
new file mode 100644
index 0000000..b6d43bc
--- /dev/null
+++ b/Dragon/src/NXTLCU.f
@@ -0,0 +1,619 @@
+*DECK NXTLCU
+      SUBROUTINE NXTLCU(IPRINT,MXSLIN,ICOMB ,
+     >                  NBCOR ,NBSINT,ISINT ,TRKLSI)
+*
+*----------
+*
+*Purpose:
+* To merge two sets of tracks and store the result in tracking
+* vector with adequate region identification.
+*
+*Copyright:
+* Copyright (C) 2005 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version.
+*
+*Author(s): G. Marleau.
+*
+*Parameters: input
+* IPRINT  print level.
+* MXSLIN  maximum number of segments in a subgeometry track.
+* ICOMB   flag for combination if two sets of tracks with:
+*         ICOMB=-1  one considers only the first
+*         track set (Cartesian); ICOMB=-2 one considers
+*         only the second track set (annular); ICOMB=1, 
+*         the Cartesian and annular regions are super imposed while
+*         preserving the outer Cartesian boundary;
+*         ICOMB=2 the Cartesian and annular regions
+*         are super imposed while
+*         preserving the outer annular boundary.
+*
+*Parameters: input/output
+* NBCOR   number of corner found for each external faces.
+* NBSINT  number of surface crossed by track.
+* ISINT   direction of plane intersected and
+*         the surfaces crossed by the track.
+* TRKLSI  the surface intersection distance.
+*
+*Reference:
+*  G. Marleau,
+*  New Geometries Processing in DRAGON: The NXT: Module,
+*  Report IGE-260, Polytechnique Montreal,
+*  Montreal, 2005.
+*
+*----------
+*
+      IMPLICIT         NONE
+*----
+*  Subroutine arguments
+*----
+      INTEGER          IPRINT,MXSLIN,ICOMB,NBCOR(2,3),NBSINT(3)
+      INTEGER          ISINT(0:5,MXSLIN,3)
+      DOUBLE PRECISION TRKLSI(MXSLIN,3)
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='NXTLCU')
+      DOUBLE PRECISION DCUTOF,DZERO,DONE,DTWO
+      PARAMETER       (DCUTOF=1.0D-7,DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0)
+*----
+*  Local variables
+*----
+      INTEGER          IT,IC,IA,IG,ISURC,JSURC,ISURA,ISURG,JSURG,
+     >                 IFIC,ILIC,ITIC,IFIA,ILIA,ITIA,ISUR,JSUR,ITB,ITE
+      INTEGER          IDIR,KDIR,JDIR,LDIR,IRADG,IFSG
+      DOUBLE PRECISION REFLOC,REFLOA,REFLOG,CURLOC,CURLOA,
+     >                 THICK,DELCUR
+*----
+*  Data
+*----
+      CHARACTER        NAMTYP(3)*12
+      SAVE             NAMTYP
+      DATA             NAMTYP
+     >               /'Cartesian   ','Annular     ','Combined    '/
+*----
+*  Print header if required
+*----
+      IF(IPRINT .GT. 1000) THEN
+        WRITE(IOUT,6000) NAMSBR
+        ITB=1
+        ITE=2
+        IF(ICOMB .EQ. -1) ITE=1
+        IF(ICOMB .EQ. -2) ITB=2
+        DO IT=ITB,ITE
+          WRITE(IOUT,6016) NAMTYP(IT)
+          WRITE(IOUT,6011) 'Initial face'
+          DO JSUR=1,NBCOR(1,IT)
+            WRITE(IOUT,6010) TRKLSI(JSUR,IT),
+     >      (ISINT(IDIR,JSUR,IT),IDIR=1,5)
+          ENDDO
+          WRITE(IOUT,6011) 'Regions '
+          DO JSUR=NBCOR(1,IT)+1,NBSINT(IT)-NBCOR(2,IT)+1
+            WRITE(IOUT,6010) TRKLSI(JSUR,IT),
+     >      (ISINT(IDIR,JSUR,IT),IDIR=1,5)
+          ENDDO
+          WRITE(IOUT,6011) 'Final face'
+          DO JSUR=NBSINT(IT)-NBCOR(2,IT)+2,NBSINT(IT)+1
+            WRITE(IOUT,6010) TRKLSI(JSUR,IT),
+     >      (ISINT(IDIR,JSUR,IT),IDIR=1,5)
+          ENDDO
+        ENDDO
+      ENDIF
+      IFIA=1
+      ILIA=1
+      LDIR=1
+      JSURC=0
+      ISURG=0
+*----
+*  Create combined tracking
+*----
+      IC=1
+      IA=2
+      IG=3
+      NBCOR(1,IG)=0
+      NBCOR(2,IG)=0
+      IF(ABS(ICOMB) .EQ. 1) THEN
+*----
+*  Save startup Cartesian face
+*----
+        ISURG=0
+        IFIC=1
+        ILIC=NBCOR(1,IC)
+        DO ISURC=IFIC,ILIC
+          ISURG=ISURG+1
+          TRKLSI(ISURG,IG)=TRKLSI(ISURC,IC)
+          DO IDIR=1,5
+            ISINT(IDIR,ISURG,IG)=ISINT(IDIR,ISURC,IC)
+          ENDDO
+        ENDDO
+        NBCOR(1,IG)=ISURG
+        ISURC=NBCOR(1,IC)
+        REFLOC=TRKLSI(ISURC,IC)
+        ISURA=NBCOR(1,2)
+        REFLOA=TRKLSI(ISURA,IA)
+        KDIR=0
+        IRADG=0
+        IFIC=NBCOR(1,IC)+1
+        ILIC=NBSINT(IC)-NBCOR(2,IC)+1
+        REFLOG=REFLOC
+        DELCUR=ABS(REFLOC-REFLOA)
+        IF(ICOMB .EQ. 1) THEN
+*----
+*  Find radial surface direction
+*----
+          DO IDIR=1,5
+            IF(ISINT(IDIR,ISURA,IA) .LT. 0) THEN
+              KDIR=IDIR
+            ENDIF
+          ENDDO
+          IF(KDIR .EQ. 4) THEN
+*----
+*  It is an inner face,
+*  Scan Cartesian regions
+*----
+            ITIC=0
+            DO ISURC=IFIC,ILIC
+              CURLOC=REFLOC+TRKLSI(ISURC,IC)
+              ISURG=ISURG+1
+              DO JDIR=1,5
+                ISINT(JDIR,ISURG,IG)=ISINT(JDIR,ISURC,IC)
+              ENDDO
+              TRKLSI(ISURG,IG)=CURLOC-REFLOG
+              ISINT(4,ISURG,IG)=IRADG
+              IF(ICOMB .EQ. 1) THEN
+                IF(REFLOA .LE. CURLOC) THEN
+                  IRADG=ISINT(4,ISURA+1,IA)
+                  THICK=REFLOA-REFLOG
+                  TRKLSI(ISURG,IG)=THICK
+                  REFLOG=REFLOA
+                  GO TO 205
+                ENDIF
+              ENDIF
+              ITIC=ITIC+1
+              REFLOC=CURLOC
+              REFLOG=REFLOC
+            ENDDO
+ 205        CONTINUE
+            IFIC=IFIC+ITIC
+          ELSE IF(DELCUR .LT. DCUTOF) THEN
+*----
+*  It is an outer face,
+*  reset radial id for Cartesian face
+*----
+            REFLOG=REFLOC
+            IRADG=ISINT(4,1,IA)
+            DO ISURC=1,NBCOR(1,IC)
+              IF(ISINT(KDIR,ISURC,IC) .LT. 0) GO TO 105
+            ENDDO
+            CALL XABORT (NAMSBR//
+     >      ': Cartesian and annular face not coherent')
+ 105        CONTINUE
+            DO JSURG=1,NBCOR(1,IG)
+              ISINT(4,JSURG,IG)=IRADG
+            ENDDO
+          ELSE
+            CALL XABORT(NAMSBR//
+     >      ': Initial Cartesian and annular faces are incoherent')
+          ENDIF
+          IFIA=NBCOR(1,IA)+1
+          ILIA=NBSINT(IA)-NBCOR(2,IA)+1
+        ELSE
+          REFLOG=REFLOC
+        ENDIF
+*----
+*  Scan Cartesian regions
+*----
+        DO ISURC=IFIC,ILIC
+          CURLOC=REFLOC+TRKLSI(ISURC,IC)
+*          IF(IPRINT .GT. 1000)
+*     >    write(6,7000) 'ISURC=',ISURC,CURLOC,
+*     >    (ISINT(JDIR,ISURC,IC),JDIR=1,5)
+          ISURG=ISURG+1
+          DO JDIR=1,5
+            ISINT(JDIR,ISURG,IG)=ISINT(JDIR,ISURC,IC)
+          ENDDO
+          TRKLSI(ISURG,IG)=CURLOC-REFLOG
+          ISINT(4,ISURG,IG)=IRADG
+          IF(ICOMB .EQ. 1) THEN
+            ITIA=0
+*----
+*  Scan annular regions for intermediate mesh
+*----
+            DO ISURA=IFIA,ILIA
+              CURLOA=REFLOA+TRKLSI(ISURA,IA)
+              DELCUR=ABS(CURLOC-CURLOA)
+*              IF(IPRINT .GT. 1000)
+*     >        write(6,7000) 'ISURA=',ISURA,CURLOA,
+*     >        (ISINT(JDIR,ISURA,IA),JDIR=0,5)
+              IF(DELCUR .LT. DCUTOF) THEN
+                ITIA=ITIA+1
+                IF(ISURA .EQ. ILIA) THEN
+                  IRADG=0
+                ELSE
+                  IRADG=ISINT(4,ISURA+1,IA)
+                ENDIF
+*                IRADG=ISINT(4,ISURA,IA)
+                DO JDIR=1,3
+                  IF(ISINT(JDIR,ISURA,IA) .NE. 0) THEN
+                    ISINT(JDIR,ISURG,IG)=ISINT(JDIR,ISURA,IA)
+                  ENDIF
+                ENDDO
+                IF(ISINT(5,ISURA,IA) .NE. 0) THEN
+                  ISINT(5,ISURG,IG)=ISINT(5,ISURA,IA)
+                ENDIF
+*                IF(IPRINT .GT. 1000)
+*     >          write(6,7001) 'CURLOA .EQ. CURLOC',ISURG,
+*     >          (ISINT(JDIR,ISURG,IG),JDIR=1,5)
+                THICK=CURLOA-REFLOG
+                REFLOA=CURLOA
+                REFLOG=REFLOA
+                GO TO 115
+              ELSE IF(CURLOA .GT. CURLOC) THEN
+*                IF(IPRINT .GT. 1000)
+*     >          write(6,7001) 'CURLOA .GT. CURLOC=',ISURG,
+*     >          (ISINT(JDIR,ISURG,IG),JDIR=1,5)
+                GO TO 115
+              ELSE
+                ITIA=ITIA+1
+                IF(ISURA .EQ. ILIA) THEN
+                  IRADG=0
+                ELSE
+                  IRADG=ISINT(4,ISURA+1,IA)
+                ENDIF
+                DO JDIR=1,3
+                  ISINT(JDIR,ISURG+1,IG)=ISINT(JDIR,ISURG,IG)
+                ENDDO
+                ISINT(5,ISURG+1,IG)=ISINT(5,ISURG,IG)
+              ENDIF
+              THICK=CURLOA-REFLOG
+              ISINT(4,ISURG+1,IG)=IRADG
+              TRKLSI(ISURG+1,IG)=TRKLSI(ISURG,IG)-THICK
+              TRKLSI(ISURG,IG)=THICK
+              ISURG=ISURG+1
+              REFLOA=CURLOA
+              REFLOG=REFLOA
+*          IF(IPRINT .GT. 1000)
+*     >    write(6,7001) 'CURLOA .LT. CURLOC',ISURG,
+*     >    (ISINT(JDIR,ISURG,IG),JDIR=1,5)
+            ENDDO
+ 115        CONTINUE
+            IFIA=IFIA+ITIA
+          ENDIF
+          REFLOC=CURLOC
+          REFLOG=REFLOC
+        ENDDO
+*----
+*  Save final Cartesian faces
+*----
+        IFIC=NBSINT(IC)-NBCOR(2,IC)+2
+        ILIC=NBSINT(IC)+1
+        IFSG=ISURG
+        DO ISURC=IFIC,ILIC
+          ISURG=ISURG+1
+          TRKLSI(ISURG,IG)=TRKLSI(ISURC,IC)
+          DO IDIR=1,5
+            ISINT(IDIR,ISURG,IG)=ISINT(IDIR,ISURC,IC)
+          ENDDO
+          ISINT(4,ISURG,IG)=IRADG
+        ENDDO
+        NBCOR(2,IG)=NBCOR(2,IC)
+*----
+*  Find radial zone associated with final Cartesian face
+*  when cylinder face coincides with final face
+*----
+        IF(ICOMB .EQ. 1) THEN
+*----
+*  Find radial surface direction
+*----
+          ISURC=NBSINT(IC)+1
+          REFLOC=TRKLSI(ISURC,IC)
+          ISURA=NBSINT(IA)+1
+          REFLOA=TRKLSI(ISURA,IA)
+          DELCUR=ABS(REFLOC-REFLOA)
+          DO IDIR=1,5
+            IF(ISINT(IDIR,ISURA,IA) .LT. 0) THEN
+              KDIR=IDIR
+            ENDIF
+          ENDDO
+          IF(KDIR .NE. 4) THEN
+*----
+*  It is an inner face,
+*  Scan Cartesian regions
+*----
+            IF(DELCUR .LT. DCUTOF) THEN
+*----
+*  It is an outer face,
+*  reset radial id for Cartesian face
+*----
+              REFLOG=REFLOC
+              IRADG=ISINT(4,ISURA,IA)
+              JSURG=IFSG
+              DO ISURC=IFIC,ILIC
+                JSURG=JSURG+1
+                ISINT(4,JSURG,IG)=IRADG
+              ENDDO
+            ELSE
+              WRITE(IOUT,9002) ISURA,KDIR,REFLOC,REFLOA
+              CALL XABORT(NAMSBR//
+     >        ': Final Cartesian and annular faces are incoherent')
+            ENDIF
+          ENDIF
+          IFIA=NBCOR(1,IA)+1
+          ILIA=NBSINT(IA)-NBCOR(2,IA)+1
+        ELSE
+          REFLOG=REFLOC
+        ENDIF
+      ELSE IF(ABS(ICOMB) .EQ. 2) THEN
+*----
+*  Find radial surface direction
+*----
+        ISURG=0
+        ISURA=NBCOR(1,IA)
+        REFLOA=TRKLSI(ISURA,IA)
+        REFLOG=REFLOA
+        IF(ICOMB .EQ. 2) THEN
+          DO IDIR=1,5
+            IF(ISINT(IDIR,ISURA,IA) .LT. 0) THEN
+              KDIR=IDIR
+            ELSE IF(ISINT(IDIR,ISURA,IA) .GT. 0) THEN
+              LDIR=IDIR
+            ENDIF
+          ENDDO
+          IF(KDIR .EQ. 4) THEN
+*----
+*  It is an outer radial face,
+*  Determine Cartesian region location for
+*  this point
+*----
+            ISURC=NBCOR(1,IC)
+            REFLOC=TRKLSI(ISURC,IC)
+            IF(REFLOC .GT. REFLOA) THEN
+              WRITE(IOUT,9000) ISURA,ISURC,REFLOA,REFLOC
+              CALL XABORT(NAMSBR//
+     >': No Cartesian region found for this annulus')
+            ENDIF
+            IFIC=NBCOR(1,IC)+1
+            ILIC=NBSINT(IC)-NBCOR(2,IC)+1
+            DO JSURC=IFIC,ILIC
+              CURLOC=REFLOC+TRKLSI(JSURC,IC)
+              IF(REFLOA .LT. CURLOC) THEN
+                ISURC=JSURC
+                GO TO 125
+              ENDIF
+              REFLOC=CURLOC
+            ENDDO
+            CALL XABORT(NAMSBR//
+     >': Impossible to find Cartesian region containing annulus')
+ 125        CONTINUE
+            IFIC=ISURC
+          ELSE
+*----
+*  It is an outer Cartesian face,
+*  Determine Cartesian region location for
+*  this face
+*----
+            ISURC=NBCOR(1,IC)
+            REFLOC=TRKLSI(ISURC,IC)
+            IF(REFLOC .NE. REFLOA) CALL XABORT(NAMSBR//
+     >': No compatible Cartesian face found for this annulus')
+            IFIC=ISURC+1
+          ENDIF
+          IFIA=1
+          ILIA=NBCOR(1,IA)
+          DO ISURA=IFIA,ILIA
+            ISURG=ISURG+1
+            TRKLSI(ISURG,IG)=TRKLSI(ISURA,IA)
+            DO IDIR=1,5
+              IF(IDIR .EQ. KDIR .OR. IDIR .EQ. LDIR ) THEN
+                ISINT(IDIR,ISURG,IG)=ISINT(IDIR,ISURA,IA)
+              ELSE
+                ISINT(IDIR,ISURG,IG)=ISINT(IDIR,ISURC,IC)
+              ENDIF
+            ENDDO
+            ISINT(4,ISURG,IG)=ISINT(4,ISURA,IA)
+          ENDDO
+          NBCOR(1,IG)=ISURG
+          IFIA=NBCOR(1,IA)+1
+          IRADG=ISINT(4,IFIA,IA)
+          ILIA=NBSINT(IA)-NBCOR(2,IA)+1
+*----
+*  Scan Cartesian regions
+*----
+          ILIC=NBSINT(IC)-NBCOR(2,IC)+1
+          DO ISURC=IFIC,ILIC
+            CURLOC=REFLOC+TRKLSI(ISURC,IC)
+            ISURG=ISURG+1
+            DO JDIR=1,5
+              ISINT(JDIR,ISURG,IG)=ISINT(JDIR,ISURC,IC)
+            ENDDO
+            ISINT(4,ISURG,IG)=IRADG
+            THICK=CURLOC-REFLOG
+            TRKLSI(ISURG,IG)=THICK
+            ITIA=0
+            JSURC=ISURC
+*----
+*  Scan annular regions for intermediate mesh
+*----
+            DO ISURA=IFIA,ILIA
+              CURLOA=REFLOA+TRKLSI(ISURA,IA)
+              IF(CURLOA .GT. CURLOC) THEN
+                GO TO 145
+              ELSE IF(CURLOA .EQ. CURLOC) THEN
+                ITIA=ITIA+1
+                IRADG=ISINT(4,ISURA,IA)
+                DO JDIR=1,5
+                  IF(ISINT(JDIR,ISURA,IA) .NE. 0) THEN
+                    ISINT(JDIR,ISURG,IG)=ISINT(JDIR,ISURA,IA)
+                  ENDIF
+                ENDDO
+                ISINT(4,ISURG,IG)=IRADG
+*                IF(IPRINT .GT. 1000)
+*     >          write(6,7001) 'CURLOA .EQ. CURLOC',ISURG,
+*     >          (ISINT(JDIR,ISURG,IG),JDIR=1,5)
+                THICK=CURLOA-REFLOG
+                REFLOA=CURLOA
+                REFLOG=REFLOA
+                GO TO 145
+              ELSE
+                IF(ISURA .EQ. ILIA) THEN
+                  THICK=CURLOA-REFLOG
+                  TRKLSI(ISURG,IG)=THICK
+                  IRADG=0
+                  GO TO 155
+                ELSE
+                  ITIA=ITIA+1
+                  IRADG=ISINT(4,ISURA+1,IA)
+                  THICK=CURLOA-REFLOG
+                  DO JDIR=1,5
+                    ISINT(JDIR,ISURG+1,IG)=ISINT(JDIR,ISURG,IG)
+                  ENDDO
+                  ISINT(4,ISURG+1,IG)=IRADG
+                  TRKLSI(ISURG+1,IG)=TRKLSI(ISURG,IG)-THICK
+                  TRKLSI(ISURG,IG)=THICK
+                  ISURG=ISURG+1
+                ENDIF
+                REFLOA=CURLOA
+                REFLOG=REFLOA
+              ENDIF
+            ENDDO
+ 145        CONTINUE
+            IFIA=IFIA+ITIA
+            REFLOC=CURLOC
+            REFLOG=REFLOC
+          ENDDO
+ 155      CONTINUE
+*----
+*  Save final annular faces
+*----
+          IFIA=NBSINT(IA)-NBCOR(2,IA)+2
+          ILIA=NBSINT(IA)+1
+          ISURA=IFIA
+          DO IDIR=1,5
+            IF(ISINT(IDIR,ISURA,IA) .LT. 0) THEN
+              KDIR=IDIR
+            ELSE IF(ISINT(IDIR,ISURA,IA) .GT. 0) THEN
+              LDIR=IDIR
+            ENDIF
+          ENDDO
+          ISURC=JSURC
+          IFIA=NBSINT(IA)-NBCOR(2,IA)+2
+          ILIA=NBSINT(IA)+1
+          DO ISURA=IFIA,ILIA
+            ISURG=ISURG+1
+            TRKLSI(ISURG,IG)=TRKLSI(ISURA,IA)
+             DO IDIR=1,5
+              IF(IDIR .EQ. KDIR .OR. IDIR .EQ. LDIR ) THEN
+                ISINT(IDIR,ISURG,IG)=ISINT(IDIR,ISURA,IA)
+              ELSE
+                ISINT(IDIR,ISURG,IG)=ISINT(IDIR,ISURC,IC)
+              ENDIF
+            ENDDO
+            ISINT(4,ISURG,IG)=ISINT(4,ISURA,IA)
+          ENDDO
+          NBCOR(2,IG)=NBCOR(2,IA)
+        ELSE
+          IFIA=1
+          ILIA=NBCOR(1,IA)
+          DO ISURA=IFIA,ILIA
+            ISURG=ISURG+1
+            TRKLSI(ISURG,IG)=TRKLSI(ISURA,IA)
+            DO IDIR=1,5
+              ISINT(IDIR,ISURG,IG)=ISINT(IDIR,ISURA,IA)
+            ENDDO
+          ENDDO
+          NBCOR(1,IG)=NBCOR(1,IA)
+*----
+*  Scan annular regions for intermediate mesh
+*----
+          IFIA=NBCOR(1,IA)+1
+          ILIA=NBSINT(IA)-NBCOR(2,IA)+1
+          DO ISURA=IFIA,ILIA
+            ISURG=ISURG+1
+            CURLOA=REFLOA+TRKLSI(ISURA,IA)
+            DO JDIR=1,5
+              ISINT(JDIR,ISURG,IG)=ISINT(JDIR,ISURA,IA)
+            ENDDO
+            TRKLSI(ISURG,IG)=TRKLSI(ISURA,IA)
+            REFLOA=CURLOA
+          ENDDO
+          IFIA=NBSINT(IA)-NBCOR(2,IA)+2
+          ILIA=NBSINT(IA)+1
+          DO ISURA=IFIA,ILIA
+            ISURG=ISURG+1
+            TRKLSI(ISURG,IG)=TRKLSI(ISURA,IA)
+            DO IDIR=1,5
+              ISINT(IDIR,ISURG,IG)=ISINT(IDIR,ISURA,IA)
+            ENDDO
+          ENDDO
+          NBCOR(2,IG)=NBCOR(2,IA)
+          NBSINT(IG)=ISURG-1
+        ENDIF
+      ENDIF
+*----
+*  Remove track segment with vanishing distance
+*----
+      IT=3
+      NBSINT(IT)=ISURG-1
+      ISUR=NBCOR(1,IT)
+      DO JSUR=NBCOR(1,IT)+1,NBSINT(IT)-NBCOR(2,IT)+1
+        IF(TRKLSI(JSUR,IT) .GT. DCUTOF) THEN
+          ISUR=ISUR+1
+          TRKLSI(ISUR,IT)=TRKLSI(JSUR,IT)
+          DO IDIR=1,5
+            ISINT(IDIR,ISUR,IT)=ISINT(IDIR,JSUR,IT)
+          ENDDO
+        ELSE
+          IF(IPRINT .GT. 1000)
+     >    WRITE(IOUT,9001) JSUR,IT,ISUR
+        ENDIF
+      ENDDO
+      IF(ISUR .LT. NBSINT(IT)-NBCOR(2,IT)+1) THEN
+        DO JSUR=NBSINT(IT)-NBCOR(2,IT)+2,NBSINT(IT)+1
+          ISUR=ISUR+1
+          TRKLSI(ISUR,IT)=TRKLSI(JSUR,IT)
+          DO IDIR=1,5
+            ISINT(IDIR,ISUR,IT)=ISINT(IDIR,JSUR,IT)
+          ENDDO
+        ENDDO
+        NBSINT(IT)=ISUR-1
+      ENDIF
+      IF(IPRINT .GT. 1000) THEN
+        WRITE(IOUT,6016) NAMTYP(IT)
+        WRITE(IOUT,6011) 'Initial face'
+        DO JSUR=1,NBCOR(1,IT)
+          WRITE(IOUT,6010) TRKLSI(JSUR,IT),
+     >    (ISINT(IDIR,JSUR,IT),IDIR=1,5)
+        ENDDO
+        WRITE(IOUT,6011) 'Regions '
+        DO JSUR=NBCOR(1,IT)+1,NBSINT(IT)-NBCOR(2,IT)+1
+          WRITE(IOUT,6010) TRKLSI(JSUR,IT),
+     >    (ISINT(IDIR,JSUR,IT),IDIR=1,5)
+        ENDDO
+        WRITE(IOUT,6011) 'Final face'
+        DO JSUR=NBSINT(IT)-NBCOR(2,IT)+2,NBSINT(IT)+1
+          WRITE(IOUT,6010) TRKLSI(JSUR,IT),
+     >    (ISINT(IDIR,JSUR,IT),IDIR=1,5)
+        ENDDO
+        WRITE(IOUT,6001) NAMSBR
+      ENDIF
+      RETURN
+*----
+*  Output formats
+*----
+ 6000 FORMAT('(* Output from --',A6,'-- follows ')
+ 6001 FORMAT('   Output from --',A6,'-- completed *)')
+ 6010 FORMAT(1X,F25.16,5I10)
+ 6011 FORMAT(A20)
+ 6016 FORMAT(1X,A12)
+* 7000 FORMAT('Verify ',A6,I10,F25.16,6I10)
+* 7001 FORMAT('Verify ',A18,I10,6I10)
+ 9000 FORMAT(' Problem with surface ',2I10,2F20.10)
+ 9001 FORMAT(' Warning : region with vanishing distance found ',
+     >       3I10)
+ 9002 FORMAT('Final Cartesian and annular faces are incoherent',
+     >       2(2X,I10),2(2X,F20.10))
+      END