1 files changed, 117 insertions, 0 deletions
diff --git a/Dragon/src/USSSEK.f b/Dragon/src/USSSEK.f
new file mode 100644
index 0000000..849de6b
--- /dev/null
+++ b/Dragon/src/USSSEK.f
@@ -0,0 +1,117 @@
+*DECK USSSEK
+      SUBROUTINE USSSEK(NBNRS,NQT,LMOD,SIGR,CONRL,WEIGH,SIGL,PIJK,DIL)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the dilution matrix preserving the non-correlated collision
+* probability matrix in each subgroup. Use a fixed point iteration.
+*
+*Copyright:
+* Copyright (C) 2003 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
+* NBNRS   number of correlated fuel regions.
+* NQT     number of subgroups in admixed resonant isotope.
+* LMOD    moderator flag (=.TRUE. if all regions are containing the
+*         resonant isotopes; =.FALSE. if a moderator region exists).
+* SIGR    macroscopic total xs of the other isotopes.
+* CONRL   number density of the admixed resonant isotope.
+* WEIGH   multiband weights for the admixed resonant isotope.
+* SIGL    microscopic total xs of the admixed resonant isotope.
+* PIJK    non-correlated collision probability matrix.
+*
+*Parameters: input/output
+* DIL     estimate and converged value of the dilution matrix.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      LOGICAL LMOD
+      INTEGER NBNRS,NQT
+      REAL    SIGR(NBNRS),CONRL(NBNRS),WEIGH(NQT),SIGL(NQT),
+     1        PIJK(0:NBNRS,0:NBNRS),DIL(0:NBNRS,0:NBNRS)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: WORK
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(WORK(0:NBNRS,0:NBNRS,3))
+*
+      DEN=0.0
+      DO 20 I=0,NBNRS
+      DO 10 J=0,NBNRS
+      WORK(I,J,3)=PIJK(I,J)
+      DEN=MAX(DEN,ABS(PIJK(I,J)))
+   10 CONTINUE
+   20 CONTINUE
+      IF(LMOD) THEN
+         CALL ALINV(NBNRS,WORK(1,1,3),NBNRS+1,IER)
+      ELSE
+         CALL ALINV(NBNRS+1,WORK(0,0,3),NBNRS+1,IER)
+      ENDIF
+      IF(IER.NE.0) CALL XABORT('USSSEK: SINGULAR MATRIX(1).')
+      ITER=0
+   30 ITER=ITER+1
+      IF(ITER.GT.50) CALL XABORT('USSSEK: MAXIMUM NB. OF ITERATIONS.')
+      DO 45 I=0,NBNRS
+      DO 40 J=0,NBNRS
+      WORK(I,J,1)=0.0
+   40 CONTINUE
+   45 CONTINUE
+      DO 72 L=1,NQT
+      DO 55 I=0,NBNRS
+      DO 50 J=0,NBNRS
+      WORK(I,J,2)=DIL(I,J)
+   50 CONTINUE
+   55 CONTINUE
+      DO 60 I=1,NBNRS
+      WORK(I,I,2)=WORK(I,I,2)+SIGR(I)+CONRL(I)*SIGL(L)
+   60 CONTINUE
+      IF(LMOD) THEN
+         CALL ALINV(NBNRS,WORK(1,1,2),NBNRS+1,IER)
+      ELSE
+         CALL ALINV(NBNRS+1,WORK(0,0,2),NBNRS+1,IER)
+      ENDIF
+      IF(IER.NE.0) CALL XABORT('USSSEK: SINGULAR MATRIX(2).')
+      DO 71 I=0,NBNRS
+      DO 70 J=0,NBNRS
+      WORK(I,J,1)=WORK(I,J,1)+WEIGH(L)*WORK(I,J,2)
+   70 CONTINUE
+   71 CONTINUE
+   72 CONTINUE
+      ERR=0.0
+      DO 85 I=0,NBNRS
+      DO 80 J=0,NBNRS
+      ERR=MAX(ERR,ABS(PIJK(I,J)-WORK(I,J,1)))
+   80 CONTINUE
+   85 CONTINUE
+      IF(ERR.LT.1.0E-4*DEN) GO TO 110
+      IF(LMOD) THEN
+         CALL ALINV(NBNRS,WORK(1,1,1),NBNRS+1,IER)
+      ELSE
+         CALL ALINV(NBNRS+1,WORK(0,0,1),NBNRS+1,IER)
+      ENDIF
+      IF(IER.NE.0) CALL XABORT('USSSEK: SINGULAR MATRIX(3).')
+      DO 105 I=0,NBNRS
+      DO 100 J=0,NBNRS
+      DIL(I,J)=DIL(I,J)+WORK(I,J,3)-WORK(I,J,1)
+  100 CONTINUE
+  105 CONTINUE
+      GO TO 30
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+  110 DEALLOCATE(WORK)
+      RETURN
+      END