summaryrefslogtreecommitdiff
path: root/Dragon/src/SYBUQ0.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SYBUQ0.f')
-rw-r--r--Dragon/src/SYBUQ0.f409
1 files changed, 409 insertions, 0 deletions
diff --git a/Dragon/src/SYBUQ0.f b/Dragon/src/SYBUQ0.f
new file mode 100644
index 0000000..905b5e5
--- /dev/null
+++ b/Dragon/src/SYBUQ0.f
@@ -0,0 +1,409 @@
+*DECK SYBUQ0
+ SUBROUTINE SYBUQ0(ZZR,ZZI,NSURF,NREG,SIGT,MNA,LGFULL,PSI,QSS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the one-group DP-1 leakage and transmission
+* probabilities in a sectorized Cartesian or hexagonal cell.
+*
+*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): A. Hebert
+*
+*Parameters: input
+* ZZR real tracking elements.
+* ZZI integer tracking elements.
+* NSURF number of surfaces.
+* NREG number of regions.
+* SIGT total macroscopic cross section.
+* MNA number of angles.
+* LGFULL voided region flag (=.TRUE. in voided regions).
+*
+*Parameters: output
+* PSI escape probabilities.
+* QSS transmission probabilities.
+*
+*-----------------------------------------------------------------------
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER ZZI(*),NSURF,NREG,MNA
+ REAL ZZR(*),SIGT(NREG),PSI(0:NSURF-1,3,NREG),
+ 1 QSS(9*NSURF*(NSURF-1)/2)
+ LOGICAL LGFULL(NREG)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (MKI3=600,MKI4=600,MKI5=600)
+ PARAMETER (ZI30=0.785398164,ZI40=0.666666667)
+ LOGICAL LGEMPT
+ REAL, ALLOCATABLE, DIMENSION(:) :: POPTX
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: COSINU
+ COMMON /BICKL3/BI3(0:MKI3),BI31(0:MKI3),BI32(0:MKI3),PAS3,XLIM3,L3
+ COMMON /BICKL4/BI4(0:MKI4),BI41(0:MKI4),BI42(0:MKI4),PAS4,XLIM4,L4
+ COMMON /BICKL5/BI5(0:MKI5),BI51(0:MKI5),BI52(0:MKI5),PAS5,XLIM5,L5
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(POPTX(2*NREG+4),COSINU(2,0:NSURF-1))
+*----
+* CHECK FOR VOIDED REGIONS
+*----
+ LGEMPT=.FALSE.
+ DO 10 IR=1,NREG
+ LGEMPT=LGEMPT.OR.(.NOT.LGFULL(IR))
+ 10 CONTINUE
+*----
+* INITIALIZATION
+*----
+ IZ0=0
+ IZR=4
+ DO 22 I=1,NREG
+ DO 21 J=1,3
+ DO 20 K=0,NSURF-1
+ PSI(K,J,I)=0.0
+ 20 CONTINUE
+ 21 CONTINUE
+ 22 CONTINUE
+ DO 30 I=1,9*NSURF*(NSURF-1)/2
+ QSS(I)=0.0
+ 30 CONTINUE
+*
+ IF(LGEMPT) THEN
+* VOIDED REGION DETECTED.
+* LOOP OVER THE ANGLE (FROM 0 TO PI/3).
+ DO 100 IA=1,MNA
+ MNT=ZZI(IZ0+1)
+ IZ0=IZ0+2
+ IZFINI=IZ0
+ IOF=IZR
+ DO 40 I=0,NSURF-1
+ IOF=IOF+1
+ COSINU(1,I)=ZZR(IOF)
+ IOF=IOF+1
+ COSINU(2,I)=ZZR(IOF)
+ 40 CONTINUE
+ IZR=IZR+NSURF*2+1
+*
+* LOOP OVER COMPLETE TRACKS.
+ DO 90 ITT=1,MNT
+ IZ0=IZ0+1
+ NH=ZZI(IZ0)
+ IZ0=IZ0+1
+ NX=ZZI(IZ0)
+*
+* WEST SIDE.
+ IZ0=IZ0+1
+ ISW=ZZI(IZ0)
+ IZDEBU=IZ0
+ IZFINI=IZ0+NH+1
+*
+ DO 80 ITX=1,NX
+ IZR=IZR+1
+ WITRAJ=ZZR(IZR)
+ IZ0=IZDEBU
+*
+* COMPUTE OPTICAL LENGTHS.
+ DO 50 IX=1,NH
+ IZ0=IZ0+1
+ IR=ZZI(IZ0)+1-NSURF
+ IZR=IZR+1
+ SEG=ZZR(IZR)
+ POPTX(IX)=SEG*SIGT(IR)
+ 50 CONTINUE
+*
+ IZ0=IZ0+1
+ ISE=ZZI(IZ0)
+*
+* EXTERNAL IS IXI=0
+ IZ0=IZDEBU
+ PKI3=ZI30*WITRAJ
+ PKI4=ZI40*WITRAJ
+ TAUX=0.
+ DO 60 IXI=1,NH
+ QIJ3=PKI3
+ QIJ4=PKI4
+ IZ0=IZ0+1
+ IRS=ZZI(IZ0)+1-NSURF
+ TAUP=TAUX
+ TAUX=TAUX+POPTX(IXI)
+ IF(TAUX.GE.XLIM3) THEN
+ PKI3=0.0
+ PKI4=0.0
+ ELSE
+ K=NINT(TAUX*PAS3)
+ PKI3=(BI3(K)+TAUX*(BI31(K)+TAUX*BI32(K)))*WITRAJ
+ PKI4=(BI4(K)+TAUX*(BI41(K)+TAUX*BI42(K)))*WITRAJ
+ ENDIF
+ SIJ3=QIJ3-PKI3
+ SIJ4=QIJ4-PKI4
+*
+* COMPUTE LEAKAGE PROBABILITIES.
+ IF(.NOT. LGFULL(IRS)) THEN
+ CALL SYB32C(SIJ3,TAUP,POPTX(IXI),2)
+ SIJ3=SIJ3*WITRAJ
+ CALL SYB32C(SIJ4,TAUP,POPTX(IXI),3)
+ SIJ4=SIJ4*WITRAJ
+ ENDIF
+ PSI(ISW,1,IRS)=PSI(ISW,1,IRS)+SIJ3
+ PSI(ISW,2,IRS)=PSI(ISW,2,IRS)+SIJ4*COSINU(2,ISW)
+ PSI(ISW,3,IRS)=PSI(ISW,3,IRS)+SIJ4*COSINU(1,ISW)
+ 60 CONTINUE
+*
+* EXTERNAL IS IXI=NH
+ IZ0=IZFINI
+ PKI3=ZI30*WITRAJ
+ PKI4=ZI40*WITRAJ
+ TAUX=0.
+ DO 70 IXI=NH,1,-1
+ QIJ3=PKI3
+ QIJ4=PKI4
+ IZ0=IZ0-1
+ IRS=ZZI(IZ0)+1-NSURF
+ TAUP=TAUX
+ TAUX=TAUX+POPTX(IXI)
+ IF(TAUX.GE.XLIM3) THEN
+ PKI3=0.0
+ PKI4=0.0
+ ELSE
+ K=NINT(TAUX*PAS3)
+ PKI3=(BI3(K)+TAUX*(BI31(K)+TAUX*BI32(K)))*WITRAJ
+ PKI4=(BI4(K)+TAUX*(BI41(K)+TAUX*BI42(K)))*WITRAJ
+ ENDIF
+ SIJ3=QIJ3-PKI3
+ SIJ4=QIJ4-PKI4
+ IF(.NOT. LGFULL(IRS)) THEN
+ CALL SYB32C(SIJ3,TAUP,POPTX(IXI),2)
+ SIJ3=SIJ3*WITRAJ
+ CALL SYB32C(SIJ4,TAUP,POPTX(IXI),3)
+ SIJ4=SIJ4*WITRAJ
+ ENDIF
+ PSI(ISE,1,IRS)=PSI(ISE,1,IRS)+SIJ3
+ PSI(ISE,2,IRS)=PSI(ISE,2,IRS)+SIJ4*COSINU(2,ISE)
+ PSI(ISE,3,IRS)=PSI(ISE,3,IRS)+SIJ4*COSINU(1,ISE)
+ 70 CONTINUE
+*
+ IF(TAUX.GE.XLIM5) THEN
+ PKI5=0.0
+ ELSE
+ K=NINT(TAUX*PAS5)
+ PKI5=(BI5(K)+TAUX*(BI51(K)+TAUX*BI52(K)))*WITRAJ
+ ENDIF
+*
+* COMPUTE TRANSMISSION PROBABILITIES.
+ IIJ=0
+ IS1=0
+ IS2=0
+ IF(ISW .LT. ISE) THEN
+ IIJ=((ISE-1)*ISE)/2+ISW
+ IS1=ISW
+ IS2=ISE
+ ELSE IF(ISE .LT. ISW) THEN
+ IIJ=((ISW-1)*ISW)/2+ISE
+ IS2=ISW
+ IS1=ISE
+ ELSE
+ CALL XABORT('SYBUQ0: IDENTICAL INCOMING AND OUTCOMING SUR'
+ 1 //'FACES(1).')
+ ENDIF
+ IIJ=IIJ*9
+ QSS(IIJ+1)=QSS(IIJ+1)+PKI3
+ QSS(IIJ+2)=QSS(IIJ+2)+PKI4*COSINU(2,IS1)
+ QSS(IIJ+3)=QSS(IIJ+3)+PKI4*COSINU(1,IS1)
+ QSS(IIJ+4)=QSS(IIJ+4)+PKI4*COSINU(2,IS2)
+ QSS(IIJ+7)=QSS(IIJ+7)+PKI4*COSINU(1,IS2)
+ QSS(IIJ+5)=QSS(IIJ+5)+PKI5*COSINU(2,IS1)*COSINU(2,IS2)
+ QSS(IIJ+6)=QSS(IIJ+6)+PKI5*COSINU(1,IS1)*COSINU(2,IS2)
+ QSS(IIJ+8)=QSS(IIJ+8)+PKI5*COSINU(2,IS1)*COSINU(1,IS2)
+ QSS(IIJ+9)=QSS(IIJ+9)+PKI5*COSINU(1,IS1)*COSINU(1,IS2)
+* END OF TRACK.
+ 80 CONTINUE
+ IZ0=IZFINI
+ 90 CONTINUE
+* END OF ANGLE
+ 100 CONTINUE
+ ELSE
+* NO VOIDED REGION DETECTED. FAST INTEGRATION METHOD.
+* LOOP OVER THE ANGLE (FROM 0 TO PI/3).
+ DO 170 IA=1,MNA
+ MNT=ZZI(IZ0+1)
+ IZ0=IZ0+2
+ IZFINI=IZ0
+ IOF=IZR
+ DO 110 I=0,NSURF-1
+ IOF=IOF+1
+ COSINU(1,I)=ZZR(IOF)
+ IOF=IOF+1
+ COSINU(2,I)=ZZR(IOF)
+ 110 CONTINUE
+ IZR=IZR+NSURF*2+1
+*
+* LOOP OVER COMPLETE TRACKS.
+ DO 160 ITT=1,MNT
+ IZ0=IZ0+1
+ NH=ZZI(IZ0)
+ IZ0=IZ0+1
+ NX=ZZI(IZ0)
+*
+* WEST SIDE.
+ IZ0=IZ0+1
+ ISW=ZZI(IZ0)
+ IZDEBU=IZ0
+ IZFINI=IZ0+NH+1
+ DO 150 ITX=1,NX
+ IZR=IZR+1
+ WITRAJ=ZZR(IZR)
+ IZ0=IZDEBU
+*
+* COMPUTE OPTICAL LENGTHS.
+ DO 120 IX=1,NH
+ IZ0=IZ0+1
+ IR=ZZI(IZ0)+1-NSURF
+ IZR=IZR+1
+ SEG=ZZR(IZR)
+ POPTX(IX)=SEG*SIGT(IR)
+ 120 CONTINUE
+*
+ IZ0=IZ0+1
+ ISE=ZZI(IZ0)
+*
+* EXTERNAL IS IXI=0
+ IZ0=IZDEBU
+ PKI3=ZI30*WITRAJ
+ PKI4=ZI40*WITRAJ
+ TAUX=0.
+ DO 130 IXI=1,NH
+ QIJ3=PKI3
+ QIJ4=PKI4
+ IZ0=IZ0+1
+ IRS=ZZI(IZ0)+1-NSURF
+ TAUX=TAUX+POPTX(IXI)
+ IF(TAUX.GE.XLIM3) THEN
+ PKI3=0.0
+ PKI4=0.0
+ ELSE
+ K=NINT(TAUX*PAS3)
+ PKI3=(BI3(K)+TAUX*(BI31(K)+TAUX*BI32(K)))*WITRAJ
+ PKI4=(BI4(K)+TAUX*(BI41(K)+TAUX*BI42(K)))*WITRAJ
+ ENDIF
+ SIJ3=QIJ3-PKI3
+ SIJ4=QIJ4-PKI4
+*
+* COMPUTE LEAKAGE PROBABILITIES.
+ PSI(ISW,1,IRS)=PSI(ISW,1,IRS)+SIJ3
+ PSI(ISW,2,IRS)=PSI(ISW,2,IRS)+SIJ4*COSINU(2,ISW)
+ PSI(ISW,3,IRS)=PSI(ISW,3,IRS)+SIJ4*COSINU(1,ISW)
+ 130 CONTINUE
+*
+* EXTERNAL IS IXI=NH
+ IZ0=IZFINI
+ PKI3=ZI30*WITRAJ
+ PKI4=ZI40*WITRAJ
+ TAUX=0.
+ DO 140 IXI=NH,1,-1
+ QIJ3=PKI3
+ QIJ4=PKI4
+ IZ0=IZ0-1
+ IRS=ZZI(IZ0)+1-NSURF
+ TAUX=TAUX+POPTX(IXI)
+ IF(TAUX.GE.XLIM3) THEN
+ PKI3=0.0
+ PKI4=0.0
+ ELSE
+ K=NINT(TAUX*PAS3)
+ PKI3=(BI3(K)+TAUX*(BI31(K)+TAUX*BI32(K)))*WITRAJ
+ PKI4=(BI4(K)+TAUX*(BI41(K)+TAUX*BI42(K)))*WITRAJ
+ ENDIF
+ SIJ3=QIJ3-PKI3
+ SIJ4=QIJ4-PKI4
+ PSI(ISE,1,IRS)=PSI(ISE,1,IRS)+SIJ3
+ PSI(ISE,2,IRS)=PSI(ISE,2,IRS)+SIJ4*COSINU(2,ISE)
+ PSI(ISE,3,IRS)=PSI(ISE,3,IRS)+SIJ4*COSINU(1,ISE)
+ 140 CONTINUE
+*
+ IF(TAUX.GE.XLIM5) THEN
+ PKI5=0.0
+ ELSE
+ K=NINT(TAUX*PAS5)
+ PKI5=(BI5(K)+TAUX*(BI51(K)+TAUX*BI52(K)))*WITRAJ
+ ENDIF
+*
+* COMPUTE TRANSMISSION PROBABILITIES.
+ IIJ=0
+ IS1=0
+ IS2=0
+ IF(ISW .LT. ISE) THEN
+ IIJ=((ISE-1)*ISE)/2+ISW
+ IS1=ISW
+ IS2=ISE
+ ELSE IF(ISE .LT. ISW) THEN
+ IIJ=((ISW-1)*ISW)/2+ISE
+ IS2=ISW
+ IS1=ISE
+ ELSE
+ CALL XABORT('SYBUQ0: IDENTICAL INCOMING AND OUTCOMING SUR'
+ 1 //'FACES(2).')
+ ENDIF
+ IIJ=IIJ*9
+ QSS(IIJ+1)=QSS(IIJ+1)+PKI3
+ QSS(IIJ+2)=QSS(IIJ+2)+PKI4*COSINU(2,IS1)
+ QSS(IIJ+3)=QSS(IIJ+3)+PKI4*COSINU(1,IS1)
+ QSS(IIJ+4)=QSS(IIJ+4)+PKI4*COSINU(2,IS2)
+ QSS(IIJ+7)=QSS(IIJ+7)+PKI4*COSINU(1,IS2)
+ QSS(IIJ+5)=QSS(IIJ+5)+PKI5*COSINU(2,IS1)*COSINU(2,IS2)
+ QSS(IIJ+6)=QSS(IIJ+6)+PKI5*COSINU(1,IS1)*COSINU(2,IS2)
+ QSS(IIJ+8)=QSS(IIJ+8)+PKI5*COSINU(2,IS1)*COSINU(1,IS2)
+ QSS(IIJ+9)=QSS(IIJ+9)+PKI5*COSINU(1,IS1)*COSINU(1,IS2)
+* END OF TRACK.
+ 150 CONTINUE
+ IZ0=IZFINI
+ 160 CONTINUE
+* END OF ANGLE.
+ 170 CONTINUE
+ ENDIF
+*----
+* NUMERICAL ORTHONORMALIZATION
+*----
+ Z1=ZZR(1)
+ Z2=ZZR(2)
+ Z3=ZZR(3)
+ Z4=ZZR(4)
+ DO 185 IRS=1,NREG
+ DO 180 ISE=0,NSURF-1
+ DEN0=PSI(ISE,1,IRS)
+ PSI(ISE,1,IRS)=0.25*Z1*Z1*DEN0
+ PSI(ISE,2,IRS)=0.25*Z1*Z4*PSI(ISE,2,IRS)
+ PSI(ISE,3,IRS)=0.25*Z1*(Z2*PSI(ISE,3,IRS)-Z3*DEN0)
+ 180 CONTINUE
+ 185 CONTINUE
+ IIJ=0
+ DO 190 IS=1,NSURF*(NSURF-1)/2
+ DEN0=QSS(IIJ+1)
+ DEN1=QSS(IIJ+3)
+ DEN2=QSS(IIJ+7)
+ QSS(IIJ+1)=Z1*Z1*DEN0
+ QSS(IIJ+3)=Z1*(Z2*DEN1-Z3*DEN0)
+ QSS(IIJ+7)=Z1*(Z2*DEN2-Z3*DEN0)
+ QSS(IIJ+9)=Z2*Z2*QSS(IIJ+9)-Z2*Z3*(DEN1+DEN2)+Z3*Z3*DEN0
+ DEN1=QSS(IIJ+2)
+ DEN2=QSS(IIJ+4)
+ QSS(IIJ+2)=Z1*Z4*DEN1
+ QSS(IIJ+4)=Z1*Z4*DEN2
+ QSS(IIJ+6)=(Z2*QSS(IIJ+6)-Z3*DEN2)*Z4
+ QSS(IIJ+8)=(Z2*QSS(IIJ+8)-Z3*DEN1)*Z4
+ QSS(IIJ+5)=Z4*Z4*QSS(IIJ+5)
+ IIJ=IIJ+9
+ 190 CONTINUE
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(COSINU,POPTX)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 17:31:32 +0000