*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