Initial commit from Polytechnique Montreal

1 files changed, 194 insertions, 0 deletions

diff --git a/Dragon/src/PIJRHL.f b/Dragon/src/PIJRHL.f
new file mode 100644
index 0000000..9e2040b
--- /dev/null
+++ b/Dragon/src/PIJRHL.f
@@ -0,0 +1,194 @@
+*DECK PIJRHL
+      SUBROUTINE PIJRHL(IPRT,NREG,NSOUT,SIGTAL,PROB)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* HELIOS type normalization of collision probs (CP).
+*
+*Copyright:
+* Copyright (C) 1994 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): R. Roy, E. Varin
+*
+*Parameters: input
+* IPRT    print level.
+* NREG    number of zones for geometry.
+* NSOUT   number of surfaces for geometry.
+* SIGTAL  albedo-sigt vector.
+*
+*Parameters: input/output
+* PROB    CP matrix for all types.
+*
+*References:
+* R. Roy and G. Marleau,
+* Normalization Techniques for CP Matrices,
+* CONF/PHYSOR-90, Marseille/France, V 2, P IX-40 (1990).
+* \\\\
+* E.A. Villarino, R.J.J. Stamm'ler, A.A. Ferri and J.J. Casal
+* HELIOS: Angularly Dependent Collision Probabilities.
+* Nucl.Sci.Eng. 112,16-31, 1992.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT   NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER   IPRT,NREG,NSOUT 
+      REAL      SIGTAL(-NSOUT:NREG)
+      DOUBLE PRECISION PROB(*)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER    IUNOUT,NITMAX,NIT,IPRINT,IR,JR,IP,IPRB,IND,I,J,CPTLB,
+     >           CPTAC,CTOT,NSURC,NSURM,NVOLC,NVOLM
+      LOGICAL    NOTCON
+      DOUBLE PRECISION NOM,DENOM,DMU,WFSPAD,WFSP,EPSCON,R1,R2,TOTCON,
+     >                 TMPCON
+      CHARACTER  HSMG*131
+      PARAMETER (IUNOUT=6, IPRINT=10, EPSCON=1.0E-6, NITMAX=20)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: CHI
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: WEIG
+*
+*----- INTRINSIC FUNCTION FOR POSITION IN CONDENSE PIJ MATRIX
+*
+      IND(I,J)=(MAX(I+NSOUT+1,J+NSOUT+1)*
+     >         (MAX(I+NSOUT+1,J+NSOUT+1)-1))/2
+     >         +MIN(I+NSOUT+1,J+NSOUT+1)
+*----
+*  SCRATCH STORAGE ALLOCATION
+*   WEIG : ADDITIVE WEIGHT
+*----
+      ALLOCATE(WEIG(-NSOUT:NREG,3),CHI(-NSOUT:NREG))
+*
+      NOTCON= .FALSE.
+      CPTLB = 3
+      CPTAC = 3
+      CTOT = CPTAC+CPTLB
+*
+*     INITIALISATION OF WEIGHTS
+      DO 60 IR=-NSOUT, NREG
+         WEIG(IR,1)=0.0D0
+         WEIG(IR,2)=0.5D0
+         WEIG(IR,3)=0.5D0
+   60 CONTINUE
+      DO 50 IR=-NSOUT, NREG
+         CHI(IR)= 1.0D0
+         IF( IR.GE.0.AND.SIGTAL(IR).EQ.0.0D0 )THEN
+            CHI(IR)= 0.0D0
+         ENDIF
+   50 CONTINUE
+*
+*     MAIN ITERATION LOOP
+      IF(IPRT.GT.2) WRITE(IUNOUT,'(A24)')
+     >       'ITER.     MU      ERROR '
+      DO 110 NIT=1,NITMAX
+*
+         DO 220 IR= -NSOUT, NREG
+            WFSPAD = PROB(IND(IR,0))
+     >                   + CHI(IR)*PROB(IND(IR,IR))*WEIG(IR,3)
+            WFSP = CHI(IR)*PROB(IND(IR,IR))
+            DO 200 JR=-NSOUT, NREG
+               WFSPAD = WFSPAD - CHI(JR)*WEIG(JR,3)*PROB(IND(IR,JR))
+               WFSP = WFSP + CHI(JR)*PROB(IND(IR,JR))
+  200       CONTINUE
+            IF(WFSP.NE.0.0) WEIG(IR,3) = WFSPAD / WFSP
+  220    CONTINUE
+*
+*        ACCELERATION TECHNIQUE
+         IF(  MOD(NIT-1,CTOT).GE.CPTAC )THEN
+            NOM   = 0.0D0
+            DENOM = 0.0D0
+            DO 10 IR=-NSOUT, NREG
+               R1= WEIG(IR,2) - WEIG(IR,1)
+               R2= WEIG(IR,3) - WEIG(IR,2)
+               NOM = NOM + R1*(R2-R1)
+               DENOM = DENOM + (R2-R1)*(R2-R1)
+   10       CONTINUE
+            IF(DENOM.EQ.0.0D0) THEN
+              DMU = 1.0D0
+            ELSE
+              DMU = - NOM / DENOM
+            ENDIF
+            IF( DMU.GT.50.0D0 .OR. DMU.LT.0.0D0 ) THEN
+              WRITE(HSMG,'(37HPIJRHL: PROBLEM OF ACCELERATION (DMU=,1P,
+     >        E11.4,2H).)') DMU
+              CALL XABORT(HSMG)
+            ENDIF
+            DO 20 IR=-NSOUT, NREG
+               WEIG(IR,3) = WEIG(IR,2) + DMU *
+     >                           (WEIG(IR,3) - WEIG(IR,2))
+               WEIG(IR,2) = WEIG(IR,1) + DMU *
+     >                           (WEIG(IR,2) - WEIG(IR,1))
+   20       CONTINUE
+         ELSE
+            DMU = 1.0D0
+         ENDIF
+*
+*        CALCULATIONS OF SQUARE DISTANCE BETWEEN 2 ITERATIONS
+*        AND UPDATING THE SOLUTION
+         TOTCON = 0.0D0
+         DO 100 IR=-NSOUT, NREG
+            TMPCON=ABS(WEIG(IR,3)-WEIG(IR,2))/WEIG(IR,3)
+            TOTCON=MAX(TMPCON,TOTCON)
+            WEIG(IR,1)= WEIG(IR,2)
+            WEIG(IR,2)= WEIG(IR,3)
+  100    CONTINUE
+         IF( IPRT.GT.2 ) WRITE(IUNOUT,'(I3,F9.5,E15.7)') NIT,DMU,TOTCON
+*
+*        CONVERGENCE TEST
+         IF( TOTCON.LT.EPSCON )GO TO 120
+*
+  110 CONTINUE
+      NOTCON=.TRUE.
+      WRITE(IUNOUT,'(35H PIJRHL: WEIGHTS NOT CONVERGED          )')
+  120 CONTINUE
+*
+*     RENORMALIZE "PIJ" SYMMETRIC MATRIX
+      IPRB = 0
+      DO 240 IR   = -NSOUT, NREG
+         DO 230 JR= -NSOUT, IR
+            IPRB= IPRB+1
+            IF( IR.NE.0.AND.JR.NE.0 )THEN
+                PROB(IPRB)=PROB(IPRB)*(WEIG(IR,1)+WEIG(JR,1))
+            ENDIF
+  230    CONTINUE
+  240 CONTINUE
+*
+*     PRINT WEIGHT FACTORS IF THERE IS A PROBLEM...
+      IF( NOTCON .OR. IPRT.GE.IPRINT )THEN
+         WRITE(IUNOUT,'(30H0 SURFACE WEIGHTS FACTORS                /)')
+         NSURC = -1
+         DO 300 IP  = 1, (9 +NSOUT) / 10
+            NSURM= MAX( -NSOUT, NSURC-9 )
+            WRITE(IUNOUT,'(10X,10( A5,    I6)/)')
+     >                     (' SUR ',-IR,IR= NSURC, NSURM, -1)
+            WRITE(IUNOUT,'(10H WEIGHT   ,10F11.5)')
+     >                   (WEIG(IR,1),IR=NSURC,NSURM,-1)
+            NSURC = NSURC - 10
+ 300     CONTINUE
+         WRITE(IUNOUT,'(30H0  VOLUME WEIGHTS FACTORS                /)')
+         NVOLC =  1
+         DO 310 IP  = 1, (9 + NREG) / 10
+            NVOLM= MIN( NREG, NVOLC+9 )
+            WRITE(IUNOUT,'(10X,10( A5 ,  I6)/)')
+     >                  (' VOL ',IR,IR=NVOLC,NVOLM, 1)
+            WRITE(IUNOUT,'(10H WEIGHT   ,10F11.5)')
+     >                   (WEIG(IR,1),IR=NVOLC,NVOLM, 1)
+            NVOLC = NVOLC + 10
+ 310     CONTINUE
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(CHI,WEIG)
+      RETURN
+      END