Initial commit from Polytechnique Montreal

1 files changed, 606 insertions, 0 deletions

diff --git a/Dragon/src/SALTCG.f b/Dragon/src/SALTCG.f
new file mode 100644
index 0000000..123a288
--- /dev/null
+++ b/Dragon/src/SALTCG.f
@@ -0,0 +1,606 @@
+*DECK SALTCG
+      SUBROUTINE SALTCG(IPTRK ,IFTRK ,IPRINT,IGTRK ,NBSLIN, GG)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* To track an assembly of cells containing clusters using the new SALT
+* tracking procedure (based on NXTTCG.f).
+*
+*Copyright:
+* Copyright (C) 2014 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):
+* A. Hebert and G. Marleau
+*
+*Parameters: input
+* IPTRK   pointer to the TRACKING data structure in
+*         update or creation mode.
+* IFTRK   pointer to the TRACKING file in creation mode.
+* IPRINT  print level.
+* IGTRK   flag to generate the tracking file. In the case where
+*         IGTRK=1, the tracking is performed and
+*         used to evaluate the track normalisation factor and the
+*         tracking file is generated. When IGTRK=0, the tracking is
+*         still performed and used to evaluate the track normalisation
+*         factor but the tracking file is not generated.
+* NBSLIN  maximum number of segments in a single tracking line.
+*         computed by default in SALTCG but limited to 100000
+*         elements. This default value can be bypassed using
+*         keyword NBSLIN.
+* GG      geometry basic information.
+*
+*----------
+*
+      USE              GANLIB
+      USE SAL_GEOMETRY_TYPES, ONLY : T_G_BASIC,ISPEC
+      USE SAL_TRACKING_TYPES, ONLY : NMAX2,ITRAC2,RTRAC2,IPART,RPART,
+     >                        NIPART,NRPART
+      IMPLICIT         NONE
+*----
+*  Subroutine arguments
+*----
+      TYPE(C_PTR)      IPTRK
+      INTEGER          IFTRK
+      INTEGER          IPRINT,IGTRK
+      INTEGER          NBSLIN
+      TYPE(T_G_BASIC)  GG
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='SALTCG')
+      INTEGER          NSTATE
+      PARAMETER       (NSTATE=40)
+      INTEGER          NMAX0
+      DOUBLE PRECISION PI,DZERO,DONE,DTWO,DSUM
+      PARAMETER       (PI=3.14159265358979, DZERO=0.0D0,DONE=1.0D0,
+     >                 DTWO=2.0D0,NMAX0=100000)
+*----
+*  Functions
+*----
+      INTEGER          KDROPN,IFTEMP,KDRCLS,ICLS
+*----
+*  Local variables
+*----
+      INTEGER          ISTATE(NSTATE),IEDIMG(NSTATE),ICODE(6)
+      REAL             RSTATT(NSTATE),ALBEDO(6)
+      INTEGER          RENO,LTRK,AZMOAQ,ISYMM,POLQUA,POLOAQ,AZMQUA,
+     >                 AZMNBA,OK
+      DOUBLE PRECISION DENUSR,RCUTOF,DENLIN,SPACLN,WEIGHT
+      DOUBLE PRECISION RADIUS,CENTER(3)
+      INTEGER          NDIM,ITYPBC,IDIRG,NBOCEL,NBUCEL,IDIAG,
+     >                 ISAXIS(3),NOCELL(3),NUCELL(3),MXMSH,MAXMSH,
+     >                 MAXREG,NBTCLS,MAXMSP,MAXRSP,NFSUR,NFREG,
+     >                 MXGSUR,NUNK,NPLANE,NPOINT,NTLINE,NBTDIR,
+     >                 MAXSUB,MAXSGL,NBDR,ILONG,ITYLCM,IPER(3)
+      INTEGER          JJ,KK,NCOR,NQUAD,NANGL,NBANGL,LINMAX
+      DOUBLE PRECISION DQUAD(4),ABSC(3,2),RCIRC,SIDEH,ANGLE
+      CHARACTER        CTRK*4,COMENT*80
+      INTEGER          IFMT,NEREG,NESUR
+*----
+*  Allocatable arrays
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: KEYMRG,MATALB
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NBSANG
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: SURVOL,DGMESH,
+     > DNSANG,DDANG,DVNOR,DSNOR
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DDENWT
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: DANGLT
+*----
+*  Processing starts:
+*  print routine opening header if required
+*  and initialize various parameters.
+*----
+      IF(IPRINT .GE. 1) WRITE(IOUT,6000) NAMSBR
+*----
+*  Open temporary tracking file if required
+*----
+      IF(IGTRK .EQ. 1) THEN
+        IFTEMP= KDROPN('DUMMYSQ',0,2,0)
+        IF(IFTEMP .LE. 0) WRITE(IOUT,9010) NAMSBR
+      ENDIF
+*----
+*  Get state vectors
+*----
+      ISTATE(:NSTATE)=0
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+      CALL LCMGET(IPTRK,'EXCELTRACKOP',RSTATT)
+      NEREG=ISTATE(1)
+      NESUR=ISTATE(5)
+      NUNK=NEREG+NESUR+1
+      RENO  =ISTATE(8)
+      LTRK  =ISTATE(9)
+      AZMOAQ=ISTATE(11)
+      ISYMM =ISTATE(12)
+      POLQUA=ISTATE(13)
+      POLOAQ=ISTATE(14)
+      AZMQUA=ISTATE(15)
+      AZMNBA=ISTATE(16)
+      IFMT=ISTATE(21)
+      DENUSR=DBLE(RSTATT(2))
+      RCUTOF=DBLE(RSTATT(3))
+      DENLIN=DBLE(RSTATT(4))
+      SPACLN=DBLE(RSTATT(5))
+      WEIGHT=RSTATT(6)
+*----
+*  Get main tracking records
+*----
+      CALL LCMLEN(IPTRK,'ICODE       ',ILONG,ITYLCM)
+      IF(ILONG.GT.6) CALL XABORT('SALTCG: ALBEDO OVERFLOW.')
+      CALL LCMGET(IPTRK,'ICODE       ',ICODE )
+      CALL LCMGET(IPTRK,'ALBEDO      ',ALBEDO)
+      CALL LCMSIX(IPTRK,'NXTRecords  ',1)
+*----
+*  Get general dimensioning vector for geometry tracking
+*----
+      IEDIMG(:NSTATE)=0
+      CALL LCMGET(IPTRK,'G00000001DIM',IEDIMG)
+      NDIM     =IEDIMG( 1)
+      ITYPBC   =IEDIMG( 2)
+      IDIRG    =IEDIMG( 3)
+      NBOCEL   =IEDIMG( 4)
+      NBUCEL   =IEDIMG( 5)
+      IDIAG    =IEDIMG( 6)
+      ISAXIS(1)=IEDIMG( 7)
+      ISAXIS(2)=IEDIMG( 8)
+      ISAXIS(3)=IEDIMG( 9)
+      NOCELL(1)=IEDIMG(10)
+      NOCELL(2)=IEDIMG(11)
+      NOCELL(3)=IEDIMG(12)
+      NUCELL(1)=IEDIMG(13)
+      NUCELL(2)=IEDIMG(14)
+      NUCELL(3)=IEDIMG(15)
+      MXMSH    =IEDIMG(16)
+      MAXREG   =IEDIMG(17)
+      NBTCLS   =IEDIMG(18)
+      MAXMSP   =IEDIMG(20)
+      MAXRSP   =IEDIMG(21)
+      NFSUR    =IEDIMG(22)
+      NFREG    =IEDIMG(23)
+      MXGSUR   =IEDIMG(24)
+      MAXMSH=MAX(1,MXMSH,MAXMSP,MAXREG)
+      IF(IPRINT .GE. 1) THEN
+        WRITE(IOUT,6011) NFREG,NEREG,NFSUR,NESUR
+      ENDIF
+      IF((ITYPBC .EQ. 0).OR.(ITYPBC .EQ. 2)) THEN
+*----
+*  Define Cell for periodicity
+*  Cartesian Boundary
+*----
+        IPER(1)=2
+        IF(ABS(ISAXIS(1)) .EQ. 3) IPER(1)=1
+        IPER(2)=2
+        IF(ABS(ISAXIS(2)) .EQ. 3) IPER(2)=1
+        IPER(3)=2
+        IF(ABS(ISAXIS(3)) .EQ. 3) IPER(3)=1
+*----
+*  Use intrinsic geometry symmetries
+*  to simplify tracking unless
+*  NOSY tracking option activated
+*----
+        IF(ISYMM .NE. 0) THEN
+          ISYMM=0
+          IF(ABS(ISAXIS(1)) .EQ. 1 .OR. ABS(ISAXIS(1)) .EQ. 2) THEN
+*----
+*  X SYMMETRY
+*----
+            ISYMM=2
+          ENDIF
+          IF(ABS(ISAXIS(2)) .EQ. 1 .OR. ABS(ISAXIS(2)) .EQ. 2) THEN
+*----
+*  Y SYMMETRY
+*----
+            ISYMM=4+2*ISYMM
+          ENDIF
+          IF(NDIM .EQ. 3) THEN
+            IF(ABS(ISAXIS(3)) .EQ. 1 .OR. ABS(ISAXIS(3)) .EQ. 2) THEN
+*----
+*  Z SYMMETRY
+*----
+              ISYMM=16+ISYMM
+            ENDIF
+          ENDIF
+          IF(ISYMM .EQ. 0) ISYMM=1
+        ENDIF
+      ENDIF
+*----
+*  Read global mesh for geometry
+*----
+      ALLOCATE(DGMESH((MAXMSH+2)*4))
+      IF(ITYPBC.EQ.0) THEN
+        CALL NXTXYZ(IPTRK ,IPRINT,NDIM  ,ITYPBC,MAXMSH,NUCELL,
+     >              ABSC,DGMESH)
+      ELSE IF(ITYPBC.GE.2) THEN
+        ! hexagonal geometry
+        CALL LCMGET(IPTRK,'G00000001SMX',DGMESH)
+        SIDEH=DGMESH(2)-DGMESH(1)
+        ABSC(:2,1)=SIDEH*SQRT(3.0)
+        ABSC(:2,2)=ABSC(:2,1)
+      ELSE
+        CALL XABORT(NAMSBR//': geometry not implemented')
+      ENDIF
+*----
+*  Verify tracking parameters and compute number of angles
+*  associated with angular order and spatial quadrature parameters
+*  1. Isotropic tracking
+*----
+      NCOR= 1
+      NPLANE=1
+      IF(LTRK .EQ. 0) THEN
+        IF(NDIM .EQ. 3) THEN
+          CALL XABORT(NAMSBR//': 3-D geometry is not allowed')
+        ENDIF
+        NQUAD=2
+        DQUAD(1)=DONE
+        DQUAD(2)=DONE
+        NANGL=AZMOAQ
+        NBANGL=NANGL
+        IF(ISYMM .EQ. 2 .OR. ISYMM .EQ. 8) THEN
+          DQUAD(1)=DONE/DTWO
+          DQUAD(2)=DZERO
+        ENDIF
+        ALLOCATE(DANGLT(NDIM,NBANGL,NQUAD),DDENWT(NBANGL,NQUAD))
+        CALL NXTQAS(IPRINT,NDIM  ,AZMQUA,NANGL ,NQUAD ,NBANGL,
+     >              DQUAD ,DANGLT,DDENWT)
+        LINMAX=NBUCEL*(4+MXGSUR+16)
+*----
+*  Select standard spatial tracking parameters
+*----
+        CALL NXTQSS(IPRINT,NDIM  ,ITYPBC,MAXMSH,NUCELL,DENUSR,
+     >              DGMESH,NPLANE,NPOINT,DENLIN,SPACLN,WEIGHT,
+     >              RADIUS,CENTER)
+        CALL LCMPUT(IPTRK,'TrackingDirc',NDIM*NBANGL*NQUAD,4,DANGLT)
+      ELSE
+*----
+*  2. Specular tracking
+*----
+        NPOINT=0
+        NQUAD=2
+        IF(NDIM .EQ. 3) CALL XABORT(NAMSBR//
+     >  ': TSPC option not valid for 3-D geometries')
+        IF(ITYPBC .EQ. 0) THEN
+          ! Cartesian geometries
+          IF(AZMOAQ .GT. 24) THEN
+            WRITE(IOUT,9002) NAMSBR,AZMOAQ,24,30
+            AZMOAQ=30
+          ELSE IF(AZMOAQ .GT. 20) THEN
+            IF(AZMOAQ .NE. 24) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,20,24,24
+              AZMOAQ=24
+            ENDIF
+          ELSE IF(AZMOAQ .GT. 18) THEN
+            IF(AZMOAQ .NE. 20) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,18,20,20
+              AZMOAQ=20
+            ENDIF
+          ELSE IF(AZMOAQ .GT. 14) THEN
+            IF(AZMOAQ .NE. 18) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,14,18,18
+              AZMOAQ=18
+            ENDIF
+          ELSE IF(AZMOAQ .GT. 12) THEN
+            IF(AZMOAQ .NE. 14) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,12,14,14
+              AZMOAQ=14
+            ENDIF
+          ELSE IF(AZMOAQ .GT. 8) THEN
+            IF(AZMOAQ .NE. 12) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,8,12,12
+              AZMOAQ=12
+            ENDIF
+          ELSE IF(AZMOAQ .GT. 6) THEN
+            IF(AZMOAQ .NE. 8) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,6,8,8
+              AZMOAQ=8
+            ENDIF
+          ELSE IF(AZMOAQ .GT. 2) THEN
+            IF(AZMOAQ .NE. 6) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,2,6,6
+              AZMOAQ=6
+            ENDIF
+          ELSE IF(AZMOAQ .GE. 0) THEN
+            IF(AZMOAQ .NE. 2) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,-1,2,2
+              AZMOAQ=2
+            ENDIF
+          ENDIF
+        ELSE IF(ITYPBC .GE.2) THEN
+          ! hexagonal geometries
+          IF(AZMOAQ .GT. 12) THEN
+            WRITE(IOUT,9002) NAMSBR,AZMOAQ,12,18
+            AZMOAQ=18
+          ELSE IF(AZMOAQ .GT. 6) THEN
+            IF(AZMOAQ .NE. 12) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,6,12,12
+              AZMOAQ=12
+            ENDIF
+          ELSE IF(AZMOAQ .GT. 3) THEN
+            IF(AZMOAQ .NE. 6) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,3,6,6
+              AZMOAQ=6
+            ENDIF
+          ELSE IF(AZMOAQ .GT. 1) THEN
+            IF(AZMOAQ .NE. 3) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,1,3,3
+              AZMOAQ=3
+            ENDIF
+          ELSE IF(AZMOAQ .GE. 0) THEN
+            IF(AZMOAQ .NE. 1) THEN
+              WRITE(IOUT,9003) NAMSBR,AZMOAQ,-1,1,1
+              AZMOAQ=1
+            ENDIF
+          ENDIF
+        ENDIF
+        NBANGL=AZMOAQ
+        NANGL =AZMOAQ
+        ALLOCATE(NBSANG(5,NBANGL))
+        ALLOCATE(DANGLT(NDIM,NBANGL,4),DDENWT(NBANGL,4),DNSANG(NBANGL),
+     >  DDANG(NBANGL))
+        LINMAX=8*NANGL*NBUCEL*(4+MXGSUR+16)
+        RCIRC=SQRT(ABSC(1,1)**2+ABSC(2,1)**2)
+        ABSC(1,1)= ABSC(1,1)/RCIRC
+        ABSC(2,1)= ABSC(2,1)/RCIRC
+        CALL NXTQAC(IPRINT,NDIM  ,NANGL ,NBANGL,ITYPBC,DENUSR,
+     >              ABSC  ,RCIRC ,AZMQUA,IPER  ,DANGLT,DDENWT,
+     >              DNSANG,NBSANG,DDANG)
+        DEALLOCATE(DDANG)
+        DO JJ=1,NBANGL
+          DANGLT(1,NBANGL-JJ+1,2)=-DANGLT(1,JJ,1)
+          DANGLT(2,NBANGL-JJ+1,2)=DANGLT(2,JJ,1)
+          DDENWT(NBANGL-JJ+1,2)=DDENWT(JJ,1)
+        ENDDO
+        DO JJ=1,NBANGL
+          DANGLT(1,NBANGL-JJ+1,4)=DANGLT(1,JJ,1)
+          DANGLT(2,NBANGL-JJ+1,4)=-DANGLT(2,JJ,1)
+          DDENWT(NBANGL-JJ+1,4)=DDENWT(JJ,1)
+          DANGLT(1,NBANGL-JJ+1,3)=DANGLT(1,JJ,2)
+          DANGLT(2,NBANGL-JJ+1,3)=-DANGLT(2,JJ,2)
+          DDENWT(NBANGL-JJ+1,3)=DDENWT(JJ,2)
+        ENDDO
+        IF(IPRINT.GT.1) THEN
+          WRITE(IOUT,'(/34H SALTCG: CYCLIC ANGULAR QUADRATURE/8X,
+     >    5HANGLE,8X,7HCOSINES,16(1h-),3X,6HWEIGHT,11X,
+     >    17HCYCLIC PARAMETERS)')
+          DSUM=0.D0
+          DO KK=1,4
+            DO JJ=1,NANGL
+              ANGLE=SIGN(ACOS(DANGLT(1,JJ,KK))/PI*180.0,DANGLT(2,JJ,KK))
+              WRITE(IOUT,'(1X,I4,1P,4E13.4,5X,2I4)') (KK-1)*NANGL+JJ,
+     >        ANGLE,DANGLT(:2,JJ,KK),0.5D0/DDENWT(JJ,KK),NBSANG(:2,JJ)
+              DSUM=DSUM+0.5D0/DDENWT(JJ,KK)
+            ENDDO
+          ENDDO
+          WRITE(IOUT,'(39X,5HDSUM=,1P,E13.4)') DSUM
+        ENDIF
+        CALL LCMPUT(IPTRK,'TrackingDirc',NDIM*NBANGL*4,4,DANGLT)
+        CALL LCMPUT(IPTRK,'TrackingTrkW',NBANGL*4,4,DDENWT)
+        CALL LCMPUT(IPTRK,'TrackingSpaD',NBANGL,4,DNSANG)
+        CALL LCMPUT(IPTRK,'TrackingNbST',5*NBANGL,1,NBSANG)
+      ENDIF
+      RSTATT(4)=REAL(DENLIN)
+      RSTATT(5)=REAL(SPACLN)
+      RSTATT(6)=REAL(WEIGHT)
+      RSTATT(7)=REAL(RADIUS)
+      RSTATT(8)=REAL(CENTER(1))
+      RSTATT(9)=REAL(CENTER(2))
+*----
+*  Allocate memory to hold tracking data
+*----
+**      tracking data buffer:
+*       integers
+*       ITRAC2(NMAX OR 2*NMAX) = integer tracking array
+*
+*       *integer descriptors in itrac2:
+*           1 = address of last data
+*           2 = total number of sub-trajectories
+*           3 = 
+*           4 = phi for trajectory (2D)
+*
+*       reals
+*       RTRAC2(nmax or 2*nmax) = real tracking array
+*
+*       *real descriptors:
+*           1 =
+*           2 = cos phi entering basic
+*           3 = sin phi left surface
+*           4 = sin phi right surface
+*           5 = cos phi left surface
+*           6 = cos phi right surface
+*           7 = total weight (DELR*WPHI)
+*           8 = radial weight (DELR)
+*
+*       IPART(NIPART,MXELEM)   = to store integer intersection data
+*       RPART(NRPART,MXELEM)   = to store real intersection data
+*
+      IF(NBSLIN <= 0)THEN
+         NMAX2=NMAX0
+      ELSE
+         NMAX2=NBSLIN
+      ENDIF
+      IF(ISPEC == 1) NMAX2=NMAX2*100
+      ALLOCATE(ITRAC2(2*NMAX2),IPART(NIPART,GG%NB_ELEM),RTRAC2(NMAX2),
+     1 RPART(NRPART,GG%NB_ELEM),STAT=OK)
+      IF(OK/=0) CALL XABORT('SALTCG: not enough memory IRD')
+*----
+*  Track
+*----
+      LINMAX=MAX(LINMAX,NBSLIN)
+      IF(IPRINT .GE. 10) WRITE(IOUT,6010) LINMAX
+      NBDR=1
+      IF(RENO .EQ. -1) THEN
+        IF(LTRK .EQ. 0) THEN
+          NBDR=NQUAD*NBANGL+1
+        ELSE IF(LTRK .EQ. 1) THEN
+          NBDR=4*NBANGL+1
+        ENDIF
+      ENDIF
+      ALLOCATE(DVNOR(NFREG*NBDR),DSNOR(NFSUR*NQUAD*NBANGL))
+      IF(LTRK .EQ. 0) THEN
+*----
+*  Standard (isotropic) tracking (white boundary conditions)
+*----
+        IF(IPRINT .GE. 1) WRITE(IOUT,6030) NBANGL*NQUAD,NPOINT
+        MAXSUB=1
+        CALL SALTLS(IFTEMP,IPRINT,IGTRK ,NFREG ,NBANGL, NQUAD ,
+     >              RENO  ,NBDR  ,IFMT  ,DENUSR,DANGLT, DDENWT,
+     >              GG    ,NBTDIR,MAXSGL,NTLINE,DVNOR )
+      ELSE
+*----
+*  Cyclic (specular) tracking (mirror like boundary conditions)
+*----
+        NBTDIR=0
+        CALL SALTLC(IFTEMP,IPRINT,IGTRK,NDIM,NFREG,NBANGL,RENO,NBDR,
+     >              IFMT,DENUSR,DANGLT,DDENWT,NBSANG,GG,MAXSUB,MAXSGL,
+     >              NTLINE,DVNOR )
+      ENDIF
+*----
+*  Release allocated memory for SALT tracking
+*----
+      DEALLOCATE(RPART,RTRAC2,IPART,ITRAC2,STAT =OK)
+      IF(OK /= 0) CALL XABORT('SALTCG: failure to deallocate storage')
+*----
+*  Save track normalisation vector
+*----
+      CALL LCMPUT(IPTRK,'NumMerge    ',NFREG,1,GG%NUM_MERGE)
+      CALL LCMPUT(IPTRK,'VolMerge    ',NFREG,4,GG%VOL_NODE)
+      CALL LCMPUT(IPTRK,'VTNormalize ',NFREG,4,DVNOR)
+      IF(NBDR .GT. 1) THEN
+        CALL LCMPUT(IPTRK,'VTNormalizeD',NFREG*(NBDR-1),4,
+     >              DVNOR(NFREG+1))
+      ENDIF
+*----
+*  Get cell description of geometry
+*----
+      ALLOCATE(KEYMRG(-NFSUR:NFREG),MATALB(-NFSUR:NFREG),
+     > SURVOL(-NFSUR:NFREG))
+      CALL LCMLEN(IPTRK,'KEYMRG      ',ILONG,ITYLCM)
+      IF(ILONG>NFREG+NFSUR+1) CALL XABORT('SALTCG: KEYMRG OVERLOW.')
+      CALL LCMGET(IPTRK,'KEYMRG      ',KEYMRG)
+      CALL LCMGET(IPTRK,'MATALB      ',MATALB)
+      CALL LCMGET(IPTRK,'SAreaRvolume',SURVOL)
+*----
+*  Build NXTRecords directory
+*----
+      CALL LCMSIX(IPTRK,'NXTRecords  ',2)
+      ISTATE(12)=ISYMM
+      ISTATE(14)=POLOAQ
+      ISTATE(17)=NPOINT
+      ISTATE(18)=LINMAX
+      ISTATE(19)=NTLINE
+      ISTATE(20)=NBTDIR
+      IF(LTRK .EQ. 0) THEN
+        ISTATE(21)=NQUAD*NBANGL
+      ELSE IF(LTRK .EQ. 1) THEN
+        ISTATE(21)=4*NBANGL
+      ENDIF
+      ISTATE(22)=NPLANE
+      CALL LCMPUT(IPTRK,'STATE-VECTOR',NSTATE,1,ISTATE)
+      CALL LCMPUT(IPTRK,'EXCELTRACKOP',NSTATE,2,RSTATT)
+*----
+*  Renormalize tracks if required and transfer to final tracking file
+*----
+      IF(IGTRK .EQ. 1) THEN
+        CTRK  = '$TRK'
+        WRITE(IFTRK) CTRK,5,NTLINE,IFMT
+        COMENT='CREATOR     : DRAGON'
+        WRITE(IFTRK) COMENT
+        COMENT='MODULE      : SALTCG'
+        WRITE(IFTRK) COMENT
+        COMENT='TYPE        : CARTESIAN'
+        WRITE(IFTRK) COMENT
+        IF(RENO .EQ. -1) THEN
+          COMENT='TRKNOR      : Directional '
+        ELSE IF(RENO .EQ. 0) THEN
+          COMENT='TRKNOR      : Global      '
+        ELSE
+          COMENT='TRKNOR      : Off         '
+        ENDIF
+        WRITE(IFTRK) COMENT
+        IF(IFMT .EQ. 1) THEN
+          COMENT='OPTION      : Extended    '
+          WRITE(IFTRK) COMENT
+        ELSE
+          COMENT='OPTION      : Short       '
+          WRITE(IFTRK) COMENT
+        ENDIF
+*----
+*  Compress VOLSUR and MATALB according to KEYMRG and save on IFTRK
+*----
+        IF(LTRK .EQ. 0) THEN
+          WRITE(IFTRK) NDIM,LTRK,NEREG,NESUR,6,NCOR,NQUAD*NBANGL,MAXSUB,
+     >    MAXSGL
+        ELSE IF(LTRK .EQ. 1) THEN
+          WRITE(IFTRK) NDIM,LTRK,NEREG,NESUR,6,NCOR,4*NBANGL,MAXSUB,
+     >    MAXSGL
+        ENDIF
+        KEYMRG(1:NFREG)=GG%NUM_MERGE(:NFREG)
+        SURVOL(1:NFREG)=GG%VOL_NODE(:NFREG)
+        MATALB(1:NFREG)=GG%MED(:NFREG)
+        CALL NXTCVM(IFTRK,IPRINT,NFREG,NFSUR,NEREG,NESUR,MATALB,SURVOL,
+     >              KEYMRG)
+        WRITE(IFTRK) ( ICODE(JJ),JJ=1,6)
+        WRITE(IFTRK) (ALBEDO(JJ),JJ=1,6)
+        IF(LTRK .EQ. 0) THEN
+          CALL NXTSQD(IFTRK,IPRINT,NDIM,NQUAD,NBANGL,DANGLT,DDENWT)
+        ELSE IF(LTRK .EQ. 1) THEN
+          WRITE(IFTRK) ((DANGLT(1,JJ,KK),DANGLT(2,JJ,KK),JJ=1,NBANGL),
+     >    KK=1,4)
+          WRITE(IFTRK) ((DDENWT(JJ,KK),JJ=1,NBANGL),KK=1,4)
+        ENDIF
+        REWIND IFTEMP
+        CALL NXTTNS(IFTRK ,IFTEMP,IPRINT,RENO  ,NFSUR ,NFREG ,NDIM  ,
+     >              MAXSUB,MAXSGL,NTLINE,NBDR  ,IFMT  ,KEYMRG,DVNOR)
+*----
+*  Close temporary tracking file if required
+*----
+        ICLS=KDRCLS(IFTEMP,2)
+        IF(ICLS .NE. 0) WRITE(IOUT,9011) NAMSBR
+      ENDIF
+*----
+*  Deallocate memory
+*----
+      DEALLOCATE(DSNOR,DVNOR)
+      IF(LTRK .EQ. 1) DEALLOCATE(NBSANG,DNSANG)
+      DEALLOCATE(DDENWT,DANGLT)
+      DEALLOCATE(DGMESH,SURVOL,MATALB,KEYMRG)
+*----
+*  Processing finished:
+*  print routine closing output header if required
+*  and return
+*----
+      IF(IPRINT .GE. 1) THEN
+        WRITE(IOUT,6012)
+        WRITE(IOUT,6001) NAMSBR
+      ENDIF
+      RETURN
+*----
+*  Output formats
+*----
+ 6000 FORMAT('(* Output from --',A6,'-- follows ')
+ 6001 FORMAT('   Output from --',A6,'-- completed *)')
+ 6010 FORMAT(' Maximum length of a line =',I10)
+ 6011 FORMAT(' Tracking of geometry begins:'/
+     >       ' Number of regions before merge =',I10/
+     >       ' Number of regions after merge  =',I10/
+     >       ' Number of surfaces before merge=',I10/
+     >       ' Number of surfaces after merge =',I10)
+ 6012 FORMAT(' Tracking of geometry completed')
+ 6030 FORMAT(' Number of directions for tracking = ',I10/
+     >       ' Number of lines per direction     = ',I10)
+ 9002 FORMAT(' ***** Warning in ',A6,' *****'/
+     >       '       Number of specular angles requested :',I10/
+     >       '       For values > ',I10,' use ',I10)
+ 9003 FORMAT(' ***** Warning in ',A6,' *****'/
+     >       '       Number of specular angles requested :',I10/
+     >       '       For values > ',I10,' and < ',I10,' use ',I10)
+ 9010 FORMAT(' ***** Warning in ',A6,' *****'/
+     >       '       Impossible to open temporary tracking file ')
+ 9011 FORMAT(' ***** Warning in ',A6,' *****'/
+     >       '       Impossible to close temporary tracking file ')
+      END