Initial commit from Polytechnique Montreal

1 files changed, 243 insertions, 0 deletions

diff --git a/Dragon/src/SYBUQV.f b/Dragon/src/SYBUQV.f
new file mode 100644
index 0000000..2da3291
--- /dev/null
+++ b/Dragon/src/SYBUQV.f
@@ -0,0 +1,243 @@
+*DECK SYBUQV
+      SUBROUTINE SYBUQV(ZZR,ZZI,NSURF,NREG,SIGT,MNA,LGFULL,PIJ)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the one-group collision, DP-0 leakage and DP-0 transmission
+* probabilities in a 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 (4 or 6).
+* NREG    number of regions.
+* SIGT    total macroscopic cross section.
+* MNA     number of angles.
+* LGFULL  void flad (=.TRUE. in voided regions).
+*
+*Parameters: output
+* PIJ     collision, DP-0 leakage and DP-0 transmission probabilities
+*         in lower triangular form.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER   ZZI(*),NSURF,NREG,MNA
+      REAL      ZZR(*),SIGT(NSURF:NSURF+NREG-1),
+     >          PIJ(0:(NREG+NSURF)*(NREG+NSURF+1)/2-1)
+      LOGICAL   LGFULL(NSURF:NSURF+NREG-1)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (MKI3=600,ZI30=0.785398164)
+      LOGICAL   LGEMPT
+      REAL, ALLOCATABLE, DIMENSION(:) :: POPTX
+      COMMON /BICKL3/BI3(0:MKI3),BI31(0:MKI3),BI32(0:MKI3),PAS3,XLIM3,L3
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(POPTX(2*NREG+4))
+*----
+*  CHECK FOR VOIDED REGIONS
+*----
+      LGEMPT=.FALSE.
+      DO 10 IR=1,NREG
+      LGEMPT=LGEMPT.OR.(.NOT.LGFULL(IR+NSURF-1))
+   10 CONTINUE
+*----
+*  INITIALIZATION
+*----
+      IZ0=0
+      IZR=4
+      PIJ(0:(NREG+NSURF)*(NREG+NSURF+1)/2-1)=0.0
+      IF(LGEMPT) THEN
+*        VOIDED REGION DETECTED.
+         DO 52 IA=1,MNA
+         MNT=ZZI(IZ0+1)
+         IZ0=IZ0+2
+         IZR=IZR+NSURF*2+1
+*
+*        LOOP OVER COMPLETE TRACKS.
+         DO 51 ITT=1,MNT
+         IZ0=IZ0+1
+         NH=ZZI(IZ0)
+         IZ0=IZ0+1
+         NX=ZZI(IZ0)
+*
+         IZ0=IZ0+1
+         ISW=ZZI(IZ0) ! WEST SIDE
+         JZI=IZ0
+         DO 50 ITX=1,NX
+         IZR=IZR+1
+         WITRAJ=ZZR(IZR)/4.0
+         IZ0=JZI
+*
+*        COMPUTE OPTICAL LENGTHS.
+         DO 20 IX=1,NH
+         IR=ZZI(IZ0+IX)
+         POPTX(IX)=ZZR(IZR+IX)*SIGT(IR)
+   20    CONTINUE
+         ISE=ZZI(IZ0+NH+1) ! EAST SIDE
+         IZR=IZR+NH
+*
+*        EXTERNAL IS IR=0
+         IRS=ISW
+         POP1=0.
+         DO 40 IXI=1,NH
+         IRT=IRS
+         IRS=ZZI(IZ0+IXI)
+*
+*        BEGINNING WITH IRS == JRS, P(i,i)=Taup  -(B3(0)-P3)
+         IIJ=IRS*(IRS+3)/2
+         CALL SYB33C(PPLUS,POPTX(IXI))
+         PIJ(IIJ)=PIJ(IIJ)+2.0*PPLUS*WITRAJ
+         TAUX=0.
+         DO 30 IXJ=IXI,NH
+         JRS=ZZI(IZ0+IXJ)
+         IIJ=MAX(IRT*(IRT+1)/2+JRS,JRS*(JRS+1)/2+IRT)
+         IF(IRT.LT.NSURF) THEN
+            CALL SYB32C(PPLUS,TAUX,POPTX(IXJ),2)
+            PIJ(IIJ)=PIJ(IIJ)+PPLUS*WITRAJ
+         ELSE
+            CALL SYB31C(PPLUS,TAUX,POP1,POPTX(IXJ))
+            IF(JRS.EQ.IRT) PPLUS=2.0*PPLUS
+            PIJ(IIJ)=PIJ(IIJ)+PPLUS*WITRAJ
+         ENDIF
+         TAUX=TAUX+POPTX(IXJ)
+   30    CONTINUE
+         IIJ=MAX(IRT*(IRT+1)/2+ISE,ISE*(ISE+1)/2+IRT)
+         IF(IRT.LT.NSURF) THEN
+            PIJ(IIJ)=PIJ(IIJ)+TABKI(3,TAUX)*WITRAJ
+         ELSE
+            CALL SYB32C(PPLUS,TAUX,POP1,2)
+            PIJ(IIJ)=PIJ(IIJ)+PPLUS*WITRAJ
+         ENDIF
+         POP1=POPTX(IXI)
+   40    CONTINUE
+         IZ0=IZ0+NH+1
+*
+*        COMPUTE REMAINING PSI FROM LAST REGION TO EAST SIDE.
+         IIJ=MAX(IRS*(IRS+1)/2+ISE,ISE*(ISE+1)/2+IRS)
+         CALL SYB32C(PPLUS,0.0,POP1,2)
+         PIJ(IIJ)=PIJ(IIJ)+PPLUS*WITRAJ
+   50    CONTINUE
+   51    CONTINUE
+   52    CONTINUE
+      ELSE
+*        NO VOIDED REGION DETECTED. FAST INTEGRATION METHOD.
+*        LOOP OVER THE ANGLE(FROM 0 TO PI/3).
+         DO 92 IA=1,MNA
+         MNT=ZZI(IZ0+1)
+         IZ0=IZ0+2
+         IZR=IZR+NSURF*2+1
+*
+*        LOOP OVER COMPLETE TRACKS.
+         DO 91 ITT=1,MNT
+         IZ0=IZ0+1
+         NH=ZZI(IZ0)
+         IZ0=IZ0+1
+         NX=ZZI(IZ0)
+*
+*        WEST SIDE
+         IZ0=IZ0+1
+         ISW=ZZI(IZ0)
+         JZI=IZ0
+         DO 90 ITX=1,NX
+         IZR=IZR+1
+         WITRAJ=ZZR(IZR)/4.0
+         JZR=IZR
+         IZ0=JZI
+         IF(WITRAJ.LT.0.0) CALL XABORT('SYBUQV: FAILURE 3.')
+         IF(NH.GT.2*NREG+4) CALL XABORT('SYBUQV: FAILURE 4.')
+*
+*        COMPUTE OPTICAL LENGTHS.
+         DO 60 IX=1,NH
+         IZ0=IZ0+1
+         IR=ZZI(IZ0)
+         IZR=IZR+1
+         SEG=ZZR(IZR)
+         POPTX(IX)=SEG*SIGT(IR)
+   60    CONTINUE
+*
+*        EXTERNAL IS IR=0
+         IZ0=JZI
+         IRS=ISW
+         JZ0=0
+         DO 80 IXI=1,NH
+         IRT=IRS
+         IZ00=IZ0
+         IZ0=IZ0+1
+         IRS=ZZI(IZ0)
+*
+*        PREVIOUS REGION J
+         JZ0=IZ00
+         PKI3=ZI30*WITRAJ
+         TAUX=0.
+*
+*        BEGINNING WITH IRS == JRS,P(i,i)=Taup  -(B3(0)-P3)
+         IIJ=IRS*(IRS+3)/2
+         PIJ(IIJ)=PIJ(IIJ)+POPTX(IXI)*WITRAJ
+         DO 70 IXJ=IXI,NH
+         JZ0=JZ0+1
+         JRS=ZZI(JZ0)
+         QIJ3=PKI3
+         TAUX=TAUX+POPTX(IXJ)
+         IF(TAUX.GE.XLIM3) THEN
+            PKI3=0.0
+         ELSE
+            K=NINT(TAUX*PAS3)
+            PKI3=(BI3(K)+TAUX*(BI31(K)+TAUX*BI32(K)))*WITRAJ
+         ENDIF
+         SIJ3=QIJ3-PKI3
+*
+*        BEGINNING OF PIJ : REGION IRS,JRS
+         IIJ=MAX(IRS*(IRS+1)/2+JRS,JRS*(JRS+1)/2+IRS)
+         PIJ(IIJ)=PIJ(IIJ)-SIJ3
+*
+*        REMAINING OF PIJ : Region IRS,JRS
+*        OR ... PSI IFF IXI=1, IRT IS THE WEST SIDE
+         IIJ=MAX(IRT*(IRT+1)/2+JRS,JRS*(JRS+1)/2+IRT)
+         PIJ(IIJ)=PIJ(IIJ)+SIJ3
+   70    CONTINUE
+*
+*        COMPUTE LEAKAGE AND TRANSMISSION PROBABILITIES PSI, PSS
+         JZ0=JZ0+1
+         ISE=ZZI(JZ0)
+         IIJ=IRS*(IRS+1)/2+ISE
+         PIJ(IIJ)=PIJ(IIJ)-PKI3
+         IIJ=MAX(IRT*(IRT+1)/2+ISE,ISE*(ISE+1)/2+IRT)
+         PIJ(IIJ)=PIJ(IIJ)+PKI3
+   80    CONTINUE
+*
+*        COMPUTE REMAINING PSI FROM LAST REGION TO EAST SIDE.
+         IZ0=IZ0+1
+         ISE=ZZI(JZ0)
+         IIJ=MAX(IRS*(IRS+1)/2+ISE,ISE*(ISE+1)/2+IRS)
+         PIJ(IIJ)=PIJ(IIJ)+ZI30*WITRAJ
+   90    CONTINUE
+   91    CONTINUE
+   92    CONTINUE
+         IIJ=NSURF*(NSURF+1)/2-1
+         DO 100 JR=1,NREG
+         IIJ=IIJ+NSURF+JR
+         PIJ(IIJ)=PIJ(IIJ)*2.0
+  100    CONTINUE
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(POPTX)
+      RETURN
+      END