Initial commit from Polytechnique Montreal

1 files changed, 217 insertions, 0 deletions

diff --git a/Dragon/src/SYBALP.f b/Dragon/src/SYBALP.f
new file mode 100644
index 0000000..c3963bd
--- /dev/null
+++ b/Dragon/src/SYBALP.f
@@ -0,0 +1,217 @@
+*DECK SYBALP
+      SUBROUTINE SYBALP(NPIJ,MAXPTS,Y,SIG,NCOD,ALB,PIJ)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Pij calculation using the method of Kavenoky in 1D slab geometry.
+*
+*Copyright:
+* Copyright (C) 2005 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
+* NPIJ    number of regions.
+* MAXPTS  first dimension of matrix PIJ.
+* Y       abscissa array.
+* SIG     total cross section array.
+* NCOD    left and right type of boundary conditions (=1 void;
+*         =2 refl; =4 tran).
+* ALB     left and right albedos.
+*
+*Parameters: output
+* PIJ     reduced collision probability matrix.
+*
+*Reference:
+* A. Kavenoky, 'Calcul et utilisation des probabilites de premiere
+* collision pour les milieux heterogenes a une dimension: Les programmes
+* ALCOLL et CORTINA', note CEA-N-1077, Commissariat a l'energie
+* atomique, mars 1969.
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER NPIJ,MAXPTS,NCOD(2)
+      REAL Y(NPIJ+1),SIG(NPIJ),ALB(2),PIJ(MAXPTS,NPIJ)
+*----
+*  LOCAL VARIABLES
+*----
+      CHARACTER BC*8
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: AUXI,F2
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(AUXI(2*NPIJ,3),F2(2*NPIJ,2*NPIJ))
+*----
+*  SET THE BOUNDARY CONDITIONS
+*----
+      IF((NCOD(1).EQ.1).AND.(NCOD(2).EQ.1)) THEN
+         BC='VOID'
+         IF((ALB(1).NE.0.0).OR.(ALB(2).NE.0.0)) BC='ALBE'
+      ELSE IF((NCOD(1).EQ.4).AND.(NCOD(2).EQ.4)) THEN
+         BC='TRAN'
+      ELSE IF((ALB(1).EQ.0.0).AND.(ALB(2).EQ.0.0)) THEN
+         BC='VOID'
+      ELSE
+         BC='ALBE'
+      ENDIF
+*----
+*  COMPUTE THE PIJ MATRIX
+*----
+      PIJ(:MAXPTS,:NPIJ)=0.0
+      IF(BC.EQ.'VOID') THEN
+         DO 10 IP=1,NPIJ
+         AUXI(IP,1)=Y(IP+1)-Y(IP)
+         AUXI(IP,3)=MAX(1.0E-10,AUXI(IP,1)*SIG(IP))
+         AUXI(IP,2)=AUXI(IP,3)/AUXI(IP,1)
+   10    CONTINUE
+         DO 25 IP=1,NPIJ
+         CALL SYBRII(RIIP,1.0,0.0,AUXI(IP,3))
+         PIJ(IP,IP)=RIIP/AUXI(IP,2)**2
+         TAU0=0.0
+         DO 20 JP=IP+1,NPIJ
+         CALL SYBRIJ(RIJP,1.0,TAU0,AUXI(IP,3),AUXI(JP,3))
+         PIJ(IP,JP)=RIJP/(AUXI(IP,2)*AUXI(JP,2))
+         TAU0=TAU0+AUXI(JP,3)
+   20    CONTINUE
+   25    CONTINUE
+         DO 35 IP=1,NPIJ
+         DO 30 JP=IP,NPIJ
+         PIJ(JP,IP)=PIJ(IP,JP)
+   30    CONTINUE
+   35    CONTINUE
+      ELSE IF(BC.EQ.'TRAN') THEN
+         DO 40 IP=1,NPIJ
+         AUXI(IP,1)=Y(IP+1)-Y(IP)
+         AUXI(IP,3)=MAX(1.0E-10,AUXI(IP,1)*SIG(IP))
+         AUXI(IP,2)=AUXI(IP,3)/AUXI(IP,1)
+   40    CONTINUE
+         TAUCEL=0.0
+         DO 50 IP=1,NPIJ
+         TAUCEL=TAUCEL+AUXI(IP,3)
+   50    CONTINUE
+         M=-1
+   60    M=M+1
+         IF(M.GT.100) CALL XABORT('SYBALP: UNABLE TO CONVERGE(1).')
+         F2(:2*NPIJ,:2*NPIJ)=0.0
+         SMALL=0.0
+         DO 75 IP=1,NPIJ
+         CALL SYBRII(RIIP,1.0,M*TAUCEL,AUXI(IP,3))
+         CALL SYBRII(RIIM,-1.0,(M+1)*TAUCEL,AUXI(IP,3))
+         F2(IP,IP)=F2(IP,IP)+(RIIP+RIIM)/AUXI(IP,2)**2
+         SMALL=MAX(SMALL,ABS(F2(IP,IP)*AUXI(IP,2)))
+         TAU0=0.0
+         DO 70 JP=IP+1,NPIJ
+         CALL SYBRIJ(RIJP,1.0,M*TAUCEL+TAU0,AUXI(IP,3),AUXI(JP,3))
+         CALL SYBRIJ(RIJM,-1.0,(M+1)*TAUCEL-TAU0,AUXI(IP,3),AUXI(JP,3))
+         F2(IP,JP)=F2(IP,JP)+(RIJP+RIJM)/(AUXI(IP,2)*AUXI(JP,2))
+         TAU0=TAU0+AUXI(JP,3)
+         SMALL=MAX(SMALL,ABS(F2(IP,JP)*AUXI(JP,2)))
+   70    CONTINUE
+   75    CONTINUE
+         DO 85 IP=1,NPIJ
+         DO 80 JP=IP,NPIJ
+         PIJ(IP,JP)=PIJ(IP,JP)+F2(IP,JP)
+   80    CONTINUE
+   85    CONTINUE
+         IF(SMALL.LE.1.0E-6) GO TO 90
+         GO TO 60
+   90    DO 105 IP=1,NPIJ
+         DO 100 JP=IP,NPIJ
+         PIJ(JP,IP)=PIJ(IP,JP)
+  100    CONTINUE
+  105    CONTINUE
+      ELSE IF(BC.EQ.'ALBE') THEN
+         TAUCEL=0.0
+         DO 110 IP=1,NPIJ
+         AUXI(IP,1)=Y(IP+1)-Y(IP)
+         AUXI(IP,3)=MAX(1.0E-10,AUXI(IP,1)*SIG(IP))
+         AUXI(IP,2)=AUXI(IP,3)/AUXI(IP,1)
+         AUXI(2*NPIJ-IP+1,1)=AUXI(IP,1)
+         AUXI(2*NPIJ-IP+1,2)=AUXI(IP,2)
+         AUXI(2*NPIJ-IP+1,3)=AUXI(IP,3)
+         TAUCEL=TAUCEL+2.0*AUXI(IP,3)
+  110    CONTINUE
+         F2(:2*NPIJ,:2*NPIJ)=0.0
+         DO 125 IP=1,2*NPIJ
+         CALL SYBRII(RIIP,1.0,0.0,AUXI(IP,3))
+         F2(IP,IP)=F2(IP,IP)+RIIP/AUXI(IP,2)**2
+         TAU0=0.0
+         DO 120 JP=IP+1,2*NPIJ
+         CALL SYBRIJ(RIJP,1.0,TAU0,AUXI(IP,3),AUXI(JP,3))
+         F2(IP,JP)=F2(IP,JP)+RIJP/(AUXI(IP,2)*AUXI(JP,2))
+         F2(JP,IP)=F2(JP,IP)+RIJP/(AUXI(IP,2)*AUXI(JP,2))
+         TAU0=TAU0+AUXI(JP,3)
+  120    CONTINUE
+  125    CONTINUE
+         DO 135 IP=1,NPIJ
+         DO 130 JP=1,NPIJ
+         PIJ(IP,JP)=PIJ(IP,JP)+F2(IP,JP)+ALB(2)*F2(2*NPIJ+1-IP,JP)
+  130    CONTINUE
+  135    CONTINUE
+         M=0
+  140    M=M+1
+         IF(M.GT.100) CALL XABORT('UNABLE TO CONVERGE(2).')
+         F2(:2*NPIJ,:2*NPIJ)=0.0
+         SMALL=0.0
+         DO 155 IP=1,2*NPIJ
+         CALL SYBRII(RIIP,1.0,M*TAUCEL,AUXI(IP,3))
+         F2(IP,IP)=F2(IP,IP)+ALB(1)**M*RIIP/AUXI(IP,2)**2
+         SMALL=MAX(SMALL,ABS(F2(IP,IP)*AUXI(IP,2)))
+         TAU0=0.0
+         DO 150 JP=IP+1,2*NPIJ
+         CALL SYBRIJ(RIJP,1.0,M*TAUCEL+TAU0,AUXI(IP,3),AUXI(JP,3))
+         F2(IP,JP)=F2(IP,JP)+ALB(1)**M*RIJP/(AUXI(IP,2)*AUXI(JP,2))
+         CALL SYBRIJ(RIJM,-1.0,M*TAUCEL-TAU0,AUXI(JP,3),AUXI(IP,3))
+         F2(JP,IP)=F2(JP,IP)+ALB(1)**M*RIJM/(AUXI(IP,2)*AUXI(JP,2))
+         TAU0=TAU0+AUXI(JP,3)
+         SMALL=MAX(SMALL,ABS(F2(IP,JP)*AUXI(JP,2)))
+         SMALL=MAX(SMALL,ABS(F2(JP,IP)*AUXI(IP,2)))
+  150    CONTINUE
+  155    CONTINUE
+         DO 165 IP=1,NPIJ
+         DO 160 JP=1,NPIJ
+         PIJ(IP,JP)=PIJ(IP,JP)+ALB(2)**M*F2(IP,JP)+ALB(2)**(M-1)
+     1   *F2(2*NPIJ+1-IP,JP)
+  160    CONTINUE
+  165    CONTINUE
+         F2(:2*NPIJ,:2*NPIJ)=0.0
+         DO 175 IP=1,NPIJ
+         CALL SYBRII(RIIM,-1.0,M*TAUCEL,AUXI(IP,3))
+         F2(IP,IP)=F2(IP,IP)+ALB(1)**M*RIIM/AUXI(IP,2)**2
+         TAU0=0.0
+         DO 170 JP=IP+1,2*NPIJ
+         CALL SYBRIJ(RIJM,-1.0,M*TAUCEL-TAU0,AUXI(IP,3),AUXI(JP,3))
+         F2(IP,JP)=F2(IP,JP)+ALB(1)**M*RIJM/(AUXI(IP,2)*AUXI(JP,2))
+         CALL SYBRIJ(RIJP,1.0,M*TAUCEL+TAU0,AUXI(JP,3),AUXI(IP,3))
+         F2(JP,IP)=F2(JP,IP)+ALB(1)**M*RIJP/(AUXI(IP,2)*AUXI(JP,2))
+         TAU0=TAU0+AUXI(JP,3)
+         SMALL=MAX(SMALL,ABS(F2(IP,JP)*AUXI(JP,2)))
+         SMALL=MAX(SMALL,ABS(F2(JP,IP)*AUXI(IP,2)))
+  170    CONTINUE
+  175    CONTINUE
+         DO 185 IP=1,NPIJ
+         DO 180 JP=1,NPIJ
+         PIJ(IP,JP)=PIJ(IP,JP)+ALB(2)**M*F2(IP,JP)+ALB(2)**(M+1)*
+     1   F2(2*NPIJ+1-IP,JP)
+  180    CONTINUE
+  185    CONTINUE
+         IF(SMALL.LE.1.0E-6) GO TO 190
+         GO TO 140
+      ENDIF
+  190 DO 210 IP=1,NPIJ
+      DO 200 JP=1,NPIJ
+      PIJ(IP,JP)=PIJ(IP,JP)/AUXI(IP,1)
+  200 CONTINUE
+  210 CONTINUE
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(F2,AUXI)
+      RETURN
+      END