Initial commit from Polytechnique Montreal

1 files changed, 784 insertions, 0 deletions

diff --git a/Dragon/src/NXTLCY.f b/Dragon/src/NXTLCY.f
new file mode 100644
index 0000000..a9f6018
--- /dev/null
+++ b/Dragon/src/NXTLCY.f
@@ -0,0 +1,784 @@
+*DECK NXTLCY
+      FUNCTION NXTLCY(IPRINT,ITST  ,NDIM  ,MXMESH,LINMAX,
+     >                MESH  ,ORITRK,DIRTRK,DCMESH,IDIRC ,
+     >                NBCOR ,NBSINT,ISINT ,TRKLSI)
+*
+*----------
+*
+*Purpose:
+* To track an annular 2-D or 3-D geometry
+* using the NXT tracking procedure.
+*
+*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.
+* ITST    type of tracking, where:
+*         =-1   only the exact geometry
+*               is considered taking into account the
+*               submesh in each direction;
+*         = 0   only the global geometry
+*               is considered without taking into account the
+*               submesh in each direction;
+*         = 1   both the global
+*               geometry (as a first step) and the exact geometry
+*               are considered taking into account the
+*               submesh in each direction.
+* NDIM    dimension of problem.
+* MXMESH  maximum number of spatial subdivision in
+*         $R$ and $X$, $Y$ or $Z$.
+* LINMAX  maximum number of segments in a track.
+* MESH    effective number of spatial subdivision in $X$
+*         $Y$, $Z$ and $R$.
+* ORITRK  a point on the track (origin).
+* DIRTRK  the track direction (director cosines).
+* DCMESH  spatial description of the cylinder.
+* IDIRC   the direction of the cylinder axis.
+*
+*Parameters: output
+* NXTLCY  number of surfaces intersected.
+* NBCOR   number of corner found for each external faces.
+* NBSINT  number of surface crossed by track.
+* ISINT   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,ITST,NDIM,MXMESH,LINMAX
+      INTEGER          MESH(5)
+      DOUBLE PRECISION ORITRK(NDIM),
+     >                 DIRTRK(NDIM),
+     >                 DCMESH(-1:MXMESH,5)
+      INTEGER          IDIRC
+      INTEGER          NBCOR(2),NBSINT
+      INTEGER          ISINT(0:5,LINMAX)
+      DOUBLE PRECISION TRKLSI(LINMAX)
+      INTEGER          NXTLCY
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='NXTLCY')
+      INTEGER          MXDIM
+      PARAMETER       (MXDIM=3)
+      DOUBLE PRECISION DZERO,DONE,DTWO
+      PARAMETER       (DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0)
+      DOUBLE PRECISION DCUTOL
+      PARAMETER       (DCUTOL=1.0D-11)
+*----
+*  Local variables
+*----
+      INTEGER          IDIR,INDIR,INEXT,IFRST,ILAST,IFACE,
+     >                 IDGR,IDG1,IDG2,IDGP,NR,NP
+      INTEGER          JFACE,JSUR,KNEXT,KSFRST,KSLAST,
+     >                 KVFRST,KVLAST,KFACE,KSUR,KSUB
+      INTEGER          IC1,IDIRS,IDIRF,IDIRB,IRFIN,IDIRFS
+      DOUBLE PRECISION DP1,DP2,TRKK
+      DOUBLE PRECISION TRKDIS,TRKOLD,TRKTMP
+      DOUBLE PRECISION TORI(MXDIM),ROTANG(MXDIM)
+      DOUBLE PRECISION XYCEN(2),RADIUS,CYLLIM(2),PROJ,PROJN,
+     >                 YLOC,XORI,XLOC(2),CYLINT,TCINTD,XXINTD
+      DOUBLE PRECISION DCUTOF
+*----
+*  Data
+*----
+      CHARACTER        CDIR(4)*1
+      SAVE             CDIR
+      DATA             CDIR /'X','Y','Z','R'/
+*----
+*  Verify ITST option and reset to default value if invalid
+*----
+      IF(ITST .LT. -1 .OR. ITST .GT. 1) THEN
+*----
+*  Reset ITST=1 (complete analysis) if the value of ITST is invalid
+*----
+        ITST=1
+      ENDIF
+*----
+*  Initialise output vectors
+*----
+      KSFRST=1
+      KSLAST=1
+      KVFRST=1
+      KVLAST=1
+      KNEXT=1
+      KSUB=1
+      TRKOLD=DZERO
+      NBCOR(1)=0
+      NBCOR(2)=0
+      NBSINT=0
+      ISINT(0:5,:LINMAX)=0
+      TRKLSI(:LINMAX)=DZERO
+*----
+*  IDG1 is first direction of plane perpendicular
+*       to main direction ($Y, $Z$ or $X$).
+*  IDG2 is second direction of plane perpendicular
+*       to main direction ($Z$, $X$ or $Y$).
+*  IDGP is main cylinder direction ($X$, $Y$ or $Z$)
+*       for 2D case IDGP=3
+*----
+      IDGR=4
+      IDGP=IDIRC
+      IDG1=MOD(IDGP,3)+1
+      IDG2=MOD(IDGP+1,3)+1
+      NR=MESH(4)
+      NP=MESH(IDGP)
+*----
+*  Print header if required
+*----
+      IF(IPRINT .GT. 1000) THEN
+        WRITE(IOUT,6000) NAMSBR
+        WRITE(IOUT,6011) 'radial'//CDIR(IDG1)//CDIR(IDG2)//
+     >  '={          '
+        WRITE(IOUT,6012)  (DCMESH(IC1,IDGR),IC1=1,NR)
+        WRITE(IOUT,6013)
+        WRITE(IOUT,6011) 'center'//CDIR(IDG1)//CDIR(IDG2)//
+     >  '={          '
+        WRITE(IOUT,6012)  DCMESH(-1,IDG1),DCMESH(-1,IDG2)
+        WRITE(IOUT,6013)
+        IF(NDIM .EQ. 3) THEN
+          WRITE(IOUT,6011) 'axial'//CDIR(IDGP)//'={            '
+          WRITE(IOUT,6012)  (DCMESH(IC1,IDGP),IC1=0,NP)
+          WRITE(IOUT,6013)
+        ENDIF
+        WRITE(IOUT,6011) 'trackorigin={       '
+        WRITE(IOUT,6012) ORITRK
+        WRITE(IOUT,6013)
+        WRITE(IOUT,6011) 'trackdirection={    '
+        WRITE(IOUT,6012) DIRTRK
+        WRITE(IOUT,6013)
+      ENDIF
+*----
+*  Select planes order in direction IDIR
+*----
+      IFRST=NR+1
+      INEXT=-1
+      ILAST=2
+*----
+*  Select planes order in direction KDIR (if required)
+*----
+      DCUTOF=DCUTOL*DCUTOL
+      IF(NDIM .EQ. 3) THEN
+        DCUTOF=DCUTOF*DCUTOL
+        IF(DIRTRK(IDGP) .LT. 0) THEN
+          KNEXT=-1
+          KSFRST=NP+1
+          KSLAST=1
+          KVFRST=NP+1
+          KVLAST=2
+          KSUB=-1
+        ELSE
+          KNEXT=1
+          KSFRST=1
+          KSLAST=NP+1
+          KVFRST=1
+          KVLAST=NP
+          KSUB=0
+        ENDIF
+      ENDIF
+*----
+*  Projection on 2-D plane
+*  and translation to origin of circle
+*----
+      XYCEN(1)=DCMESH(-1,IDG1)
+      XYCEN(2)=DCMESH(-1,IDG2)
+      TORI(1)=ORITRK(IDG1)
+      TORI(2)=ORITRK(IDG2)
+      IF(NDIM .EQ. 3) THEN
+        TORI(3)=ORITRK(IDGP)
+        ROTANG(3)=DIRTRK(IDGP)
+      ELSE
+        ROTANG(3)=DZERO
+        TORI(3)=DZERO
+      ENDIF
+      PROJ=DONE/SQRT(DONE-ROTANG(3)**2)
+      ROTANG(1)=DIRTRK(IDG1)*PROJ
+      ROTANG(2)=DIRTRK(IDG2)*PROJ
+      PROJN=DONE/SQRT(ROTANG(1)*ROTANG(1)+ROTANG(2)*ROTANG(2))
+      ROTANG(1)=ROTANG(1)*PROJN
+      ROTANG(2)=ROTANG(2)*PROJN
+*----
+*  Translate (x,y) coordinates to circle center and
+*  rotate y coordinate in such a way that tracking
+*  line parallel to the x-axis
+*----
+      XORI= (TORI(1)-XYCEN(1))*ROTANG(1)
+     >     +(TORI(2)-XYCEN(2))*ROTANG(2)
+      YLOC=-(TORI(1)-XYCEN(1))*ROTANG(2)
+     >     +(TORI(2)-XYCEN(2))*ROTANG(1)
+*      IF(IPRINT .GT. 1000) THEN
+*        write(6,*) 'XYCEN =',XYCEN
+*        write(6,*) 'TORI  =',TORI
+*        write(6,*) 'ROTANG=',ROTANG
+*        write(6,*) 'RotOri=',XORI,YLOC
+*      ENDIF
+*----
+*  For option ITST=0, 1, consider global geometry
+*  and test if line crosses this geometry.
+*----
+      IF(ITST .EQ. 0 .OR. ITST .EQ. 1) THEN
+        NBSINT=0
+        RADIUS=DCMESH(NR,IDGR)
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'YLOC, RADIUS=',YLOC,RADIUS
+        IF(NDIM .EQ. 3) THEN
+          CYLLIM(1)=DCMESH(0,IDGP)
+          CYLLIM(2)=DCMESH(NP,IDGP)
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*)'CYLLIM ',CYLLIM
+        ENDIF
+*----
+*  Test if coordinate is inside circle.
+*----
+        IF(ABS(YLOC) .LT. RADIUS) THEN
+*----
+*  Line crosses circle (infinite cylinder)
+*----
+          XLOC(2)=SQRT(RADIUS**2-YLOC**2)
+          XLOC(1)=-XLOC(2)
+          IF(NDIM .EQ. 3) THEN
+*----
+*  Find intersection points for finite cylinders in 3D
+*----
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'XLOC ',XLOC(1),XLOC(2)
+              DO 100 IFACE=1,2
+              TCINTD=YLOC*ROTANG(1)+
+     >               XLOC(IFACE)*ROTANG(2)+XYCEN(2)-TORI(2)
+              XXINTD=XLOC(IFACE)*ROTANG(1)-
+     >               YLOC*ROTANG(2)+XYCEN(1)-TORI(1)
+              TRKTMP=(TCINTD*PROJ)/ROTANG(2)
+              TRKDIS=SQRT(XXINTD*XXINTD+TCINTD*TCINTD)*PROJ
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'Changedir=',TRKDIS,TCINTD,ROTANG(2)
+              IF(ROTANG(2)*TCINTD .LT. DZERO) TRKDIS=-TRKDIS
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'TRKDIS=',TRKDIS,TRKTMP,ROTANG(2)
+              CYLINT=TRKDIS*ROTANG(3)+TORI(3)
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'CYLINT= ',CYLINT
+              IF(CYLINT .GE. CYLLIM(1) .AND. CYLINT .LE. CYLLIM(2)) THEN
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'iface r ',IFACE
+                IF(NBSINT .EQ. 0) THEN
+                  NBSINT=1
+                  TRKOLD=TRKDIS
+                ELSE IF(TRKOLD .NE. TRKDIS) THEN
+                  NBSINT=2
+                  GO TO 105
+                ENDIF
+              ENDIF
+ 100        CONTINUE
+          ELSE
+            NBSINT=2
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'Intersected in 2D'
+            GO TO 105
+          ENDIF
+        ENDIF
+*----
+*  For 3-D geometry look for intersection with bottom
+*  and top of cylinder
+*----
+        IF(NDIM .EQ. 3) THEN
+          DO 102 IFACE=1,2
+            TRKDIS=(CYLLIM(IFACE)-TORI(3))/ROTANG(3)
+            XLOC(1)=TORI(1)+TRKDIS*ROTANG(1)/PROJ-XYCEN(1)
+            XLOC(2)=TORI(2)+TRKDIS*ROTANG(2)/PROJ-XYCEN(2)
+            TCINTD=SQRT(XLOC(1)**2+XLOC(2)**2)
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'dir 3 int ', XLOC(1),XLOC(2),TRKDIS
+            IF(TCINTD .LT. RADIUS) THEN
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'iface z ',IFACE
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) TORI(1)+TRKDIS*ROTANG(1)/PROJ,
+*     >                 TORI(2)+TRKDIS*ROTANG(2)/PROJ,
+*     >                 TORI(3)+TRKDIS*ROTANG(3)
+              IF(NBSINT .EQ. 0) THEN
+                NBSINT=1
+                TRKOLD=TRKDIS
+              ELSE IF(TRKOLD .NE. TRKDIS) THEN
+                NBSINT=2
+                GO TO 105
+              ENDIF
+            ENDIF
+ 102      CONTINUE
+        ENDIF
+ 105    CONTINUE
+      ELSE
+*----
+*  If ITST=-1, assume that 2 surface are intersected
+*  and continue
+*----
+        NBSINT=2
+      ENDIF
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'Nb sur=',NBSINT
+      NXTLCY=NBSINT
+*----
+*  If line does not intersect any surfaces, exit
+*----
+      IF(NXTLCY .EQ. 0) THEN
+        IF(IPRINT .GT. 1000) THEN
+          WRITE(IOUT,6002) NAMSBR
+        ENDIF
+        RETURN
+      ENDIF
+*----
+*  If line does not intersect 2 surfaces, abort
+*----
+      IF(NBSINT .NE. 2) CALL XABORT(NAMSBR//
+     >': Invalid line since it did not intersect 0 or 2 surfaces')
+*----
+*  If ITST=0, return number of surfaces intersected by line
+*----
+      IF(ITST .EQ. 0) RETURN
+*----
+*  If ITST=1 or -1
+*  track geometry with submesh
+*----
+      NBSINT=0
+*----
+*  Loop over radial regions from maximal radius
+*  from outer to inner
+*----
+      DO 110 IFACE=IFRST,ILAST,INEXT
+        RADIUS=DCMESH(IFACE-1,IDGR)
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'YLOC, RADIUS=',YLOC,RADIUS
+        IF(ABS(YLOC) .LE. RADIUS) THEN
+          XLOC(2)=SQRT(RADIUS**2-YLOC**2)
+          XLOC(1)=-XLOC(2)
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'XLOC',XLOC
+*----
+*  Scan over two possible cylinder intersections
+*----
+          INDIR=1
+          DO 111 JFACE=1,2
+            TCINTD=YLOC*ROTANG(1)+XLOC(JFACE)*ROTANG(2)+XYCEN(2)-TORI(2)
+            XXINTD=XLOC(JFACE)*ROTANG(1)-YLOC*ROTANG(2)+XYCEN(1)-TORI(1)
+            TRKDIS=SQRT(XXINTD*XXINTD+TCINTD*TCINTD)*PROJ
+*      IF(IPRINT .GT. 1000)
+*     >  TRKTMP=(TCINTD*PROJ)/ROTANG(2)
+            IF(ROTANG(2)*TCINTD .LT. -DCUTOF) TRKDIS=-TRKDIS
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'Distances',XXINTD,TCINTD,ROTANG(2),TRKDIS,TRKTMP
+            IF(NDIM .EQ. 3) THEN
+*----
+*  Scan over planes in axial cylinder direction
+*----
+              TRKK=DBLE(KNEXT)*(TRKDIS*ROTANG(3)+TORI(3))
+              DO 112 KFACE=KVFRST,KVLAST,KNEXT
+                DP1=DBLE(KNEXT)*DCMESH(KFACE-1,IDGP)
+                DP2=DBLE(KNEXT)*DCMESH(KFACE-1+KNEXT,IDGP)
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'DP1,DP2 ',DP1,DP2,TRKK
+                IF(TRKK .GE. DP1 .AND.
+     >             TRKK .LE. DP2 ) THEN
+                  NBSINT=NBSINT+1
+                  IF(IPRINT .GT. 1000) THEN
+                    WRITE(IOUT,6011) 'Intersection at     '
+                      WRITE(IOUT,6014) IFACE,JFACE,KFACE
+                      WRITE(IOUT,6015) TRKDIS
+                  ENDIF
+                  IF(NBSINT .EQ. 1) THEN
+                    ISINT(0,NBSINT)=(IDG1+3)*INDIR
+                    ISINT(IDGR,NBSINT)=IFACE-1
+                    ISINT(IDGP,NBSINT)=KFACE+KSUB
+                    TRKLSI(NBSINT)=TRKDIS
+                  ELSE
+*----
+*  Scan previous distances and reorder if necessary
+*----
+                    DO 120 JSUR=NBSINT-1,1,-1
+                      IF(TRKDIS .GE. TRKLSI(JSUR)) THEN
+*----
+*  For more than one intersection at a given distance:
+*  a) if this is an external surface, store the distance
+*     for corner duplication
+*  b) if this is an internal surface select adequate next region
+*     to cross and reset NBSINT=NBSINT-1
+*----
+                        IF(TRKDIS .EQ. TRKLSI(JSUR) .AND.
+     >                     IFACE .NE. IFRST             ) THEN
+                          ISINT(IDGR,JSUR)=IFACE-1
+                          ISINT(IDGP,JSUR)=KFACE+KSUB
+                          NBSINT=NBSINT-1
+                          GO TO 125
+                        ENDIF
+                        DO 121 KSUR=NBSINT,JSUR+2,-1
+                          ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                          ISINT(IDGR,KSUR)=ISINT(IDGR,KSUR-1)
+                          ISINT(IDGP,KSUR)=ISINT(IDGP,KSUR-1)
+                          TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 121                    CONTINUE
+                        ISINT(0,JSUR+1)=(IDG1+3)*INDIR
+                        ISINT(IDGR,JSUR+1)=IFACE-1
+                        ISINT(IDGP,JSUR+1)=KFACE+KSUB
+                        TRKLSI(JSUR+1)=TRKDIS
+                        GO TO 125
+                      ENDIF
+ 120                CONTINUE
+                    DO 122 KSUR=NBSINT,2,-1
+                      ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                      ISINT(IDGR,KSUR)=ISINT(IDGR,KSUR-1)
+                      ISINT(IDGP,KSUR)=ISINT(IDGP,KSUR-1)
+                      TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 122                CONTINUE
+                    ISINT(0,1)=(IDG1+3)*INDIR
+                    ISINT(IDGR,1)=IFACE-1
+                    ISINT(IDGP,1)=KFACE+KSUB
+                    TRKLSI(1)=TRKDIS
+ 125                CONTINUE
+                  ENDIF
+                ENDIF
+ 112          CONTINUE
+            ELSE
+*----
+*  For 2-D, intersection with a circle
+*----
+              NBSINT=NBSINT+1
+              IF(IPRINT .GT. 1000) THEN
+                WRITE(IOUT,6011) 'Intersection at     '
+                WRITE(IOUT,6014) IFACE,JFACE
+                WRITE(IOUT,6015) TRKDIS
+              ENDIF
+              IF(NBSINT .EQ. 1) THEN
+                ISINT(0,NBSINT)=(IDG1+3)*INDIR
+                ISINT(IDGR,NBSINT)=IFACE-1
+                TRKLSI(NBSINT)=TRKDIS
+              ELSE
+*----
+*  Scan previous distances and reorder if necessary
+*----
+                DO 130 JSUR=NBSINT-1,1,-1
+                  IF(TRKDIS .GE. TRKLSI(JSUR)) THEN
+*----
+*  For more than one intersection at a given distance:
+*  a) if this is an external surface, store the distance
+*     for corner duplication
+*  b) if this is an internal surface select adequate next region
+*     to cross and reset NBSINT=NBSINT-1
+*----
+                    IF(TRKDIS .EQ. TRKLSI(JSUR) .AND.
+     >                 IFACE .NE. IFRST              ) THEN
+                      ISINT(IDGR,JSUR)=IFACE-1
+                      NBSINT=NBSINT-1
+                      GO TO 135
+                    ENDIF
+                    DO 131 KSUR=NBSINT,JSUR+2,-1
+                      ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                      ISINT(IDGR,KSUR)=ISINT(IDGR,KSUR-1)
+                      TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 131                CONTINUE
+                    ISINT(0,JSUR+1)=(IDG1+3)*INDIR
+                    ISINT(IDGR,JSUR+1)=IFACE-1
+                    TRKLSI(JSUR+1)=TRKDIS
+                    GO TO 135
+                  ENDIF
+ 130            CONTINUE
+                DO 132 KSUR=NBSINT,2,-1
+                  ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                  ISINT(IDGR,KSUR)=ISINT(IDGR,KSUR-1)
+                  TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 132            CONTINUE
+                ISINT(0,1)=(IDG1+3)*INDIR
+                ISINT(IDGR,1)=IFACE-1
+                TRKLSI(1)=TRKDIS
+ 135            CONTINUE
+              ENDIF
+            ENDIF
+            INDIR=-1
+ 111      CONTINUE
+        ENDIF
+ 110  CONTINUE
+      IF(NDIM .EQ. 3) THEN
+*----
+*  Loop over axial regions from Top to bottom or
+*  vice versa
+*----
+        DO 140 KFACE=KSFRST,KSLAST,KNEXT
+          TRKDIS=(DCMESH(KFACE-1,IDGP)-TORI(3))/ROTANG(3)
+          DP1=TORI(1)+TRKDIS*ROTANG(1)/PROJ-XYCEN(1)
+          DP2=TORI(2)+TRKDIS*ROTANG(2)/PROJ-XYCEN(2)
+          TCINTD=SQRT(DP1**2+DP2**2)
+*----
+*  Loop over radial regions from inner to outer
+*----
+          DO 141 IFACE=ILAST,IFRST,-INEXT
+            RADIUS=DCMESH(IFACE-1,IDGR)
+*      IF(IPRINT .GT. 1000)
+*     >  write(6,*) 'zz= ',DP1,DP2,TCINTD,RADIUS
+            IF(TCINTD .LE. RADIUS) THEN
+              IF(IPRINT .GT. 1000) THEN
+                WRITE(IOUT,6011) 'Intersection at     '
+                WRITE(IOUT,6014) KFACE,IFACE
+                WRITE(IOUT,6015) TRKDIS
+              ENDIF
+              NBSINT=NBSINT+1
+              IF(NBSINT .EQ. 1) THEN
+                ISINT(0,NBSINT)=IDGP
+                ISINT(IDGR,NBSINT)=IFACE-1
+                ISINT(IDGP,NBSINT)=KFACE
+                TRKLSI(NBSINT)=TRKDIS
+              ELSE
+*----
+*  Scan previous distances and reorder if necessary
+*----
+                DO 150 JSUR=NBSINT-1,1,-1
+                  IF(TRKDIS .GE. TRKLSI(JSUR)) THEN
+*----
+*  For more than one intersection at a given distance:
+*  a) if this is an external surface, store the distance
+*     for corner duplication
+*  b) if this is an internal surface select adequate next region
+*     to cross and reset NBSINT=NBSINT-1
+*----
+                    IF(TRKDIS .EQ. TRKLSI(JSUR) .AND.
+     >                 KFACE .NE. KSFRST         .AND.
+     >                 KFACE .NE. KSLAST              ) THEN
+                      ISINT(IDGR,JSUR)=IFACE-1
+                      ISINT(IDGP,JSUR)=KFACE
+                      NBSINT=NBSINT-1
+                      GO TO 155
+                    ENDIF
+                    DO 151 KSUR=NBSINT,JSUR+2,-1
+                      ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                      ISINT(IDGR,KSUR)=ISINT(IDGR,KSUR-1)
+                      ISINT(IDGP,KSUR)=ISINT(IDGP,KSUR-1)
+                      TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 151                CONTINUE
+                    ISINT(0,JSUR+1)=IDGP
+                    ISINT(IDGR,JSUR+1)=IFACE-1
+                    ISINT(IDGP,JSUR+1)=KFACE
+                    TRKLSI(JSUR+1)=TRKDIS
+                    GO TO 155
+                  ENDIF
+ 150            CONTINUE
+                DO 152 KSUR=NBSINT,2,-1
+                  ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                  ISINT(IDGR,KSUR)=ISINT(IDGR,KSUR-1)
+                  ISINT(IDGP,KSUR)=ISINT(IDGP,KSUR-1)
+                  TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 152            CONTINUE
+                ISINT(0,1)=IDGP
+                ISINT(IDGR,1)=IFACE-1
+                ISINT(IDGP,1)=KFACE
+                TRKLSI(1)=TRKDIS
+ 155            CONTINUE
+              ENDIF
+              GO TO 145
+            ENDIF
+ 141      CONTINUE
+ 145      CONTINUE
+ 140    CONTINUE
+      ENDIF
+      IF(NBSINT .GE. 2) THEN
+*----
+*  Find number of corners for initial surface intersections
+*----
+        TRKDIS=TRKLSI(1)
+        NBCOR(1)=1
+        DO 160 JSUR=2,NBSINT
+          IF(TRKLSI(JSUR) .GT. TRKDIS) GO TO 165
+          NBCOR(1)=NBCOR(1)+1
+ 160    CONTINUE
+ 165    CONTINUE
+*----
+*  Find number of corners for final surface intersections
+*----
+        TRKDIS=TRKLSI(NBSINT)
+        NBCOR(2)=1
+        DO 170 JSUR=NBSINT-1,1,-1
+          IF(TRKLSI(JSUR) .LT. TRKDIS) GO TO 175
+          NBCOR(2)=NBCOR(2)+1
+ 170    CONTINUE
+ 175    CONTINUE
+      ENDIF
+*----
+*  Print intersection points
+*----
+      IF(IPRINT .GT. 1000) THEN
+        WRITE(IOUT,6011) 'Initial planes      '
+        DO 200 JSUR=1,NBCOR(1)
+          WRITE(IOUT,6010) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5)
+ 200    CONTINUE
+        IF(NBCOR(1)+1 .LE. NBSINT-NBCOR(2)) THEN
+          WRITE(IOUT,6011) 'Intermediate planes '
+          DO 201 JSUR=NBCOR(1)+1,NBSINT-NBCOR(2)
+            WRITE(IOUT,6010) TRKLSI(JSUR),
+     >      (ISINT(IDIR,JSUR),IDIR=0,5)
+ 201      CONTINUE
+        ENDIF
+        WRITE(IOUT,6011) 'Final planes        '
+        DO 202 JSUR=NBSINT-NBCOR(2)+1,NBSINT
+          WRITE(IOUT,6010) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5)
+ 202    CONTINUE
+      ENDIF
+*----
+*  Identify final faces
+*----
+      DO JSUR=NBSINT,NBSINT-NBCOR(2)+1,-1
+        DO IDIRS=0,5
+          ISINT(IDIRS,JSUR+1)=ISINT(IDIRS,JSUR)
+        ENDDO
+        TRKLSI(JSUR+1)=TRKLSI(JSUR)
+        IDIRS=ABS(ISINT(0,JSUR+1))
+        IF(IDIRS .GE. 4) THEN
+          IDIRF=ISINT(4,JSUR+1)
+          IF(IDIRF .NE. MESH(4)) THEN
+            WRITE(IOUT,9001) NAMSBR,'final   ',
+     >                    JSUR,4,IDIRF,MESH(4)
+            CALL XABORT(NAMSBR//
+     >      ': Invalid '//CDIR(4)//' directed surface')
+          ENDIF
+          ISINT(4,JSUR+1)=-2
+        ELSE
+          IDIRF=ISINT(IDIRS,JSUR+1)
+          IF(IDIRF .EQ. 1) THEN
+            ISINT(IDIRS,JSUR+1)=-1
+          ELSE IF(IDIRF .EQ. MESH(IDIRS)+1) THEN
+            ISINT(IDIRS,JSUR+1)=-2
+          ELSE
+            WRITE(IOUT,9001) NAMSBR,'final   ',
+     >                    JSUR,IDIRS,IDIRF,MESH(IDIRS)
+            CALL XABORT(NAMSBR//': Invalid '//CDIR(IDIRS)//
+     >      ' directed surface')
+          ENDIF
+        ENDIF
+      ENDDO
+*----
+*  Identify regions
+*----
+      DO JSUR=NBSINT-NBCOR(2)+1,NBCOR(1)+1,-1
+        IDIRFS=ISINT(0,JSUR)
+        IF(IDIRFS .GE. 4) THEN
+          IDIRFS=1
+        ELSE
+          IDIRFS=0
+        ENDIF
+        IDIRB=ABS(ISINT(0,JSUR-1))
+        IF(IDIRB .GE. 4) IDIRB=4
+        IDIRF=ABS(ISINT(0,JSUR))
+        IF(IDIRF .GE. 4) IDIRF=4
+        IF(IDIRB .EQ. IDIRF) THEN
+          IF(IDIRB .EQ. 4) THEN
+            IRFIN=MAX(ISINT(4,JSUR-1),ISINT(4,JSUR))
+          ELSE
+            IRFIN=MIN(ISINT(IDIRF,JSUR-1),ISINT(IDIRF,JSUR))
+          ENDIF
+        ELSE
+          IF(IDIRF .EQ. 4) THEN
+            IF(IDIRFS .EQ. 1) THEN
+              IRFIN=ISINT(4,JSUR-1)
+            ELSE
+              IRFIN=ISINT(4,JSUR)
+            ENDIF
+          ELSE
+            IRFIN=ISINT(IDIRF,JSUR-1)
+          ENDIF
+        ENDIF
+        IF(IPRINT .GT. 1000) THEN
+          WRITE(IOUT,*) 'JSUR =',JSUR,ISINT(0,JSUR-1),ISINT(0,JSUR),
+     >                  IDIRB,IDIRF,IRFIN
+        ENDIF
+        TRKLSI(JSUR)=TRKLSI(JSUR)-TRKLSI(JSUR-1)
+        DO IDIRS=1,5
+          IF(IDIRS .EQ. IDIRF) THEN
+            ISINT(IDIRS,JSUR)=IRFIN
+          ELSE
+            ISINT(IDIRS,JSUR)=ISINT(IDIRS,JSUR)
+          ENDIF
+        ENDDO
+      ENDDO
+*----
+*  Identify initial face
+*----
+      DO JSUR=1,NBCOR(1)
+        IDIRS=ABS(ISINT(0,JSUR))
+        IF(IDIRS .GE. 4) THEN
+          IDIRB=ISINT(4,JSUR)
+          IF(IDIRB .NE. MESH(4)) THEN
+            WRITE(IOUT,9001) NAMSBR,'initial ',
+     >                    JSUR,4,IDIRB,MESH(4)
+            CALL XABORT(NAMSBR//
+     >      ': Invalid '//CDIR(4)//' directed surface')
+          ENDIF
+          ISINT(4,JSUR)=-2
+        ELSE
+          IDIRB=ISINT(IDIRS,JSUR)
+          IF(IDIRB .EQ. 1) THEN
+            ISINT(IDIRS,JSUR)=-1
+          ELSE IF(IDIRB .EQ. MESH(IDIRS)+1) THEN
+            ISINT(IDIRS,JSUR)=-2
+          ELSE
+            WRITE(IOUT,9001) NAMSBR,'initial ',
+     >                    JSUR,IDIRS,IDIRB,MESH(IDIRS)
+            CALL XABORT(NAMSBR//': Invalid '//CDIR(IDIRS)//
+     >      ' directed surface')
+          ENDIF
+        ENDIF
+      ENDDO
+*----
+*  Print final track information
+*----
+      IF(IPRINT .GT. 1000) THEN
+        WRITE(IOUT,6011) 'Initial face        '
+        DO JSUR=1,NBCOR(1)
+          WRITE(IOUT,6010) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5)
+        ENDDO
+        WRITE(IOUT,6011) 'Regions             '
+        DO JSUR=NBCOR(1)+1,NBSINT-NBCOR(2)+1
+          WRITE(IOUT,6010) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5)
+        ENDDO
+        WRITE(IOUT,6011) 'Final face          '
+        DO JSUR=NBSINT-NBCOR(2)+2,NBSINT+1
+          WRITE(IOUT,6010) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5)
+        ENDDO
+        WRITE(IOUT,6001) NAMSBR
+      ENDIF
+      RETURN
+*----
+*  Output formats
+*----
+ 6000 FORMAT('(* Output from --',A6,'-- follows ')
+ 6001 FORMAT('   Output from --',A6,'-- completed *)')
+ 6002 FORMAT('   No region crossed'/
+     >       '   Output from --',A6,'-- completed *)')
+ 6010 FORMAT(1X,F25.16,6I10)
+ 6011 FORMAT(A20)
+ 6012 FORMAT(6(1X,F25.16,:,','))
+ 6013 FORMAT('};')
+ 6014 FORMAT(6I10)
+ 6015 FORMAT(6(1X,F25.16))
+*----
+*  Warning formats
+*----
+ 9001 FORMAT(1X,'Error in --',A6,' --'/
+     >       1X,'invalid ',A8,' surface identification ',4I10)
+      END