summaryrefslogtreecommitdiff
path: root/Dragon/src/MCGPTF.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/MCGPTF.f')
-rw-r--r--Dragon/src/MCGPTF.f663
1 files changed, 663 insertions, 0 deletions
diff --git a/Dragon/src/MCGPTF.f b/Dragon/src/MCGPTF.f
new file mode 100644
index 0000000..5e88feb
--- /dev/null
+++ b/Dragon/src/MCGPTF.f
@@ -0,0 +1,663 @@
+*DECK MCGPTF
+ SUBROUTINE MCGPTF(SUBFFI,SUBFFA,SUBSCH,IFTRAK,N2BTR,N2MAX,
+ 1 KPN,K,NREG,M,NGEFF,NANGL,NMU,NANI,NFUNL,NMOD,
+ 2 KEYFLX,KEYCUR,NZON,NCONV,CAZ1,CAZ2,XMU,WZMU,PHI,
+ 3 S,SIGAL,ISGNR,NMAX,NZP,N2REG,N2SOU,DELU,INDREG,
+ 4 Z,VNORF,CMU,CMUI,SMU,SMUI,TMU,TMUI,SSYM,IDIR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Flux integration upon the tracking (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
+* SUBFFI isotropic flux integration subroutine.
+* SUBFFA anisotropic flux integration subroutine.
+* SUBSCH track coefficients calculation subroutine.
+* IFTRAK tracking file unit number.
+* N2BTR total number of 2D tracking lines.
+* N2MAX maximum number of elements in a 2D track.
+* KPN total number of unknowns in vectors PHI.
+* K total number of volumes for which specific values
+* of the neutron flux and reactions rates are required.
+* NREG number of volumes.
+* M number of material mixtures.
+* NGEFF number of groups to process.
+* NANGL number of tracking angles in the tracking file.
+* NMU order of the polar quadrature in 2D / 1 in 3D.
+* NANI scattering anisotropy (=1 for isotropic scattering).
+* NFUNL number of moments of the flux (in 2D: NFUNL=NANI*(NANI+1)/2).
+* NMOD first dimension of ISGNR.
+* KEYFLX position of flux elements in PHI vector.
+* KEYCUR position of current elements in PHI vector.
+* NZON index-number of the mixture type assigned to each volume.
+* NCONV logical array of convergence status for each group (.TRUE.
+* not converged).
+* CAZ1 first cosines of the different tracking azimuthal angles.
+* CAZ2 second cosines of the different tracking azimuthal angles.
+* XMU cosines of the different tracking polar angles.
+* (polar quadrature in 2D / tracking angles in 3D).
+* WZMU polar quadrature set in 2D.
+* S total source vector components.
+* SIGAL total cross-section and albedo array.
+* ISGNR spherical harmonic signs.
+* NMAX maximum number of segments for the 3D tracks.
+* NZP number of z-plans.
+* 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.
+* IDIR direction of fundamental current for TIBERE with MoC
+* (=0,1,2,3).
+*
+*Parameters: input/output
+* PHI vector containing the zonal scalar flux.
+*
+*-----------------------------------------------------------------------
+*
+ IMPLICIT NONE
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER NGEFF,K,KPN,M,N2MAX,NMU,NZON(K),NANI,NFUNL,NMOD,
+ 1 NREG,KEYFLX(NREG,NFUNL),KEYCUR(K-NREG),IFTRAK,N2BTR,NANGL,
+ 2 ISGNR(NMOD,NFUNL),NMAX,NZP,N2REG,N2SOU,
+ 3 INDREG(-N2SOU:N2REG,0:NZP+1),SSYM,IDIR
+ REAL WZMU(NMU),SIGAL(-6:M,NGEFF),XMU(NMU),DELU,Z(0:NZP)
+ DOUBLE PRECISION CAZ1(NANGL),CAZ2(NANGL),PHI(KPN,NGEFF),
+ 1 S(KPN,NGEFF),VNORF(NREG,NANGL,NMU,2),CMU(NMU),CMUI(NMU),
+ 2 SMU(NMU),SMUI(NMU),TMU(NMU),TMUI(NMU)
+ LOGICAL NCONV(NGEFF)
+ EXTERNAL SUBFFI,SUBFFA,SUBSCH
+*----
+* LOCAL VARIABLES
+*----
+ INTEGER MODUR,MODDL,MODDR,MODUL
+ PARAMETER(MODUR=1,MODDL=8,MODDR=5,MODUL=4)
+ INTEGER II,I2LIN,IANG,N2SEG,NR2D,NBTR,KST,IST,ILINE,N3D,I,I1,I2,
+ 1 IMU,IANG0,NOMP,INDP,NOMM,INDM,NOMI,JF,IND,TIN,N3DP,NSUB
+ DOUBLE PRECISION W2D,Q0,Q1,CPO,CPOI,SPO,SPOI,TPO,TPOI,LTOT,DELTE,
+ 1 DELZE,T,Z1,Z2,TP,Z1P,W3DPO,WPO,W3D,OMEGAX,OMEGAY,OMEGAZ
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: NOM2D,NOM3D
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: RHARM
+ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: TRHAR
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: H2D,H3D,T2D,B
+ DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:) :: STOT
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(NOM2D(N2MAX),H2D(N2MAX),NOM3D(NMAX),H3D(NMAX),B(4*NMAX),
+ 1 T2D(0:N2MAX-1))
+*
+ IF(NANI.EQ.1) THEN
+*---
+* ISOTROPIC SCATTERING
+*---
+ DO I2LIN=1,N2BTR
+ READ(IFTRAK) NSUB,N2SEG,W2D,IANG,(NOM2D(I),I=1,N2SEG),
+ 1 (H2D(I),I=1,N2SEG)
+ IF(NSUB.NE.1) CALL XABORT('MCGPTF: NSUB.NE.1.')
+ NR2D=N2SEG-2
+ T2D(0)=0.0
+ DO II=1,NR2D
+ T2D(II)=T2D(II-1)+H2D(II+1)
+ ENDDO
+ 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)=INDREG(NOM2D(K+1),0)
+ H3D(N3D)=0.5D0
+ CALL MCGPT1(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I1,K,Z1,
+ 1 T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D)
+ DO II=2,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,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,
+ 1 Z1,T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D)
+ DO II=N3DP,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,IMU,2)
+ ENDDO
+ ENDIF
+ ENDIF
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+* MCGFFIR: 'Source Term Isolation' Strategy turned on
+* MCGFFIS: 'Source Term Isolation' Strategy turned off
+* MCGFFIT: 'MOCC/MCI' Iterative Strategy
+ OMEGAX=0.0D0
+ OMEGAY=0.0D0
+ OMEGAZ=0.0D0
+ IDIR=0
+ CALL SUBFFI(SUBSCH,K,KPN,M,N3D,H3D,NOM3D,NZON,
+ 1 SIGAL(0,II),S(1,II),NREG,KEYFLX,KEYCUR,PHI(1,II),
+ 2 B,W3D,OMEGAX,OMEGAY,OMEGAZ,IDIR)
+ ENDIF
+ ENDDO
+ 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)=INDREG(NOM2D(K+1),NZP+1)
+ H3D(N3D)=0.5D0
+ CALL MCGPT2(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I2,K,Z2,
+ 1 T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D)
+ DO II=2,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,IMU,2)
+ ENDDO
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+* MCGFFIR: 'Source Term Isolation' Strategy turned on
+* MCGFFIS: 'Source Term Isolation' Strategy turned off
+* MCGFFIT: 'MOCC/MCI' Iterative Strategy
+ OMEGAX=0.0D0
+ OMEGAY=0.0D0
+ OMEGAZ=0.0D0
+ IDIR=0
+ CALL SUBFFI(SUBSCH,K,KPN,M,N3D,H3D,NOM3D,NZON,
+ 1 SIGAL(0,II),S(1,II),NREG,KEYFLX,KEYCUR,PHI(1,II),
+ 2 B,W3D,OMEGAX,OMEGAY,OMEGAZ,IDIR)
+ ENDIF
+ ENDDO
+* ---
+ 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)=INDREG(NOM2D(1),IST)
+ H3D(N3D)=0.5D0
+ 21 CONTINUE
+ 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),IANG,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,
+ 1 Z1,T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D)
+ DO II=N3DP,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,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
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+* MCGFFIR: 'Source Term Isolation' Strategy turned on
+* MCGFFIS: 'Source Term Isolation' Strategy turned off
+* MCGFFIT: 'MOCC/MCI' Iterative Strategy
+ OMEGAX=0.0D0
+ OMEGAY=0.0D0
+ OMEGAZ=0.0D0
+ IDIR=0
+ CALL SUBFFI(SUBSCH,K,KPN,M,N3D,H3D,NOM3D,NZON,
+ 1 SIGAL(0,II),S(1,II),NREG,KEYFLX,KEYCUR,PHI(1,II),
+ 2 B,W3D,OMEGAX,OMEGAY,OMEGAZ,IDIR)
+ ENDIF
+ ENDDO
+ Z1=Z1P
+ I=IST
+* ---
+* negative polar sine track
+* ---
+ K=1
+ T=T2D(K-1)
+ TIN=1
+ N3D=1
+ N3DP=2
+ NOM3D(N3D)=INDREG(NOM2D(1),IST)
+ H3D(N3D)=0.5D0
+ 22 CONTINUE
+ 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),IANG,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,
+ 1 Z1,T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D)
+ DO II=N3DP,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,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
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+* MCGFFIR: 'Source Term Isolation' Strategy turned on
+* MCGFFIS: 'Source Term Isolation' Strategy turned off
+* MCGFFIT: 'MOCC/MCI' Iterative Strategy
+ OMEGAX=0.0D0
+ OMEGAY=0.0D0
+ OMEGAZ=0.0D0
+ IDIR=0
+ CALL SUBFFI(SUBSCH,K,KPN,M,N3D,H3D,NOM3D,NZON,
+ 1 SIGAL(0,II),S(1,II),NREG,KEYFLX,KEYCUR,PHI(1,II),
+ 2 B,W3D,OMEGAX,OMEGAY,OMEGAZ,IDIR)
+ ENDIF
+ ENDDO
+* ---
+ Z1=Z1P
+ 20 CONTINUE
+ ENDDO
+ ENDDO
+ ELSE
+*---
+* ANISOTROPIC SCATTERING
+*---
+ ALLOCATE(STOT(NMAX,NMU,NGEFF,2))
+ STOT(:NMAX,:NMU,:NGEFF,:2)=0.0D0
+ ALLOCATE(RHARM(NMU,NFUNL),TRHAR(NMU,NFUNL,4))
+ IANG0=0
+ DO I2LIN=1,N2BTR
+ READ(IFTRAK) NSUB,N2SEG,W2D,IANG,(NOM2D(I),I=1,N2SEG),
+ 1 (H2D(I),I=1,N2SEG)
+ IF(NSUB.NE.1) CALL XABORT('MCGPTF: NSUB.NE.1.')
+ NR2D=N2SEG-2
+ T2D(0)=0.0
+ DO II=1,NR2D
+ T2D(II)=T2D(II-1)+H2D(II+1)
+ ENDDO
+ IF(IANG.NE.IANG0) THEN
+ IANG0=IANG
+ CALL MOCCHR(3,NANI-1,NFUNL,NMU,XMU,CAZ1(IANG),CAZ2(IANG),
+ 1 RHARM)
+ DO 26 JF=1,NFUNL
+ DO 25 IMU=1,NMU
+ ! positive polar sine track
+ TRHAR(IMU,JF,1)=ISGNR(MODUR,JF)*RHARM(IMU,JF)
+ TRHAR(IMU,JF,2)=ISGNR(MODDL,JF)*RHARM(IMU,JF)
+ ! negative polar sine track
+ TRHAR(IMU,JF,3)=ISGNR(MODDR,JF)*RHARM(IMU,JF)
+ TRHAR(IMU,JF,4)=ISGNR(MODUL,JF)*RHARM(IMU,JF)
+ 25 CONTINUE
+ 26 CONTINUE
+ 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)
+*---
+* 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 30 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)=INDREG(NOM2D(K+1),0)
+ H3D(N3D)=0.5D0
+ CALL MCGPT1(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I1,K,Z1,
+ 1 T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D)
+ DO II=2,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,IMU,1)
+ ENDDO
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+* incoming flux in + direction
+ NOMP=NOM3D(1)
+ INDP=KEYCUR(-NOMP)
+* incoming flux in - direction
+ NOMM=NOM3D(N3D)
+ INDM=KEYCUR(-NOMM)
+ STOT(1,IMU,II,1)=W3D*S(INDP,II)
+ STOT(N3D,IMU,II,2)=W3D*S(INDM,II)
+* regional sources
+ DO I=2,N3D-1
+ NOMI=NOM3D(I)
+ Q0=0.0D0
+ Q1=0.0D0
+ DO JF=1,NFUNL
+ IND=KEYFLX(NOMI,JF)
+ Q0=Q0+S(IND,II)*TRHAR(IMU,JF,1)
+ Q1=Q1+S(IND,II)*TRHAR(IMU,JF,2)
+ ENDDO
+ STOT(I,IMU,II,1)=W3D*Q0
+ STOT(I,IMU,II,2)=W3D*Q1
+ ENDDO
+* MCGFFAR: 'Source Term Isolation' Strategy turned on
+* MCGFFAS: 'Source Term Isolation' Strategy turned off
+* MCGFFAT: 'MOCC/MCI' Iterative Strategy
+ CALL SUBFFA(SUBSCH,K,KPN,M,N3D,H3D,NOM3D,NZON,
+ 1 SIGAL(0,II),STOT(1,IMU,II,1),STOT(1,IMU,II,2),
+ 2 NREG,NMU,NANI,NFUNL,TRHAR(1,1,1),KEYFLX,KEYCUR,
+ 3 IMU,PHI(1,II),B)
+ ENDIF
+ ENDDO
+ T=TP
+ K=KST
+* ---
+* negative polar sine track
+* ---
+ I2=NZP
+ Z2=Z(I2)
+ TIN=0
+ N3D=1
+ NOM3D(N3D)=INDREG(NOM2D(K+1),NZP+1)
+ H3D(N3D)=0.5D0
+ CALL MCGPT2(N2SOU,N2REG,NZP,NR2D,INDREG,Z,NOM2D,T2D,I2,K,Z2,
+ 1 T,TIN,CPOI,SPOI,TPO,TPOI,N3D,NOM3D,H3D)
+ DO II=2,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,IMU,2)
+ ENDDO
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+* incoming flux in + direction
+ NOMP=NOM3D(1)
+ INDP=KEYCUR(-NOMP)
+* incoming flux in - direction
+ NOMM=NOM3D(N3D)
+ INDM=KEYCUR(-NOMM)
+ STOT(1,IMU,II,1)=W3D*S(INDP,II)
+ STOT(N3D,IMU,II,2)=W3D*S(INDM,II)
+* regional sources
+ DO I=2,N3D-1
+ NOMI=NOM3D(I)
+ Q0=0.0D0
+ Q1=0.0D0
+ DO JF=1,NFUNL
+ IND=KEYFLX(NOMI,JF)
+ Q0=Q0+S(IND,II)*TRHAR(IMU,JF,3)
+ Q1=Q1+S(IND,II)*TRHAR(IMU,JF,4)
+ ENDDO
+ STOT(I,IMU,II,1)=W3D*Q0
+ STOT(I,IMU,II,2)=W3D*Q1
+ ENDDO
+* MCGFFAR: 'Source Term Isolation' Strategy turned on
+* MCGFFAS: 'Source Term Isolation' Strategy turned off
+* MCGFFAT: 'MOCC/MCI' Iterative Strategy
+ CALL SUBFFA(SUBSCH,K,KPN,M,N3D,H3D,NOM3D,NZON,
+ 1 SIGAL(0,II),STOT(1,IMU,II,1),STOT(1,IMU,II,2),
+ 2 NREG,NMU,NANI,NFUNL,TRHAR(1,1,3),KEYFLX,KEYCUR,
+ 3 IMU,PHI(1,II),B)
+ ENDIF
+ ENDDO
+* ---
+ T=TP
+ 30 CONTINUE
+*---
+* CONSTRUCT THE 3D TRACKS WHICH ENTER THE GEOMETRY THROUGH A LATERAL
+* SURFACE
+*---
+* length of the spatial integration interval
+ 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 40 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
+ NOM3D(N3D)=INDREG(NOM2D(1),IST)
+ H3D(N3D)=0.5D0
+ 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=2,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,IMU,1)
+ ENDDO
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+* incoming flux in + direction
+ NOMP=NOM3D(1)
+ INDP=KEYCUR(-NOMP)
+* incoming flux in - direction
+ NOMM=NOM3D(N3D)
+ INDM=KEYCUR(-NOMM)
+ STOT(1,IMU,II,1)=W3D*S(INDP,II)
+ STOT(N3D,IMU,II,2)=W3D*S(INDM,II)
+* regional sources
+ DO I=2,N3D-1
+ NOMI=NOM3D(I)
+ Q0=0.0D0
+ Q1=0.0D0
+ DO JF=1,NFUNL
+ IND=KEYFLX(NOMI,JF)
+ Q0=Q0+S(IND,II)*TRHAR(IMU,JF,1)
+ Q1=Q1+S(IND,II)*TRHAR(IMU,JF,2)
+ ENDDO
+ STOT(I,IMU,II,1)=W3D*Q0
+ STOT(I,IMU,II,2)=W3D*Q1
+ ENDDO
+* MCGFFAR: 'Source Term Isolation' Strategy turned on
+* MCGFFAS: 'Source Term Isolation' Strategy turned off
+* MCGFFAT: 'MOCC/MCI' Iterative Strategy
+ CALL SUBFFA(SUBSCH,K,KPN,M,N3D,H3D,NOM3D,NZON,
+ 1 SIGAL(0,II),STOT(1,IMU,II,1),STOT(1,IMU,II,2),
+ 2 NREG,NMU,NANI,NFUNL,TRHAR(1,1,1),KEYFLX,KEYCUR,
+ 3 IMU,PHI(1,II),B)
+ ENDIF
+ ENDDO
+ Z1=Z1P
+ I=IST
+* ---
+* negative polar sine track
+* ---
+ K=1
+ T=T2D(K-1)
+ TIN=1
+ N3D=1
+ NOM3D(N3D)=INDREG(NOM2D(1),IST)
+ H3D(N3D)=0.5D0
+ 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=2,N3D-1
+ H3D(II)=H3D(II)*VNORF(NOM3D(II),IANG,IMU,2)
+ ENDDO
+ DO II=1,NGEFF
+ IF(NCONV(II)) THEN
+* incoming flux in + direction
+ NOMP=NOM3D(1)
+ INDP=KEYCUR(-NOMP)
+* incoming flux in - direction
+ NOMM=NOM3D(N3D)
+ INDM=KEYCUR(-NOMM)
+ STOT(1,IMU,II,1)=W3D*S(INDP,II)
+ STOT(N3D,IMU,II,2)=W3D*S(INDM,II)
+* regional sources
+ DO I=2,N3D-1
+ NOMI=NOM3D(I)
+ Q0=0.0D0
+ Q1=0.0D0
+ DO JF=1,NFUNL
+ IND=KEYFLX(NOMI,JF)
+ Q0=Q0+S(IND,II)*TRHAR(IMU,JF,3)
+ Q1=Q1+S(IND,II)*TRHAR(IMU,JF,4)
+ ENDDO
+ STOT(I,IMU,II,1)=W3D*Q0
+ STOT(I,IMU,II,2)=W3D*Q1
+ ENDDO
+* MCGFFAR: 'Source Term Isolation' Strategy turned on
+* MCGFFAS: 'Source Term Isolation' Strategy turned off
+* MCGFFAT: 'MOCC/MCI' Iterative Strategy
+ CALL SUBFFA(SUBSCH,K,KPN,M,N3D,H3D,NOM3D,NZON,
+ 1 SIGAL(0,II),STOT(1,IMU,II,1),STOT(1,IMU,II,2),
+ 2 NREG,NMU,NANI,NFUNL,TRHAR(1,1,3),KEYFLX,KEYCUR,
+ 3 IMU,PHI(1,II),B)
+ ENDIF
+ ENDDO
+* ---
+ Z1=Z1P
+ 40 CONTINUE
+ ENDDO
+ ENDDO
+ DEALLOCATE(TRHAR,RHARM,STOT)
+ ENDIF
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(B,H3D,NOM3D,H2D,NOM2D,T2D)
+*
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 20:24:31 +0000