Initial commit from Polytechnique Montreal

1 files changed, 449 insertions, 0 deletions

diff --git a/Dragon/src/PIJXL3.f b/Dragon/src/PIJXL3.f
new file mode 100644
index 0000000..1d79b45
--- /dev/null
+++ b/Dragon/src/PIJXL3.f
@@ -0,0 +1,449 @@
+*DECK PIJXL3
+      SUBROUTINE PIJXL3(IPTRK,IPRT,NGRP,NANI,NBMIX,NPSYS,NRENOR,LEAKSW,
+     > XSSIGT,XSSIGW,NELPIJ,PIJ)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Calculation of the collision probabilities in EXCELL without producing
+* a tracking file. Based on subroutine XL3TRK in DRAGON 3.4.
+*
+*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): R. Roy and A. Hebert
+*
+*Parameters: input
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* IPRT    print flag (equal to zero for no print).
+* NGRP    number of energy groups.
+* NANI    number of Legendre orders (usually equal to one).
+* NBMIX   number of mixtures.
+* NPSYS   index set to zero if a group is not to be processed. Usually,
+*         NPSYS(I)=I.
+* NRENOR  normalization scheme for PIJ matrices.
+* LEAKSW  leakage flag (=.true. if neutron leakage through external
+*         boundary is present).
+* XSSIGT  total macroscopic cross sections ordered by mixture.
+* XSSIGW  P0 within-group scattering macroscopic cross sections
+*         ordered by mixture.
+* NELPIJ  number of elements in symmetrized pij matrix.
+*
+*Parameters: output
+* PIJ     reduced and symmetrized collision probabilities.
+*
+*-----------------------------------------------------------------------
+*
+      USE             GANLIB
+      IMPLICIT        NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      LOGICAL         LEAKSW
+      TYPE(C_PTR)     IPTRK
+      INTEGER         IPRT,NGRP,NANI,NBMIX,NPSYS(NGRP),NRENOR,NELPIJ
+      REAL            XSSIGT(0:NBMIX,NGRP),XSSIGW(0:NBMIX,NANI,NGRP),
+     >                PIJ(NELPIJ,NGRP)
+*----
+*  LOCAL VARIABLES
+*----
+      INTEGER         IOUT,NALB,NSTATE,ICPALL,ICPEND
+      PARAMETER      (IOUT=6,NALB=6,NSTATE=40,ICPALL=4,ICPEND=3)
+      INTEGER         NDIM  ,ISPEC ,NANGLE,NANGL ,ISYMM, NORE
+      INTEGER         NALBG ,NC    ,NTR   ,NTZ   ,ICL   ,NSOUT ,
+     >                ITGEO ,NRMV  ,NTY   ,LTRK  ,LINMAX,NUNK  ,
+     >                MAXR  ,NEXTGE,NTX   ,NCOR  ,NSUR  ,NTOTCL,
+     >                NVOL  ,NV    ,NS    ,IGRP  ,ISOUT ,ILONG ,
+     >                ITYPE ,INDPIJ,IIN   ,IBM   ,I     ,J     ,
+     >                NPIJ  ,NREG  ,NUNKMR
+      INTEGER         ISTATE(NSTATE),LCLSYM(3)
+      INTEGER         MXANGL,ICODE(NALB)
+      LOGICAL         SWVOID,SWNZBC,LSKIP
+      REAL            ALBOLD(NALB),EXTKOP(NSTATE),DENUSR, RCUTOF,
+     >                CUTOFX, FACT
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: MINDIM,MAXDIM,ICORD,INDEL,
+     > KEYMRG,MATALB,MATMRG,ICUR,INCR,NUMERO,MATRT
+      REAL, ALLOCATABLE, DIMENSION(:) :: REMESH,VOLSUR,VOLMRG,CONV,
+     > TRKBEG,TRKDIR,FFACT,ANGLES,DENSTY
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: SIGVOL,SIGTAL
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: LENGHT,VOLTRK,
+     > PSST,PSVT
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: DBLPIJ,PCSCT
+*----
+*  INTRINSIC FUNCTION FOR POSITION IN CONDENSE PIJ MATRIX
+*----
+      INTEGER INDPOS
+      INDPOS(I,J)=MAX(I,J)*(MAX(I,J)-1)/2+MIN(I,J)
+*----
+*  READ THE GEOMETRY INFORMATION STORED ON IPTRK
+*----
+      ISTATE(:NSTATE)=0
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+      NORE=ISTATE(8)
+      LTRK=ISTATE(9)+1
+      NANGLE=ISTATE(11)
+      ISYMM=ISTATE(12)
+      CALL LCMGET(IPTRK,'EXCELTRACKOP',EXTKOP)
+      CUTOFX=EXTKOP(1)
+      DENUSR=EXTKOP(2)
+      RCUTOF=EXTKOP(3)
+      CALL LCMSIX(IPTRK,'EXCELL',1)
+      ISTATE(:NSTATE)=0
+      CALL LCMGET(IPTRK,'STATE-VECTOR',ISTATE)
+      NDIM     =ISTATE(1)
+      NSUR    =-ISTATE(2)
+      NVOL     =ISTATE(3)
+      NTOTCL   =ISTATE(4)
+      MAXR     =ISTATE(5)
+      NUNK     =ISTATE(6)
+      NEXTGE   =ISTATE(7)
+*----
+*  Intrinsic symmetries used in geometry
+*  Use these to simplify tracking unless 
+*  NOSYMM tracking option activated
+*----
+      LCLSYM(1) =ISTATE(8)
+      LCLSYM(2) =ISTATE(9)
+      LCLSYM(3) =ISTATE(10)
+      IF(ISYMM .NE. 0) THEN
+        ISYMM=0
+        IF(NDIM .EQ. 2) THEN
+          IF(LCLSYM(1) .NE. 0) THEN
+*----
+*  X SYMMETRY
+*----
+            ISYMM=2
+          ENDIF
+          IF(LCLSYM(2) .NE. 0) THEN
+            IF(ISYMM .EQ. 0) THEN
+*----
+*  Y SYMMETRY
+*----
+              ISYMM=4
+            ELSE
+*----
+*  X AND Y SYMMETRY
+*----
+              ISYMM=8 
+            ENDIF
+          ENDIF
+*----
+*  X-Y DIAGONAL SYMMETRY
+*---- 
+        ELSE
+          IF(LCLSYM(1) .NE. 0) THEN
+*----
+*  X SYMMETRY
+*----
+            ISYMM=2
+          ENDIF
+          IF(LCLSYM(2) .NE. 0) THEN
+            IF(ISYMM .EQ. 0) THEN
+*----
+*  Y SYMMETRY
+*----
+              ISYMM=4
+            ELSE
+*----
+*  X AND Y SYMMETRY
+*----
+              ISYMM=8 
+            ENDIF
+          ENDIF 
+          IF(LCLSYM(3) .NE. 0) THEN
+*----
+*  Z SYMMETRY
+*----
+            ISYMM=ISYMM+16
+          ENDIF
+        ENDIF
+        IF(ISYMM .EQ. 0) ISYMM=1
+      ENDIF
+      ALLOCATE(MINDIM(NTOTCL),MAXDIM(NTOTCL),ICORD(NTOTCL),
+     > INDEL(4*NUNK),KEYMRG(NUNK),MATALB(NUNK))
+      ALLOCATE(REMESH(MAXR),VOLSUR(NUNK))
+      CALL LCMGET(IPTRK,'MINDIM      ',MINDIM)
+      CALL LCMGET(IPTRK,'MAXDIM      ',MAXDIM)
+      CALL LCMGET(IPTRK,'ICORD       ',ICORD )
+      CALL LCMGET(IPTRK,'INDEX       ',INDEL )
+      CALL LCMGET(IPTRK,'REMESH      ',REMESH)
+      CALL LCMGET(IPTRK,'KEYMRG      ',KEYMRG)
+      CALL LCMGET(IPTRK,'MATALB      ',MATALB)
+      CALL LCMGET(IPTRK,'VOLSUR      ',VOLSUR)
+      CALL LCMSIX(IPTRK,'EXCELL      ',2)
+      CALL LCMGET(IPTRK,'ICODE       ',ICODE )
+      CALL LCMGET(IPTRK,'ALBEDO      ',ALBOLD)
+*----
+*  VERIFY SYMMETRY AND STUDY TRACKING PARAMETERS. ARE THEY BASICALLY
+*  POSSIBLE ?
+*----
+      MXANGL=0
+      IF(LTRK .EQ. 1)THEN
+        NCOR= 1
+        IF(NDIM .EQ. 2) THEN
+          MXANGL=NANGLE
+          IF(ISYMM .GE. 2) THEN
+            NANGL = (NANGLE+1)/2
+          ELSE 
+            NANGL = NANGLE
+          ENDIF
+          IF( RCUTOF.GT.0.0 ) NCOR= 2
+        ELSE IF(NDIM .EQ. 3) THEN
+          IF(MOD(NANGLE,2) .EQ. 1)THEN
+            NANGLE=NANGLE+1
+            WRITE(IOUT,'(/31H MESS = ONLY EVEN # EQN ANGLES )')
+          ENDIF
+          IF(NANGLE .GT. 16)THEN
+            NANGLE=16
+            WRITE(IOUT,'(/31H MESS = 16 IS MAX # EQN ANGLES )')
+          ENDIF
+          MXANGL=(NANGLE * (NANGLE+2)) / 2 
+          IF(NEXTGE .EQ. 1) THEN
+            NANGL = (NANGLE * (NANGLE+2)) / 8
+          ELSE
+            IF(ISYMM .EQ. 8 .OR. ISYMM .EQ. 24) THEN
+              NANGL = (NANGLE * (NANGLE+2)) / 8 
+            ELSE IF(ISYMM .EQ. 2  .OR. ISYMM .EQ. 4 .OR. 
+     >              ISYMM .EQ. 18 .OR. ISYMM .EQ. 20 ) THEN
+              NANGL = (NANGLE * (NANGLE+2)) / 4 
+            ELSE
+              NANGL = (NANGLE * (NANGLE+2)) / 2 
+            ENDIF
+          ENDIF
+          IF(RCUTOF .GT. 0.0) NCOR= 4
+        ENDIF
+      ELSEIF( LTRK.EQ.2 )THEN
+        NCOR  = 1
+        MXANGL=NANGLE
+        IF( NDIM.EQ.2 )THEN
+           NANGL = NANGLE
+        ELSEIF( NDIM.EQ.3 )THEN
+           CALL XABORT('PIJXL3: *TSPC* NOT AVAILABLE FOR 3-D GEOMETRY')
+        ENDIF
+        CUTOFX= RCUTOF
+      ENDIF
+      IF( IPRT.GE.1 ) THEN
+         WRITE(IOUT,6002) NANGL,ISYMM,CUTOFX,DENUSR,RCUTOF
+      ENDIF
+      IF( IPRT.GT.1 .AND. NEXTGE.EQ.0 )THEN
+*
+*        IF PRINT REQUIRED AND OVERALL CARTESIAN GEOMETRY
+*        PRINT CARTESIAN REGION MAP
+         NTX= MAXDIM(1)-MINDIM(1)
+         NTY= MAXDIM(2)-MINDIM(2)
+         NTZ= MAXDIM(3)-MINDIM(3)
+         NTR=0
+         DO 103 ICL=4,NTOTCL
+            NTR= MAX(NTR,MAXDIM(ICL)-MINDIM(ICL)+1)
+  103    CONTINUE
+         CALL XELGPR(NDIM,NTX,NTY,NTZ,NTR,ISYMM,
+     >               NSUR,NVOL,NTOTCL,MINDIM,MAXDIM,
+     >               KEYMRG,INDEL,MATALB)
+      ENDIF
+      ALLOCATE(VOLTRK((NANGL+1)*NUNK))
+*
+      NV=  NVOL
+      NS= -NSUR
+      ALLOCATE(VOLMRG(NUNK),MATMRG(NUNK))
+      ITGEO=3
+      CALL XELCMP(   NS,    NV, VOLSUR, MATALB, KEYMRG,
+     >            NSOUT,  NREG, VOLMRG, MATMRG, ITGEO,ICODE)
+      NUNKMR= NREG+NSOUT+1
+      NPIJ= (NUNKMR*(NUNKMR+1))/2
+      IF( IPRT .GT. 1 ) WRITE(IOUT,6000) (NGRP*NPIJ/128)
+      ALLOCATE(DBLPIJ(NPIJ,NGRP))
+      IF( IPRT .GT. 1 ) WRITE(IOUT,6001)
+*
+*     ALLOCATE AND CHARGE TOTAL XS PER REGION
+      ALLOCATE(SIGTAL(NUNKMR,NGRP),SIGVOL(NREG,NGRP))
+*
+*  3) DO THE TRACKING OF THE EXACT GEOMETRY FOR *NEWT* OPTION.
+      IF( LTRK.NE.0 )THEN
+         NC= NTOTCL - 3
+         IF( IPRT.GE.1 )THEN
+            WRITE(IOUT,'(1H )')
+            IF( NC.EQ.0 )THEN
+               WRITE(IOUT,'(/38H NOW, TRACKING GEOMETRY WITH NO CYLIND,
+     >         2HER/)')
+            ELSEIF( NC.EQ.1 )THEN
+               WRITE(IOUT,'(/38H NOW, TRACKING GEOMETRY WITH ONE CYLIN,
+     >         3HDER/)')
+            ELSE
+               WRITE(IOUT,'(/28H NOW, TRACKING GEOMETRY WITH,I4,
+     >         10H CYLINDERS/)') NC
+            ENDIF
+         ENDIF
+         ALLOCATE(ICUR(NTOTCL),INCR(NTOTCL))
+         ALLOCATE(CONV(NTOTCL),TRKBEG(NTOTCL),TRKDIR(NTOTCL))
+*
+*  3.0)  WRITE FIRST RECORDS OF THE UNNORMALIZED TRACKING FILE
+         IF( LTRK.EQ.1 )THEN
+            LINMAX= 2*NVOL + 10
+         ELSE
+            LINMAX= 8*NANGL*(2*NVOL + 8)
+         ENDIF
+         ISPEC = LTRK-1
+         NALBG = 6
+         ALLOCATE(NUMERO(LINMAX))
+         ALLOCATE(LENGHT(LINMAX),ANGLES(3*MXANGL),DENSTY(MXANGL))
+*
+         NRMV=1
+         CALL XL3TI3( IPRT,   NANGLE, DENUSR, ISYMM,  ANGLES, DENSTY,
+     >                NTOTCL, NEXTGE, MAXR,   REMESH, LINMAX, RCUTOF,
+     >                NSUR,   NVOL,   INDEL,  MINDIM, MAXDIM, ICORD,
+     >                INCR,   ICUR,   TRKBEG, CONV,   TRKDIR, LENGHT, 
+     >                NUMERO, NPIJ,   NGRP,   SIGTAL, SWVOID, NORE,
+     >                NRMV,   VOLTRK, KEYMRG,-NSOUT,  NREG,   NPSYS,
+     >                DBLPIJ )
+*
+         CALL XL3NTR( IPRT, NDIM, ISPEC, NS, NV, NORE, VOLSUR, KEYMRG,
+     >                MATALB, NANGL, VOLTRK, DENSTY )
+*
+         CALL XL3SIG( NGRP, NBMIX,  XSSIGT, ALBOLD, NPSYS, NGRP, -NSOUT,
+     >                NREG, MATMRG, VOLMRG(NSOUT+2), SIGTAL, SIGVOL,
+     >                SWVOID, SWNZBC)
+*
+         NRMV=0
+         CALL XL3TI3( IPRT,   NANGLE, DENUSR, ISYMM,  ANGLES, DENSTY,
+     >                NTOTCL, NEXTGE, MAXR,   REMESH, LINMAX, RCUTOF,
+     >                NSUR,   NVOL,   INDEL,  MINDIM, MAXDIM, ICORD ,
+     >                INCR,   ICUR,   TRKBEG, CONV,   TRKDIR, LENGHT, 
+     >                NUMERO, NPIJ,   NGRP,   SIGTAL, SWVOID, NORE,
+     >                NRMV,   VOLTRK, KEYMRG, -NSOUT, NREG  , NPSYS,
+     >                DBLPIJ )
+*
+         CALL QIJCMP(NREG,-NSOUT,NPIJ,NGRP,NCOR,VOLMRG,SIGTAL,DBLPIJ,
+     >               NPSYS)
+*----
+*  RENORMALIZE ALL ISOTROPIC PROBS WITH VARIOUS OPTIONS
+*----
+         DO 2060 IGRP=1,NGRP
+         IF(NPSYS(IGRP).EQ.0) GO TO 2060
+         IF( NRENOR.EQ.1 )THEN
+*
+*           NORMALIZATION USING GELBARD SCHEME
+            CALL PIJRGL(IPRT,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP))
+         ELSEIF( NRENOR.EQ.2 )THEN
+*
+*           NORMALIZATION WORKING ON DIAGONAL COEFFICIENTS
+            CALL PIJRDG(NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP) )
+         ELSEIF( NRENOR.EQ.3 )THEN
+*
+*           NORMALIZATION WORKING ON WEIGHT FACTORS TO KEEP DIAG = 0.0
+            CALL PIJRNL(IPRT,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP))
+         ELSEIF( NRENOR .EQ. 4 )THEN  ! ATTENTION
+*
+*           NORMALIZATION WORKING ON WEIGHT FACTORS ADDITIVE (HELIOS)
+            CALL PIJRHL(IPRT,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP))
+         ENDIF
+         IF( IPRT.GE.ICPALL )THEN
+            WRITE(IOUT,'(1H )')
+            WRITE(IOUT,'(35H   COLLISION PROBABILITIES OUTPUT: ,
+     >                   35H *BEFORE* ALBEDO REDUCTION          )')
+            CALL PIJWPR(0,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP),
+     >      SIGVOL(1,IGRP),1)
+         ENDIF
+ 2060    CONTINUE
+*----
+*  ELIMINATION OF SURFACES FOR PIJ
+*----
+         IF( SWNZBC )THEN
+            ALLOCATE(PSST(NSOUT*NSOUT),PSVT(NSOUT*NREG))
+            ALLOCATE(MATRT(NSOUT))
+            CALL LCMLEN(IPTRK,'BC-REFL+TRAN',ILONG,ITYPE)
+            IF(ILONG.EQ.NSOUT) THEN
+              CALL LCMGET(IPTRK,'BC-REFL+TRAN',MATRT)
+            ELSE
+               DO 130 ISOUT=1,NSOUT
+                 MATRT(ISOUT)=ISOUT
+ 130           CONTINUE
+            ENDIF
+            DO 2080 IGRP=1,NGRP
+              IF(NPSYS(IGRP).EQ.0) GO TO 2080
+              CALL PIJABC(NREG,NSOUT,NPIJ,SIGTAL(1,IGRP),MATRT,
+     >                    DBLPIJ(1,IGRP),PSST,PSVT)
+ 2080       CONTINUE
+*
+            DEALLOCATE(MATRT)
+            DEALLOCATE(PSVT,PSST)
+         ENDIF
+*
+         ALLOCATE(FFACT(NREG))
+         DO 2090 IGRP=1,NGRP
+         IF(NPSYS(IGRP).EQ.0) GO TO 2090
+         IF( IPRT.GE.ICPEND )THEN
+            WRITE(IOUT,'(1H )')
+            WRITE(IOUT,'(35H   COLLISION PROBABILITIES OUTPUT: ,
+     >                   35H *AFTER* ALBEDO REDUCTION          )')
+            CALL PIJWPR(1,NREG,NSOUT,SIGTAL(1,IGRP),DBLPIJ(1,IGRP),
+     >      SIGVOL(1,IGRP),1)
+         ENDIF
+*----
+*  CHARGE PIJ MATRIX IN THE DRAGON SYMMETRIZED FORMAT
+*----
+         DO 160 IIN=1,NREG
+            IF(SIGTAL(NSOUT+IIN+1,IGRP).EQ.0.0) THEN
+               FFACT(IIN)=1.0
+            ELSE
+               FFACT(IIN)=1.0/SIGTAL(NSOUT+IIN+1,IGRP)
+            ENDIF
+  160    CONTINUE
+         CALL PIJD2R(NREG,NSOUT,DBLPIJ(1,IGRP),FFACT,.FALSE.,NELPIJ,
+     >               NPIJ,PIJ(1,IGRP))
+ 2090    CONTINUE
+         DEALLOCATE(FFACT)
+*
+         DEALLOCATE(DENSTY,ANGLES,LENGHT,NUMERO,TRKDIR,TRKBEG,CONV,
+     >   INCR,ICUR)
+      ENDIF
+      DEALLOCATE(INDEL,ICORD,MAXDIM,MINDIM,REMESH,DBLPIJ,SIGTAL,SIGVOL,
+     > VOLSUR,VOLTRK,KEYMRG,MATALB)
+*----
+*  CHECK IF SCATTERING REDUCTION IS REQUIRED
+*----
+      ALLOCATE(PCSCT(NREG,2*NREG))
+      DO 3000 IGRP=1,NGRP
+      IF(NPSYS(IGRP).EQ.0) GO TO 3000
+      LSKIP=.TRUE.
+      DO 200 IBM=1,NBMIX
+        LSKIP=LSKIP.AND.(XSSIGW(IBM,1,IGRP).EQ.0.0)
+  200 CONTINUE
+*----
+*  COMPUTE THE SCATTERING-REDUCED CP MATRICES
+*----
+      IF(.NOT.LSKIP) THEN
+        CALL PIJSMD(IPRT,NBMIX,NREG,MATMRG(NSOUT+2),VOLMRG(NSOUT+2),
+     >              XSSIGW(0,1,IGRP),XSSIGT(0,IGRP),LEAKSW,PIJ(1,IGRP),
+     >              PCSCT,1)
+        DO 220 I=1,NREG
+          FACT=VOLMRG(NSOUT+I+1)
+          DO 210 J=1,NREG
+            INDPIJ=INDPOS(I,J)
+            PIJ(INDPIJ,IGRP)=REAL(PCSCT(I,J))*FACT
+  210     CONTINUE
+  220   CONTINUE
+      ENDIF
+ 3000 CONTINUE
+      DEALLOCATE(PCSCT,VOLMRG,MATMRG)
+      RETURN 
+*----
+*  FORMAT
+*----
+ 6000 FORMAT(' *** SPACE REQUIRED FOR CP MATRICES = ',I10,' K ***')
+ 6001 FORMAT(' *** CP MATRICES ALLOCATED            ',10X,'   ***')
+ 6002 FORMAT(
+     > ' -----------------'/' RECOMPUTED PARAMETERS '/
+     > ' NANGL  =',I10  ,' (NUMBER OF TRACKING ANGLES)'/
+     > ' ISYMM  =',I10  ,' (TRACKING SYMMETRY FACTOR)'/
+     > ' CUTOFX =',F10.5,' (CUTOFF FOR TRACK LENGTH)'/
+     > ' DENS   =',F10.5,' (TRACK DENSITY)'/
+     > ' PCORN  =',F10.5,' (CORNER DUPLICATION DISTANCE)'/
+     > ' -----------------'/)
+      END