summaryrefslogtreecommitdiff
path: root/Dragon/src/SYBRN2.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SYBRN2.f')
-rw-r--r--Dragon/src/SYBRN2.f429
1 files changed, 429 insertions, 0 deletions
diff --git a/Dragon/src/SYBRN2.f b/Dragon/src/SYBRN2.f
new file mode 100644
index 0000000..63d5518
--- /dev/null
+++ b/Dragon/src/SYBRN2.f
@@ -0,0 +1,429 @@
+*DECK SYBRN2
+ SUBROUTINE SYBRN2 (NREG,NSURF,A,B,Z,IZ,VOL,SIGT,TRONC,PVS,PSS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the DP-1 leakage and transmission probabilities for an
+* heterogeneous non-sectorized square or rectangular cell. The tracks
+* are computed by SYBRTK.
+*
+*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
+* NREG number of regions in the cell.
+* NSURF number of surfaces.
+* A dimension of the external X sides.
+* B dimension of the external Y sides.
+* Z real tracking information.
+* IZ integer tracking information.
+* VOL volumes.
+* SIGT total macroscopic cross section.
+* TRONC voided block criterion.
+*
+*Parameters: output
+* PVS volume to surface probabilities:
+* XINF surfaces 1, 2 and 3; XSUP surfaces 4, 5 and 6;
+* YINF surfaces 7, 8 and 9; YSUP surfaces 10, 11 and 12.
+* PSS surface to surface probabilities in the following order:
+* PSS(i,j) is the probability from surface i to surface j.
+*
+*-----------------------------------------------------------------------
+*
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ INTEGER NREG,NSURF,IZ(*)
+ REAL A,B,Z(*),VOL(NREG),SIGT(NREG),TRONC,PVS(NREG,3*NSURF),
+ 1 PSS(3*NSURF,3*NSURF)
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (MKI3=600,MKI4=600,MKI5=600)
+ PARAMETER (PI=3.141592654,ZI30=0.785398164,ZI40=0.666666667)
+ LOGICAL ICARE
+ REAL KI3,KI4,KI5
+ INTEGER ISN(12,12)
+ REAL PBB(28)
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: COSINU
+*----
+* BICKLEY TABLES
+*----
+ 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
+*
+ SAVE ISN
+ DATA ISN/ 0, 0, 0, 4, 12, 0, 1, 9,-19, 1, 9, 19,
+ 1 0, 0, 0, 8, 16, 0, 5, 13,-23, 5, 13, 23,
+ 2 0, 0, 0, 0, 0, 28, 17, 21,-25,-17,-21,-25,
+ 3 4, 12, 0, 0, 0, 0, 1, 9, 19, 1, 9,-19,
+ 4 8, 16, 0, 0, 0, 0, 5, 13, 23, 5, 13,-23,
+ 5 0, 0, 28, 0, 0, 0,-17,-21,-25, 17, 21,-25,
+ 6 3, 11, 20, 3, 11,-20, 0, 0, 0, 2, 10, 0,
+ 7 7, 15, 24, 7, 15,-24, 0, 0, 0, 6, 14, 0,
+ 8 -18,-22,-27, 18, 22,-27, 0, 0, 0, 0, 0, 26,
+ 9 3, 11,-20, 3, 11, 20, 2, 10, 0, 0, 0, 0,
+ 1 7, 15,-24, 7, 15, 24, 6, 14, 0, 0, 0, 0,
+ 2 18, 22,-27,-18,-22,-27, 0, 0, 26, 0, 0, 0/
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(COSINU(2,NSURF))
+*----
+* INTEGRATION USING THE TRACKING
+*----
+ ICARE=IZ(1).EQ.2
+ ZERO=TRONC*(A+B)/(2.0*A*B)
+ AOB=A/B
+ PBB(:28)=0.0
+ PVS(:NREG,:12)=0.0
+ IZ0=2
+ IZR=4
+ DO 205 IA=1,IZ(2)
+ DO 20 I=1,NSURF
+ COSINU(1,I)=Z(IZR+1)
+ COSINU(2,I)=Z(IZR+2)
+ IZR=IZR+2
+ 20 CONTINUE
+ MNT=IZ(IZ0+1)
+ IZ0=IZ0+2
+ IZR=IZR+1
+ DO 200 IMNT=1,MNT
+ NH=IZ(IZ0+1)
+ NX=IZ(IZ0+2)
+ ISURF=IZ(IZ0+3)+1
+ JSURF=IZ(IZ0+NH+4)+1
+ DO 190 INX=1,NX
+ Z1=Z(IZR+1)
+ IZR=IZR+1
+ IF((ISURF.EQ.3).AND.(JSURF.EQ.2)) THEN
+ Z1=0.5*Z1/A
+ Z2=Z1*COSINU(1,ISURF)
+ Z3=Z1*COSINU(2,ISURF)
+ Z4=Z1*COSINU(1,ISURF)*COSINU(1,JSURF)
+ Z5=Z1*COSINU(2,ISURF)*COSINU(1,JSURF)
+ KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ DO 40 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+I)
+ IF(POP.LT.XLIM3) GO TO 30
+ IF(SIGTI.LE.ZERO) GO TO 50
+ PVS(III,7)=PVS(III,7)+2.0*KI3*Z1
+ PVS(III,8)=PVS(III,8)+2.0*KI4*Z3
+ GO TO 50
+ 30 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,7)=PVS(III,7)+2.0*TABKI(2,POP0)*Z(IZR+I)*Z1
+ PVS(III,8)=PVS(III,8)+2.0*KI3*Z(IZR+I)*Z3
+ ELSE
+ PVS(III,7)=PVS(III,7)+2.0*(KI3-WI3)*Z1
+ PVS(III,8)=PVS(III,8)+2.0*(KI4-WI4)*Z3
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 40 CONTINUE
+ 50 KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ K=0
+ DO 80 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+NH+1-I)
+ IF(POP.LT.XLIM3) GO TO 70
+ IF(SIGTI.LE.ZERO) GO TO 185
+ PVS(III,1)=PVS(III,1)+2.0*KI3*Z1*AOB
+ PVS(III,2)=PVS(III,2)+2.0*KI4*Z2*AOB
+ GO TO 185
+ 70 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+2.0*TABKI(2,POP0)*Z(IZR+NH+1-I)*Z1*AOB
+ PVS(III,2)=PVS(III,2)+2.0*KI3*Z(IZR+NH+1-I)*Z2*AOB
+ ELSE
+ PVS(III,1)=PVS(III,1)+2.0*(KI3-WI3)*Z1*AOB
+ PVS(III,2)=PVS(III,2)+2.0*(KI4-WI4)*Z2*AOB
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 80 CONTINUE
+ KI5=BI5(K)+POP*(BI51(K)+POP*BI52(K))
+ PBB(1)=PBB(1)+KI3*Z1
+ PBB(5)=PBB(5)+KI4*Z2
+ PBB(9)=PBB(9)+KI4*Z3
+ PBB(13)=PBB(13)+KI5*Z4
+ PBB(21)=PBB(21)+KI5*Z5
+ PBB(23)=PBB(23)+KI5*(Z1-Z5)
+ ELSE IF((ISURF.EQ.3).AND.(JSURF.EQ.4)) THEN
+ Z1=Z1/A
+ Z2=Z1*COSINU(2,ISURF)
+ Z3=Z1*COSINU(2,ISURF)*COSINU(2,JSURF)
+ KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ K=0
+ DO 100 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+I)
+ IF(POP.LT.XLIM3) GO TO 90
+ IF(SIGTI.LE.ZERO) GO TO 110
+ PVS(III,7)=PVS(III,7)+KI3*Z1
+ PVS(III,8)=PVS(III,8)+KI4*Z2
+ GO TO 110
+ 90 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,7)=PVS(III,7)+TABKI(2,POP0)*Z(IZR+I)*Z1
+ PVS(III,8)=PVS(III,8)+KI3*Z(IZR+I)*Z2
+ ELSE
+ PVS(III,7)=PVS(III,7)+(KI3-WI3)*Z1
+ PVS(III,8)=PVS(III,8)+(KI4-WI4)*Z2
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 100 CONTINUE
+ 110 KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ DO 130 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+NH+1-I)
+ IF(POP.LT.XLIM3) GO TO 120
+ IF(SIGTI.LE.ZERO) GO TO 185
+ PVS(III,7)=PVS(III,7)+KI3*Z1
+ PVS(III,8)=PVS(III,8)+KI4*Z2
+ GO TO 185
+ 120 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,7)=PVS(III,7)+TABKI(2,POP0)*Z(IZR+NH+1-I)*Z1
+ PVS(III,8)=PVS(III,8)+KI3*Z(IZR+NH+1-I)*Z2
+ ELSE
+ PVS(III,7)=PVS(III,7)+(KI3-WI3)*Z1
+ PVS(III,8)=PVS(III,8)+(KI4-WI4)*Z2
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 130 CONTINUE
+ KI5=BI5(K)+POP*(BI51(K)+POP*BI52(K))
+ PBB(2)=PBB(2)+2.0*KI3*Z1
+ PBB(6)=PBB(6)+2.0*KI4*Z2
+ PBB(14)=PBB(14)+2.0*KI5*Z3
+ PBB(26)=PBB(26)+2.0*KI5*(Z1-Z3)
+ ELSE IF((ISURF.EQ.1).AND.(JSURF.EQ.2)) THEN
+ Z1=Z1/B
+ Z2=Z1*COSINU(1,ISURF)
+ Z3=Z1*COSINU(1,ISURF)*COSINU(1,JSURF)
+ KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ DO 150 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+I)
+ IF(POP.LT.XLIM3) GO TO 140
+ IF(SIGTI.LE.ZERO) GO TO 160
+ PVS(III,1)=PVS(III,1)+KI3*Z1
+ PVS(III,2)=PVS(III,2)+KI4*Z2
+ GO TO 160
+ 140 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+TABKI(2,POP0)*Z(IZR+I)*Z1
+ PVS(III,2)=PVS(III,2)+KI3*Z(IZR+I)*Z2
+ ELSE
+ PVS(III,1)=PVS(III,1)+(KI3-WI3)*Z1
+ PVS(III,2)=PVS(III,2)+(KI4-WI4)*Z2
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 150 CONTINUE
+ 160 KI3=ZI30
+ KI4=ZI40
+ POP=0.0
+ K=0
+ DO 180 I=1,NH
+ III=IZ(IZ0+3+I)-NSURF+1
+ SIGTI=SIGT(III)
+ POP0=POP
+ POP=POP+SIGTI*Z(IZR+NH+1-I)
+ IF(POP.LT.XLIM3) GO TO 170
+ IF(SIGTI.LE.ZERO) GO TO 185
+ PVS(III,1)=PVS(III,1)+KI3*Z1
+ PVS(III,2)=PVS(III,2)+KI4*Z2
+ GO TO 185
+ 170 K=NINT(POP*PAS3)
+ WI3=BI3(K)+POP*(BI31(K)+POP*BI32(K))
+ WI4=BI4(K)+POP*(BI41(K)+POP*BI42(K))
+ IF(SIGTI.LE.ZERO) THEN
+ PVS(III,1)=PVS(III,1)+TABKI(2,POP0)*Z(IZR+NH+1-I)*Z1
+ PVS(III,2)=PVS(III,2)+KI3*Z(IZR+NH+1-I)*Z2
+ ELSE
+ PVS(III,1)=PVS(III,1)+(KI3-WI3)*Z1
+ PVS(III,2)=PVS(III,2)+(KI4-WI4)*Z2
+ ENDIF
+ KI3=WI3
+ KI4=WI4
+ 180 CONTINUE
+ KI5=BI5(K)+POP*(BI51(K)+POP*BI52(K))
+ PBB(4)=PBB(4)+2.0*KI3*Z1
+ PBB(8)=PBB(8)+2.0*KI4*Z2
+ PBB(16)=PBB(16)+2.0*KI5*Z3
+ PBB(28)=PBB(28)+2.0*KI5*(Z1-Z3)
+ ENDIF
+ 185 IZR=IZR+NH
+ 190 CONTINUE
+ IZ0=IZ0+NH+4
+ 200 CONTINUE
+ 205 CONTINUE
+*----
+* APPLY SYMMETRIES
+*----
+ IF(ICARE) THEN
+ PBB(1)=2.0*PBB(1)
+ PBB(5)=PBB(5)+PBB(9)
+ PBB(9)=PBB(5)
+ PBB(13)=2.0*PBB(13)
+ PBB(21)=PBB(21)+PBB(23)
+ PBB(23)=PBB(21)
+ PBB(4)=PBB(2)
+ PBB(8)=PBB(6)
+ PBB(16)=PBB(14)
+ PBB(28)=PBB(26)
+ DO 210 I=1,NREG
+ PVS(I,7)=PVS(I,7)+PVS(I,1)
+ PVS(I,8)=PVS(I,8)+PVS(I,2)
+ PVS(I,1)=PVS(I,7)
+ PVS(I,2)=PVS(I,8)
+ 210 CONTINUE
+ ENDIF
+ PBB(10)=PBB(6)
+ PBB(12)=PBB(8)
+ PBB(17)=PBB(9)
+ PBB(19)=PBB(5)
+ PBB(25)=PBB(13)
+ PBB(3)=PBB(1)*AOB
+ PBB(7)=PBB(9)*AOB
+ PBB(11)=PBB(5)*AOB
+ PBB(15)=PBB(13)*AOB
+ PBB(18)=PBB(19)*AOB
+ PBB(20)=PBB(17)*AOB
+ PBB(22)=PBB(23)*AOB
+ PBB(24)=PBB(21)*AOB
+ PBB(27)=PBB(25)*AOB
+*----
+* ORTHONORMALIZATION
+*----
+ Z1=Z(1)
+ Z2=Z(2)
+ Z3=Z(3)
+ Z4=Z(4)
+ DO 280 I=1,4
+ DEN0=PBB(I)
+ DEN1=PBB(4+I)
+ DEN2=PBB(8+I)
+ PBB(I)=Z1*Z1*DEN0
+ PBB(4+I)=Z1*(Z2*DEN1-Z3*DEN0)
+ PBB(8+I)=Z1*(Z2*DEN2-Z3*DEN0)
+ PBB(12+I)=Z2*Z2*PBB(12+I)-Z2*Z3*(DEN1+DEN2)+Z3*Z3*DEN0
+ PBB(24+I)=Z4*Z4*PBB(24+I)
+ 280 CONTINUE
+ DO 290 I=1,2
+ DEN1=PBB(16+I)
+ DEN2=PBB(18+I)
+ PBB(16+I)=Z1*Z4*DEN1
+ PBB(18+I)=Z1*Z4*DEN2
+ PBB(20+I)=(Z2*PBB(20+I)-Z3*DEN1)*Z4
+ PBB(22+I)=(Z2*PBB(22+I)-Z3*DEN2)*Z4
+ 290 CONTINUE
+ DO 300 J=1,NREG
+ DEN1=PVS(J,7)
+ DEN2=PVS(J,1)
+ PVS(J,7)=Z1*Z1*DEN1
+ PVS(J,1)=Z1*Z1*DEN2
+ PVS(J,8)=Z1*(Z2*PVS(J,8)-Z3*DEN1)
+ PVS(J,2)=Z1*(Z2*PVS(J,2)-Z3*DEN2)
+ 300 CONTINUE
+*
+ Z1=0.0
+ Z2=0.0
+ Z4=0.0
+ Z5=0.0
+ DO 510 J=1,NREG
+ X=4.0*VOL(J)
+ SIGTI=SIGT(J)
+ IF(SIGTI.GT.ZERO) THEN
+ Z1=Z1+PVS(J,7)
+ Z2=Z2+PVS(J,8)
+ Z4=Z4+PVS(J,1)
+ Z5=Z5+PVS(J,2)
+ X=X*SIGTI
+ ENDIF
+ PVS(J,7)=PVS(J,7)*A/X
+ PVS(J,8)=PVS(J,8)*A/X
+ PVS(J,1)=PVS(J,1)*B/X
+ PVS(J,2)=PVS(J,2)*B/X
+ 510 CONTINUE
+ IF((Z1.GT.TRONC).AND.(Z2.GT.TRONC).AND.(Z4.GT.TRONC).AND.(Z5.GT.
+ 1 TRONC)) THEN
+ Z1=(1.0-2.0*PBB(1)-PBB(2))/Z1
+ Z2=(-2.0*PBB(9)-PBB(10))/Z2
+ Z4=(1.0-2.0*PBB(3)-PBB(4))/Z4
+ Z5=(-2.0*PBB(11)-PBB(12))/Z5
+ DO 520 J=1,NREG
+ PVS(J,7)=PVS(J,7)*Z1
+ PVS(J,8)=PVS(J,8)*Z2
+ PVS(J,1)=PVS(J,1)*Z4
+ PVS(J,2)=PVS(J,2)*Z5
+ 520 CONTINUE
+ ENDIF
+*
+ DO 540 I=1,NREG
+ PVS(I,4)=PVS(I,1)
+ PVS(I,5)=PVS(I,2)
+ PVS(I,10)=PVS(I,7)
+ PVS(I,11)=PVS(I,8)
+ 540 CONTINUE
+ DO 560 JC=1,12
+ DO 550 IC=1,12
+ IB=ISN(IC,JC)
+ IF(IB.LT.0) THEN
+ PSS(IC,JC)=-PBB(-IB)
+ ELSE IF(IB.GT.0) THEN
+ PSS(IC,JC)=PBB(IB)
+ ELSE
+ PSS(IC,JC)=0.0
+ ENDIF
+ CONTINUE
+ 550 CONTINUE
+ 560 CONTINUE
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(COSINU)
+ RETURN
+ END