*DECK MCGPTS SUBROUTINE MCGPTS(SUBPJJ,NFI,NREG,M,NANI,NFUNL,NANGL,NMU,NMOD,LPS, 1 NPJJM,NGEFF,IANGL,NSEG,ISGNR,NZON,NOM2D,IS,JS, 2 PJJIND,W2D,XMU,CAZ1,CAZ2,ZMU,WZMU,SIGAL,T2D,PSJ, 3 PJJD,LPJJAN,NR2D,NMAX,NZP,N2REG,N2SOU,DELU, 4 INDREG,Z,VNORF,CMU,CMUI,SMU,SMUI,TMU,TMUI,SSYM) * *----------------------------------------------------------------------- * *Purpose: * Calculation of the PJJ and PSJ (3D prismatic extended tracking). * *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): R. Le Tellier * *Parameters: input * SUBPJJ PJJ calculation subroutine. * 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. * M number of material mixtures. * NANI number of Legendre orders. * NFUNL number of moments of the flux (in 2D: NFUNL=NANI*(NANI+1)/2). * NANGL number of tracking angles in the plane. * NMU order of the polar quadrature in 2D / 1 in 3D. * NMOD first dimension of ISGNR. * LPS dimension of JS. * NPJJM number of pjj modes to store for STIS option. * NGEFF number of groups to process. * IANGL direction index of this track. * NSEG number of elements in the current track. * ISGNR spherical harmonic signs. * NZON index-number of the mixture type assigned to each volume. * NOM2D vector containing the region number of the different segments * of this 2D track. * 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 LPJJAN option. * W2D track weight. * XMU polar angle cosines. * CAZ1 first cosines of the different tracking azimuthal angles. * CAZ2 second cosines of the different tracking azimuthal angles. * ZMU polar quadrature set. * WZMU polar quadrature set. * SIGAL albedos and total cross sections array. * T2D vector containing the local coordinates of the segments * boundaries for this 2D track. * LPJJAN anisotropic scattering flag. * NR2D number of segments corresponding to regions for this 2D track. * NMAX maximum number of segments for the 3D tracks. * NZP number of z-planes. * N2SOU number of external surfaces in the 2D tracking. * N2REG number of regions in the 2D tracking. * DELU input track spacing for 3D track reconstruction. * INDREG region/surface index to go from the 2D to the 3D geometry. * Z z-plan coordinates. * VNORF normalization factors per angle. * CMU polar angle cosines. * CMUI inverse of polar angle cosines. * SMU polar angle sines. * SMUI inverse of polar angle sines. * TMU polar angle tangents. * TMUI inverse of polar angle tangents. * SSYM symmetry flag. * *Parameters: input/output * PJJD collision probabilities. * PSJ leakage probabilities. * *----------------------------------------------------------------------- * IMPLICIT NONE *---- * SUBROUTINE ARGUMENTS *---- INTEGER NFI,NREG,M,NANI,NFUNL,NANGL,NMU,NMOD,LPS,NPJJM,NGEFF, 1 IANGL,NSEG,ISGNR(NMOD,NFUNL),NZON(NFI),NOM2D(NSEG),JS(LPS), 2 IS(NFI-NREG+1),PJJIND(NPJJM,2),NR2D,NMAX,NZP,N2REG,N2SOU, 3 INDREG(-N2SOU:N2REG,0:NZP+1),SSYM REAL XMU(NMU),ZMU(NMU),WZMU(NMU),SIGAL(-6:M,NGEFF),PSJ(LPS,NGEFF), 1 DELU,Z(0:NZP) DOUBLE PRECISION PJJD(NREG,NPJJM,NGEFF),VNORF(NREG,NANGL,NMU,2), 1 CMU(NMU),CMUI(NMU),SMU(NMU),SMUI(NMU),TMU(NMU),TMUI(NMU), 2 W2D,T2D(0:NR2D),CAZ1(NANGL),CAZ2(NANGL) LOGICAL LPJJAN EXTERNAL SUBPJJ *---- * LOCAL VARIABLES *---- INTEGER MODUR,MODDR PARAMETER(MODUR=1,MODDR=5) INTEGER JF,IE,IMU,NBTR,KST,IST,ILINE,I,I1,I2,K,N3D,II,TIN,N3DP, 1 NSUB,NANGL0,KANGL(1) DOUBLE PRECISION CPO,CPOI,SPO,SPOI,TPO,TPOI,LTOT,DELTE,DELZE,T, 1 Z1,Z2,TP,Z1P,WPO,W3D,W3DPO *---- * ALLOCATABLE ARRAYS *---- INTEGER, ALLOCATABLE, DIMENSION(:) :: NOM3D REAL, ALLOCATABLE, DIMENSION(:,:) :: RHARM,TRHARP,TRHARM DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: H3D *---- * SCRATCH STORAGE ALLOCATION *---- ALLOCATE(NOM3D(NMAX),RHARM(NMU,NFUNL),TRHARP(NMU,NFUNL), 1 TRHARM(NMU,NFUNL),H3D(NMAX)) * NSUB=1 NANGL0=1 KANGL(1)=1 IF (LPJJAN) THEN IF(MODDR.GT.NMOD) CALL XABORT('MCGPTS: NMOD OVERFLOW') CALL MOCCHR(3,NANI-1,NFUNL,NMU,XMU,CAZ1(IANGL),CAZ2(IANGL), 1 RHARM) DO JF=1,NFUNL DO IE=1,NMU TRHARP(IE,JF)=ISGNR(MODUR,JF)*RHARM(IE,JF) TRHARM(IE,JF)=ISGNR(MODDR,JF)*RHARM(IE,JF) ENDDO ENDDO ENDIF * DO IMU=1,NMU CPO=CMU(IMU) CPOI=CMUI(IMU) SPO=SMU(IMU) SPOI=SMUI(IMU) TPO=TMU(IMU) TPOI=TMUI(IMU) WPO=WZMU(IMU) IF (SSYM.EQ.2) GOTO 15 *--- * CONSTRUCT THE 3D TRACKS WHICH ENTER THE GEOMETRY THROUGH A BOTTOM/TOP SURFACE *--- * length of the spatial integration interval LTOT=T2D(NR2D)*CPO * number of 3D tracks generated for this x-y track and this polar direction NBTR=INT(LTOT/DELU)+1 * effective track spacing in T DELTE=T2D(NR2D)/DBLE(NBTR) W3DPO=W2D*DELTE*CPO W3D=WPO*W3DPO T=-0.5D0*DELTE KST=1 DO 10 ILINE=1,NBTR T=T+DELTE TP=T DO WHILE (T2D(KST).LT.T) KST=KST+1 ENDDO K=KST * --- * positive polar sine track * --- I1=1 Z1=Z(I1-1) TIN=0 N3D=1 NOM3D(N3D)=NREG-INDREG(NOM2D(K+1),0) H3D(N3D)=0.5D0 CALL MCGPT1(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I1,K,Z1,T, 1 TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D) DO II=2,N3D-1 H3D(II)=H3D(II)*VNORF(NOM3D(II),IANGL,IMU,1) ENDDO IF (SSYM.EQ.1) THEN * the top boundary condition is a surface symmetry IF (TIN.EQ.0) THEN * this track has encountered the top boundary -> it is reflected N3DP=N3D N3D=N3D-1 I1=I1-1 CALL MCGPT2(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I1,K,Z1, 1 T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D) DO II=N3DP,N3D-1 H3D(II)=H3D(II)*VNORF(NOM3D(II),IANGL,IMU,2) ENDDO ENDIF ENDIF NOM3D(N3D)=NREG-NOM3D(N3D) CALL MCGDS4(SUBPJJ,N3D,NSUB,NMU,LPS,NFUNL,NANGL0,NGEFF,W3D, 1 KANGL,TRHARP,H3D,ZMU,WZMU,NOM3D,NZON,NFI,NREG,3,M,IS,JS, 2 PJJD,PSJ,LPJJAN,NPJJM,PJJIND,SIGAL,IMU,1) T=TP IF (SSYM.EQ.1) GOTO 10 K=KST * --- * negative polar sine track * --- I2=NZP Z2=Z(I2) TIN=0 N3D=1 NOM3D(N3D)=NREG-INDREG(NOM2D(K+1),NZP+1) H3D(N3D)=0.5D0 CALL MCGPT2(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I2,K,Z2,T, 1 TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D) NOM3D(N3D)=NREG-NOM3D(N3D) DO II=2,N3D-1 H3D(II)=H3D(II)*VNORF(NOM3D(II),IANGL,IMU,2) ENDDO CALL MCGDS4(SUBPJJ,N3D,NSUB,NMU,LPS,NFUNL,NANGL0,NGEFF,W3D, 1 KANGL,TRHARM,H3D,ZMU,WZMU,NOM3D,NZON,NFI,NREG,3,M,IS,JS, 2 PJJD,PSJ,LPJJAN,NPJJM,PJJIND,SIGAL,IMU,1) * --- T=TP 10 CONTINUE *--- * CONSTRUCT THE 3D TRACKS WHICH ENTER THE GEOMETRY THROUGH A LATERAL SURFACE *--- * length of the spatial integration interval 15 LTOT=Z(NZP)*SPO ! LTOT=(Z(NZP)-Z(0))*SPO with Z(0)=0.0 * number of 3D tracks generated for this x-y track and this polar direction NBTR=INT(LTOT/DELU)+1 * effective track spacing in Z DELZE=Z(NZP)/DBLE(NBTR) ! DELZE=(Z(NZP)-Z(0))/DBLE(NBTR) with Z(0)=0.0 W3DPO=W2D*DELZE*SPO W3D=WPO*W3DPO Z1=-0.5D0*DELZE ! Z1=Z(0)-0.5D0*DELZE with Z(0)=0.0 IST=1 DO 20 ILINE=1,NBTR Z1=Z1+DELZE Z1P=Z1 DO WHILE (Z(IST).LT.Z1) IST=IST+1 ENDDO I=IST * --- * positive polar sine track * --- K=1 T=T2D(K-1) TIN=1 N3D=1 N3DP=2 NOM3D(N3D)=NREG-INDREG(NOM2D(1),IST) H3D(N3D)=0.5D0 21 CONTINUE CALL MCGPT1(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I,K,Z1,T, 1 TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D) DO II=N3DP,N3D-1 H3D(II)=H3D(II)*VNORF(NOM3D(II),IANGL,IMU,1) ENDDO IF (SSYM.GT.0) THEN * the top boundary condition is a surface symmetry IF (TIN.EQ.0) THEN * this track has encountered the top boundary -> it is reflected N3DP=N3D N3D=N3D-1 I=I-1 CALL MCGPT2(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I,K,Z1, 1 T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D) DO II=N3DP,N3D-1 H3D(II)=H3D(II)*VNORF(NOM3D(II),IANGL,IMU,2) ENDDO IF ((SSYM.EQ.2).AND.(TIN.EQ.0)) THEN * the bottom boundary is a surface symmetry * this track has encountered the bottom boundary -> it is reflected N3DP=N3D N3D=N3D-1 I=I+1 GOTO 21 ENDIF ENDIF ENDIF NOM3D(N3D)=NREG-NOM3D(N3D) CALL MCGDS4(SUBPJJ,N3D,NSUB,NMU,LPS,NFUNL,NANGL0,NGEFF,W3D, 1 KANGL,TRHARP,H3D,ZMU,WZMU,NOM3D,NZON,NFI,NREG,3,M,IS,JS, 2 PJJD,PSJ,LPJJAN,NPJJM,PJJIND,SIGAL,IMU,1) Z1=Z1P I=IST * --- * negative polar sine track * --- K=1 T=T2D(K-1) TIN=1 N3D=1 N3DP=2 NOM3D(N3D)=NREG-INDREG(NOM2D(1),IST) H3D(N3D)=0.5D0 22 CONTINUE CALL MCGPT2(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I,K,Z1,T, 1 TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D) DO II=N3DP,N3D-1 H3D(II)=H3D(II)*VNORF(NOM3D(II),IANGL,IMU,2) ENDDO IF (SSYM.EQ.2) THEN * the bottom boundary is a surface symmetry IF (TIN.EQ.0) THEN * this track has encountered the bottom boundary -> it is reflected N3DP=N3D N3D=N3D-1 I=I+1 CALL MCGPT1(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I,K,Z1, 1 T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D) DO II=N3DP,N3D-1 H3D(II)=H3D(II)*VNORF(NOM3D(II),IANGL,IMU,1) ENDDO IF (TIN.EQ.0) THEN * the top boundary is a surface symmetry * this track has encountered the top boundary -> it is reflected N3DP=N3D N3D=N3D-1 I=I-1 GOTO 22 ENDIF ENDIF ENDIF NOM3D(N3D)=NREG-NOM3D(N3D) CALL MCGDS4(SUBPJJ,N3D,NSUB,NMU,LPS,NFUNL,NANGL0,NGEFF,W3D, 1 KANGL,TRHARM,H3D,ZMU,WZMU,NOM3D,NZON,NFI,NREG,3,M,IS,JS, 2 PJJD,PSJ,LPJJAN,NPJJM,PJJIND,SIGAL,IMU,1) * --- Z1=Z1P 20 CONTINUE * ENDDO DEALLOCATE(H3D,TRHARM,TRHARP,RHARM,NOM3D) RETURN END