Initial commit from Polytechnique Montreal

1 files changed, 675 insertions, 0 deletions

diff --git a/Dragon/src/MCGASM.f b/Dragon/src/MCGASM.f
new file mode 100644
index 0000000..d58660d
--- /dev/null
+++ b/Dragon/src/MCGASM.f
@@ -0,0 +1,675 @@
+*DECK MCGASM
+      SUBROUTINE MCGASM(SUBPJJ,SUBDS2,SUBDSP,SUBDSC,IPTRK,KPSYS,IPRINT,
+     1                  IFTRAK,NANI,NGEFF,NFI,NREG,NLONG,M,NMU,NANGL,
+     2                  N2MAX,LC,NDIM,NGIND,CYCLIC,ISCR,CAZ0,CAZ1,CAZ2,
+     3                  CPO,LC0,PACA,LPS,LTMT,NPJJM,LACA,LPJJ,LPJJAN,
+     4                  SIGAL,LPRISM,N2REG,N2SOU,NZP,DELU,FACSYM,ISTRM)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Preconditioning matrices calculation based on the Algebraic Collapsing
+* Acceleration developed by I. R. Suslov and R. Le Tellier 
+* or Self-Collision Probabilities acceleration developed by G.J. Wu 
+* and R. Roy.
+*
+*Copyright:
+* Copyright (C) 2002 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): I. Suslov and R. Le Tellier
+*
+*Parameters: input/output
+* SUBPJJ  PJJ calculation subroutine.
+* SUBDS2  ACA coefficients summation subroutine.
+* SUBDSP  ACA coefficients position subroutine.
+* SUBDSC  ACA coefficients calculation subroutine.
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* KPSYS   pointer array for each group properties.
+* IPRINT  print parameter (equal to zero for no print).
+* IFTRAK  tracking file unit number if IOFSET=0.
+* NANI    number of Legendre orders.
+* NGEFF   number of groups to process.
+* NFI     total number of volumes and surfaces for which specific values
+*         of the neutron flux and reactions rates are required.
+* NREG    number of volumes for which specific values
+*         of the neutron flux and reactions rates are required.
+* NLONG   order of the corrective system.
+* M       number of material mixtures.
+* NMU     order of the polar quadrature in 2D / 1 in 3D.
+* NANGL   number of tracking angles in the plane.
+* N2MAX   maximum number of elements in a track.
+* LC      dimension of vector MCU.
+* NDIM    number of dimensions for the geometry.
+* NGIND   index of the groups to process.
+* CYCLIC  flag set to .true. for cyclic tracking.
+* ISCR    SCR preconditionning flag.
+* CAZ0    cosines of the tracking polar angles in 3D.
+* CAZ1    first cosines of the different tracking azimuthal angles.
+* CAZ2    second cosines of the different tracking azimuthal angles.
+* CPO     cosines of the different tracking polar angles in 2D.
+* PACA    type of preconditioner to solve the ACA corrective system.
+* LC0     used in ILU0-ACA acceleration.
+* LPS     dimension of JS.
+* LTMT    tracking merging flag.
+* NPJJM   number of pjj modes to store for STIS option.
+* LACA    ACA flag.
+* LPJJ    PJJ flag.
+* LPJJAN  anisotropic PJJ flag.
+* SIGAL   albedos and total cross sections array.
+* LPRISM  3D prismatic extended tracking flag.
+* N2REG   number of regions in the 2D tracking if LPRISM.
+* N2SOU   number of external surfaces in the 2D tracking if LPRISM.
+* NZP     number of z-plans if LPRISM.
+* DELU    input track spacing for 3D track reconstruction if LPRISM.
+* FACSYM  tracking symmetry factor for maximum track length if LPRISM.
+* ISTRM   type of streaming effect:
+*         =1 no streaming effect;
+*         =2 isotropic streaming effect;
+*         =3 anisotropic streaming effect.
+*
+*Reference:
+*  Igor R. Suslov, "An Algebraic Collapsing Acceleration in Long 
+*  Characteristics Transport Theory" Proc. of 11-th Symposium of 
+*  AER, 178/9-188, Csopak, September 2001.
+*  \\\\
+*  Igor R. Suslov, "Solution of Transport Equation in 2- and 3-
+*  dimensional Irregular Geometry by the Method of Characteristics"
+*  Int. Conf. Mathematical. Methods and Supercomputing in Nuclear
+*  Applications, Karlsruhe, 1994.
+*  \\\\
+*  G.J. Wu and R. Roy, "Acceleration Techniques for Trajectory-based
+*  Deterministic 3D Transport Solvers",
+*  Ann. Nucl. Energy, 30, 567-583 (2003).
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,KPSYS(NGEFF)
+      INTEGER IPRINT,IFTRAK,NGEFF,NFI,NREG,NLONG,M,NMU,NANGL,N2MAX,LC,
+     1 NDIM,NGIND(NGEFF),LC0,PACA,LPS,NPJJM,N2REG,N2SOU,NZP,ISTRM
+      REAL CPO(NMU),SIGAL(-6:M,NGEFF),DELU,FACSYM
+      DOUBLE PRECISION CAZ0(NANGL),CAZ1(NANGL),CAZ2(NANGL)
+      LOGICAL LTMT,LACA,LPJJ,LPJJAN,LPRISM,LFORC,CYCLIC
+      EXTERNAL SUBPJJ,SUBDS2,SUBDSP,SUBDSC
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) JPSYS
+      INTEGER NCODE(6),SSYM
+      REAL T1,T2,T3,XMUANG(1)
+      DOUBLE PRECISION WEIGHT,WEIGHT0
+      CHARACTER TEXT4*4
+      INTEGER, TARGET, SAVE, DIMENSION(1) :: IDUMMY
+      INTEGER, TARGET, SAVE, DIMENSION(1,1) :: I2DUMMY
+      REAL, TARGET, SAVE, DIMENSION(1) :: DUMMY
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, POINTER, DIMENSION(:) :: NZON,KM,IM,MCU,NZONA,IPERM,
+     1 JU,IM0,MCU0,IS,JS
+      INTEGER, POINTER, DIMENSION(:,:) :: PJJIND
+      REAL, POINTER, DIMENSION(:) :: ZMU,WZMU,V,VA
+      TYPE(C_PTR) :: ZMU_PTR,WZMU_PTR,NZON_PTR,V_PTR,KM_PTR,IM_PTR,
+     1 MCU_PTR,NZONA_PTR,VA_PTR,IPERM_PTR,JU_PTR,IM0_PTR,MCU0_PTR,
+     2 PJJIND_PTR,IS_PTR,JS_PTR
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: NOM,NOM0,PREV,NEXT,NOM3D,
+     1 NOM3D0,KANGL,KEYANI
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ISGNR
+      REAL, ALLOCATABLE, DIMENSION(:) :: ZMUA,WZMUA,PJJ,PSJ,XSW,CQ,
+     1 DIAGQ,DIAGFR,CFR,WORK,LUDF,LUCF,PJJX,PJJY,PJJZ,PJJXI,PJJYI,
+     2 PJJZI,PSJX,PSJY,PSJZ
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: RHARM
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: TRHAR
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: H,H0,HH,H3D,H3D0,
+     1 PJJD,PJJDX,PJJDY,PJJDZ,PJJDXI,PJJDYI,PJJDZI
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DIAGF,CF
+      INTEGER, POINTER, DIMENSION(:) :: INDREG
+      REAL, POINTER, DIMENSION(:) :: ZZZ
+      DOUBLE PRECISION, POINTER, DIMENSION(:) :: VNORF,CMU,CMUI,SMU,
+     1 SMUI,TMU,TMUI
+      TYPE(C_PTR) :: INDREG_PTR,ZZZ_PTR,VNORF_PTR,CMU_PTR,CMUI_PTR,
+     1 SMU_PTR,SMUI_PTR,TMU_PTR,TMUI_PTR
+*----
+*  INITIALIZE POINTERS
+*----
+      KM=>IDUMMY
+      IM=>IDUMMY
+      MCU=>IDUMMY
+      NZONA=>IDUMMY
+      VA=>DUMMY
+      IPERM=>IDUMMY
+      JU=>IDUMMY
+      IM0=>IDUMMY
+      MCU0=>IDUMMY
+      IS=>IDUMMY
+      JS=>IDUMMY
+      PJJIND=>I2DUMMY
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(DIAGF(NLONG,NGEFF),CF(LC,NGEFF))
+*----
+*  TRACKING INFORMATION PINNING AND MEMORY ALLOCATION
+*   NZON    index-number of the mixture type assigned to each volume.
+*   NZONA   index-number of the mixture type assigned to each volume for
+*           ACA.
+*   ISGNR   array of the spherical harmonics signs for the different
+*           reflections.
+*   V       volumes and surfaces.
+*   VA      volumes and surfaces reordered for ACA.
+*   ZMU     polar quadrature set in 2D.
+*   WZMU    polar quadrature set in 2D.
+*   KM      used in ACA acceleration.
+*   MCU     used in ACA acceleration.
+*   IM      used in ACA acceleration.
+*   JU      used in ACA  acceleration for ilu0.
+*   IPERM   permutation array for the unknowns of the corrective system
+*           for ilu0.
+*   LC0     used in ILU0-ACA acceleration.
+*   IM0     used in ILU0-ACA acceleration.
+*   MCU0    used in ILU0-ACA acceleration.
+*   IS      arrays for surfaces neighbors
+*   JS      JS(IS(ISOUT)+1:IS(ISOUT+1)) give the neighboring regions to
+*           surface ISOUT.
+*   PJJIND  index of the modes for STIS option.
+*---- 
+*---
+* GENERATE ALL SIGNS FOR SPHERICAL HARMONICS
+*---
+      IF(NDIM.EQ.1) THEN
+         NFUNL=NANI
+         NMOD=2
+      ELSE IF((.NOT.LPRISM).AND.(NDIM.EQ.2)) THEN
+         NFUNL=NANI*(NANI+1)/2
+         NMOD=4
+      ELSE ! NDIM.EQ.3
+         NFUNL=NANI*NANI
+         NMOD=8
+      ENDIF
+      ALLOCATE(ISGNR(NMOD,NFUNL),KEYANI(NFUNL))
+      CALL MOCIK3(NANI-1,NFUNL,NMOD,ISGNR,KEYANI) 
+      DEALLOCATE(KEYANI)
+*---
+* ASSEMBLY OF SCR AND ACA MATRICES
+*---
+*     recover polar quadrature
+      CALL LCMGPD(IPTRK,'ZMU$MCCG',ZMU_PTR)
+      CALL LCMGPD(IPTRK,'WZMU$MCCG',WZMU_PTR)
+      CALL C_F_POINTER(ZMU_PTR,ZMU,(/ NMU /))
+      CALL C_F_POINTER(WZMU_PTR,WZMU,(/ NMU /))
+*     recover MATALB and VOLSUR
+      CALL LCMGPD(IPTRK,'NZON$MCCG',NZON_PTR)
+      CALL LCMGPD(IPTRK,'V$MCCG',V_PTR)
+      CALL C_F_POINTER(NZON_PTR,NZON,(/ NFI /))
+      CALL C_F_POINTER(V_PTR,V,(/ NFI /))
+      IF(LACA) THEN
+*        recover connection matrices
+         CALL LCMGPD(IPTRK,'KM$MCCG',KM_PTR)
+         CALL LCMGPD(IPTRK,'IM$MCCG',IM_PTR)
+         CALL LCMGPD(IPTRK,'MCU$MCCG',MCU_PTR)
+         CALL C_F_POINTER(KM_PTR,KM,(/ NFI /))
+         CALL C_F_POINTER(IM_PTR,IM,(/ NLONG+1 /))
+         CALL C_F_POINTER(MCU_PTR,MCU,(/ LC /))
+*        recover modified MATALB and VOLSUR
+         CALL LCMGPD(IPTRK,'NZONA$MCCG',NZONA_PTR)
+         CALL LCMGPD(IPTRK,'VA$MCCG',VA_PTR)
+         CALL C_F_POINTER(NZONA_PTR,NZONA,(/ NFI /))
+         CALL C_F_POINTER(VA_PTR,VA,(/ NFI /))
+*        recover permutation array
+         CALL LCMGPD(IPTRK,'INVPI$MCCG',IPERM_PTR)
+         CALL C_F_POINTER(IPERM_PTR,IPERM,(/ NFI /))
+         IF(PACA.GE.2) THEN
+            CALL LCMGPD(IPTRK,'JU$MCCG',JU_PTR)
+            CALL C_F_POINTER(JU_PTR,JU,(/ NLONG /))
+            IF(PACA.EQ.3) THEN
+               CALL LCMLEN(IPTRK,'IM0$MCCG',ILONG1,ITYLCM)
+               CALL LCMLEN(IPTRK,'MCU0$MCCG',ILONG2,ITYLCM)
+               CALL LCMGPD(IPTRK,'IM0$MCCG',IM0_PTR)
+               CALL LCMGPD(IPTRK,'MCU0$MCCG',MCU0_PTR)
+               CALL C_F_POINTER(IM0_PTR,IM0,(/ ILONG1 /))
+               CALL C_F_POINTER(MCU0_PTR,MCU0,(/ ILONG2 /))
+            ENDIF
+         ENDIF
+      ENDIF
+      IF((.NOT.CYCLIC).AND.(ISCR.GT.0)) THEN
+*        recover (IS,JS) array for surfaces neighbors identification
+         CALL LCMGPD(IPTRK,'IS$MCCG',IS_PTR)
+         CALL LCMGPD(IPTRK,'JS$MCCG',JS_PTR)
+         CALL C_F_POINTER(IS_PTR,IS,(/ NFI-NREG+1 /))
+         CALL C_F_POINTER(JS_PTR,JS,(/ LPS /))
+      ENDIF
+      IF(LPJJAN) THEN
+         CALL LCMGPD(IPTRK,'PJJIND$MCCG',PJJIND_PTR)
+         CALL C_F_POINTER(PJJIND_PTR,PJJIND,(/ NPJJM,2 /))
+      ENDIF
+      IF(LPRISM) THEN
+         CALL LCMGET(IPTRK,'NCODE',NCODE)
+         IF(NCODE(6).EQ.30) THEN
+            IF(NCODE(5).EQ.30) THEN
+*           Z- and Z+ surfaces symmetry
+               SSYM=2
+            ELSE
+*           Z+ symmetry
+               SSYM=1
+            ENDIF
+         ELSE
+            SSYM=0
+         ENDIF
+         NMAX=(INT(FACSYM)+1)*N2MAX*(NZP+2)
+      ELSE
+         NMAX=N2MAX
+      ENDIF
+      ALLOCATE(NOM(N2MAX))
+      ALLOCATE(H(N2MAX),HH(N2MAX),ZMUA(NMU),WZMUA(NMU))
+      IF(LPJJ) THEN
+*     Self-Collision Probabilities
+         ALLOCATE(PJJ(NREG*NPJJM*NGEFF),PJJX(NREG*NPJJM*NGEFF),
+     >            PJJXI(NREG*NPJJM*NGEFF),PJJY(NREG*NPJJM*NGEFF),
+     >            PJJYI(NREG*NPJJM*NGEFF),PJJZ(NREG*NPJJM*NGEFF),
+     >            PJJZI(NREG*NPJJM*NGEFF))
+         IF(LPS.GT.0) THEN
+            ALLOCATE(PSJ(LPS*NGEFF),PSJX(LPS*NGEFF),PSJY(LPS*NGEFF),
+     >               PSJZ(LPS*NGEFF))
+            PSJ(:LPS*NGEFF)=0.0
+            PSJX(:LPS*NGEFF)=0.0
+            PSJY(:LPS*NGEFF)=0.0
+            PSJZ(:LPS*NGEFF)=0.0
+         ENDIF
+         ALLOCATE(PJJD(NREG*NPJJM*NGEFF),PJJDX(NREG*NPJJM*NGEFF),
+     >            PJJDXI(NREG*NPJJM*NGEFF),PJJDY(NREG*NPJJM*NGEFF),
+     >            PJJDYI(NREG*NPJJM*NGEFF),PJJDZ(NREG*NPJJM*NGEFF),
+     >            PJJDZI(NREG*NPJJM*NGEFF))
+         PJJD(:NREG*NPJJM*NGEFF)=0.0D0
+         PJJDX(:NREG*NPJJM*NGEFF)=0.0D0
+         PJJDXI(:NREG*NPJJM*NGEFF)=0.0D0
+         PJJDY(:NREG*NPJJM*NGEFF)=0.0D0
+         PJJDYI(:NREG*NPJJM*NGEFF)=0.0D0
+         PJJDZ(:NREG*NPJJM*NGEFF)=0.0D0
+         PJJDZI(:NREG*NPJJM*NGEFF)=0.0D0
+      ENDIF
+      IF(LACA) THEN
+*     Algebraic Collapsing Acceleration
+         ALLOCATE(XSW((M+1)*NANI*NGEFF))
+         IF(LTMT) ALLOCATE(NOM0(N2MAX),H0(N2MAX))
+         ALLOCATE(PREV(NMAX),NEXT(NMAX))
+         ALLOCATE(CQ(LC*NGEFF),DIAGQ(NLONG*NGEFF),DIAGFR(NLONG),CFR(LC),
+     1   WORK(6*NMAX))
+         IF(PACA.GE.2) THEN
+            ALLOCATE(LUDF(NLONG))
+            IF(LC0.GT.0) ALLOCATE(LUCF(LC0))
+         ENDIF
+         DO II=1,NGEFF
+            JPSYS=KPSYS(II)
+            CALL LCMGET(JPSYS,'DRAGON-S0XSC',XSW((II-1)*(M+1)+1))
+         ENDDO
+         CQ(:LC*NGEFF)=0.0
+         DIAGQ(:NLONG*NGEFF)=0.0
+         DIAGF(:NLONG,:NGEFF)=0.0D0
+         CF(:LC,:NGEFF)=0.0D0
+      ENDIF 
+*
+      NMUA=NMU
+      DO IE=1,NMUA
+         ZMUA(IE)=ZMU(IE)
+         WZMUA(IE)=WZMU(IE)       
+      ENDDO
+      IF((.NOT.LPRISM).AND.(NDIM.EQ.2).AND.LTMT) THEN
+         NMUA=1
+         ZMUI=ZMUA(1)
+         W=WZMUA(1)
+         TEMP=ZMUI*W
+         DO IE=2,NMU
+            ZMUI=ZMUA(IE)
+            WZMUI=WZMUA(IE)
+            W=W+WZMUI
+            TEMP=TEMP+ZMUI*WZMUI
+         ENDDO               
+         ZMUA=TEMP/W
+         WZMUA=W
+      ENDIF
+*---
+* SUMMATION UPON THE TRACKING
+*---
+      REWIND IFTRAK
+      READ(IFTRAK) TEXT4,NCOMNT,NBTR,IFMT
+      DO ICOM=1,NCOMNT
+         READ(IFTRAK)
+      ENDDO
+      READ(IFTRAK) (NITMA,II=1,7),MXSUB,NITMA
+      DO ICOM=1,6
+         READ(IFTRAK)
+      ENDDO
+      CALL KDRCPU(T1)
+      IANGL0=0
+      IPANG=1
+      NMERG=0
+      NTR=0
+      NTRTMT=0
+      NSE=0
+      NSETMT=0
+      LFORC=.FALSE.
+      NTPROC=1
+*
+      ALLOCATE(KANGL(MXSUB))
+      IF(LPRISM) THEN
+*     3D PRISMATIC GEOMETRY CONSTRUCTED FROM A 2D TRACKING
+      N3TR=0
+      N3TRTMT=0
+      N3SE=0
+      N3SETMT=0
+      ALLOCATE(NOM3D(NMAX),H3D(NMAX))
+      IF(LTMT) ALLOCATE(NOM3D0(2*NMAX),H3D0(2*NMAX))
+      CALL LCMSIX(IPTRK,'PROJECTION',1)
+      CALL LCMGPD(IPTRK,'IND2T3',INDREG_PTR)
+      CALL LCMGPD(IPTRK,'ZCOORD',ZZZ_PTR)
+      CALL LCMGPD(IPTRK,'VNORF',VNORF_PTR)
+      CALL LCMGPD(IPTRK,'CMU',CMU_PTR)
+      CALL LCMGPD(IPTRK,'CMUI',CMUI_PTR)
+      CALL LCMGPD(IPTRK,'SMU',SMU_PTR)
+      CALL LCMGPD(IPTRK,'SMUI',SMUI_PTR)
+      CALL LCMGPD(IPTRK,'TMU',TMU_PTR)
+      CALL LCMGPD(IPTRK,'TMUI',TMUI_PTR)
+      CALL C_F_POINTER(INDREG_PTR,INDREG,(/ (N2SOU+N2REG+1)*(NZP+1) /))
+      CALL C_F_POINTER(ZZZ_PTR,ZZZ,(/ NZP+1 /))
+      CALL C_F_POINTER(VNORF_PTR,VNORF,(/ NREG*NANGL*NMU*2 /))
+      CALL C_F_POINTER(CMU_PTR,CMU,(/ NMU /))
+      CALL C_F_POINTER(CMUI_PTR,CMUI,(/ NMU /))
+      CALL C_F_POINTER(SMU_PTR,SMU,(/ NMU /))
+      CALL C_F_POINTER(SMUI_PTR,SMUI,(/ NMU /))
+      CALL C_F_POINTER(TMU_PTR,TMU,(/ NMU /))
+      CALL C_F_POINTER(TMUI_PTR,TMUI,(/ NMU /))
+      CALL LCMSIX(IPTRK,'PROJECTION',2)
+      DO 10 ILINE=1,NBTR
+         READ(IFTRAK) NSUB,NSEG,WEIGHT,(KANGL(II),II=1,NSUB),
+     1        (NOM(II),II=1,NSEG),(H(II),II=1,NSEG)
+         IF(NSUB.GT.MXSUB) CALL XABORT('MCGASM: MXSUB OVERFLOW.')
+         IANGL=KANGL(1)
+         IF(LPJJ) THEN
+*        ----------------------------
+*        Self-Collision Probabilities
+*        ----------------------------
+            NR2SE=NSEG-2
+            HH=0.0
+            DO II=0,NR2SE
+               HH(II+2)=HH(II+1)+H(II+2)
+            ENDDO
+            CALL MCGPTS(SUBPJJ,NFI,NREG,M,NANI,NFUNL,NANGL,NMU,NMOD,
+     1           LPS,NPJJM,NGEFF,IANGL,NSEG,ISGNR,NZON,NOM,IS,JS,
+     2           PJJIND,WEIGHT,CPO,CAZ1,CAZ2,ZMU,WZMU,SIGAL,HH,PSJ,
+     3           PJJD,LPJJAN,NR2SE,NMAX,NZP,N2REG,N2SOU,DELU,INDREG,
+     4           ZZZ,VNORF,CMU,CMUI,SMU,SMUI,TMU,TMUI,SSYM)
+         ENDIF
+         IF(LTMT) THEN
+*        ----------------
+*        Tracking Merging (angle by angle) for ACA
+*        ----------------
+            NTR=NTR+1
+            NSE=NSE+NSEG
+            LFORC=(IPANG.NE.IANGL)
+            IF(LFORC) THEN
+               ITEMP=IANGL
+               IANGL=IPANG
+               IPANG=ITEMP
+            ENDIF
+            IF(ILINE.EQ.NBTR) LFORC=.TRUE.
+            CALL MCGTMT(NMERG,NTRTMT,NSETMT,NSEG,NSEG0,NOM,NOM0,WEIGHT,
+     1           WEIGHT0,H,H0,LFORC,NTPROC)
+            IF(NTPROC.EQ.0) GOTO 10
+         ENDIF
+         IF(LACA) THEN
+*        ---------------------------------
+*        Algebraic Collapsing Acceleration
+*        ---------------------------------
+            NR2SE=NSEG-2
+            HH=0.0
+            DO II=0,NR2SE
+               HH(II+2)=HH(II+1)+H(II+2)
+            ENDDO
+            CALL MCGPTA(NFI,NREG,NLONG,M,NANGL,NMU,LC,NGEFF,
+     1           IANGL,NSEG,NOM,NZONA,IPERM,KM,IM,MCU,PREV,NEXT,
+     2           WEIGHT,ZMU,WZMU,SIGAL,XSW,HH,DIAGQ,CQ,DIAGF,
+     3           CF,WORK,LTMT,SUBDS2,SUBDSP,SUBDSC,NR2SE,NMAX,
+     4           NZP,N2REG,N2SOU,DELU,INDREG,NOM3D,NOM3D0,H3D,
+     5           H3D0,ZZZ,VNORF,CMU,CMUI,SMU,SMUI,TMU,TMUI,N3TR,
+     6           N3TRTMT,N3SE,N3SETMT,NTPROC,SSYM)
+         ENDIF
+ 10   CONTINUE
+      IF(LTMT) THEN
+*     process last integration line for ACA TMT
+         CALL MCGTMT(NMERG,NTRTMT,NSETMT,NSEG,NSEG0,NOM,NOM0,WEIGHT,
+     1        WEIGHT0,H,H0,LFORC,NTPROC) 
+         NR2SE=NSEG-2
+         HH=0.0
+         DO II=0,NR2SE
+            HH(II+2)=HH(II+1)+H(II+2)
+         ENDDO
+         CALL MCGPTA(NFI,NREG,NLONG,M,NANGL,NMU,LC,NGEFF,IANGL,NSEG,
+     1        NOM,NZONA,IPERM,KM,IM,MCU,PREV,NEXT,WEIGHT,ZMU,WZMU,
+     2        SIGAL,XSW,HH,DIAGQ,CQ,DIAGF,CF,WORK,LTMT,SUBDS2,SUBDSP,
+     3        SUBDSC,NR2SE,NMAX,NZP,N2REG,N2SOU,DELU,INDREG,NOM3D,
+     4        NOM3D0,H3D,H3D0,ZZZ,VNORF,CMU,CMUI,SMU,SMUI,TMU,TMUI,
+     5        N3TR,N3TRTMT,N3SE,N3SETMT,NTPROC,SSYM)
+         DEALLOCATE(H3D0,NOM3D0)
+      ENDIF
+      DEALLOCATE(H3D,NOM3D)
+      ELSE
+*     REGULAR 2D OR 3D GEOMETRY
+      ALLOCATE(TRHAR(NMU,NFUNL,NANGL),RHARM(NMU,NFUNL))
+      DO IANGL=1,NANGL
+        IF(NDIM.EQ.2) THEN
+          CALL MOCCHR(NDIM,NANI-1,NFUNL,NMU,CPO(1),CAZ1(IANGL),
+     1                CAZ2(IANGL),RHARM)
+        ELSE
+          XMUANG(1)=REAL(CAZ0(IANGL))
+          CALL MOCCHR(NDIM,NANI-1,NFUNL,1,XMUANG(1),CAZ1(IANGL),
+     1                CAZ2(IANGL),RHARM)
+        ENDIF
+        DO JF=1,NFUNL
+          DO IE=1,NMU
+            TRHAR(IE,JF,IANGL)=ISGNR(1,JF)*RHARM(IE,JF)
+          ENDDO
+        ENDDO
+      ENDDO
+      DEALLOCATE(RHARM)
+      DO 20 ILINE=1,NBTR
+         READ(IFTRAK) NSUB,NSEG,WEIGHT,(KANGL(II),II=1,NSUB),
+     1        (NOM(II),II=1,NSEG),(H(II),II=1,NSEG)
+         IF(NSUB.GT.MXSUB) CALL XABORT('MCGASM: MXSUB OVERFLOW.')
+         IANGL=KANGL(1)
+         DO II=1,NSEG
+            IF(NOM(II).LT.0) THEN
+               NOM(II)=NREG-NOM(II)
+            ELSE IF(NOM(II).EQ.0) THEN
+               NOM(II)=NREG+1
+            ENDIF
+         ENDDO
+         IF(LPJJ) THEN
+*        ----------------------------
+*        Self-Collision Probabilities
+*        ----------------------------
+            IF(ISTRM.LE.2) THEN
+              CALL MCGDS4(SUBPJJ,NSEG,NSUB,NMU,LPS,NFUNL,NANGL,NGEFF,
+     1           WEIGHT,KANGL,TRHAR,H,ZMU,WZMU,NOM,NZON,NFI,NREG,NDIM,
+     2           M,IS,JS,PJJD,PSJ,LPJJAN,NPJJM,PJJIND,SIGAL,1)
+            ELSE IF(ISTRM.EQ.3) THEN
+*             TIBERE model
+              CALL MCGDSD(NSEG,NSUB,NMU,LPS,NFUNL,NANGL,NGEFF,WEIGHT,
+     1           TRHAR,H,ZMU,WZMU,KANGL,NOM,NZON,NFI,NREG,NDIM,M,IS,
+     2           JS,PJJD,PSJ,LPJJAN,NPJJM,PJJIND,SIGAL,1,CAZ1(IANGL),
+     3           CAZ2(IANGL),PJJDX,PJJDY,PJJDZ,PJJDXI,PJJDYI,PJJDZI,
+     4           CAZ0(IANGL),PSJX,PSJY,PSJZ)
+            ENDIF
+         ENDIF
+         IF(LTMT) THEN
+*        ----------------
+*        Tracking Merging (angle by angle) for ACA
+*        ----------------
+            NTR=NTR+1
+            NSE=NSE+NSEG
+            IF(ILINE.EQ.NBTR) LFORC=.TRUE.
+            CALL MCGTMT(NMERG,NTRTMT,NSETMT,NSEG,NSEG0,NOM,NOM0,WEIGHT,
+     1           WEIGHT0,H,H0,LFORC,NTPROC)
+            IF(NTPROC.EQ.0) GOTO 20
+         ENDIF
+         IF(LACA) THEN
+*        ---------------------------------
+*        Algebraic Collapsing Acceleration
+*        ---------------------------------
+            DO II=1,NSEG
+               NOM(II)=IPERM(NOM(II))
+            ENDDO
+            CALL MCGDS1(SUBDS2,SUBDSP,SUBDSC,NSEG,NMUA,NGEFF,WEIGHT,H,
+     1           ZMUA,WZMUA,NOM,NZONA,NLONG,NFI,NDIM,LC,M,KM,IM,MCU,
+     2           DIAGF,DIAGQ,CF,CQ,PREV,NEXT,SIGAL,XSW,WORK)
+         ENDIF
+ 20   CONTINUE
+      DEALLOCATE(TRHAR)
+      IF(LTMT) THEN
+*     process last integration line
+         CALL MCGTMT(NMERG,NTRTMT,NSETMT,NSEG,NSEG0,NOM,NOM0,WEIGHT,
+     1        WEIGHT0,H,H0,LFORC,NTPROC)
+         DO II=1,NSEG
+            NOM(II)=IPERM(NOM(II))
+         ENDDO
+         CALL MCGDS1(SUBDS2,SUBDSP,SUBDSC,NSEG,NMUA,NGEFF,WEIGHT,H,
+     1        ZMUA,WZMUA,NOM,NZONA,NLONG,NFI,NDIM,LC,M,KM,IM,MCU,
+     2        DIAGF,DIAGQ,CF,CQ,PREV,NEXT,SIGAL,XSW,WORK)
+      ENDIF
+      ENDIF
+*
+      IF((LTMT).AND.(IPRINT.GT.1)) THEN
+         WRITE(6,*) 'TRACKING MERGING: FROM TRACKING FILE'
+         WRITE(6,*) '                 ',
+     1              NTR,' TRACKS ->',
+     2              NTRTMT,' EQUIVALENT TRACKS'
+         WRITE(6,*) '                 ',
+     1              NSE,' SEGMENTS ->',
+     2              NSETMT, ' SEGMENTS'
+         IF((.NOT.LPRISM).AND.(NDIM.EQ.2))
+     1    WRITE(6,*) '                 ',
+     2              NMU,' POLAR ANGLES ->',
+     3              ' 1 EQUIVALENT POLAR ANGLE'
+         IF(LPRISM) THEN
+            WRITE(6,*) 'TRACKING MERGING: 3D PRISMATIC EXTENSION'
+            WRITE(6,*) '                 ',
+     1                 N3TR,' TRACKS ->',
+     2                 N3TRTMT,' EQUIVALENT TRACKS'
+            WRITE(6,*) '                 ',
+     1                 N3SE,' SEGMENTS ->',
+     2                 N3SETMT, ' SEGMENTS'
+         ENDIF
+      ENDIF
+      CALL KDRCPU(T2)
+      IF(IPRINT.GT.0) THEN
+         WRITE(6,100) 'SUMMATION UPON THE TRACKING FOR ACA/SCR ',(T2-T1)
+      ENDIF
+*---
+*  FOR ACA: (IF PACA.GE.2)
+*  CALCULATION OF ILU0 PRECONDITIONER FOR BICGSTAB ITERATIONS
+*---
+      IF(LACA) THEN
+         DO II=1,NGEFF
+           IG=NGIND(II)
+           JPSYS=KPSYS(II)
+           IF(IPRINT.GT.2) WRITE(6,200) 'GROUP',IG
+           CALL MCGDS3(NLONG,PACA,M,SIGAL(0,II),
+     1          XSW((II-1)*(M+1)+1),VA,NZONA,LC,MCU,IM,JU,LC0,IM0,
+     2          MCU0,DIAGF(1,II),CF(1,II),DIAGQ((II-1)*NLONG+1),DIAGFR,
+     3          CFR,LUDF,LUCF)
+*     
+           IF(IPRINT.GT.3) THEN
+              CALL PRINAM('DIAGF ',DIAGFR(1),NLONG)
+              CALL PRINAM('CF    ',CFR(1),LC)
+           ENDIF
+           CALL LCMPUT(JPSYS,'DIAGF$MCCG',NLONG,2,DIAGFR)
+           CALL LCMPUT(JPSYS,'CF$MCCG',LC,2,CFR)
+           IF(PACA.GE.2) THEN
+              CALL LCMPUT(JPSYS,'ILUDF$MCCG',NLONG,2,LUDF)
+              IF(LC0.GT.0) CALL LCMPUT(JPSYS,'ILUCF$MCCG',LC0,2,LUCF)
+           ENDIF
+           IF(IPRINT.GT.3) THEN
+              CALL PRINAM('DIAGQ ',DIAGQ((II-1)*NLONG+1),NLONG)
+              CALL PRINAM('CQ    ',CQ((II-1)*LC+1),LC)
+           ENDIF
+           CALL LCMPUT(JPSYS,'CQ$MCCG',LC,2,CQ((II-1)*LC+1))
+           CALL LCMPUT(JPSYS,'DIAGQ$MCCG',NLONG,2,DIAGQ((II-1)*NLONG+1))
+         ENDDO
+         CALL KDRCPU(T3)
+         IF(IPRINT.GT.0) THEN
+            WRITE(6,100) 'CALCULATION OF ACA PRECONDITIONER       ',
+     1                   (T3-T2)
+         ENDIF
+*
+         IF(PACA.GE.2) THEN
+            IF(LC0.GT.0) DEALLOCATE(LUCF)
+            DEALLOCATE(LUDF)
+         ENDIF
+         DEALLOCATE(WORK,CFR,DIAGFR,DIAGQ,CQ)
+         DEALLOCATE(NEXT,PREV)
+         IF(LTMT) DEALLOCATE(H0,NOM0)
+         DEALLOCATE(XSW)
+      ENDIF
+*----
+*  FOR SCR/STIS:
+*  VECTORS NORMALIZATION
+*----
+      IF(LPJJ) THEN
+         CALL MCGDS6(NGEFF,NPJJM,NREG,PJJD,V,PJJ)
+         CALL MCGDS6(NGEFF,NPJJM,NREG,PJJDX,V,PJJX)
+         CALL MCGDS6(NGEFF,NPJJM,NREG,PJJDY,V,PJJY)
+         CALL MCGDS6(NGEFF,NPJJM,NREG,PJJDZ,V,PJJZ)
+         CALL MCGDS6(NGEFF,NPJJM,NREG,PJJDXI,V,PJJXI)
+         CALL MCGDS6(NGEFF,NPJJM,NREG,PJJDYI,V,PJJYI)
+         CALL MCGDS6(NGEFF,NPJJM,NREG,PJJDZI,V,PJJZI)
+         DEALLOCATE(PJJD,PJJDX,PJJDY,PJJDZ,PJJDXI,PJJDYI,PJJDZI)
+*
+         DO II=1,NGEFF
+            JPSYS=KPSYS(II)
+            IF(IPRINT.GT.3) THEN
+               CALL PRINAM('PJJ   ',PJJ((II-1)*NREG*NPJJM+1),NREG*NPJJM)
+               IF(LPS.GT.0) CALL PRINAM('PSJ   ',PSJ((II-1)*LPS+1),LPS)
+            ENDIF
+            CALL LCMPUT(JPSYS,'PJJ$MCCG',NREG*NPJJM,2,
+     1                  PJJ((II-1)*NREG*NPJJM+1))
+            CALL LCMPUT(JPSYS,'PJJX$MCCG',NREG*NPJJM,2,
+     1                  PJJX((II-1)*NREG*NPJJM+1))
+            CALL LCMPUT(JPSYS,'PJJY$MCCG',NREG*NPJJM,2,
+     1                  PJJY((II-1)*NREG*NPJJM+1))
+            CALL LCMPUT(JPSYS,'PJJZ$MCCG',NREG*NPJJM,2,
+     1                  PJJZ((II-1)*NREG*NPJJM+1))
+            CALL LCMPUT(JPSYS,'PJJXI$MCCG',NREG*NPJJM,2,
+     1                  PJJXI((II-1)*NREG*NPJJM+1))
+            CALL LCMPUT(JPSYS,'PJJYI$MCCG',NREG*NPJJM,2,
+     1                  PJJYI((II-1)*NREG*NPJJM+1))
+            CALL LCMPUT(JPSYS,'PJJZI$MCCG',NREG*NPJJM,2,
+     1                  PJJZI((II-1)*NREG*NPJJM+1))
+            IF(LPS.GT.0) THEN
+               CALL LCMPUT(JPSYS,'PSJ$MCCG',LPS,2,PSJ((II-1)*LPS+1))
+               CALL LCMPUT(JPSYS,'PSJX$MCCG',LPS,2,PSJX((II-1)*LPS+1))
+               CALL LCMPUT(JPSYS,'PSJY$MCCG',LPS,2,PSJY((II-1)*LPS+1))
+               CALL LCMPUT(JPSYS,'PSJZ$MCCG',LPS,2,PSJZ((II-1)*LPS+1))
+            ENDIF
+         ENDDO
+*
+         IF(LPS.GT.0) DEALLOCATE(PSJ,PSJX,PSJY,PSJZ)
+         DEALLOCATE(PJJX,PJJXI,PJJY,PJJYI,PJJZ,PJJZI)
+         DEALLOCATE(PJJ)
+      ENDIF
+*
+      DEALLOCATE(KANGL,WZMUA,ZMUA,HH,H,NOM)
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(CF,DIAGF)
+      RETURN
+*
+ 100  FORMAT('   -->>TIME SPENT IN: ',A40,':',F13.3,' s.')
+ 200  FORMAT(1X,A6,1X,I4)
+      END