Initial commit from Polytechnique Montreal

1 files changed, 227 insertions, 0 deletions

diff --git a/Dragon/src/SYBUP1.f b/Dragon/src/SYBUP1.f
new file mode 100644
index 0000000..2173378
--- /dev/null
+++ b/Dragon/src/SYBUP1.f
@@ -0,0 +1,227 @@
+*DECK SYBUP1
+      SUBROUTINE SYBUP1(ZZR,ZZI,NSURF,NREG,SIGT,TRONC,A,B,IMPX,VOL,PIJ,
+     1 PIS,PSS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the one-group collision, DP-1 leakage and DP-1 transmission
+* probabilities in a Cartesian or hexagonal non-sectorized cell.
+*
+*Copyright:
+* Copyright (C) 2008 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.
+* TRONC   voided block criterion.
+* A       Cartesian dimension of the cell along the X axis or side of
+*         the hexagon.
+* B       Cartesian dimension of the cell along the Y axis.
+* IMPX    print flag.
+*
+*Parameters: output
+* VOL     volumes.
+* PIJ     cellwise reduced collision probability matrices.
+* PIS     cellwise reduced escape probability matrices.
+* PSS     cellwise reduced transmission probability matrices.
+*         PSS(i,j) is the probability from surface i to surface j.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER   ZZI(*),NSURF,NREG,IMPX
+      REAL      ZZR(*),SIGT(NREG),TRONC,A,B,VOL(NREG),PIJ(NREG,NREG),
+     1          PIS(NREG,3*NSURF),PSS(3*NSURF,3*NSURF)
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (SIGVID=1.0E-10,NCURR=3)
+      REAL      SURF(6)
+      REAL, ALLOCATABLE, DIMENSION(:) :: G,PIJS
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: LGFULL
+*----
+*  INLINE FUNCTIONS
+*----
+      INDPOS(I,J)=MAX(I,J)*(MAX(I,J)-1)/2+MIN(I,J)
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(G(NREG+6),LGFULL(NREG))
+*----
+*  CHECK FOR VOIDED REGIONS
+*----
+      DO 10 IR=1,NREG
+      VOL(IR)=ZZR(IR)
+      IF(VOL(IR).GT.0.) THEN
+         DR=SQRT(VOL(IR))
+      ELSE
+         DR=0.0
+      ENDIF
+      LGFULL(IR)=(SIGT(IR)*DR).GT.TRONC
+      IF(SIGT(IR).LE.SIGVID) SIGT(IR)=SIGVID
+   10 CONTINUE
+*----
+*  COMPUTE SYMMETRIZED CP MATRIX
+*----
+      IOFI=ZZI(1)
+      IOFR=ZZI(2)
+      ICARE=ZZI(3)
+      MNA=ZZI(4)
+      ALLOCATE(PIJS((NREG+NSURF)*(NREG+NSURF+1)/2))
+      CALL SYBUQV(ZZR(IOFR),ZZI(IOFI),NSURF,NREG,SIGT,MNA,LGFULL,PIJS)
+*----
+*  APPLY SYMMETRIES
+*----
+      IF(NSURF.EQ.4) THEN
+         SURF(1)=0.25*B
+         SURF(2)=0.25*B
+         SURF(3)=0.25*A
+         SURF(4)=0.25*A
+      ELSE
+         DO 20 JC=1,6
+         SURF(JC)=0.25*A
+   20    CONTINUE
+      ENDIF
+      DO 30 I=1,NSURF
+      G(I)=SURF(I)
+   30 CONTINUE
+      IF(ICARE.EQ.1) THEN
+*        RECTANGULAR CELL.
+         PIJS(2)=2.0*PIJS(2)
+         PIJS(5)=0.5*PIJS(5)
+         PIJS(9)=2.0*PIJS(9)
+         PIJS(4)=PIJS(5)
+         PIJS(7)=PIJS(5)
+         PIJS(8)=PIJS(5)
+         IOF=11
+         DO 50 I=1,NREG
+         G(4+I)=SIGT(I)*VOL(I)
+         SUM1=PIJS(IOF)+PIJS(IOF+1)
+         SUM2=PIJS(IOF+2)+PIJS(IOF+3)
+         PIJS(IOF)=SUM1
+         PIJS(IOF+1)=SUM1
+         PIJS(IOF+2)=SUM2
+         PIJS(IOF+3)=SUM2
+         DO 40 J=4,3+I
+         PIJS(IOF+J)=2.0*PIJS(IOF+J)
+   40    CONTINUE
+         IOF=IOF+4+I
+   50    CONTINUE
+      ELSE IF(ICARE.EQ.2) THEN
+*        SQUARE CELL.
+         PIJS(9)=2.0*PIJS(9)
+         PIJS(2)=PIJS(9)
+         PIJS(4)=PIJS(5)
+         PIJS(7)=PIJS(5)
+         PIJS(8)=PIJS(5)
+         IOF=11
+         DO 80 I=1,NREG
+         G(4+I)=SIGT(I)*VOL(I)
+         SUM=PIJS(IOF)+PIJS(IOF+1)+PIJS(IOF+2)+PIJS(IOF+3)
+         DO 60 J=0,3
+         PIJS(IOF+J)=SUM
+   60    CONTINUE
+         DO 70 J=4,3+I
+         PIJS(IOF+J)=4.0*PIJS(IOF+J)
+   70    CONTINUE
+         IOF=IOF+4+I
+   80    CONTINUE
+      ELSE IF(ICARE.EQ.3) THEN
+*        HEXAGONAL CELL.
+         PIJS(12)=2.0*PIJS(12)
+         PIJS(7)=PIJS(12)
+         PIJS(18)=PIJS(12)
+         PIJS(2)=PIJS(20)
+         PIJS(5)=PIJS(20)
+         PIJS(9)=PIJS(20)
+         PIJS(14)=PIJS(20)
+         PIJS(16)=PIJS(20)
+         PIJS(4)=PIJS(11)
+         PIJS(8)=PIJS(11)
+         PIJS(13)=PIJS(11)
+         PIJS(17)=PIJS(11)
+         PIJS(19)=PIJS(11)
+         IOF=22
+         DO 120 I=1,NREG
+         G(6+I)=SIGT(I)*VOL(I)
+         SUM=0.0
+         DO 90 J=0,5
+         SUM=SUM+PIJS(IOF+J)
+   90    CONTINUE
+         DO 100 J=0,5
+         PIJS(IOF+J)=SUM
+  100    CONTINUE
+         DO 110 J=6,5+I
+         PIJS(IOF+J)=6.0*PIJS(IOF+J)
+  110    CONTINUE
+         IOF=IOF+6+I
+  120    CONTINUE
+      ENDIF
+*----
+*  FIRST APPLY THE ORTHONORMALIZATION FACTOR
+*----
+      DO 130 I=1,(NSURF+NREG)*(NSURF+NREG+1)/2
+      PIJS(I)=PIJS(I)*ZZR(IOFR)*ZZR(IOFR)
+  130 CONTINUE
+*----
+*  PERFORM A DP-1 PIS AND PSS CALCULATION USING THE TRACKING
+*----
+      IF(NSURF.EQ.4) THEN
+         CALL SYBRN2(NREG,NSURF,A,B,ZZR(IOFR),ZZI(3),ZZR(1),SIGT,TRONC,
+     1   PIS,PSS)
+      ELSE IF(NSURF.EQ.6) THEN
+         CALL SYBHN2(NREG,NSURF,A,ZZR(IOFR),ZZI(3),ZZR(1),SIGT,TRONC,
+     1   PIS,PSS)
+      ENDIF
+*----
+*  VILLARINO-STAMM'LER NORMALIZATION
+*----
+      CALL SYBRHL(IMPX,NSURF,NREG,G,PIJS)
+      DO 160 IR=1,NREG
+      DO 150 IS=1,NSURF
+      ZNOR=G(IS)+G(NSURF+IR)
+      DO 140 IH=1,3
+      ISS=(IS-1)*3+IH
+      PIS(IR,ISS)=ZNOR*PIS(IR,ISS)
+  140 CONTINUE
+  150 CONTINUE
+  160 CONTINUE
+      DO 200 IS=1,NSURF
+      DO 190 JS=1,NSURF
+      ZNOR=G(IS)+G(JS)
+      DO 180 IH=1,3
+      ISS=(IS-1)*3+IH
+      DO 170 JH=1,3
+      JSS=(JS-1)*3+JH
+      PSS(ISS,JSS)=ZNOR*PSS(ISS,JSS)
+  170 CONTINUE
+  180 CONTINUE
+  190 CONTINUE
+  200 CONTINUE
+*----
+*  LOAD THE EURYDICE CP ARRAY
+*----
+      DO 220 I=1,NREG
+      DO 210 J=1,NREG
+      PIJ(I,J)=PIJS(INDPOS(NSURF+I,NSURF+J))/(VOL(I)*SIGT(I)*SIGT(J))
+  210 CONTINUE
+  220 CONTINUE
+      DEALLOCATE(PIJS)
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(LGFULL,G)
+      RETURN
+      END