Initial commit from Polytechnique Montreal

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