Initial commit from Polytechnique Montreal

1 files changed, 304 insertions, 0 deletions

diff --git a/Dragon/src/SYB003.f b/Dragon/src/SYB003.f
new file mode 100644
index 0000000..3ddfbf0
--- /dev/null
+++ b/Dragon/src/SYB003.f
@@ -0,0 +1,304 @@
+*DECK SYB003
+      SUBROUTINE SYB003 (NGEN,NPIJ,NPIS,SIGT2,SIGW2,IMPX,NCOUR,IWIGN,
+     1 IQUAD,XX,YY,NMC,RAYRE,MAIL,RZMAIL,PIJW,PISW,PSJW,PSSW)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the cellwise scattering-reduced collision, escape and
+* transmission probabilities in a 2-D Cartesian or hexagonal assembly
+* with Roth x 4 or Roth x 6 approximation.
+*
+*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
+* NGEN    total number of generating cells.
+* NPIJ    length of cellwise scattering-reduced collision probability
+*         matrices.
+* NPIS    length of cellwise scattering-reduced escape probability
+*         matrices (NPIS=NMC(NGEN+1)).
+* SIGT2   total macroscopic cross sections.
+* SIGW2   P0 within-group scattering macroscopic cross sections.
+* IMPX    print flag (equal to 0 for no print).
+* NCOUR   number of currents surrounding the cells (=4 Cartesian
+*         lattice; =6 hexagonal lattice).
+* IWIGN   type of cylinderization.
+* IQUAD   quadrature parameters.
+* XX      X-thickness of the generating cells.
+* YY      Y-thickness of the generating cells.
+* NMC     offset of the first volume in each generating cell.
+* RAYRE   radius of the tubes in each generating cell.
+* MAIL    offset of the first tracking information in each generating
+*         cell.
+* RZMAIL  real tracking information.
+*
+*Parameters: output
+* PIJW    cellwise scattering-reduced collision probability matrices.
+* PISW    cellwise scattering-reduced escape probability matrices.
+* PSJW    cellwise scattering-reduced collision probability matrices
+*         for incoming neutrons.
+* PSSW    cellwise scattering-reduced transmission probability
+*         matrices.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NGEN,NPIJ,NPIS,IMPX,NCOUR,IWIGN,IQUAD(4),NMC(NGEN+1),
+     1 MAIL(2,NGEN)
+      REAL SIGT2(NPIS),SIGW2(NPIS),XX(NGEN),YY(NGEN),RAYRE(NPIS),
+     1 RZMAIL(*),PIJW(NPIJ),PISW(NCOUR*NPIS),PSJW(NCOUR*NPIS),
+     2 PSSW(NGEN*NCOUR*NCOUR)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (PI=3.141592654)
+      LOGICAL LSKIP
+      INTEGER ISLR(4,4),ISLH(6,6)
+      REAL PBB(6),PSS(36)
+      REAL, ALLOCATABLE, DIMENSION(:) :: RAYR2,WORK
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: PIS,PSJ,PP
+*----
+*  DATA STATEMENTS
+*----
+      SAVE ISLR,ISLH
+      DATA ISLR/0,4,1,1,4,0,1,1,3,3,0,2,3,3,2,0/
+      DATA ISLH/0,1,2,3,2,1,1,0,1,2,3,2,2,1,0,1,2,3,
+     1          3,2,1,0,1,2,2,3,2,1,0,1,1,2,3,2,1,0/
+*
+      MR=IQUAD(4)
+      IPIJ=0
+      IPIS=0
+      IPSS=0
+      DO 240 JKG=1,NGEN
+      J1=NMC(JKG)
+      J2=NMC(JKG+1)-J1
+*----
+*  CYLINDERIZATION OPTIONS 
+*----
+      A=XX(JKG)
+      B=YY(JKG)
+      IB=MAIL(2,JKG)
+      RJ1=RAYRE(NMC(JKG+1))
+      SCALE1=1.0
+      SCALE2=1.0
+      ROUT=0.0
+      IF((NCOUR.EQ.4).AND.(IWIGN.EQ.1)) THEN
+*        ASKEW CYLINDERIZATION CARTESIAN.
+         RJ2=(A+B)/PI
+         SCALE1=(A*B-PI*RJ1**2)/(PI*RJ2**2-PI*RJ1**2)
+         ROUT=RJ2
+      ELSE IF((NCOUR.EQ.4).AND.(IWIGN.EQ.2)) THEN
+*        WIGNER CYLINDERIZATION CARTESIAN.
+         ROUT=SQRT(A*B/PI)
+      ELSE IF((NCOUR.EQ.4).AND.(IWIGN.EQ.3)) THEN
+*        SANCHEZ CYLINDERIZATION CARTESIAN.
+         SCALE2=SQRT(PI*A*B)/(A+B)
+         ROUT=SQRT(A*B/PI)
+      ELSE IF(IWIGN.EQ.1) THEN
+*        ASKEW CYLINDERIZATION HEXAGONAL.
+         RJ2=3.0*A/PI
+         SCALE1=(1.5*SQRT(3.0)*A*A-PI*RJ1**2)/(PI*RJ2**2-PI*RJ1**2)
+         ROUT=RJ2
+      ELSE IF(IWIGN.EQ.2) THEN
+*        WIGNER CYLINDERIZATION HEXAGONAL.
+         ROUT=SQRT(1.5*SQRT(3.0)/PI)*A
+      ELSE IF(IWIGN.EQ.3) THEN
+*        SANCHEZ CYLINDERIZATION HEXAGONAL.
+         SCALE2=SQRT(PI*SQRT(3.0)/6.0)
+         ROUT=SQRT(1.5*SQRT(3.0)/PI)*A
+      ENDIF
+      IF(ROUT.LE.RJ1) CALL XABORT('SYB003: CYLINDERIZATION ERROR.')
+*----
+*  COMPUTE THE REDUCED COLLISION PROBABILITY MATRIX
+*----
+      SURFA=0.5*PI*ROUT
+      ALLOCATE(PIS(J2,NCOUR),PSJ(NCOUR,J2),PP(J2,J2),RAYR2(J2+1))
+      DO 10 I=1,J2
+      RAYR2(I)=RAYRE(J1+I)
+   10 CONTINUE
+      RAYR2(J2+1)=ROUT
+      SIGT2(J1+J2)=SIGT2(J1+J2)*SCALE1
+      CALL SYBALC(J2,J2,RAYR2,SIGT2(J1+1),MR,0.0,RZMAIL(IB),PP)
+*
+      PSSX=0.0
+      RJ1=0.0
+      DO 50 I=1,J2
+      PISX=1.0
+      RJ2=RAYR2(I+1)**2
+      VV=PI*(RJ2-RJ1)
+      DO 20 J=1,J2
+      PISX=PISX-PP(I,J)*SIGT2(J1+J)
+   20 CONTINUE
+      PSSX=PSSX+PISX*SIGT2(J1+I)*VV/SURFA
+      IF(NCOUR.EQ.4) THEN
+         DEN1=2.0*(A+B)
+         PIS(I,1)=PISX*B/DEN1
+         PIS(I,2)=PISX*B/DEN1
+         PIS(I,3)=PISX*A/DEN1
+         PIS(I,4)=PISX*A/DEN1
+      ELSE
+         DO 30 IC=1,NCOUR
+         PIS(I,IC)=PISX/6.0
+   30    CONTINUE
+      ENDIF
+*----
+*  COMPUTE THE REDUCED COLLISION PROBABILITY MATRIX FOR INCOMING
+*  NEUTRONS
+*----
+      SURFA2=(0.5*PI*ROUT)/SCALE2
+      DO 40 IC=1,NCOUR
+      IF(I.LE.J2-1) THEN
+         PSJ(IC,I)=PISX*VV/SURFA2
+      ELSE IF(I.EQ.J2) THEN
+         PSJ(IC,I)=PISX*VV*SCALE1/SURFA2
+      ENDIF
+   40 CONTINUE
+      RJ1=RJ2
+   50 CONTINUE
+      DEALLOCATE(RAYR2)
+*----
+*  COMPUTE THE TRANSMISSION PROBABILITIES
+*----
+      PSSX=1.0-SCALE2*PSSX
+      IF(NCOUR.EQ.4) THEN
+         A=XX(JKG)
+         B=YY(JKG)
+         DEN1=MAX(2.0*B+A,2.0*A+B)
+         PBB(1)=B/DEN1
+         PBB(2)=MAX(A,2.0*A-B)/DEN1
+         PBB(3)=A/DEN1
+         PBB(4)=MAX(B,2.0*B-A)/DEN1
+      ELSE
+         PBB(1)=1.0/5.0
+         PBB(2)=1.0/5.0
+         PBB(3)=1.0/5.0
+      ENDIF
+      DO 65 JC=1,NCOUR
+      DO 60 IC=1,NCOUR
+      IF(NCOUR.EQ.4) THEN
+         IB=ISLR(IC,JC)
+      ELSE
+         IB=ISLH(IC,JC)
+      ENDIF
+      IF(IB.GT.0) THEN
+         PSS((JC-1)*NCOUR+IC)=PSSX*PBB(IB)
+      ELSE
+         PSS((JC-1)*NCOUR+IC)=0.0
+      ENDIF
+   60 CONTINUE
+   65 CONTINUE
+      SIGT2(J1+J2)=SIGT2(J1+J2)/SCALE1
+      IF(IMPX.GE.8) THEN 
+         CALL SYBPRX(1,NCOUR,J2,JKG,SIGT2(J1+1),SIGW2(J1+1),PP(1,1),
+     1   PIS(1,1),PSJ(1,1),PSS(1))
+      ENDIF
+*----
+*  CHECK IF SCATTERING REDUCTION IS REQUIRED
+*----
+      LSKIP=.TRUE.
+      DO 70 I=1,J2
+      LSKIP=LSKIP.AND.(SIGW2(J1+I).EQ.0.0)
+   70 CONTINUE
+*----
+*  SCATTERING REDUCTION IF LSKIP=.FALSE.
+*----
+      IF(LSKIP) THEN
+*        DO NOT PERFORM SCATTERING REDUCTION.
+         DO 85 I=1,J2
+         DO 80 J=1,J2-1
+         PIJW(IPIJ+(J-1)*J2+I)=PP(I,J)
+   80    CONTINUE
+         PIJW(IPIJ+(J2-1)*J2+I)=PP(I,J2)*SCALE1
+   85    CONTINUE
+         DO 95 I=1,J2
+         DO 90 JC=1,NCOUR
+         PISW(IPIS+(JC-1)*J2+I)=PIS(I,JC)
+         PSJW(IPIS+(I-1)*NCOUR+JC)=PSJ(JC,I)
+   90    CONTINUE
+   95    CONTINUE
+         DO 105 IC=1,NCOUR
+         DO 100 JC=1,NCOUR
+         PSSW(IPSS+(JC-1)*NCOUR+IC)=PSS((JC-1)*NCOUR+IC)
+  100    CONTINUE
+  105    CONTINUE
+      ELSE
+*        COMPUTE THE SCATTERING-REDUCED COLLISION AND ESCAPE MATRICES.
+         DO 120 I=1,J2
+         DO 110 J=1,J2-1
+         PIJW(IPIJ+(J-1)*J2+I)=-PP(I,J)*SIGW2(J1+J)
+  110    CONTINUE
+         PIJW(IPIJ+(J2-1)*J2+I)=-PP(I,J2)*SIGW2(J1+J2)*SCALE1
+         PIJW(IPIJ+(I-1)*J2+I)=1.0+PIJW(IPIJ+(I-1)*J2+I)
+  120    CONTINUE
+         CALL ALINV(J2,PIJW(IPIJ+1),J2,IER)
+         IF(IER.NE.0) CALL XABORT('SYB003: SINGULAR MATRIX.')
+         ALLOCATE(WORK(J2))
+         DO 175 I=1,J2
+         DO 130 K=1,J2
+         WORK(K)=PIJW(IPIJ+(K-1)*J2+I)
+  130    CONTINUE
+         DO 150 J=1,J2-1
+         WGAR=0.0
+         DO 140 K=1,J2
+         WGAR=WGAR+WORK(K)*PP(K,J)
+  140    CONTINUE
+         PIJW(IPIJ+(J-1)*J2+I)=WGAR
+  150    CONTINUE
+         WGAR=0.0
+         DO 155 K=1,J2
+         WGAR=WGAR+WORK(K)*PP(K,J2)
+  155    CONTINUE
+         PIJW(IPIJ+(J2-1)*J2+I)=WGAR*SCALE1
+         DO 170 JC=1,NCOUR
+         WGAR=0.0
+         DO 160 K=1,J2
+         WGAR=WGAR+WORK(K)*PIS(K,JC)
+  160    CONTINUE
+         PISW(IPIS+(JC-1)*J2+I)=WGAR
+  170    CONTINUE
+  175    CONTINUE
+         DEALLOCATE(WORK)
+*
+*        COMPUTE THE SCATTERING-REDUCED COLLISION PROBABILITY MATRIX
+*        FOR INCOMING NEUTRONS.
+         DO 200 IC=1,NCOUR
+         DO 190 J=1,J2
+         WGAR=PSJ(IC,J)
+         DO 180 K=1,J2
+         WGAR=WGAR+PSJ(IC,K)*SIGW2(J1+K)*PIJW(IPIJ+(J-1)*J2+K)
+  180    CONTINUE
+         PSJW(IPIS+(J-1)*NCOUR+IC)=WGAR
+  190    CONTINUE
+  200    CONTINUE
+*
+*        COMPUTE THE SCATTERING-REDUCED TRANSMISSION PROBABILITY MATRIX.
+         DO 230 IC=1,NCOUR
+         DO 220 JC=1,NCOUR
+         WGAR=PSS((JC-1)*NCOUR+IC)
+         DO 210 K=1,J2
+         WGAR=WGAR+PSJ(IC,K)*SIGW2(J1+K)*PISW(IPIS+(JC-1)*J2+K)
+  210    CONTINUE
+         PSSW(IPSS+(JC-1)*NCOUR+IC)=WGAR
+  220    CONTINUE
+  230    CONTINUE
+      ENDIF
+      DEALLOCATE(PP,PSJ,PIS)
+      IF(IMPX.GE.10) THEN 
+        CALL SYBPRX(2,NCOUR,J2,JKG,SIGT2(J1+1),SIGW2(J1+1),PIJW(IPIJ+1),
+     1  PISW(IPIS+1),PSJW(IPIS+1),PSSW(IPSS+1))
+      ENDIF
+      IPIJ=IPIJ+J2*J2
+      IPIS=IPIS+J2*NCOUR
+      IPSS=IPSS+NCOUR*NCOUR
+  240 CONTINUE
+      RETURN
+      END