Initial commit from Polytechnique Montreal

1 files changed, 809 insertions, 0 deletions

diff --git a/Dragon/src/NXTLCA.f b/Dragon/src/NXTLCA.f
new file mode 100644
index 0000000..9469b2c
--- /dev/null
+++ b/Dragon/src/NXTLCA.f
@@ -0,0 +1,809 @@
+*DECK NXTLCA
+      FUNCTION NXTLCA(IPRINT,ITST  ,NDIM  ,MXMESH,LINMAX,
+     >                MESH  ,ORITRK,DIRTRK,DCMESH,
+     >                NBCOR ,NBSINT,ISINT ,TRKLSI)
+*
+*----------
+*
+*Purpose:
+* To track a Cartesian 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
+*         $X$, $Y$ or $Z$.
+* LINMAX  maximum number of segments in a track.
+* MESH    effective number of spatial subdivision in
+*         each direction ($X$, $Y$ and $Z$).
+* ORITRK  a point on the track (origin).
+* DIRTRK  the track direction (director cosines).
+* DCMESH  spatial description of the parallepiped.
+*
+*Parameters: output
+* NXTLCA  number of surfaces intersected.
+* 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,ITST,NDIM,MXMESH,LINMAX
+      INTEGER          MESH(NDIM)
+      DOUBLE PRECISION ORITRK(NDIM),
+     >                 DIRTRK(NDIM),
+     >                 DCMESH(-1:MXMESH,5)
+      INTEGER          NBCOR(2),NBSINT
+      INTEGER          ISINT(0:5,LINMAX)
+      DOUBLE PRECISION TRKLSI(LINMAX)
+      INTEGER          NXTLCA
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='NXTLCA')
+      DOUBLE PRECISION DCUTOF,DZERO,DONE,DTWO
+      PARAMETER       (DCUTOF=1.0D-8,DZERO=0.0D0,DONE=1.0D0,DTWO=2.0D0)
+*----
+*  Local variables
+*----
+      INTEGER          IDIR,INEXT,IFRST,ILAST,IFACE,IOFF,NBTINT
+      INTEGER          JDIR,JNEXT,JFRST,JLAST,JFACE,JOFF,JSUR
+      INTEGER          KDIR,KNEXT,KFRST,KLAST,KFACE,KOFF,KSUR
+      INTEGER          IC1,IDIRS,IDIRF,IDIRB,IRFIN(5)
+      INTEGER          ISDIR,JSDIR,KSDIR
+      DOUBLE PRECISION DP1,DP2,TRKK,TRKJ
+      DOUBLE PRECISION TRKDIS,TRKOLD,DELTKD
+*----
+*  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
+      KFRST=1
+      KLAST=1
+      KNEXT=1
+      TRKK=DZERO
+      TRKOLD=DZERO
+      KOFF=1
+*----
+*  Initialise output vectors
+*----
+      NBCOR(1)=0
+      NBCOR(2)=0
+      NBSINT=0
+      ISINT(0:5,:LINMAX)=0
+      TRKLSI(:LINMAX)=DZERO
+*----
+*  Print header if required
+*----
+      IF(IPRINT .GT. 1000) THEN
+        WRITE(IOUT,6000) NAMSBR
+        WRITE(IOUT,6011) 'meshx={             '
+        WRITE(IOUT,6012) (DCMESH(IC1,1),IC1=0,MESH(1))
+        WRITE(IOUT,6013)
+        WRITE(IOUT,6011) 'meshy={             '
+        WRITE(IOUT,6012) (DCMESH(IC1,2),IC1=0,MESH(2))
+        WRITE(IOUT,6013)
+        IF(NDIM .EQ. 3) THEN
+          WRITE(IOUT,6011) 'meshz={             '
+          WRITE(IOUT,6012) (DCMESH(IC1,3),IC1=0,MESH(3))
+          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
+*----
+*  For option ITST=0, 1, consider global geometry
+*  and test if line crosses this geometry.
+*----
+      IF(ITST .EQ. 0 .OR. ITST .EQ. 1) THEN
+*----
+*  Scan over Cartesian directions
+*----
+        DO 100 IDIR=1,NDIM
+*----
+*  IDIR is main planes direction ($X$, $Y$ or $Z$)
+*  JDIR is first direction of plane perpendicular
+*       to main direction ($Y, $Z$ or $X$).
+*  KDIR is second direction of plane perpendicular
+*       to main direction ($Z$, $X$ or $Y$).
+*----
+          JDIR=MOD(IDIR,NDIM)+1
+          KDIR=MOD(IDIR+1,NDIM)+1
+*----
+*  Select planes order in direction IDIR
+*----
+          IF(DIRTRK(IDIR) .LT. 0) THEN
+            INEXT=-1
+            IFRST=MESH(IDIR)+1
+            ILAST=1
+          ELSE
+            INEXT=1
+            IFRST=1
+            ILAST=MESH(IDIR)+1
+          ENDIF
+*----
+*  Select planes order in direction JDIR
+*----
+          IF(DIRTRK(JDIR) .LT. 0) THEN
+            JNEXT=-1
+            JFRST=MESH(JDIR)+1
+            JLAST=1
+          ELSE
+            JNEXT=1
+            JFRST=1
+            JLAST=MESH(JDIR)+1
+          ENDIF
+*----
+*  Select planes order in direction KDIR (if required)
+*----
+          IF(NDIM .EQ. 3) THEN
+            IF(DIRTRK(KDIR) .LT. 0) THEN
+              KNEXT=-1
+              KFRST=MESH(KDIR)+1
+              KLAST=1
+            ELSE
+              KNEXT=1
+              KFRST=1
+              KLAST=MESH(KDIR)+1
+            ENDIF
+          ENDIF
+*----
+*  Scan over planes in direction IDIR
+*----
+          DO 101 IFACE=IFRST,ILAST,ILAST-IFRST
+*----
+*  Compute track length required to reach a face
+*  and intersection point with infinite plane
+*----
+            TRKDIS=(DCMESH(IFACE-1,IDIR)-ORITRK(IDIR))/DIRTRK(IDIR)
+            TRKJ=DBLE(JNEXT)*(TRKDIS*DIRTRK(JDIR)+ORITRK(JDIR))
+            IF(NDIM .EQ. 3) THEN
+              TRKK=DBLE(KNEXT)*(TRKDIS*DIRTRK(KDIR)+ORITRK(KDIR))
+            ENDIF
+*----
+*  Test if point is in finite Cartesian plane (j,k)
+*----
+            DP1=DBLE(JNEXT)*DCMESH(JFRST-1,JDIR)
+            DP2=DBLE(JNEXT)*DCMESH(JLAST-1,JDIR)
+            IF(TRKJ .GE. DP1-DCUTOF
+     >   .AND. TRKJ .LE. DP2+DCUTOF ) THEN
+              IF(NDIM .EQ. 3) THEN
+*----
+*  For 3-D, consider intersection with plane
+*----
+                DP1=DBLE(KNEXT)*DCMESH(KFRST-1,KDIR)
+                DP2=DBLE(KNEXT)*DCMESH(KLAST-1,KDIR)
+                IF(TRKK .GE. DP1-DCUTOF
+     >       .AND. TRKK .LE. DP2+DCUTOF ) THEN
+*----
+*  If no point, save distance set NBSINT to 1 and continue
+*  if one point already known, check if this point is
+*  at a different distance exit, otherwise, this point
+*  is already identified (intersection of 2 or more planes)
+*  just continue
+*----
+                  IF(NBSINT .EQ. 0) THEN
+                    NBSINT=1
+                    TRKOLD=TRKDIS
+                  ELSE IF(TRKOLD .NE. TRKDIS) THEN
+                    NBSINT=2
+                    GO TO 105
+                 ENDIF
+                ENDIF
+              ELSE
+*----
+*  For 2-D, intersection with a line is sufficient
+*----
+                IF(NBSINT .EQ. 0) THEN
+                  NBSINT=1
+                  TRKOLD=TRKDIS
+                ELSE IF(TRKOLD .NE. TRKDIS) THEN
+                  NBSINT=2
+                  GO TO 105
+                ENDIF
+              ENDIF
+            ENDIF
+ 101      CONTINUE
+ 100    CONTINUE
+ 105    CONTINUE
+      ELSE
+*----
+*  If ITST=-1, assume that 2 surface are intersected
+*  and continue
+*----
+        NBSINT=2
+      ENDIF
+      NXTLCA=NBSINT
+*----
+*  If line does not intersect any surfaces, exit
+*----
+      IF(NXTLCA .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) THEN
+        IF(NBSINT .EQ. 1) CALL XABORT(NAMSBR//
+     >  ': Invalid line since it intersects only 1 surface')
+        IF(IPRINT .GT. 1000)
+     >   WRITE(IOUT,9000) NBSINT
+      ENDIF
+*----
+*  If ITST=0, return number of surfaces intersected by line
+*----
+      IF(ITST .EQ. 0) RETURN
+*----
+*  If ITST=1 and NXTLCA=2 or ITST=-1
+*  track geometry with submesh
+*----
+      NBSINT=0
+*----
+*  Scan over Cartesian directions
+*----
+      DO 110 IDIR=1,NDIM
+*----
+*  IDIR is main planes direction ($X$, $Y$ or $Z$)
+*  JDIR is first direction of plane perpendicular
+*       to main direction ($Y, $Z$ or $X$).
+*  KDIR is second direction of plane perpendicular
+*       to main direction ($Z$, $X$ or $Y$).
+*----
+        JDIR=MOD(IDIR,NDIM)+1
+        KDIR=MOD(IDIR+1,NDIM)+1
+*----
+*  Select planes order in direction IDIR
+*----
+        IF(DIRTRK(IDIR) .LT. 0) THEN
+          INEXT=-1
+          IFRST=MESH(IDIR)+1
+          ILAST=1
+          IOFF=-1
+        ELSE
+          INEXT=1
+          IFRST=1
+          ILAST=MESH(IDIR)+1
+          IOFF=0
+        ENDIF
+*----
+*  Select planes order in direction JDIR
+*----
+        IF(DIRTRK(JDIR) .LT. 0) THEN
+          JNEXT=-1
+          JFRST=MESH(JDIR)+1
+          JLAST=2
+          JOFF=-1
+        ELSE
+          JNEXT=1
+          JFRST=1
+          JLAST=MESH(JDIR)
+          JOFF=0
+        ENDIF
+*----
+*  Select planes order in direction KDIR (if required)
+*----
+        IF(NDIM .EQ. 3) THEN
+          IF(DIRTRK(KDIR) .LT. 0) THEN
+            KNEXT=-1
+            KFRST=MESH(KDIR)+1
+            KLAST=2
+            KOFF=-1
+          ELSE
+            KNEXT=1
+            KFRST=1
+            KLAST=MESH(KDIR)
+            KOFF=0
+          ENDIF
+        ENDIF
+*----
+*  Scan over planes in direction IDIR
+*----
+        DO 111 IFACE=IFRST,ILAST,INEXT
+          ISDIR=IDIR
+          IF(IFACE .EQ. IFRST .OR. IFACE .EQ. ILAST)
+     >    ISDIR=-ISDIR
+*----
+*  Compute track length required to reach a face
+*  and intersection point with infinite plane
+*----
+          TRKDIS=(DCMESH(IFACE-1,IDIR)-ORITRK(IDIR))/DIRTRK(IDIR)
+          TRKJ=DBLE(JNEXT)*(TRKDIS*DIRTRK(JDIR)+ORITRK(JDIR))
+          IF(NDIM .EQ. 3) THEN
+            TRKK=DBLE(KNEXT)*(TRKDIS*DIRTRK(KDIR)+ORITRK(KDIR))
+          ENDIF
+*----
+*  Test if point is in finite Cartesian plane (j,k) and add to number
+*  of surfaces crossed if this is the case
+*  (j=1 or NPLJ+1 and k=1 or NPLK+1)
+*  For values of j=1,NPLJ and k=1,NPLK a volume was crossed
+*----
+          DO 112 JFACE=JFRST,JLAST,JNEXT
+            JSDIR=JDIR
+            IF(JFACE .EQ. JFRST .OR. JFACE .EQ. JLAST)
+     >      JSDIR=-JSDIR
+            DP1=DBLE(JNEXT)*DCMESH(JFACE-1,JDIR)
+            DP2=DBLE(JNEXT)*DCMESH(JFACE+JNEXT-1,JDIR)
+            IF(TRKJ .GE. DP1-DCUTOF
+     >   .AND. TRKJ .LE. DP2+DCUTOF ) THEN
+              IF(NDIM .EQ. 3) THEN
+*----
+*  For 3-D, consider intersection with plane
+*----
+                DO 113 KFACE=KFRST,KLAST,KNEXT
+                  KSDIR=KDIR
+                  IF(KFACE .EQ. KFRST .OR. KFACE .EQ. KLAST)
+     >            KSDIR=-KSDIR
+                  DP1=DBLE(KNEXT)*DCMESH(KFACE-1,KDIR)
+                  DP2=DBLE(KNEXT)*DCMESH(KFACE+KNEXT-1,KDIR)
+                  IF(TRKK .GE. DP1-DCUTOF
+     >         .AND. TRKK .LE. DP2+DCUTOF ) THEN
+                    NBSINT=NBSINT+1
+                    IF(IPRINT .GT. 1000) THEN
+                      WRITE(IOUT,6011) 'Intersection at     '
+                      WRITE(IOUT,6014) IDIR,IFACE,
+     >                JDIR,JFACE+JOFF,KDIR,KFACE+KOFF
+                      WRITE(IOUT,6015) TRKDIS,DIRTRK
+                    ENDIF
+                    IF(NBSINT .EQ. 1) THEN
+                      ISINT(0,NBSINT)=ISDIR
+                      ISINT(IDIR,NBSINT)=-IFACE
+                      ISINT(JDIR,NBSINT)=JFACE+JOFF
+                      ISINT(KDIR,NBSINT)=KFACE+KOFF
+                      TRKLSI(NBSINT)=TRKDIS
+                    ELSE
+*----
+*  Scan previous distances and reorder if necessary
+*----
+                      NBTINT=NBSINT
+                      DO 120 JSUR=NBTINT-1,1,-1
+                        DELTKD=TRKDIS-TRKLSI(JSUR)
+                        IF(DELTKD .GE. -DCUTOF) 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(ABS(DELTKD) .LE. DCUTOF) THEN
+                            IF(ISDIR .GT. 0 ) THEN
+                              IF(ISINT(0,JSUR) .GT. 0) THEN
+                                ISINT(0,JSUR)=ISINT(0,JSUR)+10*ISDIR
+                                ISINT(IDIR,JSUR)=-IFACE
+                              ENDIF
+                              NBSINT=NBSINT-1
+*      IF(IPRINT .GT. 1000)
+*     >write(6,*) 'SurInt A =',JSUR,IFACE,JFACE,KFACE,
+*     >ISINT(0,JSUR),ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),
+*     >ISINT(KDIR,JSUR),
+*     >TRKDIS,TRKLSI(JSUR)
+                              GO TO 125
+                            ELSE
+                              IF(ISINT(0,JSUR) .GT. 0) THEN
+                                ISINT(0,JSUR)=ISDIR
+                                ISINT(IDIR,JSUR)=-IFACE
+                                ISINT(JDIR,JSUR)=JFACE+JOFF
+                                ISINT(KDIR,JSUR)=KFACE+KOFF
+                                NBSINT=NBSINT-1
+                                GO TO 125
+                              ENDIF
+*      IF(IPRINT .GT. 1000)
+*     >write(6,*) 'SurInt B =',JSUR,IFACE,JFACE,KFACE,
+*     >ISINT(0,JSUR),ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),
+*     >ISINT(KDIR,JSUR),
+*     >TRKDIS,TRKLSI(JSUR)
+                            ENDIF
+                          ENDIF
+                          DO 121 KSUR=NBTINT,JSUR+2,-1
+                            ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                            ISINT(IDIR,KSUR)=ISINT(IDIR,KSUR-1)
+                            ISINT(JDIR,KSUR)=ISINT(JDIR,KSUR-1)
+                            ISINT(KDIR,KSUR)=ISINT(KDIR,KSUR-1)
+                            TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 121                      CONTINUE
+                          ISINT(0,JSUR+1)=ISDIR
+                          ISINT(IDIR,JSUR+1)=-IFACE
+                          ISINT(JDIR,JSUR+1)=JFACE+JOFF
+                          ISINT(KDIR,JSUR+1)=KFACE+KOFF
+                          TRKLSI(JSUR+1)=TRKDIS
+*      IF(IPRINT .GT. 1000)
+*     >write(6,*) 'SurInt C =',JSUR,IFACE,JFACE,KFACE,
+*     >ISINT(0,JSUR),ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),
+*     >ISINT(KDIR,JSUR),
+*     >TRKDIS,TRKLSI(JSUR)
+                          GO TO 125
+                        ENDIF
+ 120                  CONTINUE
+                      DO 122 KSUR=NBTINT,2,-1
+                        ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                        ISINT(IDIR,KSUR)=ISINT(IDIR,KSUR-1)
+                        ISINT(JDIR,KSUR)=ISINT(JDIR,KSUR-1)
+                        ISINT(KDIR,KSUR)=ISINT(KDIR,KSUR-1)
+                        TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 122                  CONTINUE
+                      ISINT(0,1)=ISDIR
+                      ISINT(IDIR,1)=-IFACE
+                      ISINT(JDIR,1)=JFACE+JOFF
+                      ISINT(KDIR,1)=KFACE+KOFF
+                      TRKLSI(1)=TRKDIS
+ 125                  CONTINUE
+                    ENDIF
+                  ENDIF
+ 113            CONTINUE
+              ELSE
+*----
+*  For 2-D, intersection with a line is sufficient
+*----
+                NBSINT=NBSINT+1
+                IF(IPRINT .GT. 1000) THEN
+                  WRITE(IOUT,6011) 'Intersection at     '
+                  WRITE(IOUT,6014) IDIR,IFACE,JDIR,JFACE
+                  WRITE(IOUT,6015) TRKDIS,DIRTRK
+                ENDIF
+                IF(NBSINT .EQ. 1) THEN
+                  ISINT(0,NBSINT)=ISDIR
+                  ISINT(IDIR,NBSINT)=-IFACE
+                  ISINT(JDIR,NBSINT)=JFACE+JOFF
+                  TRKLSI(NBSINT)=TRKDIS
+*      IF(IPRINT .GT. 1000)
+*     >write(6,*) 'Int sur A =',NBSINT,
+*     >ISINT(IDIR,NBSINT),ISINT(JDIR,NBSINT),TRKLSI(NBSINT)
+                ELSE
+*----
+*  Scan previous distances and reorder if necessary
+*----
+                  NBTINT=NBSINT
+                  DO 130 JSUR=NBTINT-1,1,-1
+                    DELTKD=TRKDIS-TRKLSI(JSUR)
+                    IF(DELTKD .GE. -DCUTOF) 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(ABS(DELTKD) .LE. DCUTOF) THEN
+                        IF(ISDIR .GT. 0 ) THEN
+                          IF(ISINT(0,JSUR) .GT. 0) THEN
+                            ISINT(0,JSUR)=ISINT(0,JSUR)+10*ISDIR
+                            ISINT(IDIR,JSUR)=-IFACE
+                          ENDIF
+                          NBSINT=NBSINT-1
+*      IF(IPRINT .GT. 1000)
+*     >write(6,*) 'Int sur A =',JSUR,
+*     >ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),TRKLSI(JSUR)
+                          GO TO 135
+                        ELSE
+                          IF(ISINT(0,JSUR) .GT. 0) THEN
+                            ISINT(0,JSUR)=ISDIR
+                            ISINT(IDIR,JSUR)=-IFACE
+                            ISINT(JDIR,JSUR)=JFACE+JOFF
+                            NBSINT=NBSINT-1
+*      IF(IPRINT .GT. 1000)
+*     >write(6,*) 'Int sur B =',JSUR,
+*     >ISINT(IDIR,JSUR),ISINT(JDIR,JSUR),TRKLSI(JSUR)
+                            GO TO 135
+                          ENDIF
+                        ENDIF
+                      ENDIF
+                      DO 131 KSUR=NBTINT,JSUR+2,-1
+                        ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                        ISINT(IDIR,KSUR)=ISINT(IDIR,KSUR-1)
+                        ISINT(JDIR,KSUR)=ISINT(JDIR,KSUR-1)
+                        TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 131                  CONTINUE
+                      ISINT(0,JSUR+1)=ISDIR
+                      ISINT(IDIR,JSUR+1)=-IFACE
+                      ISINT(JDIR,JSUR+1)=JFACE+JOFF
+                      TRKLSI(JSUR+1)=TRKDIS
+*      IF(IPRINT .GT. 1000)
+*     >write(6,*) 'Int sur C =',JSUR+1,
+*     >ISINT(IDIR,JSUR+1),ISINT(JDIR,JSUR+1),TRKLSI(JSUR+1)
+                      GO TO 135
+                    ENDIF
+ 130              CONTINUE
+                  DO 132 KSUR=NBTINT,2,-1
+                    ISINT(0,KSUR)=ISINT(0,KSUR-1)
+                    ISINT(IDIR,KSUR)=ISINT(IDIR,KSUR-1)
+                    ISINT(JDIR,KSUR)=ISINT(JDIR,KSUR-1)
+                    TRKLSI(KSUR)=TRKLSI(KSUR-1)
+ 132              CONTINUE
+                  ISINT(0,1)=ISDIR
+                  ISINT(IDIR,1)=-IFACE
+                  ISINT(JDIR,1)=JFACE+JOFF
+                  TRKLSI(1)=TRKDIS
+*      IF(IPRINT .GT. 1000)
+*     >write(6,*) 'Int sur D =',1,
+*     >ISINT(IDIR,1),ISINT(JDIR,1),TRKLSI(1)
+ 135              CONTINUE
+                ENDIF
+              ENDIF
+            ENDIF
+ 112      CONTINUE
+ 111    CONTINUE
+ 110  CONTINUE
+      IF(IPRINT .GT. 1000) THEN
+        WRITE(IOUT,6011) 'Intersections'
+        DO 210 JSUR=1,NBSINT
+          WRITE(IOUT,6010) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5)
+ 210    CONTINUE
+      ENDIF
+      IF(NBSINT .GE. 2) THEN
+*----
+*  Find number of corners for initial surface intersections
+*----
+        TRKDIS=TRKLSI(1)
+        NBCOR(1)=1
+        DO 140 JSUR=2,NBSINT
+          DELTKD=TRKLSI(JSUR)-TRKDIS
+          IF(DELTKD .GT. DCUTOF) GO TO 145
+          NBCOR(1)=NBCOR(1)+1
+ 140    CONTINUE
+ 145    CONTINUE
+*----
+*  Find number of corners for final surface intersections
+*----
+        TRKDIS=TRKLSI(NBSINT)
+        NBCOR(2)=1
+        DO 150 JSUR=NBSINT-1,1,-1
+          DELTKD=TRKDIS-TRKLSI(JSUR)
+          IF(DELTKD .GT. DCUTOF) GO TO 155
+          NBCOR(2)=NBCOR(2)+1
+ 150    CONTINUE
+ 155    CONTINUE
+      ENDIF
+*----
+*  Print intersection points
+*----
+      IF(IPRINT .GT. 1000) THEN
+        WRITE(IOUT,6011) 'Initial planes      '
+        DO JSUR=1,NBCOR(1)
+          WRITE(IOUT,6010) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5)
+        ENDDO
+        IF(NBCOR(1)+1 .LE. NBSINT-NBCOR(2)) THEN
+          WRITE(IOUT,6011) 'Intermediate planes '
+          DO JSUR=NBCOR(1)+1,NBSINT-NBCOR(2)
+            WRITE(IOUT,6010) TRKLSI(JSUR),
+     >      (ISINT(IDIR,JSUR),IDIR=0,5)
+          ENDDO
+        ENDIF
+        WRITE(IOUT,6011) 'Final planes        '
+        DO JSUR=NBSINT-NBCOR(2)+1,NBSINT
+          WRITE(IOUT,6010) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5)
+        ENDDO
+      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 .GT. 10) CALL XABORT(NAMSBR//
+     >': Final outer surface identified as an internal corner')
+        IDIRF=ABS(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
+*----
+* Print information for debug
+*----
+          WRITE(IOUT,6011) 'meshx={             '
+          WRITE(IOUT,6012) (DCMESH(IC1,1),IC1=0,MESH(1))
+          WRITE(IOUT,6013)
+          WRITE(IOUT,6011) 'meshy={             '
+          WRITE(IOUT,6012) (DCMESH(IC1,2),IC1=0,MESH(2))
+          WRITE(IOUT,6013)
+          IF(NDIM .EQ. 3) THEN
+            WRITE(IOUT,6011) 'meshz={             '
+            WRITE(IOUT,6012) (DCMESH(IC1,3),IC1=0,MESH(3))
+            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)
+          WRITE(IOUT,6011) 'Initial planes      '
+          DO KSUR=1,NBCOR(1)
+            WRITE(IOUT,6010) TRKLSI(KSUR),
+     >      (ISINT(IDIR,KSUR),IDIR=0,5)
+          ENDDO
+          IF(NBCOR(1)+1 .LE. NBSINT-NBCOR(2)) THEN
+            WRITE(IOUT,6011) 'Intermediate planes '
+            DO KSUR=NBCOR(1)+1,NBSINT-NBCOR(2)
+              WRITE(IOUT,6010) TRKLSI(KSUR),
+     >        (ISINT(IDIR,KSUR),IDIR=0,5)
+            ENDDO
+          ENDIF
+          WRITE(IOUT,6011) 'Final planes        '
+          DO KSUR=NBSINT-NBCOR(2)+1,NBSINT
+            WRITE(IOUT,6010) TRKLSI(KSUR),
+     >      (ISINT(IDIR,KSUR),IDIR=0,5)
+          ENDDO
+          WRITE(IOUT,9001) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5),IDIRF
+          CALL XABORT(NAMSBR//': Invalid final '//CDIR(IDIRS)//
+     >' directed surface')
+        ENDIF
+      ENDDO
+*----
+*  Identify regions
+*----
+      DO JSUR=NBSINT-NBCOR(2)+1,NBCOR(1)+1,-1
+        DO IDIRS=1,5
+          IRFIN(IDIRS)=ABS(ISINT(IDIRS,JSUR-1))
+          IF(ISINT(IDIRS,JSUR)   .LT. 0 ) THEN
+            IF(ISINT(IDIRS,JSUR-1) .LT. 0 ) THEN
+              IRFIN(IDIRS)=-MAX(ISINT(IDIRS,JSUR-1),ISINT(IDIRS,JSUR))
+            ENDIF
+          ELSE
+            IF(ISINT(IDIRS,JSUR-1) .LT. 0 ) THEN
+              IRFIN(IDIRS)=ISINT(IDIRS,JSUR)
+            ENDIF
+          ENDIF
+        ENDDO
+        TRKLSI(JSUR)=TRKLSI(JSUR)-TRKLSI(JSUR-1)
+        DO IDIRS=1,5
+          IF(ISINT(IDIRS,JSUR)   .LT. 0 ) THEN
+            ISINT(IDIRS,JSUR)=IRFIN(IDIRS)
+          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 .GT. 10) CALL XABORT(NAMSBR//
+     >': Initial outer surface identified as an internal corner')
+        IDIRB=ABS(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
+*----
+* Print information for debug
+*----
+          WRITE(IOUT,6011) 'meshx={             '
+          WRITE(IOUT,6012) (DCMESH(IC1,1),IC1=0,MESH(1))
+          WRITE(IOUT,6013)
+          WRITE(IOUT,6011) 'meshy={             '
+          WRITE(IOUT,6012) (DCMESH(IC1,2),IC1=0,MESH(2))
+          WRITE(IOUT,6013)
+          IF(NDIM .EQ. 3) THEN
+            WRITE(IOUT,6011) 'meshz={             '
+            WRITE(IOUT,6012) (DCMESH(IC1,3),IC1=0,MESH(3))
+            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)
+          WRITE(IOUT,6011) 'Initial planes      '
+          DO KSUR=1,NBCOR(1)
+            WRITE(IOUT,6010) TRKLSI(KSUR),
+     >      (ISINT(IDIR,KSUR),IDIR=0,5)
+          ENDDO
+          IF(NBCOR(1)+1 .LE. NBSINT-NBCOR(2)) THEN
+            WRITE(IOUT,6011) 'Intermediate planes '
+            DO KSUR=NBCOR(1)+1,NBSINT-NBCOR(2)
+              WRITE(IOUT,6010) TRKLSI(KSUR),
+     >        (ISINT(IDIR,KSUR),IDIR=0,5)
+            ENDDO
+          ENDIF
+          WRITE(IOUT,6011) 'Final planes        '
+          DO KSUR=NBSINT-NBCOR(2)+1,NBSINT
+            WRITE(IOUT,6010) TRKLSI(KSUR),
+     >      (ISINT(IDIR,KSUR),IDIR=0,5)
+          ENDDO
+         WRITE(IOUT,9001) TRKLSI(JSUR),
+     >    (ISINT(IDIR,JSUR),IDIR=0,5),IDIRB
+          CALL XABORT(NAMSBR//': Invalid initial '//CDIR(IDIRS)//
+     >' directed surface')
+        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(F25.16,2X))
+ 9000 FORMAT(' Warning : ',I10,' surfaces crossed')
+ 9001 FORMAT(1X,F25.16,7I10)
+      END