Initial commit from Polytechnique Montreal

1 files changed, 285 insertions, 0 deletions

diff --git a/Dragon/src/PIJRNL.f b/Dragon/src/PIJRNL.f
new file mode 100644
index 0000000..262bae8
--- /dev/null
+++ b/Dragon/src/PIJRNL.f
@@ -0,0 +1,285 @@
+*DECK PIJRNL
+      SUBROUTINE PIJRNL(IPRT,NREG,NSOUT,SIGTAL,PROB)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Non-linear 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, G. Marleau
+*
+*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).
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT   NONE
+      INTEGER    IPRT,NREG,NSOUT,IUNOUT,NITMAX,NIT,
+     >           NUNKNO,IPRB,IPRF,IVOL,IDIA,IUNK,JUNK,IR,JR,
+     >           NSURC,NSURM,NVOLC,NVOLM,IP,NPR
+      REAL       SIGTAL(-NSOUT:NREG)
+      DOUBLE PRECISION PROB(*),EPSCON,TOTCON,WFSPAD
+      PARAMETER (IUNOUT=6, EPSCON=1.0E-8, NITMAX=10)
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IDL
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: CIJ,WSPACE,WFSP,
+     > WEIG
+C----
+C  SCRATCH STORAGE ALLOCATION
+C    CIJ   : MODIFIED CP PROB MATRIX
+C    WSPACE: NON-LINEAR SYSTEM MATRIX
+C    WFSP  : NON LINEAR SYSTEM SOLUTION
+C    IDL   : WSPACE DIAGONAL LOCATION
+C    WEIG  : NON-LINEAR WEIGHT
+C----
+      NPR=(NSOUT+NREG+1)*(NSOUT+NREG+2)/2
+      ALLOCATE(IDL(NSOUT+NREG+1))
+      ALLOCATE(CIJ(NPR),WSPACE(NPR),WFSP(-NSOUT:NREG),WEIG(-NSOUT:NREG))
+C
+C     CHARGE MATRIX "CIJ"
+      NUNKNO=NREG+NSOUT+1
+      IPRB= 0
+      IUNK= 0
+      IVOL= NSOUT*(NSOUT+1)/2
+      DO 20 IR   = -NSOUT, NREG
+         IUNK= IUNK+1
+         IF( IR.LT.0.OR.SIGTAL(IR).GT.0.0 )THEN
+            DO 10 JR= -NSOUT, IR-1
+               IPRB= IPRB+1
+               IF( JR.LT.0.OR.SIGTAL(JR).GT.0.0 )THEN
+                  CIJ(IPRB)= PROB(IPRB)
+               ELSE
+                  CIJ(IPRB)= 0.0D0
+               ENDIF
+   10       CONTINUE
+         ELSE
+            DO 15 JR= -NSOUT, IR-1
+               IPRB= IPRB+1
+               CIJ(IPRB)= 0.0D0
+   15       CONTINUE
+         ENDIF
+         IPRB= IPRB+1
+         IDL(IUNK)= IPRB
+         IF( IR.LT.0 )THEN
+            IVOL= IVOL+1
+            CIJ(IPRB)= PROB(IPRB)
+         ELSEIF( IR.GT.0 )THEN
+            IVOL= IVOL+IUNK-1
+            IF( SIGTAL(IR).GT.0.0 )THEN
+               CIJ(IPRB)= PROB(IPRB)
+            ELSE
+               CIJ(IPRB)= PROB(IVOL)
+            ENDIF
+         ELSE
+            IVOL= IVOL+1
+            CIJ(IPRB)= 1.0D0
+         ENDIF
+   20 CONTINUE
+C
+C     COPY MATRIX "CIJ" IN THE "WSPACE" ARRAY FOR INVERSION
+C          AND ADD TO THE DIAGONAL ALL TERMS OF A LINE
+      IPRB= 0
+      IUNK= 0
+      IDIA= 0
+      DO 50 IR   = -NSOUT, NREG
+         IUNK= IUNK+1
+         IDIA= IDIA+IUNK
+         WSPACE(IDIA)= CIJ(IDIA) + CIJ(IDIA)
+         DO 30 JR= -NSOUT, IR-1
+            IPRB= IPRB+1
+            WSPACE(IPRB)= CIJ(IPRB)
+            WSPACE(IDIA)= WSPACE(IDIA) + CIJ(IPRB)
+   30    CONTINUE
+         IPRB= IPRB+1
+         IPRF= IPRB
+         JUNK= IUNK
+         DO 40 JR=  IR+1 , NREG
+            IPRF= IPRF+JUNK
+            JUNK= JUNK+1
+            WSPACE(IDIA)= WSPACE(IDIA) + CIJ(IPRF)
+   40    CONTINUE
+   50 CONTINUE
+      IF( IPRT.GT.100 )THEN
+         WRITE(IUNOUT,8002)
+         IPRB= 0
+         DO 55 IR= -NSOUT, NREG
+         DO 52 JR= -NSOUT, IR
+            IPRB= IPRB+1
+            WRITE(IUNOUT,8003) IR, JR, CIJ(IPRB),
+     >                            WSPACE(IPRB),PROB(IPRB)
+   52    CONTINUE
+   55    CONTINUE
+      ENDIF
+C
+C     INVERSION OF THE INITIAL SYSTEM JACOBIAN MATRIX
+      CALL ALDDLF(NUNKNO,WSPACE,IDL)
+C
+C     INITIALISATION OF WEIGHTS
+      DO 60 IR=-NSOUT, NREG
+         WEIG(IR)=1.0D0
+   60 CONTINUE
+      WEIG(0)= 0.0D0
+C
+C     THE NON-LINEAR SYSTEM FOR WEIGHTS IS:
+C      F1(W1, W2, ... WN)= W1*(W1*C11+W2*C12+ ... +WN*C1N) - TRUE1
+C      F2(W1, W2, ... WN)= W2*(W1*C21+W2*C22+ ... +WN*C2N) - TRUE2
+C      ...
+C      FN(W1, W2, ... WN)= WN*(W1*CN1+W2*CN2+ ... +WN*CNN) - TRUEN
+C     FORMING THE SYSTEM USING WEIGHTS "WEIG" & CONTRIBUTIONS "CIJ"
+C
+C     MAIN ITERATION LOOP
+      DO 110 NIT=1,NITMAX
+C
+         IPRB= 0
+         IUNK= 0
+         IVOL= NSOUT*(NSOUT+1)/2
+         DO 90 IR=-NSOUT, NREG
+            IF( IR.LE.0 )THEN
+               IVOL= IVOL+1
+            ELSE
+               IVOL= IVOL+IUNK
+            ENDIF
+            IUNK= IUNK+1
+            WFSPAD= 0.0D0
+            DO 70 JR=-NSOUT, IR
+               IPRB= IPRB+1
+               WFSPAD=WFSPAD+WEIG(JR)*CIJ(IPRB)
+   70       CONTINUE
+            IPRF= IPRB
+            JUNK= IUNK
+            DO 80 JR= IR+1 , NREG
+               IPRF= IPRF+JUNK
+               JUNK= JUNK+1
+               WFSPAD=WFSPAD+WEIG(JR)*CIJ(IPRF)
+   80       CONTINUE
+            WFSP(IR)=WEIG(IR)*WFSPAD-PROB(IVOL)
+   90    CONTINUE
+         IF( IPRT.GT.100 )THEN
+            WRITE(IUNOUT,9000)
+            DO 92 IR= -NSOUT, NREG
+               WRITE(IUNOUT,9001) IR, WFSP(IR)
+   92       CONTINUE
+         ENDIF
+         CALL ALDDLS(NUNKNO,IDL,WSPACE,WFSP)
+C
+C        CALCULATIONS OF SQUARE DISTANCE BETWEEN 2 ITERATIONS
+C        AND UPDATING THE SOLUTION
+         TOTCON = 0.0D0
+         DO 100 IR=-NSOUT, NREG
+            TOTCON= TOTCON + WFSP(IR)**2
+            WEIG(IR)= WEIG(IR) - WFSP(IR)
+  100    CONTINUE
+         IF( IPRT.GT.100 )THEN
+            WRITE(IUNOUT,9004)
+            DO 102 IR= -NSOUT, NREG
+               WRITE(IUNOUT,9005) IR, WEIG(IR)
+  102       CONTINUE
+            WRITE(IUNOUT,'( 8H TOTCON: ,E15.7)') TOTCON
+         ENDIF
+C
+C        CONVERGENCE TEST
+         IF( TOTCON.LT.EPSCON )GO TO 120
+C
+  110 CONTINUE
+      WRITE(IUNOUT,'(35H PIJRNL: WEIGHTS NOT CONVERGED          )')
+  120 CONTINUE
+C
+C     RECOMPUTE WEIGHTS FOR VOID REGIONS
+      IPRB= (NSOUT+1)*(NSOUT+2)/2
+      IVOL= IPRB
+      IUNK= NSOUT+1
+      DO 220 IR= 1, NREG
+         IVOL= IVOL+IUNK
+         IUNK= IUNK+1
+         IF( SIGTAL(IR).EQ.0.0 )THEN
+            WFSPAD= 0.0D0
+            DO 200 JR=-NSOUT, IR
+               IPRB= IPRB+1
+               IF( JR.LT.0.OR.SIGTAL(JR).GT.0.0 )THEN
+                  WFSPAD=WFSPAD+WEIG(JR)*PROB(IPRB)
+               ENDIF
+  200       CONTINUE
+            IPRF= IPRB
+            JUNK= IUNK
+            DO 210 JR= IR+1 , NREG
+               IPRF= IPRF+JUNK
+               JUNK= JUNK+1
+               IF( JR.LT.0.OR.SIGTAL(JR).GT.0.0 )THEN
+                  WFSPAD=WFSPAD+WEIG(JR)*PROB(IPRF)
+               ENDIF
+  210       CONTINUE
+            WEIG(IR)=PROB(IVOL)/WFSPAD
+         ELSE
+            IPRB= IPRB+IUNK
+         ENDIF
+  220 CONTINUE
+C
+C     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)*WEIG(JR)
+            ENDIF
+  230    CONTINUE
+  240 CONTINUE
+C
+C     PRINT WEIGHT FACTORS IF REQUESTED
+      IF( IPRT .GE. 100 )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),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),IR=NVOLC,NVOLM, 1)
+            NVOLC = NVOLC + 10
+ 310     CONTINUE
+      ENDIF
+C----
+C  SCRATCH STORAGE DEALLOCATION
+C----
+      DEALLOCATE(WEIG,WFSP,WSPACE,CIJ)
+      DEALLOCATE(IDL)
+      RETURN
+C
+ 8002 FORMAT(//' S U R F / S U R F    C O N T R I B U T I O N S'//
+     >9X,'BEGIN S',6X,'END S ',11X,'CIJ.    ',11X, 'WSPACE  ',
+     >                         11X,'PROBS.  ')
+ 8003 FORMAT(6X,I10,5X,I10,5X,1P,E15.7,5X,E15.7,5X,E15.7 )
+ 9000 FORMAT(//' F U N C T I O N     V A L U E S'//
+     >9X,'VOL/SUR',6X,'VALUE')
+ 9001 FORMAT(6X,I10,5X,F10.4)
+ 9004 FORMAT(//' W E I G H T E D     V A L U E S'//
+     >9X,'VOL/SUR',6X,'VALUE')
+ 9005 FORMAT(6X,I10,5X,F10.4)
+      END