summaryrefslogtreecommitdiff
path: root/Dragon/src/MCGPTS.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/MCGPTS.f')
-rw-r--r--Dragon/src/MCGPTS.f332
1 files changed, 332 insertions, 0 deletions
diff --git a/Dragon/src/MCGPTS.f b/Dragon/src/MCGPTS.f
new file mode 100644
index 0000000..7de7a98
--- /dev/null
+++ b/Dragon/src/MCGPTS.f
@@ -0,0 +1,332 @@
+*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
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 21:25:08 +0000